src/libengrid/gmshreader.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 "gmshreader.h"
  22 #include "correctsurfaceorientation.h"
  23
  24 #include <QFileInfo>
  25 #include "guimainwindow.h"
  26
  27 void GmshReader::readAscii1(vtkUnstructuredGrid *m_Grid)
  28 {
  29   vtkIdType Nnodes, Ncells;
  30   QFile file(getFileName());
  31   file.open(QIODevice::ReadOnly | QIODevice::Text);
  32   QTextStream f(&file);
  33   QString word;
  34   f >> word;
  35   if (word != "$NOD") EG_ERR_RETURN("$NOD expected");
  36   f >> Nnodes;
  37   EG_VTKSP(vtkUnstructuredGrid, ug);
  38   QVector<vtkIdType> idxmap(Nnodes + 1);
  39   QVector<vec3_t> x(Nnodes);
  40   for (vtkIdType i = 0; i < Nnodes; ++i) {
  41     int ir;
  42     f >> ir >> x[i][0] >> x[i][1] >> x[i][2];
  43     idxmap[ir] = i;
  44   }
  45   f >> word;
  46   if (word != "$ENDNOD") EG_ERR_RETURN("$ENDNOD expected");
  47   f >> word;
  48   if (word != "$ELM") EG_ERR_RETURN("$ELM expected");
  49   f >> Ncells;
  50   allocateGrid(ug, Ncells, Nnodes, false);
  51   for (vtkIdType i = 0; i < Nnodes; ++i) {
  52     ug->GetPoints()->SetPoint(i, x[i].data());
  53   }
  54   EG_VTKSP(vtkIntArray, cell_code);
  55   cell_code->SetName("cell_code");
  56   cell_code->SetNumberOfValues(Ncells);
  57   ug->GetCellData()->AddArray(cell_code);
  58   for (vtkIdType i = 0; i < Ncells; ++i) {
  59     f >> word;
  60     int elm_type, reg_phys;
  61     f >> elm_type;
  62     f >> reg_phys;
  63     f >> word;
  64     f >> word;
  65     cell_code->SetValue(i, reg_phys);
  66     vtkIdType pts[8];
  67     size_t node;
  68     if        (elm_type == 2) { // triangle
  69       for (vtkIdType j = 0; j < 3; ++j) {
  70         f >> node;
  71         pts[j] = idxmap[node];
  72       }
  73       ug->InsertNextCell(VTK_TRIANGLE,3,pts);
  74     } else if (elm_type == 3) { // quad
  75       for (vtkIdType j = 0; j < 4; ++j) {
  76         f >> node;
  77         pts[j] = idxmap[node];
  78       }
  79       ug->InsertNextCell(VTK_QUAD,4,pts);
  80     } else if (elm_type == 4) { // tetrahedron
  81       for (vtkIdType j = 0; j < 4; ++j) {
  82         f >> node;
  83         pts[j] = idxmap[node];
  84       }
  85       vtkIdType h = pts[0];
  86       pts[0] = pts[1];
  87       pts[1] = h;
  88       ug->InsertNextCell(VTK_TETRA,4,pts);
  89     } else if (elm_type == 5) { // hexhedron
  90       for (vtkIdType j = 0; j < 8; ++j) {
  91         f >> node;
  92         pts[j] = idxmap[node];
  93       }
  94       ug->InsertNextCell(VTK_HEXAHEDRON,8,pts);
  95     } else if (elm_type == 6) { // prism/wedge
  96       for (vtkIdType j = 0; j < 3; ++j) {
  97         f >> node;
  98         pts[j+3] = idxmap[node];
  99       }
 100       for (vtkIdType j = 0; j < 3; ++j) {
 101         f >> node;
 102         pts[j] = idxmap[node];
 103       }
 104       ug->InsertNextCell(VTK_WEDGE,6,pts);
 105     } else if (elm_type == 7) { // pyramid
 106       for (vtkIdType j = 0; j < 5; ++j) {
 107         f >> node;
 108         pts[j] = idxmap[node];
 109       }
 110       ug->InsertNextCell(VTK_PYRAMID,5,pts);
 111     }
 112   }
 113   ug->GetCellData()->AddArray(cell_code);
 114   m_Grid->DeepCopy(ug);
 115 }
 116
 117 void GmshReader::readAscii2(vtkUnstructuredGrid *m_Grid)
 118 {
 119   vtkIdType Nnodes, Ncells;
 120   QFile file(getFileName());
 121   file.open(QIODevice::ReadOnly | QIODevice::Text);
 122   QTextStream f(&file);
 123   QString word;
 124   f >> word;
 125   if (word != "$MeshFormat") EG_ERR_RETURN("$MeshFormat expected");
 126   f >> word;
 127   f >> word;
 128   f >> word;
 129   f >> word;
 130   if (word != "$EndMeshFormat") EG_ERR_RETURN("$EndMeshFormat expected");
 131   f >> word;
 132   if (word != "$Nodes") EG_ERR_RETURN("$Nodes expected");
 133   f >> Nnodes;
 134   EG_VTKSP(vtkUnstructuredGrid, ug);
 135   QVector<vtkIdType> idxmap(Nnodes + 1);
 136   QVector<vec3_t> x(Nnodes);
 137   for (vtkIdType i = 0; i < Nnodes; ++i) {
 138     int ir;
 139     f >> ir >> x[i][0] >> x[i][1] >> x[i][2];
 140     idxmap[ir] = i;
 141   }
 142   f >> word;
 143   if (word != "$EndNodes") EG_ERR_RETURN("$EndNotes expected");
 144   f >> word;
 145   if (word != "$Elements") EG_ERR_RETURN("$Elements expected");
 146   f >> Ncells;
 147   allocateGrid(ug, Ncells, Nnodes, false);
 148   for (vtkIdType i = 0; i < Nnodes; ++i) {
 149     ug->GetPoints()->SetPoint(i, x[i].data());
 150   }
 151   EG_VTKSP(vtkIntArray, cell_code);
 152   cell_code->SetName("cell_code");
 153   cell_code->SetNumberOfValues(Ncells);
 154   ug->GetCellData()->AddArray(cell_code);
 155
 156   for (vtkIdType i = 0; i < Ncells; ++i) {
 157     f >> word;
 158     int elm_type, Ntags, bc;
 159     f >> elm_type;
 160     f >> Ntags;
 161     if (Ntags == 0) {
 162       bc = 1;
 163     } else {
 164       f >> bc;
 165       if (bc <= 0) {
 166         bc = 99;
 167       }
 168       for (int j = 1; j < Ntags; ++j) f >> word;
 169     }
 170     cell_code->SetValue(i, bc);
 171     vtkIdType pts[8];
 172     size_t node;
 173     if        (elm_type == 2) { // triangle
 174       for (vtkIdType j = 0; j < 3; ++j) {
 175         f >> node;
 176         pts[j] = idxmap[node];
 177       }
 178       ug->InsertNextCell(VTK_TRIANGLE,3,pts);
 179     } else if (elm_type == 3) { // quad
 180       for (vtkIdType j = 0; j < 4; ++j) {
 181         f >> node;
 182         pts[j] = idxmap[node];
 183       }
 184       ug->InsertNextCell(VTK_QUAD,4,pts);
 185     } else if (elm_type == 4) { // tetrahedron
 186       for (vtkIdType j = 0; j < 4; ++j) {
 187         f >> node;
 188         pts[j] = idxmap[node];
 189       }
 190       ug->InsertNextCell(VTK_TETRA,4,pts);
 191     } else if (elm_type == 5) { // hexhedron
 192       for (vtkIdType j = 0; j < 8; ++j) {
 193         f >> node;
 194         pts[j] = idxmap[node];
 195       }
 196       ug->InsertNextCell(VTK_HEXAHEDRON,8,pts);
 197     } else if (elm_type == 6) { // prism/wedge
 198       for (vtkIdType j = 0; j < 3; ++j) {
 199         f >> node;
 200         pts[j+3] = idxmap[node];
 201       }
 202       for (vtkIdType j = 0; j < 3; ++j) {
 203         f >> node;
 204         pts[j] = idxmap[node];
 205       }
 206       ug->InsertNextCell(VTK_WEDGE,6,pts);
 207     } else if (elm_type == 7) { // pyramid
 208       for (vtkIdType j = 0; j < 5; ++j) {
 209         f >> node;
 210         pts[j] = idxmap[node];
 211       }
 212       ug->InsertNextCell(VTK_PYRAMID,5,pts);
 213     }
 214   }
 215   ug->GetCellData()->AddArray(cell_code);
 216   m_Grid->DeepCopy(ug);
 217 }
 218
 219 void GmshReader::operate()
 220 {
 221   try {
 222     QFileInfo file_info(GuiMainWindow::pointer()->getFilename());
 223     readInputFileName(file_info.completeBaseName() + ".msh");
 224     if (isValid()) {
 225       if (format == ascii1) {
 226         readAscii1(m_Grid);
 227       }
 228       if (format == ascii2) {
 229         readAscii2(m_Grid);
 230       }
 231       createBasicFields(m_Grid, m_Grid->GetNumberOfCells(), m_Grid->GetNumberOfPoints());
 232       updateCellIndex(m_Grid);
 233       CorrectSurfaceOrientation corr_surf;
 234       corr_surf.setGrid(m_Grid);
 235       corr_surf();
 236     }
 237   } catch (Error err) {
 238     err.display();
 239   }
 240 }