/[escript]/trunk/finley/src/Mesh_saveVTK.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_saveVTK.c

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

revision 818 by dhawcroft, Sun Aug 27 11:10:34 2006 UTC revision 1030 by phornby, Wed Mar 14 05:14:44 2007 UTC
# Line 452  void Finley_Mesh_saveVTK_MPIO(const char Line 452  void Finley_Mesh_saveVTK_MPIO(const char
452              "<UnstructuredGrid>\n" \              "<UnstructuredGrid>\n" \
453              "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \              "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
454              "<Points>\n" \              "<Points>\n" \
455              "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n"              "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n"
456              ,numPoints,numGlobalCells,MAX(3,nDim));              ,numPoints,numGlobalCells,MAX(3,nDim));
457    
458    
# Line 845  void Finley_Mesh_saveVTK_MPIO(const char Line 845  void Finley_Mesh_saveVTK_MPIO(const char
845          if( myRank == 0)          if( myRank == 0)
846          {          {
847            char header[250];            char header[250];
848            sprintf(header,"<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);            sprintf(header,"<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
849            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
850            MPI_Wait(&req,&status);            MPI_Wait(&req,&status);
851          }          }
# Line 866  void Finley_Mesh_saveVTK_MPIO(const char Line 866  void Finley_Mesh_saveVTK_MPIO(const char
866            // averaging over the number of points in the sample            // averaging over the number of points in the sample
867            for (n=0; n<nComp; n++)            for (n=0; n<nComp; n++)
868            {            {
869              rtmp = 0.;              if (isExpanded(data_pp[i_data])) {
870              for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];                 rtmp = 0.;
871              sampleAvg[k] = rtmp/numPointsPerSample;                 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];
872                   sampleAvg[n] = rtmp/numPointsPerSample;
873                } else {
874                   sampleAvg[n] = values[n];
875                }
876            }            }
877            // if the number of required components is more than the number            // if the number of required components is more than the number
878            // of actual components, pad with zeros            // of actual components, pad with zeros
# Line 1042  void Finley_Mesh_saveVTK_MPIO(const char Line 1046  void Finley_Mesh_saveVTK_MPIO(const char
1046          if( myRank == 0)          if( myRank == 0)
1047          {          {
1048            char header[250];            char header[250];
1049            sprintf(header, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);            sprintf(header, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1050            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1051            MPI_Wait(&req,&status);            MPI_Wait(&req,&status);
1052          }          }
# Line 1199  void Finley_Mesh_saveVTK_MPIO(const char Line 1203  void Finley_Mesh_saveVTK_MPIO(const char
1203    
1204  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
1205  {  {
1206      #define NCOMP_MAX 9
1207    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1208      double sampleAvg[NCOMP_MAX];
1209    /* if there is no mesh we just return */    /* if there is no mesh we just return */
   if (mesh_p==NULL) return;  
1210    
1211    int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,    int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1212    nDim, numPointsPerSample, nComp, nCompReqd;    nDim, numPointsPerSample, nComp, nCompReqd, NN;
   
1213    index_t j2;    index_t j2;
1214      int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1215      int elementtype=FINLEY_UNKNOWN;
1216    double* values, rtmp;    double* values, rtmp;
1217    char elemTypeStr[32];    char elemTypeStr[32];
1218      FILE * fileHandle_p = NULL;
1219      bool_t do_write, *isCellCentered=NULL,write_celldata=FALSE,write_pointdata=FALSE;
1220      Finley_ElementFile* elements=NULL;
1221      ElementTypeId TypeId;
1222    
1223    /* open the file and check handle */    /* open the file and check handle */
1224      if (mesh_p==NULL) return;
1225    
1226    FILE * fileHandle_p = fopen(filename_p, "w");    fileHandle_p = fopen(filename_p, "w");
1227    if (fileHandle_p==NULL)    if (fileHandle_p==NULL)
1228    {    {
1229      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
# Line 1220  void Finley_Mesh_saveVTK(const char * fi Line 1231  void Finley_Mesh_saveVTK(const char * fi
1231      return;      return;
1232    }    }
1233    /* find the mesh type to be written */    /* find the mesh type to be written */
1234    int nodetype=FINLEY_DEGREES_OF_FREEDOM;    isCellCentered=TMPMEMALLOC(num_data,bool_t);
1235    int elementtype=FINLEY_UNKNOWN;  
1236    bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;  
1237      if (Finley_checkPtr(isCellCentered)) {
1238         fclose(fileHandle_p);
1239         return;
1240      }
1241    for (i_data=0;i_data<num_data;++i_data)    for (i_data=0;i_data<num_data;++i_data)
1242    {    {
1243      if (! isEmpty(data_pp[i_data]))      if (! isEmpty(data_pp[i_data]))
# Line 1239  void Finley_Mesh_saveVTK(const char * fi Line 1254  void Finley_Mesh_saveVTK(const char * fi
1254          {          {
1255            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1256            fclose(fileHandle_p);            fclose(fileHandle_p);
1257              TMPMEMFREE(isCellCentered);
1258            return;            return;
1259          }          }
1260          isCellCentered[i_data]=FALSE;          isCellCentered[i_data]=FALSE;
# Line 1253  void Finley_Mesh_saveVTK(const char * fi Line 1269  void Finley_Mesh_saveVTK(const char * fi
1269          {          {
1270            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1271            fclose(fileHandle_p);            fclose(fileHandle_p);
1272              TMPMEMFREE(isCellCentered);
1273            return;            return;
1274          }          }
1275          isCellCentered[i_data]=FALSE;          isCellCentered[i_data]=FALSE;
# Line 1267  void Finley_Mesh_saveVTK(const char * fi Line 1284  void Finley_Mesh_saveVTK(const char * fi
1284          {          {
1285            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1286            fclose(fileHandle_p);            fclose(fileHandle_p);
1287              TMPMEMFREE(isCellCentered);
1288            return;            return;
1289          }          }
1290          isCellCentered[i_data]=FALSE;          isCellCentered[i_data]=FALSE;
# Line 1282  void Finley_Mesh_saveVTK(const char * fi Line 1300  void Finley_Mesh_saveVTK(const char * fi
1300            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1301            fclose(fileHandle_p);            fclose(fileHandle_p);
1302            return;            return;
1303              TMPMEMFREE(isCellCentered);
1304          }          }
1305          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
1306          break;          break;
# Line 1295  void Finley_Mesh_saveVTK(const char * fi Line 1314  void Finley_Mesh_saveVTK(const char * fi
1314          {          {
1315            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1316            fclose(fileHandle_p);            fclose(fileHandle_p);
1317              TMPMEMFREE(isCellCentered);
1318            return;            return;
1319          }          }
1320          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1309  void Finley_Mesh_saveVTK(const char * fi Line 1329  void Finley_Mesh_saveVTK(const char * fi
1329          {          {
1330            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1331            fclose(fileHandle_p);            fclose(fileHandle_p);
1332              TMPMEMFREE(isCellCentered);
1333            return;            return;
1334          }          }
1335          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1323  void Finley_Mesh_saveVTK(const char * fi Line 1344  void Finley_Mesh_saveVTK(const char * fi
1344          {          {
1345            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1346            fclose(fileHandle_p);            fclose(fileHandle_p);
1347              TMPMEMFREE(isCellCentered);
1348            return;            return;
1349          }          }
1350          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1337  void Finley_Mesh_saveVTK(const char * fi Line 1359  void Finley_Mesh_saveVTK(const char * fi
1359          {          {
1360            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1361            fclose(fileHandle_p);            fclose(fileHandle_p);
1362              TMPMEMFREE(isCellCentered);
1363            return;            return;
1364          }          }
1365          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1345  void Finley_Mesh_saveVTK(const char * fi Line 1368  void Finley_Mesh_saveVTK(const char * fi
1368          sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));          sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
1369          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
1370          fclose(fileHandle_p);          fclose(fileHandle_p);
1371            TMPMEMFREE(isCellCentered);
1372          return;          return;
1373        }        }
1374        if (isCellCentered[i_data])        if (isCellCentered[i_data])
# Line 1368  void Finley_Mesh_saveVTK(const char * fi Line 1392  void Finley_Mesh_saveVTK(const char * fi
1392      numPoints = mesh_p->Nodes->numNodes;      numPoints = mesh_p->Nodes->numNodes;
1393    }    }
1394    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
   Finley_ElementFile* elements=NULL;  
1395    switch(elementtype)    switch(elementtype)
1396    {    {
1397    case FINLEY_ELEMENTS:    case FINLEY_ELEMENTS:
# Line 1388  void Finley_Mesh_saveVTK(const char * fi Line 1411  void Finley_Mesh_saveVTK(const char * fi
1411    {    {
1412      Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");      Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1413      fclose(fileHandle_p);      fclose(fileHandle_p);
1414        TMPMEMFREE(isCellCentered);
1415      return;      return;
1416    }    }
1417    /* map finley element type to VTK element type */    /* map finley element type to VTK element type */
1418    numCells = elements->numElements;    numCells = elements->numElements;
   ElementTypeId TypeId;  
1419    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1420    {    {
1421      TypeId = elements->LinearReferenceElement->Type->TypeId;      TypeId = elements->LinearReferenceElement->Type->TypeId;
# Line 1501  void Finley_Mesh_saveVTK(const char * fi Line 1524  void Finley_Mesh_saveVTK(const char * fi
1524      sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);      sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
1525      Finley_setError(VALUE_ERROR,error_msg);      Finley_setError(VALUE_ERROR,error_msg);
1526      fclose(fileHandle_p);      fclose(fileHandle_p);
1527        TMPMEMFREE(isCellCentered);
1528      return;      return;
1529    }    }
1530    /* xml header */    /* xml header */
# Line 1523  void Finley_Mesh_saveVTK(const char * fi Line 1547  void Finley_Mesh_saveVTK(const char * fi
1547    * the reason for this if statement is explained in the long comment below    * the reason for this if statement is explained in the long comment below
1548    */    */
1549    nDim = mesh_p->Nodes->numDim;    nDim = mesh_p->Nodes->numDim;
1550    fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n",MAX(3,nDim));    fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n",MAX(3,nDim));
1551    /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate    /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
1552    * third dimension to handle 2D data (like a sheet of paper).  So, if    * third dimension to handle 2D data (like a sheet of paper).  So, if
1553    * nDim is 2, we have to append zeros to the array to get this third    * nDim is 2, we have to append zeros to the array to get this third
# Line 1559  void Finley_Mesh_saveVTK(const char * fi Line 1583  void Finley_Mesh_saveVTK(const char * fi
1583    
1584    /* write out the DataArray element for the connectivity */    /* write out the DataArray element for the connectivity */
1585    
1586    int NN = elements->ReferenceElement->Type->numNodes;    NN = elements->ReferenceElement->Type->numNodes;
1587    fprintf(fileHandle_p, "<Cells>\n");    fprintf(fileHandle_p, "<Cells>\n");
1588    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
1589    
# Line 1671  void Finley_Mesh_saveVTK(const char * fi Line 1695  void Finley_Mesh_saveVTK(const char * fi
1695            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1696            Finley_setError(VALUE_ERROR,error_msg);            Finley_setError(VALUE_ERROR,error_msg);
1697            fclose(fileHandle_p);            fclose(fileHandle_p);
1698              TMPMEMFREE(isCellCentered);
1699            return;            return;
1700          }          }
1701        }        }
# Line 1697  void Finley_Mesh_saveVTK(const char * fi Line 1722  void Finley_Mesh_saveVTK(const char * fi
1722            {            {
1723              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1724              fclose(fileHandle_p);              fclose(fileHandle_p);
1725                TMPMEMFREE(isCellCentered);
1726              return;              return;
1727            }            }
1728            nCompReqd = 3;            nCompReqd = 3;
# Line 1708  void Finley_Mesh_saveVTK(const char * fi Line 1734  void Finley_Mesh_saveVTK(const char * fi
1734            {            {
1735              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1736              fclose(fileHandle_p);              fclose(fileHandle_p);
1737                TMPMEMFREE(isCellCentered);
1738              return;              return;
1739            }            }
1740            nCompReqd = 9;            nCompReqd = 9;
1741          }          }
1742          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1743    
         double sampleAvg[nComp];  
