diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp index 4576dedcb67ced609d8ffd226deb6e53f081410d..27f724ff3d5509b7d7bb6c242a0c7e8615bbcb0c 100644 --- a/src/scheme/FluxingAdvectionSolver.cpp +++ b/src/scheme/FluxingAdvectionSolver.cpp @@ -222,6 +222,28 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() 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& x1 = old_xr[face_to_node[1]]; + const Rd& x2 = old_xr[face_to_node[2]]; + + const Rd& x3 = new_xr[face_to_node[0]]; + const Rd& x4 = new_xr[face_to_node[1]]; + const Rd& x5 = new_xr[face_to_node[2]]; + + const Rd& a1 = x5 - x1; + 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); + 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");