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

@@ -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)

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)

src/domain/descriptions/CMakeLists.txt

@@ -0,0 +1,8 @@
+target_sources(ablateLibrary
+        PRIVATE
+        axisymmetric.cpp
+
+        PUBLIC
+        meshDescription.hpp
+        axisymmetric.hpp
+        )

src/domain/descriptions/axisymmetric.cpp

@@ -0,0 +1,170 @@
+#include "axisymmetric.hpp"
+
+#include <utility>
+#include "utilities/vectorUtilities.hpp"
+
+ablate::domain::descriptions::Axisymmetric::Axisymmetric(const std::vector<PetscReal> &startLocation, PetscReal length, std::shared_ptr<ablate::mathFunctions::MathFunction> radiusFunction,
+                                                         PetscInt numberWedges, PetscInt numberSlices, PetscInt numberShells)
+    : startLocation(utilities::VectorUtilities::ToArray<PetscReal, 3>(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::Region> ablate::domain::descriptions::Axisymmetric::GetRegion(const std::set<PetscInt> &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<double>, "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"));

src/domain/descriptions/axisymmetric.hpp

@@ -0,0 +1,132 @@
+#ifndef ABLATELIBRARY_AXISYMMETRIC_HPP
+#define ABLATELIBRARY_AXISYMMETRIC_HPP
+
+#include <array>
+#include <memory>
+#include <vector>
+#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<PetscReal, 3> 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<ablate::mathFunctions::MathFunction> 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<ablate::domain::Region> shellBoundary = std::make_shared<ablate::domain::Region>("outerShell");
+    const static inline std::shared_ptr<ablate::domain::Region> lowerCapBoundary = std::make_shared<ablate::domain::Region>("lowerCap");
+    const static inline std::shared_ptr<ablate::domain::Region> upperCapBoundary = std::make_shared<ablate::domain::Region>("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<PetscReal>& startLocation, PetscReal length, std::shared_ptr<ablate::mathFunctions::MathFunction> 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<ablate::domain::Region> GetBoundaryRegion() const override { return std::make_shared<ablate::domain::Region>("boundary"); }
+
+    /**
+     * returns the boundary region for the end caps and outer shell
+     * @param face
+     * @return
+     */
+    [[nodiscard]] std::shared_ptr<ablate::domain::Region> GetRegion(const std::set<PetscInt>& face) const override;
+};
+}  // namespace ablate::domain::descriptions
+#endif  // ABLATELIBRARY_AXISYMMETRIC_HPP