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
  25 DrNumWriter::DrNumWriter()
  26 {
  27   m_OverlapLayers = 2;
  28 }
  29
  30 QString DrNumWriter::boundaryCode(vtkIdType id_cell, int i)
  31 {
  32   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  33   QString bcode = "0";
  34   vtkIdType id_neigh = m_Part.c2cGG(id_cell, i);
  35   if (id_neigh != -1) {
  36     if (m_CellToCartPatch[id_neigh] == -1) {
  37       if (isSurface(id_neigh, m_Grid)) {
  38         BoundaryCondition bc = GuiMainWindow::pointer()->getBC(cell_code->GetValue(id_neigh));
  39         if (bc.getName() != "unknown") {
  40           bcode = bc.getType();
  41         }
  42       }
  43     }
  44   }
  45   return bcode;
  46 }
  47
  48 void DrNumWriter::computeMeshDensity()
  49 {
  50   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings").replace("\n", " ");
  51   if (!buffer.isEmpty()) {
  52     QTextStream in(&buffer, QIODevice::ReadOnly);
  53     in >> m_MaxEdgeLength;
  54     in >> m_MinEdgeLength;
  55     in >> m_GrowthFactor;
  56   } else {
  57     m_MaxEdgeLength = 1000.0;
  58     m_MinEdgeLength = 0.0;
  59     m_GrowthFactor = 1.5;
  60   }
  61   m_ELSManager.read();
  62   m_H.fill(m_MaxEdgeLength, m_Grid->GetNumberOfCells());
  63
  64   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  65     vec3_t x = cellCentre(m_Grid, id_cell);
  66     double cl_src = m_ELSManager.minEdgeLength(x);
  67     if (cl_src > 0) {
  68       if (cl_src < m_H[id_cell]) {
  69         m_H[id_cell] = cl_src;
  70       }
  71     }
  72   }
  73
  74   double H_min = 1e99;
  75   vtkIdType id_min = -1;
  76   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  77     if (m_H[id_cell] < H_min && isVolume(id_cell, m_Grid)) {
  78       id_min = id_cell;
  79       H_min = m_H[id_cell];
  80     }
  81   }
  82   if (id_min < 0) {
  83     EG_BUG;
  84   }
  85
  86   QVector<bool> marked(m_Grid->GetNumberOfCells(), false);
  87   marked[id_min] = true;
  88   bool done = false;
  89   while (!done) {
  90     done = true;
  91     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  92       if (marked[id_cell] && m_H[id_cell] <= H_min && isVolume(id_cell, m_Grid)) {
  93         vtkIdType num_pts, *pts;
  94         m_Grid->GetCellPoints(id_cell, num_pts, pts);
  95         for (int i = 0; i < num_pts; ++i) {
  96           for (int j = 0; j < m_Part.n2cGSize(pts[i]); ++j) {
  97             vtkIdType id_neigh = m_Part.n2cGG(pts[i], j);
  98             if (!marked[id_neigh] && isVolume(id_neigh, m_Grid)) {
  99               double h = m_GrowthFactor*m_H[id_cell];
 100               if (h < 0) {
 101                 EG_BUG;
 102               }
 103               m_H[id_neigh] = min(m_H[id_neigh], h);
 104               marked[id_neigh] = true;
 105               done = false;
 106             }
 107           }
 108         }
 109       }
 110     }
 111     H_min *= m_GrowthFactor;
 112   }
 113
 114   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 115     vec3_t x = cellCentre(m_Grid, id_cell);
 116     double cl_src = m_ELSManager.minEdgeLength(x);
 117     if (cl_src > 0) {
 118       if (cl_src < m_H[id_cell]) {
 119         m_H[id_cell] = cl_src;
 120       }
 121     }
 122   }
 123
 124 }
 125
 126 void DrNumWriter::createAllCartPatches()
 127 {
 128   m_CellToCartPatch.fill(-1, m_Grid->GetNumberOfCells());
 129   int num_cart_patches = 0;
 130   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 131     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 132     if (type_cell == VTK_HEXAHEDRON) {
 133       m_CellToCartPatch[id_cell] = num_cart_patches;
 134       ++num_cart_patches;
 135     }
 136   }
 137   m_CartPatches.resize(num_cart_patches);
 138   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 139   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 140     if (m_CellToCartPatch[id_cell] != -1) {
 141       cart_patch_t patch = m_CartPatches[m_CellToCartPatch[id_cell]];
 142       patch.id_cell = id_cell;
 143       vtkIdType num_pts, *pts;
 144       vec3_t x[8];
 145       m_Grid->GetCellPoints(id_cell, num_pts, pts);
 146       if (num_pts != 8) {
 147         EG_BUG;
 148       }
 149       vec3_t xc = cellCentre(m_Grid, id_cell);
 150       for (int i = 0; i < 8; ++i) {
 151         m_Grid->GetPoint(pts[i], x[i].data());
 152       }
 153       double h = m_H[id_cell];
 154       vec3_t gi = 0.25*(x[1] + x[2] + x[6] + x[5] - x[0] - x[3] - x[7] - x[4]);
 155       vec3_t gj = 0.25*(x[4] + x[5] + x[6] + x[7] - x[0] - x[1] - x[2] - x[3]);
 156       gi.normalise();
 157       gj.normalise();
 158       vec3_t gk = gi.cross(gj);
 159       double Li = 0;
 160       double Lj = 0;
 161       double Lk = 0;
 162       for (int i = 0; i < 8; ++i) {
 163         vec3_t Dx = x[i] - xc;
 164         Li = max(2*fabs(Dx*gi), Li);
 165         Lj = max(2*fabs(Dx*gj), Lj);
 166         Lk = max(2*fabs(Dx*gk), Lk);
 167       }
 168       int Ni = max(5, int(Li/h));
 169       int Nj = max(5, int(Lj/h));
 170       int Nk = max(5, int(Lk/h));
 171       double hi = Li/Ni;
 172       double hj = Lj/Nj;
 173       double hk = Lk/Nk;
 174       vec3_t x0 = xc - 0.5*(Li*gi + Lj*gj + Lk*gk);
 175
 176       {
 177         QString bc = boundaryCode(id_cell, 4);
 178         patch.bx = bc;
 179         if (bc == "0") {
 180           patch.sl_i1 = m_OverlapLayers;
 181           Ni += m_OverlapLayers;
 182           Li += m_OverlapLayers*hi;
 183           x0 -= m_OverlapLayers*hi*gi;
 184         } else {
 185           patch.sl_i1 = 0;
 186         }
 187       }
 188       {
 189         QString bc = boundaryCode(id_cell, 5);
 190         patch.bX = bc;
 191         if (bc == "0") {
 192           patch.sl_i2 = m_OverlapLayers;
 193           Ni += m_OverlapLayers;
 194           Li += m_OverlapLayers*hi;
 195         } else {
 196           patch.sl_i2 = 0;
 197         }
 198       }
 199       {
 200         QString bc = boundaryCode(id_cell, 0);
 201         patch.by = bc;
 202         if (bc == "0") {
 203           patch.sl_j1 = m_OverlapLayers;
 204           Nj += m_OverlapLayers;
 205           Lj += m_OverlapLayers*hj;
 206           x0 -= m_OverlapLayers*hj*gj;
 207         } else {
 208           patch.sl_j1 = 0;
 209         }
 210       }
 211       {
 212         QString bc = boundaryCode(id_cell, 1);
 213         patch.bY = bc;
 214         if (bc == "0") {
 215           patch.sl_j2 = m_OverlapLayers;
 216           Nj += m_OverlapLayers;
 217           Lj += m_OverlapLayers*hj;
 218         } else {
 219           patch.sl_j2 = 0;
 220         }
 221       }
 222       {
 223         QString bc = boundaryCode(id_cell, 3);
 224         patch.bz = bc;
 225         if (bc == "0") {
 226           patch.sl_k1 = m_OverlapLayers;
 227           Nk += m_OverlapLayers;
 228           Lk += m_OverlapLayers*hk;
 229           x0 -= m_OverlapLayers*hk*gk;
 230         } else {
 231           patch.sl_k1 = 0;
 232         }
 233       }
 234       {
 235         QString bc = boundaryCode(id_cell,2);
 236         patch.bZ = bc;
 237         if (bc == "0") {
 238           patch.sl_k2 = m_OverlapLayers;
 239           Nk += m_OverlapLayers;
 240           Lk += m_OverlapLayers*hk;
 241         } else {
 242           patch.sl_k2 = 0;
 243         }
 244       }
 245
 246       patch.gi = gi;
 247       patch.gj = gj;
 248       patch.Li = Li;
 249       patch.Lj = Lj;
 250       patch.Lk = Lk;
 251       patch.Ni = Ni;
 252       patch.Nj = Nj;
 253       patch.Nk = Nk;
 254       patch.x0 = x0;
 255
 256       patch.fx = "fx";
 257       patch.fy = "fy";
 258       patch.fz = "fz";
 259       patch.s  = "0";
 260
 261       m_CartPatches[m_CellToCartPatch[id_cell]] = patch;
 262     }
 263   }
 264 }
 265
 266 void DrNumWriter::writeCartPatches(QTextStream &s)
 267 {
 268   foreach (cart_patch_t patch, m_CartPatches) {
 269     s << "1001\n";
 270     s << "{\n";
 271     s << "  " << patch.x0[0] << " " << patch.x0[1] << " " << patch.x0[2] << "\n";
 272     s << "  " << patch.gi[0] << " " << patch.gi[1] << " " << patch.gi[2] << "\n";
 273     s << "  " << patch.gj[0] << " " << patch.gj[1] << " " << patch.gj[2] << "\n";
 274     s << "  1.0\n";
 275     s << "  " << patch.Ni << " " << patch.Nj << " " << patch.Nk << "\n";
 276     s << "  " << patch.sl_i1 << " " << patch.sl_i2 << " " << patch.sl_j1 << " " << patch.sl_j2 << " " << patch.sl_k1 << " " << patch.sl_k2 << "\n";
 277     s << "  " << patch.Li << " " << patch.Lj << " " << patch.Lk << "\n";
 278     s << "  " << patch.fx << " " << patch.fy << " " << patch.fz << "\n";
 279     s << "  " << patch.bx << " " << patch.bX << " " << patch.by << " " << patch.bY << " " << patch.bz << " " << patch.bZ << "\n";
 280     s << "  " << patch.s << "\n";
 281     s << "}\n";
 282     s << "\n";
 283   }
 284 }
 285
 286 void DrNumWriter::operate()
 287 {
 288   try {
 289     readOutputDirectory();
 290     if (isValid()) {
 291       QString p1 = getFileName();
 292       QString p2 = p1 + "/patches";
 293       QDir d1(p1);
 294       QDir d2(p2);
 295       if (!d1.exists()) {
 296         EG_BUG;
 297       }
 298       if (!d2.exists()) {
 299         d1.mkdir("patches");
 300         d2 = QDir(p2);
 301       }
 302       d1 = d2;
 303       p1 = p2;
 304       p2 = p1 + "/complex";
 305       d2 = QDir(p2);
 306       if (!d2.exists()) {
 307         d1.mkdir("complex");
 308       };
 309       m_PatchFile   = getFileName() + "/patches/from_enGrid";
 310       m_ComplexPath = getFileName() + "/patches/complex";
 311       if (!QDir(m_ComplexPath).exists()) {
 312         EG_BUG;
 313       }
 314       QFile file(m_PatchFile);
 315       file.open(QIODevice::WriteOnly);
 316       QTextStream s(&file);
 317       computeMeshDensity();
 318       createAllCartPatches();
 319       writeCartPatches(s);
 320       s << "0\n";
 321     }
 322   } catch (Error err) {
 323     err.display();
 324   }
 325 }