/[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 615 by elspeth, Wed Mar 22 02:12:00 2006 UTC revision 1748 by ksteube, Wed Sep 3 06:10:39 2008 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /* $Id$ */
3   ************************************************************  
4   *          Copyright 2006 by ACcESS MNRF                   *  /*******************************************************
5   *                                                          *   *
6   *              http://www.access.edu.au                    *   *           Copyright 2003-2007 by ACceSS MNRF
7   *       Primary Business: Queensland, Australia            *   *       Copyright 2007 by University of Queensland
8   *  Licensed under the Open Software License version 3.0    *   *
9   *     http://www.opensource.org/licenses/osl-3.0.php       *   *                http://esscc.uq.edu.au
10   *                                                          *   *        Primary Business: Queensland, Australia
11   ************************************************************   *  Licensed under the Open Software License version 3.0
12  */   *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  #include "DataExpanded.h"  #include "DataExpanded.h"
17  #include "DataException.h"  #include "DataException.h"
18  #include "DataConstant.h"  #include "DataConstant.h"
19  #include "DataTagged.h"  #include "DataTagged.h"
20    #ifdef USE_NETCDF
21    #include <netcdfcpp.h>
22    #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 154  DataExpanded::~DataExpanded() Line 162  DataExpanded::~DataExpanded()
162  {  {
163  }  }
164    
 void  
 DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)  
 {  
   if (getPointDataView().getRank()!=0) {  
     stringstream temp;  
     temp << "Error - Can only reshape Data with data points of rank 0. "  
          << "This Data has data points with rank: "  
          << getPointDataView().getRank();  
     throw DataException(temp.str());  
   }  
   //  
   // create the new DataBlocks2D data array, and a corresponding DataArrayView  
   DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape));  
   DataArrayView newView(newData.getData(),shape);  
   //  
   // Copy the original data to every value for the new shape  
   int i,j;  
   int nRows=m_data.getNumRows();  
   int nCols=m_data.getNumCols();  
   #pragma omp parallel for private(i,j) schedule(static)  
   for (i=0;i<nRows;i++) {  
     for (j=0;j<nCols;j++) {  
       // NOTE: An exception may be thown from this call if  
       // DOASSERT is on which of course will play  
       // havoc with the omp threads. Run single threaded  
       // if using DOASSERT.  
       newView.copy(newData.index(i,j),m_data(i,j));  
     }  
   }  
   // swap the new data array for the original  
   m_data.Swap(newData);  
   // set the corresponding DataArrayView  
   DataArrayView temp(m_data.getData(),shape);  
   setPointDataView(temp);  
 }  
   
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 213  DataExpanded::setSlice(const DataAbstrac Line 185  DataExpanded::setSlice(const DataAbstrac
185    if (getPointDataView().getRank()!=region.size()) {    if (getPointDataView().getRank()!=region.size()) {
186      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
187    }    }
188    if (tempDataExp->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
189      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (value->getPointDataView().createShapeErrorMessage(
190          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape));
191    }    }
# Line 234  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 244  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 289  string Line 271  string
271  DataExpanded::toString() const  DataExpanded::toString() const
272  {  {
273    stringstream temp;    stringstream temp;
274      FunctionSpace fs=getFunctionSpace();
275    //    //
276    // create a temporary view as the offset will be changed    // create a temporary view as the offset will be changed
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 << "(" << i << "," << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
283        temp << tempView.toString(suffix.str());        temp << tempView.toString(suffix.str());
284        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
285          temp << endl;          temp << endl;
# Line 317  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 327  DataExpanded::getLength() const Line 310  DataExpanded::getLength() const
310    return m_data.size();    return m_data.size();
311  }  }
312    
 void  
 DataExpanded::setRefValue(int ref,  
                           const DataArray& value)  
 {  
   //  
   // Get the number of samples and data-points per sample.  
   int numSamples = getNumSamples();  
   int numDPPSample = getNumDPPSample();  
   
   //  
   // Determine the sample number which corresponds to this reference number.  
   int sampleNo = -1;  
   int tempRef = -1;  
   for (int n=0; n<numSamples; n++) {  
     tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);  
     if (tempRef == ref) {  
       sampleNo = n;  
       break;  
     }  
   }  
   if (sampleNo == -1) {  
     throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");  
   }  
   
   for (int n=0; n<numDPPSample; n++) {  
     //  
     // Get *each* data-point in the sample in turn.  
     DataArrayView pointView = getDataPoint(sampleNo, n);  
     //  
     // Assign the values in the DataArray to this data-point.  
     pointView.copy(value.getView());  
   }  
 }  
   
 void  
 DataExpanded::getRefValue(int ref,  
                           DataArray& value)  
 {  
   //  
   // Get the number of samples and data-points per sample.  
   int numSamples = getNumSamples();  
   int numDPPSample = getNumDPPSample();  
   
   //  
   // Determine the sample number which corresponds to this reference number  
   int sampleNo = -1;  
   int tempRef = -1;  
   for (int n=0; n<numSamples; n++) {  
     tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);  
     if (tempRef == ref) {  
       sampleNo = n;  
       break;  
     }  
   }  
   if (sampleNo == -1) {  
     throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");  
   }  
   
   //  
   // Get the *first* data-point associated with this sample number.  
   DataArrayView pointView = getDataPoint(sampleNo, 0);  
   
   //  
   // Load the values from this data-point into the DataArray  
   value.getView().copy(pointView);  
 }  
   
