diff --git a/src/hyteg/forms/form_hyteg_manual/P1ToP2FormDivT.hpp b/src/hyteg/forms/form_hyteg_manual/P1ToP2FormDivT.hpp index 3cf2cc76cead984f115fab5e41924a032038f842..6cd53024247f6862a0b952fc7ef7ad895fa1f896 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 ed0339ee650da6e67b5a438274c715c6db7dab28..fcbfdf2432c3bfab81d9b27a059c4b017cc29c25 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,170 @@ 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; \ + \ + 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 ## 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; \ + } + // 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( 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 +384,316 @@ 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." ); + +// 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 ## 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; \ + } + + // 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( 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; \ + \ + 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 ## 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; \ + } + + // 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( 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/src/hyteg/forms/form_hyteg_manual/ShapeFunctionMacros.hpp b/src/hyteg/forms/form_hyteg_manual/ShapeFunctionMacros.hpp index 008bc6622e1a6b9a9cd7822c4e91f73b34a214b3..9cde39b28461ff6d43b0fb64f5a0d11b535684ef 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 diff --git a/tests/hyteg/forms/HytegVsFenicsFormTest.cpp b/tests/hyteg/forms/HytegVsFenicsFormTest.cpp index 38d8016ca01f6648e408c6188c39f8dad00b6362..cef59a2e2e86e62efdf492155db076fcd084c92c 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,29 +294,34 @@ 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} )}; + 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)" ); @@ -333,8 +335,43 @@ 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, 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, 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, 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() { @@ -344,20 +381,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 @@ -365,23 +402,45 @@ 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 ) { // abort in case of common floating-point exceptions