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