src/libengrid/foamreader.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2013 enGits GmbH                                      +
   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 "foamreader.h"
  24 #include "guimainwindow.h"
  25
  26 FoamReader::FoamReader()
  27 {
  28   EG_TYPENAME;
  29 }
  30
  31 void FoamReader::operate()
  32 {
  33   try {
  34     readInputDirectory("Select OpenFOAM case directory");
  35     setCaseDir(getFileName());
  36     if (isValid()) {
  37       buildMaps();
  38       EG_VTKSP(vtkUnstructuredGrid, ug);
  39       int num_faces;
  40
  41       {
  42         readFile("constant/polyMesh/faces");
  43         QTextStream f(getBuffer());
  44         f >> num_faces;
  45       }
  46
  47       // count faces
  48       int num_triangles = 0;
  49       {
  50         readFile("constant/polyMesh/faces");
  51         QTextStream f(getBuffer());
  52         f >> num_faces;
  53         for (int i = 0; i < num_faces; ++i) {
  54           vtkIdType N_pts;
  55           f >> N_pts;
  56           QVector<vtkIdType> pts(N_pts);
  57           for (int j = 0; j < N_pts; ++j) {
  58             f >> pts[j];
  59           }
  60           if (i >= getFirstBoundaryFace()) {
  61             if (N_pts < 3) {
  62               EG_BUG;
  63             }
  64             num_triangles += (N_pts - 2);
  65           }
  66         }
  67       }
  68
  69       allocateGrid(ug, num_triangles, numSurfNodes());
  70       EG_VTKDCC(vtkIntArray, bc, ug, "cell_code");
  71       EG_VTKDCC(vtkIntArray, orgdir, ug, "cell_orgdir");
  72       EG_VTKDCC(vtkIntArray, voldir, ug, "cell_voldir");
  73       EG_VTKDCC(vtkIntArray, curdir, ug, "cell_curdir");
  74       QVector<vtkIdType> face2id(num_faces, -1);
  75       QVector<int> num_tris(num_faces, 0);
  76
  77       // insert points
  78       {
  79         readFile("constant/polyMesh/points");
  80         QTextStream f(getBuffer());
  81         int num_nodes;
  82         f >> num_nodes;
  83         for (int i = 0; i < num_nodes; ++i) {
  84           vec3_t x;
  85           f >> x[0] >> x[1] >> x[2];
  86           if (volToSurf(i) != -1) {
  87             ug->GetPoints()->SetPoint(volToSurf(i), x.data());
  88           }
  89         }
  90       }
  91
  92       qWarning() << ug->GetNumberOfPoints() << " nodes";
  93       qWarning() << num_faces - getFirstBoundaryFace() << " faces";
  94       qWarning() << num_triangles << " triangles";
  95       qWarning() << num_tris.size();
  96
  97       // insert faces
  98       {
  99         readFile("constant/polyMesh/faces");
 100         QTextStream f(getBuffer());
 101         f >> num_faces;
 102         for (int i = 0; i < num_faces; ++i) {
 103           vtkIdType N_pts;
 104           f >> N_pts;
 105           QVector<vtkIdType> pts(N_pts+2);
 106           for (int j = 0; j < N_pts; ++j) {
 107             f >> pts[j];
 108           }
 109           if (i >= getFirstBoundaryFace()) {
 110             for (int j = 0; j < N_pts; ++j) {
 111               pts[j] = volToSurf(pts[j]);
 112             }
 113             pts[N_pts]   = pts[0];
 114             pts[N_pts+1] = pts[1];
 115             QVector<vec3_t> x(N_pts+2);
 116             vec3_t xc(0,0,0);
 117             for (int j = 0; j < N_pts; ++j) {
 118               ug->GetPoint(pts[j], x[j].data());
 119               xc += x[j];
 120             }
 121             x[N_pts]   = x[0];
 122             x[N_pts+1] = x[1];
 123             xc *= 1.0/N_pts;
 124             vec3_t n(0,0,0);
 125             for (int j = 0; j < N_pts; ++j) {
 126               vec3_t u = x[j] - xc;
 127               vec3_t v = x[j+1] - xc;
 128               n += u.cross(v);
 129             }
 130             n.normalise();
 131             int j0 = 0;
 132             double scal0 = 1e99;
 133             for (int j = 1; j <= N_pts; ++j) {
 134               vec3_t u = x[j]   - x[j-1];
 135               vec3_t v = x[j+1] - x[j];
 136               u.normalise();
 137               v.normalise();
 138               vec3_t nj = u.cross(v);
 139               double scal = nj*n;
 140               if (scal < scal0) {
 141                 j0 = j % N_pts;
 142                 scal0 = scal;
 143               }
 144             }
 145             QList<int> jn;
 146             for (int j = 0; j < j0 - 1; ++j) {
 147               jn.append(j);
 148             }
 149             if (j0 > 0) {
 150               for (int j = j0 + 1; j < N_pts; ++j) {
 151                 jn.append(j);
 152               }
 153             } else {
 154               for (int j = j0 + 1; j < N_pts - 1; ++j) {
 155                 jn.append(j);
 156               }
 157             }
 158             foreach (int j, jn) {
 159               vtkIdType p[3];
 160               p[0] = pts[j];
 161               p[1] = pts[j+1];
 162               p[2] = pts[j0];
 163               vtkIdType id_cell = ug->InsertNextCell(VTK_TRIANGLE, 3, p);
 164               bc->SetValue(id_cell, 999);
 165               orgdir->SetValue(id_cell, 0);
 166               curdir->SetValue(id_cell, 0);
 167               voldir->SetValue(id_cell, 0);
 168               if (num_tris[i] == 0) {
 169                 face2id[i] = id_cell;
 170               }
 171               ++num_tris[i];
 172             }
 173           }
 174         }
 175       }
 176       GuiMainWindow::pointer()->clearBCs();
 177       {
 178         int N1 = 0;
 179         int N2 = 0;
 180         readFile("constant/polyMesh/boundary");
 181         QTextStream f(getBuffer());
 182         int num_patches;
 183         f >> num_patches;
 184         for (int i = 1; i <= num_patches; ++i) {
 185           QString name, dummy, type;
 186           f >> name >> dummy >> type;
 187           BoundaryCondition BC(name, type);
 188           GuiMainWindow::pointer()->addBC(i, BC);
 189           vtkIdType num_faces, start_face;
 190           f >> dummy >> num_faces;
 191           f >> dummy >> start_face;
 192           for (vtkIdType i_face = start_face; i_face < start_face + num_faces; ++i_face) {
 193             if (face2id[i_face] < 0) {
 194               EG_BUG;
 195             }
 196             ++N1;
 197             N2 += num_tris[i_face];
 198             for (int j = 0; j < num_tris[i_face]; ++j) {
 199               vtkIdType id_cell = face2id[i_face] + j;
 200               bc->SetValue(id_cell, i);
 201             }
 202           }
 203           qWarning() << i << ',' << N1 << ',' << N2;
 204         }
 205       }
 206       makeCopy(ug, m_Grid);
 207       createBasicFields(m_Grid, m_Grid->GetNumberOfCells(), m_Grid->GetNumberOfPoints());
 208       UpdateCellIndex(m_Grid);
 209     }
 210   } catch (Error err) {
 211     err.display();
 212   }
 213 }