1744          for (i=0; i<numCells; i++)          for (i=0; i<numCells; i++)
1745          {          {
1746            values = getSampleData(data_pp[i_data], i);            values = getSampleData(data_pp[i_data], i);
1747            /* averaging over the number of points in the sample */            /* averaging over the number of points in the sample */
1748            for (k=0; k<nComp; k++)            for (k=0; k<MIN(nComp,NCOMP_MAX); k++)
1749            {            {
1750              rtmp = 0.;              if (isExpanded(data_pp[i_data])) {
1751              for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];                 rtmp = 0.;
1752              sampleAvg[k] = rtmp/numPointsPerSample;                 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1753                   sampleAvg[k] = rtmp/numPointsPerSample;
1754                } else {
1755                   sampleAvg[k] = values[k];
1756                }
1757    
1758            }            }
1759            /* if the number of required components is more than the number            /* if the number of required components is more than the number
1760            * of actual components, pad with zeros            * of actual components, pad with zeros
# Line 1805  void Finley_Mesh_saveVTK(const char * fi Line 1836  void Finley_Mesh_saveVTK(const char * fi
1836            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1837            Finley_setError(VALUE_ERROR,error_msg);            Finley_setError(VALUE_ERROR,error_msg);
1838            fclose(fileHandle_p);            fclose(fileHandle_p);
1839              TMPMEMFREE(isCellCentered);
1840            return;            return;
1841          }          }
1842        }        }
# Line 1831  void Finley_Mesh_saveVTK(const char * fi Line 1863  void Finley_Mesh_saveVTK(const char * fi
1863            {            {
1864              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1865              fclose(fileHandle_p);              fclose(fileHandle_p);
1866                TMPMEMFREE(isCellCentered);
1867              return;              return;
1868            }            }
1869            nCompReqd = 3;            nCompReqd = 3;
# Line 1842  void Finley_Mesh_saveVTK(const char * fi Line 1875  void Finley_Mesh_saveVTK(const char * fi
1875            {            {
1876              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1877              fclose(fileHandle_p);              fclose(fileHandle_p);
1878                TMPMEMFREE(isCellCentered);
1879              return;              return;
1880            }            }
1881            nCompReqd = 9;            nCompReqd = 9;
1882          }          }
1883          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1884          /* write out the data */          /* write out the data */
1885          /* if the number of required components is more than the number          /* if the number of required components is more than the number
1886          * of actual components, pad with zeros          * of actual components, pad with zeros
1887          */          */
1888          bool_t do_write=TRUE;          do_write=TRUE;
1889          for (i=0; i<mesh_p->Nodes->numNodes; i++)          for (i=0; i<mesh_p->Nodes->numNodes; i++)
1890          {          {
1891            if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)            if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
# Line 1935  void Finley_Mesh_saveVTK(const char * fi Line 1969  void Finley_Mesh_saveVTK(const char * fi
1969    fprintf(fileHandle_p, "</VTKFile>\n");    fprintf(fileHandle_p, "</VTKFile>\n");
1970    /* close the file */    /* close the file */
1971    fclose(fileHandle_p);    fclose(fileHandle_p);
1972      TMPMEMFREE(isCellCentered);
1973    return;    return;
1974  }  }
1975  #endif  #endif

Legend:
Removed from v.818  
changed lines
  Added in v.1030

  ViewVC Help
Powered by ViewVC 1.1.26