diff --git a/src/geometry/SquareTransformation.hpp b/src/geometry/SquareTransformation.hpp
index d6474f33e4ecbaa4bcfba525dadc9fd38e95962d..b0e970c6eda673b48c5cf5c0d0a1ad7e73f8ef74 100644
--- a/src/geometry/SquareTransformation.hpp
+++ b/src/geometry/SquareTransformation.hpp
@@ -84,6 +84,25 @@ class SquareTransformation
return m_coefficients * X + m_shift;
}
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation(const TinyVector<2>& X) const
+ {
+ const auto& T = m_coefficients;
+ const double& x = X[0];
+ const double& y = X[1];
+
+ const TinyVector<Dimension> dxT{T(0, 0) + T(0, 2) * y, //
+ T(1, 0) + T(1, 2) * y, //
+ T(2, 0) + T(2, 2) * y};
+
+ const TinyVector<Dimension> dyT{T(0, 1) + T(0, 2) * x, //
+ T(1, 1) + T(1, 2) * x, //
+ T(2, 1) + T(2, 2) * x};
+
+ return crossProduct(dxT, dyT);
+ }
+
PUGS_INLINE double
areaVariationNorm(const TinyVector<2>& X) const
{
diff --git a/src/geometry/TriangleTransformation.hpp b/src/geometry/TriangleTransformation.hpp
index df8444d894a7ed56c0fe3825ffbf1cc14a742531..c49db368684d8c96aeb5061041828799e8509fac 100644
--- a/src/geometry/TriangleTransformation.hpp
+++ b/src/geometry/TriangleTransformation.hpp
@@ -66,6 +66,7 @@ class TriangleTransformation
private:
TinyMatrix<Dimension, 2> m_jacobian;
TinyVector<Dimension> m_shift;
+ TinyVector<Dimension> m_area_variation;
double m_area_variation_norm;
public:
@@ -76,6 +77,21 @@ class TriangleTransformation
return m_jacobian * x + m_shift;
}
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation() const
+ {
+ return m_area_variation;
+ }
+
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation(const TinyVector<2>&) const
+ {
+ return m_area_variation;
+ }
+
+ PUGS_INLINE
double
areaVariationNorm() const
{
@@ -92,7 +108,8 @@ class TriangleTransformation
m_shift = a;
- m_area_variation_norm = l2Norm(crossProduct(b - a, c - a));
+ m_area_variation = crossProduct(b - a, c - a);
+ m_area_variation_norm = l2Norm(m_area_variation);
}
~TriangleTransformation() = default;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index b9cfb57e151e5e8f3ba3df6b33b215e3805df345..f3d3768d9dd4a880386c66c20ca69fd5847be8e2 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -824,12 +824,7 @@ PolynomialReconstruction::_build(
return this->_build<CellCenterReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
}
case IntegrationMethodType::boundary: {
- if constexpr (MeshType::Dimension < 3) {
- return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh,
- discrete_function_variant_list);
- }
- std::clog << "falling back to element integration in 3d\n";
- [[fallthrough]];
+ return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
}
case IntegrationMethodType::element: {
return this->_build<ElementIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
index 5c7c5af3e4d39a9b903becfaeaade4a8aef4a6ca..58b8fcd5ea34a64bd8d5808c41ea95079c2fcf84 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
@@ -1,10 +1,13 @@
#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
#include <geometry/SymmetryUtils.hpp>
+#include <geometry/TriangleTransformation.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
@@ -38,7 +41,7 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh
return m_cell_face_is_reversed;
}
- SpecificDimensionalData(const Connectivity<2>& connectivity)
+ SpecificDimensionalData(const Connectivity<MeshType::Dimension>& connectivity)
: m_cell_to_face_matrix{connectivity.cellToFaceMatrix()},
m_face_to_node_matrix{connectivity.faceToNodeMatrix()},
m_cell_face_is_reversed{connectivity.cellFaceIsReversed()}
@@ -67,10 +70,68 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh
~SpecificDimensionalData() = default;
};
-template <MeshConcept MeshTypeT>
+template <>
template <typename ConformTransformationT>
void
-PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkBoundaryMean(
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<3>>::_computeEjkBoundaryMean(
+ const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double inv_Vi,
+ SmallArray<double>& mean_of_ejk)
+{
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const TinyVector<2> xi_q = quadrature.point(i_q);
+
+ const double area_variation_e1 = T.areaVariation(xi_q)[0] * inv_Vi;
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+ const double z_zj = X_Xj[2];
+
+ {
+ size_t k = 0;
+ m_tmp_array[k++] = x_xj * wq * area_variation_e1;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+ const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
+ const size_t nb_monoms = m_yz_row_size[i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_tmp_array[k] = y_yj * m_tmp_array[begin_i_y_1 + l];
+ }
+ }
+
+ for (size_t i_z = 1; i_z <= m_polynomial_degree; ++i_z) {
+ const size_t nb_y = m_yz_row_size[m_z_triangle_index[i_z]];
+ const size_t index_z = m_z_triangle_index[i_z];
+ const size_t index_z_1 = m_z_triangle_index[i_z - 1];
+ for (size_t i_y = 0; i_y < nb_y; ++i_y) {
+ const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
+ const size_t nb_monoms = m_yz_row_size[index_z + i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_tmp_array[k] = z_zj * m_tmp_array[begin_i_yz_1 + l];
+ }
+ }
+ }
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] += m_tmp_array[k];
+ }
+ }
+}
+
+template <>
+template <typename ConformTransformationT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::_computeEjkBoundaryMean(
const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
const ConformTransformationT& T,
const Rd& Xj,
@@ -97,7 +158,7 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
}
for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
- const size_t begin_i_y_1 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
+ const size_t begin_i_y_1 = m_y_row_index[i_y - 1];
for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
m_tmp_array[k++] = y_yj * m_tmp_array[begin_i_y_1 + l];
}
@@ -117,9 +178,6 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const CellId& cell_id,
SmallArray<double>& mean_of_ejk)
{
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
-
const double inv_Vi = 1. / m_Vj[cell_id];
const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
@@ -127,15 +185,60 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
mean_of_ejk.fill(0);
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = face_to_node_matrix[face_id];
- if (face_is_reversed[i_face]) {
- const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
- } else {
- const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+
+ if constexpr (MeshType::Dimension == 2) {
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ if (face_is_reversed[i_face]) {
+ const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+
+ const auto& triangle_quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(m_polynomial_degree + 1));
+ const auto& square_quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ switch (face_node_list.size()) {
+ case 3: {
+ if (face_is_reversed[i_face]) {
+ const TriangleTransformation<3> T{m_xr[face_node_list[2]], m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const TriangleTransformation<3> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]], m_xr[face_node_list[2]]};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ case 4: {
+ if (face_is_reversed[i_face]) {
+ const SquareTransformation<3> T{m_xr[face_node_list[3]], m_xr[face_node_list[2]], m_xr[face_node_list[1]],
+ m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const SquareTransformation<3> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]], m_xr[face_node_list[2]],
+ m_xr[face_node_list[3]]};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ throw NotImplementedError("invalid face type");
+ }
+ }
}
}
}
@@ -149,9 +252,6 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const CellId& cell_id,
SmallArray<double>& mean_of_ejk)
{
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
-
const double inv_Vi = 1. / m_Vj[cell_id];
const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
@@ -159,19 +259,70 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
mean_of_ejk.fill(0);
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = face_to_node_matrix[face_id];
-
- const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
- const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
-
- if (face_is_reversed[i_face]) {
- const LineTransformation<2> T{x1, x0};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
- } else {
- const LineTransformation<2> T{x0, x1};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ if constexpr (MeshType::Dimension == 2) {
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const LineTransformation<2> T{x1, x0};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{x0, x1};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+
+ const auto& triangle_quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(m_polynomial_degree + 1));
+ const auto& square_quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ switch (face_node_list.size()) {
+ case 3: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[2]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const TriangleTransformation<3> T{x2, x1, x0};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const TriangleTransformation<3> T{x0, x1, x2};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ case 4: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const SquareTransformation<3> T{x3, x2, x1, x0};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const SquareTransformation<3> T{x0, x1, x2, x3};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ throw NotImplementedError("invalid face type");
+ }
+ }
}
}
}
@@ -256,7 +407,7 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const CellId cell_j_id,
ShrinkMatrixView<SmallMatrix<double>>& A)
{
- if constexpr (MeshType::Dimension < 3) {
+ if constexpr (MeshType::Dimension <= 3) {
const auto& stencil_cell_list = m_stencil_array[cell_j_id];
const Rd& Xj = m_xj[cell_j_id];
@@ -326,6 +477,46 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
}
m_antiderivative_correction_coef = antiderivative_correction_coef;
+
+ if constexpr (MeshType::Dimension == 2) {
+ SmallArray<size_t> y_row_index(m_polynomial_degree + 1);
+
+ size_t i_y = 0;
+
+ y_row_index[i_y++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n > 1; --n, ++i_y) {
+ y_row_index[i_y] = y_row_index[i_y - 1] + n;
+ }
+
+ m_y_row_index = y_row_index;
+
+ } else if constexpr (MeshType::Dimension == 3) {
+ SmallArray<size_t> yz_row_index((m_polynomial_degree + 2) * (m_polynomial_degree + 1) / 2 + 1);
+ SmallArray<size_t> z_triangle_index(m_polynomial_degree + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_row_index[i_yz++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
+ z_triangle_index[i_z++] = i_yz - 1;
+ for (ssize_t i = n; i >= 1; --i) {
+ yz_row_index[i_yz] = yz_row_index[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ SmallArray<size_t> yz_row_size{yz_row_index.size() - 1};
+ for (size_t i = 0; i < yz_row_size.size(); ++i) {
+ yz_row_size[i] = yz_row_index[i + 1] - yz_row_index[i];
+ }
+
+ m_yz_row_index = yz_row_index;
+ m_z_triangle_index = z_triangle_index;
+ m_yz_row_size = yz_row_size;
+ }
}
template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::build(
@@ -336,6 +527,10 @@ template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuil
const CellId,
ShrinkMatrixView<SmallMatrix<double>>&);
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<3>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
Mesh<1>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
const size_t,
@@ -349,3 +544,10 @@ template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const SmallArray<const Rd>&,
const SmallArray<const Rd>&,
const CellToCellStencilArray&);
+
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<3>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
index 824aa902796cb0989c500b8b5eb05c8f0b4393fc..60707189b909e9d5f17ae2852be655cbfb877351 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -42,6 +42,14 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
const CellValue<const Rd> m_xj;
const NodeValue<const Rd> m_xr;
+ // 2D
+ SmallArray<const size_t> m_y_row_index;
+
+ // 3D
+ SmallArray<const size_t> m_yz_row_index;
+ SmallArray<const size_t> m_z_triangle_index;
+ SmallArray<const size_t> m_yz_row_size;
+
SmallArray<const double> m_antiderivative_correction_coef;
template <typename ConformTransformationT>
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
index 18306a8ecc635277f5e003d10faf53dfc41042a7..09e410811ed203610926bd605be3618e4219f396 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
@@ -77,13 +77,7 @@ PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<MeshTypeT>:
} else if constexpr (MeshType::Dimension == 3) {
static_assert(MeshType::Dimension == 3);
- const double detT = [&] {
- if constexpr (std::is_same_v<TetrahedronTransformation, std::decay_t<decltype(T)>>) {
- return T.jacobianDeterminant();
- } else {
- return T.jacobianDeterminant(xi_q);
- }
- }();
+ const double detT = T.jacobianDeterminant(xi_q);
const double x_xj = X_Xj[0];
const double y_yj = X_Xj[1];