/[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

trunk/esys2/escript/src/Data/DataExpanded.cpp revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC trunk/escript/src/DataExpanded.cpp revision 1748 by ksteube, Wed Sep 3 06:10:39 2008 UTC
# Line 1  Line 1 
 // $Id$  
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
   
 #include "escript/Data/DataException.h"  
 #include "escript/Data/DataExpanded.h"  
 #include "escript/Data/DataConstant.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/DataArrayView.h"  
1    
2  #include <boost/python/extract.hpp>  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                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"
17    #include "DataException.h"
18    #include "DataConstant.h"
19    #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 <iostream>  #include <boost/python/extract.hpp>
28    
29  using namespace std;  using namespace std;
30  using namespace boost::python;  using namespace boost::python;
# Line 35  DataExpanded::DataExpanded(const boost:: Line 38  DataExpanded::DataExpanded(const boost::
38  {  {
39    DataArrayView::ShapeType tempShape;    DataArrayView::ShapeType tempShape;
40    //    //
41    // extract the shape from the python array    // extract the shape of the python numarray
42    for (int i=0; i<value.getrank(); ++i) {    for (int i=0; i<value.getrank(); i++) {
     //cout << extract<int>(value.getshape()[i]) << endl;  
43      tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
44    }    }
45      //
46      // initialise the data array for this object
47    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
48    //    //
49    // copy the given value to every data point    // copy the given value to every data point
# Line 49  DataExpanded::DataExpanded(const boost:: Line 53  DataExpanded::DataExpanded(const boost::
53  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
54    : DataAbstract(other.getFunctionSpace()),    : DataAbstract(other.getFunctionSpace()),
55    m_data(other.m_data)    m_data(other.m_data)
56  {      {
57    //    //
58    // create the view for the data    // create the view for the data
59    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
# Line 58  DataExpanded::DataExpanded(const DataExp Line 62  DataExpanded::DataExpanded(const DataExp
62    
63  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
64    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
65  {      {
66      //
67      // initialise the data array for this object
68    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
69    //    //
70    // DataConstant only has one value    // DataConstant only has one value, copy this to every data point
71    copy(other.getPointDataView());    copy(other.getPointDataView());
72  }  }
73    
74  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
75    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
76  {      {
77      //
78      // initialise the data array for this object
79    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
80      //
81      // for each data point in this object, extract and copy the corresponding data
82      // value from the given DataTagged object
83    int i,j;    int i,j;
84    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
85    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
86  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
87    for (i=0;i<numRows;++i) {    for (i=0;i<numRows;i++) {
88      for (j=0;j<numCols;++j) {      for (j=0;j<numCols;j++) {
89        try {        try {
90          getPointDataView().copy(getPointOffset(i,j), other.getPointDataView(), other.getPointOffset(i,j));          getPointDataView().copy(getPointOffset(i,j),
91                                    other.getPointDataView(),
92                                    other.getPointOffset(i,j));
93        }        }
94        catch (std::exception& e) {        catch (std::exception& e) {
95          cout << e.what() << endl;          cout << e.what() << endl;
# Line 96  DataExpanded::DataExpanded(const DataExp Line 109  DataExpanded::DataExpanded(const DataExp
109    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
110    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(shape,other.getNumSamples(),other.getNumDPPSample());
111    //    //
   DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);  
   //  
112    // copy the data    // copy the data
113      DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
114    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
115    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
116    int i,j;    int i,j;
117  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
118    for (i=0;i<numRows;++i) {    for (i=0;i<numRows;i++) {
119      for (j=0;j<numCols;++j) {      for (j=0;j<numCols;j++) {
120        try {        try {
121          getPointDataView().copySlice(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j),region_loop_range);          getPointDataView().copySlice(getPointOffset(i,j),
122                                         other.getPointDataView(),
123                                         other.getPointOffset(i,j),
124                                         region_loop_range);
125        }        }
126        catch (std::exception& e) {        catch (std::exception& e) {
127          cout << e.what() << endl;          cout << e.what() << endl;
# Line 119  DataExpanded::DataExpanded(const DataArr Line 134  DataExpanded::DataExpanded(const DataArr
134                             const FunctionSpace& what)                             const FunctionSpace& what)
135    : DataAbstract(what)    : DataAbstract(what)
136  {  {
137      //
138      // get the shape of the given data value
139    DataArrayView::ShapeType tempShape=value.getShape();    DataArrayView::ShapeType tempShape=value.getShape();
140      //
141      // initialise this Data object to the shape of the given data value
142    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
143      //
144      // copy the given value to every data point
145    copy(value);    copy(value);
146  }  }
147    
148  DataExpanded::~DataExpanded()  DataExpanded::DataExpanded(const FunctionSpace& what,
149                               const DataArrayView::ShapeType &shape,
150                               const DataArrayView::ValueType &data)
151      : DataAbstract(what)
152  {  {
153      //
154      // create the view of the data
155      initialise(shape,what.getNumSamples(),what.getNumDPPSample());
156      //
157      // copy the data in the correct format
158      m_data.getData()=data;
159  }  }
160    
161  void  DataExpanded::~DataExpanded()
 DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)  
162  {  {
   //  
   // reshape a rank zero data point  
   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());  
   }  
   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) {  
       //  
       // Copy the data into the specified offset  
       // 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.getData()[m_data.index(i,j)]);  
     }  
   }  
   m_data.Swap(newData);  
   DataArrayView temp(m_data.getData(),shape);  
   setPointDataView(temp);  
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) {
177      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
178    }    }
179    //    //
180      // get shape of slice
181    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
182    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
183    //    //
184    // check shape:    // check shape
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    }    }
192    //    //
193    // copy the data    // copy the data from the slice into this object
194    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
195    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
196    int i, j;    int i, j;
197  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
198    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
199      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
200        getPointDataView().copySliceFrom(getPointOffset(i,j),tempDataExp->getPointDataView(),tempDataExp->getPointOffset(i,j),region_loop_range);        getPointDataView().copySliceFrom(getPointOffset(i,j),
201                                           tempDataExp->getPointDataView(),
202                                           tempDataExp->getPointOffset(i,j),
203                                           region_loop_range);
204      }      }
205    }    }
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    // copy a single value to every data point in this object
213    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
214    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
215    int i, j;    int i,j;
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
       // Copy the data into the specified offset  
       // 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      // extract the shape of the numarray
233      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    // first convert the numarray into a DataArrayView format    // check the input shape matches this shape
246    DataArray temp(value);    if (!getPointDataView().checkShape(temp_dataView.getShape())) {
   //  
   // check the input shape matches this shape, this will throw an exception  
   if (!getPointDataView().checkShape(temp.getView().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 entire data structure    // now copy over the data
253    copy(temp.getView());    copy(temp_dataView);
254  }  }
255    
256  void  void
# Line 249  DataExpanded::initialise(const DataArray Line 259  DataExpanded::initialise(const DataArray
259                           int noDataPointsPerSample)                           int noDataPointsPerSample)
260  {  {
261    //    //
262    // resize data to the required size    // resize data array to the required size
263    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
264    //    //
265    // create a point data viewer of the data    // create the data view of the data array
266    DataArrayView temp(m_data.getData(),shape);    DataArrayView temp(m_data.getData(),shape);
267    setPointDataView(temp);    setPointDataView(temp);
268  }  }
# Line 261  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 289  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    
307  DataArrayView::ValueType::size_type  DataArrayView::ValueType::size_type
308  DataExpanded::getLength() const  DataExpanded::getLength() const
309  {  {
310    return m_data.getData().size();    return m_data.size();
311  }  }
312    
313    int
314    DataExpanded::archiveData(ofstream& archiveFile,
315                              const DataArrayView::ValueType::size_type noValues) const
316    {
317      return(m_data.archiveData(archiveFile, noValues));
318    }
319    
320    int
321    DataExpanded::extractData(ifstream& archiveFile,
322                              const DataArrayView::ValueType::size_type noValues)
323    {
324      return(m_data.extractData(archiveFile, noValues));
325    }
326    
327    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) {
433      //
434      // Get the number of samples and data-points per sample.
435      int numSamples = getNumSamples();
436      int numDataPointsPerSample = getNumDPPSample();
437      int dataPointRank = getPointDataView().getRank();
438      ShapeType dataPointShape = getPointDataView().getShape();
439      //
440      // check rank:
441      if (value.getrank()!=dataPointRank+1)
442           throw DataException("Rank of numarray does not match Data object rank");
443      if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
444           throw DataException("leading dimension of numarray is too small");
445      //
446      int dataPoint = 0;
447      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
448        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
449          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
450          if (dataPointRank==0) {
451             dataPointView()=extract<double>(value[dataPoint]);
452          } else if (dataPointRank==1) {
453             for (int i=0; i<dataPointShape[0]; i++) {
454                dataPointView(i)=extract<double>(value[dataPoint][i]);
455             }
456          } else if (dataPointRank==2) {
457             for (int i=0; i<dataPointShape[0]; i++) {
458               for (int j=0; j<dataPointShape[1]; j++) {
459                 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
460               }
461             }
462           } else if (dataPointRank==3) {
463             for (int i=0; i<dataPointShape[0]; i++) {
464               for (int j=0; j<dataPointShape[1]; j++) {
465                 for (int k=0; k<dataPointShape[2]; k++) {
466                     dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
467                 }
468               }
469             }
470           } else if (dataPointRank==4) {
471             for (int i=0; i<dataPointShape[0]; i++) {
472               for (int j=0; j<dataPointShape[1]; j++) {
473                 for (int k=0; k<dataPointShape[2]; k++) {
474                   for (int l=0; l<dataPointShape[3]; l++) {
475                     dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
476                   }
477                 }
478               }
479             }
480          }
481          dataPoint++;
482        }
483      }
484    }
485    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)
589    {
590      int sampleNo,dataPointNo;
591      int numSamples = getNumSamples();
592      int numDataPointsPerSample = getNumDPPSample();
593      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
594      if (temp_ev==0) {
595        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
596      }
597      DataArrayView& thisView=getPointDataView();
598      DataArrayView& evView=ev->getPointDataView();
599      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
600      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
601        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
602             DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
603                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
604        }
605      }
606    }
607    void
608    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
609    {
610      int numSamples = getNumSamples();
611      int numDataPointsPerSample = getNumDPPSample();
612      int sampleNo,dataPointNo;
613      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
614      if (temp_ev==0) {
615        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
616      }
617      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
618      if (temp_V==0) {
619        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
620      }
621      DataArrayView& thisView=getPointDataView();
622      DataArrayView& evView=ev->getPointDataView();
623      DataArrayView& VView=V->getPointDataView();
624      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
625      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
626        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
627             DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
628                                                         evView,ev->getPointOffset(sampleNo,dataPointNo),
629                                                         VView,V->getPointOffset(sampleNo,dataPointNo),
630                                                         tol);
631        }
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.108  
changed lines
  Added in v.1748

  ViewVC Help
Powered by ViewVC 1.1.26