From: buerg Date: Mon, 13 Sep 2010 16:00:22 +0000 (+0000) Subject: Test nedelec_2 adapted and a bug in FE_Nedelec::get_subface_interpolation_matrix... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cd482535a7139c72571e502f957de0d6a606c0a0;p=dealii-svn.git Test nedelec_2 adapted and a bug in FE_Nedelec::get_subface_interpolation_matrix fixed. git-svn-id: https://svn.dealii.org/trunk@21947 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index d0bb20cd78..0a07eb57ea 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -1,16 +1,3 @@ -//--------------------------------------------------------------------------- -// $Id$ -// Version: $Name$ -// -// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors -// -// This file is subject to QPL and may not be distributed -// without copyright and license information. Please refer -// to the file deal.II/doc/license.html for the text and -// further information on this license. -// -//--------------------------------------------------------------------------- - #include #include #include @@ -1767,8 +1754,10 @@ const template void FE_Nedelec::get_subface_interpolation_matrix( - const FiniteElement&, const unsigned int, - FullMatrix&) const + const FiniteElement&, + const unsigned int, + FullMatrix&) +const { Assert (false, ExcNotImplemented ()); } @@ -1850,7 +1839,7 @@ FE_Nedelec::get_subface_interpolation_matrix // interpolation part. for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - const double weight = 2.0 * edge_quadrature.weight (q_point); + const double weight = edge_quadrature.weight (q_point); for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) interpolation_matrix (0, dof) @@ -1905,7 +1894,7 @@ FE_Nedelec::get_subface_interpolation_matrix const double tmp = 2.0 * quadrature_points[q_point] (0) + shift[subface]; const double weight - = std::sqrt (2.0 * edge_quadrature.weight (q_point)); + = std::sqrt (edge_quadrature.weight (q_point)); for (unsigned int i = 0; i < deg; ++i) assembling_matrix (i, q_point) @@ -1934,7 +1923,7 @@ FE_Nedelec::get_subface_interpolation_matrix = 2.0 * quadrature_points[q_point] (0) + shift[subface]; const double weight - = 2.0 * edge_quadrature.weight (q_point) + = edge_quadrature.weight (q_point) * (this->shape_value_component (dof, Point (0.0, quadrature_points[q_point] (0)), @@ -1978,7 +1967,7 @@ FE_Nedelec::get_subface_interpolation_matrix ++q_point) { const double - weight = 2.0 * edge_quadrature.weight (q_point); + weight = edge_quadrature.weight (q_point); for (unsigned int i = 0; i < 2; ++i) for (unsigned int dof = 0; dof < this->dofs_per_face; @@ -2072,12 +2061,9 @@ FE_Nedelec::get_subface_interpolation_matrix // face. std::vector > > edge_quadrature_points_full_dim - (GeometryInfo::lines_per_face); + (GeometryInfo::lines_per_face, + std::vector > (n_edge_points)); - for (unsigned int line = 0; - line < GeometryInfo::lines_per_face; ++line) - edge_quadrature_points_full_dim.resize (n_edge_points); - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { @@ -2097,35 +2083,32 @@ FE_Nedelec::get_subface_interpolation_matrix 0.5, 0.0); } - Vector solution (deg); - Vector system_rhs (deg); - - for (unsigned int dof = 0; dof < this->dofs_per_face; - ++dof) - { // Set up the system matrix. // This can be used for all // edges. - for (unsigned int q_point = 0; - q_point < n_edge_points; ++q_point) - { - const double tmp - = 2.0 * edge_quadrature_points[q_point] (0); - const double weight - = std::sqrt (2.0 * edge_quadrature.weight - (q_point)); - - for (unsigned int i = 0; i < deg; ++i) - assembling_matrix (i, q_point) - = weight - * lobatto_polynomials_grad[i + 1].value - (tmp); - } + for (unsigned int q_point = 0; + q_point < n_edge_points; ++q_point) + { + const double tmp + = 2.0 * edge_quadrature_points[q_point] (0); + const double weight + = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < deg; ++i) + assembling_matrix (i, q_point) + = weight * lobatto_polynomials_grad[i + 1].value + (tmp); + } - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.invert (system_matrix); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.invert (system_matrix); + Vector solution (deg); + Vector system_rhs (deg); + + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) @@ -2137,7 +2120,7 @@ FE_Nedelec::get_subface_interpolation_matrix q_point < n_edge_points; ++q_point) { const double right_hand_side_value - = std::sqrt (2.0 * edge_quadrature.weight + = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value_component (this->face_to_cell_index (dof, 4), @@ -2168,9 +2151,44 @@ FE_Nedelec::get_subface_interpolation_matrix = solution (i); } - assembling_matrix.reinit (deg * this->degree, - n_face_points); - system_rhs.reinit (assembling_matrix.m ()); + assembling_matrix.reinit (deg * this->degree, + n_face_points); + + for (unsigned int q_point = 0; q_point < n_face_points; + ++q_point) + { + const Point quadrature_point + (2.0 * face_quadrature_points[q_point] (0), + 2.0 * face_quadrature_points[q_point] (1)); + const double weight + = std::sqrt (face_quadrature.weight (q_point)); + + for (unsigned int i = 0; i <= deg; ++i) + { + const double tmp + = weight * legendre_polynomials[i].value + (quadrature_point (0)); + + for (unsigned int j = 0; j < deg; ++j) + assembling_matrix (i * deg + j, q_point) + = tmp * lobatto_polynomials[j + 2].value + (quadrature_point (1)); + } + } + + system_matrix.reinit (assembling_matrix.m (), + assembling_matrix.m ()); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + solution.reinit (system_matrix_inv.m ()); + system_rhs.reinit (system_matrix_inv.m ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) + { system_rhs = 0; // Now we project the remaining @@ -2204,42 +2222,23 @@ FE_Nedelec::get_subface_interpolation_matrix quadrature_point, 1); right_hand_side_value - *= 4.0 * face_quadrature.weight (q_point); + *= face_quadrature.weight (q_point); - const double weight - = std::sqrt (4.0 * face_quadrature.weight - (q_point)); - for (unsigned int i = 0; i <= deg; ++i) { - const double L_i - = legendre_polynomials[i].value - (quadrature_point (0)); - const double tmp1 = weight * L_i; - const double tmp2 - = right_hand_side_value * L_i; + const double tmp + = right_hand_side_value + * legendre_polynomials[i].value + (quadrature_point (0)); for (unsigned int j = 0; j < deg; ++j) - { - const double l_j - = lobatto_polynomials[j + 2].value - (quadrature_point (1)); - - assembling_matrix (i * deg + j, q_point) - = tmp1 * l_j; - system_rhs (i * deg + j) += tmp2 * l_j; - } + system_rhs (i * deg + j) + += tmp + * lobatto_polynomials[j + 2].value + (quadrature_point (1)); } } - system_matrix.reinit (assembling_matrix.m (), - assembling_matrix.m ()); - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.reinit (system_matrix.m (), - system_matrix.m ()); - system_matrix_inv.invert (system_matrix); - solution.reinit (system_matrix_inv.m ()); system_matrix_inv.vmult (solution, system_rhs); // Add the computed values @@ -2284,15 +2283,14 @@ FE_Nedelec::get_subface_interpolation_matrix quadrature_point, 0); right_hand_side_value - *= 4.0 * face_quadrature.weight (q_point); + *= face_quadrature.weight (q_point); for (unsigned int i = 0; i <= deg; ++i) { - const double L_i - = legendre_polynomials[i].value - (quadrature_point (0)); - const double tmp = right_hand_side_value - * L_i; + const double tmp + = right_hand_side_value + * legendre_polynomials[i].value + (quadrature_point (0)); for (unsigned int j = 0; j < deg; ++j) system_rhs (i * deg + j) @@ -2341,27 +2339,29 @@ FE_Nedelec::get_subface_interpolation_matrix for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { -// const double weight -// = 2.0 * edge_quadrature_x.weight (q_point); + const double weight + = edge_quadrature_x.weight (q_point); for (unsigned int i = 0; i < 2; ++i) for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { interpolation_matrix (i * source_fe.degree, dof) - += this->shape_value_component - (this->face_to_cell_index (dof, 4), - Point - (0.5 * (i + 1), - edge_quadrature_y_points[q_point] (0), 0.0), - 1); + += weight + * this->shape_value_component + (this->face_to_cell_index (dof, 4), + Point + (0.5 * (i + 1), + edge_quadrature_y_points[q_point] (0), 0.0), + 1); interpolation_matrix ((i + 2) * source_fe.degree, dof) - += this->shape_value_component - (this->face_to_cell_index (dof, 4), - Point - (edge_quadrature_x_points[q_point] (0), - 0.5 * i, 0.0), 0); + += weight + * this->shape_value_component + (this->face_to_cell_index (dof, 4), + Point + (edge_quadrature_x_points[q_point] (0), + 0.5 * i, 0.0), 0); } } @@ -2434,22 +2434,16 @@ FE_Nedelec::get_subface_interpolation_matrix // weights on the four // edges of a face. std::vector > edge_quadrature_points - (GeometryInfo::lines_per_face); + (GeometryInfo::lines_per_face, + std::vector (n_edge_points)); std::vector > edge_quadrature_weights - (GeometryInfo::lines_per_face); + (GeometryInfo::lines_per_face, + std::vector (n_edge_points)); std::vector > > edge_quadrature_points_full_dim - (GeometryInfo::lines_per_face); - - for (unsigned int line = 0; - line < GeometryInfo::lines_per_face; ++line) - { - edge_quadrature_points.resize (n_edge_points); - edge_quadrature_points_full_dim.resize - (n_edge_points); - edge_quadrature_weights.resize (n_edge_points); - } + (GeometryInfo::lines_per_face, + std::vector > (n_edge_points)); for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) @@ -2480,46 +2474,41 @@ FE_Nedelec::get_subface_interpolation_matrix (edge_quadrature_x_points[q_point] (0), 0.5, 0.0); edge_quadrature_weights[0][q_point] - = std::sqrt (2.0 * edge_quadrature_y.weight - (q_point)); + = std::sqrt (edge_quadrature_y.weight (q_point)); edge_quadrature_weights[1][q_point] = edge_quadrature_weights[0][q_point]; edge_quadrature_weights[2][q_point] - = std::sqrt (2.0 * edge_quadrature_x.weight - (q_point)); + = std::sqrt (edge_quadrature_x.weight (q_point)); edge_quadrature_weights[3][q_point] = edge_quadrature_weights[2][q_point]; } - Vector system_rhs (system_matrix.m ()); - Vector solution (system_rhs.size ()); - - for (unsigned int dof = 0; dof < this->dofs_per_face; - ++dof) - { // Set up the system matrix. // This can be used for all // edges. - for (unsigned int q_point = 0; - q_point < n_edge_points; ++q_point) - { - const double tmp - = 2.0 * edge_quadrature_y_points[q_point] (0); - const double weight - = std::sqrt (2.0 * edge_quadrature_y.weight - (q_point)); - - for (unsigned int i = 0; i < deg; ++i) - assembling_matrix (i, q_point) - = weight - * lobatto_polynomials_grad[i + 1].value - (tmp); - } + for (unsigned int q_point = 0; + q_point < n_edge_points; ++q_point) + { + const double tmp + = 2.0 * edge_quadrature_y_points[q_point] (0); + const double weight + = std::sqrt (edge_quadrature_y.weight (q_point)); + + for (unsigned int i = 0; i < deg; ++i) + assembling_matrix (i, q_point) + = weight + * lobatto_polynomials_grad[i + 1].value (tmp); + } - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.invert (system_matrix); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.invert (system_matrix); + Vector system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) for (unsigned int line = 0; line < GeometryInfo::lines_per_cell; ++line) @@ -2565,9 +2554,45 @@ FE_Nedelec::get_subface_interpolation_matrix // functions. First on the // horizontal ones, then on // the vertical ones. - assembling_matrix.reinit (deg * this->degree, - n_face_points); - system_rhs.reinit (assembling_matrix.m ()); + assembling_matrix.reinit (deg * this->degree, + n_face_points); + + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + const Point quadrature_point + (2.0 * face_quadrature_points[q_point] (0) + - 1.0, + 2.0 * face_quadrature_points[q_point] (1)); + const double weight + = std::sqrt (face_quadrature.weight (q_point)); + + for (unsigned int i = 0; i <= deg; ++i) + { + const double tmp + = weight * legendre_polynomials[i].value + (quadrature_point (0)); + + for (unsigned int j = 0; j < deg; ++j) + assembling_matrix (i * deg + j, q_point) + = tmp * lobatto_polynomials[j + 2].value + (quadrature_point (1)); + } + } + + system_matrix.reinit (assembling_matrix.m (), + assembling_matrix.m ()); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + solution.reinit (system_matrix_inv.m ()); + system_rhs.reinit (assembling_matrix.m ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) + { system_rhs = 0; for (unsigned int q_point = 0; @@ -2597,42 +2622,23 @@ FE_Nedelec::get_subface_interpolation_matrix quadrature_point, 1); right_hand_side_value - *= 4.0 * face_quadrature.weight (q_point); - - const double weight - = std::sqrt (4.0 * face_quadrature.weight - (q_point)); + *= face_quadrature.weight (q_point); for (unsigned int i = 0; i <= deg; ++i) { - const double L_i - = legendre_polynomials[i].value - (quadrature_point (0)); - const double tmp1 = weight * L_i; - const double tmp2 = right_hand_side_value - * L_i; + const double tmp + = right_hand_side_value + * legendre_polynomials[i].value + (quadrature_point (0)); for (unsigned int j = 0; j < deg; ++j) - { - const double l_j - = lobatto_polynomials[j + 2].value + system_rhs (i * deg + j) + = tmp + * lobatto_polynomials[j + 2].value (quadrature_point (1)); - - assembling_matrix (i * deg + j, q_point) - = tmp1 * l_j; - system_rhs (i * deg + j) += tmp2 * l_j; - } } } - system_matrix.reinit (assembling_matrix.m (), - assembling_matrix.m ()); - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.reinit (system_matrix.m (), - system_matrix.m ()); - system_matrix_inv.invert (system_matrix); - solution.reinit (system_matrix_inv.m ()); system_matrix_inv.vmult (solution, system_rhs); // Add the computed values @@ -2678,15 +2684,14 @@ FE_Nedelec::get_subface_interpolation_matrix quadrature_point, 0); right_hand_side_value - *= 4.0 * face_quadrature.weight (q_point); + *= face_quadrature.weight (q_point); for (unsigned int i = 0; i <= deg; ++i) { - const double L_i - = legendre_polynomials[i].value - (quadrature_point (0)); const double tmp - = right_hand_side_value * L_i; + = right_hand_side_value + * legendre_polynomials[i].value + (quadrature_point (0)); for (unsigned int j = 0; j < deg; ++j) system_rhs (i * deg + j) @@ -2735,27 +2740,29 @@ FE_Nedelec::get_subface_interpolation_matrix for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { -// const double weight -// = 2.0 * edge_quadrature_x.weight (q_point); + const double weight + = edge_quadrature_x.weight (q_point); for (unsigned int i = 0; i < 2; ++i) for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { interpolation_matrix (i * source_fe.degree, dof) - += this->shape_value_component - (this->face_to_cell_index (dof, 4), - Point - (0.5 * i, - edge_quadrature_y_points[q_point] (0), 0.0), - 1); + += weight + * this->shape_value_component + (this->face_to_cell_index (dof, 4), + Point + (0.5 * i, + edge_quadrature_y_points[q_point] (0), 0.0), + 1); interpolation_matrix ((i + 2) * source_fe.degree, dof) - += this->shape_value_component - (this->face_to_cell_index (dof, 4), - Point - (edge_quadrature_x_points[q_point] (0), - 0.5 * (i + 1), 0.0), 0); + += weight + * this->shape_value_component + (this->face_to_cell_index (dof, 4), + Point + (edge_quadrature_x_points[q_point] (0), + 0.5 * (i + 1), 0.0), 0); } } @@ -2824,22 +2831,16 @@ FE_Nedelec::get_subface_interpolation_matrix // edges of a face. std::vector > edge_quadrature_points - (GeometryInfo::lines_per_face); + (GeometryInfo::lines_per_face, + std::vector (n_edge_points)); std::vector > edge_quadrature_weights - (GeometryInfo::lines_per_face); + (GeometryInfo::lines_per_face, + std::vector (n_edge_points)); std::vector > > edge_quadrature_points_full_dim - (GeometryInfo::lines_per_face); - - for (unsigned int line = 0; - line < GeometryInfo::lines_per_face; ++line) - { - edge_quadrature_points.resize (n_edge_points); - edge_quadrature_points_full_dim.resize - (n_edge_points); - edge_quadrature_weights.resize (n_edge_points); - } + (GeometryInfo::lines_per_face, + std::vector > (n_edge_points)); for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) @@ -2870,47 +2871,44 @@ FE_Nedelec::get_subface_interpolation_matrix (edge_quadrature_x_points[q_point] (0), 1.0, 0.0); edge_quadrature_weights[0][q_point] - = std::sqrt (2.0 * edge_quadrature_y.weight - (q_point)); + = std::sqrt (edge_quadrature_y.weight (q_point)); edge_quadrature_weights[1][q_point] = edge_quadrature_weights[0][q_point]; edge_quadrature_weights[2][q_point] - = std::sqrt (2.0 * edge_quadrature_x.weight - (q_point)); + = std::sqrt (edge_quadrature_x.weight (q_point)); edge_quadrature_weights[3][q_point] = edge_quadrature_weights[2][q_point]; } - Vector system_rhs (system_matrix.m ()); - Vector solution (system_rhs.size ()); - - for (unsigned int dof = 0; dof < this->dofs_per_face; - ++dof) - { // Set up the system matrix. // This can be used for all // edges. - for (unsigned int q_point = 0; - q_point < n_edge_points; ++q_point) - { - const double weight - = std::sqrt (2.0 * edge_quadrature_y.weight - (q_point)); - const double tmp - = 2.0 * edge_quadrature_y_points[q_point] (0) - - 1.0; - - for (unsigned int i = 0; i < deg; ++i) - assembling_matrix (i, q_point) - = weight - * lobatto_polynomials_grad[i + 1].value - (tmp); - } + for (unsigned int q_point = 0; + q_point < n_edge_points; ++q_point) + { + const double weight + = std::sqrt (edge_quadrature_y.weight + (q_point)); + const double tmp + = 2.0 * edge_quadrature_y_points[q_point] (0) + - 1.0; + + for (unsigned int i = 0; i < deg; ++i) + assembling_matrix (i, q_point) + = weight + * lobatto_polynomials_grad[i + 1].value + (tmp); + } - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.invert (system_matrix); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.invert (system_matrix); + Vector system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) for (unsigned int line = 0; line < GeometryInfo::lines_per_cell; ++line) @@ -2950,16 +2948,50 @@ FE_Nedelec::get_subface_interpolation_matrix = solution (i); } - assembling_matrix.reinit (deg * this->degree, - n_face_points); - system_rhs.reinit (assembling_matrix.m ()); - system_rhs = 0; + assembling_matrix.reinit (deg * this->degree, + n_face_points); + + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + const double weight + = std::sqrt (face_quadrature.weight (q_point)); + const Point quadrature_point + (2.0 * face_quadrature_points[q_point] (0), + 2.0 * face_quadrature_points[q_point] (1) - 1.0); + + for (unsigned int i = 0; i <= deg; ++i) + { + const double tmp + = weight * legendre_polynomials[i].value + (quadrature_point (0)); + + for (unsigned int j = 0; j < deg; ++j) + assembling_matrix (i * deg + j, q_point) + = tmp * lobatto_polynomials[j + 2].value + (quadrature_point (1)); + } + } + + system_matrix.reinit (assembling_matrix.m (), + assembling_matrix.m ()); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + solution.reinit (system_matrix_inv.m ()); + system_rhs.reinit (assembling_matrix.m ()); + system_rhs = 0; // Now we project the remaining // part on the face shape // functions. First on the // horizontal ones, then on // the vertical ones. + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) + { for (unsigned int q_point = 0; q_point < n_face_points; ++q_point) { @@ -2986,42 +3018,23 @@ FE_Nedelec::get_subface_interpolation_matrix quadrature_point, 1); right_hand_side_value - *= 4.0 * face_quadrature.weight (q_point); + *= face_quadrature.weight (q_point); - const double weight - = std::sqrt (4.0 * reference_face_quadrature.weight - (q_point)); - for (unsigned int i = 0; i <= deg; ++i) { - const double L_i - = legendre_polynomials[i].value - (quadrature_point (0)); - const double tmp1 = weight * L_i; - const double tmp2 - = right_hand_side_value * L_i; + const double tmp + = right_hand_side_value + * legendre_polynomials[i].value + (quadrature_point (0)); for (unsigned int j = 0; j < deg; ++j) - { - const double l_j - = lobatto_polynomials[j + 2].value - (quadrature_point (1)); - - assembling_matrix (i * deg + j, q_point) - = tmp1 * l_j; - system_rhs (i * deg + j) += tmp2 * l_j; - } + system_rhs (i * deg + j) + += tmp + * lobatto_polynomials[j + 2].value + (quadrature_point (1)); } } - system_matrix.reinit (assembling_matrix.m (), - assembling_matrix.m ()); - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.reinit (system_matrix.m (), - system_matrix.m ()); - system_matrix_inv.invert (system_matrix); - solution.reinit (system_matrix_inv.m ()); system_matrix_inv.vmult (solution, system_rhs); // Add the computed values @@ -3065,14 +3078,15 @@ FE_Nedelec::get_subface_interpolation_matrix (i * source_fe.degree + j, quadrature_point, 0); - right_hand_side_value *= 4.0 * face_quadrature.weight (q_point); + right_hand_side_value *= face_quadrature.weight + (q_point); for (unsigned int i = 0; i <= deg; ++i) { - const double L_i - = legendre_polynomials[i].value - (quadrature_point (0)); - const double tmp = right_hand_side_value * L_i; + const double tmp + = right_hand_side_value + * legendre_polynomials[i].value + (quadrature_point (0)); for (unsigned int j = 0; j < deg; ++j) system_rhs (i * deg + j) @@ -3115,27 +3129,28 @@ FE_Nedelec::get_subface_interpolation_matrix for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { -// const double weight -// = 2.0 * edge_quadrature.weight (q_point); + const double weight = edge_quadrature.weight (q_point); for (unsigned int i = 0; i < 2; ++i) for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { interpolation_matrix (i * source_fe.degree, dof) - += this->shape_value_component - (this->face_to_cell_index (dof, 4), - Point - (0.5 * (i + 1), - edge_quadrature_points[q_point] (0), 0.0), - 1); + += weight + * this->shape_value_component + (this->face_to_cell_index (dof, 4), + Point + (0.5 * (i + 1), + edge_quadrature_points[q_point] (0), 0.0), + 1); interpolation_matrix ((i + 2) * source_fe.degree, dof) - += this->shape_value_component - (this->face_to_cell_index (dof, 4), - Point - (edge_quadrature_points[q_point] (0), - 0.5 * (i + 1), 0.0), 0); + += weight + * this->shape_value_component + (this->face_to_cell_index (dof, 4), + Point + (edge_quadrature_points[q_point] (0), + 0.5 * (i + 1), 0.0), 0); } } @@ -3209,12 +3224,8 @@ FE_Nedelec::get_subface_interpolation_matrix // face. std::vector > > edge_quadrature_points_full_dim - (GeometryInfo::lines_per_face); - - for (unsigned int line = 0; - line < GeometryInfo::lines_per_face; ++line) - edge_quadrature_points_full_dim.resize - (n_edge_points); + (GeometryInfo::lines_per_face, + std::vector > (n_edge_points)); for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) @@ -3233,36 +3244,34 @@ FE_Nedelec::get_subface_interpolation_matrix 1.0, 0.0); } - Vector system_rhs (system_matrix.m ()); - Vector solution (system_rhs.size ()); - - for (unsigned int dof = 0; dof < this->dofs_per_face; - ++dof) - { // Set up the system matrix. // This can be used for all // edges. - for (unsigned int q_point = 0; - q_point < n_edge_points; ++q_point) - { - const double tmp - = 2.0 * edge_quadrature_points[q_point] (0) - - 1.0; - const double weight - = std::sqrt (2.0 * edge_quadrature.weight - (q_point)); - - for (unsigned int i = 0; i < deg; ++i) - assembling_matrix (i, q_point) - = weight - * lobatto_polynomials_grad[i + 1].value - (tmp); - } + for (unsigned int q_point = 0; + q_point < n_edge_points; ++q_point) + { + const double tmp + = 2.0 * edge_quadrature_points[q_point] (0) + - 1.0; + const double weight + = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < deg; ++i) + assembling_matrix (i, q_point) + = weight + * lobatto_polynomials_grad[i + 1].value + (tmp); + } - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.invert (system_matrix); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.invert (system_matrix); + Vector system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) for (unsigned int line = 0; line < GeometryInfo::lines_per_cell; ++line) @@ -3274,8 +3283,8 @@ FE_Nedelec::get_subface_interpolation_matrix q_point < n_edge_points; ++q_point) { const double right_hand_side_value - = std::sqrt (2.0 * edge_quadrature.weight - (q_point)) + = std::sqrt (edge_quadrature.weight + (q_point)) * (this->shape_value_component (this->face_to_cell_index (dof, 4), edge_quadrature_points_full_dim[line][q_point], @@ -3306,11 +3315,46 @@ FE_Nedelec::get_subface_interpolation_matrix = solution (i); } - assembling_matrix.reinit (deg * this->degree, - n_face_points); - system_rhs.reinit (assembling_matrix.m ()); - system_rhs = 0; + assembling_matrix.reinit (deg * this->degree, + n_face_points); + + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + const double weight + = std::sqrt (face_quadrature.weight + (q_point)); + const Point quadrature_point + (2.0 * face_quadrature_points[q_point] (0) - 1.0, + 2.0 * face_quadrature_points[q_point] (1) - 1.0); + + for (unsigned int i = 0; i <= deg; ++i) + { + const double tmp + = weight * legendre_polynomials[i].value + (quadrature_point (0)); + + for (unsigned int j = 0; j < deg; ++j) + assembling_matrix (i * deg + j, q_point) + = tmp * lobatto_polynomials[j + 2].value + (quadrature_point (1)); + } + } + + system_matrix.reinit (assembling_matrix.m (), + assembling_matrix.m ()); + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + solution.reinit (system_matrix.m ()); + system_rhs.reinit (assembling_matrix.m ()); + system_rhs = 0; + for (unsigned int dof = 0; dof < this->dofs_per_face; + ++dof) + { // Now we project the remaining // part on the face shape // functions. First on the @@ -3343,42 +3387,23 @@ FE_Nedelec::get_subface_interpolation_matrix quadrature_point, 1); right_hand_side_value - *= 4.0 * face_quadrature.weight (q_point); - - const double weight - = std::sqrt (4.0 * face_quadrature.weight - (q_point)); + *= face_quadrature.weight (q_point); for (unsigned int i = 0; i <= deg; ++i) { - const double L_i - = legendre_polynomials[i].value - (quadrature_point (0)); - const double tmp1 = weight * L_i; - const double tmp2 = right_hand_side_value - * L_i; + const double tmp + = right_hand_side_value + * legendre_polynomials[i].value + (quadrature_point (0)); for (unsigned int j = 0; j < deg; ++j) - { - const double l_j - = lobatto_polynomials[j + 2].value - (quadrature_point (1)); - - assembling_matrix (i * deg + j, q_point) - = tmp1 * l_j; - system_rhs (i * deg + j) += tmp2 * l_j; - } + system_rhs (i * deg + j) + += tmp + * lobatto_polynomials[j + 2].value + (quadrature_point (1)); } } - system_matrix.reinit (assembling_matrix.m (), - assembling_matrix.m ()); - assembling_matrix.mTmult (system_matrix, - assembling_matrix); - system_matrix_inv.reinit (system_matrix.m (), - system_matrix.m ()); - system_matrix_inv.invert (system_matrix); - solution.reinit (system_matrix.m ()); system_matrix_inv.vmult (solution, system_rhs); // Add the computed values @@ -3424,7 +3449,7 @@ FE_Nedelec::get_subface_interpolation_matrix quadrature_point, 0); right_hand_side_value - *= 4.0 * face_quadrature.weight (q_point); + *= face_quadrature.weight (q_point); for (unsigned int i = 0; i <= deg; ++i) { @@ -3825,7 +3850,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, local_dofs[(i + 8) * this->degree] = 0.0; } - if (offset+1 < dim) + if (offset + 1 < dim) { for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j < 2; ++j) @@ -3923,7 +3948,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, system_rhs = 0; if ((((line == 0) || (line == 1) || (line == 4) || - (line == 5)) && (offset+1 < dim)) || + (line == 5)) && (offset + 1 < dim)) || (((line == 2) || (line == 3) || (line == 6) || (line == 7)) && (offset == 0)) || (line > 7)) { @@ -4010,7 +4035,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, { case 0: { - if (offset+1 < dim) + if (offset + 1 < dim) { // Set up the right hand side // for the horizontal shape @@ -4114,7 +4139,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, case 1: { - if (offset+1 < dim) + if (offset + 1 < dim) { // Set up the right hand side // for the horizontal shape @@ -4437,7 +4462,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, case 4: { - if (offset+1 < dim) + if (offset + 1 < dim) { // Set up the right hand side // for the horizontal shape @@ -4550,7 +4575,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, } default: - if (offset+1 < dim) + if (offset + 1 < dim) { // Set up the right hand side // for the horizontal shape @@ -4733,7 +4758,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, system_rhs.reinit (system_matrix_inv.m ()); solution.reinit (system_matrix.m ()); - if (offset+1 < dim) + if (offset + 1 < dim) { if (offset == 0) { @@ -5569,7 +5594,8 @@ const + GeometryInfo::lines_per_cell * n_edge_points]; - for (int i = 2;i < (int) GeometryInfo::lines_per_face; ++i) + for (unsigned int i = 2; + i < GeometryInfo::lines_per_face; ++i) for (unsigned int j = 0; j <= deg; ++j) tmp -= local_dofs[edge_indices[face][i] * this->degree + j] diff --git a/tests/fe/nedelec_2.cc b/tests/fe/nedelec_2.cc index eb089f25de..e6abc991d3 100644 --- a/tests/fe/nedelec_2.cc +++ b/tests/fe/nedelec_2.cc @@ -28,12 +28,12 @@ template void -plot (const Triangulation &tr) +plot (const Triangulation &tr, const unsigned int p) { deallog << dim << "d, " << tr.n_active_cells() << " CELLS" << std::endl; - FE_Nedelec fe_ned(1); + FE_Nedelec fe_ned (p); DoFHandler dof(tr); dof.distribute_dofs(fe_ned); @@ -87,7 +87,7 @@ plot (const Triangulation &tr) template inline void -check () +check (const unsigned int p) { Triangulation tr; GridGenerator::hyper_cube(tr, 0., 1.); @@ -95,12 +95,12 @@ check () // first everything on a // once-refined grid tr.refine_global (1); - plot (tr); + plot (tr, p); // then same with one cell refined tr.begin_active()->set_refine_flag (); tr.execute_coarsening_and_refinement (); - plot (tr); + plot (tr, p); } @@ -113,9 +113,12 @@ main() deallog.attach(logfile); deallog.depth_console(0); deallog.threshold_double(1.e-10); - - check<2>(); - check<3>(); + deallog << "Degree 0:" << std::endl; + check<2> (0); + check<3> (0); + deallog << "Degree 1:" << std::endl; + check<2> (1); + check<3> (1); return 0; } diff --git a/tests/fe/nedelec_2/cmp/generic b/tests/fe/nedelec_2/cmp/generic index 516283be33..8097c289be 100644 --- a/tests/fe/nedelec_2/cmp/generic +++ b/tests/fe/nedelec_2/cmp/generic @@ -1,4 +1,5 @@ +DEAL::Degree 0: DEAL::2d, 4 CELLS DEAL:: CELL 1.0 DEAL::, xq=0.00 0.00, f=[4.00 0] @@ -299,3 +300,304 @@ DEAL::, xq=0.50 0.25 0.50, f=[47.00 8.00 18.00] DEAL::, xq=0.25 0.50 0.50, f=[36.00 32.00 30.00] DEAL::, xq=0.50 0.50 0.50, f=[36.00 8.00 20.00] DEAL:: +DEAL::Degree 1: +DEAL::2d, 4 CELLS +DEAL:: CELL 1.0 +DEAL::, xq=0.00 0.00, f=[-9.32 -3.46] +DEAL::, xq=0.50 0.00, f=[25.32 -6.39] +DEAL::, xq=0.00 0.50, f=[-12.25 3.46] +DEAL::, xq=0.50 0.50, f=[36.25 14.39] +DEAL:: +DEAL:: CELL 1.1 +DEAL::, xq=0.50 0.00, f=[-23.96 -6.39] +DEAL::, xq=1.00 0.00, f=[79.96 -21.03] +DEAL::, xq=0.50 0.50, f=[-26.89 14.39] +DEAL::, xq=1.00 0.50, f=[90.89 69.03] +DEAL:: +DEAL:: CELL 1.2 +DEAL::, xq=0.00 0.50, f=[-12.25 -35.67] +DEAL::, xq=0.50 0.50, f=[36.25 -38.60] +DEAL::, xq=0.00 1.00, f=[-41.53 123.67] +DEAL::, xq=0.50 1.00, f=[145.53 134.60] +DEAL:: +DEAL:: CELL 1.3 +DEAL::, xq=0.50 0.50, f=[-26.89 -38.60] +DEAL::, xq=1.00 0.50, f=[90.89 -50.32] +DEAL::, xq=0.50 1.00, f=[-53.24 134.60] +DEAL::, xq=1.00 1.00, f=[189.24 178.32] +DEAL:: +DEAL::2d, 7 CELLS +DEAL:: CELL 1.1 +DEAL::, xq=0.50 0.00, f=[-9.32 -3.46] +DEAL::, xq=1.00 0.00, f=[25.32 -6.39] +DEAL::, xq=0.50 0.50, f=[-12.25 3.46] +DEAL::, xq=1.00 0.50, f=[36.25 14.39] +DEAL:: +DEAL:: CELL 1.2 +DEAL::, xq=0.00 0.50, f=[-26.89 -21.03] +DEAL::, xq=0.50 0.50, f=[90.89 -23.96] +DEAL::, xq=0.00 1.00, f=[-29.82 69.03] +DEAL::, xq=0.50 1.00, f=[101.82 79.96] +DEAL:: +DEAL:: CELL 1.3 +DEAL::, xq=0.50 0.50, f=[-12.25 -23.96] +DEAL::, xq=1.00 0.50, f=[36.25 -38.60] +DEAL::, xq=0.50 1.00, f=[-41.53 79.96] +DEAL::, xq=1.00 1.00, f=[145.53 134.60] +DEAL:: +DEAL:: CELL 2.0 +DEAL::, xq=0.00 0.00, f=[-112.34 -100.63] +DEAL::, xq=0.25 0.00, f=[400.34 -106.49] +DEAL::, xq=0.00 0.25, f=[-118.20 356.63] +DEAL::, xq=0.25 0.25, f=[422.20 378.49] +DEAL:: +DEAL:: CELL 2.1 +DEAL::, xq=0.25 0.00, f=[-141.63 -106.49] +DEAL::, xq=0.50 0.00, f=[509.63 -5.20] +DEAL::, xq=0.25 0.25, f=[-147.48 378.49] +DEAL::, xq=0.50 0.25, f=[531.48 1.73] +DEAL:: +DEAL:: CELL 2.2 +DEAL::, xq=0.00 0.25, f=[-118.20 -165.05] +DEAL::, xq=0.25 0.25, f=[422.20 -170.91] +DEAL::, xq=0.00 0.50, f=[-56.33 597.05] +DEAL::, xq=0.25 0.50, f=[61.44 618.91] +DEAL:: +DEAL:: CELL 2.3 +DEAL::, xq=0.25 0.25, f=[-147.48 -170.91] +DEAL::, xq=0.50 0.25, f=[531.48 -1.73] +DEAL::, xq=0.25 0.50, f=[2.56 618.91] +DEAL::, xq=0.50 0.50, f=[120.33 5.20] +DEAL:: +DEAL::3d, 8 CELLS +DEAL:: CELL 1.0 +DEAL::, xq=0.00 0.00 0.00, f=[-9.32 -3.46 -26.89] +DEAL::, xq=0.50 0.00 0.00, f=[25.32 -6.39 -29.82] +DEAL::, xq=0.00 0.50 0.00, f=[-12.25 3.46 -32.75] +DEAL::, xq=0.50 0.50 0.00, f=[36.25 14.39 -35.67] +DEAL::, xq=0.00 0.00 0.50, f=[-21.03 -15.18 90.89] +DEAL::, xq=0.50 0.00 0.50, f=[69.03 -18.11 101.82] +DEAL::, xq=0.00 0.50 0.50, f=[-23.96 47.18 112.75] +DEAL::, xq=0.50 0.50 0.50, f=[79.96 58.11 123.67] +DEAL:: +DEAL:: CELL 1.1 +DEAL::, xq=0.50 0.00 0.00, f=[-85.45 -6.39 -29.82] +DEAL::, xq=1.00 0.00 0.00, f=[309.45 -82.53 -100.09] +DEAL::, xq=0.50 0.50 0.00, f=[-88.38 14.39 -35.67] +DEAL::, xq=1.00 0.50 0.00, f=[320.38 298.53 -103.02] +DEAL::, xq=0.50 0.00 0.50, f=[-94.24 -18.11 101.82] +DEAL::, xq=1.00 0.00 0.50, f=[342.24 -91.31 364.09] +DEAL::, xq=0.50 0.50 0.50, f=[-97.17 58.11 123.67] +DEAL::, xq=1.00 0.50 0.50, f=[353.17 331.31 375.02] +DEAL:: +DEAL:: CELL 1.2 +DEAL::, xq=0.00 0.50 0.00, f=[-12.25 -144.02 -32.75] +DEAL::, xq=0.50 0.50 0.00, f=[36.25 -146.95 -35.67] +DEAL::, xq=0.00 1.00 0.00, f=[-149.87 528.02 -161.59] +DEAL::, xq=0.50 1.00 0.00, f=[549.87 538.95 -164.52] +DEAL::, xq=0.00 0.50 0.50, f=[-23.96 -152.80 112.75] +DEAL::, xq=0.50 0.50 0.50, f=[79.96 -155.73 123.67] +DEAL::, xq=0.00 1.00 0.50, f=[-158.66 560.80 593.59] +DEAL::, xq=0.50 1.00 0.50, f=[582.66 571.73 604.52] +DEAL:: +DEAL:: CELL 1.3 +DEAL::, xq=0.50 0.50 0.00, f=[-88.38 -146.95 -35.67] +DEAL::, xq=1.00 0.50 0.00, f=[320.38 -205.51 -103.02] +DEAL::, xq=0.50 1.00 0.00, f=[-208.44 538.95 -164.52] +DEAL::, xq=1.00 1.00 0.00, f=[768.44 757.51 -217.22] +DEAL::, xq=0.50 0.50 0.50, f=[-97.17 -155.73 123.67] +DEAL::, xq=1.00 0.50 0.50, f=[353.17 -211.37 375.02] +DEAL::, xq=0.50 1.00 0.50, f=[-214.29 571.73 604.52] +DEAL::, xq=1.00 1.00 0.50, f=[790.29 779.37 801.22] +DEAL:: +DEAL:: CELL 1.4 +DEAL::, xq=0.00 0.00 0.50, f=[-21.03 -15.18 -264.07] +DEAL::, xq=0.50 0.00 0.50, f=[69.03 -18.11 -267.00] +DEAL::, xq=0.00 0.50 0.50, f=[-23.96 47.18 -269.93] +DEAL::, xq=0.50 0.50 0.50, f=[79.96 58.11 -272.86] +DEAL::, xq=0.00 0.00 1.00, f=[-258.22 -252.36 976.07] +DEAL::, xq=0.50 0.00 1.00, f=[954.22 -255.29 987.00] +DEAL::, xq=0.00 0.50 1.00, f=[-261.15 932.36 997.93] +DEAL::, xq=0.50 0.50 1.00, f=[965.15 943.29 1008.86] +DEAL:: +DEAL:: CELL 1.5 +DEAL::, xq=0.50 0.00 0.50, f=[-94.24 -18.11 -267.00] +DEAL::, xq=1.00 0.00 0.50, f=[342.24 -91.31 -322.64] +DEAL::, xq=0.50 0.50 0.50, f=[-97.17 58.11 -272.86] +DEAL::, xq=1.00 0.50 0.50, f=[353.17 331.31 -325.57] +DEAL::, xq=0.50 0.00 1.00, f=[-316.78 -255.29 987.00] +DEAL::, xq=1.00 0.00 1.00, f=[1172.78 -313.85 1194.64] +DEAL::, xq=0.50 0.50 1.00, f=[-319.71 943.29 1008.86] +DEAL::, xq=1.00 0.50 1.00, f=[1183.71 1161.85 1205.57] +DEAL:: +DEAL:: CELL 1.6 +DEAL::, xq=0.00 0.50 0.50, f=[-23.96 -152.80 -269.93] +DEAL::, xq=0.50 0.50 0.50, f=[79.96 -155.73 -272.86] +DEAL::, xq=0.00 1.00 0.50, f=[-158.66 560.80 -369.49] +DEAL::, xq=0.50 1.00 0.50, f=[582.66 571.73 -372.42] +DEAL::, xq=0.00 0.50 1.00, f=[-261.15 -360.70 997.93] +DEAL::, xq=0.50 0.50 1.00, f=[965.15 -363.63 1008.86] +DEAL::, xq=0.00 1.00 1.00, f=[-366.56 1336.70 1369.49] +DEAL::, xq=0.50 1.00 1.00, f=[1358.56 1347.63 1380.42] +DEAL:: +DEAL:: CELL 1.7 +DEAL::, xq=0.50 0.50 0.50, f=[-97.17 -155.73 -272.86] +DEAL::, xq=1.00 0.50 0.50, f=[353.17 -211.37 -325.57] +DEAL::, xq=0.50 1.00 0.50, f=[-214.29 571.73 -372.42] +DEAL::, xq=1.00 1.00 0.50, f=[790.29 779.37 -413.41] +DEAL::, xq=0.50 0.50 1.00, f=[-319.71 -363.63 1008.86] +DEAL::, xq=1.00 0.50 1.00, f=[1183.71 -407.56 1205.57] +DEAL::, xq=0.50 1.00 1.00, f=[-410.48 1347.63 1380.42] +DEAL::, xq=1.00 1.00 1.00, f=[1522.48 1511.56 1533.41] +DEAL:: +DEAL::3d, 15 CELLS +DEAL:: CELL 1.1 +DEAL::, xq=0.50 0.00 0.00, f=[-9.32 -3.46 -26.89] +DEAL::, xq=1.00 0.00 0.00, f=[25.32 -6.39 -29.82] +DEAL::, xq=0.50 0.50 0.00, f=[-12.25 3.46 -32.75] +DEAL::, xq=1.00 0.50 0.00, f=[36.25 14.39 -35.67] +DEAL::, xq=0.50 0.00 0.50, f=[-21.03 -15.18 90.89] +DEAL::, xq=1.00 0.00 0.50, f=[69.03 -18.11 101.82] +DEAL::, xq=0.50 0.50 0.50, f=[-23.96 47.18 112.75] +DEAL::, xq=1.00 0.50 0.50, f=[79.96 58.11 123.67] +DEAL:: +DEAL:: CELL 1.2 +DEAL::, xq=0.00 0.50 0.00, f=[-88.38 -82.53 -105.95] +DEAL::, xq=0.50 0.50 0.00, f=[320.38 -85.45 -32.75] +DEAL::, xq=0.00 1.00 0.00, f=[-91.31 298.53 -108.88] +DEAL::, xq=0.50 1.00 0.00, f=[331.31 309.45 -111.81] +DEAL::, xq=0.00 0.50 0.50, f=[-100.09 -94.24 385.95] +DEAL::, xq=0.50 0.50 0.50, f=[364.09 -97.17 112.75] +DEAL::, xq=0.00 1.00 0.50, f=[-103.02 342.24 396.88] +DEAL::, xq=0.50 1.00 0.50, f=[375.02 353.17 407.81] +DEAL:: +DEAL:: CELL 1.3 +DEAL::, xq=0.50 0.50 0.00, f=[-12.25 -85.45 -32.75] +DEAL::, xq=1.00 0.50 0.00, f=[36.25 -158.66 -35.67] +DEAL::, xq=0.50 1.00 0.00, f=[-161.59 309.45 -111.81] +DEAL::, xq=1.00 1.00 0.00, f=[593.59 582.66 -170.37] +DEAL::, xq=0.50 0.50 0.50, f=[-23.96 -97.17 112.75] +DEAL::, xq=1.00 0.50 0.50, f=[79.96 -164.52 123.67] +DEAL::, xq=0.50 1.00 0.50, f=[-167.44 353.17 407.81] +DEAL::, xq=1.00 1.00 0.50, f=[615.44 604.52 626.37] +DEAL:: +DEAL:: CELL 1.4 +DEAL::, xq=0.00 0.00 0.50, f=[-208.44 -205.51 -223.08] +DEAL::, xq=0.50 0.00 0.50, f=[768.44 -15.18 -226.01] +DEAL::, xq=0.00 0.50 0.50, f=[-100.09 757.51 -228.94] +DEAL::, xq=0.50 0.50 0.50, f=[364.09 47.18 -231.86] +DEAL::, xq=0.00 0.00 1.00, f=[-217.22 -211.37 823.08] +DEAL::, xq=0.50 0.00 1.00, f=[801.22 -214.29 834.01] +DEAL::, xq=0.00 0.50 1.00, f=[-220.15 779.37 844.94] +DEAL::, xq=0.50 0.50 1.00, f=[812.15 790.29 855.86] +DEAL:: +DEAL:: CELL 1.5 +DEAL::, xq=0.50 0.00 0.50, f=[-21.03 -15.18 -226.01] +DEAL::, xq=1.00 0.00 0.50, f=[69.03 -18.11 -287.50] +DEAL::, xq=0.50 0.50 0.50, f=[-23.96 47.18 -231.86] +DEAL::, xq=1.00 0.50 0.50, f=[79.96 58.11 -290.43] +DEAL::, xq=0.50 0.00 1.00, f=[-281.64 -214.29 834.01] +DEAL::, xq=1.00 0.00 1.00, f=[1041.64 -278.72 1063.50] +DEAL::, xq=0.50 0.50 1.00, f=[-284.57 790.29 855.86] +DEAL::, xq=1.00 0.50 1.00, f=[1052.57 1030.72 1074.43] +DEAL:: +DEAL:: CELL 1.6 +DEAL::, xq=0.00 0.50 0.50, f=[-100.09 -94.24 -228.94] +DEAL::, xq=0.50 0.50 0.50, f=[364.09 -97.17 -231.86] +DEAL::, xq=0.00 1.00 0.50, f=[-103.02 342.24 -334.35] +DEAL::, xq=0.50 1.00 0.50, f=[375.02 353.17 -337.28] +DEAL::, xq=0.00 0.50 1.00, f=[-220.15 -325.57 844.94] +DEAL::, xq=0.50 0.50 1.00, f=[812.15 -328.49 855.86] +DEAL::, xq=0.00 1.00 1.00, f=[-331.42 1205.57 1238.35] +DEAL::, xq=0.50 1.00 1.00, f=[1227.42 1216.49 1249.28] +DEAL:: +DEAL:: CELL 1.7 +DEAL::, xq=0.50 0.50 0.50, f=[-23.96 -97.17 -231.86] +DEAL::, xq=1.00 0.50 0.50, f=[79.96 -164.52 -290.43] +DEAL::, xq=0.50 1.00 0.50, f=[-167.44 353.17 -337.28] +DEAL::, xq=1.00 1.00 0.50, f=[615.44 604.52 -378.27] +DEAL::, xq=0.50 0.50 1.00, f=[-284.57 -328.49 855.86] +DEAL::, xq=1.00 0.50 1.00, f=[1052.57 -372.42 1074.43] +DEAL::, xq=0.50 1.00 1.00, f=[-375.35 1216.49 1249.28] +DEAL::, xq=1.00 1.00 1.00, f=[1391.35 1380.42 1402.27] +DEAL:: +DEAL:: CELL 2.0 +DEAL::, xq=0.00 0.00 0.00, f=[-826.83 -815.11 -861.96] +DEAL::, xq=0.25 0.00 0.00, f=[3066.83 -820.97 -867.82] +DEAL::, xq=0.00 0.25 0.00, f=[-832.68 3023.11 -873.68] +DEAL::, xq=0.25 0.25 0.00, f=[3088.68 3044.97 -879.53] +DEAL::, xq=0.00 0.00 0.25, f=[-850.25 -838.54 3197.96] +DEAL::, xq=0.25 0.00 0.25, f=[3154.25 -844.39 3219.82] +DEAL::, xq=0.00 0.25 0.25, f=[-856.11 3110.54 3241.68] +DEAL::, xq=0.25 0.25 0.25, f=[3176.11 3132.39 3263.53] +DEAL:: +DEAL:: CELL 2.1 +DEAL::, xq=0.25 0.00 0.00, f=[-979.09 -820.97 -867.82] +DEAL::, xq=0.50 0.00 0.00, f=[3635.09 -31.84 -115.22] +DEAL::, xq=0.25 0.25 0.00, f=[-984.95 3044.97 -879.53] +DEAL::, xq=0.50 0.25 0.00, f=[3656.95 28.38 -49.81] +DEAL::, xq=0.25 0.00 0.25, f=[-996.66 -844.39 3219.82] +DEAL::, xq=0.50 0.00 0.25, f=[3700.66 34.34 120.33] +DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 3132.39 3263.53] +DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -39.73 51.46] +DEAL:: +DEAL:: CELL 2.2 +DEAL::, xq=0.00 0.25 0.00, f=[-832.68 -1096.22 -873.68] +DEAL::, xq=0.25 0.25 0.00, f=[3088.68 -1102.08 -879.53] +DEAL::, xq=0.00 0.50 0.00, f=[-394.95 4072.22 -83.31] +DEAL::, xq=0.25 0.50 0.00, f=[422.57 4094.08 134.29] +DEAL::, xq=0.00 0.25 0.25, f=[-856.11 -1113.79 3241.68] +DEAL::, xq=0.25 0.25 0.25, f=[3176.11 -1119.65 3263.53] +DEAL::, xq=0.00 0.50 0.25, f=[-169.82 4137.79 117.36] +DEAL::, xq=0.25 0.50 0.25, f=[181.13 4159.65 -131.63] +DEAL:: +DEAL:: CELL 2.3 +DEAL::, xq=0.25 0.25 0.00, f=[-984.95 -1102.08 -879.53] +DEAL::, xq=0.50 0.25 0.00, f=[3656.95 43.56 -49.81] +DEAL::, xq=0.25 0.50 0.00, f=[-190.57 4094.08 134.29] +DEAL::, xq=0.50 0.50 0.00, f=[626.95 -40.09 -141.87] +DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 -1119.65 3263.53] +DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -75.91 51.46] +DEAL::, xq=0.25 0.50 0.25, f=[-82.09 4159.65 -131.63] +DEAL::, xq=0.50 0.50 0.25, f=[268.87 30.16 149.12] +DEAL:: +DEAL:: CELL 2.4 +DEAL::, xq=0.00 0.00 0.25, f=[-850.25 -838.54 -1336.33] +DEAL::, xq=0.25 0.00 0.25, f=[3154.25 -844.39 -1342.19] +DEAL::, xq=0.00 0.25 0.25, f=[-856.11 3110.54 -1348.05] +DEAL::, xq=0.25 0.25 0.25, f=[3176.11 3132.39 -1353.90] +DEAL::, xq=0.00 0.00 0.50, f=[-380.71 -927.78 4968.33] +DEAL::, xq=0.25 0.00 0.50, f=[452.27 16.47 4990.19] +DEAL::, xq=0.00 0.25 0.50, f=[571.27 998.27 5012.05] +DEAL::, xq=0.25 0.25 0.50, f=[-558.41 -20.16 5033.90] +DEAL:: +DEAL:: CELL 2.5 +DEAL::, xq=0.25 0.00 0.25, f=[-996.66 -844.39 -1342.19] +DEAL::, xq=0.50 0.00 0.25, f=[3700.66 34.34 -56.33] +DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 3132.39 -1353.90] +DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -39.73 -24.49] +DEAL::, xq=0.25 0.00 0.50, f=[560.05 16.47 4990.19] +DEAL::, xq=0.50 0.00 0.50, f=[488.39 104.59 179.22] +DEAL::, xq=0.25 0.25 0.50, f=[59.35 -20.16 5033.90] +DEAL::, xq=0.50 0.25 0.50, f=[149.04 -103.77 76.78] +DEAL:: +DEAL:: CELL 2.6 +DEAL::, xq=0.00 0.25 0.25, f=[-856.11 -1113.79 -1348.05] +DEAL::, xq=0.25 0.25 0.25, f=[3176.11 -1119.65 -1353.90] +DEAL::, xq=0.00 0.50 0.25, f=[-169.82 4137.79 386.73] +DEAL::, xq=0.25 0.50 0.25, f=[181.13 4159.65 -131.78] +DEAL::, xq=0.00 0.25 0.50, f=[571.27 -446.27 5012.05] +DEAL::, xq=0.25 0.25 0.50, f=[-558.41 7.31 5033.90] +DEAL::, xq=0.00 0.50 0.50, f=[-448.24 1479.78 139.22] +DEAL::, xq=0.25 0.50 0.50, f=[480.14 -29.31 198.13] +DEAL:: +DEAL:: CELL 2.7 +DEAL::, xq=0.25 0.25 0.25, f=[-1002.52 -1119.65 -1353.90] +DEAL::, xq=0.50 0.25 0.25, f=[3722.52 -75.91 -24.49] +DEAL::, xq=0.25 0.50 0.25, f=[-82.09 4159.65 -131.78] +DEAL::, xq=0.50 0.50 0.25, f=[268.87 30.16 76.37] +DEAL::, xq=0.25 0.25 0.50, f=[59.35 7.31 5033.90] +DEAL::, xq=0.50 0.25 0.50, f=[149.04 31.59 76.78] +DEAL::, xq=0.25 0.50 0.50, f=[248.05 -29.31 198.13] +DEAL::, xq=0.50 0.50 0.50, f=[248.05 31.59 76.37] +DEAL::