From 653a5e9b7aec49ea30e263f171cb91c370401f98 Mon Sep 17 00:00:00 2001 From: Matt McGurn <75279702+mmcgurn@users.noreply.github.com> Date: Tue, 30 Jan 2024 13:10:33 -0500 Subject: [PATCH] axisymmetric cylinder mesh (#513) * working on axissymmetric * refactor for mesh generator and desciption * working on mesh * working very simple cylinder I think * restored chrest compiler flags * added some tests for simple case * working protoype code * added radius support * added example and updated tests * version bump * working on boundary labeling * problem with extrude * updated tests * fixed hex element order * fixed hex element order * fixed tri prism element order * updated tests for new axisymmetric. * updated tests * updated tests * typo fix --- CMakeLists.txt | 2 +- src/domain/CMakeLists.txt | 5 +- src/domain/descriptions/CMakeLists.txt | 8 + src/domain/descriptions/axisymmetric.cpp | 170 +++++++++++++ src/domain/descriptions/axisymmetric.hpp | 132 ++++++++++ src/domain/descriptions/meshDescription.hpp | 63 +++++ src/domain/meshGenerator.cpp | 226 +++++++++++++++++ src/domain/meshGenerator.hpp | 51 ++++ src/domain/modifiers/extrudeLabel.cpp | 55 +++- src/domain/range.cpp | 6 +- src/domain/range.hpp | 10 +- src/domain/region.cpp | 10 +- src/domain/region.hpp | 39 ++- src/utilities/vectorUtilities.hpp | 53 +++- .../axisymmetricPipeFlow.txt | 13 + .../axisymmetricPipeFlow.yaml | 190 ++++++++++++++ .../bottomBoundary_monitor.xmf | 0 .../boundaryMonitorTest2D.txt | 0 .../boundaryMonitorTest2D.yaml | 4 +- .../bottomBoundary_monitor.xmf | 0 .../boundaryMonitorTest3D.txt | 0 .../boundaryMonitorTest3D.yaml | 4 +- .../dmViewFromOptions.txt | 0 .../dmViewFromOptions.yaml | 2 +- .../extrudeBoundaryTest.txt | 0 .../extrudeBoundaryTest.yaml | 2 +- .../domain/meshGeneratorAxisymmetric.txt | 18 ++ .../domain/meshGeneratorAxisymmetric.yaml | 53 ++++ .../{machinery => domain}/meshMappingTest.txt | 0 .../meshMappingTest.yaml | 2 +- .../meshMappingTestCoordinateSpace.txt | 0 .../meshMappingTestCoordinateSpace.yaml | 2 +- .../mixedCellTypeTest/domain.xmf | 0 .../mixedCellTypeTest2D.yaml | 2 +- .../mixedCellTypeTest/mixedCells2D.msh | 0 .../{machinery => domain}/subDomainFVM.yaml | 4 +- .../subDomainFVM/fluidField.xmf | 0 .../subDomainFVM/subDomainFVM.txt | 0 tests/unitTests/domain/CMakeLists.txt | 1 + .../domain/descriptions/CMakeLists.txt | 4 + .../domain/descriptions/axisymmetricTests.cpp | 239 ++++++++++++++++++ tests/unitTests/utilities/CMakeLists.txt | 1 + .../utilities/vectorUtilitiesTests.cpp | 40 +++ 43 files changed, 1368 insertions(+), 43 deletions(-) create mode 100644 src/domain/descriptions/CMakeLists.txt create mode 100644 src/domain/descriptions/axisymmetric.cpp create mode 100644 src/domain/descriptions/axisymmetric.hpp create mode 100644 src/domain/descriptions/meshDescription.hpp create mode 100644 src/domain/meshGenerator.cpp create mode 100644 src/domain/meshGenerator.hpp create mode 100644 tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.txt create mode 100644 tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.yaml rename tests/integrationTests/inputs/{machinery => domain}/boundaryMonitorTest2D/bottomBoundary_monitor.xmf (100%) rename tests/integrationTests/inputs/{machinery => domain}/boundaryMonitorTest2D/boundaryMonitorTest2D.txt (100%) rename tests/integrationTests/inputs/{machinery => domain}/boundaryMonitorTest2D/boundaryMonitorTest2D.yaml (96%) rename tests/integrationTests/inputs/{machinery => domain}/boundaryMonitorTest3D/bottomBoundary_monitor.xmf (100%) rename tests/integrationTests/inputs/{machinery => domain}/boundaryMonitorTest3D/boundaryMonitorTest3D.txt (100%) rename tests/integrationTests/inputs/{machinery => domain}/boundaryMonitorTest3D/boundaryMonitorTest3D.yaml (96%) rename tests/integrationTests/inputs/{machinery => domain}/dmViewFromOptions.txt (100%) rename tests/integrationTests/inputs/{machinery => domain}/dmViewFromOptions.yaml (97%) rename tests/integrationTests/inputs/{machinery => domain}/extrudeBoundaryTest.txt (100%) rename tests/integrationTests/inputs/{machinery => domain}/extrudeBoundaryTest.yaml (98%) create mode 100644 tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.txt create mode 100644 tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.yaml rename tests/integrationTests/inputs/{machinery => domain}/meshMappingTest.txt (100%) rename tests/integrationTests/inputs/{machinery => domain}/meshMappingTest.yaml (97%) rename tests/integrationTests/inputs/{machinery => domain}/meshMappingTestCoordinateSpace.txt (100%) rename tests/integrationTests/inputs/{machinery => domain}/meshMappingTestCoordinateSpace.yaml (96%) rename tests/integrationTests/inputs/{machinery => domain}/mixedCellTypeTest/domain.xmf (100%) rename tests/integrationTests/inputs/{machinery => domain}/mixedCellTypeTest/mixedCellTypeTest2D.yaml (97%) rename tests/integrationTests/inputs/{machinery => domain}/mixedCellTypeTest/mixedCells2D.msh (100%) rename tests/integrationTests/inputs/{machinery => domain}/subDomainFVM.yaml (97%) rename tests/integrationTests/inputs/{machinery => domain}/subDomainFVM/fluidField.xmf (100%) rename tests/integrationTests/inputs/{machinery => domain}/subDomainFVM/subDomainFVM.txt (100%) create mode 100644 tests/unitTests/domain/descriptions/CMakeLists.txt create mode 100644 tests/unitTests/domain/descriptions/axisymmetricTests.cpp create mode 100644 tests/unitTests/utilities/vectorUtilitiesTests.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f1fc65305..1c5749937 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.18.4) include(config/petscCompilers.cmake) # Set the project details -project(ablateLibrary VERSION 0.12.21) +project(ablateLibrary VERSION 0.12.22) # Load the Required 3rd Party Libaries pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc) diff --git a/src/domain/CMakeLists.txt b/src/domain/CMakeLists.txt index 0a3ffd926..143135f25 100644 --- a/src/domain/CMakeLists.txt +++ b/src/domain/CMakeLists.txt @@ -17,6 +17,7 @@ target_sources(ablateLibrary initializer.cpp hdf5Initializer.cpp initializerList.cpp + meshGenerator.cpp PUBLIC domain.hpp @@ -39,7 +40,9 @@ target_sources(ablateLibrary initializer.hpp hdf5Initializer.hpp initializerList.hpp - ) + meshGenerator.hpp +) add_subdirectory(modifiers) add_subdirectory(RBF) +add_subdirectory(descriptions) diff --git a/src/domain/descriptions/CMakeLists.txt b/src/domain/descriptions/CMakeLists.txt new file mode 100644 index 000000000..614978a06 --- /dev/null +++ b/src/domain/descriptions/CMakeLists.txt @@ -0,0 +1,8 @@ +target_sources(ablateLibrary + PRIVATE + axisymmetric.cpp + + PUBLIC + meshDescription.hpp + axisymmetric.hpp + ) diff --git a/src/domain/descriptions/axisymmetric.cpp b/src/domain/descriptions/axisymmetric.cpp new file mode 100644 index 000000000..09d5459ee --- /dev/null +++ b/src/domain/descriptions/axisymmetric.cpp @@ -0,0 +1,170 @@ +#include "axisymmetric.hpp" + +#include +#include "utilities/vectorUtilities.hpp" + +ablate::domain::descriptions::Axisymmetric::Axisymmetric(const std::vector &startLocation, PetscReal length, std::shared_ptr radiusFunction, + PetscInt numberWedges, PetscInt numberSlices, PetscInt numberShells) + : startLocation(utilities::VectorUtilities::ToArray(startLocation)), + length(length), + radiusFunction(std::move(radiusFunction)), + numberWedges(numberWedges), + numberSlices(numberSlices), + numberShells(numberShells), + numberCellsPerSlice(numberWedges * numberShells), + numberCellsPerShell(numberWedges * numberSlices), + numberVerticesPerShell(numberWedges * (numberSlices + 1)), + numberCenterVertices(numberSlices + 1), + numberCells(numberCellsPerSlice * numberSlices), + numberVertices(numberVerticesPerShell * numberShells + numberCenterVertices), + numberTriPrismCells(numberWedges * numberSlices) { + // make sure there are at least 4 wedges and one slice + if (numberWedges < 3 || numberSlices < 1 || numberShells < 1) { + throw std::invalid_argument("Axisymmetric requires at least 3 wedges, 1 slice, and 1 shell."); + } +} + +void ablate::domain::descriptions::Axisymmetric::BuildTopology(PetscInt cell, PetscInt *cellNodes) const { + // Note that the cell/vertex ordering grows in shells so that tri prism elements are specified first and together. + // flip the cell order so that hexes appear before tri prisms + cell = CellReverser(cell); + + // determine the type of element + if (cell >= numberTriPrismCells) { + // this is a hex cell + // determine which shell this is + PetscInt cellShell = cell / numberCellsPerShell; + + // determine which slice this cell falls upon + PetscInt cellSlice = (cell - cellShell * numberCellsPerShell) / numberWedges; + + // determine the local cell index (0 - numberWedges) + PetscInt cellIndex = cell % numberWedges; + + PetscInt lowerSliceLowerShellOffset = (cellShell - 1) * numberVerticesPerShell + numberCenterVertices + cellSlice * numberWedges; + PetscInt lowerSliceUpperShellOffset = cellShell * numberVerticesPerShell + numberCenterVertices + cellSlice * numberWedges; + PetscInt upperSliceLowerShellOffset = (cellShell - 1) * numberVerticesPerShell + numberCenterVertices + (cellSlice + 1) * numberWedges; + PetscInt upperSliceUpperShellOffset = cellShell * numberVerticesPerShell + numberCenterVertices + (cellSlice + 1) * numberWedges; + + cellNodes[0] = upperSliceLowerShellOffset + cellIndex; + cellNodes[1] = lowerSliceLowerShellOffset + cellIndex; + cellNodes[2] = lowerSliceUpperShellOffset + cellIndex; + cellNodes[3] = upperSliceUpperShellOffset + cellIndex; + + cellNodes[4] = (cellIndex + 1 == numberWedges) ? upperSliceLowerShellOffset : upperSliceLowerShellOffset + cellIndex + 1 /*check for wrap around*/; + cellNodes[5] = (cellIndex + 1 == numberWedges) ? upperSliceUpperShellOffset : upperSliceUpperShellOffset + cellIndex + 1 /*check for wrap around*/; + cellNodes[6] = (cellIndex + 1 == numberWedges) ? lowerSliceUpperShellOffset : lowerSliceUpperShellOffset + cellIndex + 1 /*check for wrap around*/; + cellNodes[7] = (cellIndex + 1 == numberWedges) ? lowerSliceLowerShellOffset : lowerSliceLowerShellOffset + cellIndex + 1 /*check for wrap around*/; + } else { + // This is a tri prism + // determine which slice this is + PetscInt slice = cell / numberWedges; + + // determine the local cell/wedge number (in this slice) + PetscInt localWedge = cell % numberWedges; + + // Set the center nodes + const auto lowerCenter = slice; + const auto upperCenter = (slice + 1); + + cellNodes[0] = lowerCenter; + cellNodes[3] = upperCenter; + + // Compute the offset for this slice + const auto nodeSliceOffset = numberCenterVertices + numberWedges * slice; + const auto nextNodeSliceOffset = numberCenterVertices + numberWedges * (slice + 1); + + // set the lower nodes + cellNodes[2] = localWedge + nodeSliceOffset; + cellNodes[1] = (localWedge + 1) == numberWedges ? nodeSliceOffset : localWedge + 1 + nodeSliceOffset; // checking for wrap around + + // repeat for the upper nodes + cellNodes[5] = (localWedge + 1) == numberWedges ? nextNodeSliceOffset : localWedge + 1 + nextNodeSliceOffset; // checking for wrap around + cellNodes[4] = localWedge + nextNodeSliceOffset; + } +} +void ablate::domain::descriptions::Axisymmetric::SetCoordinate(PetscInt node, PetscReal *coordinate) const { + // start the coordinate at the start location + coordinate[0] = startLocation[0]; + coordinate[1] = startLocation[1]; + coordinate[2] = startLocation[2]; + + // determine where this node is, there is a special case for the 0 node shell or center + PetscInt nodeShell, nodeSlice, nodeRotationIndex; + if (node < numberCenterVertices) { + nodeShell = 0; + nodeSlice = node; + nodeRotationIndex = 0; + } else { + nodeShell = ((node - numberCenterVertices) / numberVerticesPerShell) + 1; + nodeSlice = (node - numberCenterVertices - (nodeShell - 1) * numberVerticesPerShell) / numberWedges; + nodeRotationIndex = (node - numberCenterVertices) % numberWedges; + } + + // offset along z + auto delta = length / numberSlices; + auto offset = delta * nodeSlice; + coordinate[2] += offset; + + // compute the maximum radius at this coordinate + auto radius = radiusFunction->Eval(coordinate, 3, NAN); + + // if we are not at the center + auto radiusFactor = ((PetscReal)nodeShell) / ((PetscReal)numberShells); + coordinate[0] += radius * radiusFactor * PetscCosReal(2.0 * nodeRotationIndex * PETSC_PI / numberWedges); + coordinate[1] += radius * radiusFactor * PetscSinReal(2.0 * nodeRotationIndex * PETSC_PI / numberWedges); +} + +DMPolytopeType ablate::domain::descriptions::Axisymmetric::GetCellType(PetscInt cell) const { + // flip the cell order so that hexes appear before tri prisms + return CellReverser(cell) < numberTriPrismCells ? DM_POLYTOPE_TRI_PRISM : DM_POLYTOPE_HEXAHEDRON; +} + +std::shared_ptr ablate::domain::descriptions::Axisymmetric::GetRegion(const std::set &face) const { + // check to see if each node in on the surface + bool onOuterShell = true; + bool onLowerEndCap = true; + bool onUpperEndCap = true; + + // compute the outer shell start + auto outerShellStart = numberCenterVertices + numberVerticesPerShell * (numberShells - 1); + + for (const auto &node : face) { + // check if we are on the outer shell + if (node < outerShellStart) { + onOuterShell = false; + } + // check if we are on the end caps + PetscInt nodeSlice; + if (node < numberCenterVertices) { + nodeSlice = node; + } else { + auto nodeShell = ((node - numberCenterVertices) / numberVerticesPerShell) + 1; + nodeSlice = (node - numberCenterVertices - (nodeShell - 1) * numberVerticesPerShell) / numberWedges; + } + + if (nodeSlice != 0) { + onLowerEndCap = false; + } + if (nodeSlice != numberSlices) { + onUpperEndCap = false; + } + } + + // determine what region to return + if (onOuterShell) { + return shellBoundary; + } else if (onLowerEndCap) { + return lowerCapBoundary; + } else if (onUpperEndCap) { + return upperCapBoundary; + } + + return nullptr; +} + +#include "registrar.hpp" +REGISTER(ablate::domain::descriptions::MeshDescription, ablate::domain::descriptions::Axisymmetric, "The Axisymmetric MeshDescription is used to create an axisymmetric mesh around the z axis", + ARG(std::vector, "start", "the start coordinate of the mesh, must be 3D"), ARG(double, "length", "the length of the domain starting at the start coordinate"), + ARG(ablate::mathFunctions::MathFunction, "radius", "a radius function that describes the radius as a function of z"), ARG(int, "numberWedges", "wedges/pie slices in the circle"), + ARG(int, "numberSlices", "slicing of the cylinder along the z axis"), ARG(int, "numberShells", "slicing of the cylinder along the radius")); diff --git a/src/domain/descriptions/axisymmetric.hpp b/src/domain/descriptions/axisymmetric.hpp new file mode 100644 index 000000000..10bd8e825 --- /dev/null +++ b/src/domain/descriptions/axisymmetric.hpp @@ -0,0 +1,132 @@ +#ifndef ABLATELIBRARY_AXISYMMETRIC_HPP +#define ABLATELIBRARY_AXISYMMETRIC_HPP + +#include +#include +#include +#include "mathFunctions/mathFunction.hpp" +#include "meshDescription.hpp" + +namespace ablate::domain::descriptions { + +/** + * produces an axisymmetric mesh around the z axis + */ +class Axisymmetric : public ablate::domain::descriptions::MeshDescription { + private: + //! Store the start and end location of the mesh + const std::array startLocation; + const PetscReal length; // this is in z + + //! function used to describe a single return value (radius) as a functino of z + const std::shared_ptr radiusFunction; + + //! hard code the assumed mesh dimension + inline static const PetscInt dim = 3; + + //! Store the number of wedges, wedges/pie slices in the circle + const PetscInt numberWedges; + + //! Store the number of slices, slicing of the cylinder along the z axis + const PetscInt numberSlices; + + //! Store the number of shells, slicing of the cylinder along the radius + const PetscInt numberShells; + + //! Store the number of cells per slice + const PetscInt numberCellsPerSlice; + + //! store the number of cells per slice + const PetscInt numberCellsPerShell; + + //! Store the number of vertices per shell, does not include the center + const PetscInt numberVerticesPerShell; + + //! Store the number of vertices in the center + const PetscInt numberCenterVertices; + + //! Compute the number of cells + const PetscInt numberCells; + + //! And the number of vertices + const PetscInt numberVertices; + + //! store the number of tri prism cells to simplify the logic + const PetscInt numberTriPrismCells; + + //! precompute the region identifier for the boundary + const static inline std::shared_ptr shellBoundary = std::make_shared("outerShell"); + const static inline std::shared_ptr lowerCapBoundary = std::make_shared("lowerCap"); + const static inline std::shared_ptr upperCapBoundary = std::make_shared("upperCap"); + + /** + * Reverse the cell order to make sure hexes appear before tri prism + * @param in + * @return + */ + inline PetscInt CellReverser(PetscInt in) const { return numberCells - in - 1; } + + public: + /** + * generate and precompute a bunch of the required parameters + * @param startLocation the start coordinate of the mesh, must be 3D + * @param length the length of the domain starting at the start coordinate + * @param radiusFunction a radius function that describes the radius as a function of z + * @param numberWedges wedges/pie slices in the circle + * @param numberSlices slicing of the cylinder along the z axis + * @param numberShells slicing of the cylinder along the radius + */ + Axisymmetric(const std::vector& startLocation, PetscReal length, std::shared_ptr radiusFunction, PetscInt numberWedges, PetscInt numberSlices, + PetscInt numberShells); + + /** + * The overall assumed dimension of the mesh + * @return + */ + [[nodiscard]] const PetscInt& GetMeshDimension() const override { return dim; } + + /** + * Add some getters, so we can virtualize this in the future if needed + * @return + */ + [[nodiscard]] const PetscInt& GetNumberCells() const override { return numberCells; } + + /** + * Add some getters, so we can virtualize this in the future if needed + * @return + */ + [[nodiscard]] const PetscInt& GetNumberVertices() const override { return numberVertices; } + + /** + * Return the cell type for this cell, based off zero offset + * @return + */ + [[nodiscard]] DMPolytopeType GetCellType(PetscInt cell) const override; + + /** + * Builds the topology based upon a zero vertex offset. The cellNodes should be sized to hold cell + * @return + */ + void BuildTopology(PetscInt cell, PetscInt* cellNodes) const override; + + /** + * Builds the topology based upon a zero vertex offset. The cellNodes should be sized to hold cell + * @return + */ + void SetCoordinate(PetscInt node, PetscReal* coordinate) const override; + + /** + * Returns a hard coded boundary region name + * @return + */ + [[nodiscard]] std::shared_ptr GetBoundaryRegion() const override { return std::make_shared("boundary"); } + + /** + * returns the boundary region for the end caps and outer shell + * @param face + * @return + */ + [[nodiscard]] std::shared_ptr GetRegion(const std::set& face) const override; +}; +} // namespace ablate::domain::descriptions +#endif // ABLATELIBRARY_AXISYMMETRIC_HPP diff --git a/src/domain/descriptions/meshDescription.hpp b/src/domain/descriptions/meshDescription.hpp new file mode 100644 index 000000000..d8d3a3d7b --- /dev/null +++ b/src/domain/descriptions/meshDescription.hpp @@ -0,0 +1,63 @@ +#ifndef ABLATELIBRARY_MESHDESCRIPTION_HPP +#define ABLATELIBRARY_MESHDESCRIPTION_HPP + +#include +#include +#include + +namespace ablate::domain::descriptions { +/** + * A simple interface that is responsible for determine cell/vertex locations in the mesh + */ +class MeshDescription { + public: + virtual ~MeshDescription() = default; + /** + * The overall assumed dimension of the mesh + * @return + */ + [[nodiscard]] virtual const PetscInt& GetMeshDimension() const = 0; + + /** + * Total number of cells in the entire mesh + * @return + */ + [[nodiscard]] virtual const PetscInt& GetNumberCells() const = 0; + + /** + * Total number of nodes/vertices in the entire mesh + * @return + */ + [[nodiscard]] virtual const PetscInt& GetNumberVertices() const = 0; + + /** + * Return the cell type for this cell, based off zero offset + * @return + */ + [[nodiscard]] virtual DMPolytopeType GetCellType(PetscInt cell) const = 0; + + /** + * Builds the topology based upon a zero vertex offset. The cellNodes should be sized to hold cell + * @return + */ + virtual void BuildTopology(PetscInt cell, PetscInt* cellNodes) const = 0; + + /** + * Builds the node coordinate for each vertex + * @return + */ + virtual void SetCoordinate(PetscInt node, PetscReal* coordinate) const = 0; + + /** + * returns the boundary region to label for the entire region + */ + [[nodiscard]] virtual std::shared_ptr GetBoundaryRegion() const = 0; + + /** + * returns a region if the nodes set (on a face) belongs to a certain region + */ + [[nodiscard]] virtual std::shared_ptr GetRegion(const std::set& face) const = 0; +}; +} // namespace ablate::domain::descriptions + +#endif // ABLATELIBRARY_MESHDESCRIPTION_HPP diff --git a/src/domain/meshGenerator.cpp b/src/domain/meshGenerator.cpp new file mode 100644 index 000000000..de521b5a9 --- /dev/null +++ b/src/domain/meshGenerator.cpp @@ -0,0 +1,226 @@ +#include "meshGenerator.hpp" +#include +#include +#include "utilities/mpiUtilities.hpp" +#include "utilities/vectorUtilities.hpp" + +ablate::domain::MeshGenerator::MeshGenerator(const std::string& name, std::vector> fieldDescriptors, + const std::shared_ptr& description, std::vector> modifiers, + const std::shared_ptr& options) + : Domain(CreateDM(name, description), name, std::move(fieldDescriptors), std::move(modifiers), options) {} + +ablate::domain::MeshGenerator::~MeshGenerator() { + if (dm) { + DMDestroy(&dm); + } +} + +void ablate::domain::MeshGenerator::ReplaceDm(DM& originalDm, DM& replaceDm) { + if (replaceDm) { + // copy over the name + const char* name; + PetscObjectName((PetscObject)originalDm) >> utilities::PetscUtilities::checkError; + PetscObjectGetName((PetscObject)originalDm, &name) >> utilities::PetscUtilities::checkError; + PetscObjectSetName((PetscObject)replaceDm, name) >> utilities::PetscUtilities::checkError; + + // Copy over the options object + PetscOptions options; + PetscObjectGetOptions((PetscObject)originalDm, &options) >> utilities::PetscUtilities::checkError; + PetscObjectSetOptions((PetscObject)replaceDm, options) >> utilities::PetscUtilities::checkError; + ((DM_Plex*)(replaceDm)->data)->useHashLocation = ((DM_Plex*)originalDm->data)->useHashLocation; + + DMDestroy(&originalDm) >> utilities::PetscUtilities::checkError; + originalDm = replaceDm; + } +} + +DM ablate::domain::MeshGenerator::CreateDM(const std::string& name, const std::shared_ptr& description) { + // Create the new dm and set it to be a dmplex + DM dm; + DMCreate(PETSC_COMM_WORLD, &dm) >> utilities::MpiUtilities::checkError; + DMSetType(dm, DMPLEX) >> utilities::MpiUtilities::checkError; + + // Get the rank + PetscInt dim = description->GetMeshDimension(); + PetscMPIInt rank; + MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank) >> utilities::MpiUtilities::checkError; + DMSetDimension(dm, dim) >> utilities::PetscUtilities::checkError; + + /* Must create the celltype label here so that we do not automatically try to compute the types */ + DMCreateLabel(dm, "celltype") >> utilities::PetscUtilities::checkError; + + /* Create topology */ + // Get the number of cells and numVertices from the + auto numVertices = description->GetNumberVertices(); + auto numCells = description->GetNumberCells(); + + // set the values to zero if not on the first rank + numCells = rank == 0 ? numCells : 0; + numVertices = rank == 0 ? numVertices : 0; + DMPlexSetChart(dm, 0, numCells + numVertices) >> utilities::PetscUtilities::checkError; + + // Determine the max cone size and set the value for each cell. Cells come before points + PetscInt maxConeSize = 0; + for (PetscInt c = 0; c < numCells; ++c) { + auto cellType = description->GetCellType(c); + + // Determine the cone size and dim + auto coneSize = DMPolytopeTypeGetNumVertices(cellType); + maxConeSize = PetscMax(maxConeSize, coneSize); + + // Set the cone size + DMPlexSetConeSize(dm, c, coneSize) >> utilities::PetscUtilities::checkError; + } + + // With all the cones set, set up the dm + DMSetUp(dm) >> utilities::PetscUtilities::checkError; + + // Size up a buffer for the cone + PetscInt cone[maxConeSize]; + + // Compute and set the cone for each cell + for (PetscInt c = 0; c < numCells; ++c) { + description->BuildTopology(c, cone); + + // Offset the cone from the number of numVertices + for (auto& node : cone) { + node += numCells; + } + + auto cellType = description->GetCellType(c); + DMPlexSetCone(dm, c, cone) >> utilities::PetscUtilities::checkError; + DMPlexSetCellType(dm, c, cellType) >> utilities::PetscUtilities::checkError; + } + DMPlexSymmetrize(dm) >> utilities::PetscUtilities::checkError; + DMPlexStratify(dm) >> utilities::PetscUtilities::checkError; + + // Now compute the location for each of the cells + for (PetscInt v = numCells; v < numCells + numVertices; ++v) { + DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT) >> utilities::PetscUtilities::checkError; + } + /* Create cylinder geometry */ + { + Vec coordinates; + PetscSection coordSection; + PetscScalar* coords; + + /* Build coordinates */ + DMGetCoordinateSection(dm, &coordSection) >> utilities::PetscUtilities::checkError; + PetscSectionSetNumFields(coordSection, 1) >> utilities::PetscUtilities::checkError; + PetscSectionSetFieldComponents(coordSection, 0, dim) >> utilities::PetscUtilities::checkError; + PetscSectionSetChart(coordSection, numCells, numCells + numVertices) >> utilities::PetscUtilities::checkError; + for (PetscInt v = numCells; v < numCells + numVertices; ++v) { + PetscSectionSetDof(coordSection, v, dim) >> utilities::PetscUtilities::checkError; + PetscSectionSetFieldDof(coordSection, v, 0, dim) >> utilities::PetscUtilities::checkError; + } + PetscSectionSetUp(coordSection) >> utilities::PetscUtilities::checkError; + PetscInt coordSize; + PetscSectionGetStorageSize(coordSection, &coordSize) >> utilities::PetscUtilities::checkError; + VecCreate(PETSC_COMM_SELF, &coordinates) >> utilities::PetscUtilities::checkError; + PetscObjectSetName((PetscObject)coordinates, "coordinates") >> utilities::PetscUtilities::checkError; + VecSetSizes(coordinates, coordSize, PETSC_DETERMINE) >> utilities::PetscUtilities::checkError; + VecSetBlockSize(coordinates, dim) >> utilities::PetscUtilities::checkError; + VecSetType(coordinates, VECSTANDARD) >> utilities::PetscUtilities::checkError; + VecGetArray(coordinates, &coords) >> utilities::PetscUtilities::checkError; + + // March over each vertex + for (PetscInt v = 0; v < numVertices; ++v) { + // Get the pointer for this coordinate + auto vertex = coords + (v * dim); + + // Set it to zero for safety + for (PetscInt d = 0; d < dim; ++d) { + vertex[d] = 0.0; + } + + // Compute the value in the description + description->SetCoordinate(v, vertex); + } + + VecRestoreArray(coordinates, &coords) >> utilities::PetscUtilities::checkError; + DMSetCoordinatesLocal(dm, coordinates) >> utilities::PetscUtilities::checkError; + VecDestroy(&coordinates) >> utilities::PetscUtilities::checkError; + } + /* Interpolate */ + DM idm; + DMPlexInterpolate(dm, &idm) >> utilities::PetscUtilities::checkError; + DMPlexCopyCoordinates(dm, idm) >> utilities::PetscUtilities::checkError; + ReplaceDm(dm, idm); + + // set the name + PetscObjectSetName((PetscObject)dm, name.c_str()) >> utilities::PetscUtilities::checkError; + + // label of the boundaries + LabelBoundaries(description, dm); + + return dm; +} + +void ablate::domain::MeshGenerator::LabelBoundaries(const std::shared_ptr& description, DM& dm) { + // start by marking all the boundary faces + auto boundaryRegion = description->GetBoundaryRegion(); + DMLabel boundaryLabel; + PetscInt boundaryValue; + boundaryRegion->CreateLabel(dm, boundaryLabel, boundaryValue); + DMPlexMarkBoundaryFaces(dm, boundaryValue, boundaryLabel) >> utilities::PetscUtilities::checkError; + DMPlexLabelComplete(dm, boundaryLabel) >> utilities::PetscUtilities::checkError; + + // determine the new vertex/nodes bounds + PetscInt vStart, vEnd; + DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd) >> utilities::PetscUtilities::checkError; + ; // Range of vertices + + // determine the new face bounds + PetscInt fStart, fEnd; + DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd) >> utilities::PetscUtilities::checkError; + ; // Range of faces + + // now march over each face to see if it is part of another label + std::set faceSet; + std::map labels; + for (PetscInt f = fStart; f < fEnd; ++f) { + // check to see if this face is in the boundary label, if not skip for now + if (!boundaryRegion->InRegion(dm, f)) { + continue; + } + + // reset the faceSet + faceSet.clear(); + + // This returns everything associated with the cell in the correct ordering + PetscInt cl, nClosure, *closure = nullptr; + DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &nClosure, &closure) >> utilities::PetscUtilities::checkError; + + // we look at every other point in the closure because we do not need point orientations + for (cl = 0; cl < nClosure * 2; cl += 2) { + if (closure[cl] >= vStart && closure[cl] < vEnd) { // Only use the points corresponding to a vertex + faceSet.insert(closure[cl] - vStart); // we subtract away vStart to reset the node numbering + } + } + + DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &nClosure, &closure) >> utilities::PetscUtilities::checkError; // Restore the points + + // see if the face gets a label + if (auto region = description->GetRegion(faceSet)) { + auto& label = labels[*region]; + PetscInt labelValue = region->GetValue(); + // create and get the label if needed + if (!label) { + region->CreateLabel(dm, label, labelValue); + } + + // Set the value + DMLabelSetValue(label, f, labelValue) >> utilities::PetscUtilities::checkError; + } + } + // complete each label + for (auto& [region, label] : labels) { + DMPlexLabelComplete(dm, label) >> utilities::PetscUtilities::checkError; + } +} +#include "registrar.hpp" +REGISTER(ablate::domain::Domain, ablate::domain::MeshGenerator, "The MeshGenerator will use a mesh description to generate an arbitrary mesh based upon supplied cells/nodes from the description.", + ARG(std::string, "name", "the name of the domain/mesh object"), OPT(std::vector, "fields", "a list of fields/field descriptors"), + ARG(ablate::domain::descriptions::MeshDescription, "description", "the mesh description used to describe the cell/nodes used to create the new mesh"), + OPT(std::vector, "modifiers", "a list of domain modifier"), + OPT(ablate::parameters::Parameters, "options", "PETSc options specific to this dm. Default value allows the dm to access global options.")); diff --git a/src/domain/meshGenerator.hpp b/src/domain/meshGenerator.hpp new file mode 100644 index 000000000..a04d9c317 --- /dev/null +++ b/src/domain/meshGenerator.hpp @@ -0,0 +1,51 @@ +#ifndef ABLATELIBRARY_MESHGENERATOR_HPP +#define ABLATELIBRARY_MESHGENERATOR_HPP + +#include +#include +#include +#include "domain.hpp" +#include "domain/descriptions/meshDescription.hpp" +#include "parameters/parameters.hpp" + +namespace ablate::domain { + +class MeshGenerator : public Domain { + private: + /** + * Private function to read the description and generate the mesh + * @param name + * @return + */ + static DM CreateDM(const std::string& name, const std::shared_ptr& description); + + /** + * Helper function to replace the dm with a new dm + * @param originalDm + * @param replaceDm + */ + static void ReplaceDm(DM& originalDm, DM& replaceDm); + + /** + * Use the mesh description to label the boundaries + * @param description + * @param dm + */ + static void LabelBoundaries(const std::shared_ptr& description, DM& dm); + + public: + /** + * Create a mesh using a simple description of nodes/vertices + * @param name + * @param fieldDescriptors + * @param description + * @param modifiers + * @param options + */ + MeshGenerator(const std::string& name, std::vector> fieldDescriptors, const std::shared_ptr& description, + std::vector> modifiers, const std::shared_ptr& options = {}); + + ~MeshGenerator() override; +}; +} // namespace ablate::domain +#endif // ABLATELIBRARY_MESHGENERATOR_HPP diff --git a/src/domain/modifiers/extrudeLabel.cpp b/src/domain/modifiers/extrudeLabel.cpp index 9bebae89f..def57adf9 100644 --- a/src/domain/modifiers/extrudeLabel.cpp +++ b/src/domain/modifiers/extrudeLabel.cpp @@ -1,5 +1,6 @@ #include "extrudeLabel.hpp" #include +#include #include #include #include "tagLabelInterface.hpp" @@ -84,18 +85,49 @@ void ablate::domain::modifiers::ExtrudeLabel::Modify(DM &dm) { originalRegion->CreateLabel(dmAdapt, originalRegionLabel, originalRegionValue); extrudedRegion->CreateLabel(dmAdapt, extrudedRegionLabel, extrudedRegionValue); - // Determine the current max cell int - PetscInt originalMaxCell; - DMPlexGetHeightStratum(dm, 0, nullptr, &originalMaxCell) >> utilities::PetscUtilities::checkError; - // March over each cell in this rank and determine if it is original or not - PetscInt cStart, cEnd; - DMPlexGetHeightStratum(dmAdapt, 0, &cStart, &cEnd) >> utilities::PetscUtilities::checkError; - for (PetscInt c = cStart; c < cEnd; ++c) { - if (c < originalMaxCell) { - DMLabelSetValue(originalRegionLabel, c, originalRegionValue) >> utilities::PetscUtilities::checkError; - } else { - DMLabelSetValue(extrudedRegionLabel, c, extrudedRegionValue) >> utilities::PetscUtilities::checkError; + PetscInt cAdaptStart, cAdaptEnd; + DMPlexGetHeightStratum(dmAdapt, 0, &cAdaptStart, &cAdaptEnd) >> utilities::PetscUtilities::checkError; + + // cell the depths of the cell layer + PetscInt cellDepth; + DMPlexGetDepth(dm, &cellDepth) >> utilities::PetscUtilities::checkError; + DMLabel depthAdaptLabel, ctAdaptLabel, ctLabel; + DMPlexGetDepthLabel(dmAdapt, &depthAdaptLabel) >> utilities::PetscUtilities::checkError; + DMPlexGetCellTypeLabel(dmAdapt, &ctAdaptLabel) >> utilities::PetscUtilities::checkError; + DMPlexGetCellTypeLabel(dm, &ctLabel) >> utilities::PetscUtilities::checkError; + + // because the new cells can be intertwined with the old cells for mixed use we need to do this cell type by cell type + for (PetscInt cellType = 0; cellType < DM_NUM_POLYTOPES; ++cellType) { + auto ict = (DMPolytopeType)cellType; + + // get the new range for this cell type + PetscInt tAdaptStart, tAdaptEnd; + DMLabelGetStratumBounds(ctAdaptLabel, ict, &tAdaptStart, &tAdaptEnd) >> utilities::PetscUtilities::checkError; + + // only check if there are cell of this type + if (tAdaptStart < 0) { + continue; + } + // determine the depth of this cell type + PetscInt cellTypeDepth; + DMLabelGetValue(depthAdaptLabel, tAdaptStart, &cellTypeDepth) >> utilities::PetscUtilities::checkError; + if (cellTypeDepth != cellDepth) { + continue; + } + + // Get the original range for this cell type + PetscInt tStart, tEnd; + DMLabelGetStratumBounds(ctLabel, ict, &tStart, &tEnd) >> utilities::PetscUtilities::checkError; + PetscInt numberOldCells = tStart >= 0 ? tEnd - tStart : 0; + + // march over each new cell + for (PetscInt c = tAdaptStart; c < tAdaptEnd; ++c) { + if ((c - tAdaptStart) < numberOldCells) { + DMLabelSetValue(originalRegionLabel, c, originalRegionValue) >> utilities::PetscUtilities::checkError; + } else { + DMLabelSetValue(extrudedRegionLabel, c, extrudedRegionValue) >> utilities::PetscUtilities::checkError; + } } } @@ -134,6 +166,7 @@ PetscErrorCode ablate::domain::modifiers::ExtrudeLabel::DMPlexTransformAdaptLabe PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); PetscCall(DMCopyDisc(dm, *rdm)); PetscCall(DMPlexTransformDestroy(&tr)); + ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; PetscFunctionReturn(0); } diff --git a/src/domain/range.cpp b/src/domain/range.cpp index 8a7466098..60c4faa11 100644 --- a/src/domain/range.cpp +++ b/src/domain/range.cpp @@ -2,7 +2,7 @@ #include // For ISIntersect_Caching_Internal, used in ablate::domain::SubDomain::GetRange #include "utilities/petscUtilities.hpp" -void ablate::domain::GetRange(DM dm, const std::shared_ptr region, PetscInt depth, ablate::domain::Range &range) { +void ablate::domain::GetRange(DM dm, const std::shared_ptr ®ion, PetscInt depth, ablate::domain::Range &range) { // Start out getting all the points IS allPointIS; DMGetStratumIS(dm, "dim", depth, &allPointIS) >> utilities::PetscUtilities::checkError; @@ -39,14 +39,14 @@ void ablate::domain::GetRange(DM dm, const std::shared_ptr> utilities::PetscUtilities::checkError; } -void ablate::domain::GetCellRange(DM dm, const std::shared_ptr region, ablate::domain::Range &cellRange) { +void ablate::domain::GetCellRange(DM dm, const std::shared_ptr ®ion, ablate::domain::Range &cellRange) { // Start out getting all the cells PetscInt depth; DMPlexGetDepth(dm, &depth) >> utilities::PetscUtilities::checkError; ablate::domain::GetRange(dm, region, depth, cellRange); } -void ablate::domain::GetFaceRange(DM dm, const std::shared_ptr region, ablate::domain::Range &faceRange) { +void ablate::domain::GetFaceRange(DM dm, const std::shared_ptr ®ion, ablate::domain::Range &faceRange) { // Start out getting all the faces PetscInt depth; DMPlexGetDepth(dm, &depth) >> utilities::PetscUtilities::checkError; diff --git a/src/domain/range.hpp b/src/domain/range.hpp index 357f25dc3..b94fb59e2 100644 --- a/src/domain/range.hpp +++ b/src/domain/range.hpp @@ -10,8 +10,8 @@ namespace ablate::domain { */ struct Range { IS is = nullptr; - PetscInt start; - PetscInt end; + PetscInt start{}; + PetscInt end{}; const PetscInt *points = nullptr; /** @@ -29,7 +29,7 @@ struct Range { * @param depth * @param range */ -void GetRange(DM dm, const std::shared_ptr region, PetscInt depth, Range &range); +void GetRange(DM dm, const std::shared_ptr ®ion, PetscInt depth, Range &range); /** * Get the range of cells defined over the region for this solver. @@ -37,7 +37,7 @@ void GetRange(DM dm, const std::shared_ptr region, PetscInt depth, Range * @param region * @param cellRange */ -void GetCellRange(DM dm, const std::shared_ptr region, Range &cellRange); +void GetCellRange(DM dm, const std::shared_ptr ®ion, Range &cellRange); /** * Get the range of faces/edges defined over the region for this solver. @@ -45,7 +45,7 @@ void GetCellRange(DM dm, const std::shared_ptr region, Range &cellRange) * @param region * @param faceRange */ -void GetFaceRange(DM dm, const std::shared_ptr region, Range &faceRange); +void GetFaceRange(DM dm, const std::shared_ptr ®ion, Range &faceRange); /** * Restores the range diff --git a/src/domain/region.cpp b/src/domain/region.cpp index 04e5a0796..a385d512e 100644 --- a/src/domain/region.cpp +++ b/src/domain/region.cpp @@ -3,13 +3,13 @@ #include "utilities/mpiUtilities.hpp" #include "utilities/petscUtilities.hpp" -ablate::domain::Region::Region(std::string name, int valueIn) : name(name), value(valueIn == 0 ? 1 : valueIn) { +ablate::domain::Region::Region(const std::string& name, int valueIn) : name(name), value(valueIn == 0 ? 1 : valueIn) { // Create a unique string auto hashString = name + ":" + std::to_string(value); id = std::hash()(hashString); } -void ablate::domain::Region::CreateLabel(DM dm, DMLabel& regionLabel, PetscInt& regionValue) { +void ablate::domain::Region::CreateLabel(DM dm, DMLabel& regionLabel, PetscInt& regionValue) const { DMCreateLabel(dm, GetName().c_str()) >> utilities::PetscUtilities::checkError; DMGetLabel(dm, GetName().c_str(), ®ionLabel) >> utilities::PetscUtilities::checkError; regionValue = GetValue(); @@ -37,7 +37,7 @@ void ablate::domain::Region::CheckForLabel(DM dm, MPI_Comm comm) const { DMLabel label = nullptr; DMGetLabel(dm, GetName().c_str(), &label) >> utilities::PetscUtilities::checkError; - PetscMPIInt found = (PetscMPIInt)(label != nullptr); + auto found = (PetscMPIInt)(label != nullptr); PetscMPIInt anyFound; MPI_Allreduce(&found, &anyFound, 1, MPI_INT, MPI_MAX, comm) >> utilities::MpiUtilities::checkError; @@ -51,9 +51,7 @@ bool ablate::domain::Region::InRegion(const std::shared_ptr& region, DM if (!region) { return true; } - PetscInt ptValue; - DMGetLabelValue(dm, region->name.c_str(), point, &ptValue) >> utilities::PetscUtilities::checkError; - return ptValue == region->value; + return region->InRegion(dm, point); } std::ostream& ablate::domain::operator<<(std::ostream& os, const ablate::domain::Region& region) { diff --git a/src/domain/region.hpp b/src/domain/region.hpp index ef05cfb2d..b7f621225 100644 --- a/src/domain/region.hpp +++ b/src/domain/region.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include namespace ablate::domain { @@ -22,7 +23,7 @@ class Region { * @param name label name * @param value label value */ - explicit Region(std::string name = {}, int value = 1); + explicit Region(const std::string& name = {}, int value = 1); [[nodiscard]] inline const std::size_t& GetId() const { return id; } @@ -30,7 +31,7 @@ class Region { [[nodiscard]] inline const PetscInt& GetValue() const { return value; } - [[nodiscard]] inline const std::string ToString() const { return name + ":" + std::to_string(value); }; + [[nodiscard]] inline std::string ToString() const { return name + ":" + std::to_string(value); }; /** * create and returns a label/region value @@ -39,12 +40,31 @@ class Region { * @param regionLabel * @param regionValue */ - void CreateLabel(DM dm, DMLabel& regionLabel, PetscInt& regionValue); + void CreateLabel(DM dm, DMLabel& regionLabel, PetscInt& regionValue) const; static void GetLabel(const std::shared_ptr& region, DM dm, DMLabel& regionLabel, PetscInt& regionValue); + /** + * static call to see if a point is in region/or null + * @param region + * @param dm + * @param point + * @return + */ static bool InRegion(const std::shared_ptr& region, DM dm, PetscInt point); + /** + * Non static call to see if this point is in a given region + * @param dm + * @param point + * @return + */ + inline bool InRegion(DM dm, PetscInt point) const { + PetscInt ptValue; + DMGetLabelValue(dm, name.c_str(), point, &ptValue) >> utilities::PetscUtilities::checkError; + return ptValue == value; + } + /** * throws exception if the label is not in the dm * @param region @@ -58,6 +78,19 @@ class Region { * @param dm */ void CheckForLabel(DM dm, MPI_Comm comm) const; + + /** + * provide simple comparsion == operator + * @return + */ + inline bool operator==(const Region& otherRegion) const { return id == otherRegion.id; } + + /** + * simple comparsion operator allowing regions to be used in a map + * @param otherRegion + * @return + */ + inline bool operator<(const Region& otherRegion) const { return id < otherRegion.id; } }; std::ostream& operator<<(std::ostream& os, const Region& region); diff --git a/src/utilities/vectorUtilities.hpp b/src/utilities/vectorUtilities.hpp index 91f469d44..706fab274 100644 --- a/src/utilities/vectorUtilities.hpp +++ b/src/utilities/vectorUtilities.hpp @@ -65,7 +65,7 @@ class VectorUtilities { } /** - * Fills an array based upon a key vector and map + * Fills an vector based upon a key vector and map * @tparam T * @param list * @return @@ -171,7 +171,56 @@ class VectorUtilities { return transformed; } - private: + /** + * Helper function to convert a vector to array and fill in missing values to zero. It will throw an exception if the size is to large + * @tparam T the array/vector type + * @tparam S the desired size of the array + * @param vector the input vector + * @param defaultValue the fill value + * @return + */ + template + static inline std::array ToArray(const std::vector& vector, const T defaultValue) { + std::array array; + + // check to make sure there is room + if (array.size() < vector.size()) { + throw std::invalid_argument("The supplied vector {" + std::to_string(vector.size()) + "} is to large for the array {" + std::to_string(array.size()) + "}"); + } + + std::size_t i = 0; + for (; i < vector.size(); ++i) { + array[i] = vector[i]; + } + + for (; i < array.size(); ++i) { + array[i] = defaultValue; + } + + return array; + } + /** + * Helper function to convert a vector to array and fill in missing values to zero. It will throw an exception if the size of the vector and array are not the same + * @tparam T the array/vector type + * @tparam S the desired size of the array + * @param vector the input vector + * @param defaultValue the fill value + * @return + */ + template + static inline std::array ToArray(const std::vector& vector) { + std::array array; + + // check to make sure there is room + if (array.size() != vector.size()) { + throw std::invalid_argument("The supplied vector {" + std::to_string(vector.size()) + "} is not sized correctly for the array {" + std::to_string(array.size()) + "}"); + } + + std::copy(vector.begin(), vector.end(), array.begin()); + + return array; + } + VectorUtilities() = delete; }; diff --git a/tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.txt b/tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.txt new file mode 100644 index 000000000..7be555113 --- /dev/null +++ b/tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.txt @@ -0,0 +1,13 @@ +Timestep: 0000 time = 0 dt = 1e-10 +Timestep: 0010 time = (.*) dt = (.*) ~ ~ +Timestep: 0020 time = (.*) dt = (.*) ~ ~ +Timestep: 0030 time = (.*) dt = (.*) ~ ~ +Timestep: 0040 time = (.*) dt = (.*) ~ ~ +Timestep: 0050 time = (.*) dt = (.*) ~ ~ +ResultFiles: +domain.hdf5 +domain.xmf +flowField.hdf5 +flowField.xmf +pgsLog +restart.rst diff --git a/tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.yaml b/tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.yaml new file mode 100644 index 000000000..5add5692b --- /dev/null +++ b/tests/integrationTests/inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.yaml @@ -0,0 +1,190 @@ +# Simple flow through an axisymmetric pipe +--- +test: + # a unique test name for this integration tests + name: axisymmetricPipeFlow + # create a default assert that compares the log file + assert: "inputs/compressibleFlow/axisymmetricPipeFlow/axisymmetricPipeFlow.txt" + +# metadata for the simulation +environment: + title: _axisymmetricPipeFlow + tagDirectory: false +# global arguments that can be used by petsc +arguments: [ ] + +# set up the time stepper responsible for marching in time +timestepper: + # write the output to show the mesh + io: + interval: 0 + # time stepper specific input arguments + arguments: + ts_type: rk + ts_max_time: 100000 + ts_max_steps: 50 + ts_dt: 1.0E-10 + ts_adapt_safety: 0.9 + ts_adapt_type: physicsConstrained + # create a simple box mesh for simulation + domain: !ablate::domain::MeshGenerator + name: exampleAxisymmetricMesh + # specify the axisymmetric mesh description + description: !ablate::domain::descriptions::Axisymmetric + start: [ 0.0, 0.0, 0.0 ] + length: 0.5 + radius: ".05 + (z-.25)*(z-.25)" + numberWedges: 15 + numberSlices: 29 + numberShells: 10 + # output the result label for interior cells + options: + # output the flow region to visually check the result + dm_label_view: flowRegion + # force the new mesh to check everything + dm_plex_check_all: true + + # specify any modifications to be performed to the mesh/domain + modifiers: + - # use the newly labels to extrude the boundary. Do not extrude the cell + !ablate::domain::modifiers::ExtrudeLabel + regions: + - name: boundary + # mark all the resulting boundary faces with boundaryFaces label + boundaryRegion: + name: boundaryFaces + # tag the original mesh as the flow region + originalRegion: + name: flowRegion + # tag the new boundary cells for easy boundary condition specifications + extrudedRegion: + name: boundaryCells + + # if using mpi, this modifier distributes cells + - !ablate::domain::modifiers::DistributeWithGhostCells + ghostCellDepth: 2 + + # setup some dummy fields + fields: + # all fields must be defined before solvers. The ablate::finiteVolume::CompressibleFlowFields is a helper + # class that creates the required fields for the compressible flow solver (rho, rhoE, rhoU, ...) + - !ablate::finiteVolume::CompressibleFlowFields + eos: !ablate::eos::PerfectGas &eos + parameters: + gamma: 1.4 + Rgas: 287.0 + # species are added to the flow through the eos. This allows testing of the species transport equations + species: [ N2, H2O, O2 ] + # by adding a pressure field the code will compute and output pressure + - name: pressure + location: AUX + type: FVM + + # initialize the dummy field + initialization: + # ablate::finiteVolume::CompressibleFlowFields is a helper + # class that creates the required fields for the compressible flow solver (rho, rhoE, rhoU, ...) + - !ablate::finiteVolume::fieldFunctions::Euler + state: + &flowFieldState + eos: *eos + pressure: 101325.0 + temperature: 300 + velocity: "0.0, 0.0, 0.0" + # individual mass fractions must be passed to the flow field state to compute density, energy, etc. + other: !ablate::finiteVolume::fieldFunctions::MassFractions + eos: *eos + values: + - fieldName: N2 + field: "z > 0.25 ? .2 : 1.0" + - fieldName: H2O + field: "z> 0.25 ? .3 :0" + - fieldName: O2 + field: " z > 0.25 ? .5 : 0" + + # the same state can be used to internalize the DensityMassFractions field from density and mass fractions + - !ablate::finiteVolume::fieldFunctions::DensityMassFractions + state: *flowFieldState + +# this is a test input file with no solvers +solvers: + # The compressible flow solver will solve the compressible flow equations over the interiorCells + - !ablate::finiteVolume::CompressibleFlowSolver + id: flowField + # only apply this solver to the flowRegion, area without faces + region: + name: flowRegion + additionalProcesses: + - !ablate::finiteVolume::processes::PressureGradientScaling + &pgs + eos: *eos + alphaInit: 100.0 + maxAlphaAllowed: 100.0 + domainLength: 0.165354 + log: !ablate::monitors::logs::CsvLog + name: pgsLog + + # a flux calculator must be specified to so solver for advection + fluxCalculator: !ablate::finiteVolume::fluxCalculator::AusmpUp + pgs: *pgs + + # the default transport object assumes constant values for k, mu, diff + transport: + k: .2 + mu: .1 + diff: 1E-4 + + # cfl is used to compute the physics time step + parameters: + cfl: 0.5 + + # share the existing eos with the compressible flow solver + eos: *eos + + monitors: + # output the timestep and dt at each time step + - !ablate::monitors::TimeStepMonitor + interval: 10 + + # use a boundary solver to update the cells in the gMsh inlet region to represent an inlet + - !ablate::boundarySolver::BoundarySolver + id: inlet + region: + name: lowerCap + fieldBoundary: + name: boundaryFaces + mergeFaces: false + processes: + - !ablate::boundarySolver::lodi::Inlet + eos: *eos + pgs: *pgs + velocity: "0.0, 0.0, min(10, 10*t)" # for stability, increase the velocity slowly + + # use a boundary solver to update the cells in the gMsh outlet region to represent an open pipe + - !ablate::boundarySolver::BoundarySolver + id: openBoundary + region: + name: upperCap + fieldBoundary: + name: boundaryFaces + mergeFaces: true + processes: + - !ablate::boundarySolver::lodi::OpenBoundary + eos: *eos + reflectFactor: 0.0 + referencePressure: 101325.0 + maxAcousticsLength: 1 + pgs: *pgs + + # use a boundary solver to update the cells in the wall region to represent standard wall + - !ablate::boundarySolver::BoundarySolver + id: wall + region: + name: outerShell + fieldBoundary: + name: boundaryFaces + mergeFaces: true + processes: + - !ablate::boundarySolver::lodi::IsothermalWall + eos: *eos + pgs: *pgs diff --git a/tests/integrationTests/inputs/machinery/boundaryMonitorTest2D/bottomBoundary_monitor.xmf b/tests/integrationTests/inputs/domain/boundaryMonitorTest2D/bottomBoundary_monitor.xmf similarity index 100% rename from tests/integrationTests/inputs/machinery/boundaryMonitorTest2D/bottomBoundary_monitor.xmf rename to tests/integrationTests/inputs/domain/boundaryMonitorTest2D/bottomBoundary_monitor.xmf diff --git a/tests/integrationTests/inputs/machinery/boundaryMonitorTest2D/boundaryMonitorTest2D.txt b/tests/integrationTests/inputs/domain/boundaryMonitorTest2D/boundaryMonitorTest2D.txt similarity index 100% rename from tests/integrationTests/inputs/machinery/boundaryMonitorTest2D/boundaryMonitorTest2D.txt rename to tests/integrationTests/inputs/domain/boundaryMonitorTest2D/boundaryMonitorTest2D.txt diff --git a/tests/integrationTests/inputs/machinery/boundaryMonitorTest2D/boundaryMonitorTest2D.yaml b/tests/integrationTests/inputs/domain/boundaryMonitorTest2D/boundaryMonitorTest2D.yaml similarity index 96% rename from tests/integrationTests/inputs/machinery/boundaryMonitorTest2D/boundaryMonitorTest2D.yaml rename to tests/integrationTests/inputs/domain/boundaryMonitorTest2D/boundaryMonitorTest2D.yaml index a151e8b1c..33c2aab71 100644 --- a/tests/integrationTests/inputs/machinery/boundaryMonitorTest2D/boundaryMonitorTest2D.yaml +++ b/tests/integrationTests/inputs/domain/boundaryMonitorTest2D/boundaryMonitorTest2D.yaml @@ -6,9 +6,9 @@ test: ranks: 2 asserts: # create a default assert that compares the log file - - "inputs/machinery/boundaryMonitorTest2D/boundaryMonitorTest2D.txt" + - "inputs/domain/boundaryMonitorTest2D/boundaryMonitorTest2D.txt" - !testingResources::asserts::TextFileAssert - expected: "inputs/machinery/boundaryMonitorTest2D/bottomBoundary_monitor.xmf" + expected: "inputs/domain/boundaryMonitorTest2D/bottomBoundary_monitor.xmf" actual: "bottomBoundary_monitor.xmf" # metadata for the simulation diff --git a/tests/integrationTests/inputs/machinery/boundaryMonitorTest3D/bottomBoundary_monitor.xmf b/tests/integrationTests/inputs/domain/boundaryMonitorTest3D/bottomBoundary_monitor.xmf similarity index 100% rename from tests/integrationTests/inputs/machinery/boundaryMonitorTest3D/bottomBoundary_monitor.xmf rename to tests/integrationTests/inputs/domain/boundaryMonitorTest3D/bottomBoundary_monitor.xmf diff --git a/tests/integrationTests/inputs/machinery/boundaryMonitorTest3D/boundaryMonitorTest3D.txt b/tests/integrationTests/inputs/domain/boundaryMonitorTest3D/boundaryMonitorTest3D.txt similarity index 100% rename from tests/integrationTests/inputs/machinery/boundaryMonitorTest3D/boundaryMonitorTest3D.txt rename to tests/integrationTests/inputs/domain/boundaryMonitorTest3D/boundaryMonitorTest3D.txt diff --git a/tests/integrationTests/inputs/machinery/boundaryMonitorTest3D/boundaryMonitorTest3D.yaml b/tests/integrationTests/inputs/domain/boundaryMonitorTest3D/boundaryMonitorTest3D.yaml similarity index 96% rename from tests/integrationTests/inputs/machinery/boundaryMonitorTest3D/boundaryMonitorTest3D.yaml rename to tests/integrationTests/inputs/domain/boundaryMonitorTest3D/boundaryMonitorTest3D.yaml index cc4a8310f..2755c63b1 100644 --- a/tests/integrationTests/inputs/machinery/boundaryMonitorTest3D/boundaryMonitorTest3D.yaml +++ b/tests/integrationTests/inputs/domain/boundaryMonitorTest3D/boundaryMonitorTest3D.yaml @@ -5,9 +5,9 @@ test: name: boundaryMonitorTest3D asserts: # create a default assert that compares the log file - - "inputs/machinery/boundaryMonitorTest3D/boundaryMonitorTest3D.txt" + - "inputs/domain/boundaryMonitorTest3D/boundaryMonitorTest3D.txt" - !testingResources::asserts::TextFileAssert - expected: "inputs/machinery/boundaryMonitorTest3D/bottomBoundary_monitor.xmf" + expected: "inputs/domain/boundaryMonitorTest3D/bottomBoundary_monitor.xmf" actual: "bottomBoundary_monitor.xmf" # metadata for the simulation diff --git a/tests/integrationTests/inputs/machinery/dmViewFromOptions.txt b/tests/integrationTests/inputs/domain/dmViewFromOptions.txt similarity index 100% rename from tests/integrationTests/inputs/machinery/dmViewFromOptions.txt rename to tests/integrationTests/inputs/domain/dmViewFromOptions.txt diff --git a/tests/integrationTests/inputs/machinery/dmViewFromOptions.yaml b/tests/integrationTests/inputs/domain/dmViewFromOptions.yaml similarity index 97% rename from tests/integrationTests/inputs/machinery/dmViewFromOptions.yaml rename to tests/integrationTests/inputs/domain/dmViewFromOptions.yaml index d1453d150..8ad062b3c 100644 --- a/tests/integrationTests/inputs/machinery/dmViewFromOptions.yaml +++ b/tests/integrationTests/inputs/domain/dmViewFromOptions.yaml @@ -5,7 +5,7 @@ test: # a unique test name for this integration tests name: dmViewFromOptions # create a default assert that compares the log file - assert: "inputs/machinery/dmViewFromOptions.txt" + assert: "inputs/domain/dmViewFromOptions.txt" environment: title: incompessibleFlow diff --git a/tests/integrationTests/inputs/machinery/extrudeBoundaryTest.txt b/tests/integrationTests/inputs/domain/extrudeBoundaryTest.txt similarity index 100% rename from tests/integrationTests/inputs/machinery/extrudeBoundaryTest.txt rename to tests/integrationTests/inputs/domain/extrudeBoundaryTest.txt diff --git a/tests/integrationTests/inputs/machinery/extrudeBoundaryTest.yaml b/tests/integrationTests/inputs/domain/extrudeBoundaryTest.yaml similarity index 98% rename from tests/integrationTests/inputs/machinery/extrudeBoundaryTest.yaml rename to tests/integrationTests/inputs/domain/extrudeBoundaryTest.yaml index 93662bf4b..22d6aa8a5 100644 --- a/tests/integrationTests/inputs/machinery/extrudeBoundaryTest.yaml +++ b/tests/integrationTests/inputs/domain/extrudeBoundaryTest.yaml @@ -4,7 +4,7 @@ test: # a unique test name for this integration tests name: extrudeBoundaryTest # create a default assert that compares the log file - assert: "inputs/machinery/extrudeBoundaryTest.txt" + assert: "inputs/domain/extrudeBoundaryTest.txt" # metadata for the simulation environment: diff --git a/tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.txt b/tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.txt new file mode 100644 index 000000000..7f7872373 --- /dev/null +++ b/tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.txt @@ -0,0 +1,18 @@ +DM Object: exampleAxisymmetricMesh 1 MPI process + type: plex +exampleAxisymmetricMesh in 3 dimensions: + Number of 0-cells per rank: 8241 + Number of 1-cells per rank: 24440 + Number of 2-cells per rank: 24200 (23380) + Number of 3-cells per rank: 8000 (800) +Labels: + celltype: 6 strata with value/size (0 (8241), 7 (7200), 8 (800), 3 (820), 1 (24440), 4 (23380)) + depth: 4 strata with value/size (0 (8241), 1 (24440), 2 (24200), 3 (8000)) + boundary: 1 strata with value/size (1 (4762)) + upperCap: 1 strata with value/size (1 (801)) + lowerCap: 1 strata with value/size (1 (801)) + outerShell: 1 strata with value/size (1 (3240)) +ResultFiles: +domain.hdf5 +domain.xmf +restart.rst diff --git a/tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.yaml b/tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.yaml new file mode 100644 index 000000000..26df2daac --- /dev/null +++ b/tests/integrationTests/inputs/domain/meshGeneratorAxisymmetric.yaml @@ -0,0 +1,53 @@ +# Simple test that checks the functionality of the Axisymmetric mesh generator +--- +test: + # a unique test name for this integration tests + name: meshGeneratorAxisymmetric + # create a default assert that compares the log file + assert: "inputs/domain/meshGeneratorAxisymmetric.txt" + +# metadata for the simulation +environment: + title: _meshGeneratorAxisymmetric + tagDirectory: false +# global arguments that can be used by petsc +arguments: + # force the new mesh to check everything + dm_plex_check_all: "" + +# set up the time stepper responsible for marching in time +timestepper: + # write the output to show the mesh + io: + interval: 0 + # for this example there are no time stepper arguments (empty simulation) + arguments: { } + # create a simple box mesh for simulation + domain: !ablate::domain::MeshGenerator + name: exampleAxisymmetricMesh + # specify the axisymmetric mesh description + description: !ablate::domain::descriptions::Axisymmetric + start: [0.0, 0.0, 0.0] + length: 2.0 + radius: ".1*z*z+ 0.5" + numberWedges: 20 + numberSlices: 40 + numberShells: 10 + + # setup some dummy fields + fields: + - name: exampleFVField + components: [ "xx", "yy", "zz" ] + type: FVM + modifiers: + + # the DmViewFromOptions "modifier" does not modify the dm but outputs. See [PetscOptionsGetViewer](https://petsc.org/release/docs/manualpages/Viewer/PetscOptionsGetViewer.html) for more details + - !ablate::monitors::DmViewFromOptions + options: ascii + # initialize the dummy field + initialization: + - fieldName: "exampleFVField" + field: "x, y, z" + +# this is a test input file with no solvers +solvers: [ ] \ No newline at end of file diff --git a/tests/integrationTests/inputs/machinery/meshMappingTest.txt b/tests/integrationTests/inputs/domain/meshMappingTest.txt similarity index 100% rename from tests/integrationTests/inputs/machinery/meshMappingTest.txt rename to tests/integrationTests/inputs/domain/meshMappingTest.txt diff --git a/tests/integrationTests/inputs/machinery/meshMappingTest.yaml b/tests/integrationTests/inputs/domain/meshMappingTest.yaml similarity index 97% rename from tests/integrationTests/inputs/machinery/meshMappingTest.yaml rename to tests/integrationTests/inputs/domain/meshMappingTest.yaml index 747940e3f..b3ccdf949 100644 --- a/tests/integrationTests/inputs/machinery/meshMappingTest.yaml +++ b/tests/integrationTests/inputs/domain/meshMappingTest.yaml @@ -4,7 +4,7 @@ test: # a unique test name for this integration tests name: meshMappingTest # create a default assert that compares the log file - assert: "inputs/machinery/meshMappingTest.txt" + assert: "inputs/domain/meshMappingTest.txt" # metadata for the simulation environment: diff --git a/tests/integrationTests/inputs/machinery/meshMappingTestCoordinateSpace.txt b/tests/integrationTests/inputs/domain/meshMappingTestCoordinateSpace.txt similarity index 100% rename from tests/integrationTests/inputs/machinery/meshMappingTestCoordinateSpace.txt rename to tests/integrationTests/inputs/domain/meshMappingTestCoordinateSpace.txt diff --git a/tests/integrationTests/inputs/machinery/meshMappingTestCoordinateSpace.yaml b/tests/integrationTests/inputs/domain/meshMappingTestCoordinateSpace.yaml similarity index 96% rename from tests/integrationTests/inputs/machinery/meshMappingTestCoordinateSpace.yaml rename to tests/integrationTests/inputs/domain/meshMappingTestCoordinateSpace.yaml index 3c226892c..acdfa237c 100644 --- a/tests/integrationTests/inputs/machinery/meshMappingTestCoordinateSpace.yaml +++ b/tests/integrationTests/inputs/domain/meshMappingTestCoordinateSpace.yaml @@ -4,7 +4,7 @@ test: # a unique test name for this integration tests name: meshMappingTestCoordinateSpace # create a default assert that compares the log file - assert: "inputs/machinery/meshMappingTestCoordinateSpace.txt" + assert: "inputs/domain/meshMappingTestCoordinateSpace.txt" # metadata for the simulation environment: diff --git a/tests/integrationTests/inputs/machinery/mixedCellTypeTest/domain.xmf b/tests/integrationTests/inputs/domain/mixedCellTypeTest/domain.xmf similarity index 100% rename from tests/integrationTests/inputs/machinery/mixedCellTypeTest/domain.xmf rename to tests/integrationTests/inputs/domain/mixedCellTypeTest/domain.xmf diff --git a/tests/integrationTests/inputs/machinery/mixedCellTypeTest/mixedCellTypeTest2D.yaml b/tests/integrationTests/inputs/domain/mixedCellTypeTest/mixedCellTypeTest2D.yaml similarity index 97% rename from tests/integrationTests/inputs/machinery/mixedCellTypeTest/mixedCellTypeTest2D.yaml rename to tests/integrationTests/inputs/domain/mixedCellTypeTest/mixedCellTypeTest2D.yaml index 6029937d9..9b27333af 100644 --- a/tests/integrationTests/inputs/machinery/mixedCellTypeTest/mixedCellTypeTest2D.yaml +++ b/tests/integrationTests/inputs/domain/mixedCellTypeTest/mixedCellTypeTest2D.yaml @@ -8,7 +8,7 @@ test: ranks: 1 # compare the generated xmf file with the expected assert: !testingResources::asserts::TextFileAssert - expected: "inputs/machinery/mixedCellTypeTest/domain.xmf" + expected: "inputs/domain/mixedCellTypeTest/domain.xmf" actual: "domain.xmf" # metadata for the simulation diff --git a/tests/integrationTests/inputs/machinery/mixedCellTypeTest/mixedCells2D.msh b/tests/integrationTests/inputs/domain/mixedCellTypeTest/mixedCells2D.msh similarity index 100% rename from tests/integrationTests/inputs/machinery/mixedCellTypeTest/mixedCells2D.msh rename to tests/integrationTests/inputs/domain/mixedCellTypeTest/mixedCells2D.msh diff --git a/tests/integrationTests/inputs/machinery/subDomainFVM.yaml b/tests/integrationTests/inputs/domain/subDomainFVM.yaml similarity index 97% rename from tests/integrationTests/inputs/machinery/subDomainFVM.yaml rename to tests/integrationTests/inputs/domain/subDomainFVM.yaml index f5c873c5e..b7730b8e6 100644 --- a/tests/integrationTests/inputs/machinery/subDomainFVM.yaml +++ b/tests/integrationTests/inputs/domain/subDomainFVM.yaml @@ -5,9 +5,9 @@ test: name: subDomainFVM asserts: # create a default assert that compares the log file - - "inputs/machinery/subDomainFVM/subDomainFVM.txt" + - "inputs/domain/subDomainFVM/subDomainFVM.txt" - !testingResources::asserts::TextFileAssert - expected: "inputs/machinery/subDomainFVM/fluidField.xmf" + expected: "inputs/domain/subDomainFVM/fluidField.xmf" actual: "fluidField.xmf" # metadata for the simulation diff --git a/tests/integrationTests/inputs/machinery/subDomainFVM/fluidField.xmf b/tests/integrationTests/inputs/domain/subDomainFVM/fluidField.xmf similarity index 100% rename from tests/integrationTests/inputs/machinery/subDomainFVM/fluidField.xmf rename to tests/integrationTests/inputs/domain/subDomainFVM/fluidField.xmf diff --git a/tests/integrationTests/inputs/machinery/subDomainFVM/subDomainFVM.txt b/tests/integrationTests/inputs/domain/subDomainFVM/subDomainFVM.txt similarity index 100% rename from tests/integrationTests/inputs/machinery/subDomainFVM/subDomainFVM.txt rename to tests/integrationTests/inputs/domain/subDomainFVM/subDomainFVM.txt diff --git a/tests/unitTests/domain/CMakeLists.txt b/tests/unitTests/domain/CMakeLists.txt index 18806083e..02c35e84a 100644 --- a/tests/unitTests/domain/CMakeLists.txt +++ b/tests/unitTests/domain/CMakeLists.txt @@ -12,3 +12,4 @@ target_sources(ablateUnitTestLibrary add_subdirectory(modifiers) add_subdirectory(RBF) +add_subdirectory(descriptions) diff --git a/tests/unitTests/domain/descriptions/CMakeLists.txt b/tests/unitTests/domain/descriptions/CMakeLists.txt new file mode 100644 index 000000000..1fc51950e --- /dev/null +++ b/tests/unitTests/domain/descriptions/CMakeLists.txt @@ -0,0 +1,4 @@ +target_sources(ablateUnitTestLibrary + PRIVATE + axisymmetricTests.cpp + ) diff --git a/tests/unitTests/domain/descriptions/axisymmetricTests.cpp b/tests/unitTests/domain/descriptions/axisymmetricTests.cpp new file mode 100644 index 000000000..dc4387ad4 --- /dev/null +++ b/tests/unitTests/domain/descriptions/axisymmetricTests.cpp @@ -0,0 +1,239 @@ +#include +#include "domain/descriptions/axisymmetric.hpp" +#include "gtest/gtest.h" +#include "mathFunctions/functionFactory.hpp" + +struct AxisymmetricMeshDescriptionParameters { + // Hold an mesh generator + std::shared_ptr description; + + // Expected values + PetscInt expectedMeshDimension; + PetscInt expectedNumberCells; + PetscInt expectedNumberVertices; + std::map expectedCellTypes; + std::map> expectedTopology; + std::map> expectedVertices; + std::map, std::vector>> expectedBoundaries; +}; + +class AxisymmetricMeshDescriptionFixture : public ::testing::TestWithParam {}; + +TEST_P(AxisymmetricMeshDescriptionFixture, ShouldComputeCorrectNumberOfCellsAndVertices) { + // arrange + auto description = GetParam().description; + // act + // assert + ASSERT_EQ(GetParam().expectedNumberCells, description->GetNumberCells()); + ASSERT_EQ(GetParam().expectedNumberVertices, description->GetNumberVertices()); + ASSERT_EQ(GetParam().expectedMeshDimension, description->GetMeshDimension()); +} + +TEST_P(AxisymmetricMeshDescriptionFixture, ShouldProduceCorrectCellTypes) { + // arrange + auto description = GetParam().description; + // act + // assert + for (const auto& test : GetParam().expectedCellTypes) { + ASSERT_EQ(test.second, description->GetCellType(test.first)) << "For Cell " << test.first << ". Check https://petsc.org/release/manualpages/DM/DMPolytopeType/ for types." << std::endl; + } +} + +TEST_P(AxisymmetricMeshDescriptionFixture, ShouldBuildCorrectTopology) { + // arrange + auto description = GetParam().description; + // act + // assert + for (const auto& test : GetParam().expectedTopology) { + // Build the testCellNodes + std::vector testNodes(test.second.size()); + description->BuildTopology(test.first, testNodes.data()); + + ASSERT_EQ(test.second, testNodes) << "For Cell " << test.first << std::endl; + } +} + +TEST_P(AxisymmetricMeshDescriptionFixture, ShouldBuildCorrectVertices) { + // arrange + auto description = GetParam().description; + // act + // assert + for (const auto& test : GetParam().expectedVertices) { + // Build the testCellNodes + std::vector testVertex(test.second.size()); + description->SetCoordinate(test.first, testVertex.data()); + + for (std::size_t i = 0; i < testVertex.size(); ++i) { + ASSERT_NEAR(test.second[i], testVertex[i], 1E-8) << "For Node " << test.first << " index " << i << std::endl; + } + } +} + +TEST_P(AxisymmetricMeshDescriptionFixture, ShouldIdentifyCorrectBoundaries) { + // arrange + auto description = GetParam().description; + // act + // assert + for (const auto& test : GetParam().expectedBoundaries) { + // march over each node set + for (const auto& nodeSet : test.second) { + auto region = description->GetRegion(nodeSet); + + // build the error nodes if needed + std::stringstream nodes; + for (const auto& n : nodeSet) { + nodes << n << " "; + } + if (test.first) { + if (!region) { + FAIL() << "Should be region for nodes: " << nodes.str() << std::endl; + } + ASSERT_EQ(*test.first, *region) << "Should be equal for nodes: " << nodes.str() << std::endl; + } else { + ASSERT_EQ(region, nullptr) << "Should be null " << nodes.str() << std::endl; + } + } + } +} + +INSTANTIATE_TEST_SUITE_P(AxisymmetricTests, AxisymmetricMeshDescriptionFixture, + testing::Values( + // 1 slice, 1 shell + AxisymmetricMeshDescriptionParameters{ + .description = std::make_shared(std::vector{0.0, 0.0, 0.0}, 0.5, ablate::mathFunctions::Create(1.0), 8, 1, 1), + .expectedMeshDimension = 3, + .expectedNumberCells = 8, + .expectedNumberVertices = 18, + .expectedCellTypes = {{7, DM_POLYTOPE_TRI_PRISM}, {3, DM_POLYTOPE_TRI_PRISM}, {0, DM_POLYTOPE_TRI_PRISM}}, + .expectedTopology = {{7, {0, 3, 2, 1, 10, 11}}, {0, {0, 2, 9, 1, 17, 10}}}, + .expectedVertices = {{0, {0.0, 0.0, 0.0}}, {2, {1.0, 0.0, 0.0}}, {4, {0.0, 1.0, 0.0}}, {1, {0.0, 0.0, 0.5}}, {10, {1.0, 0.0, 0.5}}, {12, {0.0, 1.0, 0.5}}}, + .expectedBoundaries = {{ablate::domain::Region::ENTIREDOMAIN, {{0, 3, 1, 11}}}, + {std::make_shared("outerShell"), {{2, 9, 10, 17}, {5, 4, 12, 13}}}, + {std::make_shared("lowerCap"), {{0, 2, 3}, {0, 5, 6}}}, + {std::make_shared("upperCap"), {{1, 16, 17}, {1, 17, 10}}}}}, + // 2 slices, 1 shell + AxisymmetricMeshDescriptionParameters{ + .description = std::make_shared(std::vector{0.0, 0.0, 0.0}, 0.5, ablate::mathFunctions::Create(1.0), 8, 2, 1), + .expectedMeshDimension = 3, + .expectedNumberCells = 16, + .expectedNumberVertices = 27, + .expectedCellTypes = {{15, DM_POLYTOPE_TRI_PRISM}, {8, DM_POLYTOPE_TRI_PRISM}, {7, DM_POLYTOPE_TRI_PRISM}, {0, DM_POLYTOPE_TRI_PRISM}}, + .expectedTopology = {{15, {0, 4, 3, 1, 11, 12}}, {8, {0, 3, 10, 1, 18, 11}}, {7, {1, 12, 11, 2, 19, 20}}, {0, {1, 11, 18, 2, 26, 19}}}, + .expectedVertices = {{0, {0.0, 0.0, 0.0}}, + {3, {1.0, 0.0, 0.0}}, + {5, {0.0, 1.0, 0.0}}, + {1, {0.0, 0.0, 0.25}}, + {11, {1.0, 0.0, 0.25}}, + {13, {0.0, 1.0, 0.25}}, + {2, {0.0, 0.0, 0.5}}, + {19, {1.0, 0.0, 0.5}}, + {21, {0.0, 1.0, 0.5}}}, + .expectedBoundaries = {{ablate::domain::Region::ENTIREDOMAIN, {{1, 12, 13}}}, + {std::make_shared("outerShell"), {{3, 4, 11, 12}, {19, 25, 11, 18}}}, + {std::make_shared("lowerCap"), {{0, 3, 4}, {0, 10, 3}}}, + {std::make_shared("upperCap"), {{2, 19, 20}, {2, 19, 26}}}}}, + // 1 slice, 2 shells + AxisymmetricMeshDescriptionParameters{ + .description = std::make_shared(std::vector{0.0, 0.0, 0.0}, 0.5, ablate::mathFunctions::Create(1.0), 8, 1, 2), + .expectedMeshDimension = 3, + .expectedNumberCells = 16, + .expectedNumberVertices = 34, + .expectedCellTypes = {{15, DM_POLYTOPE_TRI_PRISM}, {8, DM_POLYTOPE_TRI_PRISM}, {7, DM_POLYTOPE_HEXAHEDRON}, {0, DM_POLYTOPE_HEXAHEDRON}}, + .expectedTopology = {{15, {0, 3, 2, 1, 10, 11}}, {8, {0, 2, 9, 1, 17, 10}}, {7, {10, 2, 18, 26, 11, 27, 19, 3}}, {0, {17, 9, 25, 33, 10, 26, 18, 2}}}, + .expectedVertices = {{0, {0.0, 0.0, 0.0}}, + {2, {0.5, 0.0, 0.0}}, + {4, {0.0, 0.5, 0.0}}, + {18, {1.0, 0.0, 0.0}}, + {20, {0.0, 1.0, 0.0}}, + {1, {0.0, 0.0, 0.5}}, + {10, {0.5, 0.0, 0.5}}, + {12, {0.0, 0.5, 0.5}}, + {26, {1.0, 0.0, 0.5}}, + {28, {0.0, 1.0, 0.5}}}, + .expectedBoundaries = {{ablate::domain::Region::ENTIREDOMAIN, {{0, 3, 11}}}, + {std::make_shared("outerShell"), {{18, 19, 26, 27}, {26, 33, 18, 25}, {23, 24, 31, 32}}}, + {std::make_shared("lowerCap"), {{0, 2, 3}, {0, 9, 2}, {2, 19, 19, 3}, {2, 18, 9, 25}}}, + {std::make_shared("upperCap"), {{1, 10, 11}, {1, 17, 10}, {10, 26, 11, 27}, {10, 26, 17, 33}}}}}, + // 2 slice, 2 shells + AxisymmetricMeshDescriptionParameters{ + .description = std::make_shared(std::vector{0.0, 0.0, 0.0}, 1.0, ablate::mathFunctions::Create(1.0), 8, 2, 2), + .expectedMeshDimension = 3, + .expectedNumberCells = 32, + .expectedNumberVertices = 51, + .expectedCellTypes = {{31, DM_POLYTOPE_TRI_PRISM}, + {24, DM_POLYTOPE_TRI_PRISM}, + {23, DM_POLYTOPE_TRI_PRISM}, + {16, DM_POLYTOPE_TRI_PRISM}, + {15, DM_POLYTOPE_HEXAHEDRON}, + {8, DM_POLYTOPE_HEXAHEDRON}, + {7, DM_POLYTOPE_HEXAHEDRON}, + {0, DM_POLYTOPE_HEXAHEDRON}}, + .expectedTopology = {{31, {0, 4, 3, 1, 11, 12}}, + {24, {0, 3, 10, 1, 18, 11}}, + {15, {11, 3, 27, 35, 12, 36, 28, 4}}, + {8, {18, 10, 34, 42, 11, 35, 27, 3}}, + {23, {1, 12, 11, 2, 19, 20}}, + {16, {1, 11, 18, 2, 26, 19}}, + {7, {19, 11, 35, 43, 20, 44, 36, 12}}, + {0, {26, 18, 42, 50, 19, 43, 35, 11}}}, + .expectedVertices = {{0, {0.0, 0.0, 0.0}}, + {3, {0.5, 0.0, 0.0}}, + {5, {0.0, 0.5, 0.0}}, + {27, {1.0, 0.0, 0.0}}, + {29, {0.0, 1.0, 0.0}}, + {1, {0.0, 0.0, 0.5}}, + {11, {0.5, 0.0, 0.5}}, + {13, {0.0, 0.5, 0.5}}, + {35, {1.0, 0.0, 0.5}}, + {37, {0.0, 1.0, 0.5}}, + {2, {0.0, 0.0, 1.0}}, + {19, {0.5, 0.0, 1.0}}, + {21, {0.0, 0.5, 1.0}}, + {43, {1.0, 0.0, 1.0}}, + {45, {0.0, 1.0, 1.0}}}, + .expectedBoundaries = {{ablate::domain::Region::ENTIREDOMAIN, {{1, 11, 12}}}, + {std::make_shared("outerShell"), {{28, 27, 36, 35}, {44, 43, 36, 35}, {48, 49, 30, 41}}}, + {std::make_shared("lowerCap"), {{0, 3, 4}, {0, 10, 3}, {3, 4, 28, 27}, {10, 3, 27, 34}}}, + {std::make_shared("upperCap"), {{2, 19, 20}, {2, 19, 26}, {26, 50, 43, 19}}}}}, + // 2 slice, 2 shells, variable radius + AxisymmetricMeshDescriptionParameters{ + .description = std::make_shared(std::vector{0.0, 0.0, 0.0}, 2.0, ablate::mathFunctions::Create(".1*z+ 0.5"), 8, + 2, 2), + .expectedMeshDimension = 3, + .expectedNumberCells = 32, + .expectedNumberVertices = 51, + .expectedCellTypes = {{31, DM_POLYTOPE_TRI_PRISM}, + {24, DM_POLYTOPE_TRI_PRISM}, + {23, DM_POLYTOPE_TRI_PRISM}, + {16, DM_POLYTOPE_TRI_PRISM}, + {15, DM_POLYTOPE_HEXAHEDRON}, + {8, DM_POLYTOPE_HEXAHEDRON}, + {7, DM_POLYTOPE_HEXAHEDRON}, + {0, DM_POLYTOPE_HEXAHEDRON}}, + .expectedTopology = {{31, {0, 4, 3, 1, 11, 12}}, + {24, {0, 3, 10, 1, 18, 11}}, + {15, {11, 3, 27, 35, 12, 36, 28, 4}}, + {8, {18, 10, 34, 42, 11, 35, 27, 3}}, + {23, {1, 12, 11, 2, 19, 20}}, + {16, {1, 11, 18, 2, 26, 19}}, + {7, {19, 11, 35, 43, 20, 44, 36, 12}}, + {0, {26, 18, 42, 50, 19, 43, 35, 11}}}, + .expectedVertices = {{0, {0.0, 0.0, 0.0}}, + {3, {0.25, 0.0, 0.0}}, + {5, {0.0, 0.25, 0.0}}, + {27, {0.5, 0.0, 0.0}}, + {29, {0.0, 0.5, 0.0}}, + {1, {0.0, 0.0, 1.0}}, + {11, {0.3, 0.0, 1.0}}, + {13, {0.0, 0.3, 1.0}}, + {35, {0.6, 0.0, 1.0}}, + {37, {0.0, 0.6, 1.0}}, + {2, {0.0, 0.0, 2.0}}, + {19, {0.35, 0.0, 2.0}}, + {21, {0.0, 0.35, 2.0}}, + {43, {0.7, 0.0, 2.0}}, + {45, {0.0, 0.7, 2.0}}}, + .expectedBoundaries = {{ablate::domain::Region::ENTIREDOMAIN, {{1, 11, 12}}}, + {std::make_shared("outerShell"), {{28, 27, 36, 35}, {44, 43, 36, 35}, {48, 49, 30, 41}}}, + {std::make_shared("lowerCap"), {{0, 3, 4}, {0, 10, 3}, {3, 4, 28, 27}, {10, 3, 27, 34}}}, + {std::make_shared("upperCap"), {{2, 19, 20}, {2, 19, 26}, {26, 50, 43, 19}}}}})); diff --git a/tests/unitTests/utilities/CMakeLists.txt b/tests/unitTests/utilities/CMakeLists.txt index 793baa3df..61487aa0d 100644 --- a/tests/unitTests/utilities/CMakeLists.txt +++ b/tests/unitTests/utilities/CMakeLists.txt @@ -4,4 +4,5 @@ target_sources(ablateUnitTestLibrary petscUtilitiesTests.cpp petscSupportTests.cpp stringUtilitiesTests.cpp + vectorUtilitiesTests.cpp ) diff --git a/tests/unitTests/utilities/vectorUtilitiesTests.cpp b/tests/unitTests/utilities/vectorUtilitiesTests.cpp new file mode 100644 index 000000000..b1e5090d2 --- /dev/null +++ b/tests/unitTests/utilities/vectorUtilitiesTests.cpp @@ -0,0 +1,40 @@ +#include "gtest/gtest.h" +#include "utilities/vectorUtilities.hpp" + +TEST(VectorUtilitiesTests, ShouldConvertVectorsToArraysWithDefaults) { + { + auto actual = ablate::utilities::VectorUtilities::ToArray({1.0, 2.0, 3.0}, 0.0); + auto expected = std::array{1.0, 2.0, 3.0}; + ASSERT_EQ(expected, actual); + } + { + auto actual = ablate::utilities::VectorUtilities::ToArray({1.0, 2.0}, 0.0); + auto expected = std::array{1.0, 2.0, 0.0}; + ASSERT_EQ(expected, actual); + } + { + auto actual = ablate::utilities::VectorUtilities::ToArray({1.0, 2.0}, 100.0); + auto expected = std::array{1.0, 2.0, 100.0}; + ASSERT_EQ(expected, actual); + } +} + +TEST(VectorUtilitiesTests, ShouldConvertVectorsToArrays) { + { + auto actual = ablate::utilities::VectorUtilities::ToArray({1.0, 2.0, 3.0}); + auto expected = std::array{1.0, 2.0, 3.0}; + ASSERT_EQ(expected, actual); + } + { + auto actual = ablate::utilities::VectorUtilities::ToArray({1.0, 2.0}); + auto expected = std::array{1.0, 2.0}; + ASSERT_EQ(expected, actual); + } +} + +TEST(VectorUtilitiesTests, ShouldConvertThrowExceptionWhenTooLarge) { ASSERT_ANY_THROW((ablate::utilities::VectorUtilities::ToArray({1.0, 2.0, 3.0, 4.0}, 0.0))); } + +TEST(VectorUtilitiesTests, ShouldConvertThrowExceptionWhenWrongSize) { + ASSERT_ANY_THROW((ablate::utilities::VectorUtilities::ToArray({1.0, 2.0, 3.0, 4.0}))) << "When the vector to too large"; + ASSERT_ANY_THROW((ablate::utilities::VectorUtilities::ToArray({1.0, 2.0, 3.0, 4.0}))) << "When the vector to too small"; +} \ No newline at end of file