limited volume meshing to boundary layer only

diff --git a/src/libengrid/boundarylayeroperation.cpp b/src/libengrid/boundarylayeroperation.cpp
index ea74e2d..ddebfd3 100644 (file)
--- a/src/libengrid/boundarylayeroperation.cpp
+++ b/src/libengrid/boundarylayeroperation.cpp
@@ -27,6 +27,7 @@
 #include "geometrysmoother.h"
 #include "vtkEgPolyDataToUnstructuredGridFilter.h"
 
+//#include <boost/multiprecision/cpp_int/cpp_int_config.hpp>
 #include <vtkDataSetSurfaceFilter.h>
 #include <vtkIdList.h>
 #include <vtkLoopSubdivisionFilter.h>
@@ -106,6 +107,7 @@ void BoundaryLayerOperation::readSettings()
     m_UseGrouping = use_grouping;
     in >> m_GroupingAngle;
     m_GroupingAngle = deg2rad(m_GroupingAngle);
+    in >> m_NumBufferLayers;
   }
   m_ELSManagerBLayer.clear();
   m_ELSManagerBLayer.readBoundaryLayerRules(m_Grid);
@@ -216,6 +218,11 @@ void BoundaryLayerOperation::computeBoundaryLayerVectors()
   computeLimitPlaneNormals();
   writeBoundaryLayerVectors("blayer", 1);
   smoothBoundaryLayerVectors(m_NumBoundaryLayerVectorRelaxations);
+  for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
+    if (m_BoundaryLayerNode[id_node]) {
+      m_BoundaryLayerVectors[id_node].normalise();
+    }
+  }
   writeBoundaryLayerVectors("blayer", 2);
 
   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
@@ -233,14 +240,16 @@ void BoundaryLayerOperation::computeBoundaryLayerVectors()
         n.normalise();
         m_BoundaryLayerVectors[id_node] = n;
       }
-      if (num_bcs[id_node] > 1) {
-        m_BoundaryLayerVectors[id_node] = n_opt[id_node](m_BoundaryLayerVectors[id_node]);
-      }
-      if (!checkVector(m_BoundaryLayerVectors[id_node])) {
-        EG_ERR_RETURN("invalid layer vector computed");
-      }
     }
+    if (num_bcs[id_node] > 1) {
+      m_BoundaryLayerVectors[id_node] = n_opt[id_node](m_BoundaryLayerVectors[id_node]);
+    }
+    if (!checkVector(m_BoundaryLayerVectors[id_node])) {
+      EG_ERR_RETURN("invalid layer vector computed");
+    }
+    m_BoundaryLayerVectors[id_node].normalise();
   }
+  writeBoundaryLayerVectors("blayer", 3);
 
   computeHeights();
   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
@@ -251,7 +260,7 @@ void BoundaryLayerOperation::computeBoundaryLayerVectors()
       }
     }
   }
-  writeBoundaryLayerVectors("blayer", 3);
+  writeBoundaryLayerVectors("blayer", 4);
 
 }
 
