src/libengrid/drnumwriter.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 "drnumwriter.h"
  23 #include "guimainwindow.h"
  24 #include "vtkEgNormalExtrusion.h"
  25 #include "padsurface.h"
  26
  27 DrNumWriter::DrNumWriter()
  28 {
  29   m_BackupGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
  30 }
  31
  32 QList<BoundaryCondition> DrNumWriter::getBcsOfType(QString type)
  33 {
  34   QList<BoundaryCondition> bcs_of_type;
  35   QSet<int> bcs = GuiMainWindow::pointer()->getAllBoundaryCodes();
  36   foreach (int bc_code, bcs) {
  37     BoundaryCondition bc = GuiMainWindow::pointer()->getBC(bc_code);
  38     PhysicalBoundaryCondition pbc = GuiMainWindow::pointer()->getPhysicalBoundaryCondition(bc.getType());
  39     if (pbc.getType() == type) {
  40       bcs_of_type << bc;
  41     }
  42   }
  43   qSort(bcs_of_type);
  44   return bcs_of_type;
  45 }
  46
  47 void DrNumWriter::readSettings()
  48 {
  49   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings");
  50   if(!buffer.isEmpty()) {
  51     QTextStream in(&buffer, QIODevice::ReadOnly);
  52     in >> m_MaximalEdgeLength;
  53     in >> m_MinimalEdgeLength;
  54     in >> m_GrowthFactor;
  55   }
  56 }
  57
  58 double DrNumWriter::edgeLength(QString bc_name)
  59 {
  60   QString rules_txt = GuiMainWindow::pointer()->getXmlSection("engrid/surface/rules");
  61   rules_txt = rules_txt.replace("\n", " ");
  62   rules_txt = rules_txt.trimmed();
  63   QStringList rules = rules_txt.split(";", QString::SkipEmptyParts);
  64   double h_min = EG_LARGE_REAL;
  65   foreach (QString rule, rules) {
  66     QStringList parts = rule.split("=");
  67     if (parts.size() > 1) {
  68       QString left  = parts[0].trimmed();
  69       double h = parts[1].trimmed().toDouble();
  70       QStringList or_parts = left.split("<OR>");
  71       foreach (QString or_part, or_parts) {
  72         or_part = or_part.trimmed();
  73         QStringList and_parts = or_part.split("<AND>");
  74         if (and_parts.size() == 1) {
  75           QString and_part = and_parts[0].trimmed();
  76           if (and_part == bc_name) {
  77             h_min = min(h, h_min);
  78           }
  79         }
  80       }
  81     }
  82   }
  83   return h_min;
  84 }
  85
  86 void DrNumWriter::prepareLevelSets(QList<BoundaryCondition> bcs, double distance)
  87 {
  88   bool extrude = distance > 1e-3;
  89   foreach (BoundaryCondition bc, bcs) {
  90     BoundaryCondition new_bc = GuiMainWindow::pointer()->getBC(BoundaryCondition(bc.getName(), "level-set"));
  91     BoundaryCondition tmp_bc = GuiMainWindow::pointer()->getBC(BoundaryCondition(bc.getName(), "auxilary"));
  92     MeshPartition part(m_Grid);
  93     part.setBC(bc.getCode());
  94     if (!part.isPlanar()) {
  95       PhysicalBoundaryCondition pbc = GuiMainWindow::pointer()->getPhysicalBoundaryCondition(bc.getType());
  96       QString msg =  "only planar surfaces allowed for boundaries of type \"" + pbc.getType() + "\"";
  97       EG_ERR_RETURN(msg);
  98     }
  99     double A, P, Dh;
 100     vec3_t x, n;
 101     part.calcPlanarSurfaceMetrics(Dh, A, P, x, n);
 102     if (extrude) {
 103       QList<double> h;
 104       h << distance*Dh + 10*m_MaximalEdgeLength;
 105       x += h.last()*n;
 106       part.extrude(n, h, bc, false, false, true, BoundaryCondition(), BoundaryCondition(), tmp_bc);
 107       part.setBC(tmp_bc.getCode());
 108     }
 109     part.duplicate();
 110     part.resetBC(new_bc.getName(), new_bc.getType());
 111
 112     //part.setBC(new_bc.getCode());
 113     //part.scale((Dh + 10*m_MaximalEdgeLength)/Dh, x);
 114
 115     PadSurface pad;
 116     pad.setGrid(m_Grid);
 117     pad.addBC(new_bc.getCode());
 118     pad.relativeOff();
 119     pad.setDistance(10*m_MaximalEdgeLength);
 120     pad.setNewBC(new_bc.getCode());
 121     pad();
 122     part.setBC(new_bc.getCode());
 123
 124     if (extrude) {
 125       part.translate(-10*m_MaximalEdgeLength*n);
 126     }
 127     part.writeSTL(getFileName() + "/engrid/" + new_bc.getName() + ".stl");
 128     QFile info_file(getFileName() + "/engrid/" + new_bc.getName() + ".dnc");
 129     info_file.open(QIODevice::WriteOnly);
 130     QTextStream info(&info_file);
 131     PhysicalBoundaryCondition pbc = GuiMainWindow::pointer()->getPhysicalBoundaryCondition(bc.getType());
 132     int max_size = 6;
 133     for (int i = 0; i < pbc.getNumVars(); ++i) {
 134       max_size = max(max_size, pbc.getVarName(i).size());
 135     }
 136     info << "string " + QString("name").leftJustified(max_size) << " = " << bc.getName() << ";\n";
 137     info << "string " + QString("type").leftJustified(max_size) << " = " << pbc.getType() << ";\n";
 138     info << "real   " + QString("A").leftJustified(max_size) << " = " << A << ";\n";
 139     info << "real   " + QString("P").leftJustified(max_size) << " = " << P << ";\n";
 140     info << "real   " + QString("Dh").leftJustified(max_size) << " = " << Dh << ";\n";
 141     info << "vec3_t " + QString("centre").leftJustified(max_size) << " = (" << x[0] << ", " << x[1] << ", " << x[2] << ");\n";
 142     info << "vec3_t " + QString("normal").leftJustified(max_size) << " = (" << n[0] << ", " << n[1] << ", " << n[2] << ");\n";
 143     for (int i = 0; i < pbc.getNumVars(); ++i) {
 144       info << pbc.getVarType(i).leftJustified(7) << pbc.getVarName(i).leftJustified(max_size) << " = " << pbc.getVarValueAsString(i) << ";\n";
 145     }
 146     double h = edgeLength(bc.getName());
 147     if (h < m_MaximalEdgeLength) {
 148       info << "real   " + QString("h").leftJustified(max_size) << " = " << h << ";\n";
 149     }
 150   }
 151 }
 152
 153 void DrNumWriter::prepareWallLevelSets(QList<BoundaryCondition> bcs)
 154 {
 155   foreach (BoundaryCondition bc, bcs) {
 156     MeshPartition part(m_Grid);
 157     part.setBC(bc.getCode());
 158     part.writeSTL(getFileName() + "/engrid/" + bc.getName() + ".stl");
 159     QFile info_file(getFileName() + "/engrid/" + bc.getName() + ".dnc");
 160     info_file.open(QIODevice::WriteOnly);
 161     QTextStream info(&info_file);
 162     PhysicalBoundaryCondition pbc = GuiMainWindow::pointer()->getPhysicalBoundaryCondition(bc.getType());
 163     int max_size = 4;
 164     for (int i = 0; i < pbc.getNumVars(); ++i) {
 165       max_size = max(max_size, pbc.getVarName(i).size());
 166     }
 167     info << "string " + QString("name").leftJustified(max_size) << " = " << bc.getName() << ";\n";
 168     info << "string " + QString("type").leftJustified(max_size) << " = " << pbc.getType() << ";\n";
 169     for (int i = 0; i < pbc.getNumVars(); ++i) {
 170       info << pbc.getVarType(i).leftJustified(7) << pbc.getVarName(i).leftJustified(max_size) << " = " << pbc.getVarValueAsString(i) << ";\n";
 171     }
 172     double h = edgeLength(bc.getName());
 173     if (h < m_MaximalEdgeLength) {
 174       info << "real   " + QString("h").leftJustified(max_size) << " = " << h << ";\n";
 175     }
 176   }
 177 }
 178
 179 void DrNumWriter::computeBBox()
 180 {
 181   m_X1 = vec3_t(EG_LARGE_REAL, EG_LARGE_REAL, EG_LARGE_REAL);
 182   m_X2 = -1*m_X1;
 183   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 184     vec3_t x;
 185     m_Grid->GetPoint(id_node, x.data());
 186     for (int i = 0; i < 3; ++i) {
 187       m_X1[i] = min(m_X1[i], x[i]);
 188       m_X2[i] = max(m_X2[i], x[i]);
 189     }
 190   }
 191 }
 192
 193 void DrNumWriter::writeGlobals()
 194 {
 195   QFile info_file(getFileName() + "/engrid/global.dnc");
 196   info_file.open(QIODevice::WriteOnly);
 197   QTextStream info(&info_file);
 198
 199   info << "string name  = global;\n";
 200   info << "real   h_max = " << m_MaximalEdgeLength << ";\n";
 201   info << "real   gf    = " << m_GrowthFactor << ";\n";
 202   info << "vector x1    = (" << m_X1[0] << ", " << m_X1[1] << ", " << m_X1[2] << ");\n";
 203   info << "vector x2    = (" << m_X2[0] << ", " << m_X2[1] << ", " << m_X2[2] << ");\n";
 204 }
 205
 206 void DrNumWriter::backup()
 207 {
 208   makeCopy(m_Grid, m_BackupGrid);
 209 }
 210
 211 void DrNumWriter::restore()
 212 {
 213   makeCopy(m_BackupGrid, m_Grid);
 214   GuiMainWindow::pointer()->updateBoundaryCodes(true);
 215   updateNodeIndex(m_Grid);
 216   updateCellIndex(m_Grid);
 217 }
 218
 219 void DrNumWriter::operate()
 220 {
 221   readOutputDirectory();
 222   if (isValid()) {
 223     backup();
 224     readSettings();
 225
 226     QList<BoundaryCondition> turb_duct_in_bcs  = getBcsOfType("turbulent-duct-inlet");
 227     QList<BoundaryCondition> lam_duct_in_bcs   = getBcsOfType("laminar-duct-inlet");
 228     QList<BoundaryCondition> lam_ext_in_bcs    = getBcsOfType("inflow");
 229     QList<BoundaryCondition> out_bcs           = getBcsOfType("outflow");
 230     QList<BoundaryCondition> cyclic_bcs        = getBcsOfType("cyclic");
 231     QList<BoundaryCondition> sym_bcs           = getBcsOfType("symmetry");
 232     QList<BoundaryCondition> turb_wall_bcs     = getBcsOfType("DrNUM-turbulent-wall");
 233     QList<BoundaryCondition> lam_wall_bcs      = getBcsOfType("wall");
 234     QList<BoundaryCondition> slip_wall_bcs     = getBcsOfType("inviscid-wall");
 235
 236     QString root_path = getFileName();
 237     QDir    root_dir(root_path);
 238     QString engrid_path = "engrid";
 239     QDir    engrid_dir(root_path + "/" + engrid_path);
 240     if (!engrid_dir.exists()) {
 241       root_dir.mkdir(engrid_path);
 242     } else {
 243       // delete all files in the levelset directory
 244       QStringList files = engrid_dir.entryList(QDir::Files);
 245       foreach (QString file, files) {
 246         engrid_dir.remove(file);
 247       }
 248     }
 249
 250     prepareLevelSets(turb_duct_in_bcs,  3.0);
 251     prepareLevelSets(lam_duct_in_bcs,   0.0);
 252     prepareLevelSets(lam_ext_in_bcs,    0.0);
 253     prepareLevelSets(out_bcs,           0.0);
 254     prepareLevelSets(cyclic_bcs,        0.0);
 255     prepareLevelSets(sym_bcs,           0.0);
 256
 257     prepareWallLevelSets(turb_wall_bcs);
 258     prepareWallLevelSets(lam_wall_bcs);
 259     prepareWallLevelSets(slip_wall_bcs);
 260
 261     computeBBox();
 262     writeGlobals();
 263     writeSolverParameters(getFileName());
 264     restore();
 265   }
 266 }