From 05a2f2ed3783ff0f34acbc4162eb28e46c652a30 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Wed, 9 Feb 2022 10:42:33 +0100
Subject: [PATCH] [MeL] Bulk element mapping uses non-template points

---
 MeshLib/Elements/MapBulkElementPoint.cpp | 102 +++++++++++++++++++----
 MeshLib/Elements/MapBulkElementPoint.h   |  47 +++++------
 2 files changed, 103 insertions(+), 46 deletions(-)

diff --git a/MeshLib/Elements/MapBulkElementPoint.cpp b/MeshLib/Elements/MapBulkElementPoint.cpp
index e996414b05c..5bf83b89b01 100644
--- a/MeshLib/Elements/MapBulkElementPoint.cpp
+++ b/MeshLib/Elements/MapBulkElementPoint.cpp
@@ -14,9 +14,24 @@
 
 namespace MeshLib
 {
+MathLib::Point3d getBulkElementPoint(Line const& /*line*/,
+                                     std::size_t const face_id,
+                                     MathLib::WeightedPoint const& /*wp*/)
+{
+    switch (face_id)
+    {
+        case 0:
+            return MathLib::Point3d{std::array<double, 3>{{-1.0, 0.0, 0.0}}};
+        case 1:
+            return MathLib::Point3d{std::array<double, 3>{{1.0, 0.0, 0.0}}};
+        default:
+            OGS_FATAL("Invalid face id '{:d}' for a line.", face_id);
+    }
+}
+
 MathLib::Point3d getBulkElementPoint(Tri const& /*tri*/,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint1D const& wp)
+                                     MathLib::WeightedPoint const& wp)
 {
     switch (face_id)
     {
@@ -35,7 +50,7 @@ MathLib::Point3d getBulkElementPoint(Tri const& /*tri*/,
 
 MathLib::Point3d getBulkElementPoint(Quad const& /*quad*/,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint1D const& wp)
+                                     MathLib::WeightedPoint const& wp)
 {
     switch (face_id)
     {
@@ -56,7 +71,7 @@ MathLib::Point3d getBulkElementPoint(Quad const& /*quad*/,
 
 MathLib::Point3d getBulkElementPoint(Hex const& /*hex*/,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp)
+                                     MathLib::WeightedPoint const& wp)
 {
     switch (face_id)
     {
@@ -83,7 +98,7 @@ MathLib::Point3d getBulkElementPoint(Hex const& /*hex*/,
 
 MathLib::Point3d getBulkElementPoint(Prism const& /*prism*/,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp)
+                                     MathLib::WeightedPoint const& wp)
 {
     switch (face_id)
     {
@@ -108,7 +123,7 @@ MathLib::Point3d getBulkElementPoint(Prism const& /*prism*/,
 
 MathLib::Point3d getBulkElementPoint(Pyramid const& /*pyramid*/,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp)
+                                     MathLib::WeightedPoint const& wp)
 {
     switch (face_id)
     {
@@ -134,7 +149,7 @@ MathLib::Point3d getBulkElementPoint(Pyramid const& /*pyramid*/,
 
 MathLib::Point3d getBulkElementPoint(Tet const& /*tet*/,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp)
+                                     MathLib::WeightedPoint const& wp)
 {
     switch (face_id)
     {
@@ -152,11 +167,32 @@ MathLib::Point3d getBulkElementPoint(Tet const& /*tet*/,
     }
 }
 
-MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
-                                     std::size_t const bulk_element_id,
-                                     std::size_t const bulk_face_id,
-                                     MathLib::WeightedPoint1D const& wp)
+}  // namespace MeshLib
+
+static MathLib::Point3d getBulkElementPoint0D(MeshLib::Mesh const& mesh,
+                                              std::size_t const bulk_element_id,
+                                              std::size_t const bulk_face_id,
+                                              MathLib::WeightedPoint const& wp)
+{
+    using namespace MeshLib;
+
+    auto const* element = mesh.getElement(bulk_element_id);
+    if (element->getCellType() == CellType::LINE2)
+    {
+        Line const& line(*dynamic_cast<Line const*>(element));
+        return getBulkElementPoint(line, bulk_face_id, wp);
+    }
+    OGS_FATAL("Wrong cell type '{:s}' or functionality not yet implemented.",
+              CellType2String(element->getCellType()));
+}
+
+static MathLib::Point3d getBulkElementPoint1D(MeshLib::Mesh const& mesh,
+                                              std::size_t const bulk_element_id,
+                                              std::size_t const bulk_face_id,
+                                              MathLib::WeightedPoint const& wp)
 {
+    using namespace MeshLib;
+
     auto const* element = mesh.getElement(bulk_element_id);
     if (element->getCellType() == CellType::QUAD4)
     {
@@ -172,11 +208,13 @@ MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
               CellType2String(element->getCellType()));
 }
 
-MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
-                                     std::size_t bulk_element_id,
-                                     std::size_t bulk_face_id,
-                                     MathLib::WeightedPoint2D const& wp)
+static MathLib::Point3d getBulkElementPoint2D(MeshLib::Mesh const& mesh,
+                                              std::size_t bulk_element_id,
+                                              std::size_t bulk_face_id,
+                                              MathLib::WeightedPoint const& wp)
 {
+    using namespace MeshLib;
+
     auto const* element = mesh.getElement(bulk_element_id);
     if (element->getCellType() == CellType::HEX8)
     {
@@ -203,11 +241,39 @@ MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
 }
 
 // TODO disable the 3d elements in the local assembler creator
-MathLib::Point3d getBulkElementPoint(Mesh const& /*mesh*/,
-                                     std::size_t /*bulk_element_id*/,
-                                     std::size_t /*bulk_face_id*/,
-                                     MathLib::WeightedPoint3D const& /*wp*/)
+static MathLib::Point3d getBulkElementPoint3D(
+    MeshLib::Mesh const& /*mesh*/,
+    std::size_t /*bulk_element_id*/,
+    std::size_t /*bulk_face_id*/,
+    MathLib::WeightedPoint const& /*wp*/)
 {
     return MathLib::ORIGIN;
 }
+
+namespace MeshLib
+{
+MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
+                                     std::size_t const bulk_element_id,
+                                     std::size_t const bulk_face_id,
+                                     MathLib::WeightedPoint const& wp)
+{
+    switch (wp.getDimension())
+    {
+        case 0:
+            return getBulkElementPoint0D(
+                mesh, bulk_element_id, bulk_face_id, wp);
+        case 1:
+            return getBulkElementPoint1D(
+                mesh, bulk_element_id, bulk_face_id, wp);
+        case 2:
+            return getBulkElementPoint2D(
+                mesh, bulk_element_id, bulk_face_id, wp);
+        case 3:
+            return getBulkElementPoint3D(
+                mesh, bulk_element_id, bulk_face_id, wp);
+    }
+
+    OGS_FATAL("Weighted point dimension of {} not supported.",
+              wp.getDimension());
+}
 }  // namespace MeshLib
diff --git a/MeshLib/Elements/MapBulkElementPoint.h b/MeshLib/Elements/MapBulkElementPoint.h
index a28925ce3bc..0a0d6a6e55d 100644
--- a/MeshLib/Elements/MapBulkElementPoint.h
+++ b/MeshLib/Elements/MapBulkElementPoint.h
@@ -10,18 +10,23 @@
 
 #pragma once
 
-#include "MathLib/TemplateWeightedPoint.h"
-
-#include "MeshLib/Elements/FaceRule.h"
+#include "MathLib/WeightedPoint.h"
 #include "MeshLib/Elements/Elements.h"
 
 namespace MeshLib
 {
+
 /// \page BulkMappingDocuPage
 /// [Documentation](https://www.opengeosys.org/pdf/BulkElementMappings.pdf) for
 /// the mapping of a point given in local coordinates of a boundary face/element
 /// to the corresponding bulk element point.
 
+/// Maps the given 0d end point \c wp of a line, to the corresponding 1d point
+/// on the line. The result is given in local coordinates of the line element.
+MathLib::Point3d getBulkElementPoint(Line const& line,
+                                     std::size_t const face_id,
+                                     MathLib::WeightedPoint const& wp);
+
 /// Maps the given lower dimensional boundary point \c wp of a line, i.e. the 1d
 /// integration point given in local coordinates of a line, to higher
 /// dimensional point of the triangle face (defined by the triangle element and
@@ -32,7 +37,7 @@ namespace MeshLib
 /// \return the mapped point
 MathLib::Point3d getBulkElementPoint(Tri const& tri,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint1D const& wp);
+                                     MathLib::WeightedPoint const& wp);
 
 /// Maps the given lower dimensional boundary point \c wp of a line, i.e. the 1d
 /// integration point given in local coordinates of a line, to higher dimensional
@@ -44,55 +49,41 @@ MathLib::Point3d getBulkElementPoint(Tri const& tri,
 /// \return the mapped point
 MathLib::Point3d getBulkElementPoint(Quad const& quad,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint1D const& wp);
+                                     MathLib::WeightedPoint const& wp);
 
 /// Maps the given 2d boundary point \c wp of a quad face, given in local
 /// coordinates of the quad, to 3d point existing on a hexahedron face also in
 /// local coordinates.
 MathLib::Point3d getBulkElementPoint(Hex const& hex,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp);
+                                     MathLib::WeightedPoint const& wp);
 
-/// Maps the given lower dimensional 2d boundary point \c wp of a quad or
-/// triangle element, given in local coordinates of the quad or triangle, to a
-/// 3d point existing on a tet face also in local coordinates.
+/// Maps the given lower dimensional 2d boundary point \c wp of a triangle
+/// element, given in local coordinates of the quad or triangle, to a 3d point
+/// existing on a tet face also in local coordinates.
 MathLib::Point3d getBulkElementPoint(Tet const& tet,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp);
+                                     MathLib::WeightedPoint const& wp);
 
 /// Maps the given lower dimensional 2d boundary point \c wp of a quad or
 /// triangle element, given in local coordinates of the quad or triangle, to a
 /// 3d point existing on a prism face also in local coordinates.
 MathLib::Point3d getBulkElementPoint(Prism const& prism,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp);
+                                     MathLib::WeightedPoint const& wp);
 
 /// Maps the given lower dimensional 2d boundary point \c wp of a quad or
 /// triangle element, given in local coordinates of the quad or triangle, to a
 /// 3d point existing on a pyramid face also in local coordinates.
 MathLib::Point3d getBulkElementPoint(Pyramid const& pyramid,
                                      std::size_t const face_id,
-                                     MathLib::WeightedPoint2D const& wp);
+                                     MathLib::WeightedPoint const& wp);
 
-/// Maps the given 1d boundary point \c wp of a boundary element, given in local
+/// Maps the given boundary point \c wp of a boundary element, given in local
 /// coordinates of the boundary element, to 3d point existing on a bulk element
 /// also in local coordinates.
 MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
                                      std::size_t const bulk_element_id,
                                      std::size_t const bulk_face_id,
-                                     MathLib::WeightedPoint1D const& wp);
-
-/// Maps the given 2d boundary point \c wp of a boundary element, given in local
-/// coordinates of the boundary element, to 3d point existing on a bulk element
-/// also in local coordinates.
-MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
-                                     std::size_t bulk_element_id,
-                                     std::size_t bulk_face_id,
-                                     MathLib::WeightedPoint2D const& wp);
-
-// TODO disable the 3d elements in the local assembler creator
-MathLib::Point3d getBulkElementPoint(Mesh const& mesh,
-                                     std::size_t bulk_element_id,
-                                     std::size_t bulk_face_id,
-                                     MathLib::WeightedPoint3D const& wp);
+                                     MathLib::WeightedPoint const& wp);
 }  // namespace MeshLib
-- 
GitLab