/[escript]/trunk/bruce/src/Bruce/Bruce.cpp
ViewVC logotype

Diff of /trunk/bruce/src/Bruce/Bruce.cpp

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

revision 153 by jgs, Tue Oct 25 01:51:20 2005 UTC revision 154 by jgs, Mon Nov 7 05:51:17 2005 UTC
# Line 16  Line 16 
16  #include "bruce/Bruce/Bruce.h"  #include "bruce/Bruce/Bruce.h"
17  #include "bruce/Bruce/BruceException.h"  #include "bruce/Bruce/BruceException.h"
18    
19    #include "vtkCellType.h"
20    #include <boost/python/extract.hpp>
21    
22  using namespace std;  using namespace std;
23  using namespace escript;  using namespace escript;
24    
# Line 364  Bruce::getDataShape(int functionSpaceCod Line 367  Bruce::getDataShape(int functionSpaceCod
367    return pair<int,int>(numDataPointsPerSample,numSamples);    return pair<int,int>(numDataPointsPerSample,numSamples);
368  }  }
369    
370    int
371    Bruce::getNumSamples(int functionSpaceCode) const
372    {
373      std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
374      return domainShape.second;
375    }
376    
377  Data  Data
378  Bruce::getX() const  Bruce::getX() const
379  {  {
# Line 813  Bruce::setToSize(escript::Data& out) con Line 823  Bruce::setToSize(escript::Data& out) con
823    
824  }  }
825    
826    void
827    Bruce::setToGradient(escript::Data& grad,
828                         const escript::Data& arg) const
829    {
830      // pre-conditions: arg must be rank 0->3 only, rank 4 cannot be accomodated
831      //                 grad will be rank of arg, +1
832      // grad is calculated by, for each data-point in this Bruce domain, retrieving
833      // the associated value from arg, calculating the grad, and loading this value
834      // into the corresponding data-point in grad
835    }
836    
837  bool  bool
838  Bruce::operator==(const AbstractDomain& other) const  Bruce::operator==(const AbstractDomain& other) const
839  {  {
# Line 840  Bruce::operator!=(const AbstractDomain& Line 861  Bruce::operator!=(const AbstractDomain&
861  }  }
862    
863  int  int
864  Bruce::getTagFromSampleNo(int functionSpaceCode, int sampleNo) const  Bruce::getReferenceNoFromSampleNo(int functionSpaceCode,
865                                      int sampleNo) const
866  {  {
867    return 0;    // ensure supplied sampleNo is valid
868      std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
869      int numSamples = domainShape.second;
870      if ( (sampleNo>=numSamples) || (sampleNo < 0)) {
871        stringstream temp;
872        temp << "Bruce::getReferenceNoFromSampleNo: Error - invalid sampleNo supplied";
873        throw BruceException(temp.str());
874      }
875      // we just set the reference number to be the sample number for all samples
876      return sampleNo;
877  }  }
878    
879  int  void
880  Bruce::getReferenceNoFromSampleNo(int functionSpaceCode, int sampleNo) const  Bruce::saveVTK(const std::string& filename,
881                   const boost::python::dict& dataDict) const
882  {  {
883    return 0;  
884      //
885      // extract Data objects and associated names from dictionary object
886      const int num_data=boost::python::extract<int>(dataDict.attr("__len__")());
887    
888      int MAX_namelength=256;
889      char names[num_data][MAX_namelength];
890      char* c_names[num_data];
891    
892      escriptDataC data[num_data];
893      escriptDataC* ptr_data[num_data];
894    
895      boost::python::list keys=dataDict.keys();
896      for (int i=0; i<num_data; i++) {
897        Data& d=boost::python::extract<Data&>(dataDict[keys[i]]);
898        if (dynamic_cast<const Bruce&>(d.getFunctionSpace().getDomain()) != *this) {
899          throw BruceException("Error: Bruce::saveVTK: Data must be defined on same Domain");
900        }
901        data[i]=d.getDataC();
902        ptr_data[i]=&(data[i]);
903        std::string n=boost::python::extract<std::string>(keys[i]);
904        c_names[i]=&(names[i][0]);
905        if (n.length()>MAX_namelength-1) {
906          strncpy(c_names[i],n.c_str(),MAX_namelength-1);
907          c_names[i][MAX_namelength-1]='\0';
908        } else {
909          strcpy(c_names[i],n.c_str());
910        }
911      }
912    
913      //
914      // open archive file
915      ofstream archiveFile;
916      archiveFile.open(filename.data(), ios::out);
917      if (!archiveFile.good()) {
918          throw BruceException("Error in Bruce::saveVTK: File could not be opened for writing");
919      }
920    
921      //
922      // determine mesh type to be written
923      bool isCellCentered[num_data], write_celldata=false, write_pointdata=false;
924      for (int i_data=0; i_data<num_data; i_data++) {
925        if (!isEmpty(ptr_data[i_data])) {
926          switch(getFunctionSpaceType(ptr_data[i_data])) {
927            case ContinuousFunction:
928              isCellCentered[i_data]=false;
929              break;
930            case Function:
931              isCellCentered[i_data]=true;
932              break;
933          }
934          if (isCellCentered[i_data]) {
935            write_celldata=true;
936          } else {
937            write_pointdata=true;
938          }
939        }
940      }
941    
942      //
943      // determine number of points and cells
944      int numPoints=getDataShape(ContinuousFunction).second;
945      int numCells=getDataShape(Function).second;
946    
947      //
948      // determine VTK element type
949      int cellType;
950      char elemTypeStr[32];
951      int numVTKNodesPerElement;
952    
953      int nDim = getDim();
954      switch (nDim) {
955    
956        case 0:
957          cellType = VTK_VERTEX;
958          numVTKNodesPerElement = 1;
959          strcpy(elemTypeStr, "VTK_VERTEX");
960          break;
961    
962        case 1:
963          cellType = VTK_LINE;
964          numVTKNodesPerElement = 2;
965          strcpy(elemTypeStr, "VTK_LINE");
966          break;
967    
968        case 2:
969          cellType = VTK_QUAD;
970          numVTKNodesPerElement = 4;
971          strcpy(elemTypeStr, "VTK_QUAD");
972          break;
973    
974        case 3:
975          cellType = VTK_HEXAHEDRON;
976          numVTKNodesPerElement = 8;
977          strcpy(elemTypeStr, "VTK_HEXAHEDRON");
978          break;
979    
980      }
981    
982      //
983      // write xml header
984      archiveFile << "<?xml version=\"1.0\"?>" << endl;
985    
986      //
987      // determine grid extent
988    
989      // ??
990    
991      //
992      // write grid type and extent
993      archiveFile << "\t<VTKFile type=\"StructuredGrid\" version=\"0.1\">" << endl;
994      archiveFile << "\t\t<StructuredGrid WholeExtent=\"x1 x2 y1 y2 z1 z2\">" << endl;
995      archiveFile << "\t\t\t<Piece Extent=\"x1 x2 y1 y2 z1 z2\">" << endl;
996    
997      //
998      // start to write out point definitions
999      archiveFile << "\t\t\t\t<Points>" << endl;
1000    
1001      //
1002      // determine grid cooordinates
1003    
1004      // ??
1005    
1006      // vtk/mayavi doesn't like 2D data, it likes 3D data with
1007      // a degenerate third dimension to handle 2D data.
1008      // So if nDim is less than 3, must pad all empty dimensions,
1009      // so that the total number of dims is 3.
1010    
1011      archiveFile << "\t\t\t\t\t<DataArray NumberOfComponents=" << max(3,nDim) << " type=\"Float32\" format=\"ascii\">" << endl;
1012    /*
1013      for (int i=0; i<numPoints; i++) {
1014        for (int j=0; j<nDim; j++)
1015          archiveFile << "%e "; //mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)])
1016        for (int k=0; !(k>=3-nDim); k++)
1017          archiveFile << "0. ";
1018        archiveFile << endl;
1019      }
1020    */
1021      archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1022      archiveFile << "\t\t\t\t</Points>" << endl;
1023    
1024      //
1025      // cell data
1026      if (write_celldata) {
1027    
1028        //
1029        // mark the active cell-data data arrays
1030        bool set_scalar=false, set_vector=false, set_tensor=false;
1031        archiveFile << "\t\t\t\t<CellData";
1032        for (int i_data=0; i_data<num_data; i_data++) {
1033          if (!isEmpty(ptr_data[i_data]) && isCellCentered[i_data]) {
1034    
1035            switch(getDataPointRank(ptr_data[i_data])) {
1036    
1037              case 0:
1038                if (!set_scalar) {
1039                  archiveFile << " Scalars=" << names[i_data];
1040                  set_scalar=true;
1041                }
1042                break;
1043    
1044              case 1:
1045                if (!set_vector) {
1046                  archiveFile << " Vectors=" << names[i_data];
1047                  set_vector=true;
1048                }
1049                break;
1050    
1051              case 2:
1052                if (!set_tensor) {
1053                  archiveFile << " Tensors=" << names[i_data];
1054                  set_tensor=true;
1055                }
1056                break;
1057    
1058              default:
1059                archiveFile.close();
1060                throw BruceException("saveVTK: VTK can't handle objects with rank greater than 2.");
1061    
1062            }
1063          }
1064        }
1065        archiveFile << ">" << endl;
1066    
1067        //
1068        // write the cell-data data arrays
1069        for (int i_data=0; i_data<num_data; i_data++) {
1070          if (!isEmpty(ptr_data[i_data]) && isCellCentered[i_data]) {
1071    
1072            int numPointsPerSample=1;
1073            int rank=getDataPointRank(ptr_data[i_data]);
1074            int nComp=getDataPointSize(ptr_data[i_data]);
1075            int nCompReqd; // the number of components required by vtk
1076            int shape=0;
1077    
1078            switch (rank) {
1079    
1080              case 0:
1081                nCompReqd=1;
1082                break;
1083    
1084              case 1:
1085                shape=getDataPointShape(ptr_data[i_data], 0);
1086                if (shape>3) {
1087                  archiveFile.close();
1088                  throw BruceException("saveVTK: rank 1 object must have less than 4 components");
1089                }
1090                nCompReqd=3;
1091                break;
1092    
1093              case 2:
1094                shape=getDataPointShape(ptr_data[i_data], 0);
1095                if (shape>3 || shape!=getDataPointShape(ptr_data[i_data], 1)) {
1096                  archiveFile.close();
1097                  throw BruceException("saveVTK: rank 2 object must have less than 4x4 components and must have a square shape");
1098                }
1099                nCompReqd=9;
1100                break;
1101    
1102            }
1103    
1104            archiveFile << "\t\t\t\t\t<DataArray Name=" << names[i_data] << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1105    
1106    /*
1107            //
1108            // write out the cell data
1109    
1110            // if the number of required components is more than the number
1111            // of actual components, pad with zeros
1112            // probably only need to get shape of first element
1113            // write the data different ways for scalar, vector and tensor
1114    
1115            for (int i=0; i<numCells; i++) {
1116    
1117              double *values=getSampleData(ptr_data[i_data], i);
1118              double sampleAvg[nComp];
1119    
1120              // average over the number of points in the sample (ie: 1)
1121              for (int k=0; k<nComp; k++) {
1122                sampleAvg[k] = values[INDEX2(k,0,nComp)];
1123              }
1124    
1125              if (nCompReqd == 1) {
1126    
1127                archiveFile << sampleAvg[0] << endl;
1128    
1129              } else if (nCompReqd == 3) {
1130    
1131                for (int m=0; m<shape; m++) {
1132                  archiveFile << sampleAvg[m] << endl;
1133                }
1134                for (int m=0; m<nCompReqd-shape; m++) {
1135                  archiveFile << "0." << endl;
1136                }
1137    
1138              } else if (nCompReqd == 9) {
1139    
1140                int count=0;
1141                for (int m=0; m<shape; m++) {
1142                  for (int n=0; n<shape; n++) {
1143                    archiveFile << sampleAvg[count] << endl;
1144                    count++;
1145                  }
1146                  for (int n=0; n<3-shape; n++) {
1147                    archiveFile << "0." << endl;
1148                  }
1149                }
1150                for (int m=0; m<3-shape; m++) {
1151                  for (int n=0; n<3; n++) {
1152                    archiveFile << "0." << endl;
1153                  }
1154                }
1155    
1156              }
1157            }
1158    */
1159            archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1160    
1161          }
1162        }
1163    
1164        archiveFile << "\t\t\t\t</CellData>" << endl;
1165      }
1166    
1167      //
1168      // point data
1169      if (write_pointdata) {
1170    
1171        //
1172        // mark the active point-data data arrays
1173        bool set_scalar=false, set_vector=false, set_tensor=false;
1174        archiveFile << "\t\t\t\t<PointData";
1175        for (int i_data=0; i_data<num_data; i_data++) {
1176          if (!isEmpty(ptr_data[i_data]) && !isCellCentered[i_data]) {
1177    
1178            switch(getDataPointRank(ptr_data[i_data])) {
1179    
1180              case 0:
1181                if (!set_scalar) {
1182                  archiveFile << " Scalars=" << names[i_data];
1183                  set_scalar=true;
1184                }
1185                break;
1186    
1187              case 1:
1188                if (!set_vector) {
1189                  archiveFile << " Vectors=" << names[i_data];
1190                  set_vector=true;
1191                }
1192                break;
1193    
1194              case 2:
1195                if (!set_tensor) {
1196                  archiveFile << " Tensors=" << names[i_data];
1197                  set_tensor=true;
1198                }
1199                break;
1200    
1201              default:
1202                archiveFile.close();
1203                throw BruceException("saveVTK: Vtk can't handle objects with rank greater than 2");
1204            }
1205    
1206          }
1207        }
1208        archiveFile << ">" << endl;
1209    
1210        //
1211        // write the point-data data arrays
1212        for (int i_data=0; i_data<num_data; i_data++) {
1213          if (!isEmpty(ptr_data[i_data]) && !isCellCentered[i_data]) {
1214    
1215            int numPointsPerSample=1;
1216            int rank=getDataPointRank(ptr_data[i_data]);
1217            int nComp=getDataPointSize(ptr_data[i_data]);
1218            int nCompReqd; // the number of components required by vtk
1219            int shape=0;
1220    
1221            switch (rank) {
1222    
1223            case 0:
1224              nCompReqd=1;
1225    
1226            case 1:
1227              shape=getDataPointShape(ptr_data[i_data], 0);
1228              if (shape>3) {
1229                archiveFile.close();
1230                throw BruceException("saveVTK: rank 1 object must have less than 4 components");
1231              }
1232              nCompReqd=3;
1233    
1234            case 2:
1235              shape=getDataPointShape(ptr_data[i_data], 0);
1236              if (shape>3 || shape!=getDataPointShape(ptr_data[i_data], 1)) {
1237                archiveFile.close();
1238                throw BruceException("saveVTK: rank 2 object must have less than 4x4 components and must have a square shape");
1239              }
1240              nCompReqd=9;
1241    
1242            }
1243    
1244            archiveFile << "\t\t\t\t\t<DataArray Name=" << names[i_data] << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1245    
1246    /*
1247            //
1248            // write out the point data
1249    
1250            // if the number of required components is more than
1251            // the number of actual components, pad with zeros
1252            // write the data different ways for scalar, vector and tensor
1253    
1254            for (int i=0; i<numNodes; i++) {
1255    
1256              values=getSampleData(data[i_data], i);
1257    
1258              if (nCompReqd==1) {
1259    
1260                archiveFile << values[0];
1261    
1262              } else if (nCompReqd==3) {
1263    
1264                for (int m=0; m<shape; m++) {
1265                  archiveFile << values[m];
1266                }
1267                for (int m=0; m<nCompReqd-shape; m++) {
1268                  archiveFile << "0.";
1269                }
1270    
1271              } else if (nCompReqd==9) {
1272    
1273                 int count=0;
1274                 for (int m=0; m<shape; m++) {
1275                   for (int n=0; n<shape; n++) {
1276                     archiveFile << values[count];
1277                     count++;
1278                   }
1279                   for (int n=0; n<3-shape; n++) {
1280                     archiveFile << "0.";
1281                   }
1282                 }
1283                 for (int m=0; m<3-shape; m++) {
1284                   for (int n=0; n<3; n++) {
1285                     archiveFile << "0.";
1286                   }
1287                 }
1288    
1289              }
1290    
1291              archiveFile << endl;
1292            }
1293    */
1294    
1295            archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1296          }
1297        }
1298    
1299        archiveFile << "\t\t\t\t</PointData>" << endl;
1300      }
1301    
1302      //
1303      // finish off the grid definition
1304      archiveFile << "\t\t\t</Piece>" << endl;
1305      archiveFile << "\t\t</StructuredGrid>" << endl;
1306    
1307      //
1308      // write the xml footer
1309      archiveFile << "\t</VTKFile>" << endl;
1310    
1311      //
1312      // Close archive file
1313      archiveFile.close();
1314    
1315      if (!archiveFile.good()) {
1316        throw BruceException("Error in Bruce::saveVTK: problem closing file");
1317      }
1318    
1319    }
1320    
1321    void
1322    Bruce::interpolateOnDomain(escript::Data& target,
1323                               const escript::Data& source) const
1324    {
1325    }
1326    
1327    bool
1328    Bruce::probeInterpolationOnDomain(int functionSpaceType_source,
1329                                      int functionSpaceType_target) const
1330    {
1331      return true;
1332    }
1333    
1334    void
1335    Bruce::interpolateACross(escript::Data& target,
1336                             const escript::Data& source) const
1337    {
1338    }
1339    
1340    bool
1341    Bruce::probeInterpolationACross(int functionSpaceType_source,
1342                                    const AbstractDomain& targetDomain,
1343                                    int functionSpaceType_target) const
1344    {
1345      return true;
1346  }  }
1347    
1348  bool  bool

Legend:
Removed from v.153  
changed lines
  Added in v.154

  ViewVC Help
Powered by ViewVC 1.1.26