src/cgnswriter.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "cgnswriter.h"
  24
  25 #include <QFileInfo>
  26 #include "guimainwindow.h"
  27
  28 CgnsWriter::CgnsWriter()
  29 {
  30   EG_TYPENAME;
  31 #ifdef CGNS_SUPPORT
  32   setFormat("CGNS files(*.cgns)");
  33 #endif
  34 }
  35
  36
  37 void CgnsWriter::writeGrid()
  38 {
  39 #ifdef CGNS_SUPPORT
  40   setAllCells();
  41   l2g_t cells = m_Part.getCells();
  42   l2g_t nodes = m_Part.getNodes();
  43   eg2cgns.fill(-1, cells.size());
  44
  45   // create the base node
  46   if (cg_base_write(fn, "Base", 3, 3, &B)) {
  47     EG_ERR_RETURN("error creating CGNS base node");
  48   }
  49
  50   int Nvcells = 0;
  51   foreach (vtkIdType id_cell, cells) {
  52     if (isVolume(id_cell, m_Grid)) {
  53       ++Nvcells;
  54     }
  55   }
  56
  57   // create the zone node for the main grid
  58   int size[3];
  59   size[0] = nodes.size();
  60   size[1] = Nvcells;
  61   size[2] = 0;
  62   if (cg_zone_write(fn,B,"main_grid",size,Unstructured,&Z)) {
  63     EG_ERR_RETURN("error creating volume zone");
  64   }
  65
  66   // create a node for the grid coordinates
  67   int G;
  68   if (cg_grid_write(fn,B,Z,"GridCoordinates",&G)) {
  69     EG_ERR_RETURN("error creating GridCoordinates node");
  70   }
  71
  72   // create the arrays for the x,y, and z coordinates
  73   {
  74     int C;
  75     double coord_array[nodes.size()];
  76     vec3_t x;
  77     { // x coordinates
  78       for (int i = 0; i < nodes.size(); ++i) {
  79         if (nodes[i] != i) {
  80           EG_ERR_RETURN("unexpected hole in the point list");
  81         }
  82         m_Grid->GetPoint(nodes[i], x.data());
  83         coord_array[i] = x[0];
  84       }
  85       if (cg_coord_write(fn, B, Z, RealDouble, "CoordinateX", coord_array, &C)) {
  86         EG_ERR_RETURN("error writing x-coordinates");
  87       }
  88     }
  89     { // y coordinates
  90       for (int i = 0; i < nodes.size(); ++i) {
  91         m_Grid->GetPoint(nodes[i], x.data());
  92         coord_array[i] = x[1];
  93       }
  94       if (cg_coord_write(fn, B, Z, RealDouble, "CoordinateY", coord_array, &C)) {
  95         EG_ERR_RETURN("error writing y-coordinates");
  96       }
  97     }
  98     { // z coordinates
  99       for (int i = 0; i < nodes.size(); ++i) {
 100         m_Grid->GetPoint(nodes[i], x.data());
 101         coord_array[i] = x[2];
 102       }
 103       if (cg_coord_write(fn, B, Z, RealDouble, "CoordinateZ", coord_array, &C)) {
 104         EG_ERR_RETURN("error writing z-coordinates");
 105       }
 106     }
 107   }
 108
 109   // count occurrences of different element types
 110   int ntet = 0;
 111   int npyr = 0;
 112   int npri = 0;
 113   int nhex = 0;
 114   int ntri = 0;
 115   int nqua = 0;
 116   foreach (vtkIdType id_cell, cells) {
 117     if      (m_Grid->GetCellType(id_cell) == VTK_TETRA)      ++ntet;
 118     else if (m_Grid->GetCellType(id_cell) == VTK_PYRAMID)    ++npyr;
 119     else if (m_Grid->GetCellType(id_cell) == VTK_WEDGE)      ++npri;
 120     else if (m_Grid->GetCellType(id_cell) == VTK_HEXAHEDRON) ++nhex;
 121     else if (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE)   ++ntri;
 122     else if (m_Grid->GetCellType(id_cell) == VTK_QUAD)       ++nqua;
 123   };
 124
 125   // write element connectivity for tetras
 126   int start  = 0;
 127   int end    = 0;
 128   int i_cgns = 0;
 129   if (ntet) {
 130     int S;
 131     int elements[ntet*4];
 132     size_t j = 0;
 133     foreach (vtkIdType id_cell, cells) {
 134       if (m_Grid->GetCellType(id_cell) == VTK_TETRA) {
 135         eg2cgns[id_cell] = i_cgns;
 136         ++i_cgns;
 137         vtkIdType Npts, *pts;
 138         m_Grid->GetCellPoints(id_cell, Npts, pts);
 139         elements[j++] = pts[0]+1;
 140         elements[j++] = pts[1]+1;
 141         elements[j++] = pts[2]+1;
 142         elements[j++] = pts[3]+1;
 143       }
 144     }
 145     start = 1;
 146     end   = start+ntet-1;
 147     if (cg_section_write(fn, B, Z, "Tetras", TETRA_4, start, end, 0, elements, &S)) {
 148       EG_ERR_RETURN("error writing tetras");
 149     }
 150   }
 151
 152   // write element connectivity for pyramids
 153   if (npyr) {
 154     int S;
 155     int elements[npyr*5];
 156     size_t j = 0;
 157     foreach (vtkIdType id_cell, cells) {
 158       if (m_Grid->GetCellType(id_cell) == VTK_PYRAMID) {
 159         eg2cgns[id_cell] = i_cgns;
 160         ++i_cgns;
 161         vtkIdType Npts, *pts;
 162         m_Grid->GetCellPoints(id_cell, Npts, pts);
 163         elements[j++] = pts[0]+1;
 164         elements[j++] = pts[1]+1;
 165         elements[j++] = pts[2]+1;
 166         elements[j++] = pts[3]+1;
 167         elements[j++] = pts[4]+1;
 168       }
 169     }
 170     start = end+1;
 171     end   = start+npyr-1;
 172     if (cg_section_write(fn, B, Z, "Pyramids", PYRA_5, start, end, 0, elements, &S)) {
 173       EG_ERR_RETURN("error writing pyramids");
 174     }
 175   }
 176
 177   // write element connectivity for prisms
 178   if (npri) {
 179     int S;
 180     int elements[npri*6];
 181     size_t j = 0;
 182     foreach (vtkIdType id_cell, cells) {
 183       if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 184         eg2cgns[id_cell] = i_cgns;
 185         ++i_cgns;
 186         vtkIdType Npts, *pts;
 187         m_Grid->GetCellPoints(id_cell, Npts, pts);
 188         elements[j++] = pts[3]+1;
 189         elements[j++] = pts[4]+1;
 190         elements[j++] = pts[5]+1;
 191         elements[j++] = pts[0]+1;
 192         elements[j++] = pts[1]+1;
 193         elements[j++] = pts[2]+1;
 194       }
 195     }
 196     start = end+1;
 197     end   = start+npri-1;
 198     if (cg_section_write(fn, B, Z, "Prisms", PENTA_6, start, end, 0, elements, &S)) {
 199       EG_ERR_RETURN("error writing prisms");
 200     }
 201   }
 202
 203   // write element connectivity for hexas
 204   if (nhex) {
 205     int S;
 206     int elements[nhex*8];
 207     size_t j = 0;
 208     foreach (vtkIdType id_cell, cells) {
 209       if (m_Grid->GetCellType(id_cell) == VTK_HEXAHEDRON) {
 210         eg2cgns[id_cell] = i_cgns;
 211         ++i_cgns;
 212         vtkIdType Npts, *pts;
 213         m_Grid->GetCellPoints(id_cell, Npts, pts);
 214         elements[j++] = pts[0]+1;
 215         elements[j++] = pts[1]+1;
 216         elements[j++] = pts[3]+1;
 217         elements[j++] = pts[2]+1;
 218         elements[j++] = pts[4]+1;
 219         elements[j++] = pts[5]+1;
 220         elements[j++] = pts[7]+1;
 221         elements[j++] = pts[6]+1;
 222       }
 223     }
 224     start = end+1;
 225     end   = start+nhex-1;
 226     if (cg_section_write(fn, B, Z, "Hexes", HEXA_8, start, end, 0, elements, &S)) {
 227       EG_ERR_RETURN("error writing hexes");
 228     }
 229   }
 230
 231   // write element connectivity for triangles
 232   if (ntri) {
 233     int S;
 234     int elements[ntri*3];
 235     size_t j = 0;
 236     foreach (vtkIdType id_cell, cells) {
 237       if (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) {
 238         eg2cgns[id_cell] = i_cgns;
 239         ++i_cgns;
 240         vtkIdType Npts, *pts;
 241         m_Grid->GetCellPoints(id_cell, Npts, pts);
 242         elements[j++] = pts[2]+1;
 243         elements[j++] = pts[1]+1;
 244         elements[j++] = pts[0]+1;
 245       }
 246     }
 247     start = end+1;
 248     end   = start+ntri-1;
 249     if (cg_section_write(fn, B, Z, "Triangles", TRI_3, start, end, 0, elements, &S)) {
 250       EG_ERR_RETURN("error writing triangles");
 251     }
 252   }
 253
 254   // write element connectivity for quads
 255   if (nqua) {
 256     int S;
 257     int elements[nqua*4];
 258     size_t j = 0;
 259     foreach (vtkIdType id_cell, cells) {
 260       if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
 261         eg2cgns[id_cell] = i_cgns;
 262         ++i_cgns;
 263         vtkIdType Npts, *pts;
 264         m_Grid->GetCellPoints(id_cell, Npts, pts);
 265         elements[j++] = pts[3]+1;
 266         elements[j++] = pts[2]+1;
 267         elements[j++] = pts[1]+1;
 268         elements[j++] = pts[0]+1;
 269       }
 270     }
 271     start = end+1;
 272     end   = start+nqua-1;
 273     if (cg_section_write(fn, B, Z, "Quads", QUAD_4, start, end, 0, elements, &S)) {
 274       EG_ERR_RETURN("error writing quads");
 275     }
 276   }
 277 #endif
 278 }
 279
 280
 281 void CgnsWriter::writeBcs()
 282 {
 283 #ifdef CGNS_SUPPORT
 284   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 285   QSet<int> bcs;
 286   {
 287     QVector<vtkIdType> faces;
 288     getAllSurfaceCells(faces, m_Grid);
 289     foreach (vtkIdType id_face, faces) {
 290       bcs.insert(cell_code->GetValue(id_face));
 291     }
 292   }
 293   int BC_cgns;
 294   foreach (int bc, bcs) {
 295     QVector<vtkIdType> bc_faces;
 296     QSet<int> tmp_bcs;
 297     tmp_bcs.insert(bc);
 298     getSurfaceCells(tmp_bcs, bc_faces, m_Grid);
 299     int data[bc_faces.size()];
 300     for (int i = 0; i < bc_faces.size(); ++i) {
 301       data[i] = eg2cgns[bc_faces[i]]+1;
 302     };
 303     if (cg_boco_write(fn, B, Z, qPrintable(getBC(bc).getName()), BCTypeNull, ElementList, bc_faces.size(), data, &BC_cgns)) {
 304       cout << cg_get_error() << endl;
 305       EG_ERR_RETURN("error writing boundary condition");
 306     }
 307   }
 308 #endif
 309 }
 310
 311
 312 void CgnsWriter::operate()
 313 {
 314 #ifdef CGNS_SUPPORT
 315   try {
 316     QFileInfo file_info(GuiMainWindow::pointer()->getFilename());
 317     readOutputFileName(file_info.completeBaseName() + ".cgns");
 318     if (isValid()) {
 319       QString file_name = getFileName();
 320       if (cg_open(qPrintable(file_name), MODE_WRITE, &fn)) {
 321         EG_ERR_RETURN("error while opening CGNS file for writing");
 322       }
 323       writeGrid();
 324       writeBcs();
 325       if (cg_close(fn)) {
 326         EG_ERR_RETURN("error while closing CGNS file");
 327       }
 328     }
 329   } catch (Error err) {
 330     err.display();
 331   }
 332 #else
 333   EG_ERR_RETURN("CGNS support has not been compiled");
 334 #endif
 335 }
 336