/[escript]/trunk/escript/src/DataExpanded.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataExpanded.cpp

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

revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1748 by ksteube, Wed Sep 3 06:10:39 2008 UTC
# Line 20  Line 20 
20  #ifdef USE_NETCDF  #ifdef USE_NETCDF
21  #include <netcdfcpp.h>  #include <netcdfcpp.h>
22  #endif  #endif
23    #ifdef PASO_MPI
24    #include <mpi.h>
25    #endif
26    
27  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
28    
# Line 160  DataExpanded::~DataExpanded() Line 163  DataExpanded::~DataExpanded()
163  }  }
164    
165  DataAbstract*  DataAbstract*
166  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::getSlice(const DataArrayView::RegionType& region) const
167  {  {
168    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
169  }  }
170    
171  void  void
172  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
173                         const DataArrayView::RegionType& region)                         const DataArrayView::RegionType& region)
174  {  {
175    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
176    if (tempDataExp==0) {    if (tempDataExp==0) {
# Line 203  DataExpanded::setSlice(const DataAbstrac Line 206  DataExpanded::setSlice(const DataAbstrac
206  }  }
207    
208  void  void
209  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataArrayView& value)
210  {  {
211    //    //
212    // copy a single value to every data point in this object    // copy a single value to every data point in this object
# Line 213  DataExpanded::copy(const DataArrayView& Line 216  DataExpanded::copy(const DataArrayView&
216    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
217    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
218      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
219        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
220        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
221        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
222        // if using DOASSERT.        // if using DOASSERT.
223        getPointDataView().copy(m_data.index(i,j),value);        getPointDataView().copy(getPointOffset(i,j),value);
224      }      }
225    }    }
226  }  }
227    
228  void  void
229  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
230  {  {
231    //  
232    // first convert the numarray into a DataArray object    // extract the shape of the numarray
233    DataArray temp(value);    DataArrayView::ShapeType tempShape;
234      for (int i=0; i < value.getrank(); i++) {
235        tempShape.push_back(extract<int>(value.getshape()[i]));
236      }
237    
238      // get the space for the data vector
239      int len = DataArrayView::noValues(tempShape);
240      DataVector temp_data(len, 0.0, len);
241      DataArrayView temp_dataView(temp_data, tempShape);
242      temp_dataView.copy(value);
243    
244    //    //
245    // check the input shape matches this shape    // check the input shape matches this shape
246    if (!getPointDataView().checkShape(temp.getView().getShape())) {    if (!getPointDataView().checkShape(temp_dataView.getShape())) {
247      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(getPointDataView().createShapeErrorMessage(
248                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
249                          temp.getView().getShape()));                          temp_dataView.getShape()));
250    }    }
251    //    //
252    // now copy over the data    // now copy over the data
253    copy(temp.getView());    copy(temp_dataView);
254  }  }
255    
256  void  void
# Line 264  DataExpanded::toString() const Line 277  DataExpanded::toString() const
277    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
278    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
279      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
280        tempView.setOffset(m_data.index(i,j));        tempView.setOffset(getPointOffset(i,j));
281        stringstream suffix;        stringstream suffix;
282        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
283        temp << tempView.toString(suffix.str());        temp << tempView.toString(suffix.str());
# Line 287  DataArrayView Line 300  DataArrayView
300  DataExpanded::getDataPoint(int sampleNo,  DataExpanded::getDataPoint(int sampleNo,
301                             int dataPointNo)                             int dataPointNo)
302  {  {
303    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));
304    return temp;    return temp;
305  }  }
306    
# Line 311  DataExpanded::extractData(ifstream& arch Line 324  DataExpanded::extractData(ifstream& arch
324    return(m_data.extractData(archiveFile, noValues));    return(m_data.extractData(archiveFile, noValues));
325  }  }
326    
327  void  void
328  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
329    //    //
330    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
# Line 358  DataExpanded::copyToDataPoint(const int Line 371  DataExpanded::copyToDataPoint(const int
371               }               }
372             }             }
373           }           }
374       }       }
375    }    }
376  }  }
377  void  void
378  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
379    //    //
380    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
# Line 412  DataExpanded::copyToDataPoint(const int Line 425  DataExpanded::copyToDataPoint(const int
425               }               }
426             }             }
427           }           }
428       }       }
429    }    }
430  }  }
431  void  void
# Line 631  DataExpanded::setToZero(){ Line 644  DataExpanded::setToZero(){
644    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
645      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
646          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
647          for (int i=0; i<n ;++i) p[i]=0.;          for (i=0; i<n ;++i) p[i]=0.;
648      }      }
649    }    }
650  }  }
651    
652    /* Append MPI rank to file name if multiple MPI processes */
653    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
654      /* Make plenty of room for the mpi_rank number and terminating '\0' */
655      char *newFileName = (char *)malloc(strlen(fileName)+20);
656      strncpy(newFileName, fileName, strlen(fileName)+1);
657      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
658      return(newFileName);
659    }
660    
661  void  void
662  DataExpanded::dump(const std::string fileName) const  DataExpanded::dump(const std::string fileName) const
663  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");  
    #endif  
