/[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 1796 by jfenwick, Wed Sep 17 01:45:46 2008 UTC revision 1827 by ksteube, Thu Oct 2 04:28:07 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15    #include "Data.h"
16  #include "DataExpanded.h"  #include "DataExpanded.h"
17  #include "DataException.h"  #include "DataException.h"
18  #include "DataConstant.h"  #include "DataConstant.h"
# Line 38  DataExpanded::DataExpanded(const boost:: Line 38  DataExpanded::DataExpanded(const boost::
38                             const FunctionSpace& what)                             const FunctionSpace& what)
39    : DataAbstract(what,DataTypes::shapeFromNumArray(value))    : DataAbstract(what,DataTypes::shapeFromNumArray(value))
40  {  {
 /*  DataTypes::ShapeType tempShape;  
   //  
   // extract the shape of the python numarray  
   for (int i=0; i<value.getrank(); i++) {  
     tempShape.push_back(extract<int>(value.getshape()[i]));  
   }*/  
41    //    //
42    // initialise the data array for this object    // initialise the data array for this object
43    initialise(what.getNumSamples(),what.getNumDPPSample());    initialise(what.getNumSamples(),what.getNumDPPSample());
# Line 56  DataExpanded::DataExpanded(const DataExp Line 50  DataExpanded::DataExpanded(const DataExp
50    : DataAbstract(other.getFunctionSpace(), other.getShape()),    : DataAbstract(other.getFunctionSpace(), other.getShape()),
51    m_data(other.m_data)    m_data(other.m_data)
52  {  {
 /*  
   //  
   // create the view for the data  
   DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());  
   setPointDataView(temp);  
   
   // keep local shape in sync  
   setShape(other.getPointDataView().getShape());*/  
53  }  }
54    
55  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
# Line 93  DataExpanded::DataExpanded(const DataTag Line 79  DataExpanded::DataExpanded(const DataTag
79    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
80      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
81        try {        try {
 //         getPointDataView().copy(getPointOffset(i,j),  
 //                                 other.getPointDataView(),  
 //                                 other.getPointOffset(i,j));  
82             DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),             DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
83                                  other.getVector(),                                  other.getVector(),
84                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
# Line 197  DataExpanded::~DataExpanded() Line 180  DataExpanded::~DataExpanded()
180  }  }
181    
182  DataAbstract*  DataAbstract*
183    DataExpanded::deepCopy()
184    {
185      return new DataExpanded(*this);
186    }
187    
188    
189    DataAbstract*
190  DataExpanded::getSlice(const DataTypes::RegionType& region) const  DataExpanded::getSlice(const DataTypes::RegionType& region) const
191  {  {
192    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
# Line 328  void Line 318  void
318  DataExpanded::initialise(int noSamples,  DataExpanded::initialise(int noSamples,
319                           int noDataPointsPerSample)                           int noDataPointsPerSample)
320  {  {
321      if (noSamples==0)     //retain the default empty object
322      {
323         return;
324      }
325    //    //
326    // resize data array to the required size    // resize data array to the required size
327    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
   
   //  
 //   // create the data view of the data array  
 //   DataArrayView temp(m_data.getData(),shape);  
 //   setPointDataView(temp);  
   
 //   // keep shape in sync  
 //   setShape(shape);  
328  }  }
329    
330  string  string
# Line 346  DataExpanded::toString() const Line 332  DataExpanded::toString() const
332  {  {
333    stringstream temp;    stringstream temp;
334    FunctionSpace fs=getFunctionSpace();    FunctionSpace fs=getFunctionSpace();
   //  
   // create a temporary view as the offset will be changed  
 //  DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());  
335    
336    int offset=0;    int offset=0;
337    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
# Line 362  DataExpanded::toString() const Line 345  DataExpanded::toString() const
345        }        }
346      }      }
347    }    }
348      string result=temp.str();
349      if (result.empty())
350      {
351        return "(data contains no samples)\n";
352      }
353    return temp.str();    return temp.str();
354  }  }
355    
# Line 372  DataExpanded::getPointOffset(int sampleN Line 360  DataExpanded::getPointOffset(int sampleN
360    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
361  }  }
362    
 // DataArrayView  
 // DataExpanded::getDataPoint(int sampleNo,  
 //                            int dataPointNo)  
 // {  
 //   DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));  
 //   return temp;  
 // }  
   
363  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
364  DataExpanded::getLength() const  DataExpanded::getLength() const
365  {  {
# Line 459  DataExpanded::copyToDataPoint(const int Line 439  DataExpanded::copyToDataPoint(const int
439       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
440             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
441       }       }
  //    DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