@@ -546,54 +555,64 @@ void BoundaryLayerOperation::writeBoundaryLayerVectors(QString file_name, int co
 void BoundaryLayerOperation::computeDesiredHeights()
 {
   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
-
-  // first pass (intial height)
-  m_Height.fill(0, m_Grid->GetNumberOfPoints());
-  int k = 1;
+  //
+  // As a hack we will take the minimal height as the height for the first prism for the whole mesh
+  // This can be improved if (and when) we decide to go ahead with improving enGrid for use with DrNUM
+  //
+  double h_min = 1e10;
   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
     if (m_BoundaryLayerNode[id_node]) {
       vec3_t x;
       m_Grid->GetPoint(id_node, x.data());
-      double h0 = m_ELSManagerBLayer.minEdgeLength(x);
-      double h1 = cl->GetValue(id_node)*m_FarfieldRatio;
-      k = max(k, int(logarithm(m_StretchingRatio, h1/h0)));
+      h_min = std::min(m_ELSManagerBLayer.minEdgeLength(x), h_min);
     }
   }
+  //
+  // compute the heights (valid for the whole boundary layer)
+  //
+  std::vector<double> y_layer;
+  {
+    double dy     = h_min;
+    double dy_max = m_FarfieldRatio * h_min;
+    bool   done   = false;
+    //
+    y_layer.push_back(0.0);
+    y_layer.push_back(h_min);
+    //
+    while (!done) {
+      dy *= m_StretchingRatio;
+      if (dy > dy_max) {
+        dy   = dy_max;
+        done = true;
+      }
+      y_layer.push_back(y_layer.back() + dy);
+    }
+    //
+    // add buffer layers
+    //
+    for (int i = 0; i < m_NumBufferLayers; ++i) {
+      y_layer.push_back(y_layer.back() + dy);
+    }
+  }
+  m_RelativeHeights.resize(y_layer.size());
+  for (int i = 0; i < y_layer.size(); ++i) {
+    m_RelativeHeights[i] = y_layer[i]/y_layer.back();
+  }
+  //
+  // first pass (initial height)
+  //
+  m_Height.fill(0, m_Grid->GetNumberOfPoints());
+  int k = 1;
   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
     if (m_BoundaryLayerNode[id_node]) {
-      vec3_t x;
-      m_Grid->GetPoint(id_node, x.data());
-      double h0 = m_ELSManagerBLayer.minEdgeLength(x);
-      double h1 = cl->GetValue(id_node)*m_FarfieldRatio;
-      double s  = pow(h1/h0, 1.0/k);
-      double H  = h0*(1 - pow(s, k))/(1 - s);
-      if (!checkReal(H)) {
-        EG_ERR_RETURN("floating point error while computing heights");
-      }
-      if (H < 0) {
-        EG_ERR_RETURN("negative height computed");
-      }
-      if (H > 1000*h1) {
-        cout << H << ", " << h1 << endl;
-        EG_ERR_RETURN("unrealistically large height computed");
-      }
-      if (H < 1e-3*h0) {
-        EG_ERR_RETURN("unrealistically small height computed");
-      }
-      if (h1 < h0) {
-        QString h0_txt, h1_txt, id_txt;
-        h0_txt.setNum(h0);
-        h1_txt.setNum(h1);
-        id_txt.setNum(id_node);
-        EG_ERR_RETURN("h1 < h0 (" + h1_txt + " < " + h0_txt + ", for node " + id_txt + ")");
-      }
-      m_Height[id_node] = H;
+      m_Height[id_node] = y_layer.back();
     }
   }
-
-  m_NumLayers = k;
-
+  //
+  m_NumLayers = y_layer.size() - 1;
+  //
   // correct with angle between face normal and propagation direction (node normals)
