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