442       ValueType& vec=getVector();       ValueType& vec=getVector();
443       if (dataPointRank==0) {       if (dataPointRank==0) {
444           vec[0]=extract<double>(value[0]);           vec[0]=extract<double>(value[0]);
# Line 513  DataExpanded::copyAll(const boost::pytho Line 492  DataExpanded::copyAll(const boost::pytho
492    int dataPoint = 0;    int dataPoint = 0;
493    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
494      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
 //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
495        ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);        ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
496        if (dataPointRank==0) {        if (dataPointRank==0) {
497           vec[offset]=extract<double>(value[dataPoint]);           vec[offset]=extract<double>(value[dataPoint]);
# Line 525  DataExpanded::copyAll(const boost::pytho Line 503  DataExpanded::copyAll(const boost::pytho
503           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
504             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
505           vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);           vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
 //             dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);  
506             }             }
507           }           }
508         } else if (dataPointRank==3) {         } else if (dataPointRank==3) {
# Line 561  DataExpanded::symmetric(DataAbstract* ev Line 538  DataExpanded::symmetric(DataAbstract* ev
538    if (temp_ev==0) {    if (temp_ev==0) {
539      throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
540    }    }
 //  DataArrayView& thisView=getPointDataView();  
 //  DataArrayView& evView=ev->getPointDataView();  
541    ValueType& vec=getVector();    ValueType& vec=getVector();
542    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
543    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVector();
# Line 585  DataExpanded::nonsymmetric(DataAbstract* Line 560  DataExpanded::nonsymmetric(DataAbstract*
560    if (temp_ev==0) {    if (temp_ev==0) {
561      throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
562    }    }
 //   DataArrayView& thisView=getPointDataView();  
 //   DataArrayView& evView=ev->getPointDataView();  
563    ValueType& vec=getVector();    ValueType& vec=getVector();
564    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
565    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVector();
# Line 594  DataExpanded::nonsymmetric(DataAbstract* Line 567  DataExpanded::nonsymmetric(DataAbstract*
567    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
568    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
569      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
 //          DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),  
 //                                     evView,ev->getPointOffset(sampleNo,dataPointNo));  
570           DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),           DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
571                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
572      }      }
# Line 618  DataExpanded::trace(DataAbstract* ev, in Line 589  DataExpanded::trace(DataAbstract* ev, in
589    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
590    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
591      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
 /*         DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),  
                                     evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);*/  
592           DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),           DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
593                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
594      }      }
# Line 763  DataExpanded::dump(const std::string fil Line 732  DataExpanded::dump(const std::string fil
732     long dims[ldims];     long dims[ldims];
733     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
734     const DataTypes::ShapeType& shape = getShape();     const DataTypes::ShapeType& shape = getShape();
735     int mpi_iam=0, mpi_num=1;     int mpi_iam=getFunctionSpace().getDomain().getMPIRank();
736       int mpi_num=getFunctionSpace().getDomain().getMPISize();
737    #ifdef PASO_MPI
738       MPI_Status status;
739    #endif
740    
741    #ifdef PASO_MPI
742       /* Serialize NetCDF I/O */
743       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
744    #endif
745    
746     // netCDF error handler     // netCDF error handler
747     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
748     // Create the file.     // Create the file.
 #ifdef PASO_MPI  
    MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);  
    MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);  
 #endif  
749     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
750     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
751     // check if writing was successful     // check if writing was successful
# Line 787  DataExpanded::dump(const std::string fil Line 761  DataExpanded::dump(const std::string fil
761     if ( rank >0 ) {     if ( rank >0 ) {
762         dims[0]=shape[0];         dims[0]=shape[0];
763         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
764              throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
765     }     }
766     if ( rank >1 ) {     if ( rank >1 ) {
767         dims[1]=shape[1];         dims[1]=shape[1];
768         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
769              throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
770     }     }
771     if ( rank >2 ) {     if ( rank >2 ) {
772         dims[2]=shape[2];         dims[2]=shape[2];
773         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
774              throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
775     }     }
776     if ( rank >3 ) {     if ( rank >3 ) {
777         dims[3]=shape[3];         dims[3]=shape[3];
778         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
779              throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
780     }     }
781     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
782     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
# Line 811  DataExpanded::dump(const std::string fil Line 785  DataExpanded::dump(const std::string fil
785     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
786              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
787    
788     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )     if (getFunctionSpace().getNumSamples()>0)
789       {
790    
791         if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
792          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
793     const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();       const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
794     if (! (ids->put(ids_p,dims[rank+1])) )       if (! (ids->put(ids_p,dims[rank+1])) )
795          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
796         if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
    if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )  
797          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
798     if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
799          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
800       }
801    #ifdef PASO_MPI
802       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
803    #endif
804     #else     #else
805     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.");
806     #endif     #endif
807  }  }
808    
 // void  
 // DataExpanded::setTaggedValue(int tagKey,  
 //                              const DataArrayView& value)  
 // {  
 //   int numSamples = getNumSamples();  
 //   int numDataPointsPerSample = getNumDPPSample();  
 //   int sampleNo,dataPointNo, i;  
 //   DataTypes::ValueType::size_type n = getNoValues();  
 //   double* p,*in=&(value.getData()[0]);  
 //    
 //   if (value.noValues() != n) {  
 //     throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");  
 //   }  
 //  
 //   #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)  
 //   for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {  
 //     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {  
 //         for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {  
 //             p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);  
 //             for (i=0; i<n ;++i) p[i]=in[i];  
 //         }  
 //     }  
 //   }  
 // }  
   
809  void    void  
810  DataExpanded::setTaggedValue(int tagKey,  DataExpanded::setTaggedValue(int tagKey,
811             const DataTypes::ShapeType& pointshape,             const DataTypes::ShapeType& pointshape,

Legend:
Removed from v.1796  
changed lines
  Added in v.1827

  ViewVC Help
Powered by ViewVC 1.1.26