From c73b86d10c6eb03267f7e7c1c60d4827d2307590 Mon Sep 17 00:00:00 2001 From: buerg Date: Mon, 13 Sep 2010 22:12:11 +0000 Subject: [PATCH] Test fe/check_derivatives adapted. Style changes. git-svn-id: https://svn.dealii.org/trunk@21955 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/source/polynomials_nedelec.cc | 1720 +++++++++++++++----- tests/fe/check_derivatives.cc | 2 +- tests/fe/check_derivatives/cmp/generic | 6 +- 3 files changed, 1292 insertions(+), 436 deletions(-) diff --git a/deal.II/base/source/polynomials_nedelec.cc b/deal.II/base/source/polynomials_nedelec.cc index de27924579..b95b251447 100644 --- a/deal.II/base/source/polynomials_nedelec.cc +++ b/deal.II/base/source/polynomials_nedelec.cc @@ -1,16 +1,3 @@ -//--------------------------------------------------------------------------- -// $Id: polynomial.cc 21932 2010-09-12 09:03:07Z buerg $ -// Version: $Name$ -// -// Copyright (C) 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 @@ -22,482 +9,1349 @@ DEAL_II_NAMESPACE_OPEN template -PolynomialsNedelec::PolynomialsNedelec (const unsigned int k) : my_degree (k), polynomial_space (create_polynomials (k)), n_pols (compute_n_pols (k)) { +PolynomialsNedelec::PolynomialsNedelec (const unsigned int k) : + my_degree (k), polynomial_space (create_polynomials (k)), + n_pols (compute_n_pols (k)) +{ } template -std::vector > > PolynomialsNedelec::create_polynomials (const unsigned int k) { - std::vector > > pols (dim); +std::vector > > +PolynomialsNedelec::create_polynomials (const unsigned int k) +{ + std::vector > > pols (dim); - pols[0] = Polynomials::Legendre::generate_complete_basis (k); + pols[0] = Polynomials::Legendre::generate_complete_basis (k); - for (unsigned int i = 1; i < dim; ++i) - pols[i] = Polynomials::Lobatto::generate_complete_basis (k + 1); + for (unsigned int i = 1; i < dim; ++i) + pols[i] = Polynomials::Lobatto::generate_complete_basis (k + 1); - return pols; + return pols; } + // Compute the values, gradients + // and double gradients of the + // polynomial at the given point. template -void PolynomialsNedelec::compute (const Point &unit_point, std::vector > &values, std::vector > &grads, std::vector > &grad_grads) const { - Assert(values.size () == n_pols || values.size () == 0, ExcDimensionMismatch(values.size (), n_pols)); - Assert(grads.size () == n_pols || grads.size () == 0, ExcDimensionMismatch(grads.size (), n_pols)); - Assert(grad_grads.size () == n_pols || grad_grads.size () == 0, ExcDimensionMismatch(grad_grads.size (), n_pols)); - // Declare the values, derivatives and second derivatives vectors of polynomial_space at unit_point - const unsigned int n_basis = polynomial_space.n (); - std::vector unit_point_values ((values.size () == 0) ? 0 : n_basis); - std::vector > unit_point_grads ((grads.size () == 0) ? 0 : n_basis); - std::vector > unit_point_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); - - switch (dim) { - case 1: { - polynomial_space.compute (unit_point, unit_point_values, unit_point_grads, unit_point_grad_grads); - // Assign the correct values to the corresponding shape functions. - if (values.size () > 0) +void +PolynomialsNedelec::compute (const Point &unit_point, + std::vector > &values, + std::vector > &grads, + std::vector > &grad_grads) +const +{ + Assert(values.size () == n_pols || values.size () == 0, + ExcDimensionMismatch(values.size (), n_pols)); + Assert(grads.size () == n_pols || grads.size () == 0, + ExcDimensionMismatch(grads.size (), n_pols)); + Assert(grad_grads.size () == n_pols || grad_grads.size () == 0, + ExcDimensionMismatch(grad_grads.size (), n_pols)); + + // Declare the values, derivatives + // and second derivatives vectors of + // polynomial_space at + // unit_point + const unsigned int& n_basis = polynomial_space.n (); + std::vector unit_point_values ((values.size () == 0) ? 0 : n_basis); + std::vector > + unit_point_grads ((grads.size () == 0) ? 0 : n_basis); + std::vector > + unit_point_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); + + switch (dim) + { + case 1: + { + polynomial_space.compute (unit_point, unit_point_values, + unit_point_grads, unit_point_grad_grads); + + // Assign the correct values to the + // corresponding shape functions. + if (values.size () > 0) for (unsigned int i = 0; i < unit_point_values.size (); ++i) - values[i][0] = unit_point_values[i]; + values[i][0] = unit_point_values[i]; - if (grads.size () > 0) + if (grads.size () > 0) for (unsigned int i = 0; i < unit_point_grads.size (); ++i) - grads[i][0][0] = unit_point_grads[i][0]; + grads[i][0][0] = unit_point_grads[i][0]; - if (grad_grads.size () > 0) + if (grad_grads.size () > 0) for (unsigned int i = 0; i < unit_point_grad_grads.size (); ++i) - grad_grads[i][0][0][0] = unit_point_grad_grads[i][0][0]; - - break; - } - - case 2: { - polynomial_space.compute (unit_point, unit_point_values, unit_point_grads, unit_point_grad_grads); - // Declare the values, derivatives and second derivatives vectors of polynomial_space at unit_point with coordinates - // shifted one step in positive direction - Point p; - std::vector p_values ((values.size () == 0) ? 0 : n_basis); - std::vector > p_grads ((grads.size () == 0) ? 0 : n_basis); - std::vector > p_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); - - p (0) = unit_point (1); - p (1) = unit_point (0); - polynomial_space.compute (p, p_values, p_grads, p_grad_grads); - // Assign the correct values to the corresponding shape functions. - if (values.size () > 0) { - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < 2; ++j) { - values[i + j * (my_degree + 1)][0] = 0.0; - values[i + j * (my_degree + 1)][1] = p_values[i + j * (my_degree + 1)]; - values[i + (j + 2) * (my_degree + 1)][0] = unit_point_values[i + j * (my_degree + 1)]; - values[i + (j + 2) * (my_degree + 1)][1] = 0.0; - } - - if (my_degree > 0) - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < my_degree; ++j) { - values[(i + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0] = unit_point_values[i + (j + 2) * (my_degree + 1)]; - values[(i + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1] = 0.0; - values[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][0] = 0.0; - values[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1] = p_values[i + (j + 2) * (my_degree + 1)]; - } - } - - if (grads.size () > 0) { - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < 2; ++j) { - for (unsigned int k = 0; k < dim; ++k) { - grads[i + j * (my_degree + 1)][0][k] = 0.0; - grads[i + (j + 2) * (my_degree + 1)][0][k] = unit_point_grads[i + j * (my_degree + 1)][k]; - grads[i + (j + 2) * (my_degree + 1)][1][k] = 0.0; + grad_grads[i][0][0][0] = unit_point_grad_grads[i][0][0]; + + break; + } + + case 2: + { + polynomial_space.compute (unit_point, unit_point_values, + unit_point_grads, unit_point_grad_grads); + + // Declare the values, derivatives and + // second derivatives vectors of + // polynomial_space at + // unit_point with coordinates + // shifted one step in positive direction + Point p; + + p (0) = unit_point (1); + p (1) = unit_point (0); + + std::vector p_values ((values.size () == 0) ? 0 : n_basis); + std::vector > + p_grads ((grads.size () == 0) ? 0 : n_basis); + std::vector > + p_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); + + polynomial_space.compute (p, p_values, p_grads, p_grad_grads); + + // Assign the correct values to the + // corresponding shape functions. + if (values.size () > 0) + { + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < 2; ++j) + { + values[i + j * (my_degree + 1)][0] = 0.0; + values[i + j * (my_degree + 1)][1] + = p_values[i + j * (my_degree + 1)]; + values[i + (j + 2) * (my_degree + 1)][0] + = unit_point_values[i + j * (my_degree + 1)]; + values[i + (j + 2) * (my_degree + 1)][1] = 0.0; } + + if (my_degree > 0) + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < my_degree; ++j) + { + values[(i + GeometryInfo::lines_per_cell) * my_degree + + j + GeometryInfo::lines_per_cell][0] + = unit_point_values[i + (j + 2) * (my_degree + 1)]; + values[(i + GeometryInfo::lines_per_cell) * my_degree + + j + GeometryInfo::lines_per_cell][1] = 0.0; + values[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][0] = 0.0; + values[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1] + = p_values[i + (j + 2) * (my_degree + 1)]; + } + } + + if (grads.size () > 0) + { + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int k = 0; k < dim; ++k) + { + grads[i + j * (my_degree + 1)][0][k] = 0.0; + grads[i + (j + 2) * (my_degree + 1)][0][k] + = unit_point_grads[i + j * (my_degree + 1)][k]; + grads[i + (j + 2) * (my_degree + 1)][1][k] = 0.0; + } - grads[i + j * (my_degree + 1)][1][0] = p_grads[i + j * (my_degree + 1)][1]; - grads[i + j * (my_degree + 1)][1][1] = p_grads[i + j * (my_degree + 1)][0]; - } - - if (my_degree > 0) - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < my_degree; ++j) { - for (unsigned int k = 0; k < dim; ++k) { - grads[(i + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0][k] = unit_point_grads[i + (j + 2) * (my_degree + 1)][k]; - grads[(i + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][k] = 0.0; - grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][0][k] = 0.0; - } - - grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][0] = p_grads[i + (j + 2) * (my_degree + 1)][1]; - grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][1] = p_grads[i + (j + 2) * (my_degree + 1)][0]; + grads[i + j * (my_degree + 1)][1][0] + = p_grads[i + j * (my_degree + 1)][1]; + grads[i + j * (my_degree + 1)][1][1] + = p_grads[i + j * (my_degree + 1)][0]; } - } - - if (grad_grads.size () > 0) { - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < 2; ++j) { - for (unsigned int k = 0; k < dim; ++k) - for (unsigned int l = 0; l < dim; ++l) { - grad_grads[i + j * (my_degree + 1)][0][k][l] = 0.0; - grad_grads[i + (j + 2) * (my_degree + 1)][0][k][l] = unit_point_grad_grads[i + l * (my_degree + 1)][k][l]; - grad_grads[i + (j + 2) * (my_degree + 1)][1][k][l] = 0.0; - } - - grad_grads[i + j * (my_degree + 1)][1][0][0] = p_grad_grads[i + j * (my_degree + 1)][1][1]; - grad_grads[i + j * (my_degree + 1)][1][0][1] = p_grad_grads[i + j * (my_degree + 1)][1][0]; - grad_grads[i + j * (my_degree + 1)][1][1][0] = p_grad_grads[i + j * (my_degree + 1)][0][1]; - grad_grads[i + j * (my_degree + 1)][1][1][1] = p_grad_grads[i + j * (my_degree + 1)][0][0]; - } - - if (my_degree > 0) - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < my_degree; ++j) { - for (unsigned int k = 0; k < dim; ++k) - for (unsigned int l = 0; l < dim; ++l) { - grad_grads[(i + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0][k][l] = unit_point_grad_grads[i + (j + 2) * (my_degree + 1)][k][l]; - grad_grads[(i + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][k][l] = 0.0; - grad_grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][0][k][l] = 0.0; + + if (my_degree > 0) + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < my_degree; ++j) + { + for (unsigned int k = 0; k < dim; ++k) + { + grads[(i + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][0][k] + = unit_point_grads[i + (j + 2) * (my_degree + 1)] + [k]; + grads[(i + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][k] + = 0.0; + grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][0][k] = 0.0; } - - grad_grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][0][0] = p_grad_grads[i + (j + 2) * (my_degree + 1)][1][1]; - grad_grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][0][1] = p_grad_grads[i + (j + 2) * (my_degree + 1)][1][0]; - grad_grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][1][0] = p_grad_grads[i + (j + 2) * (my_degree + 1)][0][1]; - grad_grads[i + (j + my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][1][1] = p_grad_grads[i + (j + 2) * (my_degree + 1)][0][0]; - } - } - - break; - } - - case 3: { - polynomial_space.compute (unit_point, unit_point_values, unit_point_grads, unit_point_grad_grads); - // Declare the values, derivatives and second derivatives vectors of polynomial_space at unit_point with coordinates - // shifted two steps in positive direction - Point p1, p2; - std::vector p1_values ((values.size () == 0) ? 0 : n_basis); - std::vector > p1_grads ((grads.size () == 0) ? 0 : n_basis); - std::vector > p1_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); - std::vector p2_values ((values.size () == 0) ? 0 : n_basis); - std::vector > p2_grads ((grads.size () == 0) ? 0 : n_basis); - std::vector > p2_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); - - p1 (0) = unit_point (1); - p1 (1) = unit_point (2); - p1 (2) = unit_point (0); - polynomial_space.compute (p1, p1_values, p1_grads, p1_grad_grads); - p2 (0) = unit_point (2); - p2 (1) = unit_point (0); - p2 (2) = unit_point (1); - polynomial_space.compute (p2, p2_values, p2_grads, p2_grad_grads); - // Assign the correct values to the corresponding shape functions. - if (values.size () > 0) { - for (unsigned int i = 0; i <= my_degree; ++i) { - for (unsigned int j = 0; j < 2; ++j) { - for (unsigned int k = 0; k < 2; ++k) { - for (unsigned int l = 0; l < 2; ++l) { - values[i + (j + 4 * k) * (my_degree + 1)][2 * l] = 0; - values[i + (j + 4 * k + 2) * (my_degree + 1)][l + 1] = 0; - values[i + (j + 2 * (k + 4)) * (my_degree + 1)][l] = 0; - } - - values[i + (j + 4 * k + 2) * (my_degree + 1)][0] = unit_point_values[i + (j + k * (my_degree + 2)) * (my_degree + 1)]; - values[i + (j + 2 * (k + 4)) * (my_degree + 1)][2] = p2_values[i + (j + k * (my_degree + 2)) * (my_degree + 1)]; - } + + grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][0] + = p_grads[i + (j + 2) * (my_degree + 1)][1]; + grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][1] + = p_grads[i + (j + 2) * (my_degree + 1)][0]; + } + } - values[i + j * (my_degree + 1)][1] = p1_values[i + j * (my_degree + 1) * (my_degree + 2)]; - } + if (grad_grads.size () > 0) + { + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + { + grad_grads[i + j * (my_degree + 1)][0][k][l] = 0.0; + grad_grads[i + (j + 2) * (my_degree + 1)][0][k][l] + = unit_point_grad_grads[i + l * (my_degree + 1)][k] + [l]; + grad_grads[i + (j + 2) * (my_degree + 1)][1][k][l] + = 0.0; + } + + grad_grads[i + j * (my_degree + 1)][1][0][0] + = p_grad_grads[i + j * (my_degree + 1)][1][1]; + grad_grads[i + j * (my_degree + 1)][1][0][1] + = p_grad_grads[i + j * (my_degree + 1)][1][0]; + grad_grads[i + j * (my_degree + 1)][1][1][0] + = p_grad_grads[i + j * (my_degree + 1)][0][1]; + grad_grads[i + j * (my_degree + 1)][1][1][1] + = p_grad_grads[i + j * (my_degree + 1)][0][0]; + } - values[i + 4 * (my_degree + 1)][1] = p1_values[i + my_degree + 1]; - values[i + 5 * (my_degree + 1)][1] = p1_values[i + (my_degree + 1) * (my_degree + 3)]; + if (my_degree > 0) + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < my_degree; ++j) + { + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + { + grad_grads[(i + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][0] + [k][l] + = unit_point_grad_grads[i + (j + 2) + * (my_degree + 1)][k][l]; + grad_grads[(i + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1] + [k][l] = 0.0; + grad_grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][0][k][l] = 0.0; + } + + grad_grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][0][0] + = p_grad_grads[i + (j + 2) * (my_degree + 1)][1][1]; + grad_grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][0][1] + = p_grad_grads[i + (j + 2) * (my_degree + 1)][1][0]; + grad_grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][1][0] + = p_grad_grads[i + (j + 2) * (my_degree + 1)][0][1]; + grad_grads[i + (j + my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][1][1] + = p_grad_grads[i + (j + 2) * (my_degree + 1)][0][0]; + } } - if (my_degree > 0) - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < my_degree; ++j) { - for (unsigned int k = 0; k < my_degree; ++k) { - for (unsigned int l = 0; l < 2; ++l) { - values[((i + 2 * GeometryInfo::faces_per_cell) * my_degree + j + GeometryInfo::lines_per_cell + 2 * GeometryInfo::faces_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][l + 1] = 0; - values[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][2 * l] = 0; - values[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][l] = 0; - } + break; + } - values[((i + 2 * GeometryInfo::faces_per_cell) * my_degree + j + GeometryInfo::lines_per_cell + 2 * GeometryInfo::faces_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][0] = unit_point_values[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)]; - values[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1] = p1_values[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)]; - values[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2] = p2_values[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)]; - } + case 3: + { + polynomial_space.compute (unit_point, unit_point_values, + unit_point_grads, unit_point_grad_grads); + + // Declare the values, derivatives + // and second derivatives vectors of + // polynomial_space at + // unit_point with coordinates + // shifted two steps in positive + // direction + Point p1, p2; + std::vector p1_values ((values.size () == 0) ? 0 : n_basis); + std::vector > + p1_grads ((grads.size () == 0) ? 0 : n_basis); + std::vector > + p1_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); + std::vector p2_values ((values.size () == 0) ? 0 : n_basis); + std::vector > + p2_grads ((grads.size () == 0) ? 0 : n_basis); + std::vector > + p2_grad_grads ((grad_grads.size () == 0) ? 0 : n_basis); - for (unsigned int k = 0; k < 2; ++k) { - for (unsigned int l = 0; l < 2; ++l) { - for (unsigned int m = 0; m < 2; ++m) { - values[i + (j + (2 * (k + 2 * l) + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][m] = 0; - values[(i + 2 * (k + 2 * (l + 1)) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][m + 1] = 0; - } + p1 (0) = unit_point (1); + p1 (1) = unit_point (2); + p1 (2) = unit_point (0); + polynomial_space.compute (p1, p1_values, p1_grads, p1_grad_grads); + p2 (0) = unit_point (2); + p2 (1) = unit_point (0); + p2 (2) = unit_point (1); + polynomial_space.compute (p2, p2_values, p2_grads, p2_grad_grads); + + // Assign the correct values to the + // corresponding shape functions. + if (values.size () > 0) + { + for (unsigned int i = 0; i <= my_degree; ++i) + { + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int k = 0; k < 2; ++k) + { + for (unsigned int l = 0; l < 2; ++l) + { + values[i + (j + 4 * k) * (my_degree + 1)][2 * l] + = 0.0; + values[i + (j + 4 * k + 2) * (my_degree + 1)] + [l + 1] = 0.0; + values[i + (j + 2 * (k + 4)) * (my_degree + 1)][l] + = 0.0; + } - values[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][2 * l] = 0; - values[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2 * l] = 0; + values[i + (j + 4 * k + 2) * (my_degree + 1)][0] + = unit_point_values[i + (j + k * (my_degree + 2)) + * (my_degree + 1)]; + values[i + (j + 2 * (k + 4)) * (my_degree + 1)][2] + = p2_values[i + (j + k * (my_degree + 2)) + * (my_degree + 1)]; } - values[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1] = p1_values[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)]; - values[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2] = p2_values[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)]; - values[(i + 2 * (k + 2) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0] = unit_point_values[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)]; - values[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2] = p2_values[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)]; - values[(i + 2 * (k + 4) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0] = unit_point_values[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)]; - values[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1] = p1_values[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)]; - } - } - } - - if (grads.size () > 0) { - for (unsigned int i = 0; i <= my_degree; ++i) { - for (unsigned int j = 0; j < 2; ++j) { - for (unsigned int k = 0; k < 2; ++k) { - for (unsigned int l = 0; l < 2; ++l) - for (unsigned int m = 0; m < dim; ++m) { - grads[i + (j + 4 * k) * (my_degree + 1)][2 * l][m] = 0; - grads[i + (j + 4 * k + 2) * (my_degree + 1)][l + 1][m] = 0; - grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][l][m] = 0; + values[i + j * (my_degree + 1)][1] + = p1_values[i + j * (my_degree + 1) * (my_degree + 2)]; + } + + values[i + 4 * (my_degree + 1)][1] + = p1_values[i + my_degree + 1]; + values[i + 5 * (my_degree + 1)][1] + = p1_values[i + (my_degree + 1) * (my_degree + 3)]; + } + + if (my_degree > 0) + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < my_degree; ++j) + { + for (unsigned int k = 0; k < my_degree; ++k) + { + for (unsigned int l = 0; l < 2; ++l) + { + values[((i + 2 + * GeometryInfo::faces_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell + 2 + * GeometryInfo::faces_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][l + 1] + = 0.0; + values[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) + * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][2 * l] + = 0.0; + values[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][l] = 0.0; + } + + values[((i + 2 * GeometryInfo::faces_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell + 2 + * GeometryInfo::faces_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][0] + = unit_point_values[i + (j + (k + 2) * (my_degree + + 2) + 2) + * (my_degree + 1)]; + values[(i + (j + 2 * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1] + = p1_values[i + ((j + 2) * (my_degree + 2) + k + 2) + * (my_degree + 1)]; + values[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2] + = p2_values[i + (j + (k + 2) * (my_degree + 2) + 2) + * (my_degree + 1)]; } - for (unsigned int l = 0; l < dim; ++l) - grads[i + (j + 4 * k + 2) * (my_degree + 1)][0][l] = unit_point_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][l]; - - grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][0] = p2_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][1]; - grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][1] = p2_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][2]; - grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][2] = p2_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][0]; - } + for (unsigned int k = 0; k < 2; ++k) + { + for (unsigned int l = 0; l < 2; ++l) + { + for (unsigned int m = 0; m < 2; ++m) + { + values[i + (j + (2 * (k + 2 * l) + 1) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][m] = 0.0; + values[(i + 2 * (k + 2 * (l + 1)) * (my_degree + + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell] + [m + 1] = 0.0; + } + + values[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][2 * l] + = 0.0; + values[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2 * l] = 0.0; + } - grads[i + j * (my_degree + 1)][1][0] = p1_grads[i + j * (my_degree + 1) * (my_degree + 2)][2]; - grads[i + j * (my_degree + 1)][1][1] = p1_grads[i + j * (my_degree + 1) * (my_degree + 2)][0]; - grads[i + j * (my_degree + 1)][1][2] = p1_grads[i + j * (my_degree + 1) * (my_degree + 2)][1]; - } - - grads[i + 4 * (my_degree + 1)][1][0] = p1_grads[i + my_degree + 1][2]; - grads[i + 4 * (my_degree + 1)][1][1] = p1_grads[i + my_degree + 1][0]; - grads[i + 4 * (my_degree + 1)][1][2] = p1_grads[i + my_degree + 1][1]; - grads[i + 5 * (my_degree + 1)][1][0] = p1_grads[i + (my_degree + 1) * (my_degree + 3)][2]; - grads[i + 5 * (my_degree + 1)][1][1] = p1_grads[i + (my_degree + 1) * (my_degree + 3)][0]; - grads[i + 5 * (my_degree + 1)][1][2] = p1_grads[i + (my_degree + 1) * (my_degree + 3)][1]; + values[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1] + = p1_values[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)]; + values[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2] + = p2_values[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)]; + values[(i + 2 * (k + 2) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][0] + = unit_point_values[i + ((j + 2) * (my_degree + 2) + + k) * (my_degree + 1)]; + values[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2] + = p2_values[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)]; + values[(i + 2 * (k + 4) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][0] + = unit_point_values[i + (j + k * (my_degree + 2) + + 2) * (my_degree + 1)]; + values[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1] + = p1_values[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)]; + } + } } - if (my_degree > 0) - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < my_degree; ++j) { - for (unsigned int k = 0; k < my_degree; ++k) { - for (unsigned int l = 0; l < dim; ++l) { - for (unsigned int m = 0; m < 2; ++m) { - grads[((i + 2 * GeometryInfo::faces_per_cell) * my_degree + j + GeometryInfo::lines_per_cell + 2 * GeometryInfo::faces_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][m + 1][l] = 0; - grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][2 * m][l] = 0; - grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][m][l] = 0; - } - - grads[((i + 2 * GeometryInfo::faces_per_cell) * my_degree + j + GeometryInfo::lines_per_cell + 2 * GeometryInfo::faces_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][0][l] = unit_point_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][l]; + if (grads.size () > 0) + { + for (unsigned int i = 0; i <= my_degree; ++i) + { + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int k = 0; k < 2; ++k) + { + for (unsigned int l = 0; l < 2; ++l) + for (unsigned int m = 0; m < dim; ++m) + { + grads[i + (j + 4 * k) * (my_degree + 1)][2 * l] + [m] = 0.0; + grads[i + (j + 4 * k + 2) * (my_degree + 1)] + [l + 1][m] = 0.0; + grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [l][m] = 0.0; + } + + for (unsigned int l = 0; l < dim; ++l) + grads[i + (j + 4 * k + 2) * (my_degree + 1)][0][l] + = unit_point_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][l]; + + grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][0] + = p2_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][1]; + grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][1] + = p2_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][2]; + grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][2] + = p2_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][0]; } + + grads[i + j * (my_degree + 1)][1][0] + = p1_grads[i + j * (my_degree + 1) * (my_degree + 2)] + [2]; + grads[i + j * (my_degree + 1)][1][1] + = p1_grads[i + j * (my_degree + 1) * (my_degree + 2)] + [0]; + grads[i + j * (my_degree + 1)][1][2] + = p1_grads[i + j * (my_degree + 1) * (my_degree + 2)] + [1]; + } + + grads[i + 4 * (my_degree + 1)][1][0] + = p1_grads[i + my_degree + 1][2]; + grads[i + 4 * (my_degree + 1)][1][1] + = p1_grads[i + my_degree + 1][0]; + grads[i + 4 * (my_degree + 1)][1][2] + = p1_grads[i + my_degree + 1][1]; + grads[i + 5 * (my_degree + 1)][1][0] + = p1_grads[i + (my_degree + 1) * (my_degree + 3)][2]; + grads[i + 5 * (my_degree + 1)][1][1] + = p1_grads[i + (my_degree + 1) * (my_degree + 3)][0]; + grads[i + 5 * (my_degree + 1)][1][2] + = p1_grads[i + (my_degree + 1) * (my_degree + 3)][1]; + } + + if (my_degree > 0) + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < my_degree; ++j) + { + for (unsigned int k = 0; k < my_degree; ++k) + { + for (unsigned int l = 0; l < dim; ++l) + { + for (unsigned int m = 0; m < 2; ++m) + { + grads[((i + 2 + * GeometryInfo::faces_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell + + 2 + * GeometryInfo::faces_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell] + [m + 1][l] = 0.0; + grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell] + [2 * m][l] = 0.0; + grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][m][l] = 0.0; + } + + grads[((i + 2 * GeometryInfo::faces_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell + 2 + * GeometryInfo::faces_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][0][l] + = unit_point_grads[i + (j + (k + 2) + * (my_degree + 2) + 2) + * (my_degree + 1)][l]; + } - grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][0] = p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][2]; - grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][1] = p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][0]; - grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][2] = p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][1]; - grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0] = p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][1]; - grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1] = p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][2]; - grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2] = p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][0]; - } - - for (unsigned int k = 0; k < 2; ++k) { - for (unsigned int l = 0; l < 2; ++l) - for (unsigned int m = 0; m < dim; ++m) { - for (unsigned int n = 0; n < 2; ++n) { - grads[i + (j + (2 * (k + 2 * l) + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][n][m] = 0; - grads[(i + 2 * (k + 2 * (l + 1)) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][n + 1][m] = 0; + grads[(i + (j + 2 * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][0] + = p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) + * (my_degree + 1)][2]; + grads[(i + (j + 2 * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][1] + = p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) + * (my_degree + 1)][0]; + grads[(i + (j + 2 * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][2] + = p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) + * (my_degree + 1)][1]; + grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0] + = p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) + * (my_degree + 1)][1]; + grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1] + = p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) + * (my_degree + 1)][2]; + grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2] + = p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) + * (my_degree + 1)][0]; + } + + for (unsigned int k = 0; k < 2; ++k) + { + for (unsigned int l = 0; l < 2; ++l) + for (unsigned int m = 0; m < dim; ++m) + { + for (unsigned int n = 0; n < 2; ++n) + { + grads[i + (j + (2 * (k + 2 * l) + 1) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][n][m] = 0.0; + grads[(i + 2 * (k + 2 * (l + 1)) + * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell] + [n + 1][m] = 0.0; + } + + grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell] + [2 * l][m] = 0.0; + grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2 * l][m] = 0.0; } - grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][2 * l][m] = 0; - grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2 * l][m] = 0; - } + for (unsigned int l = 0; l < dim; ++l) + { + grads[(i + 2 * (k + 2) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][0][l] + = unit_point_grads[i + ((j + 2) * (my_degree + + 2) + k) + * (my_degree + 1)][l]; + grads[(i + 2 * (k + 4) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][0][l] + = unit_point_grads[i + (j + k * (my_degree + 2) + + 2) * (my_degree + 1)] + [l]; + } - for (unsigned int l = 0; l < dim; ++l) { - grads[(i + 2 * (k + 2) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0][l] = unit_point_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][l]; - grads[(i + 2 * (k + 4) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0][l] = unit_point_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][l]; + grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][0] + = p1_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2]; + grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][1] + = p1_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0]; + grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][2] + = p1_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1]; + grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0] + = p2_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1]; + grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1] + = p2_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2]; + grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2] + = p2_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0]; + grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0] + = p2_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1]; + grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1] + = p2_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2]; + grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2] + = p2_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0]; + grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][0] + = p1_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2]; + grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][1] + = p1_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0]; + grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][2] + = p1_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1]; } + } + } - grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][0] = p1_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2]; - grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][1] = p1_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0]; - grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][2] = p1_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1]; - grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0] = p2_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1]; - grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1] = p2_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2]; - grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2] = p2_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0]; - grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0] = p2_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1]; - grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1] = p2_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2]; - grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2] = p2_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0]; - grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][0] = p1_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1]; - grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][1] = p1_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2]; - grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][2] = p1_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0]; - } - } - } - - if (grad_grads.size () > 0) { - for (unsigned int i = 0; i <= my_degree; ++i) { - for (unsigned int j = 0; j < 2; ++j) { - for (unsigned int k = 0; k < 2; ++k) { - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) { - for (unsigned int n = 0; n < 2; ++n) { - grad_grads[i + (j + 4 * k) * (my_degree + 1)][2 * n][l][m] = 0; - grad_grads[i + (j + 4 * k + 2) * (my_degree + 1)][n + 1][l][m] = 0; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][n][l][m] = 0; - } - - grad_grads[i + (j + 4 * k + 2) * (my_degree + 1)][0][l][m] = unit_point_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][l][m]; - } + if (grad_grads.size () > 0) + { + for (unsigned int i = 0; i <= my_degree; ++i) + { + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int k = 0; k < 2; ++k) + { + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + { + for (unsigned int n = 0; n < 2; ++n) + { + grad_grads[i + (j + 4 * k) * (my_degree + + 1)][2 * n] + [l][m] = 0.0; + grad_grads[i + (j + 4 * k + 2) * (my_degree + + 1)] + [n + 1][l][m] = 0.0; + grad_grads[i + (j + 2 * (k + 4)) + * (my_degree + 1)][n][l][m] + = 0.0; + } + + grad_grads[i + (j + 4 * k + 2) * (my_degree + + 1)][0][l][m] + = unit_point_grad_grads[i + (j + k + * (my_degree + + 2)) + * (my_degree + 1)][l] + [m]; + } - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][0][0] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][1][1]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][0][1] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][1][2]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][0][2] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][1][0]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][1][0] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][2][1]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][1][1] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][2][2]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][1][2] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][2][0]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][2][0] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][0][1]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][2][1] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][0][2]; - grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][2][2] = p2_grad_grads[i + (j + k * (my_degree + 2)) * (my_degree + 1)][0][0]; - } + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][0][0] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][1][1]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][0][1] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][1][2]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][0][2] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][1][0]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][1][0] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][2][1]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][1][1] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][2][2]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][1][2] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][2][0]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][2][0] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][0][1]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][2][1] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][0][2]; + grad_grads[i + (j + 2 * (k + 4)) * (my_degree + 1)] + [2][2][2] + = p2_grad_grads[i + (j + k * (my_degree + 2)) + * (my_degree + 1)][0][0]; + } - grad_grads[i + j * (my_degree + 1)][1][0][0] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][2][2]; - grad_grads[i + j * (my_degree + 1)][1][0][1] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][2][0]; - grad_grads[i + j * (my_degree + 1)][1][0][2] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][2][1]; - grad_grads[i + j * (my_degree + 1)][1][1][0] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][0][2]; - grad_grads[i + j * (my_degree + 1)][1][1][1] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][0][0]; - grad_grads[i + j * (my_degree + 1)][1][1][2] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][0][1]; - grad_grads[i + j * (my_degree + 1)][1][2][0] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][1][2]; - grad_grads[i + j * (my_degree + 1)][1][2][1] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][1][0]; - grad_grads[i + j * (my_degree + 1)][1][2][2] = p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)][1][1]; - } - - grad_grads[i + 4 * (my_degree + 1)][1][0][0] = p1_grad_grads[i + my_degree + 1][2][2]; - grad_grads[i + 4 * (my_degree + 1)][1][0][1] = p1_grad_grads[i + my_degree + 1][2][0]; - grad_grads[i + 4 * (my_degree + 1)][1][0][2] = p1_grad_grads[i + my_degree + 1][2][1]; - grad_grads[i + 4 * (my_degree + 1)][1][1][0] = p1_grad_grads[i + my_degree + 1][0][2]; - grad_grads[i + 4 * (my_degree + 1)][1][1][1] = p1_grad_grads[i + my_degree + 1][0][0]; - grad_grads[i + 4 * (my_degree + 1)][1][1][2] = p1_grad_grads[i + my_degree + 1][0][1]; - grad_grads[i + 4 * (my_degree + 1)][1][2][0] = p1_grad_grads[i + my_degree + 1][1][2]; - grad_grads[i + 4 * (my_degree + 1)][1][2][1] = p1_grad_grads[i + my_degree + 1][1][0]; - grad_grads[i + 4 * (my_degree + 1)][1][2][2] = p1_grad_grads[i + my_degree + 1][1][1]; - grad_grads[i + 5 * (my_degree + 1)][1][0][0] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2][2]; - grad_grads[i + 5 * (my_degree + 1)][1][0][1] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2][0]; - grad_grads[i + 5 * (my_degree + 1)][1][0][2] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2][1]; - grad_grads[i + 5 * (my_degree + 1)][1][1][0] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0][2]; - grad_grads[i + 5 * (my_degree + 1)][1][1][1] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0][0]; - grad_grads[i + 5 * (my_degree + 1)][1][1][2] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0][1]; - grad_grads[i + 5 * (my_degree + 1)][1][2][0] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1][2]; - grad_grads[i + 5 * (my_degree + 1)][1][2][1] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1][0]; - grad_grads[i + 5 * (my_degree + 1)][1][2][2] = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1][1]; - } + grad_grads[i + j * (my_degree + 1)][1][0][0] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][2][2]; + grad_grads[i + j * (my_degree + 1)][1][0][1] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][2][0]; + grad_grads[i + j * (my_degree + 1)][1][0][2] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][2][1]; + grad_grads[i + j * (my_degree + 1)][1][1][0] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][0][2]; + grad_grads[i + j * (my_degree + 1)][1][1][1] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][0][0]; + grad_grads[i + j * (my_degree + 1)][1][1][2] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][0][1]; + grad_grads[i + j * (my_degree + 1)][1][2][0] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][1][2]; + grad_grads[i + j * (my_degree + 1)][1][2][1] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][1][0]; + grad_grads[i + j * (my_degree + 1)][1][2][2] + = p1_grad_grads[i + j * (my_degree + 1) + * (my_degree + 2)][1][1]; + } - if (my_degree > 0) - for (unsigned int i = 0; i <= my_degree; ++i) - for (unsigned int j = 0; j < my_degree; ++j) { - for (unsigned int k = 0; k < my_degree; ++k) { - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) { - for (unsigned int n = 0; n < 2; ++n) { - grad_grads[((i + 2 * GeometryInfo::faces_per_cell) * my_degree + j + GeometryInfo::lines_per_cell + 2 * GeometryInfo::faces_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][n + 1][l][m] = 0; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][2 * n][l][m] = 0; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][n][l][m] = 0; - } + grad_grads[i + 4 * (my_degree + 1)][1][0][0] + = p1_grad_grads[i + my_degree + 1][2][2]; + grad_grads[i + 4 * (my_degree + 1)][1][0][1] + = p1_grad_grads[i + my_degree + 1][2][0]; + grad_grads[i + 4 * (my_degree + 1)][1][0][2] + = p1_grad_grads[i + my_degree + 1][2][1]; + grad_grads[i + 4 * (my_degree + 1)][1][1][0] + = p1_grad_grads[i + my_degree + 1][0][2]; + grad_grads[i + 4 * (my_degree + 1)][1][1][1] + = p1_grad_grads[i + my_degree + 1][0][0]; + grad_grads[i + 4 * (my_degree + 1)][1][1][2] + = p1_grad_grads[i + my_degree + 1][0][1]; + grad_grads[i + 4 * (my_degree + 1)][1][2][0] + = p1_grad_grads[i + my_degree + 1][1][2]; + grad_grads[i + 4 * (my_degree + 1)][1][2][1] + = p1_grad_grads[i + my_degree + 1][1][0]; + grad_grads[i + 4 * (my_degree + 1)][1][2][2] + = p1_grad_grads[i + my_degree + 1][1][1]; + grad_grads[i + 5 * (my_degree + 1)][1][0][0] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2] + [2]; + grad_grads[i + 5 * (my_degree + 1)][1][0][1] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2] + [0]; + grad_grads[i + 5 * (my_degree + 1)][1][0][2] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2] + [1]; + grad_grads[i + 5 * (my_degree + 1)][1][1][0] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0] + [2]; + grad_grads[i + 5 * (my_degree + 1)][1][1][1] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0] + [0]; + grad_grads[i + 5 * (my_degree + 1)][1][1][2] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0] + [1]; + grad_grads[i + 5 * (my_degree + 1)][1][2][0] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1] + [2]; + grad_grads[i + 5 * (my_degree + 1)][1][2][1] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1] + [0]; + grad_grads[i + 5 * (my_degree + 1)][1][2][2] + = p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1] + [1]; + } - grad_grads[((i + 2 * GeometryInfo::faces_per_cell) * my_degree + j + GeometryInfo::lines_per_cell + 2 * GeometryInfo::faces_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][0][l][m] = unit_point_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][l][m]; - } - - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][0][0] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][2][2]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][0][1] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][2][0]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][0][2] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][2][1]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][1][0] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][0][2]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][1][1] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][0][0]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][1][2] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][0][1]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][2][0] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][1][2]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][2][1] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][1][0]; - grad_grads[(i + (j + 2 * GeometryInfo::faces_per_cell + my_degree) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + k + GeometryInfo::lines_per_cell][1][2][2] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + 2) * (my_degree + 1)][1][1]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][0] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][1][1]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][1] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][1][2]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][2] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][1][0]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][0] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][2][1]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][1] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][2][2]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][2] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][2][0]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][0] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][0][1]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][1] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][0][2]; - grad_grads[i + (j + (k + 2 * (GeometryInfo::faces_per_cell + my_degree)) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][2] = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + 2) * (my_degree + 1)][0][0]; - } - - for (unsigned int k = 0; k < 2; ++k) { - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) { - for (unsigned int n = 0; n < 2; ++n) { - for (unsigned int o = 0; o < 2; ++o) { - grad_grads[i + (j + (2 * (k + 2 * n) + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][o][l][m] = 0; - grad_grads[(i + 2 * (k + 2 * (n + 1)) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][o + 1][l][m] = 0; - } - - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][2 * n][l][m] = 0; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2 * n][l][m] = 0; + if (my_degree > 0) + for (unsigned int i = 0; i <= my_degree; ++i) + for (unsigned int j = 0; j < my_degree; ++j) + { + for (unsigned int k = 0; k < my_degree; ++k) + { + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + { + for (unsigned int n = 0; n < 2; ++n) + { + grad_grads[((i + 2 + * GeometryInfo::faces_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell + + 2 + * GeometryInfo::faces_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell] + [n + 1][l][m] = 0.0; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell] + [2 * n][l][m] = 0.0; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][n][l][m] + = 0.0; + } + + grad_grads[((i + 2 + * GeometryInfo::faces_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell + + 2 + * GeometryInfo::faces_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell] + [0][l][m] + = unit_point_grad_grads[i + (j + (k + 2) + * (my_degree + 2) + + 2) * (my_degree + + 1)][l][m]; } + + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) + * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][0] + [0] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][2][2]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][0] + [1] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][2][0]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][0] + [2] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][2][1]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][1] + [0] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][0][2]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][1] + [1] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][0][0]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][1] + [2] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][0][1]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][2] + [0] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][1][2]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][2] + [1] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][1][0]; + grad_grads[(i + (j + 2 + * GeometryInfo::faces_per_cell + + my_degree) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + k + + GeometryInfo::lines_per_cell][1][2] + [2] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k + + 2) * (my_degree + 1)][1][1]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][0] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][1][1]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][1] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][1][2]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][2] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][1][0]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][0] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][2][1]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][1] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][2][2]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][2] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][2][0]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][0] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][0][1]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][1] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][0][2]; + grad_grads[i + (j + (k + 2 + * (GeometryInfo::faces_per_cell + + my_degree)) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][2] + = p2_grad_grads[i + (j + (k + 2) * (my_degree + 2) + + 2) * (my_degree + 1)][0][0]; + } + + for (unsigned int k = 0; k < 2; ++k) + { + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + { + for (unsigned int n = 0; n < 2; ++n) + { + for (unsigned int o = 0; o < 2; ++o) + { + grad_grads[i + (j + (2 * (k + 2 * n) + + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][o][l][m] + = 0.0; + grad_grads[(i + 2 * (k + 2 * (n + 1)) + * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell] + [o + 1][l][m] = 0.0; + } + + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell] + [2 * n][l][m] = 0.0; + grad_grads[i + (j + (2 * k + 9) + * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2 * n][l][m] + = 0.0; + } - grad_grads[(i + 2 * (k + 2) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0][l][m] = unit_point_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][l][m]; - grad_grads[(i + 2 * (k + 4) * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][0][l][m] = unit_point_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][l][m]; - } - - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][0][0] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2][2]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][0][1] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2][0]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][0][2] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2][1]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][1][0] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0][2]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][1][1] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0][0]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][1][2] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0][1]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][2][0] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1][2]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][2][1] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1][0]; - grad_grads[(i + 2 * k * (my_degree + 1) + GeometryInfo::lines_per_cell) * my_degree + j + GeometryInfo::lines_per_cell][1][2][2] = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1][1]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][0] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1][1]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][1] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1][2]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][2] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1][0]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][0] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2][1]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][1] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2][2]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][2] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2][0]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][0] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0][1]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][1] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0][2]; - grad_grads[i + (j + (2 * k + 1) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][2] = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0][0]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][0] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1][1]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][1] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1][2]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][0][2] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][1][0]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][0] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2][1]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][1] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2][2]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][1][2] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][2][0]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][0] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0][1]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][1] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0][2]; - grad_grads[i + (j + (2 * k + 5) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][2][2][2] = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) * (my_degree + 1)][0][0]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][0][0] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2][2]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][0][1] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2][0]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][0][2] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][2][1]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][1][0] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0][2]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][1][1] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0][0]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][1][2] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][0][1]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][2][0] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1][2]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][2][1] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1][0]; - grad_grads[i + (j + (2 * k + 9) * my_degree + GeometryInfo::lines_per_cell) * (my_degree + 1)][1][2][2] = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) * (my_degree + 1)][1][1]; - } - } - } - } - } + grad_grads[(i + 2 * (k + 2) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell] + [0][l][m] + = unit_point_grad_grads[i + ((j + 2) + * (my_degree + + 2) + k) + * (my_degree + 1)][l] + [m]; + grad_grads[(i + 2 * (k + 4) * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell] + [0][l][m] + = unit_point_grad_grads[i + (j + k + * (my_degree + + 2) + 2) + * (my_degree + 1)][l] + [m]; + } + + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][0] + [0] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2][2]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][0] + [1] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2][0]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][0] + [2] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2][1]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][1] + [0] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0][2]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][1] + [1] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0][0]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][1] + [2] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0][1]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][2] + [0] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1][2]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][2] + [1] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1][0]; + grad_grads[(i + 2 * k * (my_degree + 1) + + GeometryInfo::lines_per_cell) + * my_degree + j + + GeometryInfo::lines_per_cell][1][2] + [2] + = p1_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1][1]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][0] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1][1]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][1] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1][2]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][2] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1][0]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][0] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2][1]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][1] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2][2]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][2] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2][0]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][0] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0][1]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][1] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0][2]; + grad_grads[i + (j + (2 * k + 1) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][2] + = p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0][0]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][0] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1][1]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][1] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1][2]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][0][2] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][1][0]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][0] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2][1]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][1] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2][2]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][1][2] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][2][0]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][0] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0][1]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][1] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0][2]; + grad_grads[i + (j + (2 * k + 5) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][2][2][2] + = p2_grad_grads[i + (j + k * (my_degree + 2) + 2) + * (my_degree + 1)][0][0]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][0][0] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2][2]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][0][1] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2][0]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][0][2] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][2][1]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][1][0] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0][2]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][1][1] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0][0]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][1][2] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][0][1]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][2][0] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1][2]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][2][1] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1][0]; + grad_grads[i + (j + (2 * k + 9) * my_degree + + GeometryInfo::lines_per_cell) + * (my_degree + 1)][1][2][2] + = p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) + * (my_degree + 1)][1][1]; + } + } + } + + break; + } + + default: + Assert (false, ExcNotImplemented ()); + } } template -unsigned int PolynomialsNedelec::compute_n_pols (unsigned int k) { - switch (dim) { +unsigned int +PolynomialsNedelec::compute_n_pols (unsigned int k) +{ + switch (dim) + { case 1: - return k + 1; + return k + 1; case 2: - return 2 * (k + 1) * (k + 2); + return 2 * (k + 1) * (k + 2); case 3: - return 3 * (k + 1) * (k + 2) * (k + 2); + return 3 * (k + 1) * (k + 2) * (k + 2); - default: { - Assert (false, ExcNotImplemented ()); - return 0; - } - } + default: + { + Assert (false, ExcNotImplemented ()); + return 0; + } + } } diff --git a/tests/fe/check_derivatives.cc b/tests/fe/check_derivatives.cc index 637cc8e85d..e062974ef8 100644 --- a/tests/fe/check_derivatives.cc +++ b/tests/fe/check_derivatives.cc @@ -153,7 +153,7 @@ main() check (0,4); check1 (0,3); - check1 (1,1); + check1 (0,1); check1 (0,4); check1 (0,2); diff --git a/tests/fe/check_derivatives/cmp/generic b/tests/fe/check_derivatives/cmp/generic index a2b9eb7531..e3429963d8 100644 --- a/tests/fe/check_derivatives/cmp/generic +++ b/tests/fe/check_derivatives/cmp/generic @@ -80,8 +80,10 @@ DEAL::FE_DGPMonomial<3>(2) 10 OK DEAL::FE_DGPMonomial<1>(3) 4 OK DEAL::FE_DGPMonomial<2>(3) 10 OK DEAL::FE_DGPMonomial<3>(3) 20 OK -DEAL::FE_Nedelec<2>(1) 4 OK -DEAL::FE_Nedelec<3>(1) 12 OK +DEAL::FE_Nedelec<2>(0) 4 OK +DEAL::FE_Nedelec<3>(0) 12 OK +DEAL::FE_Nedelec<2>(1) 12 OK +DEAL::FE_Nedelec<3>(1) 54 OK DEAL::FE_RaviartThomas<2>(0) 4 OK DEAL::FE_RaviartThomas<2>(1) 12 OK DEAL::FE_RaviartThomas<2>(2) 24 OK -- 2.39.5