+  //
   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
     if (m_BoundaryLayerNode[id_node]) {
       int N = 0;
@@ -745,7 +764,7 @@ void BoundaryLayerOperation::createSmoothShell()
   EG_VTKSP(vtkWindowedSincPolyDataFilter, smooth);
   smooth->SetInputConnection(subdiv->GetOutputPort());
   smooth->BoundarySmoothingOn();
-  smooth->FeatureEdgeSmoothingOn();
+  smooth->FeatureEdgeSmoothingOff();
   smooth->SetFeatureAngle(180);
   smooth->SetEdgeAngle(180);
   smooth->SetNumberOfIterations(100);
@@ -755,8 +774,20 @@ void BoundaryLayerOperation::createSmoothShell()
   smooth->SetPassBand(pb);
   smooth->Update();
 
+  /*
+  EG_VTKSP(vtkSmoothPolyDataFilter, smooth);
+  smooth->SetInputConnection(subdiv->GetOutputPort());
+  smooth->SetNumberOfIterations(10);
+  smooth->SetRelaxationFactor(0.1);
+  smooth->FeatureEdgeSmoothingOff();
+  smooth->BoundarySmoothingOn();
+  smooth->Update();
+  */
+
+
   EG_VTKSP(vtkEgPolyDataToUnstructuredGridFilter, pdata_to_grid);
-  pdata_to_grid->SetInputConnection(smooth->GetOutputPort());
+  //pdata_to_grid->SetInputConnection(smooth->GetOutputPort());
+  pdata_to_grid->SetInputConnection(subdiv->GetOutputPort());
   pdata_to_grid->Update();
   makeCopy(pdata_to_grid->GetOutput(), m_ShellGrid);
 
@@ -836,8 +867,8 @@ void BoundaryLayerOperation::smoothUsingBLVectors()
 {
   // create shell
   createSmoothShell();
-
   newHeightFromShellIntersect(1.0);
+  return;
 
   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 

diff --git a/src/libengrid/boundarylayeroperation.h b/src/libengrid/boundarylayeroperation.h
index c23c93d..8935eac 100644 (file)
--- a/src/libengrid/boundarylayeroperation.h
+++ b/src/libengrid/boundarylayeroperation.h
@@ -61,9 +61,13 @@ protected: // attributes
   int                       m_NumBoundaryLayerHeightRelaxations;
   double                    m_ShellPassBand;
   int                       m_NumLayers;
+  int                       m_NumBufferLayers;
   EdgeLengthSourceManager   m_ELSManagerBLayer;
   EdgeLengthSourceManager   m_ELSManagerSurface;
   vtkUnstructuredGrid*      m_ShellGrid;
+  std::vector<double>       m_RelativeHeights;
+
+  int m_DrnumBuffer = 4;
 
 
 protected: // methods

diff --git a/src/libengrid/createboundarylayershell.cpp b/src/libengrid/createboundarylayershell.cpp
index 74d5522..25f6e72 100644 (file)
--- a/src/libengrid/createboundarylayershell.cpp
+++ b/src/libengrid/createboundarylayershell.cpp
@@ -22,6 +22,7 @@
 #include "createboundarylayershell.h"
 
 #include "createvolumemesh.h"
+#include "math/mathvector.h"
 #include "swaptriangles.h"
 #include "deletetetras.h"
 #include "deletecells.h"
@@ -336,14 +337,10 @@ void CreateBoundaryLayerShell::createLayerNodes(vtkIdType id_node)
   m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
   vec3_t x2 = x1 + m_BoundaryLayerVectors[id_node];
 
-  double H  = m_BoundaryLayerVectors[id_node].abs();
-  double h  = H*(1.0 - m_StretchingRatio)/(1.0 - pow(m_StretchingRatio, m_NumLayers));
-  vec3_t dx = (1.0/H)*m_BoundaryLayerVectors[id_node];
-  vec3_t x  = x1;
   m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
+  double h1 = 0;
   for (int i = 1; i < m_NumLayers; ++i) {
-    x += h*dx;
-    h *= m_StretchingRatio;
+    vec3_t x = x1 + m_RelativeHeights[i]*(x2 - x1);
     m_PrismaticGrid->GetPoints()->SetPoint(i*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x.data());
   }
   m_PrismaticGrid->GetPoints()->SetPoint(m_NumLayers*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x2.data());

diff --git a/src/libengrid/createvolumemesh.cpp b/src/libengrid/createvolumemesh.cpp
index 401bc15..49e25af 100755 (executable)
--- a/src/libengrid/createvolumemesh.cpp
+++ b/src/libengrid/createvolumemesh.cpp
@@ -22,10 +22,13 @@
 #include "createvolumemesh.h"
 #include "deletetetras.h"
 #include "deletecells.h"
+#include "engrid.h"
 #include "guimainwindow.h"
+#include "meshpartition.h"
 #include "updatedesiredmeshdensity.h"
 #include "createboundarylayershell.h"
 
+#include <vtkUnstructuredGrid.h>
 #include <vtkXMLUnstructuredGridWriter.h>
 
 CreateVolumeMesh::CreateVolumeMesh()
@@ -106,16 +109,43 @@ void CreateVolumeMesh::operate()
     blayer();
     if (!m_Debug) {
       if (blayer.success()) {
-        createTetMesh(2, true);
+        //createTetMesh(2, true);
+        //
         vtkUnstructuredGrid *prismatic_grid = blayer.getPrismaticGrid();
         MeshPartition prismatic_part(prismatic_grid, true);
         QVector<vtkIdType> shell_cells;
         getSurfaceCells(blayer.getBoundaryLayerCodes(), shell_cells, m_Grid);
+        MeshPartition shell_part(m_Grid);
+        shell_part.setCells(shell_cells);
+        EG_VTKSP(vtkUnstructuredGrid, shell_grid);
+        shell_part.extractToVtkGrid(shell_grid);
+        //
+        auto new_bc = GuiMainWindow::pointer()->getBC(BoundaryCondition("transfer", "patch"));
+        {
+          EG_VTKDCC(vtkIntArray, cell_code, shell_grid, "cell_code");
+          for (vtkIdType id_cell = 0; id_cell < shell_grid->GetNumberOfCells(); ++id_cell) {
+            //
+            // set boundary code for shell cells (outside end of boundary layer mesh)
+            //
+            cell_code->SetValue(id_cell, new_bc.getCode());
+            //
+            // invert the direction of the normals of the shell cells
+            //
+            shell_grid->GetCells()->ReverseCellAtId(id_cell);
+          }
+        }
+        //
         DeleteCells delete_cells;
+        MeshPartition all(m_Grid, true);
         delete_cells.setGrid(m_Grid);
-        delete_cells.setCellsToDelete(shell_cells);
+        delete_cells.setCellsToDelete(all.getCells());
         delete_cells();
+        m_Part.setAllCells();
         m_Part.addPartition(prismatic_part);
+        //
+        MeshPartition new_shell_part(shell_grid, true);
+        m_Part.addPartition(new_shell_part);
+        //
       } else {
         cout << "An error occurred while creating the prismatic layers!" << endl;
       }
@@ -123,5 +153,6 @@ void CreateVolumeMesh::operate()
   } else {
     createTetMesh(2, true);
   }
+  m_Grid->Modified();
 }
 

diff --git a/src/libengrid/guicreatesurfacemesh.cpp b/src/libengrid/guicreatesurfacemesh.cpp
index 4d8794c..f6da3ce 100644 (file)
--- a/src/libengrid/guicreatesurfacemesh.cpp
+++ b/src/libengrid/guicreatesurfacemesh.cpp
@@ -97,7 +97,7 @@ GuiCreateSurfaceMesh::GuiCreateSurfaceMesh()
     }
   }
   setTextFromTable();
