Skip to content
Snippets Groups Projects
Commit ba3e2f14 authored by chantrait's avatar chantrait
Browse files

Merge branch 'develop' into feature/local-dt-fsi

parents 87aad460 1fe8f9ee
No related branches found
No related tags found
No related merge requests found
#include <scheme/FluxingAdvectionSolver.hpp> #include <scheme/FluxingAdvectionSolver.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/PrismTransformation.hpp>
#include <language/utils/EvaluateAtPoints.hpp> #include <language/utils/EvaluateAtPoints.hpp>
#include <language/utils/InterpolateItemArray.hpp> #include <language/utils/InterpolateItemArray.hpp>
#include <language/utils/InterpolateItemValue.hpp> #include <language/utils/InterpolateItemValue.hpp>
...@@ -322,6 +325,10 @@ template <> ...@@ -322,6 +325,10 @@ template <>
FaceValue<double> FaceValue<double>
FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
{ {
// due to the z component of the jacobian determinant, degree 3
// polynomials must be exactly integrated
const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(3));
const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix(); const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity()); FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
auto old_xr = m_old_mesh->xr(); auto old_xr = m_old_mesh->xr();
...@@ -329,31 +336,10 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() ...@@ -329,31 +336,10 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
parallel_for( parallel_for(
m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
const auto& face_to_node = face_to_node_matrix[face_id]; const auto& face_to_node = face_to_node_matrix[face_id];
if (face_to_node.size() == 4) {
const Rd& x0 = old_xr[face_to_node[0]];
const Rd& x1 = old_xr[face_to_node[1]];
const Rd& x2 = old_xr[face_to_node[2]];
const Rd& x3 = old_xr[face_to_node[3]];
const Rd& x4 = new_xr[face_to_node[0]]; const size_t face_nb_node = face_to_node.size();
const Rd& x5 = new_xr[face_to_node[1]];
const Rd& x6 = new_xr[face_to_node[2]];
const Rd& x7 = new_xr[face_to_node[3]];
const Rd& a1 = x6 - x1; if (face_nb_node == 3) {
const Rd& a2 = x6 - x3;
const Rd& a3 = x6 - x4;
const Rd& b1 = x7 - x0;
const Rd& b2 = x5 - x0;
const Rd& b3 = x2 - x0;
TinyMatrix<3> M1(a1 + b1, a2, a3);
TinyMatrix<3> M2(b1, a2 + b2, a3);
TinyMatrix<3> M3(a1, b2, a3 + b3);
algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12;
} else if (face_to_node.size() == 3) {
const Rd& x0 = old_xr[face_to_node[0]]; const Rd& x0 = old_xr[face_to_node[0]];
const Rd& x1 = old_xr[face_to_node[1]]; const Rd& x1 = old_xr[face_to_node[1]];
const Rd& x2 = old_xr[face_to_node[2]]; const Rd& x2 = old_xr[face_to_node[2]];
...@@ -362,21 +348,40 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() ...@@ -362,21 +348,40 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
const Rd& x4 = new_xr[face_to_node[1]]; const Rd& x4 = new_xr[face_to_node[1]];
const Rd& x5 = new_xr[face_to_node[2]]; const Rd& x5 = new_xr[face_to_node[2]];
const Rd& a1 = x5 - x1; double volume = 0;
const Rd& a2 = x5 - x2;
const Rd& a3 = x5 - x3;
const Rd& b1 = x5 - x0;
const Rd& b2 = x4 - x0;
const Rd& b3 = x2 - x0;
TinyMatrix<3> M1(a1 + b1, a2, a3); PrismTransformation T(x0, x1, x2, x3, x4, x5);
TinyMatrix<3> M2(b1, a2 + b2, a3); for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
TinyMatrix<3> M3(a1, b2, a3 + b3); volume += gauss.weight(i) * T.jacobianDeterminant(gauss.point(i));
}
algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12; algebraic_fluxing_volume[face_id] = volume;
} else { } else {
throw NotImplementedError("Not implemented for non quad faces"); Rd xg_old = zero;
for (size_t i_node = 0; i_node < face_nb_node; ++i_node) {
xg_old += old_xr[face_to_node[i_node]];
}
xg_old *= 1. / face_nb_node;
Rd xg_new = zero;
for (size_t i_node = 0; i_node < face_nb_node; ++i_node) {
xg_new += new_xr[face_to_node[i_node]];
}
xg_new *= 1. / face_nb_node;
double volume = 0;
for (size_t i_node = 0; i_node < face_nb_node; ++i_node) {
PrismTransformation T(xg_old, //
old_xr[face_to_node[i_node]], //
old_xr[face_to_node[(i_node + 1) % face_nb_node]], //
xg_new, //
new_xr[face_to_node[i_node]], //
new_xr[face_to_node[(i_node + 1) % face_nb_node]]);
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
volume += gauss.weight(i) * T.jacobianDeterminant(gauss.point(i));
}
}
algebraic_fluxing_volume[face_id] = volume;
} }
}); });
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment