From ee1a9a158d84cc13d91b1e0f31b6b61945472df5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 2 Dec 2021 00:03:29 +0100
Subject: [PATCH] Fix Quadrature formula for diamonds
now decompose into 2 pyramids or 2 tetrahedra in the case of "direct" formulae
---
src/language/utils/IntegrateOnCells.hpp | 14 +++++---------
src/scheme/CellIntegrator.hpp | 14 +++++---------
2 files changed, 10 insertions(+), 18 deletions(-)
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index 0290b3692..8b6008d73 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -488,12 +488,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
}
} else if (cell_to_node.size() == 6) {
- const auto qf = QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
+ const auto qf = QuadratureManager::instance().getPyramidFormula(quadrature_descriptor);
{ // top pyramid
- const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
- xr[cell_to_node[4]], //
- xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
- xr[cell_to_node[5]]);
+ const PyramidTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], xr[cell_to_node[5]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -503,10 +501,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
}
{ // bottom pyramid
- const CubeTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
- xr[cell_to_node[1]], //
- xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
- xr[cell_to_node[0]]);
+ const PyramidTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], xr[cell_to_node[0]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
diff --git a/src/scheme/CellIntegrator.hpp b/src/scheme/CellIntegrator.hpp
index 51758560a..e3a73a89c 100644
--- a/src/scheme/CellIntegrator.hpp
+++ b/src/scheme/CellIntegrator.hpp
@@ -425,12 +425,10 @@ class CellIntegrator
}
}
} else if (cell_to_node.size() == 6) {
- const auto qf = QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
+ const auto qf = QuadratureManager::instance().getPyramidFormula(quadrature_descriptor);
{ // top pyramid
- const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
- xr[cell_to_node[4]], //
- xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
- xr[cell_to_node[5]]);
+ const PyramidTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], xr[cell_to_node[5]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -438,10 +436,8 @@ class CellIntegrator
}
}
{ // bottom pyramid
- const CubeTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
- xr[cell_to_node[1]], //
- xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
- xr[cell_to_node[0]]);
+ const PyramidTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], xr[cell_to_node[0]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
--
GitLab