/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Mesh_saveVTK.c

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

revision 2744 by caltinay, Mon Nov 16 23:45:55 2009 UTC revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 95  int nodeInQuadrant(const double *coords, Line 95  int nodeInQuadrant(const double *coords,
95  #define INSIDE_3D(_X_,_Y_,_Z_,_CX_,_CY_,_CZ_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) && INSIDE_1D(_Y_,_CY_,_R_) && INSIDE_1D(_Z_,_CZ_,_R_) )  #define INSIDE_3D(_X_,_Y_,_Z_,_CX_,_CY_,_CZ_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) && INSIDE_1D(_Y_,_CY_,_R_) && INSIDE_1D(_Z_,_CZ_,_R_) )
96    
97      int ret;      int ret;
98      if (type == Rec9) {      if ( (type == Rec9) || (type == Rec9Macro) ) {
99          if (q==0)          if (q==0)
100              ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.25, 0.25, 0.25);              ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.25, 0.25, 0.25);
101          else if (q==1)          else if (q==1)
# Line 106  int nodeInQuadrant(const double *coords, Line 106  int nodeInQuadrant(const double *coords,
106              ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.75, 0.75, 0.25);              ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.75, 0.75, 0.25);
107          else          else
108              ret = 0;              ret = 0;
109      } else if (type == Hex27) {      } else if ((type == Hex27) || (type == Hex27Macro) ){
110          if (q==0)          if (q==0)
111              ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],              ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
112                      0.25, 0.25, 0.25, 0.25);                      0.25, 0.25, 0.25, 0.25);
# Line 170  void Finley_Mesh_saveVTK(const char *fil Line 170  void Finley_Mesh_saveVTK(const char *fil
170      int mpi_size, i, j, l;      int mpi_size, i, j, l;
171      int cellType=0, nodeType=FINLEY_NODES, elementType=FINLEY_UNKNOWN;      int cellType=0, nodeType=FINLEY_NODES, elementType=FINLEY_UNKNOWN;
172      Finley_ElementFile *elements = NULL;      Finley_ElementFile *elements = NULL;
173      ElementTypeId typeId = NoType;      ElementTypeId typeId = NoRef;
174    
175      const char *vtkHeader = \      const char *vtkHeader = \
176        "<?xml version=\"1.0\"?>\n" \        "<?xml version=\"1.0\"?>\n" \
# Line 375  void Finley_Mesh_saveVTK(const char *fil Line 375  void Finley_Mesh_saveVTK(const char *fil
375              myFirstCell = Finley_ElementFile_getFirstElement(elements);              myFirstCell = Finley_ElementFile_getFirstElement(elements);
376              NN = elements->numNodes;              NN = elements->numNodes;
377              if (nodeType==FINLEY_REDUCED_NODES) {              if (nodeType==FINLEY_REDUCED_NODES) {
378                  typeId = elements->LinearReferenceElement->Type->TypeId;          typeId = elements->referenceElementSet->referenceElement->LinearType->TypeId;
379                  if (hasReducedElements) {                  if (hasReducedElements) {
380                      quadNodes_p=elements->LinearReferenceElementReducedOrder->QuadNodes;                      quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadNodes;
381                  } else {                  } else {
382                      quadNodes_p=elements->LinearReferenceElement->QuadNodes;                      quadNodes_p=elements->referenceElementSet->referenceElement->Parametrization->QuadNodes;
383                  }                  }
384              } else {              } else {
385                  typeId = elements->ReferenceElement->Type->TypeId;                  typeId = elements->referenceElementSet->referenceElement->Type->TypeId;
386                  if (hasReducedElements)                  if (hasReducedElements) {
387                      quadNodes_p=elements->ReferenceElementReducedOrder->QuadNodes;                      quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadNodes;
388                  else                  } else {
389                      quadNodes_p=elements->ReferenceElement->QuadNodes;                      quadNodes_p=elements->referenceElementSet->referenceElement->Parametrization->QuadNodes;
390                    }
391              }              }
392              switch (typeId) {              switch (typeId) {
393                  case Point1:                  case Point1:
# Line 425  void Finley_Mesh_saveVTK(const char *fil Line 426  void Finley_Mesh_saveVTK(const char *fil
426                      numVTKNodesPerElement = 4;                      numVTKNodesPerElement = 4;
427                  break;                  break;
428    
429                  case Rec9:                  case Rec9Macro:
430            case Rec9:
431                      numCellFactor = 4;                      numCellFactor = 4;
432                      cellType = VTK_QUAD;                      cellType = VTK_QUAD;
433                      numVTKNodesPerElement = 4;                      numVTKNodesPerElement = 4;
# Line 442  void Finley_Mesh_saveVTK(const char *fil Line 444  void Finley_Mesh_saveVTK(const char *fil
444                  break;                  break;
445    
446                  case Line3:                  case Line3:
447                    case Line3Macro:
448                  case Tri6Face:                  case Tri6Face:
449                  case Rec8Face:                  case Rec8Face:
450                  case Line3_Contact:                  case Line3_Contact:
# Line 452  void Finley_Mesh_saveVTK(const char *fil Line 455  void Finley_Mesh_saveVTK(const char *fil
455                  break;                  break;
456    
457                  case Tri6:                  case Tri6:
458            case Tri6Macro:    
459                  case Tet10Face:                  case Tet10Face:
460                  case Tri6_Contact:                  case Tri6_Contact:
461                  case Tet10Face_Contact:                  case Tet10Face_Contact:
# Line 467  void Finley_Mesh_saveVTK(const char *fil Line 471  void Finley_Mesh_saveVTK(const char *fil
471                      numVTKNodesPerElement = 8;                      numVTKNodesPerElement = 8;
472                  break;                  break;
473    
474            case Tet10Macro:
475                  case Tet10:                  case Tet10:
476                      cellType = VTK_QUADRATIC_TETRA;                      cellType = VTK_QUADRATIC_TETRA;
477                      numVTKNodesPerElement = 10;                      numVTKNodesPerElement = 10;
# Line 477  void Finley_Mesh_saveVTK(const char *fil Line 482  void Finley_Mesh_saveVTK(const char *fil
482                      numVTKNodesPerElement = 20;                      numVTKNodesPerElement = 20;
483                  break;                  break;
484    
485            case Hex27Macro:
486                  case Hex27:                  case Hex27:
487                      numCellFactor = 8;                      numCellFactor = 8;
488                      cellType = VTK_HEXAHEDRON;                      cellType = VTK_HEXAHEDRON;
# Line 484  void Finley_Mesh_saveVTK(const char *fil Line 490  void Finley_Mesh_saveVTK(const char *fil
490                  break;                  break;
491    
492                  default:                  default:
493                      sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->ReferenceElement->Type->Name);                      sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->referenceElementSet->referenceElement->Type->Name);
494                      Finley_setError(VALUE_ERROR, errorMsg);                      Finley_setError(VALUE_ERROR, errorMsg);
495              }              }
496          }          }
# Line 512  void Finley_Mesh_saveVTK(const char *fil Line 518  void Finley_Mesh_saveVTK(const char *fil
518      if (Finley_noError()) {      if (Finley_noError()) {
519          const index_t *nodeIndex;          const index_t *nodeIndex;
520          if (FINLEY_REDUCED_NODES == nodeType) {          if (FINLEY_REDUCED_NODES == nodeType) {
521              nodeIndex = elements->ReferenceElement->Type->linearNodes;              nodeIndex = elements->referenceElementSet->referenceElement->Type->linearNodes;
522          } else if (Rec9 == typeId) {          } else if ( (Rec9 == typeId) || (Rec9Macro == typeId) ) {
523              nodeIndex = VTK_REC9_INDEX;              nodeIndex = VTK_REC9_INDEX;
524          } else if (Hex20 == typeId) {          } else if (Hex20 == typeId) {
525              nodeIndex = VTK_HEX20_INDEX;              nodeIndex = VTK_HEX20_INDEX;
526          } else if (Hex27 == typeId) {          } else if ( (Hex27 == typeId) || (Hex27Macro == typeId) ){
527              nodeIndex = VTK_HEX27_INDEX;              nodeIndex = VTK_HEX27_INDEX;
528          } else if (numVTKNodesPerElement != NN) {          } else if (numVTKNodesPerElement  !=  elements->referenceElementSet->referenceElement->Type->numNodes) {
529              nodeIndex = elements->ReferenceElement->Type->geoNodes;              nodeIndex = elements->referenceElementSet->referenceElement->Type->relevantGeoNodes;
530          } else {          } else {
531              nodeIndex = NULL;              nodeIndex = NULL;
532          }          }
# Line 787  void Finley_Mesh_saveVTK(const char *fil Line 793  void Finley_Mesh_saveVTK(const char *fil
793    
794                              /* average over number of points in the sample */                              /* average over number of points in the sample */
795                              if (isExpanded(data_pp[dataIdx])) {                              if (isExpanded(data_pp[dataIdx])) {
796                                  dim_t hits=0, hits_old;                                  dim_t hits=0;
797                                  for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;                                  for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
798                                  for (j=0; j<numPointsPerSample; j++) {                                  for (j=0; j<numPointsPerSample; j++) {
                                     hits_old=hits;  
799                                      if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {                                      if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
800                                          hits++;                                          hits++;
801                                          for (k=0; k<nCompUsed; k++) {                                          for (k=0; k<nCompUsed; k++) {

Legend:
Removed from v.2744  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26