src/libengrid/cgnswriter.cpp

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