From 7e3bb32b2b859f140f9eafa35f8fe24bb4a8b4ca Mon Sep 17 00:00:00 2001 From: Nicola Giuliani Date: Thu, 9 Jan 2020 09:36:10 +0100 Subject: [PATCH] using cross product but with problems still --- source/grid/tria_accessor.cc | 51 ++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 62919885f6..a7edc8283d 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1291,30 +1291,37 @@ namespace measure(const dealii::TriaAccessor<2, dim, 3> &accessor) { // In general the area can be computed as - // 0.25*(v_0+v_1-v_2-v_3)*(v_0-v_1+v_2-v_3) - - const Tensor<1, 3> piece_1 = accessor.vertex(0) + accessor.vertex(1) - - accessor.vertex(2) - accessor.vertex(3); - const Tensor<1, 3> piece_2 = accessor.vertex(0) - accessor.vertex(1) + - accessor.vertex(2) - accessor.vertex(3); - - return 0.25 * cross_product_3d(piece_1, piece_2).norm(); + // the integral of the cross product of the two tangential vectors + + // If we assume a bilinear patch parametrized in u and v we get that + // t_u = (v_1 - v_0) + v (v_3 - v_2 - v_1 + v_0) + // t_v = (v_2 - v_0) + u (v_3 - v_2 - v_1 + v_0) + // So t_u x t_v = (v_1 - v_0) x (v_2 - v_0) + u (v_1 - v_0) x (v_3 - v_2 - + // v_1 + v_0) + v (v_3 - v_2 - v_1 + v_0) x (v_2 - v_0) t_u x t_v = w_1 + u + // w_2 + v w_3 we can integrate the square norm (t_u x t_v) * (t_u x t_v) = + // w_1*w_1 + u^2 w_2*w_2 + v^2 w_3*w_3 + 2u w_1*w_2 + 2v w_1*w_3 + 2uv + // w_2*w_3 in u and v getting (between zero and one) w_1*w_1 + 1/3 w_2*w_2 + + // 1/3 w_3*w_3 + w_1*w_2 + w_1*w_3 + 1/2 w_2*w_3 + + const Tensor<1, 3> w_1 = + cross_product_3d(accessor.vertex(1) - accessor.vertex(0), + accessor.vertex(2) - accessor.vertex(0)); + const Tensor<1, 3> w_2 = + cross_product_3d(accessor.vertex(1) - accessor.vertex(0), + accessor.vertex(3) - accessor.vertex(2) - + accessor.vertex(1) + accessor.vertex(0)); + const Tensor<1, 3> w_3 = + cross_product_3d(accessor.vertex(3) - accessor.vertex(2) - + accessor.vertex(1) + accessor.vertex(0), + accessor.vertex(2) - accessor.vertex(0)); + + + return std::sqrt(scalar_product(w_1, w_1) + scalar_product(w_1, w_2) + + scalar_product(w_1, w_3) + 0.5 * scalar_product(w_2, w_3) + + 1. / 3 * scalar_product(w_2, w_2) + + 1. / 3 * scalar_product(w_3, w_3)); } - // // a 2d cell in 3d space - // double - // measure(const dealii::TriaAccessor<2, 2, 3> &accessor) - // { - // // In general the area can be computed as - // // 0.25*(v_0+v_1-v_2-v_3)*(v_0-v_1+v_2-v_3) - - // const Tensor<1, 3> piece_1 = accessor.vertex(0) + accessor.vertex(1) - - // accessor.vertex(2) - accessor.vertex(3); - // const Tensor<1, 3> piece_2 = accessor.vertex(0) - accessor.vertex(1) + - // accessor.vertex(2) - accessor.vertex(3); - // return 0.25 * cross_product_3d(piece_1, piece_2).norm(); - // } - template -- 2.39.5