/[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 117 by jgs, Fri Apr 1 05:48:57 2005 UTC trunk/escript/src/DataExpanded.cpp revision 1962 by gross, Tue Nov 4 07:17:31 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>  /*******************************************************
3    *
4    * Copyright (c) 2003-2008 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15    #include "Data.h"
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    #include "DataMaths.h"
29    
30  using namespace std;  using namespace std;
31  using namespace boost::python;  using namespace boost::python;
32  using namespace boost;  using namespace boost;
33    using namespace escript::DataTypes;
34    
35  namespace escript {  namespace escript {
36    
37  DataExpanded::DataExpanded(const boost::python::numeric::array& value,  DataExpanded::DataExpanded(const boost::python::numeric::array& value,
38                             const FunctionSpace& what)                             const FunctionSpace& what)
39    : DataAbstract(what)    : DataAbstract(what,DataTypes::shapeFromNumArray(value))
40  {  {
   DataArrayView::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(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(what.getNumSamples(),what.getNumDPPSample());
44    //    //
45    // copy the given value to every data point    // copy the given value to every data point
46    copy(value);    copy(value);
47  }  }
48    
49  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
50    : DataAbstract(other.getFunctionSpace()),    : 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);  
53  }  }
54    
55  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
56    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(), other.getShape())
57  {  {
58    //    //
59    // initialise the data array for this object    // initialise the data array for this object
60    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
61    //    //
62    // DataConstant only has one value, copy this to every data point    // DataConstant only has one value, copy this to every data point
63    copy(other.getPointDataView());    copy(other);
64  }  }
65    
66  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
67    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(), other.getShape())
68  {  {
69    //    //
70    // initialise the data for this object    // initialise the data array for this object
71    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
72    //    //
73    // for each data point in this object, extract and copy the corresponding data    // for each data point in this object, extract and copy the corresponding data
74    // value from the given DataTagged object    // value from the given DataTagged object
75    int i,j;    int i,j;
76    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
77    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
78  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
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 {
82          getPointDataView().copy(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j));             DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
83                                    other.getVector(),
84                                    other.getPointOffset(i,j));
85        }        }
86        catch (std::exception& e) {        catch (std::exception& e) {
87          cout << e.what() << endl;          cout << e.what() << endl;
# Line 94  DataExpanded::DataExpanded(const DataTag Line 91  DataExpanded::DataExpanded(const DataTag
91  }  }
92    
93  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
94                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
95    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
96  {  {
97    //    //
98    // get the shape of the slice    // get the shape of the slice
99    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
100    //    //
101    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
102    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
103    //    //
104    // copy the data    // copy the data
105    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
106    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
107    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
108    int i,j;    int i,j;
109  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
110    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
111      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
112        try {        try {
113          getPointDataView().copySlice(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j),region_loop_range);  //         getPointDataView().copySlice(getPointOffset(i,j),
114    //                                      other.getPointDataView(),
115    //                                      other.getPointOffset(i,j),
116    //                                      region_loop_range);
117            DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
118                                         other.getVector(),
119                         other.getShape(),
120                                         other.getPointOffset(i,j),
121                                         region_loop_range);
122        }        }
123        catch (std::exception& e) {        catch (std::exception& e) {
124          cout << e.what() << endl;          cout << e.what() << endl;
# Line 122  DataExpanded::DataExpanded(const DataExp Line 127  DataExpanded::DataExpanded(const DataExp
127    }    }
128  }  }
129    
130  DataExpanded::DataExpanded(const DataArrayView& value,  // DataExpanded::DataExpanded(const DataArrayView& value,
131                             const FunctionSpace& what)  //                            const FunctionSpace& what)
132    : DataAbstract(what)  //   : DataAbstract(what)
133  {  // {
134    //  //   //
135    // get the shape of the given data value  //   // get the shape of the given data value
136    DataArrayView::ShapeType tempShape=value.getShape();  //   DataTypes::ShapeType tempShape=value.getShape();
137    //  //   //
138    // initialise this Data object to the shape of the given data value  //   // initialise this Data object to the shape of the given data value
139    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());  //   initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
140    //  //   //
141    // copy the given value to every data point  //   // copy the given value to every data point
142    copy(value);  //   copy(value);
143    // }
144    
145    DataExpanded::DataExpanded(const FunctionSpace& what,
146                               const DataTypes::ShapeType &shape,
147                               const DataTypes::ValueType &data)
148      : DataAbstract(what,shape)
149    {
150      EsysAssert(data.size()%getNoValues()==0,
151                     "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
152    
153      if (data.size()==getNoValues())
154      {
155         ValueType& vec=m_data.getData();
156         //
157         // create the view of the data
158         initialise(what.getNumSamples(),what.getNumDPPSample());
159         // now we copy this value to all elements
160         for (int i=0;i<getLength();)
161         {
162        for (int j=0;j<getNoValues();++j,++i)
163        {
164            vec[i]=data[j];
165        }
166         }
167      }
168      else
169      {
170         //
171         // copy the data in the correct format
172         m_data.getData()=data;
173      }
174    
175    
176  }  }
177    
178  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
179  {  {
180  }  }
181    
182  void  DataAbstract*
183  DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)  DataExpanded::deepCopy()
184  {  {
185    if (getPointDataView().getRank()!=0) {    return new DataExpanded(*this);
     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);  
186  }  }
187    
188    
189  DataAbstract*  DataAbstract*
190  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::getSlice(const DataTypes::RegionType& region) const
191  {  {
192    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
193  }  }
194    
195  void  void
196  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
197                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
198  {  {
199    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
200    if (tempDataExp==0) {    if (tempDataExp==0) {
# Line 193  DataExpanded::setSlice(const DataAbstrac Line 202  DataExpanded::setSlice(const DataAbstrac
202    }    }
203    //    //
204    // get shape of slice    // get shape of slice
205    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
206    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
207    //    //
208    // check shape    // check shape
209    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
210      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
211    }    }
212    if (tempDataExp->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
213      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
214          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
215    }    }
216    //    //
217    // copy the data from the slice into this object    // copy the data from the slice into this object
218    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
219    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
220    int i, j;    int i, j;
221  #pragma omp parallel for private(i,j) schedule(static)    ValueType& vec=getVector();
222      const ShapeType& mshape=getShape();
223      const ValueType& tVec=tempDataExp->getVector();
224      const ShapeType& tShape=tempDataExp->getShape();
225      #pragma omp parallel for private(i,j) schedule(static)
226    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
227      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
228        getPointDataView().copySliceFrom(getPointOffset(i,j),tempDataExp->getPointDataView(),tempDataExp->getPointOffset(i,j),region_loop_range);  /*      getPointDataView().copySliceFrom(getPointOffset(i,j),
229                                           tempDataExp->getPointDataView(),
230                                           tempDataExp->getPointOffset(i,j),
231                                           region_loop_range);*/
232            DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
233                                           tVec,
234                           tShape,
235                                           tempDataExp->getPointOffset(i,j),
236                                           region_loop_range);
237    
238      }      }
239    }    }
240  }  }
241    
242  void  void
243  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataConstant& value)
244  {  {
245      EsysAssert((checkShape(getShape(), value.getShape())),
246                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
247    
248    //    //
249    // copy a single value to every data point in this object    // copy a single value to every data point in this object
250    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
251    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
252    int i,j;    int i,j;
253  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
254    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
255      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
256        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
257        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
258        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
259        // if using DOASSERT.        // if using DOASSERT.
260        getPointDataView().copy(m_data.index(i,j),value);        //getPointDataView().copy(getPointOffset(i,j),value);
261          DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
262      }      }
263    }    }
264  }  }
265    
266    
267    // void
268    // DataExpanded::copy(const DataArrayView& value)
269    // {
270    //   //
271    //   // copy a single value to every data point in this object
272    //   int nRows=m_data.getNumRows();
273    //   int nCols=m_data.getNumCols();
274    //   int i,j;
275    //   #pragma omp parallel for private(i,j) schedule(static)
276    //   for (i=0;i<nRows;i++) {
277    //     for (j=0;j<nCols;j++) {
278    //       // NOTE: An exception may be thown from this call if
279    //       // DOASSERT is on which of course will play
280    //       // havoc with the omp threads. Run single threaded
281    //       // if using DOASSERT.
282    //       getPointDataView().copy(getPointOffset(i,j),value);
283    //     }
284    //   }
285    // }
286    
287  void  void
288  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
289  {  {
290    //  
291    // first convert the numarray into a DataArray object    // extract the shape of the numarray
292    DataArray temp(value);    DataTypes::ShapeType tempShape;
293      for (int i=0; i < value.getrank(); i++) {
294        tempShape.push_back(extract<int>(value.getshape()[i]));
295      }
296    
297      // get the space for the data vector
298    //   int len = DataTypes::noValues(tempShape);
299    //   DataVector temp_data(len, 0.0, len);
300    //   DataArrayView temp_dataView(temp_data, tempShape);
301    //   temp_dataView.copy(value);
302    
303    //    //
304    // check the input shape matches this shape    // check the input shape matches this shape
305    if (!getPointDataView().checkShape(temp.getView().getShape())) {    if (!DataTypes::checkShape(getShape(),tempShape)) {
306      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
307                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
308                          temp.getView().getShape()));                          tempShape,getShape()));
309    }    }
310    //    //
311    // now copy over the data    // now copy over the data
312    copy(temp.getView());    //copy(temp_dataView);
313      getVector().copyFromNumArray(value);
314  }  }
315    
316    
317  void  void
318  DataExpanded::initialise(const DataArrayView::ShapeType& shape,  DataExpanded::initialise(int noSamples,
                          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,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
   //  
   // create the data view of the data array  
   DataArrayView temp(m_data.getData(),shape);  
   setPointDataView(temp);  
328  }  }
329    
330  string  string
331  DataExpanded::toString() const  DataExpanded::toString() const
332  {  {
333    stringstream temp;    stringstream temp;
334    //    FunctionSpace fs=getFunctionSpace();
335    // create a temporary view as the offset will be changed  
336    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());    int offset=0;
337    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
338      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
339        tempView.setOffset(m_data.index(i,j));        offset=getPointOffset(i,j);
340        stringstream suffix;        stringstream suffix;
341        suffix << "(" << i << "," << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
342        temp << tempView.toString(suffix.str());        temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
343        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
344          temp << endl;          temp << endl;
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    
356  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
357  DataExpanded::getPointOffset(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
358                               int dataPointNo) const                               int dataPointNo) const
359  {  {
360    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
361  }  }
362    
363  DataArrayView  DataTypes::ValueType::size_type
 DataExpanded::getDataPoint(int sampleNo,  
                            int dataPointNo)  
 {  
   DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));  
   return temp;  
 }  
   
 DataArrayView::ValueType::size_type  