-  m_Ui.m_TextEditPrismaticLayers->setText(GuiMainWindow::pointer()->getXmlSection("engrid/blayer/rules"));
+  m_Ui.m_TextEditPrismaticLayers->setText(GuiMainWindow::pointer()->getXmlSection("engrid/blayer/rules").trimmed());
   
   connect(m_Ui.pushButton_SelectAll_BC, SIGNAL(clicked()), this, SLOT(SelectAll_BC()));
   connect(m_Ui.pushButton_ClearAll_BC, SIGNAL(clicked()), this, SLOT(ClearAll_BC()));
@@ -205,6 +205,7 @@ int GuiCreateSurfaceMesh::readSettings()
     in >> dv; m_Ui.m_DoubleSpinBoxBoundaryLayerRadarAngle->setValue(dv);
     in >> iv; m_Ui.m_CheckBoxBoundaryLayerNormalVectorGrouping->setChecked(iv);
     in >> dv; m_Ui.m_DoubleSpinBoxBoundaryLayerGroupingAngle->setValue(dv);
+    in >> iv; m_Ui.m_SpinBoxBufferLayers->setValue(iv); 
   }
 
   return(0);
@@ -259,6 +260,7 @@ int GuiCreateSurfaceMesh::writeSettings()
     if (m_Ui.m_CheckBoxBoundaryLayerNormalVectorGrouping->checkState() == Qt::Checked) out << "1\n";
     else                                                                               out << "0\n";
     out << m_Ui.m_DoubleSpinBoxBoundaryLayerGroupingAngle->value() << "\n";
+    out << m_Ui.m_SpinBoxBufferLayers->value() << "\n";
   }
   GuiMainWindow::pointer()->setXmlSection("engrid/blayer/settings", buffer);
 
