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

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

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

revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC revision 119 by jgs, Tue Apr 12 04:45:05 2005 UTC
# Line 18  Line 18 
18  #include "escript/Data/Data.h"  #include "escript/Data/Data.h"
19    
20  #include <iostream>  #include <iostream>
21    #include <fstream>
22  #include <algorithm>  #include <algorithm>
23  #include <vector>  #include <vector>
24  #include <exception>  #include <exception>
# Line 29  Line 30 
30  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
31    
32  #include "escript/Data/DataException.h"  #include "escript/Data/DataException.h"
   
33  #include "escript/Data/DataExpanded.h"  #include "escript/Data/DataExpanded.h"
34  #include "escript/Data/DataConstant.h"  #include "escript/Data/DataConstant.h"
35  #include "escript/Data/DataTagged.h"  #include "escript/Data/DataTagged.h"
# Line 444  Data::getDataPointShape() const Line 444  Data::getDataPointShape() const
444  }  }
445    
446  boost::python::numeric::array  boost::python::numeric::array
447    Data::convertToNumArray()
448    {
449      //
450      // determine the current number of data points
451      int numSamples = getNumSamples();
452      int numDataPointsPerSample = getNumDataPointsPerSample();
453      int numDataPoints = numSamples * numDataPointsPerSample;
454    
455      //
456      // determine the rank and shape of each data point
457      int dataPointRank = getDataPointRank();
458      DataArrayView::ShapeType dataPointShape = getDataPointShape();
459    
460      //
461      // create the numeric array to be returned
462      boost::python::numeric::array numArray(0.0);
463    
464      //
465      // the rank of the returned numeric array will be the rank of
466      // the data points, plus one. Where the rank of the array is n,
467      // the last n-1 dimensions will be equal to the shape of the
468      // data points, whilst the first dimension will be equal to the
469      // total number of data points. Thus the array will consist of
470      // a serial vector of the data points.
471      int arrayRank = dataPointRank + 1;
472      DataArrayView::ShapeType arrayShape;
473      arrayShape.push_back(numDataPoints);
474      for (int d=0; d<dataPointRank; d++) {
475         arrayShape.push_back(dataPointShape[d]);
476      }
477    
478      //
479      // resize the numeric array to the shape just calculated
480      if (arrayRank==1) {
481        numArray.resize(arrayShape[0]);
482      }
483      if (arrayRank==2) {
484        numArray.resize(arrayShape[0],arrayShape[1]);
485      }
486      if (arrayRank==3) {
487        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
488      }
489      if (arrayRank==4) {
490        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
491      }
492      if (arrayRank==5) {
493        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
494      }
495    
496      //
497      // loop through each data point in turn, loading the values for that data point
498      // into the numeric array.
499      int dataPoint = 0;
500      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
501        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
502    
503          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
504    
505          if (dataPointRank==0) {
506            numArray[dataPoint]=dataPointView();
507          }
508    
509          if (dataPointRank==1) {
510            for (int i=0; i<dataPointShape[0]; i++) {
511              numArray[dataPoint][i]=dataPointView(i);
512            }
513          }
514    
515          if (dataPointRank==2) {
516            for (int i=0; i<dataPointShape[0]; i++) {
517              for (int j=0; j<dataPointShape[1]; j++) {
518                numArray[dataPoint][i][j] = dataPointView(i,j);
519              }
520            }
521          }
522    
523          if (dataPointRank==3) {
524            for (int i=0; i<dataPointShape[0]; i++) {
525              for (int j=0; j<dataPointShape[1]; j++) {
526                for (int k=0; k<dataPointShape[2]; k++) {
527                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
528                }
529              }
530            }
531          }
532    
533          if (dataPointRank==4) {
534            for (int i=0; i<dataPointShape[0]; i++) {
535              for (int j=0; j<dataPointShape[1]; j++) {
536                for (int k=0; k<dataPointShape[2]; k++) {
537                  for (int l=0; l<dataPointShape[3]; l++) {
538                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
539                  }
540                }
541              }
542            }
543          }
544    
545          dataPoint++;
546    
547        }
548      }
549    
550      //
551      // return the loaded array
552      return numArray;
553    }
554    
555    boost::python::numeric::array
556  Data::integrate() const  Data::integrate() const
557  {  {
558    int index;    int index;
# Line 585  Data::Lsup() const Line 694  Data::Lsup() const
694  }  }
695    
696  double  double
697    Data::Linf() const
698    {
699      //
700      // set the initial absolute minimum value to max double
701      return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));
702    }
703    
704    double
705  Data::sup() const  Data::sup() const
706  {  {
707    //    //
# Line 603  Data::inf() const Line 720  Data::inf() const
720  Data  Data
721  Data::maxval() const  Data::maxval() const
722  {  {
723      //
724      // set the initial maximum value to min possible double
725    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
726  }  }
727    
728  Data  Data
729  Data::minval() const  Data::minval() const
730  {  {
731      //
732      // set the initial minimum value to max possible double
733    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
734  }  }
735    
# Line 639  Data::saveDX(std::string fileName) const Line 760  Data::saveDX(std::string fileName) const
760    return;    return;
761  }  }
762    
763    void
764    Data::saveVTK(std::string fileName) const
765    {
766      getDomain().saveVTK(fileName,*this);
767      return;
768    }
769    
770  Data&  Data&
771  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
772  {  {
# Line 1030  Data::setTaggedValue(int tagKey, Line 1158  Data::setTaggedValue(int tagKey,
1158    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1159  }  }
1160    
1161    void
1162    Data::setRefValue(int ref,
1163                      const boost::python::numeric::array& value)
1164    {
1165      //
1166      // Construct DataArray from boost::python::object input value
1167      DataArray valueDataArray(value);
1168    
1169      //
1170      // Call DataAbstract::setRefValue
1171      m_data->setRefValue(ref,valueDataArray);
1172    }
1173    
1174    void
1175    Data::getRefValue(int ref,
1176                      boost::python::numeric::array& value)
1177    {
1178      //
1179      // Construct DataArray for boost::python::object return value
1180      DataArray valueDataArray(value);
1181    
1182      //
1183      // Load DataArray with values from data-points specified by ref
1184      m_data->getRefValue(ref,valueDataArray);
1185    
1186      //
1187      // Load values from valueDataArray into return numarray
1188    
1189      // extract the shape of the numarray
1190      int rank = value.getrank();
1191      DataArrayView::ShapeType shape;
1192      for (int i=0; i < rank; i++) {
1193        shape.push_back(extract<int>(value.getshape()[i]));
1194      }
1195    
1196      // and load the numarray with the data from the DataArray
1197      DataArrayView valueView = valueDataArray.getView();
1198    
1199      if (rank==0) {
1200        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1201      }
1202      if (rank==1) {
1203        for (int i=0; i < shape[0]; i++) {
1204          value[i] = valueView(i);
1205        }
1206      }
1207      if (rank==2) {
1208        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1209      }
1210      if (rank==3) {
1211        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1212      }
1213      if (rank==4) {
1214        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1215      }
1216    
1217    }
1218    
1219  /*  /*
1220  Note: this version removed for now. Not needed, and breaks escript.cpp  Note: this version removed for now. Not needed, and breaks escript.cpp
1221  void  void
# Line 1050  Data::setTaggedValue(int tagKey, Line 1236  Data::setTaggedValue(int tagKey,
1236  }  }
1237  */  */
1238    
1239    void
1240    Data::archiveData(const std::string fileName)
1241    {
1242      cout << "Archiving Data object to: " << fileName << endl;
1243    
1244      //
1245      // Determine type of this Data object
1246      int dataType = -1;
1247    
1248      if (isEmpty()) {
1249        dataType = 0;
1250        cout << "\tdataType: DataEmpty" << endl;
1251      }
1252      if (isConstant()) {
1253        dataType = 1;
1254        cout << "\tdataType: DataConstant" << endl;
1255      }
1256      if (isTagged()) {
1257        dataType = 2;
1258        cout << "\tdataType: DataTagged" << endl;
1259      }
1260      if (isExpanded()) {
1261        dataType = 3;
1262        cout << "\tdataType: DataExpanded" << endl;
1263      }
1264      if (dataType == -1) {
1265        throw DataException("archiveData Error: undefined dataType");
1266      }
1267    
1268      //
1269      // Collect data items common to all Data types
1270      int noSamples = getNumSamples();
1271      int noDPPSample = getNumDataPointsPerSample();
1272      int functionSpaceType = getFunctionSpace().getTypeCode();
1273      int dataPointRank = getDataPointRank();
1274      int dataPointSize = getDataPointSize();
1275      int dataLength = getLength();
1276      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1277      int referenceNumbers[noSamples];
1278      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1279        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1280      }
1281      int tagNumbers[noSamples];
1282      if (isTagged()) {
1283        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1284          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1285        }
1286      }
1287    
1288      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1289      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1290      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1291    
1292      //
1293      // Flatten Shape to an array of integers suitable for writing to file
1294      int flatShape[4] = {0,0,0,0};
1295      cout << "\tshape: < ";
1296      for (int dim=0; dim<dataPointRank; dim++) {
1297        flatShape[dim] = dataPointShape[dim];
1298        cout << dataPointShape[dim] << " ";
1299      }
1300      cout << ">" << endl;
1301    
1302      //
1303      // Write common data items to archive file
1304      ofstream archiveFile;
1305      archiveFile.open(fileName.data(), ios::out);
1306    
1307      if (!archiveFile.good()) {
1308        throw DataException("archiveData Error: problem opening archive file");
1309      }
1310    
1311      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1312      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1313      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1314      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1315      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1316      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1317      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1318      for (int dim = 0; dim < 4; dim++) {
1319        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1320      }
1321      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1322        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1323      }
1324      if (isTagged()) {
1325        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1326          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1327        }
1328      }
1329    
1330      if (!archiveFile.good()) {
1331        throw DataException("archiveData Error: problem writing to archive file");
1332      }
1333    
1334      archiveFile.close();
1335    
1336      if (!archiveFile.good()) {
1337        throw DataException("archiveData Error: problem closing archive file");
1338      }
1339    
1340      //
1341      // Collect and archive underlying data values for each Data type
1342      switch (dataType) {
1343        case 0:
1344          // DataEmpty
1345          break;
1346        case 1:
1347          // DataConstant
1348          break;
1349        case 2:
1350          // DataTagged
1351          break;
1352        case 3:
1353          // DataExpanded
1354          break;
1355      }
1356    
1357    }
1358    
1359    void
1360    Data::extractData(const std::string fileName,
1361                      const FunctionSpace& fspace)
1362    {
1363      //
1364      // Can only extract Data to an object which is initially DataEmpty
1365      if (!isEmpty()) {
1366        throw DataException("extractData Error: can only extract to DataEmpty object");
1367      }
1368    
1369      cout << "Extracting Data object from: " << fileName << endl;
1370    
1371      int dataType;
1372      int noSamples;
1373      int noDPPSample;
1374      int functionSpaceType;
1375      int dataPointRank;
1376      int dataPointSize;
1377      int dataLength;
1378      DataArrayView::ShapeType dataPointShape;
1379      int flatShape[4];
1380    
1381      //
1382      // Open the archive file and read common data items
1383      ifstream archiveFile;
1384      archiveFile.open(fileName.data(), ios::in);
1385    
1386      if (!archiveFile.good()) {
1387        throw DataException("extractData Error: problem opening archive file");
1388      }
1389    
1390      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1391      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1392      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1393      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1394      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1395      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1396      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1397      for (int dim = 0; dim < 4; dim++) {
1398        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1399        if (flatShape[dim]>0) {
1400          dataPointShape.push_back(flatShape[dim]);
1401        }
1402      }
1403      int referenceNumbers[noSamples];
1404      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1405        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1406      }
1407      int tagNumbers[noSamples];
1408      if (dataType==2) {
1409        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1410          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1411        }
1412      }
1413    
1414      if (!archiveFile.good()) {
1415        throw DataException("extractData Error: problem reading from archive file");
1416      }
1417    
1418      archiveFile.close();
1419    
1420      if (!archiveFile.good()) {
1421        throw DataException("extractData Error: problem closing archive file");
1422      }
1423    
1424      switch (dataType) {
1425        case 0:
1426          cout << "\tdataType: DataEmpty" << endl;
1427          break;
1428        case 1:
1429          cout << "\tdataType: DataConstant" << endl;
1430          break;
1431        case 2:
1432          cout << "\tdataType: DataTagged" << endl;
1433          break;
1434        case 3:
1435          cout << "\tdataType: DataExpanded" << endl;
1436          break;
1437        default:
1438          throw DataException("extractData Error: undefined dataType read from archive file");
1439          break;
1440      }
1441    
1442      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1443      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1444      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1445      cout << "\tshape: < ";
1446      for (int dim = 0; dim < dataPointRank; dim++) {
1447        cout << dataPointShape[dim] << " ";
1448      }
1449      cout << ">" << endl;
1450    
1451      //
1452      // Verify that supplied FunctionSpace object is compatible with this Data object.
1453      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1454           (fspace.getNumSamples()!=noSamples) ||
1455           (fspace.getNumDPPSample()!=noDPPSample)
1456         ) {
1457        throw DataException("extractData Error: incompatible FunctionSpace");
1458      }
1459      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1460        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1461          throw DataException("extractData Error: incompatible FunctionSpace");
1462        }
1463      }
1464      if (dataType==2) {
1465        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1466          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1467            throw DataException("extractData Error: incompatible FunctionSpace");
1468          }
1469        }
1470      }
1471    
1472      //
1473      // Construct a DataVector to hold underlying data values
1474      DataVector dataVec(dataLength);
1475    
1476      //
1477      // Load this DataVector with the appropriate values
1478      switch (dataType) {
1479        case 0:
1480          // DataEmpty
1481          break;
1482        case 1:
1483          // DataConstant
1484          break;
1485        case 2:
1486          // DataTagged
1487          break;
1488        case 3:
1489          // DataExpanded
1490          break;
1491      }
1492    
1493      //
1494      // Construct an appropriate Data object
1495      DataAbstract* tempData;
1496      switch (dataType) {
1497        case 0:
1498          // DataEmpty
1499          tempData=new DataEmpty();
1500          break;
1501        case 1:
1502          // DataConstant
1503          tempData=new DataConstant(fspace,dataPointShape,dataVec);
1504          break;
1505        case 2:
1506          // DataTagged
1507          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1508          break;
1509        case 3:
1510          // DataExpanded
1511          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1512          break;
1513      }
1514      shared_ptr<DataAbstract> temp_data(tempData);
1515      m_data=temp_data;
1516    
1517    }
1518    
1519  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
1520  {  {
1521    o << data.toString();    o << data.toString();

Legend:
Removed from v.108  
changed lines
  Added in v.119

  ViewVC Help
Powered by ViewVC 1.1.26