From 8612d75b099c035ee9712cbba1f3dca2802ce426 Mon Sep 17 00:00:00 2001 From: James Edgeley <james.edgeley@ukaea.uk> Date: Fri, 12 Jul 2024 13:36:35 +0000 Subject: [PATCH] Feature/NekMesh-Revolve --- CHANGELOG.md | 3 + docs/user-guide/utilities/nekmesh.tex | 24 + library/NekMesh/CMakeLists.txt | 2 + .../Module/ProcessModules/ProcessExtrude.cpp | 7 +- .../Module/ProcessModules/ProcessRevolve.cpp | 487 ++++++++++++++++++ .../Module/ProcessModules/ProcessRevolve.h | 68 +++ utilities/NekMesh/CMakeLists.txt | 1 + utilities/NekMesh/Tests/Nektar++/revolve.tst | 58 +++ utilities/NekMesh/Tests/Nektar++/revolve.xml | 54 ++ 9 files changed, 699 insertions(+), 5 deletions(-) create mode 100644 library/NekMesh/Module/ProcessModules/ProcessRevolve.cpp create mode 100644 library/NekMesh/Module/ProcessModules/ProcessRevolve.h create mode 100644 utilities/NekMesh/Tests/Nektar++/revolve.tst create mode 100644 utilities/NekMesh/Tests/Nektar++/revolve.xml diff --git a/CHANGELOG.md b/CHANGELOG.md index 5cedfdf949..f6f4a85335 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,9 @@ v5.7.0 **ShallowWaterSolver** - Implement implicit time-discritization (!1784) +**NekMesh** +- Added revolve module (!1825) + **Miscellaneous** - Use std::stod instead of boost::lexical_cast<NekDouble> (!1819) diff --git a/docs/user-guide/utilities/nekmesh.tex b/docs/user-guide/utilities/nekmesh.tex index 60434e4aa8..3fbde4ed5a 100644 --- a/docs/user-guide/utilities/nekmesh.tex +++ b/docs/user-guide/utilities/nekmesh.tex @@ -888,6 +888,30 @@ requires two parameters: length which determines how long the $z$ extrusion will be and layers, the number of elements in the $z$ direction. +\subsection{2D mesh revolution} + +This module allows a 2D mesh, quads, triangles or both, to be revolved around +the $y$ axis to make a simple 3D mesh made of prisms and hexahedra. It is also +capable of revolving the high-order curvature within the 2D mesh. The module +requires two parameters: + +\begin{lstlisting}[style=BashInputStyle] + NekMesh -m revolve:layers=n:angle=a 2D.xml 3D.xml +\end{lstlisting} +angle which determines the angle of revolution ($\theta$) around the $y$ axis and +layers, the number of elements in the $\theta$ direction. + +The angle can also be set to "full" for a full revolution of $2\pi$, although +this is the default setting: + +\begin{lstlisting}[style=BashInputStyle] +NekMesh -m revolve:layers=n:angle=full 2D.xml 3D.xml +\end{lstlisting} + +The angle must be between $0$ and $2\pi$. +The 2D mesh must not cross or touch the $y$ axis i.e. $x>0$ for all points. +By default the mesh is revolved through the full $2\pi$. + \subsection{Variational Optimisation} \label{sec:varopti} This module can correct invalid and improve the quality of elements in diff --git a/library/NekMesh/CMakeLists.txt b/library/NekMesh/CMakeLists.txt index 15fae92f9c..3cc4e5162e 100644 --- a/library/NekMesh/CMakeLists.txt +++ b/library/NekMesh/CMakeLists.txt @@ -30,6 +30,7 @@ SET(NEKMESH_SOURCES Module/ProcessModules/ProcessOptiExtract.cpp Module/ProcessModules/ProcessInsertSurface.cpp Module/ProcessModules/ProcessExtrude.cpp + Module/ProcessModules/ProcessRevolve.cpp Module/ProcessModules/ProcessVarOpti/ProcessVarOpti.cpp Module/ProcessModules/ProcessVarOpti/PreProcessing.cpp Module/ProcessModules/ProcessVarOpti/NodeOpti.cpp @@ -89,6 +90,7 @@ SET(NEKMESH_HEADERS Module/ProcessModules/ProcessOptiExtract.h Module/ProcessModules/ProcessInsertSurface.h Module/ProcessModules/ProcessExtrude.h + Module/ProcessModules/ProcessRevolve.h Module/ProcessModules/ProcessVarOpti/ProcessVarOpti.h Module/ProcessModules/ProcessVarOpti/NodeOpti.h Module/ProcessModules/ProcessVarOpti/ElUtil.h diff --git a/library/NekMesh/Module/ProcessModules/ProcessExtrude.cpp b/library/NekMesh/Module/ProcessModules/ProcessExtrude.cpp index 5eda6568e9..9a921581ea 100644 --- a/library/NekMesh/Module/ProcessModules/ProcessExtrude.cpp +++ b/library/NekMesh/Module/ProcessModules/ProcessExtrude.cpp @@ -129,7 +129,7 @@ void ProcessExtrude::Process() } } - EdgeSet es = m_mesh->m_edgeSet; // copy edges for curvature + EdgeSet esOld = m_mesh->m_edgeSet; // copy edges for curvature for (int j = 0; j < nLayers; ++j) { @@ -186,7 +186,7 @@ void ProcessExtrude::Process() ProcessComposites(); // Copy edge information - for (auto &edge : es) + for (auto &edge : esOld) { if (edge->m_edgeNodes.size() > 0) { @@ -381,9 +381,6 @@ void ProcessExtrude::Process() // corresponding composite if (inCommon == 4) { - m_log(VERBOSE) - << "Face " << itQ->GetId() << "\t" - << "Composite " << itOc.second->m_id << "\t" << endl; auto newC = m_mesh->m_composite.find(itOc.second->m_id); // Quad ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, diff --git a/library/NekMesh/Module/ProcessModules/ProcessRevolve.cpp b/library/NekMesh/Module/ProcessModules/ProcessRevolve.cpp new file mode 100644 index 0000000000..d5a8ac0716 --- /dev/null +++ b/library/NekMesh/Module/ProcessModules/ProcessRevolve.cpp @@ -0,0 +1,487 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File: ProcessRevolve.cpp +// +// For more information, please see: http://www.nektar.info/ +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: Revolve a two-dimensional mesh around the y-axis to form a +// three-dimensional mesh. Arguments are "angle" and "layers" +// +//////////////////////////////////////////////////////////////////////////////// + +#include "ProcessRevolve.h" +#include <NekMesh/MeshElements/Element.h> + +using namespace std; +using namespace Nektar::NekMesh; + +namespace Nektar::NekMesh +{ +ModuleKey ProcessRevolve::className = + GetModuleFactory().RegisterCreatorFunction( + ModuleKey(eProcessModule, "revolve"), ProcessRevolve::create); + +ProcessRevolve::ProcessRevolve(MeshSharedPtr m) : ProcessModule(m) +{ + m_config["layers"] = ConfigOption(false, "24", "Number of layers"); + m_config["angle"] = + ConfigOption(false, "full", "Angle of revolution (in radians)"); +} + +ProcessRevolve::~ProcessRevolve() +{ +} + +void ProcessRevolve::Process() +{ + m_log(VERBOSE) << "Revolving grid." << endl; + + if (m_mesh->m_spaceDim != 2) + { + m_log(FATAL) << "Revolve should only be called for a two dimensional " + << "mesh" << endl; + } + + int nLayers = m_config["layers"].as<int>(); + // Special case needed for full revolutions + bool full = m_config["angle"].as<std::string>() == "full"; + NekDouble angle = full ? 2 * M_PI : m_config["angle"].as<NekDouble>(); + if (angle < 0 || angle > 2 * M_PI) + { + m_log(FATAL) << "Angle should be between 0 and 2pi" << endl; + } + NekDouble dphi = angle / nLayers; + // Increment space and expansion dimensions. + m_mesh->m_spaceDim++; + m_mesh->m_expDim++; + + // Grab a copy of the existing two-dimensional elements. + vector<ElementSharedPtr> el = m_mesh->m_element[2]; + + // Grab a copy of existing composites. + CompositeMap oldComp = m_mesh->m_composite; + m_log(VERBOSE) << "Boundary composites" << endl; + for (auto &it : oldComp) + { + if (it.second->m_tag != "E") + { + continue; + } + m_log(VERBOSE) << it.first << "\t" << it.second->m_tag; + for (int i = 0; i < it.second->m_items.size(); ++i) + { + m_log(VERBOSE) << "\t" << it.second->m_items[i]->GetId() << " (" + << it.second->m_items[i]->GetVertex(0) << ", " + << it.second->m_items[i]->GetVertex(1) << ")"; + vector<NodeSharedPtr> vv = it.second->m_items[i]->GetVertexList(); + m_log(VERBOSE) << "\t(" << vv[0]->GetID() << ", " << vv[1]->GetID() + << ")"; + } + m_log(VERBOSE) << endl; + } + + // Reset mesh. + for (int d = 0; d <= 3; ++d) + { + m_mesh->m_element[d].clear(); + } + + NodeSet nodes = m_mesh->m_vertexSet; + + map<int, NodeSharedPtr> id2node; + + for (auto &n : nodes) + { + id2node[n->m_id] = n; + } + // Save z plane coordinate + NekDouble phi0 = 0; + + // Create vertices for subsequent layers. + for (int i = 1; i < nLayers + 1; ++i) + { + if (full && i == nLayers) + { + // For last layer we will connect up to first layer + break; + } + for (auto &n : nodes) + { + if (n->m_x <= 0) + { + // Disallow x=0 for now + m_log(FATAL) + << "All vertices require positive (non-zero) x coordinates" + << endl; + } + NekDouble r0 = sqrt(n->m_x * n->m_x + n->m_z * n->m_z); + NodeSharedPtr newNode(new Node(i * nodes.size() + n->m_id, + r0 * cos(i * dphi), n->m_y, + r0 * sin(i * dphi))); + + m_mesh->m_vertexSet.insert(newNode); + id2node[i * nodes.size() + n->m_id] = newNode; + } + } + + EdgeSet esOld = m_mesh->m_edgeSet; // copy edges for curvature + + for (int j = 0; j < nLayers; ++j) + { + int next = full ? (j + 1) % nLayers : j + 1; + for (int i = 0; i < el.size(); ++i) + { + vector<NodeSharedPtr> verts = el[i]->GetVertexList(); + if (verts.size() == 4) + { + vector<NodeSharedPtr> nodeList(8); + + nodeList[0] = id2node[verts[0]->m_id + j * nodes.size()]; + nodeList[1] = id2node[verts[1]->m_id + j * nodes.size()]; + nodeList[2] = id2node[verts[2]->m_id + j * nodes.size()]; + nodeList[3] = id2node[verts[3]->m_id + j * nodes.size()]; + nodeList[4] = id2node[verts[0]->m_id + next * nodes.size()]; + nodeList[5] = id2node[verts[1]->m_id + next * nodes.size()]; + nodeList[6] = id2node[verts[2]->m_id + next * nodes.size()]; + nodeList[7] = id2node[verts[3]->m_id + next * nodes.size()]; + + vector<int> tags(1); + tags[0] = 0; + + ElmtConfig conf(LibUtilities::eHexahedron, 1, false, false, + false); + ElementSharedPtr E = GetElementFactory().CreateInstance( + LibUtilities::eHexahedron, conf, nodeList, tags); + + m_mesh->m_element[3].push_back(E); + } + else + { + vector<NodeSharedPtr> nodeList(6); + nodeList[0] = id2node[verts[0]->m_id + next * nodes.size()]; + nodeList[1] = id2node[verts[1]->m_id + next * nodes.size()]; + nodeList[2] = id2node[verts[1]->m_id + j * nodes.size()]; + nodeList[3] = id2node[verts[0]->m_id + j * nodes.size()]; + nodeList[4] = id2node[verts[2]->m_id + next * nodes.size()]; + nodeList[5] = id2node[verts[2]->m_id + j * nodes.size()]; + + vector<int> tags(1); + tags[0] = 1; + + ElmtConfig conf(LibUtilities::ePrism, 1, false, false, false); + ElementSharedPtr E = GetElementFactory().CreateInstance( + LibUtilities::ePrism, conf, nodeList, tags); + + m_mesh->m_element[3].push_back(E); + } + } + } + + ProcessEdges(); + ProcessFaces(); + ProcessElements(); + ProcessComposites(); + + // Copy edge information + for (auto &edge : esOld) + { + if (edge->m_edgeNodes.size() > 0) + { + for (int j = 0; j < nLayers + 1; ++j) + { + if (full && j == nLayers) + { + // For last layer we will connect up to first layer + break; + } + vector<NodeSharedPtr> ns(edge->m_edgeNodes.size()); + for (int i = 0; i < ns.size(); i++) + { + NodeSharedPtr n = edge->m_edgeNodes[i]; + NekDouble r0 = sqrt(n->m_x * n->m_x + n->m_z * n->m_z); + + ns[i] = std::shared_ptr<Node>(new Node( + 0, r0 * cos(j * dphi), n->m_y, r0 * sin(j * dphi))); + } + + EdgeSharedPtr e = std::shared_ptr<Edge>( + new Edge(id2node[edge->m_n1->m_id + j * nodes.size()], + id2node[edge->m_n2->m_id + j * nodes.size()])); + + auto f = m_mesh->m_edgeSet.find(e); + ASSERTL1(f != m_mesh->m_edgeSet.end(), "could not find edge"); + + // Copy edge type + (*f)->m_curveType = edge->m_curveType; + // Copy points + if ((*f)->m_n1 == e->m_n1) + { + (*f)->m_edgeNodes = ns; + } + else + { + reverse(ns.begin(), ns.end()); + (*f)->m_edgeNodes = ns; + } + } + } + } + + // Get composites max id + unsigned int maxCompId = 0; + for (auto &it : oldComp) + { + if (it.second->m_id >= maxCompId) + { + maxCompId = it.second->m_id; + } + } + + // First rename surface to volume composites to out of range + int outCompId = maxCompId + 1; + std::vector<int> toErase; + for (auto &it2 : m_mesh->m_composite) + { + if (it2.second->m_id > maxCompId) + { + // done! + break; + } + if (it2.second->m_tag == "H" || it2.second->m_tag == "R") + { + it2.second->m_id = outCompId; + m_mesh->m_composite.insert(std::make_pair(outCompId, it2.second)); + toErase.push_back(it2.first); + outCompId += 1; + } + } + + for (auto &e : toErase) + { + m_mesh->m_composite.erase(e); + } + + toErase.clear(); + + // Then copy surface to volume composites names + for (auto &it2 : m_mesh->m_composite) + { + if (it2.second->m_tag == "H" || it2.second->m_tag == "R") + { + for (auto &it1 : oldComp) + { + if (it2.second->m_tag == "H" && it1.second->m_tag == "Q") + { + it2.second->m_id = it1.second->m_id; + m_mesh->m_composite.insert( + std::make_pair(it1.second->m_id, it2.second)); + toErase.push_back(it2.first); + oldComp.erase(it1.first); + break; + } + else if (it2.second->m_tag == "R" && it1.second->m_tag == "T") + { + it2.second->m_id = it1.second->m_id; + m_mesh->m_composite.insert( + std::make_pair(it1.second->m_id, it2.second)); + toErase.push_back(it2.first); + oldComp.erase(it1.first); + break; + } + } + } + } + + for (auto &e : toErase) + { + m_mesh->m_composite.erase(e); + } + + // Add new composite to be filled with all boundary faces + CompositeSharedPtr comp(new Composite()); + comp->m_id = ++maxCompId; + unsigned int compAllFaceId = + maxCompId; // save it so we can remove it later on + comp->m_tag = "F"; + m_mesh->m_composite.insert(std::make_pair(maxCompId, comp)); + + // Add all boundary faces to the composite + auto allFaceC = m_mesh->m_composite.find(maxCompId); + m_log(VERBOSE) << "Faces boundary list" << endl; + + for (auto &it : m_mesh->m_faceSet) + { + // Add to composite if boundary face + if (it->m_elLink.size() < 2) + { + if (it->m_vertexList.size() == 3) + { + // Triangle + ElmtConfig conf(LibUtilities::eTriangle, 1, false, false); + vector<int> tags(1); + tags[0] = 1; + ElementSharedPtr E = GetElementFactory().CreateInstance( + LibUtilities::eTriangle, conf, it->m_vertexList, tags); + E->SetId(it->m_id); + allFaceC->second->m_items.push_back(E); + } + else if (it->m_vertexList.size() == 4) + { + // Quad + ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false); + vector<int> tags(1); + tags[0] = 0; + ElementSharedPtr E = GetElementFactory().CreateInstance( + LibUtilities::eQuadrilateral, conf, it->m_vertexList, tags); + E->SetId(it->m_id); + allFaceC->second->m_items.push_back(E); + } + } + } + + // Create boundary composites + for (auto &itOc : oldComp) + { + CompositeSharedPtr comp(new Composite()); + comp->m_id = itOc.second->m_id; + comp->m_tag = "F"; + m_mesh->m_composite.insert(std::make_pair(itOc.second->m_id, comp)); + } + // Create periodic composites + if (!full) + { + for (int i = 0; i < 2; i++) + { + CompositeSharedPtr comp(new Composite()); + comp->m_id = ++maxCompId; + comp->m_tag = "F"; + m_mesh->m_composite.insert(std::make_pair(maxCompId, comp)); + } + } + + // Populates boundary composites + for (auto &itQ : allFaceC->second->m_items) + { + // Check if this quad belongs to previous boundary + for (auto &itOc : oldComp) + { + for (int iEd = 0; iEd < itOc.second->m_items.size(); ++iEd) + { + int inCommon = 0; + for (int iV = 0; iV < itQ->GetVertexList().size(); iV++) + { + for (int j = 0; j < 2; j++) + { + NekDouble oldr2 = + itOc.second->m_items[iEd]->GetVertex(j)->m_x * + itOc.second->m_items[iEd]->GetVertex(j)->m_x + + itOc.second->m_items[iEd]->GetVertex(j)->m_z * + itOc.second->m_items[iEd]->GetVertex(j)->m_z; + NekDouble newr2 = + itQ->GetVertex(iV)->m_x * itQ->GetVertex(iV)->m_x + + itQ->GetVertex(iV)->m_z * itQ->GetVertex(iV)->m_z; + + if (LibUtilities::IsRealEqual(oldr2, newr2) && + LibUtilities::IsRealEqual( + itQ->GetVertex(iV)->m_y, + itOc.second->m_items[iEd]->GetVertex(j)->m_y)) + { + ++inCommon; + } + } + } + // If the face contains 4 xy pairs in common with 1 edge it + // must be an extruded edge and it should be added to the + // corresponding composite + if (inCommon == 4) + { + auto newC = m_mesh->m_composite.find(itOc.second->m_id); + // Quad + ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, + false); + vector<int> tags(1); + tags[0] = 0; + ElementSharedPtr E = GetElementFactory().CreateInstance( + LibUtilities::eQuadrilateral, conf, + itQ->GetVertexList(), tags); + E->SetId(itQ->GetId()); + newC->second->m_items.push_back(E); + } + } + } + + // Populates periodic composites + if (!full) + { + NekDouble phidist = 0.0; + for (int iV = 0; iV < itQ->GetVertexList().size(); iV++) + { + phidist += + (atan2(-itQ->GetVertex(iV)->m_z, -itQ->GetVertex(iV)->m_x) + + M_PI); + } + phidist = phidist / itQ->GetVertexList().size(); + unsigned int compPerId = 0; + if (LibUtilities::IsRealEqual(phidist, phi0)) + { + compPerId = maxCompId - 1; + } + else if (LibUtilities::IsRealEqual(phidist - phi0, angle)) + { + compPerId = maxCompId; + } + if (compPerId > 0 && itQ->GetVertexList().size() == 3) + { + // Triangle + auto perC = m_mesh->m_composite.find(compPerId); + ElmtConfig conf(LibUtilities::eTriangle, 1, false, false); + vector<int> tags(1); + tags[0] = 1; + ElementSharedPtr E = GetElementFactory().CreateInstance( + LibUtilities::eTriangle, conf, itQ->GetVertexList(), tags); + E->SetId(itQ->GetId()); + perC->second->m_items.push_back(E); + } + else if (compPerId > 0 && itQ->GetVertexList().size() == 4) + { + // Quad + auto perC = m_mesh->m_composite.find(compPerId); + ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false); + vector<int> tags(1); + tags[0] = 0; + ElementSharedPtr E = GetElementFactory().CreateInstance( + LibUtilities::eQuadrilateral, conf, itQ->GetVertexList(), + tags); + E->SetId(itQ->GetId()); + perC->second->m_items.push_back(E); + } + } + } + // Remove all faces composite + m_mesh->m_composite.erase(compAllFaceId); +} +} // namespace Nektar::NekMesh diff --git a/library/NekMesh/Module/ProcessModules/ProcessRevolve.h b/library/NekMesh/Module/ProcessModules/ProcessRevolve.h new file mode 100644 index 0000000000..e8a8fe0bd5 --- /dev/null +++ b/library/NekMesh/Module/ProcessModules/ProcessRevolve.h @@ -0,0 +1,68 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// File: ProcessRevolve.h +// +// For more information, please see: http://www.nektar.info/ +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: Revolve a two-dimensional mesh around the y-axis to form a +// three-dimensional mesh. Arguments are "angle" and "layers" +// +//////////////////////////////////////////////////////////////////////////////// + +#ifndef UTILITIES_NEKMESH_PROCESSREVOLVE +#define UTILITIES_NEKMESH_PROCESSREVOLVE + +#include <NekMesh/Module/Module.h> + +namespace Nektar::NekMesh +{ +/** + * @brief This processing module revolves a 2d mesh around the y-axis + */ +class ProcessRevolve : public NekMesh::ProcessModule +{ +public: + /// Creates an instance of this class + static std::shared_ptr<Module> create(NekMesh::MeshSharedPtr m) + { + return MemoryManager<ProcessRevolve>::AllocateSharedPtr(m); + } + static NekMesh::ModuleKey className; + + ProcessRevolve(NekMesh::MeshSharedPtr m); + ~ProcessRevolve() override; + + void Process() override; + + std::string GetModuleName() override + { + return "ProcessRevolve"; + } +}; +} // namespace Nektar::NekMesh + +#endif diff --git a/utilities/NekMesh/CMakeLists.txt b/utilities/NekMesh/CMakeLists.txt index 0115a493e9..748fc89102 100644 --- a/utilities/NekMesh/CMakeLists.txt +++ b/utilities/NekMesh/CMakeLists.txt @@ -11,6 +11,7 @@ ADD_NEKTAR_TEST(Nektar++/extract_curved_edge) ADD_NEKTAR_TEST(Nektar++/extract_curved_face) ADD_NEKTAR_TEST(Nektar++/extract_detectbnd_tube) ADD_NEKTAR_TEST(Nektar++/extrude) +ADD_NEKTAR_TEST(Nektar++/revolve) ADD_NEKTAR_TEST(Nektar++/jac_list_tet_face) ADD_NEKTAR_TEST(Nektar++/linearise_invalid_quad) ADD_NEKTAR_TEST(Nektar++/linearise_invalid_tet) diff --git a/utilities/NekMesh/Tests/Nektar++/revolve.tst b/utilities/NekMesh/Tests/Nektar++/revolve.tst new file mode 100644 index 0000000000..7c3a514b2f --- /dev/null +++ b/utilities/NekMesh/Tests/Nektar++/revolve.tst @@ -0,0 +1,58 @@ +<?xml version="1.0" encoding="utf-8" ?> +<test> + <description>Revolution of 2D mesh</description> + <executable>NekMesh</executable> + <parameters>-m jac:list -m revolve:layers=1:angle=0.1 revolve.xml + revolve-out.xml:xml:test:stats</parameters> + <files> + <file description="Input File">revolve.xml</file> + </files> + <metrics> + <metric type="regex" id="1"> + <regex>.*Total negative Jacobians\s*:\s*(\d+)</regex> + <matches> + <match> + <field id="0">0</field> + </match> + </matches> + </metric> + <metric type="regex" id="2"> + <regex>^.*Node count\s*:\s*(\d+)</regex> + <matches> + <match> + <field id="0">12</field> + </match> + </matches> + </metric> + <metric type="regex" id="3"> + <regex>^.*Number of composites\s*:\s*(\d+)</regex> + <matches> + <match> + <field id="0">8</field> + </match> + </matches> + </metric> + <metric type="regex" id="4"> + <regex>^.*Lower mesh extent\s*:\s*(-?\d+\.?\d*)\s*(-?\d+\.?\d*)\s* + (-?\d+\.?\d*)</regex> + <matches> + <match> + <field id="0">0.995004</field> + <field id="1">-0.1</field> + <field id="2">0</field> + </match> + </matches> + </metric> + <metric type="regex" id="5"> + <regex>^.*Upper mesh extent\s*:\s*(-?\d+\.?\d*)\s*(-?\d+\.?\d*)\s* + (-?\d+\.?\d*)</regex> + <matches> + <match> + <field id="0">2</field> + <field id="1">1</field> + <field id="2">0.199667</field> + </match> + </matches> + </metric> + </metrics> +</test> diff --git a/utilities/NekMesh/Tests/Nektar++/revolve.xml b/utilities/NekMesh/Tests/Nektar++/revolve.xml new file mode 100644 index 0000000000..965088844e --- /dev/null +++ b/utilities/NekMesh/Tests/Nektar++/revolve.xml @@ -0,0 +1,54 @@ +<?xml version="1.0" encoding="utf-8" ?> +<NEKTAR> + <GEOMETRY DIM="2" SPACE="2"> + <VERTEX> + <V ID="0">2.00000000e+00 0.50000000e+00 0.00000000e+00</V> + <V ID="1">1.00000000e+00 0.00000000e+00 0.00000000e+00</V> + <V ID="2">2.00000000e+00 0.00000000e+00 0.00000000e+00</V> + <V ID="3">1.00000000e+00 0.50000000e+00 0.00000000e+00</V> + <V ID="4">2.00000000e+00 1.00000000e+00 0.00000000e+00</V> + <V ID="5">1.00000000e+00 1.00000000e+00 0.00000000e+00</V> + </VERTEX> + <EDGE> + <E ID="0">2 1</E> + <E ID="1">0 2</E> + <E ID="2">3 0</E> + <E ID="3">1 3</E> + <E ID="4">0 5</E> + <E ID="5">3 5</E> + <E ID="6">0 4</E> + <E ID="7">4 5</E> + </EDGE> + <ELEMENT> + <Q ID="0">0 1 2 3</Q> + <T ID="1">2 4 5</T> + <T ID="2">4 6 7</T> + </ELEMENT> + <CURVED> + <E ID="0" EDGEID="0" NUMPOINTS="5" TYPE="GaussLobattoLegendre">2.00000000e+00 0.00000000e+00 0.00000000e+00 1.87326835e+00 -5.71428571e-02 0.00000000e+00 1.500000000e+00 -1.00000000e-01 0.00000000e+00 1.172673165e+00 -5.71428571e-02 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 </E> + </CURVED> + <COMPOSITE> + <C ID="1"> E[7] </C> + <C ID="2"> E[1,6] </C> + <C ID="3"> E[0] </C> + <C ID="4"> E[3,5] </C> + <C ID="100"> Q[0] </C> + <C ID="111"> T[1,2] </C> + </COMPOSITE> + <DOMAIN> C[100,111] </DOMAIN> + </GEOMETRY> + <Metadata> + <Provenance> + <GitBranch>refs/heads/giacomo-current</GitBranch> + <GitSHA1>4e9c52d07baedb4dba5d971a9ca604f912ac85c4</GitSHA1> + <Hostname>gcastigl</Hostname> + <NektarVersion>4.5.0</NektarVersion> + <Timestamp>31-Jan-2019 12:30:52</Timestamp> + </Provenance> + <NekMeshCommandLine>square_curved.xml square_curved.xml:xml:uncompress </NekMeshCommandLine> + </Metadata> + <EXPANSIONS> + <E COMPOSITE="C[100]" NUMMODES="4" TYPE="MODIFIED" FIELDS="u" /> + <E COMPOSITE="C[111]" NUMMODES="4" TYPE="MODIFIED" FIELDS="u" /> + </EXPANSIONS> +</NEKTAR> -- GitLab