313  int  int
314  DataExpanded::archiveData(ofstream& archiveFile,  DataExpanded::archiveData(ofstream& archiveFile,
315                            const DataArrayView::ValueType::size_type noValues) const                            const DataArrayView::ValueType::size_type noValues) const
# Line 409  DataExpanded::extractData(ifstream& arch Line 325  DataExpanded::extractData(ifstream& arch
325  }  }
326    
327  void  void
328    DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
329      //
330      // Get the number of samples and data-points per sample.
331      int numSamples = getNumSamples();
332      int numDataPointsPerSample = getNumDPPSample();
333      int dataPointRank = getPointDataView().getRank();
334      ShapeType dataPointShape = getPointDataView().getShape();
335      if (numSamples*numDataPointsPerSample > 0) {
336         //TODO: global error handling
337         if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
338              throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
339         }
340         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
341               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
342         }
343         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
344         if (dataPointRank==0) {
345             dataPointView()=value;
346         } else if (dataPointRank==1) {
347            for (int i=0; i<dataPointShape[0]; i++) {
348                dataPointView(i)=value;
349            }
350         } else if (dataPointRank==2) {
351            for (int i=0; i<dataPointShape[0]; i++) {
352               for (int j=0; j<dataPointShape[1]; j++) {
353                  dataPointView(i,j)=value;
354               }
355            }
356         } else if (dataPointRank==3) {
357            for (int i=0; i<dataPointShape[0]; i++) {
358               for (int j=0; j<dataPointShape[1]; j++) {
359                  for (int k=0; k<dataPointShape[2]; k++) {
360                     dataPointView(i,j,k)=value;
361                  }
362               }
363            }
364         } else if (dataPointRank==4) {
365             for (int i=0; i<dataPointShape[0]; i++) {
366               for (int j=0; j<dataPointShape[1]; j++) {
367                 for (int k=0; k<dataPointShape[2]; k++) {
368                   for (int l=0; l<dataPointShape[3]; l++) {
369                      dataPointView(i,j,k,l)=value;
370                   }
371                 }
372               }
373             }
374         }
375      }
376    }
377    void
378    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.
381      int numSamples = getNumSamples();
382      int numDataPointsPerSample = getNumDPPSample();
383      int dataPointRank = getPointDataView().getRank();
384      ShapeType dataPointShape = getPointDataView().getShape();
385      //
386      // check rank:
387      if (value.getrank()!=dataPointRank)
388           throw DataException("Rank of numarray does not match Data object rank");
389      if (numSamples*numDataPointsPerSample > 0) {
390         //TODO: global error handling
391         if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
392              throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
393         }
394         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
395               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
396         }
397         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
398         if (dataPointRank==0) {
399             dataPointView()=extract<double>(value[0]);
400         } else if (dataPointRank==1) {
401            for (int i=0; i<dataPointShape[0]; i++) {
402                dataPointView(i)=extract<double>(value[i]);
403            }
404         } else if (dataPointRank==2) {
405            for (int i=0; i<dataPointShape[0]; i++) {
406               for (int j=0; j<dataPointShape[1]; j++) {
407                  dataPointView(i,j)=extract<double>(value[i][j]);
408               }
409            }
410         } else if (dataPointRank==3) {
411            for (int i=0; i<dataPointShape[0]; i++) {
412               for (int j=0; j<dataPointShape[1]; j++) {
413                  for (int k=0; k<dataPointShape[2]; k++) {
414                     dataPointView(i,j,k)=extract<double>(value[i][j][k]);
415                  }
416               }
417            }
418         } else if (dataPointRank==4) {
419             for (int i=0; i<dataPointShape[0]; i++) {
420               for (int j=0; j<dataPointShape[1]; j++) {
421                 for (int k=0; k<dataPointShape[2]; k++) {
422                   for (int l=0; l<dataPointShape[3]; l++) {
423                      dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]);
424                   }
425                 }
426               }
427             }
428         }
429      }
430    }
431    void
432  DataExpanded::copyAll(const boost::python::numeric::array& value) {  DataExpanded::copyAll(const boost::python::numeric::array& value) {
433    //    //
434    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
# Line 463  DataExpanded::copyAll(const boost::pytho Line 483  DataExpanded::copyAll(const boost::pytho
483    }    }
484  }  }
485  void  void
486    DataExpanded::symmetric(DataAbstract* ev)
487    {
488      int sampleNo,dataPointNo;
489      int numSamples = getNumSamples();
490      int numDataPointsPerSample = getNumDPPSample();
491      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
492      if (temp_ev==0) {
493        throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
494      }
495      DataArrayView& thisView=getPointDataView();
496      DataArrayView& evView=ev->getPointDataView();
497      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
498      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
499        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
500             DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
501                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
502        }
503      }
504    }
505    void
506    DataExpanded::nonsymmetric(DataAbstract* ev)
507    {
508      int sampleNo,dataPointNo;
509      int numSamples = getNumSamples();
510      int numDataPointsPerSample = getNumDPPSample();
511      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
512      if (temp_ev==0) {
513        throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
514      }
515      DataArrayView& thisView=getPointDataView();
516      DataArrayView& evView=ev->getPointDataView();
517      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
518      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
519        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
520             DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
521                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
522        }
523      }
524    }
525    void
526    DataExpanded::trace(DataAbstract* ev, int axis_offset)
527    {
528      int sampleNo,dataPointNo;
529      int numSamples = getNumSamples();
530      int numDataPointsPerSample = getNumDPPSample();
531      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
532      if (temp_ev==0) {
533        throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
534      }
535      DataArrayView& thisView=getPointDataView();
536      DataArrayView& evView=ev->getPointDataView();
537      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
538      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
539        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
540             DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
541                                        evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
542        }
543      }
544    }
545    
546    void
547    DataExpanded::transpose(DataAbstract* ev, int axis_offset)
548    {
549      int sampleNo,dataPointNo;
550      int numSamples = getNumSamples();
551      int numDataPointsPerSample = getNumDPPSample();
552      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
553      if (temp_ev==0) {
554        throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
555      }
556      DataArrayView& thisView=getPointDataView();
557      DataArrayView& evView=ev->getPointDataView();
558      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
559      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
560        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
561             DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
562                                        evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
563        }
564      }
565    }
566    
567    void
568    DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
569    {
570      int sampleNo,dataPointNo;
571      int numSamples = getNumSamples();
572      int numDataPointsPerSample = getNumDPPSample();
573      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
574      if (temp_ev==0) {
575        throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
576      }
577      DataArrayView& thisView=getPointDataView();
578      DataArrayView& evView=ev->getPointDataView();
579      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
580      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
581        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
582             DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
583                                        evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
584        }
585      }
586    }
587    void
588  DataExpanded::eigenvalues(DataAbstract* ev)  DataExpanded::eigenvalues(DataAbstract* ev)
589  {  {
590    int sampleNo,dataPointNo;    int sampleNo,dataPointNo;
# Line 510  DataExpanded::eigenvalues_and_eigenvecto Line 632  DataExpanded::eigenvalues_and_eigenvecto
632    }    }
633  }  }
634    
635    void
636    DataExpanded::setToZero(){
637      int numSamples = getNumSamples();
638      int numDataPointsPerSample = getNumDPPSample();
639      DataArrayView& thisView=getPointDataView();
640      DataArrayView::ValueType::size_type n = thisView.noValues();
641      double* p;
642      int  sampleNo,dataPointNo, i;
643      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
644      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
645        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
646            p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
647            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
662    DataExpanded::dump(const std::string fileName) const
663    {
664       #ifdef USE_NETCDF
665       const int ldims=2+DataArrayView::maxRank;
666       const NcDim* ncdims[ldims];
667       NcVar *var, *ids;
668       int rank = getPointDataView().getRank();
669       int type=  getFunctionSpace().getTypeCode();
670       int ndims =0;
671       long dims[ldims];
672       const double* d_ptr=&(m_data[0]);
673       DataArrayView::ShapeType shape = getPointDataView().getShape();
674       int mpi_iam=0, mpi_num=1;
675    
676       // netCDF error handler
677       NcError err(NcError::verbose_nonfatal);
678       // Create the file.
679    #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
686       if (!dataFile.is_valid())
687            throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
688       if (!dataFile.add_att("type_id",2) )
689            throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
690       if (!dataFile.add_att("rank",rank) )
691            throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
692       if (!dataFile.add_att("function_space_type",type))
693            throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
694       ndims=rank+2;
695       if ( rank >0 ) {
696           dims[0]=shape[0];
697           if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
698                throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");
699       }
700       if ( rank >1 ) {
701           dims[1]=shape[1];
702           if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
703                throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");
704       }
705       if ( rank >2 ) {
706           dims[2]=shape[2];
707           if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
708                throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");
709       }
710       if ( rank >3 ) {
711           dims[3]=shape[3];
712           if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
713                throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");
714       }
715       dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
716       if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
717                throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
718       dims[rank+1]=getFunctionSpace().getNumSamples();
719       if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
720                throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
721    
722       if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
723            throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
724       const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
725       if (! (ids->put(ids_p,dims[rank+1])) )
726            throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
727    
728       if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
729            throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
730       if (! (var->put(d_ptr,dims)) )
731            throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
732       #else
733       throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
734       #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.615  
changed lines
  Added in v.1748

  ViewVC Help
Powered by ViewVC 1.1.26