Add min/max functions: Vh*Vh->Vh, R*Vh->Vh, Vh*R->Vh, Vh->R

It is required that Vh has scalar values.

src/language/modules/MathFunctionRegisterForVh.cpp

@@ -197,4 +197,58 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
                            const IDiscreteFunction>(const TinyVector<3>&, std::shared_ptr<const IDiscreteFunction>)>>(
                            [](const TinyVector<3>& a, std::shared_ptr<const IDiscreteFunction> b)
                              -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a) -> double { return min(a); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>,
+                                                                    std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(double, std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](double a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, double)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a,
+                              double b) -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a) -> double { return max(a); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>,
+                                                                    std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(double, std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](double a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, double)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a,
+                              double b) -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); }));
 }

src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp

@@ -111,7 +111,7 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<c
     }
   } else {
     std::stringstream os;
-    os << "incompatible operands type " << operand_type_name(f) << " and " << operand_type_name(g);
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g);
     throw NormalError(os.str());
   }
 }
@@ -139,7 +139,7 @@ atan2(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
     }
   } else {
     std::stringstream os;
-    os << "incompatible operands type " << operand_type_name(a) << " and " << operand_type_name(f);
+    os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f);
     throw NormalError(os.str());
   }
 }
@@ -167,7 +167,7 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
     }
   } else {
     std::stringstream os;
-    os << "incompatible operands type " << operand_type_name(f) << " and " << operand_type_name(a);
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a);
     throw NormalError(os.str());
   }
 }
@@ -479,3 +479,243 @@ template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<2>&,
 
 template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<3>&,
                                                       const std::shared_ptr<const IDiscreteFunction>&);
+
+double
+min(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return min(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return min(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return min(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw NormalError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+min(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
+{
+  if ((f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+      (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
+    std::shared_ptr mesh = getCommonMesh({f, g});
+
+    if (mesh.use_count() == 0) {
+      throw NormalError("operands are defined on different meshes");
+    }
+
+    switch (mesh->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g);
+    throw NormalError(os.str());
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+min(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f);
+    throw NormalError(os.str());
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+min(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a);
+    throw NormalError(os.str());
+  }
+}
+
+double
+max(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return max(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return max(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return max(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw NormalError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+max(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
+{
+  if ((f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+      (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
+    std::shared_ptr mesh = getCommonMesh({f, g});
+
+    if (mesh.use_count() == 0) {
+      throw NormalError("operands are defined on different meshes");
+    }
+
+    switch (mesh->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g);
+    throw NormalError(os.str());
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+max(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f);
+    throw NormalError(os.str());
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+max(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a);
+    throw NormalError(os.str());
+  }
+}

src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp

@@ -65,4 +65,22 @@ template <size_t VectorDimension>
 std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&,
                                              const std::shared_ptr<const IDiscreteFunction>&);
 
+double min(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> min(const std::shared_ptr<const IDiscreteFunction>&,
+                                             const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> min(const double, const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> min(const std::shared_ptr<const IDiscreteFunction>&, const double);
+
+double max(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> max(const std::shared_ptr<const IDiscreteFunction>&,
+                                             const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> max(const double, const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> max(const std::shared_ptr<const IDiscreteFunction>&, const double);
+
 #endif   // EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP

src/scheme/DiscreteFunctionP0.hpp

@@ -4,6 +4,7 @@
 #include <scheme/IDiscreteFunction.hpp>
 
 #include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -559,6 +560,98 @@ class DiscreteFunctionP0 : public IDiscreteFunction
     return result;
   }
 
+  PUGS_INLINE friend double
+  min(const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+
+    return min(f.m_cell_values);
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  min(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
+    Assert(f.m_mesh == g.m_mesh);
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(f[cell_id], g[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  min(const double a, const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(a, f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  min(const DiscreteFunctionP0& f, const double a)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(f[cell_id], a); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend double
+  max(const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+
+    return max(f.m_cell_values);
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  max(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
+    Assert(f.m_mesh == g.m_mesh);
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(f[cell_id], g[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  max(const double a, const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(a, f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  max(const DiscreteFunctionP0& f, const double a)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(f[cell_id], a); });
+
+    return result;
+  }
+
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()} {}
 
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value)