Generalized Axis specification for AxiSymmetric Mesh (#519)

* abstracted out the 1D mesh used for axisymmetric

* added a fixed point input

* cleaned up the mesh mapper

* renamed the axis generator stuff

* working on regions and mesh mapping

* bug fix

* version bump

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.24)
+project(ablateLibrary VERSION 0.12.25)
 
 # Load the Required 3rd Party Libaries
 pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc)

src/domain/descriptions/CMakeLists.txt

@@ -1,8 +1,13 @@
 target_sources(ablateLibrary
         PRIVATE
         axisymmetric.cpp
+        generatedAxis.cpp
+        specifiedAxis.cpp
 
         PUBLIC
         meshDescription.hpp
         axisymmetric.hpp
+        axisDescription.hpp
+        generatedAxis.hpp
+        specifiedAxis.hpp
         )

src/domain/descriptions/axisDescription.hpp

@@ -0,0 +1,29 @@
+#ifndef ABLATELIBRARY_AXISDESCRIPTION_HPP
+#define ABLATELIBRARY_AXISDESCRIPTION_HPP
+
+#include <petsc.h>
+#include <domain/region.hpp>
+#include <set>
+
+namespace ablate::domain::descriptions {
+/**
+ * A simple interface that describes points a long a line on axis
+ */
+class AxisDescription {
+   public:
+    virtual ~AxisDescription() = default;
+    /**
+     * Total number of nodes/vertices in the entire mesh
+     * @return
+     */
+    [[nodiscard]] virtual const PetscInt& GetNumberVertices() const = 0;
+
+    /**
+     * Builds the node coordinate for each vertex
+     * @return
+     */
+    virtual void SetCoordinate(PetscInt node, PetscReal coordinate[3]) const = 0;
+};
+}  // namespace ablate::domain::descriptions
+
+#endif  // ABLATELIBRARY_AXISDESCRIPTION_HPP

src/domain/descriptions/axisymmetric.cpp

@@ -3,13 +3,12 @@
 #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),
+ablate::domain::descriptions::Axisymmetric::Axisymmetric(std::shared_ptr<ablate::domain::descriptions::AxisDescription> axis, std::shared_ptr<ablate::mathFunctions::MathFunction> radiusFunction,
+                                                         PetscInt numberWedges, PetscInt numberShells)
+    : axisDescription(std::move(axis)),
       radiusFunction(std::move(radiusFunction)),
       numberWedges(numberWedges),
-      numberSlices(numberSlices),
+      numberSlices(axisDescription->GetNumberVertices() - 1),
       numberShells(numberShells),
       numberCellsPerSlice(numberWedges * numberShells),
       numberCellsPerShell(numberWedges * numberSlices),
@@ -84,11 +83,6 @@ void ablate::domain::descriptions::Axisymmetric::BuildTopology(PetscInt cell, Pe
     }
 }
 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) {
@@ -101,10 +95,8 @@ void ablate::domain::descriptions::Axisymmetric::SetCoordinate(PetscInt node, Pe
         nodeRotationIndex = (node - numberCenterVertices) % numberWedges;
     }
 
-    // offset along z
-    auto delta = length / numberSlices;
-    auto offset = delta * nodeSlice;
-    coordinate[2] += offset;
+    // Get the coordinate along the axis.
+    axisDescription->SetCoordinate(nodeSlice, coordinate);
 
     // compute the maximum radius at this coordinate
     auto radius = radiusFunction->Eval(coordinate, 3, NAN);
@@ -165,6 +157,6 @@ std::shared_ptr<ablate::domain::Region> ablate::domain::descriptions::Axisymmetr
 
 #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::domain::descriptions::AxisDescription, "axis", "describes the nodes along the z axis"),
          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"));
+         ARG(int, "numberShells", "slicing of the cylinder along the radius"));

src/domain/descriptions/axisymmetric.hpp

@@ -4,6 +4,7 @@
 #include <array>
 #include <memory>
 #include <vector>
+#include "axisDescription.hpp"
 #include "mathFunctions/mathFunction.hpp"
 #include "meshDescription.hpp"
 
