From ec7e69538d0944d16b87537caff6c15b3cca8b16 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 12 Jun 2020 14:50:23 +0200
Subject: [PATCH] [test] Add test for cylinder integration
---
test/common/geometry/CMakeLists.txt | 1 +
.../geometry/test_cylinderintegration.cc | 168 ++++++++++++++++++
2 files changed, 169 insertions(+)
create mode 100644 test/common/geometry/test_cylinderintegration.cc
diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index 92ea708139..f1a6d093ce 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -9,6 +9,7 @@ dumux_add_test(SOURCES test_2d3d_intersection.cc LABELS unit)
dumux_add_test(SOURCES test_graham_convex_hull.cc LABELS unit)
dumux_add_test(SOURCES test_intersectingentity_cartesiangrid.cc LABELS unit)
dumux_add_test(SOURCES test_circlepoints.cc LABELS unit)
+dumux_add_test(SOURCES test_cylinderintegration.cc LABELS unit)
dune_symlink_to_source_files(FILES ball.msh)
dumux_add_test(SOURCES test_intersectionentityset.cc
diff --git a/test/common/geometry/test_cylinderintegration.cc b/test/common/geometry/test_cylinderintegration.cc
new file mode 100644
index 0000000000..f756c24dbb
--- /dev/null
+++ b/test/common/geometry/test_cylinderintegration.cc
@@ -0,0 +1,168 @@
+#include <config.h>
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/multidomain/embedded/cylinderintegration.hh>
+#include "transformation.hh"
+
+namespace Dumux {
+
+template<class Scalar>
+void checkVolumeIntegral(const Scalar volIntegralComputed, const Scalar volIntegralExact)
+{
+ std::cout << std::setprecision(7) << std::scientific
+ << "Volume integral: " << volIntegralComputed << " (computed), "
+ << volIntegralExact << " (exact)" << "\n";
+ if (Dune::FloatCmp::ne(volIntegralComputed, volIntegralExact, 1e-7))
+ DUNE_THROW(Dune::Exception, "Volume integral is not exact!");
+}
+
+template<class Scalar>
+void checkAreaIntegral(const Scalar areaIntegralComputed, const Scalar areaIntegralExact)
+{
+ std::cout << std::setprecision(7) << std::scientific
+ << "Area integral: " << areaIntegralComputed << " (computed), "
+ << areaIntegralExact << " (exact)" << "\n";
+ if (Dune::FloatCmp::ne(areaIntegralComputed, areaIntegralExact, 1e-7))
+ DUNE_THROW(Dune::Exception, "Area integral is not exact!");
+}
+
+template<class GlobalPosition>
+void checkCentroid(const GlobalPosition& centroidComputed, const GlobalPosition& centroidExact)
+{
+ std::cout << std::setprecision(7) << std::scientific
+ << "Centroid: " << centroidComputed << " (computed), "
+ << centroidExact << " (exact)" << "\n";
+ if (Dune::FloatCmp::ne(centroidComputed, centroidExact, 1e-7))
+ DUNE_THROW(Dune::Exception, "Centroid is wrong!");
+}
+
+template<class TestData>
+void checkCylinderIntegration(const TestData& data, const double characteristicLength)
+{
+ using Point = typename TestData::Point;
+
+ std::cout << std::endl;
+ EmbeddedCoupling::CylinderIntegration<double> cylinderIntegration(characteristicLength);
+ cylinderIntegration.setGeometry(data.bottomCenter, data.topCenter, data.radius, /*verbosity=*/1);
+ double volumeIntegral = 0.0;
+ Point centroid({0.0, 0.0, 0.0});
+ for (std::size_t i = 0; i < cylinderIntegration.size(); ++i)
+ {
+ volumeIntegral += cylinderIntegration.integrationElement(i);
+ centroid.axpy(cylinderIntegration.integrationElement(i), cylinderIntegration.integrationPoint(i));
+ }
+ centroid /= volumeIntegral;
+
+ checkVolumeIntegral(volumeIntegral, data.exactVolumeIntegral);
+ checkCentroid(centroid, data.exactCentroid);
+}
+
+template<class TestData, class EllipseData>
+void checkEllipticCylinderIntegration(const TestData& data, const EllipseData& ell, const double characteristicLength)
+{
+ using Point = typename TestData::Point;
+
+ std::cout << std::endl;
+ EmbeddedCoupling::EllipticCylinderIntegration<double> ellCylIntegration(characteristicLength);
+ ellCylIntegration.setGeometry(data.bottomCenter, data.topCenter, ell.firstAxis, ell.secondAxis, /*verbosity=*/1);
+ double volumeIntegral = 0.0;
+ Point centroid({0.0, 0.0, 0.0});
+ for (std::size_t i = 0; i < ellCylIntegration.size(); ++i)
+ {
+ volumeIntegral += ellCylIntegration.integrationElement(i);
+ centroid.axpy(ellCylIntegration.integrationElement(i), ellCylIntegration.integrationPoint(i));
+ }
+ centroid /= volumeIntegral;
+
+ checkVolumeIntegral(volumeIntegral, data.exactVolumeIntegral);
+ checkCentroid(centroid, data.exactCentroid);
+
+ // check top cap ellipse
+ std::cout << std::endl;
+ EmbeddedCoupling::EllipseIntegration<double> ellIntegration(characteristicLength);
+ ellIntegration.setGeometry(data.topCenter, ell.firstAxis, ell.secondAxis, /*verbosity=*/1);
+ double areaIntegral = 0.0;
+ Point topCentroid({0.0, 0.0, 0.0});
+ for (std::size_t i = 0; i < ellIntegration.size(); ++i)
+ {
+ areaIntegral += ellIntegration.integrationElement(i);
+ topCentroid.axpy(ellIntegration.integrationElement(i), ellIntegration.integrationPoint(i));
+ }
+ topCentroid /= areaIntegral;
+
+ checkAreaIntegral(areaIntegral, ell.firstAxis.two_norm()*ell.secondAxis.two_norm()*M_PI);
+ checkCentroid(topCentroid, data.topCenter);
+}
+
+} // end namespace Dumux
+
+int main(int argc, char* argv[])
+{
+ using namespace Dumux;
+ using Point = Dune::FieldVector<double, 3>;
+
+ struct TestData {
+ using Point = Dune::FieldVector<double, 3>;
+ Point bottomCenter, topCenter;
+ double radius;
+ double exactVolumeIntegral;
+ Point exactCentroid;
+ };
+
+ TestData data;
+ data.bottomCenter = Point({1.2, 2.3, 4.5});
+ data.topCenter = Point({8.2, 4.3, 6.5});
+ data.radius = 4.0;
+ data.exactVolumeIntegral = M_PI*data.radius*data.radius*(data.topCenter-data.bottomCenter).two_norm();
+ data.exactCentroid = data.bottomCenter; data.exactCentroid.axpy(0.5, data.topCenter-data.bottomCenter);
+
+ /////////////////////////////////////////////////////////////////////
+ // Determine cylinder volume and centroid by integration
+ ////////////////////////////////////////////////////////////////////
+ checkCylinderIntegration(data, 0.05);
+
+ /////////////////////////////////////////////////////////////////////
+ // Determine elliptic cylinder volume and centroid by integration
+ // Construct an perfect cylinder (a == b == radius and orthogonal caps)
+ ////////////////////////////////////////////////////////////////////
+ struct EllipseData { Point firstAxis, secondAxis; };
+ auto topNormal = data.topCenter-data.bottomCenter; topNormal /= topNormal.two_norm();
+ auto firstAxis = crossProduct(topNormal, Point({1.0, 0.0, 0.0}));
+ firstAxis /= firstAxis.two_norm();
+ auto secondAxis = crossProduct(topNormal, firstAxis);
+ // scale both axis with the same radius -> circle
+ firstAxis *= data.radius; secondAxis *= data.radius;
+ checkEllipticCylinderIntegration(data, EllipseData{firstAxis, secondAxis}, 0.05);
+
+ /////////////////////////////////////////////////////////////////////
+ // Determine elliptic cylinder volume and centroid by integration
+ // Construct an elliptic cylinder (a == 4*b and orthogonal caps)
+ ////////////////////////////////////////////////////////////////////
+
+ firstAxis *= 2; secondAxis *= 0.5;
+ checkEllipticCylinderIntegration(data, EllipseData{firstAxis, secondAxis}, 0.025);
+
+ /////////////////////////////////////////////////////////////////////
+ // Determine elliptic cylinder volume and centroid by integration
+ // Construct an elliptic cylinder (a == 2*b and tilted caps)
+ ////////////////////////////////////////////////////////////////////
+
+ // rotate about the first major axis
+ const auto rotationAngle = M_PI*0.25;
+ auto rotationAxis = firstAxis; rotationAxis /= rotationAxis.two_norm();
+ const auto rotate = make3DTransformation(1.0, Point(0.0), rotationAxis, rotationAngle, /*verbose*/false);
+ secondAxis = rotate(secondAxis);
+ // scale second axis so that the volume stays the same as for the orthonal cap case
+ secondAxis /= std::cos(rotationAngle);
+ checkEllipticCylinderIntegration(data, EllipseData{firstAxis, secondAxis}, 0.025);
+
+ return 0;
+}
--
GitLab