364  DataExpanded::getLength() const  DataExpanded::getLength() const
365  {  {
366    return m_data.size();    return m_data.size();
367  }  }
368    
369    
370    
371  void  void
372  DataExpanded::setRefValue(int ref,  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
                           const DataArray& value)  
 {  
373    //    //
374    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
375    int numSamples = getNumSamples();    int numSamples = getNumSamples();
376    int numDPPSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
377      int dataPointRank = getRank();
378      ShapeType dataPointShape = getShape();
379      if (numSamples*numDataPointsPerSample > 0) {
380         //TODO: global error handling
381         if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
382              throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
383         }
384         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
385               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
386         }
387         ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
388         ValueType& vec=getVector();
389         if (dataPointRank==0) {
390             vec[offset]=value;
391         } else if (dataPointRank==1) {
392            for (int i=0; i<dataPointShape[0]; i++) {
393                vec[offset+i]=value;
394            }
395         } else if (dataPointRank==2) {
396            for (int i=0; i<dataPointShape[0]; i++) {
397               for (int j=0; j<dataPointShape[1]; j++) {
398                  vec[offset+getRelIndex(dataPointShape,i,j)]=value;
399               }
400            }
401         } else if (dataPointRank==3) {
402            for (int i=0; i<dataPointShape[0]; i++) {
403               for (int j=0; j<dataPointShape[1]; j++) {
404                  for (int k=0; k<dataPointShape[2]; k++) {
405                     vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
406                  }
407               }
408            }
409         } else if (dataPointRank==4) {
410             for (int i=0; i<dataPointShape[0]; i++) {
411               for (int j=0; j<dataPointShape[1]; j++) {
412                 for (int k=0; k<dataPointShape[2]; k++) {
413                   for (int l=0; l<dataPointShape[3]; l++) {
414                      vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
415                   }
416                 }
417               }
418             }
419         }
420      }
421    }
422    void
423    DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
424      //
425      // Get the number of samples and data-points per sample.
426      int numSamples = getNumSamples();
427      int numDataPointsPerSample = getNumDPPSample();
428      int dataPointRank = getRank();
429      const ShapeType& shape = getShape();
430      //
431      // check rank:
432      if (value.getrank()!=dataPointRank)
433           throw DataException("Rank of numarray does not match Data object rank");
434      if (numSamples*numDataPointsPerSample > 0) {
435         //TODO: global error handling
436         if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
437              throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
438         }
439         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
440               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
441         }
442         ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
443         ValueType& vec=getVector();
444         if (dataPointRank==0) {
445             vec[offset]=extract<double>(value[0]);
446         } else if (dataPointRank==1) {
447            for (int i=0; i<shape[0]; i++) {
448                vec[offset+i]=extract<double>(value[i]);
449            }
450         } else if (dataPointRank==2) {
451            for (int i=0; i<shape[0]; i++) {
452               for (int j=0; j<shape[1]; j++) {
453                  vec[offset+getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
454               }
455            }
456         } else if (dataPointRank==3) {
457            for (int i=0; i<shape[0]; i++) {
458               for (int j=0; j<shape[1]; j++) {
459                  for (int k=0; k<shape[2]; k++) {
460                     vec[offset+getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
461                  }
462               }
463            }
464         } else if (dataPointRank==4) {
465             for (int i=0; i<shape[0]; i++) {
466               for (int j=0; j<shape[1]; j++) {
467                 for (int k=0; k<shape[2]; k++) {
468                   for (int l=0; l<shape[3]; l++) {
469                      vec[offset+getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
470                   }
471                 }
472               }
473             }
474         }
475      }
476    }
477    void
478    DataExpanded::copyAll(const boost::python::numeric::array& value) {
479    //    //
480    // Determine the sample number which corresponds to this reference number.    // Get the number of samples and data-points per sample.
481    int sampleNo = -1;    int numSamples = getNumSamples();
482    int tempRef = -1;    int numDataPointsPerSample = getNumDPPSample();
483    for (int n=0; n<numSamples; n++) {    int dataPointRank = getRank();
484      tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);    const ShapeType& dataPointShape = getShape();
485      if (tempRef == ref) {    //
486        sampleNo = n;    // check rank:
487        break;    if (value.getrank()!=dataPointRank+1)
488           throw DataException("Rank of numarray does not match Data object rank");
489      if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
490           throw DataException("leading dimension of numarray is too small");
491      //
492      ValueType& vec=getVector();
493      int dataPoint = 0;
494      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
495        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
496          ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
497          if (dataPointRank==0) {
498             vec[offset]=extract<double>(value[dataPoint]);
499          } else if (dataPointRank==1) {
500             for (int i=0; i<dataPointShape[0]; i++) {
501                vec[offset+i]=extract<double>(value[dataPoint][i]);
502             }
503          } else if (dataPointRank==2) {
504             for (int i=0; i<dataPointShape[0]; i++) {
505               for (int j=0; j<dataPointShape[1]; j++) {
506             vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
507               }
508             }
509           } else if (dataPointRank==3) {
510             for (int i=0; i<dataPointShape[0]; i++) {
511               for (int j=0; j<dataPointShape[1]; j++) {
512                 for (int k=0; k<dataPointShape[2]; k++) {
513             vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
514                 }
515               }
516             }
517           } else if (dataPointRank==4) {
518             for (int i=0; i<dataPointShape[0]; i++) {
519               for (int j=0; j<dataPointShape[1]; j++) {
520                 for (int k=0; k<dataPointShape[2]; k++) {
521                   for (int l=0; l<dataPointShape[3]; l++) {
522                     vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
523                   }
524                 }
525               }
526             }
527          }
528          dataPoint++;
529        }
530      }
531    }
532    void
533    DataExpanded::symmetric(DataAbstract* ev)
534    {
535      int sampleNo,dataPointNo;
536      int numSamples = getNumSamples();
537      int numDataPointsPerSample = getNumDPPSample();
538      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
539      if (temp_ev==0) {
540        throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
541      }
542      ValueType& vec=getVector();
543      const ShapeType& shape=getShape();
544      ValueType& evVec=temp_ev->getVector();
545      const ShapeType& evShape=temp_ev->getShape();
546      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
547      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
548        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
549             DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
550                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
551        }
552      }
553    }
554    void
555    DataExpanded::nonsymmetric(DataAbstract* ev)
556    {
557      int sampleNo,dataPointNo;
558      int numSamples = getNumSamples();
559      int numDataPointsPerSample = getNumDPPSample();
560      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
561      if (temp_ev==0) {
562        throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
563      }
564      ValueType& vec=getVector();
565      const ShapeType& shape=getShape();
566      ValueType& evVec=temp_ev->getVector();
567      const ShapeType& evShape=temp_ev->getShape();
568      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
569      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
570        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
571             DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
572                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
573      }      }
574    }    }
575    if (sampleNo == -1) {  }
576      throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");  void
577    DataExpanded::trace(DataAbstract* ev, int axis_offset)
578    {
579      int sampleNo,dataPointNo;
580      int numSamples = getNumSamples();
581      int numDataPointsPerSample = getNumDPPSample();
582      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
583      if (temp_ev==0) {
584        throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
585      }
586      ValueType& vec=getVector();
587      const ShapeType& shape=getShape();
588      ValueType& evVec=temp_ev->getVector();
589      const ShapeType& evShape=temp_ev->getShape();
590      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
591      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
592        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
593             DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
594                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
595        }
596    }    }
597    }
598    
599    for (int n=0; n<numDPPSample; n++) {  void
600      //  DataExpanded::transpose(DataAbstract* ev, int axis_offset)
601      // Get each data-point in the sample in turn.  {
602      DataArrayView pointView = getDataPoint(sampleNo, n);    int sampleNo,dataPointNo;
603      //    int numSamples = getNumSamples();
604      // Assign the values in the DataArray to this data-point.    int numDataPointsPerSample = getNumDPPSample();
605      pointView.copy(value.getView());    DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
606      if (temp_ev==0) {
607        throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
608      }
609      ValueType& vec=getVector();
610      const ShapeType& shape=getShape();
611      ValueType& evVec=temp_ev->getVector();
612      const ShapeType& evShape=temp_ev->getShape();
613      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
614      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
615        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
616             DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
617                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
618        }
619    }    }
620  }  }
621    
622  void  void
623  DataExpanded::getRefValue(int ref,  DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
624                            DataArray& value)  {
625      int sampleNo,dataPointNo;
626      int numSamples = getNumSamples();
627      int numDataPointsPerSample = getNumDPPSample();
628      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
629      if (temp_ev==0) {
630        throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
631      }
632      ValueType& vec=getVector();
633      const ShapeType& shape=getShape();
634      ValueType& evVec=temp_ev->getVector();
635      const ShapeType& evShape=temp_ev->getShape();
636      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
637      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
638        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
639             DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
640                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
641        }
642      }
643    }
644    void
645    DataExpanded::eigenvalues(DataAbstract* ev)
646    {
647      int sampleNo,dataPointNo;
648      int numSamples = getNumSamples();
649      int numDataPointsPerSample = getNumDPPSample();
650      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
651      if (temp_ev==0) {
652        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
653      }
654      ValueType& vec=getVector();
655      const ShapeType& shape=getShape();
656      ValueType& evVec=temp_ev->getVector();
657      const ShapeType& evShape=temp_ev->getShape();
658      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
659      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
660        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
661             DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
662                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
663        }
664      }
665    }
666    void
667    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
668  {  {
   //  
   // Get the number of samples and data-points per sample.  
669    int numSamples = getNumSamples();    int numSamples = getNumSamples();
670    int numDPPSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
671      int sampleNo,dataPointNo;
672      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
673      if (temp_ev==0) {
674        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
675      }
676      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
677      if (temp_V==0) {
678        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
679      }
680      ValueType& vec=getVector();
681      const ShapeType& shape=getShape();
682      ValueType& evVec=temp_ev->getVector();
683      const ShapeType& evShape=temp_ev->getShape();
684      ValueType& VVec=temp_V->getVector();
685      const ShapeType& VShape=temp_V->getShape();
686      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
687      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
688        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
689             DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
690                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
691                                        VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
692        }
693      }
694    }
695    
696    //  void
697    // Determine the sample number which corresponds to this reference number  DataExpanded::setToZero(){
698    int sampleNo = -1;  // TODO: Surely there is a more efficient way to do this????
699    int tempRef = -1;  // Why is there no memset here? Parallel issues?
700    for (int n=0; n<numSamples; n++) {    int numSamples = getNumSamples();
701      tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);    int numDataPointsPerSample = getNumDPPSample();
702      if (tempRef == ref) {    DataTypes::ValueType::size_type n = getNoValues();
703        sampleNo = n;    double* p;
704        break;    int  sampleNo,dataPointNo, i;
705      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
706      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
707        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
708            p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
709            for (i=0; i<n ;++i) p[i]=0.;
710      }      }
711    }    }
712    if (sampleNo == -1) {  }
713      throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");  
714    /* Append MPI rank to file name if multiple MPI processes */
715    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
716      /* Make plenty of room for the mpi_rank number and terminating '\0' */
717      char *newFileName = (char *)malloc(strlen(fileName)+20);
718      strncpy(newFileName, fileName, strlen(fileName)+1);
719      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
720      return(newFileName);
721    }
722    
723    void
724    DataExpanded::dump(const std::string fileName) const
725    {
726       #ifdef USE_NETCDF
727       const int ldims=2+DataTypes::maxRank;
728       const NcDim* ncdims[ldims];
729       NcVar *var, *ids;
730       int rank = getRank();
731       int type=  getFunctionSpace().getTypeCode();
732       int ndims =0;
733       long dims[ldims];
734       const double* d_ptr=&(m_data[0]);
735       const DataTypes::ShapeType& shape = getShape();
736       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
737       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
738    #ifdef PASO_MPI
739       MPI_Status status;
740    #endif
741    
742    #ifdef PASO_MPI
743       /* Serialize NetCDF I/O */
744       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
745    #endif
746    
747       // netCDF error handler
748       NcError err(NcError::verbose_nonfatal);
749       // Create the file.
750       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
751       NcFile dataFile(newFileName, NcFile::Replace);
752       // check if writing was successful
753       if (!dataFile.is_valid())
754            throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
755       if (!dataFile.add_att("type_id",2) )
756            throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
757       if (!dataFile.add_att("rank",rank) )
758            throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
759       if (!dataFile.add_att("function_space_type",type))
760            throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
761       ndims=rank+2;
762       if ( rank >0 ) {
763           dims[0]=shape[0];
764           if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
765                throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
766       }
767       if ( rank >1 ) {
768           dims[1]=shape[1];
769           if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
770                throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
771       }
772       if ( rank >2 ) {
773           dims[2]=shape[2];
774           if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
775                throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
776       }
777       if ( rank >3 ) {
778           dims[3]=shape[3];
779           if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
780                throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
781       }
782       dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
783       if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
784                throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
785       dims[rank+1]=getFunctionSpace().getNumSamples();
786       if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
787                throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
788    
789       if (getFunctionSpace().getNumSamples()>0)
790       {
791    
792         if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
793            throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
794         const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
795         if (! (ids->put(ids_p,dims[rank+1])) )
796            throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
797         if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
798            throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
799         if (! (var->put(d_ptr,dims)) )
800            throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
801       }
802    #ifdef PASO_MPI
803       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
804    #endif
805       #else
806       throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
807       #endif
808    }
809    
810    void  
811    DataExpanded::setTaggedValue(int tagKey,
812               const DataTypes::ShapeType& pointshape,
813                   const DataTypes::ValueType& value,
814               int dataOffset)
815    {
816      int numSamples = getNumSamples();
817      int numDataPointsPerSample = getNumDPPSample();
818      int sampleNo,dataPointNo, i;
819      DataTypes::ValueType::size_type n = getNoValues();
820      double* p;
821      const double* in=&value[0+dataOffset];
822      
823      if (value.size() != n) {
824        throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
825      }
826    
827      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
828      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
829        if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
830            for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
831                p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
832                for (i=0; i<n ;++i) p[i]=in[i];
833            }
834        }
835    }    }
836    }
837    
   //  
   // Get the first data-point associated with this sample number.  
   DataArrayView pointView = getDataPoint(sampleNo, 0);  