@@ -14,9 +15,8 @@ namespace ablate::domain::descriptions {
  */
 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
+    //! Store the start, end, and nodes in the axis of the mesh
+    const std::shared_ptr<ablate::domain::descriptions::AxisDescription> axisDescription;
 
     //! function used to describe a single return value (radius) as a functino of z
     const std::shared_ptr<ablate::mathFunctions::MathFunction> radiusFunction;
@@ -64,19 +64,17 @@ class Axisymmetric : public ablate::domain::descriptions::MeshDescription {
      * @param in
      * @return
      */
-    inline PetscInt CellReverser(PetscInt in) const { return numberCells - in - 1; }
+    [[nodiscard]] 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 axis describes the mesh along the z axis, must be 3D
      * @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,
+    Axisymmetric(std::shared_ptr<ablate::domain::descriptions::AxisDescription> axis, std::shared_ptr<ablate::mathFunctions::MathFunction> radiusFunction, PetscInt numberWedges,
                  PetscInt numberShells);
 
     /**

src/domain/descriptions/generatedAxis.cpp

@@ -0,0 +1,28 @@
+#include "generatedAxis.hpp"
+
+#include "utilities/vectorUtilities.hpp"
+
+ablate::domain::descriptions::GeneratedAxis::GeneratedAxis(const std::vector<PetscReal> &startLocation, PetscReal length, PetscInt numberNodes)
+    : startLocation(utilities::VectorUtilities::ToArray<PetscReal, 3>(startLocation)), length(length), numberNodes(numberNodes) {
+    // make sure there are at least 2 nodes
+    if (this->numberNodes < 2) {
+        throw std::invalid_argument("GeneratedAxis requires at least 2 nodes");
+    }
+}
+
+void ablate::domain::descriptions::GeneratedAxis::SetCoordinate(PetscInt node, PetscReal coordinate[3]) const {
+    // start the coordinate at the start location
+    coordinate[0] = startLocation[0];
+    coordinate[1] = startLocation[1];
+    coordinate[2] = startLocation[2];
+
+    // offset along z
+    auto delta = length / (numberNodes - 1);
+    auto offset = delta * node;
+    coordinate[2] += offset;
+}
+
+#include "registrar.hpp"
+REGISTER_DEFAULT(ablate::domain::descriptions::AxisDescription, ablate::domain::descriptions::GeneratedAxis, "Creates a simple fixed spacing along 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(int, "nodes", "number of nodes along the z-axis"));

src/domain/descriptions/generatedAxis.hpp

@@ -0,0 +1,48 @@
+#ifndef ABLATELIBRARY_GENERATEDAXIS_HPP
+#define ABLATELIBRARY_GENERATEDAXIS_HPP
+
+#include <array>
+#include <memory>
+#include <vector>
+#include "axisDescription.hpp"
+#include "mathFunctions/mathFunction.hpp"
+
+namespace ablate::domain::descriptions {
+
+/**
+ * Describes a simple strait line along the z axis
+ */
+class GeneratedAxis : public ablate::domain::descriptions::AxisDescription {
+   private:
+    //! Store the start and end location of the mesh
+    const std::array<PetscReal, 3> startLocation;
+
+    //! the length of the domain
+    const PetscReal length;  // this is in z
+
+    //! the number of nodes in the axis
+    const PetscInt numberNodes;
+
+   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 numberNodes the number of nodes along the axis
+     */
+    GeneratedAxis(const std::vector<PetscReal>& startLocation, PetscReal length, PetscInt numberNodes);
+
+    /**
+     * Total number of nodes/vertices in the entire mesh
+     * @return
+     */
+    [[nodiscard]] const PetscInt& GetNumberVertices() const override { return numberNodes; }
+
+    /**
+     * Sets the axis coordinate for this node
+     * @return
+     */
+    void SetCoordinate(PetscInt node, PetscReal coordinate[3]) const override;
+};
+}  // namespace ablate::domain::descriptions
+#endif  // ABLATELIBRARY_AXISYMMETRIC_HPP

src/domain/descriptions/specifiedAxis.cpp

@@ -0,0 +1,31 @@
+#include "specifiedAxis.hpp"
+
+#include <algorithm>
+#include "utilities/vectorUtilities.hpp"
+
+ablate::domain::descriptions::SpecifiedAxis::SpecifiedAxis(const std::vector<PetscReal>& nodeOffsets, const std::vector<PetscReal>& startLocation)
+    : startLocation(startLocation.empty() ? std::array<PetscReal, 3>{0.0, 0.0, 0.0} : utilities::VectorUtilities::ToArray<PetscReal, 3>(startLocation)),
+      zOffset(nodeOffsets),
+      numberNodes((PetscInt)nodeOffsets.size()) {
+    // make sure there are at least 2 nodes
+    if (this->numberNodes < 2) {
+        throw std::invalid_argument("SpecifiedAxis requires at least 2 nodes");
+    }
+    if (!std::is_sorted(zOffset.begin(), zOffset.end())) {
+        throw std::invalid_argument("SpecifiedAxis requires zOffset vector be sorted from smallest to largest");
+    }
+}
+
+void ablate::domain::descriptions::SpecifiedAxis::SetCoordinate(PetscInt node, PetscReal coordinate[3]) const {
+    // start the coordinate at the start location
+    coordinate[0] = startLocation[0];
+    coordinate[1] = startLocation[1];
+    coordinate[2] = startLocation[2];
+
+    // offset along z
+    coordinate[2] += zOffset[node];
+}
+
+#include "registrar.hpp"
+REGISTER(ablate::domain::descriptions::AxisDescription, ablate::domain::descriptions::SpecifiedAxis, "Describes a simple strait line along the z axis with specified offsets",
+         ARG(std::vector<double>, "offsets", "an ordered list of offsets in z"), OPT(std::vector<double>, "start", "the start coordinate of the mesh, must be 3D"));

src/domain/descriptions/specifiedAxis.hpp

@@ -0,0 +1,47 @@
+#ifndef ABLATELIBRARY_SPECIFIEDAXIS_HPP
+#define ABLATELIBRARY_SPECIFIEDAXIS_HPP
+
+#include <array>
+#include <memory>
+#include <vector>
+#include "axisDescription.hpp"
+#include "mathFunctions/mathFunction.hpp"
+
+namespace ablate::domain::descriptions {
+
+/**
+ * Describes a simple strait line along the z axis with specified offsets
+ */
+class SpecifiedAxis : public ablate::domain::descriptions::AxisDescription {
+   private:
+    //! Store the start and end location of the mesh
+    const std::array<PetscReal, 3> startLocation;
+
+    //! Store the z offset from each node
+    const std::vector<PetscReal> zOffset;
+
+    //! store the number of vertices for quick lookup
+    const PetscInt numberNodes;
+
+   public:
+    /**
+     * read in and compute the axis from a list of z offsets and optional start location
+     * @param nodeOffsets z offsets from the start location
+     * @param startLocation optional start location (default is {0, 0, 0})
+     */
+    explicit SpecifiedAxis(const std::vector<PetscReal>& nodeOffsets, const std::vector<PetscReal>& startLocation = {});
+
+    /**
+     * Total number of nodes/vertices in the entire mesh
+     * @return
+     */
+    [[nodiscard]] const PetscInt& GetNumberVertices() const override { return numberNodes; }
+
+    /**
+     * Sets the axis coordinate for this node
+     * @return
+     */
+    void SetCoordinate(PetscInt node, PetscReal coordinate[3]) const override;
+};
+}  // namespace ablate::domain::descriptions
+#endif  // ABLATELIBRARY_SPECIFIEDAXIS_HPP

src/domain/modifiers/edgeClusteringMapper.cpp

@@ -1,8 +1,10 @@
 #include "edgeClusteringMapper.hpp"
+
+#include <utility>
 #include "mathFunctions/functionFactory.hpp"
 
-ablate::domain::modifiers::EdgeClusteringMapper::EdgeClusteringMapper(int direction, double startIn, double end, double beta)
-    : ablate::domain::modifiers::MeshMapper(mathFunctions::Create(MappingFunction, this)), direction(direction), start(startIn), size(end - startIn), beta(beta) {
+ablate::domain::modifiers::EdgeClusteringMapper::EdgeClusteringMapper(int direction, double startIn, double end, double beta, std::shared_ptr<ablate::domain::Region> mappingRegion)
+    : ablate::domain::modifiers::MeshMapper(mathFunctions::Create(MappingFunction, this), std::move(mappingRegion)), direction(direction), start(startIn), size(end - startIn), beta(beta) {
     // make sure that the direction is valid
     if (direction < 0 || direction > 2) {
         throw std::invalid_argument("The direction must be 0, 1, or 2. Direction " + std::to_string(direction) + " is invalid.");
@@ -39,4 +41,5 @@ REGISTER(ablate::domain::modifiers::Modifier, ablate::domain::modifiers::EdgeClu
          "Performs clustering mapping using an algebraic relationship at the edges of the domain using Equation 9-42 from Hoffmann, Klaus A., and Steve T. Chiang. \"Computational fluid dynamics "
          "volume I. Forth Edition\" Engineering education system (2000).",
          ARG(int, "direction", "The direction (0, 1, 2) to perform the mapping"), ARG(double, "start", "The start of the domain in direction"),
-         ARG(double, "end", "The end of the domain in direction"), ARG(double, "beta", "The clustering factor (1 -> infinity, where 1 is more clustering)"));
+         ARG(double, "end", "The end of the domain in direction"), ARG(double, "beta", "The clustering factor (1 -> infinity, where 1 is more clustering)"),
+         OPT(ablate::domain::Region, "region", "optional region to apply this mapper (Default is everywhere)"));

src/domain/modifiers/edgeClusteringMapper.hpp

@@ -28,14 +28,15 @@ class EdgeClusteringMapper : public MeshMapper {
      * @param start The start of the domain in direction
      * @param end The end of the domain in direction
      * @param beta The clustering factor
+     * @param mappingRegion optional region to apply this mapper.  Default is everywhere
      */
-    explicit EdgeClusteringMapper(int direction, double start, double end, double beta);
+    explicit EdgeClusteringMapper(int direction, double start, double end, double beta, std::shared_ptr<ablate::domain::Region> mappingRegion = {});
 
     /**
      * Provide name of modifier for debug/output
      * @return
      */
-    std::string ToString() const override;
+    [[nodiscard]] std::string ToString() const override;
 
    private:
     static PetscErrorCode MappingFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar* u, void* ctx);