/[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 97 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 122 by jgs, Thu Jun 9 05:38: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"
36  #include "escript/Data/DataEmpty.h"  #include "escript/Data/DataEmpty.h"
37  #include "escript/Data/DataArray.h"  #include "escript/Data/DataArray.h"
 #include "escript/Data/DataAlgorithm.h"  
38  #include "escript/Data/FunctionSpaceFactory.h"  #include "escript/Data/FunctionSpaceFactory.h"
39  #include "escript/Data/AbstractContinuousDomain.h"  #include "escript/Data/AbstractContinuousDomain.h"
40  #include "escript/Data/UnaryFuncs.h"  #include "escript/Data/UnaryFuncs.h"
# Line 185  Data::getDataC() const Line 184  Data::getDataC() const
184    return temp;    return temp;
185  }  }
186    
187  tuple  const boost::python::tuple
188  Data::getShapeTuple() const  Data::getShapeTuple() const
189  {  {
190    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 443  Data::getDataPointShape() const Line 442  Data::getDataPointShape() const
442    return getPointDataView().getShape();    return getPointDataView().getShape();
443  }  }
444    
445    const
446    boost::python::numeric::array
447    Data::convertToNumArray()
448    {
449      //
450      // determine the total 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          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
503          if (dataPointRank==0) {
504            numArray[dataPoint]=dataPointView();
505          }
506          if (dataPointRank==1) {
507            for (int i=0; i<dataPointShape[0]; i++) {
508              numArray[dataPoint][i]=dataPointView(i);
509            }
510          }
511          if (dataPointRank==2) {
512            for (int i=0; i<dataPointShape[0]; i++) {
513              for (int j=0; j<dataPointShape[1]; j++) {
514                numArray[dataPoint][i][j] = dataPointView(i,j);
515              }
516            }
517          }
518          if (dataPointRank==3) {
519            for (int i=0; i<dataPointShape[0]; i++) {
520              for (int j=0; j<dataPointShape[1]; j++) {
521                for (int k=0; k<dataPointShape[2]; k++) {
522                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
523                }
524              }
525            }
526          }
527          if (dataPointRank==4) {
528            for (int i=0; i<dataPointShape[0]; i++) {
529              for (int j=0; j<dataPointShape[1]; j++) {
530                for (int k=0; k<dataPointShape[2]; k++) {
531                  for (int l=0; l<dataPointShape[3]; l++) {
532                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
533                  }
534                }
535              }
536            }
537          }
538          dataPoint++;
539        }
540      }
541    
542      //
543      // return the loaded array
544      return numArray;
545    }
546    
547    const
548    boost::python::numeric::array
549    Data::convertToNumArrayFromSampleNo(int sampleNo)
550    {
551      //
552      // Check a valid sample number has been supplied
553      if (sampleNo >= getNumSamples()) {
554        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
555      }
556    
557      //
558      // determine the number of data points per sample
559      int numDataPointsPerSample = getNumDataPointsPerSample();
560    
561      //
562      // determine the rank and shape of each data point
563      int dataPointRank = getDataPointRank();
564      DataArrayView::ShapeType dataPointShape = getDataPointShape();
565    
566      //
567      // create the numeric array to be returned
568      boost::python::numeric::array numArray(0.0);
569    
570      //
571      // the rank of the returned numeric array will be the rank of
572      // the data points, plus one. Where the rank of the array is n,
573      // the last n-1 dimensions will be equal to the shape of the
574      // data points, whilst the first dimension will be equal to the
575      // total number of data points. Thus the array will consist of
576      // a serial vector of the data points.
577      int arrayRank = dataPointRank + 1;
578      DataArrayView::ShapeType arrayShape;
579      arrayShape.push_back(numDataPointsPerSample);
580      for (int d=0; d<dataPointRank; d++) {
581         arrayShape.push_back(dataPointShape[d]);
582      }
583    
584      //
585      // resize the numeric array to the shape just calculated
586      if (arrayRank==1) {
587        numArray.resize(arrayShape[0]);
588      }
589      if (arrayRank==2) {
590        numArray.resize(arrayShape[0],arrayShape[1]);
591      }
592      if (arrayRank==3) {
593        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
594      }
595      if (arrayRank==4) {
596        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
597      }
598      if (arrayRank==5) {
599        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
600      }
601    
602      //
603      // loop through each data point in turn, loading the values for that data point
604      // into the numeric array.
605      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
606        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
607        if (dataPointRank==0) {
608          numArray[dataPoint]=dataPointView();
609        }
610        if (dataPointRank==1) {
611          for (int i=0; i<dataPointShape[0]; i++) {
612            numArray[dataPoint][i]=dataPointView(i);
613          }
614        }
615        if (dataPointRank==2) {
616          for (int i=0; i<dataPointShape[0]; i++) {
617            for (int j=0; j<dataPointShape[1]; j++) {
618              numArray[dataPoint][i][j] = dataPointView(i,j);
619            }
620          }
621        }
622        if (dataPointRank==3) {
623          for (int i=0; i<dataPointShape[0]; i++) {
624            for (int j=0; j<dataPointShape[1]; j++) {
625              for (int k=0; k<dataPointShape[2]; k++) {
626                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
627              }
628            }
629          }
630        }
631        if (dataPointRank==4) {
632          for (int i=0; i<dataPointShape[0]; i++) {
633            for (int j=0; j<dataPointShape[1]; j++) {
634              for (int k=0; k<dataPointShape[2]; k++) {
635                for (int l=0; l<dataPointShape[3]; l++) {
636                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
637                }
638              }
639            }
640          }
641        }
642      }
643    
644      //
645      // return the loaded array
646      return numArray;
647    }
648    
649    const
650    boost::python::numeric::array
651    Data::convertToNumArrayFromDPNo(int sampleNo,
652                                    int dataPointNo)
653    {
654      //
655      // Check a valid sample number has been supplied
656      if (sampleNo >= getNumSamples()) {
657        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
658      }
659    
660      //
661      // Check a valid data point number has been supplied
662      if (dataPointNo >= getNumDataPointsPerSample()) {
663        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
664      }
665    
666      //
667      // determine the rank and shape of each data point
668      int dataPointRank = getDataPointRank();
669      DataArrayView::ShapeType dataPointShape = getDataPointShape();
670    
671      //
672      // create the numeric array to be returned
673      boost::python::numeric::array numArray(0.0);
674    
675      //
676      // the shape of the returned numeric array will be the same
677      // as that of the data point
678      int arrayRank = dataPointRank;
679      DataArrayView::ShapeType arrayShape = dataPointShape;
680    
681      //
682      // resize the numeric array to the shape just calculated
683      if (arrayRank==0) {
684        numArray.resize(1);
685      }
686      if (arrayRank==1) {
687        numArray.resize(arrayShape[0]);
688      }
689      if (arrayRank==2) {
690        numArray.resize(arrayShape[0],arrayShape[1]);
691      }
692      if (arrayRank==3) {
693        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
694      }
695      if (arrayRank==4) {
696        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
697      }
698    
699      //
700      // load the values for the data point into the numeric array.
701      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
702      if (dataPointRank==0) {
703        numArray[0]=dataPointView();
704      }
705      if (dataPointRank==1) {
706        for (int i=0; i<dataPointShape[0]; i++) {
707          numArray[i]=dataPointView(i);
708        }
709      }
710      if (dataPointRank==2) {
711        for (int i=0; i<dataPointShape[0]; i++) {
712          for (int j=0; j<dataPointShape[1]; j++) {
713            numArray[i][j] = dataPointView(i,j);
714          }
715        }
716      }
717      if (dataPointRank==3) {
718        for (int i=0; i<dataPointShape[0]; i++) {
719          for (int j=0; j<dataPointShape[1]; j++) {
720            for (int k=0; k<dataPointShape[2]; k++) {
721              numArray[i][j][k]=dataPointView(i,j,k);
722            }
723          }
724        }
725      }
726      if (dataPointRank==4) {
727        for (int i=0; i<dataPointShape[0]; i++) {
728          for (int j=0; j<dataPointShape[1]; j++) {
729            for (int k=0; k<dataPointShape[2]; k++) {
730              for (int l=0; l<dataPointShape[3]; l++) {
731                numArray[i][j][k][l]=dataPointView(i,j,k,l);
732              }
733            }
734          }
735        }
736      }
737    
738      //
739      // return the loaded array
740      return numArray;
741    }
742    
743  boost::python::numeric::array  boost::python::numeric::array
744  Data::integrate() const  Data::integrate() const
745  {  {
# Line 460  Data::integrate() const Line 757  Data::integrate() const
757    // and load the array with the integral values    // and load the array with the integral values
758    boost::python::numeric::array bp_array(1.0);    boost::python::numeric::array bp_array(1.0);
759    if (rank==0) {    if (rank==0) {
760        bp_array.resize(1);
761      index = 0;      index = 0;
762      bp_array[0] = integrals[index];      bp_array[0] = integrals[index];
763    }    }
# Line 539  Data::ln() const Line 837  Data::ln() const
837    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
838  }  }
839    
840    Data
841    Data::sign() const
842    {
843      return escript::unaryOp(*this,escript::fsign);
844    }
845    
846    Data
847    Data::abs() const
848    {
849      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
850    }
851    
852    Data
853    Data::neg() const
854    {
855      return escript::unaryOp(*this,negate<double>());
856    }
857    
858    Data
859    Data::pos() const
860    {
861      return (*this);
862    }
863    
864    Data
865    Data::exp() const
866    {
867      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
868    }
869    
870    Data
871    Data::sqrt() const
872    {
873      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
874    }
875    
876  double  double
877  Data::Lsup() const  Data::Lsup() const
878  {  {
# Line 548  Data::Lsup() const Line 882  Data::Lsup() const
882  }  }
883    
884  double  double
885    Data::Linf() const
886    {
887      //
888      // set the initial absolute minimum value to max double
889      return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));
890    }
891    
892    double
893  Data::sup() const  Data::sup() const
894  {  {
895    //    //
896    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
897    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::min()));    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
898  }  }
899    
900  double  double
# Line 566  Data::inf() const Line 908  Data::inf() const
908  Data  Data
909  Data::maxval() const  Data::maxval() const
910  {  {
911    // not implemented - will use dp_algorithm    //
912    return (*this);    // set the initial maximum value to min possible double
913      return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
914  }  }
915    
916  Data  Data
917  Data::minval() const  Data::minval() const
918  {  {
919    // not implemented - will use dp_algorithm    //
920    return (*this);    // set the initial minimum value to max possible double
921      return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
922  }  }
923    
924  Data  const boost::python::tuple
925  Data::length() const  Data::mindp() const
926  {  {
927    // not implemented - will use dp_algorithm    Data temp=minval();
   return (*this);  
 }  
928    
929  Data    int numSamples=temp.getNumSamples();
930  Data::trace() const    int numDPPSample=temp.getNumDataPointsPerSample();
 {  
   // not implemented - will use dp_algorithm  
   return (*this);  
 }  
931    
932  Data    int i,j,lowi=0,lowj=0;
933  Data::transpose(int axis) const    double min=numeric_limits<double>::max();
 {  
   // not implemented  
   return (*this);  
 }  
934    
935  Data    for (i=0; i<numSamples; i++) {
936  Data::sign() const      for (j=0; j<numDPPSample; j++) {
937  {        double next=temp.getDataPoint(i,j)();
938    return escript::unaryOp(*this,escript::fsign);        if (next<min) {
939            min=next;
940            lowi=i;
941            lowj=j;
942          }
943        }
944      }
945    
946      return make_tuple(lowi,lowj);
947  }  }
948    
949  Data  Data
950  Data::abs() const  Data::length() const
951  {  {
952    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return dp_algorithm(DataAlgorithmAdapter<Length>(0));
953  }  }
954    
955  Data  Data
956  Data::neg() const  Data::trace() const
957  {  {
958    return escript::unaryOp(*this,negate<double>());    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));
959  }  }
960    
961  Data  Data
962  Data::pos() const  Data::transpose(int axis) const
963  {  {
964    return (*this);    // not implemented
965      throw DataException("Error - Data::transpose not implemented yet.");
966      return Data();
967  }  }
968    
969  Data  void
970  Data::exp() const  Data::saveDX(std::string fileName) const
971  {  {
972    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    getDomain().saveDX(fileName,*this);
973      return;
974  }  }
975    
976  Data  void
977  Data::sqrt() const  Data::saveVTK(std::string fileName) const
978  {  {
979    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    getDomain().saveVTK(fileName,*this);
980      return;
981  }  }
982    
983  Data&  Data&
# Line 959  Data::setItemD(const boost::python::obje Line 1305  Data::setItemD(const boost::python::obje
1305    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1306      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1307    }    }
1308    setSlice(value,slice_region);    if (getFunctionSpace()!=value.getFunctionSpace()) {
1309         setSlice(Data(value,getFunctionSpace()),slice_region);
1310      } else {
1311         setSlice(value,slice_region);
1312      }
1313  }  }
1314    
1315  void  void
# Line 1021  Data::setTaggedValue(int tagKey, Line 1371  Data::setTaggedValue(int tagKey,
1371    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1372  }  }
1373    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1374  void  void
1375  Data::setTaggedValue(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1376                       const DataArrayView& value)                              const DataArrayView& value)
1377  {  {
1378    //    //
1379    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
# Line 1039  Data::setTaggedValue(int tagKey, Line 1387  Data::setTaggedValue(int tagKey,
1387    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1388    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1389  }  }
1390  */  
1391    void
1392    Data::setRefValue(int ref,
1393                      const boost::python::numeric::array& value)
1394    {
1395      //
1396      // Construct DataArray from boost::python::object input value
1397      DataArray valueDataArray(value);
1398    
1399      //
1400      // Call DataAbstract::setRefValue
1401      m_data->setRefValue(ref,valueDataArray);
1402    }
1403    
1404    void
1405    Data::getRefValue(int ref,
1406                      boost::python::numeric::array& value)
1407    {
1408      //
1409      // Construct DataArray for boost::python::object return value
1410      DataArray valueDataArray(value);
1411    
1412      //
1413      // Load DataArray with values from data-points specified by ref
1414      m_data->getRefValue(ref,valueDataArray);
1415    
1416      //
1417      // Load values from valueDataArray into return numarray
1418    
1419      // extract the shape of the numarray
1420      int rank = value.getrank();
1421      DataArrayView::ShapeType shape;
1422      for (int i=0; i < rank; i++) {
1423        shape.push_back(extract<int>(value.getshape()[i]));
1424      }
1425    
1426      // and load the numarray with the data from the DataArray
1427      DataArrayView valueView = valueDataArray.getView();
1428    
1429      if (rank==0) {
1430        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1431      }
1432      if (rank==1) {
1433        for (int i=0; i < shape[0]; i++) {
1434          value[i] = valueView(i);
1435        }
1436      }
1437      if (rank==2) {
1438        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1439      }
1440      if (rank==3) {
1441        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1442      }
1443      if (rank==4) {
1444        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1445      }
1446    
1447    }
1448    
1449    void
1450    Data::archiveData(const std::string fileName)
1451    {
1452      cout << "Archiving Data object to: " << fileName << endl;
1453    
1454      //
1455      // Determine type of this Data object
1456      int dataType = -1;
1457    
1458      if (isEmpty()) {
1459        dataType = 0;
1460        cout << "\tdataType: DataEmpty" << endl;
1461      }
1462      if (isConstant()) {
1463        dataType = 1;
1464        cout << "\tdataType: DataConstant" << endl;
1465      }
1466      if (isTagged()) {
1467        dataType = 2;
1468        cout << "\tdataType: DataTagged" << endl;
1469      }
1470      if (isExpanded()) {
1471        dataType = 3;
1472        cout << "\tdataType: DataExpanded" << endl;
1473      }
1474      if (dataType == -1) {
1475        throw DataException("archiveData Error: undefined dataType");
1476      }
1477    
1478      //
1479      // Collect data items common to all Data types
1480      int noSamples = getNumSamples();
1481      int noDPPSample = getNumDataPointsPerSample();
1482      int functionSpaceType = getFunctionSpace().getTypeCode();
1483      int dataPointRank = getDataPointRank();
1484      int dataPointSize = getDataPointSize();
1485      int dataLength = getLength();
1486      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1487      int referenceNumbers[noSamples];
1488      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1489        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1490      }
1491      int tagNumbers[noSamples];
1492      if (isTagged()) {
1493        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1494          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1495        }
1496      }
1497    
1498      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1499      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1500      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1501    
1502      //
1503      // Flatten Shape to an array of integers suitable for writing to file
1504      int flatShape[4] = {0,0,0,0};
1505      cout << "\tshape: < ";
1506      for (int dim=0; dim<dataPointRank; dim++) {
1507        flatShape[dim] = dataPointShape[dim];
1508        cout << dataPointShape[dim] << " ";
1509      }
1510      cout << ">" << endl;
1511    
1512      //
1513      // Write common data items to archive file
1514      ofstream archiveFile;
1515      archiveFile.open(fileName.data(), ios::out);
1516    
1517      if (!archiveFile.good()) {
1518        throw DataException("archiveData Error: problem opening archive file");
1519      }
1520    
1521      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1522      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1523      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1524      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1525      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1526      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1527      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1528      for (int dim = 0; dim < 4; dim++) {
1529        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1530      }
1531      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1532        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1533      }
1534      if (isTagged()) {
1535        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1536          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1537        }
1538      }
1539    
1540      if (!archiveFile.good()) {
1541        throw DataException("archiveData Error: problem writing to archive file");
1542      }
1543    
1544      archiveFile.close();
1545    
1546      if (!archiveFile.good()) {
1547        throw DataException("archiveData Error: problem closing archive file");
1548      }
1549    
1550      //
1551      // Collect and archive underlying data values for each Data type
1552      switch (dataType) {
1553        case 0:
1554          // DataEmpty
1555          break;
1556        case 1:
1557          // DataConstant
1558          break;
1559        case 2:
1560          // DataTagged
1561          break;
1562        case 3:
1563          // DataExpanded
1564          break;
1565      }
1566    
1567    }
1568    
1569    void
1570    Data::extractData(const std::string fileName,
1571                      const FunctionSpace& fspace)
1572    {
1573      //
1574      // Can only extract Data to an object which is initially DataEmpty
1575      if (!isEmpty()) {
1576        throw DataException("extractData Error: can only extract to DataEmpty object");
1577      }
1578    
1579      cout << "Extracting Data object from: " << fileName << endl;
1580    
1581      int dataType;
1582      int noSamples;
1583      int noDPPSample;
1584      int functionSpaceType;
1585      int dataPointRank;
1586      int dataPointSize;
1587      int dataLength;
1588      DataArrayView::ShapeType dataPointShape;
1589      int flatShape[4];
1590    
1591      //
1592      // Open the archive file and read common data items
1593      ifstream archiveFile;
1594      archiveFile.open(fileName.data(), ios::in);
1595    
1596      if (!archiveFile.good()) {
1597        throw DataException("extractData Error: problem opening archive file");
1598      }
1599    
1600      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1601      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1602      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1603      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1604      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1605      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1606      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1607      for (int dim = 0; dim < 4; dim++) {
1608        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1609        if (flatShape[dim]>0) {
1610          dataPointShape.push_back(flatShape[dim]);
1611        }
1612      }
1613      int referenceNumbers[noSamples];
1614      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1615        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1616      }
1617      int tagNumbers[noSamples];
1618      if (dataType==2) {
1619        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1620          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1621        }
1622      }
1623    
1624      if (!archiveFile.good()) {
1625        throw DataException("extractData Error: problem reading from archive file");
1626      }
1627    
1628      archiveFile.close();
1629    
1630      if (!archiveFile.good()) {
1631        throw DataException("extractData Error: problem closing archive file");
1632      }
1633    
1634      switch (dataType) {
1635        case 0:
1636          cout << "\tdataType: DataEmpty" << endl;
1637          break;
1638        case 1:
1639          cout << "\tdataType: DataConstant" << endl;
1640          break;
1641        case 2:
1642          cout << "\tdataType: DataTagged" << endl;
1643          break;
1644        case 3:
1645          cout << "\tdataType: DataExpanded" << endl;
1646          break;
1647        default:
1648          throw DataException("extractData Error: undefined dataType read from archive file");
1649          break;
1650      }
1651    
1652      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1653      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1654      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1655      cout << "\tshape: < ";
1656      for (int dim = 0; dim < dataPointRank; dim++) {
1657        cout << dataPointShape[dim] << " ";
1658      }
1659      cout << ">" << endl;
1660    
1661      //
1662      // Verify that supplied FunctionSpace object is compatible with this Data object.
1663      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1664           (fspace.getNumSamples()!=noSamples) ||
1665           (fspace.getNumDPPSample()!=noDPPSample)
1666         ) {
1667        throw DataException("extractData Error: incompatible FunctionSpace");
1668      }
1669      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1670        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1671          throw DataException("extractData Error: incompatible FunctionSpace");
1672        }
1673      }
1674      if (dataType==2) {
1675        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1676          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1677            throw DataException("extractData Error: incompatible FunctionSpace");
1678          }
1679        }
1680      }
1681    
1682      //
1683      // Construct a DataVector to hold underlying data values
1684      DataVector dataVec(dataLength);
1685    
1686      //
1687      // Load this DataVector with the appropriate values
1688      switch (dataType) {
1689        case 0:
1690          // DataEmpty
1691          break;
1692        case 1:
1693          // DataConstant
1694          break;
1695        case 2:
1696          // DataTagged
1697          break;
1698        case 3:
1699          // DataExpanded
1700          break;
1701      }
1702    
1703      //
1704      // Construct an appropriate Data object
1705      DataAbstract* tempData;
1706      switch (dataType) {
1707        case 0:
1708          // DataEmpty
1709          tempData=new DataEmpty();
1710          break;
1711        case 1:
1712          // DataConstant
1713          tempData=new DataConstant(fspace,dataPointShape,dataVec);
1714          break;
1715        case 2:
1716          // DataTagged
1717          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1718          break;
1719        case 3:
1720          // DataExpanded
1721          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1722          break;
1723      }
1724      shared_ptr<DataAbstract> temp_data(tempData);
1725      m_data=temp_data;
1726    
1727    }
1728    
1729  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
1730  {  {

Legend:
Removed from v.97  
changed lines
  Added in v.122

  ViewVC Help
Powered by ViewVC 1.1.26