cgnswriter.cpp

   1 #include "cgnswriter.h"
   2
   3 CgnsWriter::CgnsWriter()
   4 {
   5   setFormat("CGNS files(*.cgns)");
   6 }
   7
   8 void CgnsWriter::writeGrid()
   9 {
  10   setAllCells();
  11   eg2cgns.fill(-1, cells.size());
  12
  13   // create the base node
  14   if (cg_base_write(fn, "Base", 3, 3, &B)) {
  15     EG_ERR_RETURN("error creating CGNS base node");
  16   }
  17
  18   int Nvcells = 0;
  19   foreach (vtkIdType id_cell, cells) {
  20     if (isVolume(id_cell, grid)) {
  21       ++Nvcells;
  22     }
  23   }
  24
  25   // create the zone node for the main grid
  26   int size[3];
  27   size[0] = nodes.size();
  28   size[1] = Nvcells;
  29   size[2] = 0;
  30   if (cg_zone_write(fn,B,"main_grid",size,Unstructured,&Z)) {
  31     EG_ERR_RETURN("error creating volume zone");
  32   }
  33
  34   // create a node for the grid coordinates
  35   int G;
  36   if (cg_grid_write(fn,B,Z,"GridCoordinates",&G)) {
  37     EG_ERR_RETURN("error creating GridCoordinates node");
  38   }
  39
  40   // create the arrays for the x,y, and z coordinates
  41   {
  42     int C;
  43     double coord_array[nodes.size()];
  44     vec3_t x;
  45     { // x coordinates
  46       for (int i = 0; i < nodes.size(); ++i) {
  47         if (nodes[i] != i) {
  48           EG_ERR_RETURN("unexpected hole in the point list");
  49         }
  50         grid->GetPoint(nodes[i], x.data());
  51         coord_array[i] = x[0];
  52       }
  53       if (cg_coord_write(fn, B, Z, RealDouble, "CoordinateX", coord_array, &C)) {
  54         EG_ERR_RETURN("error writing x-coordinates");
  55       }
  56     }
  57     { // y coordinates
  58       for (int i = 0; i < nodes.size(); ++i) {
  59         grid->GetPoint(nodes[i], x.data());
  60         coord_array[i] = x[1];
  61       }
  62       if (cg_coord_write(fn, B, Z, RealDouble, "CoordinateY", coord_array, &C)) {
  63         EG_ERR_RETURN("error writing y-coordinates");
  64       }
  65     }
  66     { // z coordinates
  67       for (int i = 0; i < nodes.size(); ++i) {
  68         grid->GetPoint(nodes[i], x.data());
  69         coord_array[i] = x[2];
  70       }
  71       if (cg_coord_write(fn, B, Z, RealDouble, "CoordinateZ", coord_array, &C)) {
  72         EG_ERR_RETURN("error writing z-coordinates");
  73       }
  74     }
  75   }
  76
  77   // count occurrences of different element types
  78   int ntet = 0;
  79   int npyr = 0;
  80   int npri = 0;
  81   int nhex = 0;
  82   int ntri = 0;
  83   int nqua = 0;
  84   foreach (vtkIdType id_cell, cells) {
  85     if      (grid->GetCellType(id_cell) == VTK_TETRA)      ++ntet;
  86     else if (grid->GetCellType(id_cell) == VTK_PYRAMID)    ++npyr;
  87     else if (grid->GetCellType(id_cell) == VTK_WEDGE)      ++npri;
  88     else if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) ++nhex;
  89     else if (grid->GetCellType(id_cell) == VTK_TRIANGLE)   ++ntri;
  90     else if (grid->GetCellType(id_cell) == VTK_QUAD)       ++nqua;
  91   };
  92
  93   // write element connectivity for tetras
  94   int start  = 0;
  95   int end    = 0;
  96   int i_cgns = 0;
  97   if (ntet) {
  98     int S;
  99     int elements[ntet*4];
 100     size_t j = 0;
 101     foreach (vtkIdType id_cell, cells) {
 102       if (grid->GetCellType(id_cell) == VTK_TETRA) {
 103         eg2cgns[id_cell] = i_cgns;
 104         ++i_cgns;
 105         vtkIdType Npts, *pts;
 106         grid->GetCellPoints(id_cell, Npts, pts);
 107         elements[j++] = pts[0]+1;
 108         elements[j++] = pts[1]+1;
 109         elements[j++] = pts[2]+1;
 110         elements[j++] = pts[3]+1;
 111       }
 112     }
 113     start = 1;
 114     end   = start+ntet-1;
 115     if (cg_section_write(fn, B, Z, "Tetras", TETRA_4, start, end, 0, elements, &S)) {
 116       EG_ERR_RETURN("error writing tetras");
 117     }
 118   }
 119
 120   // write element connectivity for pyramids
 121   if (npyr) {
 122     int S;
 123     int elements[npyr*5];
 124     size_t j = 0;
 125     foreach (vtkIdType id_cell, cells) {
 126       if (grid->GetCellType(id_cell) == VTK_PYRAMID) {
 127         eg2cgns[id_cell] = i_cgns;
 128         ++i_cgns;
 129         vtkIdType Npts, *pts;
 130         grid->GetCellPoints(id_cell, Npts, pts);
 131         elements[j++] = pts[0]+1;
 132         elements[j++] = pts[1]+1;
 133         elements[j++] = pts[2]+1;
 134         elements[j++] = pts[3]+1;
 135         elements[j++] = pts[4]+1;
 136       }
 137     }
 138     start = end+1;
 139     end   = start+npyr-1;
 140     if (cg_section_write(fn, B, Z, "Pyramids", PYRA_5, start, end, 0, elements, &S)) {
 141       EG_ERR_RETURN("error writing pyramids");
 142     }
 143   }
 144
 145   // write element connectivity for prisms
 146   if (npri) {
 147     int S;
 148     int elements[npri*6];
 149     size_t j = 0;
 150     foreach (vtkIdType id_cell, cells) {
 151       if (grid->GetCellType(id_cell) == VTK_WEDGE) {
 152         eg2cgns[id_cell] = i_cgns;
 153         ++i_cgns;
 154         vtkIdType Npts, *pts;
 155         grid->GetCellPoints(id_cell, Npts, pts);
 156         elements[j++] = pts[3]+1;
 157         elements[j++] = pts[4]+1;
 158         elements[j++] = pts[5]+1;
 159         elements[j++] = pts[0]+1;
 160         elements[j++] = pts[1]+1;
 161         elements[j++] = pts[2]+1;
 162       }
 163     }
 164     start = end+1;
 165     end   = start+npri-1;
 166     if (cg_section_write(fn, B, Z, "Prisms", PENTA_6, start, end, 0, elements, &S)) {
 167       EG_ERR_RETURN("error writing prisms");
 168     }
 169   }
 170
 171   // write element connectivity for hexas
 172   if (nhex) {
 173     int S;
 174     int elements[nhex*8];
 175     size_t j = 0;
 176     foreach (vtkIdType id_cell, cells) {
 177       if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) {
 178         eg2cgns[id_cell] = i_cgns;
 179         ++i_cgns;
 180         vtkIdType Npts, *pts;
 181         grid->GetCellPoints(id_cell, Npts, pts);
 182         elements[j++] = pts[0]+1;
 183         elements[j++] = pts[1]+1;
 184         elements[j++] = pts[3]+1;
 185         elements[j++] = pts[2]+1;
 186         elements[j++] = pts[4]+1;
 187         elements[j++] = pts[5]+1;
 188         elements[j++] = pts[7]+1;
 189         elements[j++] = pts[6]+1;
 190       }
 191     }
 192     start = end+1;
 193     end   = start+nhex-1;
 194     if (cg_section_write(fn, B, Z, "Hexes", HEXA_8, start, end, 0, elements, &S)) {
 195       EG_ERR_RETURN("error writing hexes");
 196     }
 197   }
 198
 199   // write element connectivity for triangles
 200   if (ntri) {
 201     int S;
 202     int elements[ntri*3];
 203     size_t j = 0;
 204     foreach (vtkIdType id_cell, cells) {
 205       if (grid->GetCellType(id_cell) == VTK_TRIANGLE) {
 206         eg2cgns[id_cell] = i_cgns;
 207         ++i_cgns;
 208         vtkIdType Npts, *pts;
 209         grid->GetCellPoints(id_cell, Npts, pts);
 210         elements[j++] = pts[2]+1;
 211         elements[j++] = pts[1]+1;
 212         elements[j++] = pts[0]+1;
 213       }
 214     }
 215     start = end+1;
 216     end   = start+ntri-1;
 217     if (cg_section_write(fn, B, Z, "Triangles", TRI_3, start, end, 0, elements, &S)) {
 218       EG_ERR_RETURN("error writing triangles");
 219     }
 220   }
 221
 222   // write element connectivity for quads
 223   if (nqua) {
 224     int S;
 225     int elements[nqua*4];
 226     size_t j = 0;
 227     foreach (vtkIdType id_cell, cells) {
 228       if (grid->GetCellType(id_cell) == VTK_QUAD) {
 229         eg2cgns[id_cell] = i_cgns;
 230         ++i_cgns;
 231         vtkIdType Npts, *pts;
 232         grid->GetCellPoints(id_cell, Npts, pts);
 233         elements[j++] = pts[3]+1;
 234         elements[j++] = pts[2]+1;
 235         elements[j++] = pts[1]+1;
 236         elements[j++] = pts[0]+1;
 237       }
 238     }
 239     start = end+1;
 240     end   = start+nqua-1;
 241     if (cg_section_write(fn, B, Z, "Quads", QUAD_4, start, end, 0, elements, &S)) {
 242       EG_ERR_RETURN("error writing quads");
 243     }
 244   }
 245 }
 246
 247 void CgnsWriter::writeBcs()
 248 {
 249   EG_VTKDCC(vtkIntArray, cell_code,   grid, "cell_code");
 250   QSet<int> bcs;
 251   {
 252     QVector<vtkIdType> faces;
 253     getAllSurfaceCells(faces, grid);
 254     foreach (vtkIdType id_face, faces) {
 255       bcs.insert(cell_code->GetValue(id_face));
 256     }
 257   }
 258   int BC_cgns;
 259   foreach (int bc, bcs) {
 260     QVector<vtkIdType> bc_faces;
 261     QSet<int> tmp_bcs;
 262     tmp_bcs.insert(bc);
 263     getSurfaceCells(tmp_bcs, bc_faces, grid);
 264     int data[bc_faces.size()];
 265     for (int i = 0; i < bc_faces.size(); ++i) {
 266       data[i] = eg2cgns[bc_faces[i]]+1;
 267     };
 268     if (cg_boco_write(fn, B, Z, getBC(bc).getName().toAscii().data(), BCTypeNull, ElementList, bc_faces.size(), data, &BC_cgns)) {
 269       cout << cg_get_error() << endl;
 270       EG_ERR_RETURN("error writing boundary condition");
 271     }
 272   }
 273 }
 274
 275 void CgnsWriter::operate()
 276 {
 277
 278 #ifndef CGNS_SUPPORT
 279   EG_ERR_RETURN("CGNS support has not been compiled");
 280 #endif
 281
 282   try {
 283     readOutputFileName();
 284     if (isValid()) {
 285       QString file_name = getFileName();
 286       if (cg_open(file_name.toAscii().data(), MODE_WRITE, &fn)) {
 287         EG_ERR_RETURN("error while opening CGNS file  for writing");
 288       }
 289       writeGrid();
 290       writeBcs();
 291       if (cg_close(fn)) {
 292         EG_ERR_RETURN("error while closing CGNS file");
 293       }
 294     }
 295   } catch (Error err) {
 296     err.display();
 297   }
 298 }
 299