838    
839    //  void
840    // Load the values from this data-point into the DataArray  DataExpanded::reorderByReferenceIDs(int *reference_ids)
841    value.getView().copy(pointView);  {
842      int numSamples = getNumSamples();
843      DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
844      int sampleNo, sampleNo2,i;
845      double* p,*p2;
846      register double rtmp;
847      FunctionSpace fs=getFunctionSpace();
848    
849      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
850         const int id_in=reference_ids[sampleNo];
851         const int id=fs.getReferenceIDOfSample(sampleNo);
852         if (id!=id_in) {
853             bool matched=false;
854             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
855                  if (id == reference_ids[sampleNo2]) {
856                     p=&(m_data[getPointOffset(sampleNo,0)]);
857                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
858                     for (i=0; i<n ;i++) {
859                             rtmp=p[i];
860                             p[i]=p2[i];
861                             p2[i]=rtmp;
862                     }
863                     reference_ids[sampleNo]=id;
864                     reference_ids[sampleNo2]=id_in;
865                     matched=true;
866                     break;
867                  }
868             }
869             if (! matched) {
870                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
871             }
872         }
873       }
874    }
875    
876    DataTypes::ValueType&
877    DataExpanded::getVector()
878    {
879        return m_data.getData();
880    }
881    
882    const DataTypes::ValueType&
883    DataExpanded::getVector() const
884    {
885        return m_data.getData();
886  }  }
887    
888  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.117  
changed lines
  Added in v.1962

  ViewVC Help
Powered by ViewVC 1.1.26