src/libengrid/simplefoamwriter.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
  22 #include "simplefoamwriter.h"
  23 #include "guimainwindow.h"
  24
  25 #include <QFileInfo>
  26 #include <QDir>
  27
  28 bool SimpleFoamWriter::face_t::operator<(const face_t &F) const
  29 {
  30   bool less = false;
  31   if (bc < F.bc) {
  32     less = true;
  33   } else if (bc == F.bc) {
  34     if (owner < F.owner) {
  35       less = true;
  36     } else if (owner == F.owner) {
  37       if (neighbour < F.neighbour) {
  38         less = true;
  39       }
  40     }
  41   }
  42   return less;
  43 }
  44
  45 vec3_t SimpleFoamWriter::face_t::normal(vtkUnstructuredGrid *m_Grid)
  46 {
  47   if (node.size() < 3) EG_BUG;
  48   vec3_t xc(0,0,0);
  49   QVector<vec3_t> x(node.size());
  50   for (int i = 0; i < node.size(); ++i) {
  51     m_Grid->GetPoint(node[i],x[i].data());
  52     xc += x[i];
  53   }
  54   xc *= 1.0/node.size();
  55   vec3_t n(0,0,0);
  56   for (int i = 0; i < node.size()-1; ++i) {
  57     vec3_t a = x[i]   - xc;
  58     vec3_t b = x[i+1] - xc;
  59     n += 0.5*(a.cross(b));
  60   }
  61   vec3_t a = x[node.size()-1] - xc;
  62   vec3_t b = x[0]             - xc;
  63   n += 0.5*(a.cross(b));
  64   return n;
  65 }
  66
  67
  68 SimpleFoamWriter::SimpleFoamWriter()
  69 {
  70   setFormat("Foam boundary files(boundary)");
  71   setExtension("");
  72 }
  73
  74 vtkIdType SimpleFoamWriter::getNeigh(int i_cells, int i_neigh)
  75 {
  76   l2g_t cells = getPartCells();
  77   l2l_t c2c   = getPartC2C();
  78   int n = c2c[i_cells][i_neigh];
  79   if (n >= 0) return cells[n];
  80   EG_ERR_RETURN("The grid is not suitable for OpenFOAM export (e.g. missing volume mesh).");
  81   return -1;
  82 }
  83
  84 void SimpleFoamWriter::addFace(face_t F)
  85 {
  86   if (isVolume(F.neighbour,m_Grid)) {
  87     if (F.neighbour > F.owner) {
  88       F.bc = 0;
  89       m_LFaces.append(F);
  90     }
  91   } else {
  92     F.bc = m_BC->GetValue(F.neighbour);
  93     F.neighbour = -1;
  94     m_LFaces.append(F);
  95   }
  96 }
  97
  98 void SimpleFoamWriter::createFaces()
  99 {
 100   l2g_t cells = getPartCells();
 101
 102   m_LFaces.clear();
 103   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 104   m_BC = cell_code;
 105   m_Eg2Of.fill(-1,cells.size());
 106   int Nvol = 0;
 107   foreach(int i_cells, cells) {
 108     vtkIdType *pts;
 109     vtkIdType  Npts;
 110     m_Grid->GetCellPoints(cells[i_cells], Npts, pts);
 111     vtkIdType type_cell = m_Grid->GetCellType(cells[i_cells]);
 112     vtkIdType id_cell = cells[i_cells];
 113
 114     // tetras
 115     //
 116     if (type_cell == VTK_TETRA) {
 117       m_Eg2Of[id_cell] = Nvol++;
 118       {
 119         face_t F(3,id_cell,getNeigh(i_cells,0));
 120         F.node[0] = pts[2]; F.node[1] = pts[1]; F.node[2] = pts[0];
 121         addFace(F);
 122       }
 123       {
 124         face_t F(3,id_cell,getNeigh(i_cells,1));
 125         F.node[0] = pts[0]; F.node[1] = pts[1]; F.node[2] = pts[3];
 126         addFace(F);
 127       }
 128       {
 129         face_t F(3,id_cell,getNeigh(i_cells,2));
 130         F.node[0] = pts[0]; F.node[1] = pts[3]; F.node[2] = pts[2];
 131         addFace(F);
 132       }
 133       {
 134         face_t F(3,id_cell,getNeigh(i_cells,3));
 135         F.node[0] = pts[1]; F.node[1] = pts[2]; F.node[2] = pts[3];
 136         addFace(F);
 137       }
 138     }
 139
 140     // prisms
 141     //
 142     if (type_cell == VTK_WEDGE) {
 143       m_Eg2Of[id_cell] = Nvol++;
 144       {
 145         face_t F(3,id_cell,getNeigh(i_cells,0));
 146         F.node[0] = pts[0]; F.node[1] = pts[1]; F.node[2] = pts[2];
 147         addFace(F);
 148       }
 149       {
 150         face_t F(3,id_cell,getNeigh(i_cells,1));
 151         F.node[0] = pts[3]; F.node[1] = pts[5]; F.node[2] = pts[4];
 152         addFace(F);
 153       }
 154       {
 155         face_t F(4,id_cell,getNeigh(i_cells,2));
 156         F.node[0] = pts[3]; F.node[1] = pts[4]; F.node[2] = pts[1]; F.node[3] = pts[0];
 157         addFace(F);
 158       }
 159       {
 160         face_t F(4,id_cell,getNeigh(i_cells,3));
 161         F.node[0] = pts[1]; F.node[1] = pts[4]; F.node[2] = pts[5]; F.node[3] = pts[2];
 162         addFace(F);
 163       }
 164       {
 165         face_t F(4,id_cell,getNeigh(i_cells,4));
 166         F.node[0] = pts[0]; F.node[1] = pts[2]; F.node[2] = pts[5]; F.node[3] = pts[3];
 167         addFace(F);
 168       }
 169     }
 170
 171     // hexes
 172     //
 173     if (type_cell == VTK_HEXAHEDRON) {
 174       m_Eg2Of[id_cell] = Nvol++;
 175       {
 176         face_t F(4,id_cell,getNeigh(i_cells,0),0);
 177         F.node[0] = pts[3]; F.node[1] = pts[2]; F.node[2] = pts[1]; F.node[3] = pts[0];
 178         addFace(F);
 179       }
 180       {
 181         face_t F(4,id_cell,getNeigh(i_cells,1),0);
 182         F.node[0] = pts[4]; F.node[1] = pts[5]; F.node[2] = pts[6]; F.node[3] = pts[7];
 183         addFace(F);
 184       }
 185       {
 186         face_t F(4,id_cell,getNeigh(i_cells,2),0);
 187         F.node[0] = pts[0]; F.node[1] = pts[1]; F.node[2] = pts[5]; F.node[3] = pts[4];
 188         addFace(F);
 189       }
 190       {
 191         face_t F(4,id_cell,getNeigh(i_cells,3),0);
 192         F.node[0] = pts[3]; F.node[1] = pts[7]; F.node[2] = pts[6]; F.node[3] = pts[2];
 193         addFace(F);
 194       }
 195       {
 196         face_t F(4,id_cell,getNeigh(i_cells,4),0);
 197         F.node[0] = pts[0]; F.node[1] = pts[4]; F.node[2] = pts[7]; F.node[3] = pts[3];
 198         addFace(F);
 199       }
 200       {
 201         face_t F(4,id_cell,getNeigh(i_cells,5),0);
 202         F.node[0] = pts[1]; F.node[1] = pts[2]; F.node[2] = pts[6]; F.node[3] = pts[5];
 203         addFace(F);
 204       }
 205     }
 206   }
 207
 208   m_Faces.resize(m_LFaces.size());
 209   qCopy(m_LFaces.begin(), m_LFaces.end(), m_Faces.begin());
 210   qSort(m_Faces);
 211
 212 }
 213
 214
 215 void SimpleFoamWriter::writePoints()
 216 {
 217   l2g_t nodes = getPartNodes();
 218
 219   QString filename = m_Path + "points";
 220   QFile file(filename);
 221   file.open(QIODevice::WriteOnly);
 222   QTextStream f(&file);
 223   f << "/*--------------------------------*- C++ -*----------------------------------*\\\n";
 224   f << "| =========                 |                                                 |\n";
 225   f << "| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n";
 226   f << "|  \\    /   O peration     | Version:  1.5                                   |\n";
 227   f << "|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |\n";
 228   f << "|    \\/     M anipulation  |                                                 |\n";
 229   f << "\\*---------------------------------------------------------------------------*/\n\n";
 230   f << "FoamFile\n";
 231   f << "{\n";
 232   f << "    version     2.0;\n";
 233   f << "    format      ascii;\n";
 234   f << "    class       vectorField;\n";
 235   f << "    location    \"constant/polyMesh\";\n";
 236   f << "    object      points;\n";
 237   f << "}\n\n";
 238   f << "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n\n";
 239   f << nodes.size() << "\n(\n";
 240   foreach(int i_nodes, nodes) {
 241     vtkIdType id_node = nodes[i_nodes];
 242     vec3_t x;
 243     m_Grid->GetPoint(id_node,x.data());
 244     f.setRealNumberPrecision(16);
 245     f << "(" << x[0] << " " << x[1] << " " << x[2] << ")\n";
 246   }
 247   f << ")\n\n";
 248   f << "// ************************************************************************* //\n\n\n";
 249 }
 250
 251 void SimpleFoamWriter::writeFaces()
 252 {
 253   QString filename = m_Path + "faces";
 254   QFile file(filename);
 255   file.open(QIODevice::WriteOnly);
 256   QTextStream f(&file);
 257   f << "/*--------------------------------*- C++ -*----------------------------------*\\\n";
 258   f << "| =========                 |                                                 |\n";
 259   f << "| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n";
 260   f << "|  \\    /   O peration     | Version:  1.5                                   |\n";
 261   f << "|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |\n";
 262   f << "|    \\/     M anipulation  |                                                 |\n";
 263   f << "\\*---------------------------------------------------------------------------*/\n\n";
 264   f << "FoamFile\n";
 265   f << "{\n";
 266   f << "    version     2.0;\n";
 267   f << "    format      ascii;\n";
 268   f << "    class       faceList;\n";
 269   f << "    location    \"constant/polyMesh\";\n";
 270   f << "    object      faces;\n";
 271   f << "}\n\n";
 272   f << "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n\n";
 273   f << m_Faces.size() << "\n(\n";
 274   foreach (face_t F, m_Faces) {
 275     f << F.node.size() << "(";
 276     for (int i = 0; i < F.node.size(); ++i) {
 277       f << F.node[i];
 278       if (i == F.node.size()-1) {
 279         f << ")\n";
 280       } else {
 281         f << " ";
 282       }
 283     }
 284   }
 285   f << ")\n\n";
 286   f << "// ************************************************************************* //\n\n\n";
 287 }
 288
 289 void SimpleFoamWriter::writeOwner()
 290 {
 291   QString filename = m_Path + "owner";
 292   QFile file(filename);
 293   file.open(QIODevice::WriteOnly);
 294   QTextStream f(&file);
 295   f << "/*--------------------------------*- C++ -*----------------------------------*\\\n";
 296   f << "| =========                 |                                                 |\n";
 297   f << "| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n";
 298   f << "|  \\    /   O peration     | Version:  1.5                                   |\n";
 299   f << "|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |\n";
 300   f << "|    \\/     M anipulation  |                                                 |\n";
 301   f << "\\*---------------------------------------------------------------------------*/\n\n";
 302   f << "FoamFile\n";
 303   f << "{\n";
 304   f << "    version     2.0;\n";
 305   f << "    format      ascii;\n";
 306   f << "    class       labelList;\n";
 307   f << "    location    \"constant/polyMesh\";\n";
 308   f << "    object      owner;\n";
 309   f << "}\n\n";
 310   f << "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n\n";
 311   f << m_Faces.size() << "\n(\n";
 312   foreach (face_t F, m_Faces) {
 313     f << m_Eg2Of[F.owner] << "\n";
 314   }
 315   f << ")\n\n";
 316   f << "// ************************************************************************* //\n\n\n";
 317 }
 318
 319 void SimpleFoamWriter::writeNeighbour()
 320 {
 321   QString filename = m_Path + "neighbour";
 322   QFile file(filename);
 323   file.open(QIODevice::WriteOnly);
 324   QTextStream f(&file);
 325   f << "/*--------------------------------*- C++ -*----------------------------------*\\\n";
 326   f << "| =========                 |                                                 |\n";
 327   f << "| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n";
 328   f << "|  \\    /   O peration     | Version:  1.5                                   |\n";
 329   f << "|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |\n";
 330   f << "|    \\/     M anipulation  |                                                 |\n";
 331   f << "\\*---------------------------------------------------------------------------*/\n\n";
 332   f << "FoamFile\n";
 333   f << "{\n";
 334   f << "    version     2.0;\n";
 335   f << "    format      ascii;\n";
 336   f << "    class       labelList;\n";
 337   f << "    location    \"constant/polyMesh\";\n";
 338   f << "    object      neighbour;\n";
 339   f << "}\n\n";
 340   f << "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n\n";
 341   f << m_Faces.size() << "\n(\n";
 342   foreach (face_t F, m_Faces) {
 343     if (F.neighbour == -1) {
 344       f << "-1\n";
 345     } else {
 346       f << m_Eg2Of[F.neighbour] << "\n";
 347     }
 348   }
 349   f << ")\n\n";
 350   f << "// ************************************************************************* //\n\n\n";
 351 }
 352
 353 void SimpleFoamWriter::writeBoundary(int faces_offset)
 354 {
 355   QString filename = m_Path + "boundary";
 356   QFile file(filename);
 357   file.open(QIODevice::WriteOnly);
 358   QTextStream f(&file);
 359   f << "/*--------------------------------*- C++ -*----------------------------------*\\\n";
 360   f << "| =========                 |                                                 |\n";
 361   f << "| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n";
 362   f << "|  \\    /   O peration     | Version:  1.5                                   |\n";
 363   f << "|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |\n";
 364   f << "|    \\/     M anipulation  |                                                 |\n";
 365   f << "\\*---------------------------------------------------------------------------*/\n\n";
 366   f << "FoamFile\n";
 367   f << "{\n";
 368   f << "    version     2.0;\n";
 369   f << "    format      ascii;\n";
 370   f << "    class       polyBoundaryMesh;\n";
 371   f << "    location    \"constant/polyMesh\";\n";
 372   f << "    object      boundary;\n";
 373   f << "}\n\n";
 374   f << "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n\n";
 375   QSet<int> bcs;
 376   foreach (face_t F, m_Faces) {
 377     if (F.bc != 0) {
 378       bcs.insert(F.bc);
 379     }
 380   }
 381   f << bcs.size() << "\n(\n";
 382   QVector<patch_t> patch(bcs.size());
 383   int N_bc = 0;
 384   foreach (int bc, bcs) {
 385     int nFaces = 0;
 386     int startFace = -1;
 387     for (int i = 0; i < m_Faces.size(); ++i) {
 388       if (m_Faces[i].bc == bc) {
 389         ++nFaces;
 390         if (startFace == -1) {
 391           startFace = i;
 392         }
 393       }
 394     }
 395     if (startFace == -1) {
 396       EG_BUG;
 397     }
 398     patch_t P;
 399     P.bc = bc;
 400     P.startFace = startFace;
 401     P.nFaces = nFaces;
 402     patch[N_bc] = P;
 403     ++N_bc;
 404   }
 405   qSort(patch);
 406   foreach (patch_t P, patch) {
 407     BoundaryCondition BC = getBC(P.bc);
 408     QString num;
 409     num.setNum(P.bc);
 410     QString name = BC.getName();
 411     QString type = BC.getType();
 412     if (GuiMainWindow::pointer()->physicalTypeDefined(type)) {
 413       PhysicalBoundaryCondition PBC = GuiMainWindow::pointer()->getPhysicalBoundaryCondition(type);
 414       type = PBC.getFoamType();
 415     }
 416     if (name == "unknown") {
 417       name = "BC_" + num;
 418     }
 419     f << "    " << name << "\n";
 420     f << "    {\n";
 421     f << "        type        " << type << ";\n";
 422     f << "        nFaces      " << P.nFaces << ";\n";
 423     f << "        startFace   " << P.startFace + faces_offset << ";\n";
 424     f << "    }\n";
 425   }
 426   f << ")\n\n";
 427   f << "// ************************************************************************* //\n\n\n";
 428 }
 429
 430 void SimpleFoamWriter::operateOnGivenFileName()
 431 {
 432   if (isValid()) {
 433     QString p1 = getFileName();
 434     QString p2 = p1 + "/constant";
 435     QDir d1(p1);
 436     QDir d2(p2);
 437     if (!d1.exists()) {
 438       EG_BUG;
 439     }
 440     if (!d2.exists()) {
 441       d1.mkdir("constant");
 442       d2 = QDir(p2);
 443     }
 444     d1 = d2;
 445     p1 = p2;
 446     p2 = p1 + "/polyMesh";
 447     d2 = QDir(p2);
 448     if (!d2.exists()) {
 449       d1.mkdir("polyMesh");
 450     }
 451     m_Path = getFileName() + "/constant/polyMesh/";
 452     if (!QDir(m_Path).exists()) {
 453       EG_BUG;
 454     }
 455     createFaces();
 456     writePoints();
 457     writeFaces();
 458     writeOwner();
 459     writeNeighbour();
 460     writeBoundary();
 461   }
 462 }
 463
 464 void SimpleFoamWriter::operate()
 465 {
 466   try {
 467     readOutputDirectory();
 468     if (isValid()) {
 469       QString p1 = getFileName();
 470       QString p2 = p1 + "/constant";
 471       QDir d1(p1);
 472       QDir d2(p2);
 473       if (!d1.exists()) {
 474         EG_BUG;
 475       }
 476       if (!d2.exists()) {
 477         d1.mkdir("constant");
 478         d2 = QDir(p2);
 479       }
 480       d1 = d2;
 481       p1 = p2;
 482       p2 = p1 + "/polyMesh";
 483       d2 = QDir(p2);
 484       if (!d2.exists()) {
 485         d1.mkdir("polyMesh");
 486       }
 487       m_Path = getFileName() + "/constant/polyMesh/";
 488       if (!QDir(m_Path).exists()) {
 489         EG_BUG;
 490       }
 491       createFaces();
 492       writePoints();
 493       writeFaces();
 494       writeOwner();
 495       writeNeighbour();
 496       writeBoundary();
 497     }
 498   } catch (Error err) {
 499     err.display();
 500   }
 501 }