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