From 577baaea19c6193b49b0ac20d3b4da29e939b4f7 Mon Sep 17 00:00:00 2001 From: Marcus Mohr Date: Fri, 24 Jul 2020 18:04:53 +0200 Subject: [PATCH 1/4] Fixes UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET --- .../form_hyteg_manual/ShapeFunctionMacros.hpp | 60 +++++++++---------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/src/hyteg/forms/form_hyteg_manual/ShapeFunctionMacros.hpp b/src/hyteg/forms/form_hyteg_manual/ShapeFunctionMacros.hpp index 008bc6622..9cde39b28 100644 --- a/src/hyteg/forms/form_hyteg_manual/ShapeFunctionMacros.hpp +++ b/src/hyteg/forms/form_hyteg_manual/ShapeFunctionMacros.hpp @@ -171,44 +171,44 @@ #ifdef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET -#define DX1N0 -#define DX2N0 -#define DX3N0 +#undef DX1N0 +#undef DX2N0 +#undef DX3N0 -#define DX1N1 -#define DX2N1 -#define DX3N1 +#undef DX1N1 +#undef DX2N1 +#undef DX3N1 -#define DX1N2 -#define DX2N2 -#define DX3N2 +#undef DX1N2 +#undef DX2N2 +#undef DX3N2 -#define DX1N3 -#define DX2N3 -#define DX3N3 +#undef DX1N3 +#undef DX2N3 +#undef DX3N3 -#define DX1N4 -#define DX2N4 -#define DX3N4 +#undef DX1N4 +#undef DX2N4 +#undef DX3N4 -#define DX1N5 -#define DX2N5 -#define DX3N5 +#undef DX1N5 +#undef DX2N5 +#undef DX3N5 -#define DX1N6 -#define DX2N6 -#define DX3N6 +#undef DX1N6 +#undef DX2N6 +#undef DX3N6 -#define DX1N7 -#define DX2N7 -#define DX3N7 +#undef DX1N7 +#undef DX2N7 +#undef DX3N7 -#define DX1N8 -#define DX2N8 -#define DX3N8 +#undef DX1N8 +#undef DX2N8 +#undef DX3N8 -#define DX1N9 -#define DX2N9 -#define DX3N9 +#undef DX1N9 +#undef DX2N9 +#undef DX3N9 #endif -- GitLab From 76e2c1b92fe3b941e42f921c2ac1649cff69e801 Mon Sep 17 00:00:00 2001 From: Marcus Mohr Date: Mon, 27 Jul 2020 16:46:17 +0200 Subject: [PATCH 2/4] Implements P2ToP1Form_div without blending for 3D. This is another step towards the implementation including blending. We test that the results we get from our HyTeG form match those from the FEniCS form. Next step will be to include the blending. --- .../forms/form_hyteg_manual/P2ToP1FormDiv.hpp | 388 ++++++++++++++++-- tests/hyteg/forms/HytegVsFenicsFormTest.cpp | 21 +- 2 files changed, 383 insertions(+), 26 deletions(-) diff --git a/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp b/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp index ed0339ee6..366110339 100644 --- a/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp +++ b/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp @@ -39,24 +39,26 @@ class P2ToP1Form_div : public P2ToP1FormHyTeG }; }; -// -------- + +// ======== // div_x -// -------- +// ======== template <> class P2ToP1Form_div< 0 > : public P2ToP1FormHyTeG { public: + + // --------- + // 2D Case + // --------- void integrateAll( const std::array< Point3D, 3 >& coords, Matrixr< 3, 6 >& elMat ) const final { -// Derivatives of P2 shape functions on unit triangle +// Set shape functions and their derivatives #define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE -#include "ShapeFunctionMacros.hpp" -#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE - -// P1 shape functions on unit triangle #define DEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #include "ShapeFunctionMacros.hpp" #undef DEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE // Select quadrature rule #define QUADPOINTS quadrature::D5_points @@ -132,48 +134,143 @@ class P2ToP1Form_div< 0 > : public P2ToP1FormHyTeG INTEGRATE2D( 2, 5 ); #define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE -#include "ShapeFunctionMacros.hpp" -#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE - #define UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #include "ShapeFunctionMacros.hpp" #undef UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE #undef INTEGRATE2D }; + // --------- + // 3D Case + // --------- void integrateAll( const std::array< Point3D, 4 >& coords, Matrixr< 4, 10 >& elMat ) const final { -#define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET -#include "ShapeFunctionMacros.hpp" -#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +// Set shape functions and their derivatives +#define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET #define DEFINE_P1_SHAPE_FUNCTIONS_TET #include "ShapeFunctionMacros.hpp" #undef DEFINE_P1_SHAPE_FUNCTIONS_TET +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET - WALBERLA_ABORT( "P2ToP1Form_div not implemented for 3D, yet." ); +// Select cubature rule +#define CUBAPOINTS cubature::T4_points +#define CUBAWEIGHTS cubature::T4_weights -#define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET -#include "ShapeFunctionMacros.hpp" -#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +// Executing quadrature rule +#define INTEGRATE3D( i, j ) \ + elMat( i, j ) = 0.0; \ + for ( uint_t k = 0; k < CUBAWEIGHTS.size(); k++ ) \ + { \ + real_t L2 = CUBAPOINTS[k][0]; \ + real_t L3 = CUBAPOINTS[k][1]; \ + real_t L4 = CUBAPOINTS[k][2]; \ + real_t L1 = 1.0 - L2 - L3 - L4; \ + \ + real_t dy = DX1N ## j*(tmp12 + tmp13 - tmp14 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2]) + DX2N ## j*(-tmp13 + tmp14 + tmp15 - tmp16 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2]) + DX3N ## j*(-tmp12 - tmp15 + tmp16 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2]); \ + \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * dy * SF_M##i; \ + } + // mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + // mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + // mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + + // Quadrature point mapped to computational triangle + Point3D mappedPt; + + // pre-compute values including Jacobian determinant of inverse pull-back mapping + real_t tmp0 = -coords[0][0]; + real_t tmp1 = tmp0 + coords[1][0]; + real_t tmp2 = -coords[0][1]; + real_t tmp3 = tmp2 + coords[2][1]; + real_t tmp4 = -coords[0][2]; + real_t tmp5 = tmp4 + coords[3][2]; + real_t tmp6 = tmp0 + coords[2][0]; + real_t tmp7 = tmp2 + coords[3][1]; + real_t tmp8 = tmp4 + coords[1][2]; + real_t tmp9 = tmp0 + coords[3][0]; + real_t tmp10 = tmp2 + coords[1][1]; + real_t tmp11 = tmp4 + coords[2][2]; + real_t tmp12 = coords[0][1]*coords[2][2]; + real_t tmp13 = coords[3][1]*coords[0][2]; + real_t tmp14 = coords[0][1]*coords[3][2]; + real_t tmp15 = coords[1][1]*coords[0][2]; + real_t tmp16 = coords[0][1]*coords[1][2]; + + real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; + real_t detPfac = std::abs( detP ) / detP; + + INTEGRATE3D( 0, 0 ); + INTEGRATE3D( 0, 1 ); + INTEGRATE3D( 0, 2 ); + INTEGRATE3D( 0, 3 ); + INTEGRATE3D( 0, 4 ); + INTEGRATE3D( 0, 5 ); + INTEGRATE3D( 0, 6 ); + INTEGRATE3D( 0, 7 ); + INTEGRATE3D( 0, 8 ); + INTEGRATE3D( 0, 9 ); + + INTEGRATE3D( 1, 0 ); + INTEGRATE3D( 1, 1 ); + INTEGRATE3D( 1, 2 ); + INTEGRATE3D( 1, 3 ); + INTEGRATE3D( 1, 4 ); + INTEGRATE3D( 1, 5 ); + INTEGRATE3D( 1, 6 ); + INTEGRATE3D( 1, 7 ); + INTEGRATE3D( 1, 8 ); + INTEGRATE3D( 1, 9 ); + + INTEGRATE3D( 2, 0 ); + INTEGRATE3D( 2, 1 ); + INTEGRATE3D( 2, 2 ); + INTEGRATE3D( 2, 3 ); + INTEGRATE3D( 2, 4 ); + INTEGRATE3D( 2, 5 ); + INTEGRATE3D( 2, 6 ); + INTEGRATE3D( 2, 7 ); + INTEGRATE3D( 2, 8 ); + INTEGRATE3D( 2, 9 ); + + INTEGRATE3D( 3, 0 ); + INTEGRATE3D( 3, 1 ); + INTEGRATE3D( 3, 2 ); + INTEGRATE3D( 3, 3 ); + INTEGRATE3D( 3, 4 ); + INTEGRATE3D( 3, 5 ); + INTEGRATE3D( 3, 6 ); + INTEGRATE3D( 3, 7 ); + INTEGRATE3D( 3, 8 ); + INTEGRATE3D( 3, 9 ); + +// Undefine shape functions +#define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET #define UNDEFINE_P1_SHAPE_FUNCTIONS_TET #include "ShapeFunctionMacros.hpp" #undef UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET }; }; -// -------- + +// ======== // div_y -// -------- +// ======== template <> class P2ToP1Form_div< 1 > : public P2ToP1FormHyTeG { public: + + // --------- + // 2D Case + // --------- void integrateAll( const std::array< Point3D, 3 >& coords, Matrixr< 3, 6 >& elMat ) const final { @@ -260,20 +357,265 @@ class P2ToP1Form_div< 1 > : public P2ToP1FormHyTeG INTEGRATE2D( 2, 4 ); INTEGRATE2D( 2, 5 ); +#define UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE #include "ShapeFunctionMacros.hpp" #undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE - -#define UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE -#include "ShapeFunctionMacros.hpp" #undef UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #undef INTEGRATE2D }; + + // --------- + // 3D Case + // --------- void integrateAll( const std::array< Point3D, 4 >& coords, Matrixr< 4, 10 >& elMat ) const final { - WALBERLA_ABORT( "P2ToP1Form_div not implemented for 3D, yet." ); + // WALBERLA_ABORT( "P2ToP1Form_div<1> not implemented for 3D, yet." ); + +// Set shape functions and their derivatives +#define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define DEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef DEFINE_P1_SHAPE_FUNCTIONS_TET +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + +// Select cubature rule +#define CUBAPOINTS cubature::T4_points +#define CUBAWEIGHTS cubature::T4_weights + +// Executing quadrature rule +#define INTEGRATE3D( i, j ) \ + elMat( i, j ) = 0.0; \ + for ( uint_t k = 0; k < CUBAWEIGHTS.size(); k++ ) \ + { \ + real_t L2 = CUBAPOINTS[k][0]; \ + real_t L3 = CUBAPOINTS[k][1]; \ + real_t L4 = CUBAPOINTS[k][2]; \ + real_t L1 = 1.0 - L2 - L3 - L4; \ + \ + real_t dy = DX1N ## j*(tmp12 + tmp13 - tmp14 - tmp15 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + DX2N ## j*(-tmp12 + tmp15 + tmp16 - tmp17 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + DX3N ## j*(-tmp13 + tmp14 - tmp16 + tmp17 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]); \ + \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * dy * SF_M##i; \ + } + + // mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + // mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + // mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + + // Quadrature point mapped to computational triangle + Point3D mappedPt; + + // pre-compute values including Jacobian determinant of inverse pull-back mapping + real_t tmp0 = -coords[0][0]; + real_t tmp1 = tmp0 + coords[1][0]; + real_t tmp2 = -coords[0][1]; + real_t tmp3 = tmp2 + coords[2][1]; + real_t tmp4 = -coords[0][2]; + real_t tmp5 = tmp4 + coords[3][2]; + real_t tmp6 = tmp0 + coords[2][0]; + real_t tmp7 = tmp2 + coords[3][1]; + real_t tmp8 = tmp4 + coords[1][2]; + real_t tmp9 = tmp0 + coords[3][0]; + real_t tmp10 = tmp2 + coords[1][1]; + real_t tmp11 = tmp4 + coords[2][2]; + real_t tmp12 = coords[0][0]*coords[3][2]; + real_t tmp13 = coords[2][0]*coords[0][2]; + real_t tmp14 = coords[0][0]*coords[2][2]; + real_t tmp15 = coords[3][0]*coords[0][2]; + real_t tmp16 = coords[0][0]*coords[1][2]; + real_t tmp17 = coords[1][0]*coords[0][2]; + + real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; + real_t detPfac = std::abs( detP ) / detP; + + INTEGRATE3D( 0, 0 ); + INTEGRATE3D( 0, 1 ); + INTEGRATE3D( 0, 2 ); + INTEGRATE3D( 0, 3 ); + INTEGRATE3D( 0, 4 ); + INTEGRATE3D( 0, 5 ); + INTEGRATE3D( 0, 6 ); + INTEGRATE3D( 0, 7 ); + INTEGRATE3D( 0, 8 ); + INTEGRATE3D( 0, 9 ); + + INTEGRATE3D( 1, 0 ); + INTEGRATE3D( 1, 1 ); + INTEGRATE3D( 1, 2 ); + INTEGRATE3D( 1, 3 ); + INTEGRATE3D( 1, 4 ); + INTEGRATE3D( 1, 5 ); + INTEGRATE3D( 1, 6 ); + INTEGRATE3D( 1, 7 ); + INTEGRATE3D( 1, 8 ); + INTEGRATE3D( 1, 9 ); + + INTEGRATE3D( 2, 0 ); + INTEGRATE3D( 2, 1 ); + INTEGRATE3D( 2, 2 ); + INTEGRATE3D( 2, 3 ); + INTEGRATE3D( 2, 4 ); + INTEGRATE3D( 2, 5 ); + INTEGRATE3D( 2, 6 ); + INTEGRATE3D( 2, 7 ); + INTEGRATE3D( 2, 8 ); + INTEGRATE3D( 2, 9 ); + + INTEGRATE3D( 3, 0 ); + INTEGRATE3D( 3, 1 ); + INTEGRATE3D( 3, 2 ); + INTEGRATE3D( 3, 3 ); + INTEGRATE3D( 3, 4 ); + INTEGRATE3D( 3, 5 ); + INTEGRATE3D( 3, 6 ); + INTEGRATE3D( 3, 7 ); + INTEGRATE3D( 3, 8 ); + INTEGRATE3D( 3, 9 ); + +// Undefine shape functions +#define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + + }; +}; + + + +// ======== +// div_z +// ======== +template <> +class P2ToP1Form_div< 2 > : public P2ToP1FormHyTeG +{ + public: + + // --------- + // 2D Case + // --------- + void integrateAll( const std::array< Point3D, 3 >& coords, Matrixr< 3, 6 >& elMat ) const final + { + WALBERLA_ABORT( "P2ToP1Form_div< 2 >.integrateAll() not available for 3D!" ); + }; + + + // --------- + // 3D Case + // --------- + void integrateAll( const std::array< Point3D, 4 >& coords, Matrixr< 4, 10 >& elMat ) const final + { + +// Set shape functions and their derivatives +#define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define DEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef DEFINE_P1_SHAPE_FUNCTIONS_TET +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + +// Select cubature rule +#define CUBAPOINTS cubature::T4_points +#define CUBAWEIGHTS cubature::T4_weights + +// Executing quadrature rule +#define INTEGRATE3D( i, j ) \ + elMat( i, j ) = 0.0; \ + for ( uint_t k = 0; k < CUBAWEIGHTS.size(); k++ ) \ + { \ + real_t L2 = CUBAPOINTS[k][0]; \ + real_t L3 = CUBAPOINTS[k][1]; \ + real_t L4 = CUBAPOINTS[k][2]; \ + real_t L1 = 1.0 - L2 - L3 - L4; \ + \ + real_t dy = DX1N ## j*(tmp12 + tmp13 - tmp14 - tmp15 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + DX2N ## j*(-tmp13 + tmp14 + tmp16 - tmp17 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + DX3N ## j*(-tmp12 + tmp15 - tmp16 + tmp17 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]); \ + \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * dy * SF_M##i; \ + } + + // mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + // mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + // mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + + // Quadrature point mapped to computational triangle + Point3D mappedPt; + + // pre-compute values including Jacobian determinant of inverse pull-back mapping + real_t tmp0 = -coords[0][0]; + real_t tmp1 = tmp0 + coords[1][0]; + real_t tmp2 = -coords[0][1]; + real_t tmp3 = tmp2 + coords[2][1]; + real_t tmp4 = -coords[0][2]; + real_t tmp5 = tmp4 + coords[3][2]; + real_t tmp6 = tmp0 + coords[2][0]; + real_t tmp7 = tmp2 + coords[3][1]; + real_t tmp8 = tmp4 + coords[1][2]; + real_t tmp9 = tmp0 + coords[3][0]; + real_t tmp10 = tmp2 + coords[1][1]; + real_t tmp11 = tmp4 + coords[2][2]; + real_t tmp12 = coords[0][0]*coords[2][1]; + real_t tmp13 = coords[3][0]*coords[0][1]; + real_t tmp14 = coords[0][0]*coords[3][1]; + real_t tmp15 = coords[2][0]*coords[0][1]; + real_t tmp16 = coords[1][0]*coords[0][1]; + real_t tmp17 = coords[0][0]*coords[1][1]; + + real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; + real_t detPfac = std::abs( detP ) / detP; + + INTEGRATE3D( 0, 0 ); + INTEGRATE3D( 0, 1 ); + INTEGRATE3D( 0, 2 ); + INTEGRATE3D( 0, 3 ); + INTEGRATE3D( 0, 4 ); + INTEGRATE3D( 0, 5 ); + INTEGRATE3D( 0, 6 ); + INTEGRATE3D( 0, 7 ); + INTEGRATE3D( 0, 8 ); + INTEGRATE3D( 0, 9 ); + + INTEGRATE3D( 1, 0 ); + INTEGRATE3D( 1, 1 ); + INTEGRATE3D( 1, 2 ); + INTEGRATE3D( 1, 3 ); + INTEGRATE3D( 1, 4 ); + INTEGRATE3D( 1, 5 ); + INTEGRATE3D( 1, 6 ); + INTEGRATE3D( 1, 7 ); + INTEGRATE3D( 1, 8 ); + INTEGRATE3D( 1, 9 ); + + INTEGRATE3D( 2, 0 ); + INTEGRATE3D( 2, 1 ); + INTEGRATE3D( 2, 2 ); + INTEGRATE3D( 2, 3 ); + INTEGRATE3D( 2, 4 ); + INTEGRATE3D( 2, 5 ); + INTEGRATE3D( 2, 6 ); + INTEGRATE3D( 2, 7 ); + INTEGRATE3D( 2, 8 ); + INTEGRATE3D( 2, 9 ); + + INTEGRATE3D( 3, 0 ); + INTEGRATE3D( 3, 1 ); + INTEGRATE3D( 3, 2 ); + INTEGRATE3D( 3, 3 ); + INTEGRATE3D( 3, 4 ); + INTEGRATE3D( 3, 5 ); + INTEGRATE3D( 3, 6 ); + INTEGRATE3D( 3, 7 ); + INTEGRATE3D( 3, 8 ); + INTEGRATE3D( 3, 9 ); + +// Undefine shape functions +#define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + }; }; diff --git a/tests/hyteg/forms/HytegVsFenicsFormTest.cpp b/tests/hyteg/forms/HytegVsFenicsFormTest.cpp index 38d8016ca..42699faff 100644 --- a/tests/hyteg/forms/HytegVsFenicsFormTest.cpp +++ b/tests/hyteg/forms/HytegVsFenicsFormTest.cpp @@ -318,9 +318,9 @@ void run2DTestsWithAffineMap() void run3DTestsWithoutBlending() { // define our test tetrahedron - std::array< Point3D, 4 > theTet{ - Point3D( {0.0, 0.0, 0.0} ), Point3D( {1.0, 1.0, 0.0} ), Point3D( {-1.0, 0.5, 0.0} ), Point3D( {0.3, 0.21, -1.2} )}; - // std::array theTet{ Point3D({0.0, 0.0, 0.0}), Point3D({1.0, 0.0, 0.0}), Point3D({0.0, 1.0, 0.0}), Point3D({0.0, 0.0, 1.0}) }; + // std::array< Point3D, 4 > theTet{ + // Point3D( {0.0, 0.0, 0.0} ), Point3D( {1.0, 1.0, 0.0} ), Point3D( {-1.0, 0.5, 0.0} ), Point3D( {0.3, 0.21, -1.2} )}; + std::array theTet{ Point3D({0.0, 0.0, 0.0}), Point3D({1.0, 0.0, 0.0}), Point3D({0.0, 1.0, 0.0}), Point3D({0.0, 0.0, 1.0}) }; logSectionHeader( "P1 Mass Forms (3D)" ); compareForms< P1FenicsForm< fenics::NoAssemble, p1_tet_mass_cell_integral_0_otherwise >, P1Form_mass3D, Matrix4r, 3 >( theTet, @@ -333,6 +333,21 @@ void run3DTestsWithoutBlending() logSectionHeader( "P2 Laplace Forms (3D)" ); compareForms< P2FenicsForm< fenics::NoAssemble, p2_tet_diffusion_cell_integral_0_otherwise >, P2Form_laplace, Matrix10r, 3 >( theTet, 3e-14 ); + + logSectionHeader( "P2ToP1 DivX Forms (3D)" ); + compareForms< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_0_otherwise >, + P2ToP1Form_div< 0 >, + Matrixr< 4, 10 >, 3 >( theTet, 1.2e-14 ); + + logSectionHeader( "P2ToP1 DivY Forms (3D)" ); + compareForms< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_1_otherwise >, + P2ToP1Form_div< 1 >, + Matrixr< 4, 10 >, 3 >( theTet, 1.2e-14 ); + + logSectionHeader( "P2ToP1 DivZ Forms (3D)" ); + compareForms< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_2_otherwise >, + P2ToP1Form_div< 2 >, + Matrixr< 4, 10 >, 3 >( theTet, 1.2e-14 ); } -- GitLab From fddf2333f643233f3e994514f5459ee788e1c8dc Mon Sep 17 00:00:00 2001 From: Marcus Mohr Date: Mon, 27 Jul 2020 17:18:31 +0200 Subject: [PATCH 3/4] Adds blending to P2ToP1Form_div Comparison of local element matrix between HyTeG and FEniCS with "challenging" affine blending looks good. --- .../forms/form_hyteg_manual/P2ToP1FormDiv.hpp | 151 +++++++++++++----- tests/hyteg/forms/HytegVsFenicsFormTest.cpp | 144 ++++++++++------- 2 files changed, 200 insertions(+), 95 deletions(-) diff --git a/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp b/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp index 366110339..f6f6e0d82 100644 --- a/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp +++ b/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp @@ -171,15 +171,30 @@ class P2ToP1Form_div< 0 > : public P2ToP1FormHyTeG real_t L4 = CUBAPOINTS[k][2]; \ real_t L1 = 1.0 - L2 - L3 - L4; \ \ - real_t dy = DX1N ## j*(tmp12 + tmp13 - tmp14 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2]) + DX2N ## j*(-tmp13 + tmp14 + tmp15 - tmp16 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2]) + DX3N ## j*(-tmp12 - tmp15 + tmp16 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2]); \ + mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ \ - elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * dy * SF_M##i; \ + Matrix3r DPsi; \ + geometryMap_->evalDF( mappedPt, DPsi ); \ + \ + real_t tmp12 = DPsi(1,1)*DPsi(2,2); \ + real_t tmp13 = DPsi(1,2)*DPsi(2,0); \ + real_t tmp14 = DPsi(1,0)*DPsi(2,1); \ + real_t tmp15 = DPsi(1,2)*DPsi(2,1); \ + real_t tmp16 = DPsi(1,0)*DPsi(2,2); \ + real_t tmp17 = DPsi(1,1)*DPsi(2,0); \ + real_t tmp18 = tmp14 - tmp17; \ + real_t tmp23 = tmp13 - tmp16; \ + real_t tmp28 = tmp12 - tmp15; \ + \ + real_t diff = DX1N ## j*(tmp18*(tmp19 + tmp20 - tmp21 - tmp22 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + tmp23*(tmp24 + tmp25 - tmp26 - tmp27 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + tmp28*(tmp29 + tmp30 - tmp31 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2])) + DX2N ## j*(tmp18*(-tmp20 + tmp21 + tmp32 - tmp33 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + tmp23*(-tmp24 + tmp27 + tmp34 - tmp35 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + tmp28*(-tmp30 + tmp31 + tmp36 - tmp37 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2])) + DX3N ## j*(tmp18*(-tmp19 + tmp22 - tmp32 + tmp33 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]) + tmp23*(-tmp25 + tmp26 - tmp34 + tmp35 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]) + tmp28*(-tmp29 - tmp36 + tmp37 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2])); \ + \ + real_t detQ = DPsi(0,0)*tmp12 - DPsi(0,0)*tmp15 + DPsi(0,1)*tmp13 - DPsi(0,1)*tmp16 + DPsi(0,2)*tmp14 - DPsi(0,2)*tmp17; \ + real_t detQfac = std::abs( detQ ) / detQ; \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * detQfac * diff * SF_M##i; \ } - // mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ - // mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ - // mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ - // Quadrature point mapped to computational triangle Point3D mappedPt; @@ -196,13 +211,25 @@ class P2ToP1Form_div< 0 > : public P2ToP1FormHyTeG real_t tmp9 = tmp0 + coords[3][0]; real_t tmp10 = tmp2 + coords[1][1]; real_t tmp11 = tmp4 + coords[2][2]; - real_t tmp12 = coords[0][1]*coords[2][2]; - real_t tmp13 = coords[3][1]*coords[0][2]; - real_t tmp14 = coords[0][1]*coords[3][2]; - real_t tmp15 = coords[1][1]*coords[0][2]; - real_t tmp16 = coords[0][1]*coords[1][2]; + real_t tmp19 = coords[0][0]*coords[2][1]; + real_t tmp20 = coords[3][0]*coords[0][1]; + real_t tmp21 = coords[0][0]*coords[3][1]; + real_t tmp22 = coords[2][0]*coords[0][1]; + real_t tmp24 = coords[0][0]*coords[3][2]; + real_t tmp25 = coords[2][0]*coords[0][2]; + real_t tmp26 = coords[0][0]*coords[2][2]; + real_t tmp27 = coords[3][0]*coords[0][2]; + real_t tmp29 = coords[0][1]*coords[2][2]; + real_t tmp30 = coords[3][1]*coords[0][2]; + real_t tmp31 = coords[0][1]*coords[3][2]; + real_t tmp32 = coords[1][0]*coords[0][1]; + real_t tmp33 = coords[0][0]*coords[1][1]; + real_t tmp34 = coords[0][0]*coords[1][2]; + real_t tmp35 = coords[1][0]*coords[0][2]; + real_t tmp36 = coords[1][1]*coords[0][2]; + real_t tmp37 = coords[0][1]*coords[1][2]; - real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; + real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; real_t detPfac = std::abs( detP ) / detP; INTEGRATE3D( 0, 0 ); @@ -395,15 +422,30 @@ class P2ToP1Form_div< 1 > : public P2ToP1FormHyTeG real_t L4 = CUBAPOINTS[k][2]; \ real_t L1 = 1.0 - L2 - L3 - L4; \ \ - real_t dy = DX1N ## j*(tmp12 + tmp13 - tmp14 - tmp15 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + DX2N ## j*(-tmp12 + tmp15 + tmp16 - tmp17 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + DX3N ## j*(-tmp13 + tmp14 - tmp16 + tmp17 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]); \ + mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + \ + Matrix3r DPsi; \ + geometryMap_->evalDF( mappedPt, DPsi ); \ + \ + real_t tmp12 = DPsi(0,0)*DPsi(2,2); \ + real_t tmp13 = DPsi(0,1)*DPsi(2,0); \ + real_t tmp14 = DPsi(0,2)*DPsi(2,1); \ + real_t tmp15 = DPsi(0,0)*DPsi(2,1); \ + real_t tmp16 = DPsi(0,1)*DPsi(2,2); \ + real_t tmp17 = DPsi(0,2)*DPsi(2,0); \ + real_t tmp18 = tmp13 - tmp15; \ + real_t tmp23 = tmp12 - tmp17; \ + real_t tmp28 = tmp14 - tmp16; \ \ - elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * dy * SF_M##i; \ + real_t diff = DX1N ## j*(tmp18*(tmp19 + tmp20 - tmp21 - tmp22 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + tmp23*(tmp24 + tmp25 - tmp26 - tmp27 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + tmp28*(tmp29 + tmp30 - tmp31 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2])) + DX2N ## j*(tmp18*(-tmp20 + tmp21 + tmp32 - tmp33 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + tmp23*(-tmp24 + tmp27 + tmp34 - tmp35 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + tmp28*(-tmp30 + tmp31 + tmp36 - tmp37 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2])) + DX3N ## j*(tmp18*(-tmp19 + tmp22 - tmp32 + tmp33 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]) + tmp23*(-tmp25 + tmp26 - tmp34 + tmp35 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]) + tmp28*(-tmp29 - tmp36 + tmp37 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2])); \ + \ + real_t detQ = DPsi(1,0)*tmp14 - DPsi(1,0)*tmp16 + DPsi(1,1)*tmp12 - DPsi(1,1)*tmp17 + DPsi(1,2)*tmp13 - DPsi(1,2)*tmp15; \ + real_t detQfac = std::abs( detQ ) / detQ; \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * detQfac * diff * SF_M##i; \ } - // mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ - // mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ - // mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ - // Quadrature point mapped to computational triangle Point3D mappedPt; @@ -420,12 +462,23 @@ class P2ToP1Form_div< 1 > : public P2ToP1FormHyTeG real_t tmp9 = tmp0 + coords[3][0]; real_t tmp10 = tmp2 + coords[1][1]; real_t tmp11 = tmp4 + coords[2][2]; - real_t tmp12 = coords[0][0]*coords[3][2]; - real_t tmp13 = coords[2][0]*coords[0][2]; - real_t tmp14 = coords[0][0]*coords[2][2]; - real_t tmp15 = coords[3][0]*coords[0][2]; - real_t tmp16 = coords[0][0]*coords[1][2]; - real_t tmp17 = coords[1][0]*coords[0][2]; + real_t tmp19 = coords[0][0]*coords[2][1]; + real_t tmp20 = coords[3][0]*coords[0][1]; + real_t tmp21 = coords[0][0]*coords[3][1]; + real_t tmp22 = coords[2][0]*coords[0][1]; + real_t tmp24 = coords[0][0]*coords[3][2]; + real_t tmp25 = coords[2][0]*coords[0][2]; + real_t tmp26 = coords[0][0]*coords[2][2]; + real_t tmp27 = coords[3][0]*coords[0][2]; + real_t tmp29 = coords[0][1]*coords[2][2]; + real_t tmp30 = coords[3][1]*coords[0][2]; + real_t tmp31 = coords[0][1]*coords[3][2]; + real_t tmp32 = coords[1][0]*coords[0][1]; + real_t tmp33 = coords[0][0]*coords[1][1]; + real_t tmp34 = coords[0][0]*coords[1][2]; + real_t tmp35 = coords[1][0]*coords[0][2]; + real_t tmp36 = coords[1][1]*coords[0][2]; + real_t tmp37 = coords[0][1]*coords[1][2]; real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; real_t detPfac = std::abs( detP ) / detP; @@ -530,15 +583,30 @@ class P2ToP1Form_div< 2 > : public P2ToP1FormHyTeG real_t L4 = CUBAPOINTS[k][2]; \ real_t L1 = 1.0 - L2 - L3 - L4; \ \ - real_t dy = DX1N ## j*(tmp12 + tmp13 - tmp14 - tmp15 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + DX2N ## j*(-tmp13 + tmp14 + tmp16 - tmp17 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + DX3N ## j*(-tmp12 + tmp15 - tmp16 + tmp17 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]); \ + mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + \ + Matrix3r DPsi; \ + geometryMap_->evalDF( mappedPt, DPsi ); \ \ - elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * dy * SF_M##i; \ + real_t tmp12 = DPsi(0,0)*DPsi(1,1); \ + real_t tmp13 = DPsi(0,1)*DPsi(1,2); \ + real_t tmp14 = DPsi(0,2)*DPsi(1,0); \ + real_t tmp15 = DPsi(0,0)*DPsi(1,2); \ + real_t tmp16 = DPsi(0,1)*DPsi(1,0); \ + real_t tmp17 = DPsi(0,2)*DPsi(1,1); \ + real_t tmp18 = tmp12 - tmp16; \ + real_t tmp23 = tmp14 - tmp15; \ + real_t tmp28 = tmp13 - tmp17; \ + \ + real_t diff = DX1N ## j*(tmp18*(tmp19 + tmp20 - tmp21 - tmp22 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + tmp23*(tmp24 + tmp25 - tmp26 - tmp27 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + tmp28*(tmp29 + tmp30 - tmp31 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2])) + DX2N ## j*(tmp18*(-tmp20 + tmp21 + tmp32 - tmp33 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + tmp23*(-tmp24 + tmp27 + tmp34 - tmp35 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + tmp28*(-tmp30 + tmp31 + tmp36 - tmp37 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2])) + DX3N ## j*(tmp18*(-tmp19 + tmp22 - tmp32 + tmp33 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]) + tmp23*(-tmp25 + tmp26 - tmp34 + tmp35 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]) + tmp28*(-tmp29 - tmp36 + tmp37 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2])); \ + \ + real_t detQ = DPsi(2,0)*tmp13 - DPsi(2,0)*tmp17 + DPsi(2,1)*tmp14 - DPsi(2,1)*tmp15 + DPsi(2,2)*tmp12 - DPsi(2,2)*tmp16; \ + real_t detQfac = std::abs( detQ ) / detQ; \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * detQfac * diff * SF_M##i; \ } - // mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ - // mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ - // mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ - // Quadrature point mapped to computational triangle Point3D mappedPt; @@ -555,12 +623,23 @@ class P2ToP1Form_div< 2 > : public P2ToP1FormHyTeG real_t tmp9 = tmp0 + coords[3][0]; real_t tmp10 = tmp2 + coords[1][1]; real_t tmp11 = tmp4 + coords[2][2]; - real_t tmp12 = coords[0][0]*coords[2][1]; - real_t tmp13 = coords[3][0]*coords[0][1]; - real_t tmp14 = coords[0][0]*coords[3][1]; - real_t tmp15 = coords[2][0]*coords[0][1]; - real_t tmp16 = coords[1][0]*coords[0][1]; - real_t tmp17 = coords[0][0]*coords[1][1]; + real_t tmp19 = coords[0][0]*coords[2][1]; + real_t tmp20 = coords[3][0]*coords[0][1]; + real_t tmp21 = coords[0][0]*coords[3][1]; + real_t tmp22 = coords[2][0]*coords[0][1]; + real_t tmp24 = coords[0][0]*coords[3][2]; + real_t tmp25 = coords[2][0]*coords[0][2]; + real_t tmp26 = coords[0][0]*coords[2][2]; + real_t tmp27 = coords[3][0]*coords[0][2]; + real_t tmp29 = coords[0][1]*coords[2][2]; + real_t tmp30 = coords[3][1]*coords[0][2]; + real_t tmp31 = coords[0][1]*coords[3][2]; + real_t tmp32 = coords[1][0]*coords[0][1]; + real_t tmp33 = coords[0][0]*coords[1][1]; + real_t tmp34 = coords[0][0]*coords[1][2]; + real_t tmp35 = coords[1][0]*coords[0][2]; + real_t tmp36 = coords[1][1]*coords[0][2]; + real_t tmp37 = coords[0][1]*coords[1][2]; real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; real_t detPfac = std::abs( detP ) / detP; diff --git a/tests/hyteg/forms/HytegVsFenicsFormTest.cpp b/tests/hyteg/forms/HytegVsFenicsFormTest.cpp index 42699faff..17df6d3ec 100644 --- a/tests/hyteg/forms/HytegVsFenicsFormTest.cpp +++ b/tests/hyteg/forms/HytegVsFenicsFormTest.cpp @@ -78,9 +78,9 @@ real_t template < class FormFEniCS, class FormHyTeG, typename matType, uint_t dim > void compareForms( const std::array< Point3D, dim + 1 >& element, real_t tol, - std::shared_ptr< GeometryMap > map = nullptr, - real_t scaleFactor = real_c( 1 ), - bool applyMapToElement = false ) + std::shared_ptr< GeometryMap > map = nullptr, + real_t scaleFactor = real_c( 1 ), + bool applyMapToElement = false ) { bool wBlending = true; @@ -98,10 +98,12 @@ void compareForms( const std::array< Point3D, dim + 1 >& element, // apply blending to element fed to FEniCS (useful for the tests with // an affine mapping) std::array< Point3D, dim + 1 > elementForFenics = element; - if( applyMapToElement ) { - for( uint_t k = 0; k <= dim; ++k ) { - map->evalF( element[k], elementForFenics[k] ); - } + if ( applyMapToElement ) + { + for ( uint_t k = 0; k <= dim; ++k ) + { + map->evalF( element[k], elementForFenics[k] ); + } } // assemble element matrices @@ -124,7 +126,6 @@ void compareForms( const std::array< Point3D, dim + 1 >& element, WALBERLA_CHECK_LESS_EQUAL( error, tol ); } - template < class FormFEniCS, class FormHyTeG, typename matType, uint_t dim > void compareScaled( const std::array< Point3D, dim + 1 >& element, real_t tol, @@ -134,16 +135,12 @@ void compareScaled( const std::array< Point3D, dim + 1 >& element, compareForms< FormFEniCS, FormHyTeG, matType, dim >( element, tol, map, scaleFactor, false ); } - template < class FormFEniCS, class FormHyTeG, typename matType, uint_t dim > -void compareUsingAffineMap( const std::array< Point3D, dim + 1 >& element, - real_t tol, - std::shared_ptr< GeometryMap > map ) +void compareUsingAffineMap( const std::array< Point3D, dim + 1 >& element, real_t tol, std::shared_ptr< GeometryMap > map ) { compareForms< FormFEniCS, FormHyTeG, matType, dim >( element, tol, map, 0.0, true ); } - void run2DTestsWithoutBlending() { // define our test triangle @@ -231,11 +228,11 @@ void run2DTestsWithAffineMap() { // define our affine map (rotation + scaling + shift) Matrix2r mat; - real_t phi = 2.0 / 9.0 * pi; - mat(0,0) = +std::cos( phi ); - mat(0,1) = -std::sin( phi ); - mat(1,0) = +std::sin( phi ) * 2.25; - mat(1,1) = +std::cos( phi ) * 2.25; + real_t phi = 2.0 / 9.0 * pi; + mat( 0, 0 ) = +std::cos( phi ); + mat( 0, 1 ) = -std::sin( phi ); + mat( 1, 0 ) = +std::sin( phi ) * 2.25; + mat( 1, 1 ) = +std::cos( phi ) * 2.25; Point2D vec( {-7.0, 3.0} ); auto map = std::make_shared< AffineMap2D >( mat, vec ); @@ -253,15 +250,15 @@ void run2DTestsWithAffineMap() logSectionHeader( "P2ToP1 DivX Forms" ); compareUsingAffineMap< P2ToP1FenicsForm< p2_to_p1_div_cell_integral_0_otherwise, fenics::NoAssemble >, - P2ToP1Form_div< 0 >, - Matrixr< 3, 6 >, - 2 >( triangle, 3.5e-14, map ); + P2ToP1Form_div< 0 >, + Matrixr< 3, 6 >, + 2 >( triangle, 3.5e-14, map ); logSectionHeader( "P2ToP1 DivY Forms" ); compareUsingAffineMap< P2ToP1FenicsForm< p2_to_p1_div_cell_integral_1_otherwise, fenics::NoAssemble >, - P2ToP1Form_div< 1 >, - Matrixr< 3, 6 >, - 2 >( triangle, 5e-14, map ); + P2ToP1Form_div< 1 >, + Matrixr< 3, 6 >, + 2 >( triangle, 5e-14, map ); logSectionHeader( "P1 DivX Forms" ); compareUsingAffineMap< P1FenicsForm< p1_div_cell_integral_0_otherwise, fenics::NoAssemble >, P1Form_div_1, Matrix3r, 2 >( @@ -297,30 +294,35 @@ void run2DTestsWithAffineMap() logSectionHeader( "P1ToP2 DivX^T Forms" ); compareUsingAffineMap< P1ToP2FenicsForm< p1_to_p2_divt_cell_integral_0_otherwise, fenics::NoAssemble >, - P1ToP2Form_divt< 0 >, - Matrixr< 6, 3 >, - 2 >( triangle, 5e-14, map ); + P1ToP2Form_divt< 0 >, + Matrixr< 6, 3 >, + 2 >( triangle, 5e-14, map ); logSectionHeader( "P1ToP2 DivY^T Forms" ); compareUsingAffineMap< P1ToP2FenicsForm< p1_to_p2_divt_cell_integral_1_otherwise, fenics::NoAssemble >, - P1ToP2Form_divt< 1 >, - Matrixr< 6, 3 >, - 2 >( triangle, 5e-14, map ); + P1ToP2Form_divt< 1 >, + Matrixr< 6, 3 >, + 2 >( triangle, 5e-14, map ); logSectionHeader( "P1 Diffusion Forms" ); - compareUsingAffineMap< P1FenicsForm< p1_diffusion_cell_integral_0_otherwise, fenics::NoAssemble >, P1Form_laplace, Matrix3r, 2 >( triangle, 3e-15, map ); + compareUsingAffineMap< P1FenicsForm< p1_diffusion_cell_integral_0_otherwise, fenics::NoAssemble >, + P1Form_laplace, + Matrix3r, + 2 >( triangle, 3e-15, map ); logSectionHeader( "P2 Laplace Form" ); - compareUsingAffineMap< P2FenicsForm< p2_diffusion_cell_integral_0_otherwise, fenics::NoAssemble >, P2Form_laplace, Matrix6r, 2 >( triangle, 1e-13, map ); + compareUsingAffineMap< P2FenicsForm< p2_diffusion_cell_integral_0_otherwise, fenics::NoAssemble >, + P2Form_laplace, + Matrix6r, + 2 >( triangle, 1e-13, map ); } - void run3DTestsWithoutBlending() { // define our test tetrahedron - // std::array< Point3D, 4 > theTet{ - // Point3D( {0.0, 0.0, 0.0} ), Point3D( {1.0, 1.0, 0.0} ), Point3D( {-1.0, 0.5, 0.0} ), Point3D( {0.3, 0.21, -1.2} )}; - std::array theTet{ Point3D({0.0, 0.0, 0.0}), Point3D({1.0, 0.0, 0.0}), Point3D({0.0, 1.0, 0.0}), Point3D({0.0, 0.0, 1.0}) }; + std::array< Point3D, 4 > theTet{ + Point3D( {0.0, 0.0, 0.0} ), Point3D( {1.0, 1.0, 0.0} ), Point3D( {-1.0, 0.5, 0.0} ), Point3D( {0.3, 0.21, -1.2} )}; + // std::array theTet{ Point3D({0.0, 0.0, 0.0}), Point3D({1.0, 0.0, 0.0}), Point3D({0.0, 1.0, 0.0}), Point3D({0.0, 0.0, 1.0}) }; logSectionHeader( "P1 Mass Forms (3D)" ); compareForms< P1FenicsForm< fenics::NoAssemble, p1_tet_mass_cell_integral_0_otherwise >, P1Form_mass3D, Matrix4r, 3 >( theTet, @@ -337,20 +339,22 @@ void run3DTestsWithoutBlending() logSectionHeader( "P2ToP1 DivX Forms (3D)" ); compareForms< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_0_otherwise >, P2ToP1Form_div< 0 >, - Matrixr< 4, 10 >, 3 >( theTet, 1.2e-14 ); + Matrixr< 4, 10 >, + 3 >( theTet, 1e-14 ); logSectionHeader( "P2ToP1 DivY Forms (3D)" ); compareForms< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_1_otherwise >, P2ToP1Form_div< 1 >, - Matrixr< 4, 10 >, 3 >( theTet, 1.2e-14 ); + Matrixr< 4, 10 >, + 3 >( theTet, 1e-14 ); logSectionHeader( "P2ToP1 DivZ Forms (3D)" ); compareForms< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_2_otherwise >, P2ToP1Form_div< 2 >, - Matrixr< 4, 10 >, 3 >( theTet, 1.2e-14 ); + Matrixr< 4, 10 >, + 3 >( theTet, 1e-14 ); } - void run3DTestsWithAffineMap() { // define our affine map @@ -359,20 +363,20 @@ void run3DTestsWithAffineMap() // simple scaling mat( 0, 0 ) = real_c( 2.0 ); mat( 1, 1 ) = real_c( 1.0 ); - mat( 2, 2 ) = real_c( std::exp(1.0) ); + mat( 2, 2 ) = real_c( std::exp( 1.0 ) ); // more interesting variant (rotation around x-axis, y-axis and scaling of z-component) #define CHALLENGING #ifdef CHALLENGING - mat(0,0) = +8.660254037844387e-01; - mat(0,1) = -1.545084971874737e-01; - mat(0,2) = -9.510565162951534e-01; - mat(1,0) = +0.000000000000000e+00; - mat(1,1) = +9.510565162951535e-01; - mat(1,2) = -6.180339887498948e-01; - mat(2,0) = +4.999999999999999e-01; - mat(2,1) = +2.676165673298174e-01; - mat(2,2) = +1.647278207092664e+00; + mat( 0, 0 ) = +8.660254037844387e-01; + mat( 0, 1 ) = -1.545084971874737e-01; + mat( 0, 2 ) = -9.510565162951534e-01; + mat( 1, 0 ) = +0.000000000000000e+00; + mat( 1, 1 ) = +9.510565162951535e-01; + mat( 1, 2 ) = -6.180339887498948e-01; + mat( 2, 0 ) = +4.999999999999999e-01; + mat( 2, 1 ) = +2.676165673298174e-01; + mat( 2, 2 ) = +1.647278207092664e+00; #endif #undef CHALLENGING @@ -380,22 +384,44 @@ void run3DTestsWithAffineMap() auto map = std::make_shared< AffineMap3D >( mat, vec ); // define our test tetrahedrons - Point3D v1( { 0.0, 0.00, 0.0} ); - Point3D v2( { 1.0, 1.00, 0.0} ); - Point3D v3( {-1.0, 0.50, 0.0} ); - Point3D v4( { 0.3, 0.21, -1.2} ); - std::array< Point3D, 4 > theTet{ v1, v2, v3, v4 }; + Point3D v1( {0.0, 0.00, 0.0} ); + Point3D v2( {1.0, 1.00, 0.0} ); + Point3D v3( {-1.0, 0.50, 0.0} ); + Point3D v4( {0.3, 0.21, -1.2} ); + std::array< Point3D, 4 > theTet{v1, v2, v3, v4}; logSectionHeader( "P1 Mass Forms (3D)" ); - compareUsingAffineMap< P1FenicsForm< fenics::NoAssemble, p1_tet_mass_cell_integral_0_otherwise >, P1Form_mass3D, Matrix4r, 3 >( theTet, 5e-15, map ); + compareUsingAffineMap< P1FenicsForm< fenics::NoAssemble, p1_tet_mass_cell_integral_0_otherwise >, P1Form_mass3D, Matrix4r, 3 >( + theTet, 5e-15, map ); logSectionHeader( "P2 Mass Forms (3D)" ); - compareUsingAffineMap< P2FenicsForm< fenics::NoAssemble, p2_tet_mass_cell_integral_0_otherwise >, P2Form_mass, Matrix10r, 3 >( theTet, 5e-7, map ); // need to improve our cubature !!! + compareUsingAffineMap< P2FenicsForm< fenics::NoAssemble, p2_tet_mass_cell_integral_0_otherwise >, P2Form_mass, Matrix10r, 3 >( + theTet, 5e-7, map ); // need to improve our cubature !!! logSectionHeader( "P2 Laplace Forms (3D)" ); - compareUsingAffineMap< P2FenicsForm< fenics::NoAssemble, p2_tet_diffusion_cell_integral_0_otherwise >, P2Form_laplace, Matrix10r, 3 >( theTet, 5e-14, map ); -} + compareUsingAffineMap< P2FenicsForm< fenics::NoAssemble, p2_tet_diffusion_cell_integral_0_otherwise >, + P2Form_laplace, + Matrix10r, + 3 >( theTet, 5e-14, map ); + logSectionHeader( "P2ToP1 DivX Forms (3D)" ); + compareUsingAffineMap< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_0_otherwise >, + P2ToP1Form_div< 0 >, + Matrixr< 4, 10 >, + 3 >( theTet, 1e-14, map ); + + logSectionHeader( "P2ToP1 DivY Forms (3D)" ); + compareUsingAffineMap< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_1_otherwise >, + P2ToP1Form_div< 1 >, + Matrixr< 4, 10 >, + 3 >( theTet, 1e-14, map ); + + logSectionHeader( "P2ToP1 DivZ Forms (3D)" ); + compareUsingAffineMap< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_2_otherwise >, + P2ToP1Form_div< 2 >, + Matrixr< 4, 10 >, + 3 >( theTet, 1e-14, map ); +} int main( int argc, char** argv ) { -- GitLab From 3d86745186a9ba49ded1426ec2bbbf49e11ac882 Mon Sep 17 00:00:00 2001 From: Marcus Mohr Date: Mon, 27 Jul 2020 17:46:13 +0200 Subject: [PATCH 4/4] Implements P1ToP2FormDivT with blending. --- .../form_hyteg_manual/P1ToP2FormDivT.hpp | 492 +++++++++++++++++- .../forms/form_hyteg_manual/P2ToP1FormDiv.hpp | 1 - tests/hyteg/forms/HytegVsFenicsFormTest.cpp | 18 + 3 files changed, 488 insertions(+), 23 deletions(-) diff --git a/src/hyteg/forms/form_hyteg_manual/P1ToP2FormDivT.hpp b/src/hyteg/forms/form_hyteg_manual/P1ToP2FormDivT.hpp index 3cf2cc76c..6cd530242 100644 --- a/src/hyteg/forms/form_hyteg_manual/P1ToP2FormDivT.hpp +++ b/src/hyteg/forms/form_hyteg_manual/P1ToP2FormDivT.hpp @@ -39,25 +39,26 @@ class P1ToP2Form_divt : public P1ToP2FormHyTeG }; }; -// --------- + +// ========= // div_x^T -// --------- +// ========= template <> class P1ToP2Form_divt< 0 > : public P1ToP2FormHyTeG { public: + + // --------- + // 2D Case + // --------- void integrateAll( const std::array< Point3D, 3 >& coords, Matrixr< 6, 3 >& elMat ) const final { -// Derivatives of P2 shape functions on unit triangle #define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE -#include "ShapeFunctionMacros.hpp" -#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE - -// P1 shape functions on unit triangle #define DEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #include "ShapeFunctionMacros.hpp" #undef DEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE // Select quadrature rule #define QUADPOINTS quadrature::D5_points @@ -127,41 +128,183 @@ class P1ToP2Form_divt< 0 > : public P1ToP2FormHyTeG INTEGRATE2D( 5, 2 ); #define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE -#include "ShapeFunctionMacros.hpp" -#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE - #define UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #include "ShapeFunctionMacros.hpp" #undef UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE #undef INTEGRATE2D }; + + // --------- + // 3D Case + // --------- void integrateAll( const std::array< Point3D, 4 >& coords, Matrixr< 10, 4 >& elMat ) const final { - WALBERLA_ABORT( "P1ToP2Form_divt not implemented for 3D, yet." ); + +// Set shape functions and their derivatives +#define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define DEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef DEFINE_P1_SHAPE_FUNCTIONS_TET +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + +// Select cubature rule +#define CUBAPOINTS cubature::T4_points +#define CUBAWEIGHTS cubature::T4_weights + +// Executing quadrature rule +#define INTEGRATE3D( i, j ) \ + elMat( i, j ) = 0.0; \ + for ( uint_t k = 0; k < CUBAWEIGHTS.size(); k++ ) \ + { \ + real_t L2 = CUBAPOINTS[k][0]; \ + real_t L3 = CUBAPOINTS[k][1]; \ + real_t L4 = CUBAPOINTS[k][2]; \ + real_t L1 = 1.0 - L2 - L3 - L4; \ + \ + mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + \ + Matrix3r DPsi; \ + geometryMap_->evalDF( mappedPt, DPsi ); \ + \ + real_t tmp12 = DPsi(1,1)*DPsi(2,2); \ + real_t tmp13 = DPsi(1,2)*DPsi(2,0); \ + real_t tmp14 = DPsi(1,0)*DPsi(2,1); \ + real_t tmp15 = DPsi(1,2)*DPsi(2,1); \ + real_t tmp16 = DPsi(1,0)*DPsi(2,2); \ + real_t tmp17 = DPsi(1,1)*DPsi(2,0); \ + real_t tmp18 = tmp14 - tmp17; \ + real_t tmp23 = tmp13 - tmp16; \ + real_t tmp28 = tmp12 - tmp15; \ + \ + real_t diff = DX1N ## i*(tmp18*(tmp19 + tmp20 - tmp21 - tmp22 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + tmp23*(tmp24 + tmp25 - tmp26 - tmp27 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + tmp28*(tmp29 + tmp30 - tmp31 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2])) + DX2N ## i*(tmp18*(-tmp20 + tmp21 + tmp32 - tmp33 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + tmp23*(-tmp24 + tmp27 + tmp34 - tmp35 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + tmp28*(-tmp30 + tmp31 + tmp36 - tmp37 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2])) + DX3N ## i*(tmp18*(-tmp19 + tmp22 - tmp32 + tmp33 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]) + tmp23*(-tmp25 + tmp26 - tmp34 + tmp35 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]) + tmp28*(-tmp29 - tmp36 + tmp37 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2])); \ + \ + real_t detQ = DPsi(0,0)*tmp12 - DPsi(0,0)*tmp15 + DPsi(0,1)*tmp13 - DPsi(0,1)*tmp16 + DPsi(0,2)*tmp14 - DPsi(0,2)*tmp17; \ + real_t detQfac = std::abs( detQ ) / detQ; \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * detQfac * diff * SF_M##j; \ + } + + // Quadrature point mapped to computational triangle + Point3D mappedPt; + + // pre-compute values including Jacobian determinant of inverse pull-back mapping + real_t tmp0 = -coords[0][0]; + real_t tmp1 = tmp0 + coords[1][0]; + real_t tmp2 = -coords[0][1]; + real_t tmp3 = tmp2 + coords[2][1]; + real_t tmp4 = -coords[0][2]; + real_t tmp5 = tmp4 + coords[3][2]; + real_t tmp6 = tmp0 + coords[2][0]; + real_t tmp7 = tmp2 + coords[3][1]; + real_t tmp8 = tmp4 + coords[1][2]; + real_t tmp9 = tmp0 + coords[3][0]; + real_t tmp10 = tmp2 + coords[1][1]; + real_t tmp11 = tmp4 + coords[2][2]; + real_t tmp19 = coords[0][0]*coords[2][1]; + real_t tmp20 = coords[3][0]*coords[0][1]; + real_t tmp21 = coords[0][0]*coords[3][1]; + real_t tmp22 = coords[2][0]*coords[0][1]; + real_t tmp24 = coords[0][0]*coords[3][2]; + real_t tmp25 = coords[2][0]*coords[0][2]; + real_t tmp26 = coords[0][0]*coords[2][2]; + real_t tmp27 = coords[3][0]*coords[0][2]; + real_t tmp29 = coords[0][1]*coords[2][2]; + real_t tmp30 = coords[3][1]*coords[0][2]; + real_t tmp31 = coords[0][1]*coords[3][2]; + real_t tmp32 = coords[1][0]*coords[0][1]; + real_t tmp33 = coords[0][0]*coords[1][1]; + real_t tmp34 = coords[0][0]*coords[1][2]; + real_t tmp35 = coords[1][0]*coords[0][2]; + real_t tmp36 = coords[1][1]*coords[0][2]; + real_t tmp37 = coords[0][1]*coords[1][2]; + + real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; + real_t detPfac = std::abs( detP ) / detP; + + INTEGRATE3D( 0, 0 ); + INTEGRATE3D( 0, 1 ); + INTEGRATE3D( 0, 2 ); + INTEGRATE3D( 0, 3 ); + + INTEGRATE3D( 1, 0 ); + INTEGRATE3D( 1, 1 ); + INTEGRATE3D( 1, 2 ); + INTEGRATE3D( 1, 3 ); + + INTEGRATE3D( 2, 0 ); + INTEGRATE3D( 2, 1 ); + INTEGRATE3D( 2, 2 ); + INTEGRATE3D( 2, 3 ); + + INTEGRATE3D( 3, 0 ); + INTEGRATE3D( 3, 1 ); + INTEGRATE3D( 3, 2 ); + INTEGRATE3D( 3, 3 ); + + INTEGRATE3D( 4, 0 ); + INTEGRATE3D( 4, 1 ); + INTEGRATE3D( 4, 2 ); + INTEGRATE3D( 4, 3 ); + + INTEGRATE3D( 5, 0 ); + INTEGRATE3D( 5, 1 ); + INTEGRATE3D( 5, 2 ); + INTEGRATE3D( 5, 3 ); + + INTEGRATE3D( 6, 0 ); + INTEGRATE3D( 6, 1 ); + INTEGRATE3D( 6, 2 ); + INTEGRATE3D( 6, 3 ); + + INTEGRATE3D( 7, 0 ); + INTEGRATE3D( 7, 1 ); + INTEGRATE3D( 7, 2 ); + INTEGRATE3D( 7, 3 ); + + INTEGRATE3D( 8, 0 ); + INTEGRATE3D( 8, 1 ); + INTEGRATE3D( 8, 2 ); + INTEGRATE3D( 8, 3 ); + + INTEGRATE3D( 9, 0 ); + INTEGRATE3D( 9, 1 ); + INTEGRATE3D( 9, 2 ); + INTEGRATE3D( 9, 3 ); + +// Undefine shape functions +#define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + }; }; -// --------- + +// ========= // div_y^T -// --------- +// ========= template <> class P1ToP2Form_divt< 1 > : public P1ToP2FormHyTeG { public: + + // --------- + // 2D Case + // --------- void integrateAll( const std::array< Point3D, 3 >& coords, Matrixr< 6, 3 >& elMat ) const final { -// Derivatives of P2 shape functions on unit triangle #define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE -#include "ShapeFunctionMacros.hpp" -#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE - -// P1 shape functions on unit triangle #define DEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #include "ShapeFunctionMacros.hpp" #undef DEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE // Select quadrature rule #define QUADPOINTS quadrature::D5_points @@ -231,19 +374,324 @@ class P1ToP2Form_divt< 1 > : public P1ToP2FormHyTeG INTEGRATE2D( 5, 2 ); #define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE -#include "ShapeFunctionMacros.hpp" -#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE - #define UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE #include "ShapeFunctionMacros.hpp" #undef UNDEFINE_P1_SHAPE_FUNCTIONS_TRIANGLE +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TRIANGLE #undef INTEGRATE2D }; + + // --------- + // 3D Case + // --------- + void integrateAll( const std::array< Point3D, 4 >& coords, Matrixr< 10, 4 >& elMat ) const final + { + +// Set shape functions and their derivatives +#define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define DEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef DEFINE_P1_SHAPE_FUNCTIONS_TET +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + +// Select cubature rule +#define CUBAPOINTS cubature::T4_points +#define CUBAWEIGHTS cubature::T4_weights + +// Executing quadrature rule +#define INTEGRATE3D( i, j ) \ + elMat( i, j ) = 0.0; \ + for ( uint_t k = 0; k < CUBAWEIGHTS.size(); k++ ) \ + { \ + real_t L2 = CUBAPOINTS[k][0]; \ + real_t L3 = CUBAPOINTS[k][1]; \ + real_t L4 = CUBAPOINTS[k][2]; \ + real_t L1 = 1.0 - L2 - L3 - L4; \ + \ + mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + \ + Matrix3r DPsi; \ + geometryMap_->evalDF( mappedPt, DPsi ); \ + \ + real_t tmp12 = DPsi(0,0)*DPsi(2,2); \ + real_t tmp13 = DPsi(0,1)*DPsi(2,0); \ + real_t tmp14 = DPsi(0,2)*DPsi(2,1); \ + real_t tmp15 = DPsi(0,0)*DPsi(2,1); \ + real_t tmp16 = DPsi(0,1)*DPsi(2,2); \ + real_t tmp17 = DPsi(0,2)*DPsi(2,0); \ + real_t tmp18 = tmp13 - tmp15; \ + real_t tmp23 = tmp12 - tmp17; \ + real_t tmp28 = tmp14 - tmp16; \ + \ + real_t diff = DX1N ## i*(tmp18*(tmp19 + tmp20 - tmp21 - tmp22 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + tmp23*(tmp24 + tmp25 - tmp26 - tmp27 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + tmp28*(tmp29 + tmp30 - tmp31 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2])) + DX2N ## i*(tmp18*(-tmp20 + tmp21 + tmp32 - tmp33 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + tmp23*(-tmp24 + tmp27 + tmp34 - tmp35 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + tmp28*(-tmp30 + tmp31 + tmp36 - tmp37 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2])) + DX3N ## i*(tmp18*(-tmp19 + tmp22 - tmp32 + tmp33 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]) + tmp23*(-tmp25 + tmp26 - tmp34 + tmp35 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]) + tmp28*(-tmp29 - tmp36 + tmp37 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2])); \ + \ + real_t detQ = DPsi(1,0)*tmp14 - DPsi(1,0)*tmp16 + DPsi(1,1)*tmp12 - DPsi(1,1)*tmp17 + DPsi(1,2)*tmp13 - DPsi(1,2)*tmp15; \ + real_t detQfac = std::abs( detQ ) / detQ; \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * detQfac * diff * SF_M##j; \ + } + + // Quadrature point mapped to computational triangle + Point3D mappedPt; + + // pre-compute values including Jacobian determinant of inverse pull-back mapping + real_t tmp0 = -coords[0][0]; + real_t tmp1 = tmp0 + coords[1][0]; + real_t tmp2 = -coords[0][1]; + real_t tmp3 = tmp2 + coords[2][1]; + real_t tmp4 = -coords[0][2]; + real_t tmp5 = tmp4 + coords[3][2]; + real_t tmp6 = tmp0 + coords[2][0]; + real_t tmp7 = tmp2 + coords[3][1]; + real_t tmp8 = tmp4 + coords[1][2]; + real_t tmp9 = tmp0 + coords[3][0]; + real_t tmp10 = tmp2 + coords[1][1]; + real_t tmp11 = tmp4 + coords[2][2]; + real_t tmp19 = coords[0][0]*coords[2][1]; + real_t tmp20 = coords[3][0]*coords[0][1]; + real_t tmp21 = coords[0][0]*coords[3][1]; + real_t tmp22 = coords[2][0]*coords[0][1]; + real_t tmp24 = coords[0][0]*coords[3][2]; + real_t tmp25 = coords[2][0]*coords[0][2]; + real_t tmp26 = coords[0][0]*coords[2][2]; + real_t tmp27 = coords[3][0]*coords[0][2]; + real_t tmp29 = coords[0][1]*coords[2][2]; + real_t tmp30 = coords[3][1]*coords[0][2]; + real_t tmp31 = coords[0][1]*coords[3][2]; + real_t tmp32 = coords[1][0]*coords[0][1]; + real_t tmp33 = coords[0][0]*coords[1][1]; + real_t tmp34 = coords[0][0]*coords[1][2]; + real_t tmp35 = coords[1][0]*coords[0][2]; + real_t tmp36 = coords[1][1]*coords[0][2]; + real_t tmp37 = coords[0][1]*coords[1][2]; + + real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; + real_t detPfac = std::abs( detP ) / detP; + + INTEGRATE3D( 0, 0 ); + INTEGRATE3D( 0, 1 ); + INTEGRATE3D( 0, 2 ); + INTEGRATE3D( 0, 3 ); + + INTEGRATE3D( 1, 0 ); + INTEGRATE3D( 1, 1 ); + INTEGRATE3D( 1, 2 ); + INTEGRATE3D( 1, 3 ); + + INTEGRATE3D( 2, 0 ); + INTEGRATE3D( 2, 1 ); + INTEGRATE3D( 2, 2 ); + INTEGRATE3D( 2, 3 ); + + INTEGRATE3D( 3, 0 ); + INTEGRATE3D( 3, 1 ); + INTEGRATE3D( 3, 2 ); + INTEGRATE3D( 3, 3 ); + + INTEGRATE3D( 4, 0 ); + INTEGRATE3D( 4, 1 ); + INTEGRATE3D( 4, 2 ); + INTEGRATE3D( 4, 3 ); + + INTEGRATE3D( 5, 0 ); + INTEGRATE3D( 5, 1 ); + INTEGRATE3D( 5, 2 ); + INTEGRATE3D( 5, 3 ); + + INTEGRATE3D( 6, 0 ); + INTEGRATE3D( 6, 1 ); + INTEGRATE3D( 6, 2 ); + INTEGRATE3D( 6, 3 ); + + INTEGRATE3D( 7, 0 ); + INTEGRATE3D( 7, 1 ); + INTEGRATE3D( 7, 2 ); + INTEGRATE3D( 7, 3 ); + + INTEGRATE3D( 8, 0 ); + INTEGRATE3D( 8, 1 ); + INTEGRATE3D( 8, 2 ); + INTEGRATE3D( 8, 3 ); + + INTEGRATE3D( 9, 0 ); + INTEGRATE3D( 9, 1 ); + INTEGRATE3D( 9, 2 ); + INTEGRATE3D( 9, 3 ); + +// Undefine shape functions +#define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + }; +}; + + +// ========= +// div_z^T +// ========= +template <> +class P1ToP2Form_divt< 2 > : public P1ToP2FormHyTeG +{ + public: + + // --------- + // 2D Case + // --------- + void integrateAll( const std::array< Point3D, 3 >& coords, Matrixr< 6, 3 >& elMat ) const final + { + WALBERLA_ABORT( "P1ToP2Form_divt<2>.integrateAll not available in 2D!" ); + }; + + + // --------- + // 3D Case + // --------- void integrateAll( const std::array< Point3D, 4 >& coords, Matrixr< 10, 4 >& elMat ) const final { - WALBERLA_ABORT( "P1ToP2Form_divt not implemented for 3D, yet." ); + +// Set shape functions and their derivatives +#define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define DEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef DEFINE_P1_SHAPE_FUNCTIONS_TET +#undef DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET + +// Select cubature rule +#define CUBAPOINTS cubature::T4_points +#define CUBAWEIGHTS cubature::T4_weights + +// Executing quadrature rule +#define INTEGRATE3D( i, j ) \ + elMat( i, j ) = 0.0; \ + for ( uint_t k = 0; k < CUBAWEIGHTS.size(); k++ ) \ + { \ + real_t L2 = CUBAPOINTS[k][0]; \ + real_t L3 = CUBAPOINTS[k][1]; \ + real_t L4 = CUBAPOINTS[k][2]; \ + real_t L1 = 1.0 - L2 - L3 - L4; \ + \ + mappedPt[0] = ( coords[1][0] - coords[0][0] ) * L2 + ( coords[2][0] - coords[0][0] ) * L3 + ( coords[3][0] - coords[0][0] ) * L4 + coords[0][0]; \ + mappedPt[1] = ( coords[1][1] - coords[0][1] ) * L2 + ( coords[2][1] - coords[0][1] ) * L3 + ( coords[3][1] - coords[0][1] ) * L4 + coords[0][1]; \ + mappedPt[2] = ( coords[1][2] - coords[0][2] ) * L2 + ( coords[2][2] - coords[0][2] ) * L3 + ( coords[3][2] - coords[0][2] ) * L4 + coords[0][2]; \ + \ + Matrix3r DPsi; \ + geometryMap_->evalDF( mappedPt, DPsi ); \ + \ + real_t tmp12 = DPsi(0,0)*DPsi(1,1); \ + real_t tmp13 = DPsi(0,1)*DPsi(1,2); \ + real_t tmp14 = DPsi(0,2)*DPsi(1,0); \ + real_t tmp15 = DPsi(0,0)*DPsi(1,2); \ + real_t tmp16 = DPsi(0,1)*DPsi(1,0); \ + real_t tmp17 = DPsi(0,2)*DPsi(1,1); \ + real_t tmp18 = tmp12 - tmp16; \ + real_t tmp23 = tmp14 - tmp15; \ + real_t tmp28 = tmp13 - tmp17; \ + \ + real_t diff = DX1N ## i*(tmp18*(tmp19 + tmp20 - tmp21 - tmp22 + coords[2][0]*coords[3][1] - coords[3][0]*coords[2][1]) + tmp23*(tmp24 + tmp25 - tmp26 - tmp27 - coords[2][0]*coords[3][2] + coords[3][0]*coords[2][2]) + tmp28*(tmp29 + tmp30 - tmp31 - coords[2][1]*coords[0][2] + coords[2][1]*coords[3][2] - coords[3][1]*coords[2][2])) + DX2N ## i*(tmp18*(-tmp20 + tmp21 + tmp32 - tmp33 - coords[1][0]*coords[3][1] + coords[3][0]*coords[1][1]) + tmp23*(-tmp24 + tmp27 + tmp34 - tmp35 + coords[1][0]*coords[3][2] - coords[3][0]*coords[1][2]) + tmp28*(-tmp30 + tmp31 + tmp36 - tmp37 - coords[1][1]*coords[3][2] + coords[3][1]*coords[1][2])) + DX3N ## i*(tmp18*(-tmp19 + tmp22 - tmp32 + tmp33 + coords[1][0]*coords[2][1] - coords[2][0]*coords[1][1]) + tmp23*(-tmp25 + tmp26 - tmp34 + tmp35 - coords[1][0]*coords[2][2] + coords[2][0]*coords[1][2]) + tmp28*(-tmp29 - tmp36 + tmp37 + coords[1][1]*coords[2][2] + coords[2][1]*coords[0][2] - coords[2][1]*coords[1][2])); \ + \ + real_t detQ = DPsi(2,0)*tmp13 - DPsi(2,0)*tmp17 + DPsi(2,1)*tmp14 - DPsi(2,1)*tmp15 + DPsi(2,2)*tmp12 - DPsi(2,2)*tmp16; \ + real_t detQfac = std::abs( detQ ) / detQ; \ + elMat( i, j ) -= CUBAWEIGHTS[k] * detPfac * detQfac * diff * SF_M##j; \ + } + + // Quadrature point mapped to computational triangle + Point3D mappedPt; + + // pre-compute values including Jacobian determinant of inverse pull-back mapping + real_t tmp0 = -coords[0][0]; + real_t tmp1 = tmp0 + coords[1][0]; + real_t tmp2 = -coords[0][1]; + real_t tmp3 = tmp2 + coords[2][1]; + real_t tmp4 = -coords[0][2]; + real_t tmp5 = tmp4 + coords[3][2]; + real_t tmp6 = tmp0 + coords[2][0]; + real_t tmp7 = tmp2 + coords[3][1]; + real_t tmp8 = tmp4 + coords[1][2]; + real_t tmp9 = tmp0 + coords[3][0]; + real_t tmp10 = tmp2 + coords[1][1]; + real_t tmp11 = tmp4 + coords[2][2]; + real_t tmp19 = coords[0][0]*coords[2][1]; + real_t tmp20 = coords[3][0]*coords[0][1]; + real_t tmp21 = coords[0][0]*coords[3][1]; + real_t tmp22 = coords[2][0]*coords[0][1]; + real_t tmp24 = coords[0][0]*coords[3][2]; + real_t tmp25 = coords[2][0]*coords[0][2]; + real_t tmp26 = coords[0][0]*coords[2][2]; + real_t tmp27 = coords[3][0]*coords[0][2]; + real_t tmp29 = coords[0][1]*coords[2][2]; + real_t tmp30 = coords[3][1]*coords[0][2]; + real_t tmp31 = coords[0][1]*coords[3][2]; + real_t tmp32 = coords[1][0]*coords[0][1]; + real_t tmp33 = coords[0][0]*coords[1][1]; + real_t tmp34 = coords[0][0]*coords[1][2]; + real_t tmp35 = coords[1][0]*coords[0][2]; + real_t tmp36 = coords[1][1]*coords[0][2]; + real_t tmp37 = coords[0][1]*coords[1][2]; + + real_t detP = -tmp1*tmp11*tmp7 + tmp1*tmp3*tmp5 + tmp10*tmp11*tmp9 - tmp10*tmp5*tmp6 - tmp3*tmp8*tmp9 + tmp6*tmp7*tmp8; + real_t detPfac = std::abs( detP ) / detP; + + INTEGRATE3D( 0, 0 ); + INTEGRATE3D( 0, 1 ); + INTEGRATE3D( 0, 2 ); + INTEGRATE3D( 0, 3 ); + + INTEGRATE3D( 1, 0 ); + INTEGRATE3D( 1, 1 ); + INTEGRATE3D( 1, 2 ); + INTEGRATE3D( 1, 3 ); + + INTEGRATE3D( 2, 0 ); + INTEGRATE3D( 2, 1 ); + INTEGRATE3D( 2, 2 ); + INTEGRATE3D( 2, 3 ); + + INTEGRATE3D( 3, 0 ); + INTEGRATE3D( 3, 1 ); + INTEGRATE3D( 3, 2 ); + INTEGRATE3D( 3, 3 ); + + INTEGRATE3D( 4, 0 ); + INTEGRATE3D( 4, 1 ); + INTEGRATE3D( 4, 2 ); + INTEGRATE3D( 4, 3 ); + + INTEGRATE3D( 5, 0 ); + INTEGRATE3D( 5, 1 ); + INTEGRATE3D( 5, 2 ); + INTEGRATE3D( 5, 3 ); + + INTEGRATE3D( 6, 0 ); + INTEGRATE3D( 6, 1 ); + INTEGRATE3D( 6, 2 ); + INTEGRATE3D( 6, 3 ); + + INTEGRATE3D( 7, 0 ); + INTEGRATE3D( 7, 1 ); + INTEGRATE3D( 7, 2 ); + INTEGRATE3D( 7, 3 ); + + INTEGRATE3D( 8, 0 ); + INTEGRATE3D( 8, 1 ); + INTEGRATE3D( 8, 2 ); + INTEGRATE3D( 8, 3 ); + + INTEGRATE3D( 9, 0 ); + INTEGRATE3D( 9, 1 ); + INTEGRATE3D( 9, 2 ); + INTEGRATE3D( 9, 3 ); + +// Undefine shape functions +#define UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET +#define UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#include "ShapeFunctionMacros.hpp" +#undef UNDEFINE_P1_SHAPE_FUNCTIONS_TET +#undef UNDEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET }; }; diff --git a/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp b/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp index f6f6e0d82..fcbfdf243 100644 --- a/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp +++ b/src/hyteg/forms/form_hyteg_manual/P2ToP1FormDiv.hpp @@ -399,7 +399,6 @@ class P2ToP1Form_div< 1 > : public P2ToP1FormHyTeG // --------- void integrateAll( const std::array< Point3D, 4 >& coords, Matrixr< 4, 10 >& elMat ) const final { - // WALBERLA_ABORT( "P2ToP1Form_div<1> not implemented for 3D, yet." ); // Set shape functions and their derivatives #define DEFINE_P2_SHAPE_FUNCTION_DERIVATIVES_TET diff --git a/tests/hyteg/forms/HytegVsFenicsFormTest.cpp b/tests/hyteg/forms/HytegVsFenicsFormTest.cpp index 17df6d3ec..cef59a2e2 100644 --- a/tests/hyteg/forms/HytegVsFenicsFormTest.cpp +++ b/tests/hyteg/forms/HytegVsFenicsFormTest.cpp @@ -353,6 +353,24 @@ void run3DTestsWithoutBlending() P2ToP1Form_div< 2 >, Matrixr< 4, 10 >, 3 >( theTet, 1e-14 ); + + logSectionHeader( "P1ToP2 DivX^T Forms (3D)" ); + compareForms< P1ToP2FenicsForm< fenics::NoAssemble, p1_to_p2_tet_divt_tet_cell_integral_0_otherwise >, + P1ToP2Form_divt< 0 >, + Matrixr< 10, 4 >, + 3 >( theTet, 1e-14 ); + + logSectionHeader( "P1ToP2 DivY^T Forms (3D)" ); + compareForms< P1ToP2FenicsForm< fenics::NoAssemble, p1_to_p2_tet_divt_tet_cell_integral_1_otherwise >, + P1ToP2Form_divt< 1 >, + Matrixr< 10, 4 >, + 3 >( theTet, 1e-14 ); + + logSectionHeader( "P1ToP2 DivZ^T Forms (3D)" ); + compareForms< P1ToP2FenicsForm< fenics::NoAssemble, p1_to_p2_tet_divt_tet_cell_integral_2_otherwise >, + P1ToP2Form_divt< 2 >, + Matrixr< 10, 4 >, + 3 >( theTet, 1e-14 ); } void run3DTestsWithAffineMap() -- GitLab