src/libengrid/tauwriter.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 #ifdef TAU_SUPPORT
  22
  23 #include "tauwriter.h"
  24 #include "guimainwindow.h"
  25 #include <netcdfcpp.h>
  26
  27
  28 TauWriter::TauWriter()
  29 {
  30   setFormat("Tau files(*.grid)");
  31   m_MenuText = "export to Tau format";
  32 }
  33
  34 void TauWriter::operate()
  35 {
  36   try {
  37     QFileInfo file_info(GuiMainWindow::pointer()->getFilename());
  38     readOutputFileName(file_info.completeBaseName() + ".grid");
  39     if (isValid()) {
  40       QString file_name = getFileName();
  41       NcFile *nc_file = new NcFile(file_name.toAscii(), NcFile::Replace);
  42       if (!nc_file->is_valid()) {
  43         EG_ERR_RETURN("unable to open NetCFD file for writing");
  44       }
  45
  46       // point coordinates
  47       size_t Np = m_Grid->GetNumberOfPoints();
  48       vector<double> x(Np),y(Np),z(Np);
  49       for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  50         vec3_t xv;
  51         m_Grid->GetPoint(id_node, xv.data());
  52         x[id_node] = xv[0];
  53         y[id_node] = xv[1];
  54         z[id_node] = xv[2];
  55       }
  56       NcDim *no_of_points = nc_file->add_dim("no_of_points", Np);
  57       NcVar *points_xc = nc_file->add_var("points_xc", ncDouble, no_of_points);
  58       NcVar *points_yc = nc_file->add_var("points_yc", ncDouble, no_of_points);
  59       NcVar *points_zc = nc_file->add_var("points_zc", ncDouble, no_of_points);
  60       points_xc->put(&x[0],Np);
  61       points_yc->put(&y[0],Np);
  62       points_zc->put(&z[0],Np);
  63
  64       // boundary faces
  65       size_t Nbe   = 0;;
  66       size_t Nquad = 0;
  67       size_t Ntri  = 0;
  68       for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  69         if (isSurface(id_cell, m_Grid)) {
  70           ++Nbe;
  71           if (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) {
  72             ++Ntri;
  73           } else if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
  74             ++Nquad;
  75           } else {
  76             EG_ERR_RETURN("unsupported boundary element type encountered");
  77           }
  78         }
  79       }
  80       NcDim *no_of_surfaceelements = nc_file->add_dim("no_of_surfaceelements", Nbe);
  81       NcVar *boundarymarker_of_surfaces = nc_file->add_var("boundarymarker_of_surfaces", ncInt, no_of_surfaceelements);
  82       vector<int> bm(Nbe,0), tri(3*Ntri), quad(4*Nquad);
  83       int i_bm = 0;
  84       QSet<int> bc_set = getAllBoundaryCodes(m_Grid);
  85       QVector<int> bcs(bc_set.size());
  86       qCopy(bc_set.begin(), bc_set.end(), bcs.begin());
  87       EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  88       if (Ntri) {
  89         NcDim *no_of_surfacetriangles = nc_file->add_dim("no_of_surfacetriangles", Ntri);
  90         NcDim *points_per_surfacetriangle = nc_file->add_dim("points_per_surfacetriangle", 3);
  91         NcVar *points_of_surfacetriangles = nc_file->add_var("points_of_surfacetriangles", ncInt, no_of_surfacetriangles, points_per_surfacetriangle);
  92         int i_tri = 0;
  93         for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  94           if (isSurface(id_cell, m_Grid)) {
  95             if (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) {
  96               for (int i_bc = 0; i_bc < bcs.size(); ++i_bc) {
  97                 if (bcs[i_bc] == cell_code->GetValue(id_cell)) {
  98                   bm[i_bm] = i_bc+1;
  99                   break;
 100                 }
 101               }
 102               vtkIdType N_pts, *pts;
 103               m_Grid->GetCellPoints(id_cell, N_pts, pts);
 104               tri[i_tri + 0] = pts[0];
 105               tri[i_tri + 1] = pts[1];
 106               tri[i_tri + 2] = pts[2];
 107               ++i_bm;
 108               i_tri += 3;
 109             }
 110           }
 111         }
 112         points_of_surfacetriangles->put(&tri[0],Ntri,3);
 113       }
 114       if (Nquad) {
 115         NcDim *no_of_surfacequadrilaterals     = nc_file->add_dim("no_of_surfacequadrilaterals",Nquad);
 116         NcDim *points_per_surfacequadrilateral = nc_file->add_dim("points_per_surfacequadrilateral",4);
 117         NcVar *points_of_surfacequadrilaterals = nc_file->add_var("points_of_surfacequadrilaterals",ncInt, no_of_surfacequadrilaterals, points_per_surfacequadrilateral);
 118         int i_quad = 0;
 119         for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 120           if (isSurface(id_cell, m_Grid)) {
 121             if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
 122               for (int i_bc = 0; i_bc < bcs.size(); ++i_bc) {
 123                 if (bcs[i_bc] == cell_code->GetValue(id_cell)) {
 124                   bm[i_bm] = i_bc+1;
 125                   break;
 126                 }
 127               }
 128               vtkIdType N_pts, *pts;
 129               m_Grid->GetCellPoints(id_cell, N_pts, pts);
 130               quad[i_quad + 0] = pts[0];
 131               quad[i_quad + 1] = pts[1];
 132               quad[i_quad + 2] = pts[2];
 133               quad[i_quad + 3] = pts[3];
 134               ++i_bm;
 135               i_quad += 4;
 136             }
 137           }
 138         }
 139         points_of_surfacequadrilaterals->put(&quad[0],Nquad,4);
 140       }
 141       boundarymarker_of_surfaces->put(&bm[0],Nbe);
 142       NcDim *no_of_markers = nc_file->add_dim("no_of_markers", bcs.size());
 143
 144
 145       int Ntet = 0;
 146       int Npri = 0;
 147       int Nhex = 0;
 148       int Npyr = 0;
 149       for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 150         if (m_Grid->GetCellType(id_cell) == VTK_TETRA)      ++Ntet;
 151         if (m_Grid->GetCellType(id_cell) == VTK_PYRAMID)    ++Npyr;
 152         if (m_Grid->GetCellType(id_cell) == VTK_WEDGE)      ++Npri;
 153         if (m_Grid->GetCellType(id_cell) == VTK_HEXAHEDRON) ++Nhex;
 154       }
 155
 156       vector<int> tet(Ntet*4),pyr(Npyr*5),pri(Npri*6),hex(Nhex*8);
 157       int i_tet = 0;
 158       int i_pyr = 0;
 159       int i_pri = 0;
 160       int i_hex = 0;
 161       if (Ntet) {
 162         NcDim *no_of_tetraeders = nc_file->add_dim("no_of_tetraeders",Ntet);
 163         NcDim *points_per_tetraeder = nc_file->add_dim("points_per_tetraeder",4);
 164         NcVar *points_of_tetraeders = nc_file->add_var("points_of_tetraeders",ncInt, no_of_tetraeders, points_per_tetraeder);
 165         for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 166           if (m_Grid->GetCellType(id_cell) == VTK_TETRA) {
 167             vtkIdType *pts, N_pts;
 168             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 169             tet[i_tet + 0] = pts[0];
 170             tet[i_tet + 1] = pts[1];
 171             tet[i_tet + 2] = pts[2];
 172             tet[i_tet + 3] = pts[3];
 173             i_tet += 4;
 174           }
 175         }
 176         points_of_tetraeders->put(&tet[0],Ntet,4);
 177       }
 178       if (Npyr) {
 179         NcDim *no_of_pyramids = nc_file->add_dim("no_of_pyramids",Npyr);
 180         NcDim *points_per_pyramid = nc_file->add_dim("points_per_pyramid",5);
 181         NcVar *points_of_pyramids = nc_file->add_var("points_of_pyramids",ncInt, no_of_pyramids, points_per_pyramid);
 182         for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 183           if (m_Grid->GetCellType(id_cell) == VTK_PYRAMID) {
 184             vtkIdType *pts, N_pts;
 185             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 186             pyr[i_pyr + 0] = pts[0];
 187             pyr[i_pyr + 1] = pts[1];
 188             pyr[i_pyr + 2] = pts[2];
 189             pyr[i_pyr + 3] = pts[3];
 190             pyr[i_pyr + 4] = pts[4];
 191             i_pyr += 5;
 192           }
 193         }
 194         points_of_pyramids->put(&pyr[0],Npyr,5);
 195       }
 196       if (Npri) {
 197         NcDim *no_of_prisms = nc_file->add_dim("no_of_prisms",Npri);
 198         NcDim *points_per_prism = nc_file->add_dim("points_per_prism",6);
 199         NcVar *points_of_prisms = nc_file->add_var("points_of_prisms",ncInt, no_of_prisms,  points_per_prism);
 200         for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 201           if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 202             vtkIdType *pts, N_pts;
 203             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 204             pri[i_pri + 0] = pts[0];
 205             pri[i_pri + 1] = pts[2];
 206             pri[i_pri + 2] = pts[1];
 207             pri[i_pri + 3] = pts[3];
 208             pri[i_pri + 4] = pts[5];
 209             pri[i_pri + 5] = pts[4];
 210             i_pri += 6;
 211           }
 212         }
 213         points_of_prisms->put(&pri[0],Npri,6);
 214       }
 215       if (Nhex) {
 216         NcDim *no_of_hexaeders = nc_file->add_dim("no_of_hexaeders",Nhex);
 217         NcDim *points_per_hexaeder = nc_file->add_dim("points_per_hexaeder",8);
 218         NcVar *points_of_hexaeders = nc_file->add_var("points_of_hexaeders",ncInt, no_of_hexaeders, points_per_hexaeder);
 219         for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 220           if (m_Grid->GetCellType(id_cell) == VTK_HEXAHEDRON) {
 221             vtkIdType *pts, N_pts;
 222             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 223             hex[i_hex + 0] = pts[4];
 224             hex[i_hex + 1] = pts[7];
 225             hex[i_hex + 2] = pts[6];
 226             hex[i_hex + 3] = pts[5];
 227             hex[i_hex + 4] = pts[0];
 228             hex[i_hex + 5] = pts[3];
 229             hex[i_hex + 6] = pts[2];
 230             hex[i_hex + 7] = pts[1];
 231             i_hex += 8;
 232           }
 233         }
 234         points_of_hexaeders->put(&hex[0],Nhex,8);
 235       }
 236
 237       delete nc_file;
 238     }
 239   } catch (Error err) {
 240     err.display();
 241   }
 242 }
 243
 244 #endif // TAU_SUPPORT