664     #ifdef USE_NETCDF     #ifdef USE_NETCDF
665     const int ldims=2+DataArrayView::maxRank;     const int ldims=2+DataArrayView::maxRank;
666     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
# Line 653  DataExpanded::dump(const std::string fil Line 671  DataExpanded::dump(const std::string fil
671     long dims[ldims];     long dims[ldims];
672     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
673     DataArrayView::ShapeType shape = getPointDataView().getShape();     DataArrayView::ShapeType shape = getPointDataView().getShape();
674       int mpi_iam=0, mpi_num=1;
675    
676     // netCDF error handler     // netCDF error handler
677     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
678     // Create the file.     // Create the file.
679     NcFile dataFile(fileName.c_str(), NcFile::Replace);  #ifdef PASO_MPI
680       MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
681       MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
682    #endif
683       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
684       NcFile dataFile(newFileName, NcFile::Replace);
685     // check if writing was successful     // check if writing was successful
686     if (!dataFile.is_valid())     if (!dataFile.is_valid())
687          throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");          throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
# Line 709  DataExpanded::dump(const std::string fil Line 733  DataExpanded::dump(const std::string fil
733     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
734     #endif     #endif
735  }  }
736    
737    void
738    DataExpanded::setTaggedValue(int tagKey,
739                                 const DataArrayView& value)
740    {
741      int numSamples = getNumSamples();
742      int numDataPointsPerSample = getNumDPPSample();
743      int sampleNo,dataPointNo, i;
744      DataArrayView& thisView=getPointDataView();
745      DataArrayView::ValueType::size_type n = thisView.noValues();
746      double* p,*in=&(value.getData()[0]);
747      
748      if (value.noValues() != n) {
749        throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
750      }
751    
752      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
753      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
754        if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
755            for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
756                p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
757                for (i=0; i<n ;++i) p[i]=in[i];
758            }
759        }
760      }
761    }
762    
763    void
764    DataExpanded::reorderByReferenceIDs(int *reference_ids)
765    {
766      int numSamples = getNumSamples();
767      DataArrayView& thisView=getPointDataView();
768      DataArrayView::ValueType::size_type n = thisView.noValues() * getNumDPPSample();
769      int sampleNo, sampleNo2,i;
770      double* p,*p2;
771      register double rtmp;
772      FunctionSpace fs=getFunctionSpace();
773    
774      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
775         const int id_in=reference_ids[sampleNo];
776         const int id=fs.getReferenceIDOfSample(sampleNo);
777         if (id!=id_in) {
778             bool matched=false;
779             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
780                  if (id == reference_ids[sampleNo2]) {
781                     p=&(m_data[getPointOffset(sampleNo,0)]);
782                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
783                     for (i=0; i<n ;i++) {
784                             rtmp=p[i];
785                             p[i]=p2[i];
786                             p2[i]=rtmp;
787                     }
788                     reference_ids[sampleNo]=id;
789                     reference_ids[sampleNo2]=id_in;
790                     matched=true;
791                     break;
792                  }
793             }
794             if (! matched) {
795                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
796             }
797         }
798       }
799    }
800    
801    
802  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1312  
changed lines
  Added in v.1748

  ViewVC Help
Powered by ViewVC 1.1.26