diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp index d2f7cd45c96dd3ed299f9cf3fe7ca62a3aee3755..eba4c3278151291dea9a411e0e528cf6d19e08c2 100644 --- a/src/scheme/FluxingAdvectionSolver.cpp +++ b/src/scheme/FluxingAdvectionSolver.cpp @@ -192,7 +192,43 @@ template <> FaceValue<double> FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() { - throw NotImplementedError("ViensViensViens"); + const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix(); + FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity()); + auto old_xr = m_old_mesh->xr(); + auto new_xr = m_new_mesh->xr(); + parallel_for( + m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId 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[1]]; + const Rd& x5 = new_xr[face_to_node[0]]; + const Rd& x6 = new_xr[face_to_node[1]]; + const Rd& x7 = new_xr[face_to_node[0]]; + + const Rd& a1 = x6 - x1; + 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 { + throw NotImplementedError("Not implemented for non quad faces"); + } + }); + + return algebraic_fluxing_volume; } template <size_t Dimension>