From af8cd99fc893e04274d4c58b07b0a64be68c954c Mon Sep 17 00:00:00 2001 From: Oliver Gloth Date: Fri, 19 Apr 2024 09:14:01 +0200 Subject: [PATCH] limited volume meshing to boundary layer only --- src/libengrid/boundarylayeroperation.cpp | 125 ++++++++++++++++++----------- src/libengrid/boundarylayeroperation.h | 4 + src/libengrid/createboundarylayershell.cpp | 9 +-- src/libengrid/createvolumemesh.cpp | 35 +++++++- src/libengrid/guicreatesurfacemesh.cpp | 8 +- src/libengrid/guicreatesurfacemesh.ui | 55 ++++++++++--- src/libengrid/guicreatevolumemesh.cpp | 1 + src/libengrid/meshpartition.cpp | 10 ++- 8 files changed, 177 insertions(+), 70 deletions(-) diff --git a/src/libengrid/boundarylayeroperation.cpp b/src/libengrid/boundarylayeroperation.cpp index ea74e2d0..ddebfd3b 100644 --- a/src/libengrid/boundarylayeroperation.cpp +++ b/src/libengrid/boundarylayeroperation.cpp @@ -27,6 +27,7 @@ #include "geometrysmoother.h" #include "vtkEgPolyDataToUnstructuredGridFilter.h" +//#include #include #include #include @@ -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 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 c23c93d8..8935eac7 100644 --- 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 m_RelativeHeights; + + int m_DrnumBuffer = 4; protected: // methods diff --git a/src/libengrid/createboundarylayershell.cpp b/src/libengrid/createboundarylayershell.cpp index 74d5522a..25f6e72e 100644 --- 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 401bc157..49e25af3 100755 --- 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 #include 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 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 4d8794c1..f6da3cea 100644 --- 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 7d279709..58831d1d 100644 --- a/src/libengrid/guicreatesurfacemesh.ui +++ b/src/libengrid/guicreatesurfacemesh.ui @@ -6,8 +6,8 @@ 0 0 - 905 - 875 + 1063 + 878 @@ -24,7 +24,7 @@ true - 0 + 3 @@ -398,7 +398,7 @@ - stretching ration between layers + stretching ratio between layers @@ -438,7 +438,7 @@ - ratio between last layer and farfield + ratio between last layer and first layer @@ -471,7 +471,7 @@ 0.050000000000000 - 1.200000000000000 + 1.500000000000000 @@ -558,9 +558,18 @@ + + + 150 + 0 + + grouping angle + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + @@ -641,16 +650,16 @@ - 0.050000000000000 + 1.000000000000000 - 1.000000000000000 + 1000.000000000000000 - 0.050000000000000 + 0.100000000000000 - 0.500000000000000 + 10.000000000000000 @@ -786,6 +795,32 @@ + + + + 0 + + + 10 + + + 3 + + + + + + + Qt::LeftToRight + + + buffer layers + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + diff --git a/src/libengrid/guicreatevolumemesh.cpp b/src/libengrid/guicreatevolumemesh.cpp index 08766d5a..3c7b5fd0 100644 --- 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 52f99c09..986c2f45 100644 --- 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; -- 2.11.4.GIT