@@ -289,8 +291,8 @@ void GuiCreateSurfaceMesh::ClearAll_BC()
 void GuiCreateSurfaceMesh::setTextFromTable()
 {
   //EG_STOPDATE("2015-06-01");
-  m_Ui.textEdit->setText(GuiMainWindow::pointer()->getXmlSection("engrid/surface/rules"));
-  m_Ui.m_TextEditPrismaticLayers->setText(GuiMainWindow::pointer()->getXmlSection("engrid/blayer/rules"));
+  m_Ui.textEdit->setText(GuiMainWindow::pointer()->getXmlSection("engrid/surface/rules").trimmed());
+  m_Ui.m_TextEditPrismaticLayers->setText(GuiMainWindow::pointer()->getXmlSection("engrid/blayer/rules").trimmed());
   return;
 
   QString text = "";

diff --git a/src/libengrid/guicreatesurfacemesh.ui b/src/libengrid/guicreatesurfacemesh.ui
index 7d27970..58831d1 100644 (file)
--- a/src/libengrid/guicreatesurfacemesh.ui
+++ b/src/libengrid/guicreatesurfacemesh.ui
@@ -6,8 +6,8 @@
    <rect>
     <x>0</x>
     <y>0</y>
-    <width>905</width>
-    <height>875</height>
+    <width>1063</width>
+    <height>878</height>
    </rect>
   </property>
   <property name="windowTitle">
@@ -24,7 +24,7 @@
       <bool>true</bool>
      </property>
      <property name="currentIndex">
-      <number>0</number>
+      <number>3</number>
      </property>
      <widget class="QWidget" name="tab_4">
       <attribute name="title">
@@ -398,7 +398,7 @@
          <item row="1" column="0">
           <widget class="QLabel" name="label_13">
            <property name="text">
-            <string>stretching ration between layers</string>
+            <string>stretching ratio between layers</string>
            </property>
           </widget>
          </item>
@@ -438,7 +438,7 @@
          <item row="2" column="0">
           <widget class="QLabel" name="label_15">
            <property name="text">
-            <string>ratio between last layer and farfield</string>
+            <string>ratio between last layer and first layer</string>
            </property>
           </widget>
          </item>
@@ -471,7 +471,7 @@
             <double>0.050000000000000</double>
            </property>
            <property name="value">
-            <double>1.200000000000000</double>
+            <double>1.500000000000000</double>
            </property>
           </widget>
          </item>
@@ -558,9 +558,18 @@
          </item>
          <item row="11" column="2">
           <widget class="QLabel" name="label_22">
+           <property name="minimumSize">
+            <size>
+             <width>150</width>
+             <height>0</height>
+            </size>
+           </property>
            <property name="text">
             <string>   grouping angle</string>
            </property>
+           <property name="alignment">
+            <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+           </property>
           </widget>
          </item>
          <item row="11" column="3">
@@ -641,16 +650,16 @@
          <item row="2" column="1">
           <widget class="QDoubleSpinBox" name="m_DoubleSpinBoxBoundaryLayerFarfieldRatio">
            <property name="minimum">
-            <double>0.050000000000000</double>
+            <double>1.000000000000000</double>
            </property>
            <property name="maximum">
-            <double>1.000000000000000</double>
+            <double>1000.000000000000000</double>
            </property>
            <property name="singleStep">
-            <double>0.050000000000000</double>
+            <double>0.100000000000000</double>
            </property>
            <property name="value">
-            <double>0.500000000000000</double>
+            <double>10.000000000000000</double>
            </property>
           </widget>
          </item>
@@ -786,6 +795,32 @@
            </property>
           </widget>
          </item>
+         <item row="2" column="3">
+          <widget class="QSpinBox" name="m_SpinBoxBufferLayers">
+           <property name="minimum">
+            <number>0</number>
+           </property>
+           <property name="maximum">
+            <number>10</number>
+           </property>
+           <property name="value">
+            <number>3</number>
+           </property>
+          </widget>
+         </item>
+         <item row="2" column="2">
+          <widget class="QLabel" name="label_28">
+           <property name="layoutDirection">
+            <enum>Qt::LeftToRight</enum>
+           </property>
+           <property name="text">
+            <string>buffer layers</string>
+           </property>
+           <property name="alignment">
+            <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+           </property>
+          </widget>
+         </item>
         </layout>
        </item>
       </layout>

diff --git a/src/libengrid/guicreatevolumemesh.cpp b/src/libengrid/guicreatevolumemesh.cpp
index 08766d5..3c7b5fd 100644 (file)
--- a/src/libengrid/guicreatevolumemesh.cpp
+++ b/src/libengrid/guicreatevolumemesh.cpp
@@ -69,6 +69,7 @@ void GuiCreateVolumeMesh::operate()
     MeshPartition rest(rest_grid, true);
     part.addPartition(rest);
   }
+  GuiMainWindow::pointer()->updateBoundaryCodes(true);
   resetOrientation(m_Grid);
   createIndices(m_Grid);
 }

diff --git a/src/libengrid/meshpartition.cpp b/src/libengrid/meshpartition.cpp
index 52f99c0..986c2f4 100644 (file)
--- a/src/libengrid/meshpartition.cpp
+++ b/src/libengrid/meshpartition.cpp
@@ -207,9 +207,15 @@ void MeshPartition::addPartition(const MeshPartition& part, double tol)
     for (vtkIdType id_pnode = 0; id_pnode < part.m_Grid->GetNumberOfPoints(); ++id_pnode) {
       vec3_t xp, x;
       part.m_Grid->GetPoint(id_pnode, xp.data());
+      bool found_match = false;
       vtkIdType id_node = loc->FindClosestPoint(xp.data());
-      m_Grid->GetPoint(id_node, x.data());
-      if ((x - xp).abs() < tol) {
+      if (id_node >= 0) {
+        m_Grid->GetPoint(id_node, x.data());
+        if ((x - xp).abs() < tol) {
+          found_match = true;
+        }
+      }
+      if (found_match) {
         pnode2node[id_pnode] = id_node;
       } else {
         pnode2node[id_pnode] = N;