#include <sstream>
#include <iostream>
-//TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
-//adjust_line_dof_index_for_line_orientation_table fields, and write tests
-//similar to bits/face_orientation_and_fe_q_*
-
DEAL_II_NAMESPACE_OPEN
default:
Assert (false, ExcNotImplemented ());
}
-
}
+
+
// Restriction operator
template <int dim>
void
n_edge_quadrature_points = edge_quadrature.size ();
const unsigned int
index = RefinementCase<dim>::isotropic_refinement - 1;
+ const unsigned int deg = this->degree-1;
+ const std::vector<Polynomials::Polynomial<double> >&
+ legendre_polynomials = Polynomials::Legendre::generate_complete_basis (deg);
switch (dim)
{
{
// First interpolate the shape
// functions of the child cells
- // to the lowest order shape
- // functions of the parent cell.
+ // to the face shape functions
+ // of the parent cell.
for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
++q_point)
{
- const double weight = 2.0 * edge_quadrature.weight (q_point);
-
if (edge_quadrature_points[q_point] (0) < 0.5)
- {
- Point<dim> quadrature_point (0.0,
- 2.0 * edge_quadrature_points[q_point] (0));
-
- this->restriction[index][0] (0, dof) += weight
- * this->shape_value_component
- (dof,
- quadrature_point,
- 1);
- quadrature_point (0) = 1.0;
- this->restriction[index][1] (this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 1);
- quadrature_point (0) = quadrature_point (1);
- quadrature_point (1) = 0.0;
- this->restriction[index][0] (2 * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 0);
- quadrature_point (1) = 1.0;
- this->restriction[index][2] (3 * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 0);
- }
+ for (unsigned int i = 0; i <= deg; ++i)
+ {
+ const double weight
+ = 2.0 * edge_quadrature.weight (q_point)
+ * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+ Point<dim> quadrature_point (0.0,
+ 2.0 * edge_quadrature_points[q_point] (0));
+
+ this->restriction[index][0] (i, dof) += weight
+ * this->shape_value_component
+ (dof,
+ quadrature_point,
+ 1);
+ quadrature_point (0) = 1.0;
+ this->restriction[index][1] (i + this->degree, dof)
+ += weight * this->shape_value_component (dof,
+ quadrature_point,
+ 1);
+ quadrature_point (0) = quadrature_point (1);
+ quadrature_point (1) = 0.0;
+ this->restriction[index][0] (i + 2 * this->degree, dof)
+ += weight * this->shape_value_component (dof,
+ quadrature_point,
+ 0);
+ quadrature_point (1) = 1.0;
+ this->restriction[index][2] (i + 3 * this->degree, dof)
+ += weight * this->shape_value_component (dof,
+ quadrature_point,
+ 0);
+ }
else
- {
- Point<dim> quadrature_point (0.0,
- 2.0 * edge_quadrature_points[q_point] (0)
- - 1.0);
-
- this->restriction[index][2] (0, dof) += weight
- * this->shape_value_component
- (dof,
- quadrature_point,
- 1);
- quadrature_point (0) = 1.0;
- this->restriction[index][3] (this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 1);
- quadrature_point (0) = quadrature_point (1);
- quadrature_point (1) = 0.0;
- this->restriction[index][1] (2 * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 0);
- quadrature_point (1) = 1.0;
- this->restriction[index][3] (3 * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 0);
- }
+ for (unsigned int i = 0; i <= deg; ++i)
+ {
+ const double weight
+ = 2.0 * edge_quadrature.weight (q_point)
+ * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+ Point<dim> quadrature_point (0.0,
+ 2.0 * edge_quadrature_points[q_point] (0)
+ - 1.0);
+
+ this->restriction[index][2] (i, dof) += weight
+ * this->shape_value_component
+ (dof,
+ quadrature_point,
+ 1);
+ quadrature_point (0) = 1.0;
+ this->restriction[index][3] (i + this->degree, dof)
+ += weight * this->shape_value_component (dof,
+ quadrature_point,
+ 1);
+ quadrature_point (0) = quadrature_point (1);
+ quadrature_point (1) = 0.0;
+ this->restriction[index][1] (i + 2 * this->degree, dof)
+ += weight * this->shape_value_component (dof,
+ quadrature_point,
+ 0);
+ quadrature_point (1) = 1.0;
+ this->restriction[index][3] (i + 3 * this->degree, dof)
+ += weight * this->shape_value_component (dof,
+ quadrature_point,
+ 0);
+ }
}
- // Then project the shape functions
- // of the child cells to the higher
- // order shape functions of the
- // parent cell.
+ // Then interpolate the shape functions
+ // of the child cells to the interior
+ // shape functions of the parent cell.
if (this->degree > 1)
{
- const unsigned int deg = this->degree-1;
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (deg);
- FullMatrix<double> system_matrix_inv (deg, deg);
-
- {
- FullMatrix<double> assembling_matrix (deg,
- n_edge_quadrature_points);
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
- {
- const double weight
- = std::sqrt (edge_quadrature.weight (q_point));
-
- for (unsigned int i = 0; i < deg; ++i)
- assembling_matrix (i, q_point) = weight
- * legendre_polynomials[i + 1].value
- (edge_quadrature_points[q_point] (0));
- }
-
- FullMatrix<double> system_matrix (deg, deg);
-
- assembling_matrix.mTmult (system_matrix, assembling_matrix);
- system_matrix_inv.invert (system_matrix);
- }
-
- FullMatrix<double> solution (this->degree-1, 4);
- FullMatrix<double> system_rhs (this->degree-1, 4);
- Vector<double> tmp (4);
-
- for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
- for (unsigned int i = 0; i < 2; ++i)
- {
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
- {
- const double weight
- = edge_quadrature.weight (q_point);
- const Point<dim> quadrature_point_0 (i,
- edge_quadrature_points[q_point] (0));
- const Point<dim> quadrature_point_1
- (edge_quadrature_points[q_point] (0),
- i);
-
- if (edge_quadrature_points[q_point] (0) < 0.5)
- {
- Point<dim> quadrature_point_2 (i,
- 2.0 * edge_quadrature_points[q_point] (0));
-
- tmp (0) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_2, 1)
- - this->restriction[index][i]
- (i * this->degree, dof)
- * this->shape_value_component
- (i * this->degree,
- quadrature_point_0, 1));
- tmp (1) = -1.0 * weight
- * this->restriction[index][i + 2]
- (i * this->degree, dof)
- * this->shape_value_component
- (i * this->degree,
- quadrature_point_0, 1);
- quadrature_point_2
- = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
- i);
- tmp (2) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_2, 0)
- - this->restriction[index][2 * i]
- ((i + 2) * this->degree, dof)
- * this->shape_value_component
- ((i + 2) * this->degree,
- quadrature_point_1, 0));
- tmp (3) = -1.0 * weight
- * this->restriction[index][2 * i + 1]
- ((i + 2) * this->degree, dof)
- * this->shape_value_component
- ((i + 2) * this->degree,
- quadrature_point_1, 0);
- }
-
- else
- {
- tmp (0) = -1.0 * weight
- * this->restriction[index][i]
- (i * this->degree, dof)
- * this->shape_value_component
- (i * this->degree,
- quadrature_point_0, 1);
-
- Point<dim> quadrature_point_2 (i,
- 2.0 * edge_quadrature_points[q_point] (0)
- - 1.0);
-
- tmp (1) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_2, 1)
- - this->restriction[index][i + 2]
- (i * this->degree, dof)
- * this->shape_value_component
- (i * this->degree,
- quadrature_point_0, 1));
- tmp (2) = -1.0 * weight
- * this->restriction[index][2 * i]
- ((i + 2) * this->degree, dof)
- * this->shape_value_component
- ((i + 2) * this->degree,
- quadrature_point_1, 0);
- quadrature_point_2
- = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
- - 1.0, i);
- tmp (3) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_2, 0)
- - this->restriction[index][2 * i + 1]
- ((i + 2) * this->degree, dof)
- * this->shape_value_component
- ((i + 2) * this->degree,
- quadrature_point_1, 0));
- }
-
- for (unsigned int j = 0; j < this->degree-1; ++j)
- {
- const double L_j
- = legendre_polynomials[j + 1].value
- (edge_quadrature_points[q_point] (0));
-
- for (unsigned int k = 0; k < tmp.size (); ++k)
- system_rhs (j, k) += tmp (k) * L_j;
- }
- }
-
- system_matrix_inv.mmult (solution, system_rhs);
-
- for (unsigned int j = 0; j < this->degree-1; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- {
- if (std::abs (solution (j, k)) > 1e-14)
- this->restriction[index][i + 2 * k]
- (i * this->degree + j + 1, dof)
- = solution (j, k);
-
- if (std::abs (solution (j, k + 2)) > 1e-14)
- this->restriction[index][2 * i + k]
- ((i + 2) * this->degree + j + 1, dof)
- = solution (j, k + 2);
- }
- }
-
const QGauss<dim> quadrature (2 * this->degree);
const std::vector<Point<dim> >&
quadrature_points = quadrature.get_points ();
(this->degree);
const unsigned int n_boundary_dofs
= GeometryInfo<dim>::faces_per_cell * this->degree;
- const unsigned int& n_quadrature_points = quadrature.size ();
- {
- FullMatrix<double> assembling_matrix ((this->degree-1) * this->degree,
- n_quadrature_points);
-
- for (unsigned int q_point = 0; q_point < n_quadrature_points;
+ for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
+ for (unsigned int q_point = 0; q_point < quadrature.size ();
++q_point)
{
- const double weight
- = std::sqrt (quadrature.weight (q_point));
+ const double weight = 2.0 * quadrature.weight (q_point);
+ Point<dim> quadrature_point (2.0 * quadrature_points[q_point] (0),
+ 2.0 * quadrature_points[q_point] (1));
+ unsigned int k;
- for (unsigned int i = 0; i < this->degree; ++i)
+ if (quadrature_points[q_point] (0) < 0.5)
{
- const double L_i = weight
- * legendre_polynomials[i].value
- (quadrature_points[q_point] (0));
-
- for (unsigned int j = 0; j < this->degree-1; ++j)
- assembling_matrix (i * (this->degree-1) + j, q_point)
- = L_i * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (1));
- }
- }
-
- FullMatrix<double> system_matrix (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 (), 8);
- system_rhs.reinit (system_matrix_inv.m (), 8);
- tmp.reinit (8);
-
- for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
- {
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_quadrature_points; ++q_point)
- {
- tmp = 0.0;
-
- if (quadrature_points[q_point] (0) < 0.5)
- {
- if (quadrature_points[q_point] (1) < 0.5)
- {
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0),
- 2.0 * quadrature_points[q_point] (1));
-
- tmp (0) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (1) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- }
-
- else
- {
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0),
- 2.0 * quadrature_points[q_point] (1)
- - 1.0);
-
- tmp (4) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (5) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- }
- }
-
- else
if (quadrature_points[q_point] (1) < 0.5)
- {
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * quadrature_points[q_point] (1));
-
- tmp (2) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (3) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- }
-
+ k = 0;
+
else
{
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * quadrature_points[q_point] (1)
- - 1.0);
-
- tmp (6) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (7) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- }
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- {
- tmp (2 * i) -= this->restriction[index][i]
- (j + 2 * this->degree, dof)
- * this->shape_value_component
- (j + 2 * this->degree,
- quadrature_points[q_point], 0);
- tmp (2 * i + 1) -= this->restriction[index][i]
- (i * this->degree + j, dof)
- * this->shape_value_component
- (i * this->degree + j,
- quadrature_points[q_point], 1);
- tmp (2 * (i + 2)) -= this->restriction[index][i + 2]
- (j + 3 * this->degree, dof)
- * this->shape_value_component
- (j + 3 * this->degree,
- quadrature_points[q_point],
- 0);
- tmp (2 * i + 5) -= this->restriction[index][i + 2]
- (i * this->degree + j, dof)
- * this->shape_value_component
- (i * this->degree + j,
- quadrature_points[q_point], 1);
+ quadrature_point (1) -= 1.0;
+ k = 1;
}
-
- tmp *= quadrature.weight (q_point);
-
- for (unsigned int i = 0; i < this->degree; ++i)
+ }
+
+ else
+ if (quadrature_points[q_point] (1) < 0.5)
{
- const double L_i_0
- = legendre_polynomials[i].value
- (quadrature_points[q_point] (0));
- const double L_i_1
- = legendre_polynomials[i].value
- (quadrature_points[q_point] (1));
-
- for (unsigned int j = 0; j < this->degree-1; ++j)
- {
- const double l_j_0
- = L_i_0 * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (1));
- const double l_j_1
- = L_i_1 * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (0));
-
- for (unsigned int k = 0; k < 4; ++k)
- {
- system_rhs (i * (this->degree-1) + j, 2 * k)
- += tmp (2 * k) * l_j_0;
- system_rhs (i * (this->degree-1) + j, 2 * k + 1)
- += tmp (2 * k + 1) * l_j_1;
- }
- }
+ quadrature_point (0) -= 1.0;
+ k = 2;
}
- }
-
- system_matrix_inv.mmult (solution, system_rhs);
-
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- for (unsigned int k = 0; k < 4; ++k)
+
+ else
{
- if (std::abs (solution (i * (this->degree-1) + j, 2 * k))
- > 1e-14)
- this->restriction[index][k]
- (i * (this->degree-1) + j + n_boundary_dofs, dof)
- = solution (i * (this->degree-1) + j, 2 * k);
-
- if (std::abs (solution (i * (this->degree-1) + j, 2 * k + 1))
- > 1e-14)
- this->restriction[index][k]
- (i + (this->degree-1 + j) * this->degree + n_boundary_dofs,
- dof)
- = solution (i * (this->degree-1) + j, 2 * k + 1);
+ quadrature_point (0) -= 1.0;
+ quadrature_point (1) -= 1.0;
+ k = 3;
}
- }
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ {
+ const double L_i_0 = weight
+ * legendre_polynomials[i].value
+ (quadrature_points[q_point] (0));
+ const double L_i_1 = weight
+ * legendre_polynomials[i].value
+ (quadrature_points[q_point] (1));
+
+ for (unsigned int j = 0; j < deg; ++j)
+ {
+ this->restriction[index][k] (i * deg + j + n_boundary_dofs, dof)
+ += L_i_0 * (this->shape_grad_component (dof, quadrature_point, 0)[1]
+ * legendre_polynomials[j + 1].value
+ (quadrature_points[q_point] (1))
+ + this->shape_value_component (dof, quadrature_point, 0)
+ * lobatto_polynomials[j + 2].value
+ (quadrature_points[q_point] (1)));
+ this->restriction[index][k] (i + (deg + j) * this->degree
+ + n_boundary_dofs, dof)
+ += L_i_1 * (this->shape_grad_component (dof, quadrature_point, 1)[0]
+ * legendre_polynomials[j + 1].value
+ (quadrature_points[q_point] (0))
+ + this->shape_value_component (dof, quadrature_point, 1)
+ * lobatto_polynomials[j + 2].value
+ (quadrature_points[q_point] (0)));
+ }
+ }
+ }
}
break;
{
// First interpolate the shape
// functions of the child cells
- // to the lowest order shape
- // functions of the parent cell.
+ // to the edge shape functions
+ // of the parent cell.
for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
++q_point)
{
const double weight = 2.0 * edge_quadrature.weight (q_point);
-
+
if (edge_quadrature_points[q_point] (0) < 0.5)
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- {
- Point<dim> quadrature_point (i,
- 2.0 * edge_quadrature_points[q_point] (0),
- j);
-
- this->restriction[index][i + 4 * j]
- ((i + 4 * j) * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 1);
- quadrature_point
- = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
- i, j);
- this->restriction[index][2 * (i + 2 * j)]
- ((i + 4 * j + 2) * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 0);
- quadrature_point = Point<dim> (i, j,
- 2.0 * edge_quadrature_points[q_point] (0));
- this->restriction[index][i + 2 * j]
- ((i + 2 * (j + 4)) * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 2);
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ {
+ const double L_i
+ = weight * legendre_polynomials[i].value
+ (edge_quadrature_points[q_point] (0));
- else
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int k = 0; k < 2; ++k)
+ {
+ Point<dim> quadrature_point (j,
+ 2.0 * edge_quadrature_points[q_point] (0),
+ k);
+
+ this->restriction[index][j + 4 * k]
+ (i + (j + 4 * k) * this->degree, dof)
+ += L_i * this->shape_value_component (dof,
+ quadrature_point,
+ 1);
+ quadrature_point
+ = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
+ j, k);
+ this->restriction[index][2 * (j + 2 * k)]
+ (i + (j + 4 * k + 2) * this->degree, dof)
+ += L_i * this->shape_value_component (dof,
+ quadrature_point,
+ 0);
+ quadrature_point = Point<dim> (j, k,
+ 2.0 * edge_quadrature_points[q_point] (0));
+ this->restriction[index][j + 2 * k]
+ (i + (j + 2 * (k + 4)) * this->degree, dof)
+ += L_i * this->shape_value_component (dof,
+ quadrature_point,
+ 2);
+ }
+ }
+
+ else
+ for (unsigned int i = 0; i < this->degree; ++i)
+ {
+ const double L_i
+ = weight * legendre_polynomials[i].value
+ (edge_quadrature_points[q_point] (0));
+
+ for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int k = 0; k < 2; ++k)
{
- Point<dim> quadrature_point (i,
+ Point<dim> quadrature_point (j,
2.0 * edge_quadrature_points[q_point] (0)
- - 1.0, j);
+ - 1.0, k);
- this->restriction[index][i + 4 * j + 2]
- ((i + 4 * j) * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 1);
+ this->restriction[index][j + 4 * k + 2]
+ (i + (j + 4 * k) * this->degree, dof)
+ += L_i * this->shape_value_component (dof,
+ quadrature_point,
+ 1);
quadrature_point
- = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
- - 1.0, i, j);
- this->restriction[index][2 * (i + 2 * j) + 1]
- ((i + 4 * j + 2) * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 0);
- quadrature_point = Point<dim> (i, j,
+ = Point<dim> (2.0 * edge_quadrature_points[q_point] (0) - 1.0,
+ j, k);
+ this->restriction[index][2 * (j + 2 * k) + 1]
+ (i + (j + 4 * k + 2) * this->degree, dof)
+ += L_i * this->shape_value_component (dof,
+ quadrature_point,
+ 0);
+ quadrature_point = Point<dim> (j, k,
2.0 * edge_quadrature_points[q_point] (0)
- 1.0);
- this->restriction[index][i + 2 * (j + 2)]
- ((i + 2 * (j + 4)) * this->degree, dof)
- += weight * this->shape_value_component (dof,
- quadrature_point,
- 2);
+ this->restriction[index][j + 2 * (k + 2)]
+ (i + (j + 2 * (k + 4)) * this->degree, dof)
+ += L_i * this->shape_value_component (dof,
+ quadrature_point,
+ 2);
}
+ }
}
- // Then project the shape functions
- // of the child cells to the higher
- // order shape functions of the
- // parent cell.
+ // Then interpolate the shape functions
+ // of the child cells to the face
+ // and interior shape functions of
+ // the parent cell.
if (this->degree > 1)
{
- const unsigned int deg = this->degree-1;
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (deg);
- FullMatrix<double> system_matrix_inv (deg, deg);
-
- {
- FullMatrix<double> assembling_matrix (deg,
- n_edge_quadrature_points);
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
- {
- const double weight = std::sqrt (edge_quadrature.weight
- (q_point));
-
- for (unsigned int i = 0; i < deg; ++i)
- assembling_matrix (i, q_point) = weight
- * legendre_polynomials[i + 1].value
- (edge_quadrature_points[q_point] (0));
- }
-
- FullMatrix<double> system_matrix (deg, deg);
-
- assembling_matrix.mTmult (system_matrix, assembling_matrix);
- system_matrix_inv.invert (system_matrix);
- }
-
- FullMatrix<double> solution (deg, 6);
- FullMatrix<double> system_rhs (deg, 6);
- Vector<double> tmp (6);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
- {
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
- {
- const double weight = edge_quadrature.weight
- (q_point);
- const Point<dim> quadrature_point_0 (i,
- edge_quadrature_points[q_point] (0),
- j);
- const Point<dim>
- quadrature_point_1
- (edge_quadrature_points[q_point] (0), i, j);
- const Point<dim> quadrature_point_2 (i, j,
- edge_quadrature_points[q_point] (0));
-
- if (edge_quadrature_points[q_point] (0) < 0.5)
- {
- Point<dim> quadrature_point_3 (i,
- 2.0 * edge_quadrature_points[q_point] (0),
- j);
-
- tmp (0) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_3, 1)
- - this->restriction[index][i + 4 * j]
- ((i + 4 * j) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j) * this->degree,
- quadrature_point_0, 1));
- tmp (1) = -1.0 * weight
- * this->restriction[index][i + 4 * j + 2]
- ((i + 4 * j) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j) * this->degree,
- quadrature_point_0, 1);
- quadrature_point_3
- = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
- i, j);
- tmp (2) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_3, 0)
- - this->restriction[index][2 * (i + 2 * j)]
- ((i + 4 * j + 2) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j + 2) * this->degree,
- quadrature_point_1, 0));
- tmp (3) = -1.0 * weight
- * this->restriction[index][2 * (i + 2 * j) + 1]
- ((i + 4 * j + 2) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j + 2) * this->degree,
- quadrature_point_1, 0);
- quadrature_point_3 = Point<dim> (i, j,
- 2.0 * edge_quadrature_points[q_point] (0));
- tmp (4) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_3, 2)
- - this->restriction[index][i + 2 * j]
- ((i + 2 * (j + 4)) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 2 * (j + 4)) * this->degree,
- quadrature_point_2, 2));
- tmp (5) = -1.0 * weight
- * this->restriction[index][i + 2 * (j + 2)]
- ((i + 2 * (j + 4)) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 2 * (j + 4)) * this->degree,
- quadrature_point_2, 2);
- }
-
- else
- {
- tmp (0) = -1.0 * weight
- * this->restriction[index][i + 4 * j]
- ((i + 4 * j) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j) * this->degree,
- quadrature_point_0, 1);
-
- Point<dim> quadrature_point_3 (i,
- 2.0 * edge_quadrature_points[q_point] (0)
- - 1.0, j);
-
- tmp (1) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_3, 1)
- - this->restriction[index][i + 4 * j + 2]
- ((i + 4 * j) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j) * this->degree,
- quadrature_point_0, 1));
- tmp (2) = -1.0 * weight
- * this->restriction[index][2 * (i + 2 * j)]
- ((i + 4 * j + 2) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j + 2) * this->degree,
- quadrature_point_1, 0);
- quadrature_point_3
- = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
- - 1.0, i, j);
- tmp (3) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_3, 0)
- - this->restriction[index][2 * (i + 2 * j) + 1]
- ((i + 4 * j + 2) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 4 * j + 2) * this->degree,
- quadrature_point_1, 0));
- tmp (4) = -1.0 * weight
- * this->restriction[index][i + 2 * j]
- ((i + 2 * (j + 4)) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 2 * (j + 4)) * this->degree,
- quadrature_point_2, 2);
- quadrature_point_3 = Point<dim> (i, j,
- 2.0 * edge_quadrature_points[q_point] (0)
- - 1.0);
- tmp (5) = weight
- * (2.0 * this->shape_value_component
- (dof, quadrature_point_3, 2)
- - this->restriction[index][i + 2 * (j + 2)]
- ((i + 2 * (j + 4)) * this->degree,
- dof)
- * this->shape_value_component
- ((i + 2 * (j + 4)) * this->degree,
- quadrature_point_2, 2));
- }
-
- for (unsigned int k = 0; k < deg; ++k)
- {
- const double L_k
- = legendre_polynomials[k + 1].value
- (edge_quadrature_points[q_point] (0));
-
- for (unsigned int l = 0; l < tmp.size (); ++l)
- system_rhs (k, l) += tmp (l) * L_k;
- }
- }
-
- system_matrix_inv.mmult (solution, system_rhs);
-
- for (unsigned int k = 0; k < 2; ++k)
- for (unsigned int l = 0; l < deg; ++l)
- {
- if (std::abs (solution (l, k)) > 1e-14)
- this->restriction[index][i + 2 * (2 * j + k)]
- ((i + 4 * j) * this->degree + l + 1, dof)
- = solution (l, k);
-
- if (std::abs (solution (l, k + 2)) > 1e-14)
- this->restriction[index][2 * (i + 2 * j) + k]
- ((i + 4 * j + 2) * this->degree + l + 1, dof)
- = solution (l, k + 2);
-
- if (std::abs (solution (l, k + 4)) > 1e-14)
- this->restriction[index][i + 2 * (j + 2 * k)]
- ((i + 2 * (j + 4)) * this->degree + l + 1,
- dof)
- = solution (l, k + 4);
- }
- }
-
const QGauss<2> face_quadrature (2 * this->degree);
const std::vector<Point<2> >& face_quadrature_points
= face_quadrature.get_points ();
lobatto_polynomials
= Polynomials::Lobatto::generate_complete_basis
(this->degree);
+ const QGauss<dim> quadrature (2 * this->degree);
+ const std::vector<Point<dim> >&
+ quadrature_points = quadrature.get_points ();
const unsigned int n_edge_dofs
= GeometryInfo<dim>::lines_per_cell * this->degree;
+ const unsigned int n_boundary_dofs
+ = 2 * GeometryInfo<dim>::faces_per_cell * deg * this->degree
+ + n_edge_dofs;
const unsigned int& n_face_quadrature_points
= face_quadrature.size ();
-
- {
- FullMatrix<double> assembling_matrix
- (deg * this->degree,
- n_face_quadrature_points);
-
- for (unsigned int q_point = 0;
- q_point < n_face_quadrature_points; ++q_point)
- {
- const double weight
- = std::sqrt (face_quadrature.weight (q_point));
-
- for (unsigned int i = 0; i <= deg; ++i)
- {
- const double L_i = weight
- * legendre_polynomials[i].value
- (face_quadrature_points[q_point] (0));
-
- for (unsigned int j = 0; j < deg; ++j)
- assembling_matrix (i * deg + j, q_point)
- = L_i * lobatto_polynomials[j + 2].value
- (face_quadrature_points[q_point] (1));
- }
- }
-
- FullMatrix<double> system_matrix (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 (), 24);
- system_rhs.reinit (system_matrix_inv.m (), 24);
- tmp.reinit (24);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
- {
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_face_quadrature_points; ++q_point)
- {
- tmp = 0.0;
-
- if (face_quadrature_points[q_point] (0) < 0.5)
- {
- if (face_quadrature_points[q_point] (1) < 0.5)
- {
- Point<dim> quadrature_point_0 (i,
- 2.0 * face_quadrature_points[q_point] (0),
- 2.0 * face_quadrature_points[q_point] (1));
-
- tmp (0) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
- tmp (1) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
- i,
- 2.0 * face_quadrature_points[q_point] (1));
- tmp (8) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- tmp (9) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
- 2.0 * face_quadrature_points[q_point] (1),
- i);
- tmp (16) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- tmp (17) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
- }
-
- else
- {
- Point<dim> quadrature_point_0 (i,
- 2.0 * face_quadrature_points[q_point] (0),
- 2.0 * face_quadrature_points[q_point] (1)
- - 1.0);
-
- tmp (2) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
- tmp (3) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
- i,
- 2.0 * face_quadrature_points[q_point] (1)
- - 1.0);
- tmp (10) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- tmp (11) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
- 2.0 * face_quadrature_points[q_point] (1)
- - 1.0, i);
- tmp (18) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- tmp (19) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
- }
- }
-
- else
- if (face_quadrature_points[q_point] (1) < 0.5)
+ const unsigned int& n_quadrature_points = quadrature.size ();
+
+ // First, the interpolation to
+ // the face shape functions.
+ for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
+ {
+ for (unsigned int q_point = 0; q_point < n_face_quadrature_points;
+ ++q_point)
+ {
+ const double weight = 2.0 * face_quadrature.weight (q_point);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ {
+ Point<dim> quadrature_point (i,
+ 2.0 * face_quadrature_points[q_point] (0),
+ 2.0 * face_quadrature_points[q_point] (1));
+ unsigned int l;
+ unsigned int m;
+
+ if (face_quadrature_points[q_point] (0) < 0.5)
{
- Point<dim> quadrature_point_0 (i,
- 2.0 * face_quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * face_quadrature_points[q_point] (1));
-
- tmp (4) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
- tmp (5) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
- - 1.0, i,
- 2.0 * face_quadrature_points[q_point] (1));
- tmp (12) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- tmp (13) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * face_quadrature_points[q_point] (1),
- i);
- tmp (20) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- tmp (21) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
+ m = 0;
+
+ if (face_quadrature_points[q_point] (1) < 0.5)
+ l = 0;
+
+ else
+ {
+ quadrature_point (2) -= 1.0;
+ l = 1;
+ }
}
-
+
else
{
- Point<dim> quadrature_point_0 (i,
- 2.0 * face_quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * face_quadrature_points[q_point] (1)
- - 1.0);
-
- tmp (6) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
- tmp (7) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
- - 1.0, i,
- 2.0 * face_quadrature_points[q_point] (1)
- - 1.0);
- tmp (14) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 2);
- tmp (15) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- quadrature_point_0
- = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * face_quadrature_points[q_point] (1)
- - 1.0, i);
- tmp (22) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 0);
- tmp (23) += 2.0 * this->shape_value_component
- (dof, quadrature_point_0, 1);
+ quadrature_point (1) -= 1.0;
+ m = 1;
+
+ if (face_quadrature_points[q_point] (1) < 0.5)
+ l = 0;
+
+ else
+ {
+ quadrature_point (2) -= 1.0;
+ l = 1;
+ }
}
-
- const Point<dim> quadrature_point_0 (i,
- face_quadrature_points[q_point] (0),
- face_quadrature_points[q_point] (1));
- const Point<dim> quadrature_point_1
- (face_quadrature_points[q_point] (0),
- i,
- face_quadrature_points[q_point] (1));
- const Point<dim> quadrature_point_2
- (face_quadrature_points[q_point] (0),
- face_quadrature_points[q_point] (1),
- i);
-
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- for (unsigned int l = 0; l <= deg; ++l)
- {
- tmp (2 * (j + 2 * k))
- -= this->restriction[index][i + 2 * (2 * j + k)]
- ((i + 4 * j) * this->degree + l, dof)
- * this->shape_value_component
- ((i + 4 * j) * this->degree + l,
- quadrature_point_0, 1);
- tmp (2 * (j + 2 * k) + 1)
- -= this->restriction[index][i + 2 * (2 * j + k)]
- ((i + 2 * (k + 4)) * this->degree + l,
- dof)
- * this->shape_value_component
- ((i + 2 * (k + 4)) * this->degree + l,
- quadrature_point_0, 2);
- tmp (2 * (j + 2 * (k + 2)))
- -= this->restriction[index][2 * (i + 2 * j) + k]
- ((2 * (i + 4) + k) * this->degree + l,
- dof)
- * this->shape_value_component
- ((2 * (i + 4) + k) * this->degree + l,
- quadrature_point_1, 2);
- tmp (2 * (j + 2 * k) + 9)
- -= this->restriction[index][2 * (i + 2 * j) + k]
- ((i + 4 * j + 2) * this->degree + l,
- dof)
- * this->shape_value_component
- ((i + 4 * j + 2) * this->degree + l,
- quadrature_point_1, 0);
- tmp (2 * (j + 2 * (k + 4)))
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((4 * i + j + 2) * this->degree + l,
- dof)
- * this->shape_value_component
- ((4 * i + j + 2) * this->degree + l,
- quadrature_point_2, 0);
- tmp (2 * (j + 2 * k) + 17)
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((4 * i + k) * this->degree + l, dof)
- * this->shape_value_component
- ((4 * i + k) * this->degree + l,
- quadrature_point_2, 1);
- }
-
- tmp *= face_quadrature.weight (q_point);
-
- for (unsigned int j = 0; j <= deg; ++j)
- {
- const double L_j_0
- = legendre_polynomials[j].value
- (face_quadrature_points[q_point] (0));
- const double L_j_1
- = legendre_polynomials[j].value
- (face_quadrature_points[q_point] (1));
-
- for (unsigned int k = 0; k < deg; ++k)
- {
- const double l_k_0
- = L_j_0 * lobatto_polynomials[k + 2].value
- (face_quadrature_points[q_point] (1));
- const double l_k_1
- = L_j_1 * lobatto_polynomials[k + 2].value
- (face_quadrature_points[q_point] (0));
-
- for (unsigned int l = 0; l < 4; ++l)
- {
- system_rhs (j * deg + k, 2 * l)
- += tmp (2 * l) * l_k_0;
- system_rhs (j * deg + k, 2 * l + 1)
- += tmp (2 * l + 1) * l_k_1;
- system_rhs (j * deg + k, 2 * (l + 4))
- += tmp (2 * (l + 4)) * l_k_1;
- system_rhs (j * deg + k, 2 * l + 9)
- += tmp (2 * l + 9) * l_k_0;
- system_rhs (j * deg + k, 2 * (l + 8))
- += tmp (2 * (l + 8)) * l_k_0;
- system_rhs (j * deg + k, 2 * l + 17)
- += tmp (2 * l + 17) * l_k_1;
- }
- }
- }
- }
-
- system_matrix_inv.mmult (solution, system_rhs);
-
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- for (unsigned int l = 0; l <= deg; ++l)
- for (unsigned int m = 0; m < deg; ++m)
+
+ for (unsigned int j = 0; j < this->degree; ++j)
{
- if (std::abs (solution (l * deg + m,
- 2 * (j + 2 * k)))
- > 1e-14)
- this->restriction[index][i + 2 * (2 * j + k)]
- ((2 * i * this->degree + l) * deg + m
- + n_edge_dofs,
- dof) = solution (l * deg + m,
- 2 * (j + 2 * k));
-
- if (std::abs (solution (l * deg + m,
- 2 * (j + 2 * k) + 1))
- > 1e-14)
- this->restriction[index][i + 2 * (2 * j + k)]
- (((2 * i + 1) * deg + m) * this->degree + l
- + n_edge_dofs, dof)
- = solution (l * deg + m,
- 2 * (j + 2 * k) + 1);
-
- if (std::abs (solution (l * deg + m,
- 2 * (j + 2 * (k + 2))))
- > 1e-14)
- this->restriction[index][2 * (i + 2 * j) + k]
- ((2 * (i + 2) * this->degree + l) * deg + m
- + n_edge_dofs,
- dof) = solution (l * deg + m,
- 2 * (j + 2 * (k + 2)));
-
- if (std::abs (solution (l * deg + m,
- 2 * (j + 2 * k) + 9))
- > 1e-14)
- this->restriction[index][2 * (i + 2 * j) + k]
- (((2 * i + 5) * deg + m) * this->degree + l
- + n_edge_dofs, dof)
- = solution (l * deg + m,
- 2 * (j + 2 * k) + 9);
-
- if (std::abs (solution (l * deg + m,
- 2 * (j + 2 * (k + 4))))
- > 1e-14)
- this->restriction[index][2 * (2 * i + j) + k]
- ((2 * (i + 4) * this->degree + l) * deg + m
- + n_edge_dofs,
- dof) = solution (l * deg + m,
- 2 * (j + 2 * (k + 4)));
-
- if (std::abs (solution (l * deg + m,
- 2 * (j + 2 * k) + 17))
- > 1e-14)
- this->restriction[index][2 * (2 * i + j) + k]
- (((2 * i + 9) * deg + m) * this->degree + l
- + n_edge_dofs, dof)
- = solution (l * deg + m,
- 2 * (j + 2 * k) + 17);
- }
- }
-
- const QGauss<dim> quadrature (2 * this->degree);
- const std::vector<Point<dim> >&
- quadrature_points = quadrature.get_points ();
- const unsigned int n_boundary_dofs
- = 2 * GeometryInfo<dim>::faces_per_cell * deg * this->degree
- + n_edge_dofs;
- const unsigned int& n_quadrature_points = quadrature.size ();
-
- {
- FullMatrix<double>
- assembling_matrix (deg * deg * this->degree,
- n_quadrature_points);
-
- for (unsigned int q_point = 0; q_point < n_quadrature_points;
- ++q_point)
- {
- const double weight = std::sqrt (quadrature.weight
- (q_point));
-
- for (unsigned int i = 0; i <= deg; ++i)
- {
- const double L_i = weight
- * legendre_polynomials[i].value
- (quadrature_points[q_point] (0));
-
- for (unsigned int j = 0; j < deg; ++j)
- {
- const double l_j
- = L_i * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (1));
-
- for (unsigned int k = 0; k < deg; ++k)
- assembling_matrix ((i * deg + j) * deg + k,
- q_point)
- = l_j * lobatto_polynomials[k + 2].value
- (quadrature_points[q_point] (2));
- }
- }
- }
-
- FullMatrix<double> system_matrix (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 (), 24);
- system_rhs.reinit (system_matrix_inv.m (), 24);
- tmp.reinit (24);
-
- for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
- {
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_quadrature_points; ++q_point)
+ const double L_j_0
+ = weight * legendre_polynomials[j].value
+ (face_quadrature_points[q_point] (0));
+ const double L_j_1
+ = weight * legendre_polynomials[j].value
+ (face_quadrature_points[q_point] (1));
+
+ for (unsigned int k = 0; k < deg; ++k)
+ {
+ const double Le_k_0
+ = legendre_polynomials[k + 1].value
+ (face_quadrature_points[q_point] (0));
+ const double Le_k_1
+ = legendre_polynomials[k + 1].value
+ (face_quadrature_points[q_point] (1));
+ const double lo_k_0
+ = lobatto_polynomials[k + 2].value
+ (face_quadrature_points[q_point] (0));
+ const double lo_k_1
+ = lobatto_polynomials[k + 2].value
+ (face_quadrature_points[q_point] (1));
+
+ this->restriction[index][i + 2 * (2 * l + m)]
+ ((2 * i * this->degree + j) * deg + k + n_edge_dofs, dof)
+ += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 1)[2]
+ * Le_k_1
+ + this->shape_value_component (dof, quadrature_point, 1)
+ * lo_k_1);
+ this->restriction[index][i + 2 * (2 * l + m)]
+ (((2 * i + 1) * deg + k) * this->degree + j + n_edge_dofs, dof)
+ += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 2)[1]
+ * Le_k_0
+ + this->shape_value_component (dof, quadrature_point, 2)
+ * lo_k_0);
+ this->restriction[index][2 * (i + 2 * l) + m]
+ ((2 * (i + 2) * this->degree + j) * deg + k + n_edge_dofs, dof)
+ += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 2)[0]
+ * Le_k_1
+ + this->shape_value_component (dof, quadrature_point, 2)
+ * lo_k_1);
+ this->restriction[index][2 * (i + 2 * l) + m]
+ (((2 * i + 5) * deg + k) * this->degree + j + n_edge_dofs, dof)
+ += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 0)[2]
+ * Le_k_0
+ + this->shape_value_component (dof, quadrature_point, 0)
+ * lo_k_0);
+ this->restriction[index][2 * (2 * i + l) + m]
+ ((2 * (i + 4) * this->degree + j) * deg + k + n_edge_dofs, dof)
+ += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 0)[1]
+ * Le_k_1
+ + this->shape_value_component (dof, quadrature_point, 0)
+ * lo_k_1);
+ this->restriction[index][2 * (2 * i + l) + m]
+ (((2 * i + 9) * deg + k) * this->degree + j + n_edge_dofs, dof)
+ += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 1)[0]
+ * Le_k_0
+ + this->shape_value_component (dof, quadrature_point, 1)
+ * lo_k_0);
+ }
+ }
+ }
+ }
+
+ // Then, the interpolation to
+ // the interior shape functions.
+ for (unsigned int q_point = 0; q_point < n_quadrature_points;
+ ++q_point)
{
- tmp = 0.0;
-
+ Point<dim> quadrature_point = 2.0 * quadrature_points[q_point];
+ unsigned int l;
+ unsigned int m;
+ unsigned int n;
+
if (quadrature_points[q_point] (0) < 0.5)
{
+ n = 0;
+
if (quadrature_points[q_point] (1) < 0.5)
{
+ m = 0;
+
if (quadrature_points[q_point] (2) < 0.5)
+ l = 0;
+
+ else
{
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0),
- 2.0 * quadrature_points[q_point] (1),
- 2.0 * quadrature_points[q_point] (2));
-
- tmp (0) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (1) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (2) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
+ quadrature_point (2) -= 1.0;
+ l = 1;
}
-
+ }
+
+ else
+ {
+ quadrature_point (1) -= 1.0;
+ m = 1;
+
+ if (quadrature_points[q_point] (2) < 0.5)
+ l = 0;
+
else
{
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0),
- 2.0 * quadrature_points[q_point] (1),
- 2.0 * quadrature_points[q_point] (2)
- - 1.0);
-
- tmp (3) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (4) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (5) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
+ quadrature_point (2) -= 1.0;
+ l = 1;
}
}
-
- else
- if (quadrature_points[q_point] (2) < 0.5)
- {
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0),
- 2.0 * quadrature_points[q_point] (1)
- - 1.0,
- 2.0 * quadrature_points[q_point] (2));
-
- tmp (6) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (7) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (8) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
- }
-
- else
- {
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0),
- 2.0 * quadrature_points[q_point] (1)
- - 1.0,
- 2.0 * quadrature_points[q_point] (2)
- - 1.0);
-
- tmp (9) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (10) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (11) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
- }
}
-
+
else
- if (quadrature_points[q_point] (1) < 0.5)
- {
- if (quadrature_points[q_point] (2) < 0.5)
- {
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * quadrature_points[q_point] (1),
- 2.0 * quadrature_points[q_point] (2));
-
- tmp (12) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (13) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (14) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
- }
-
- else
- {
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * quadrature_points[q_point] (1),
- 2.0 * quadrature_points[q_point] (2)
- - 1.0);
-
- tmp (15) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (16) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (17) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
- }
- }
-
- else
- if (quadrature_points[q_point] (2) < 0.5)
+ {
+ quadrature_point (0) -= 1.0;
+ n = 1;
+
+ if (quadrature_points[q_point] (1) < 0.5)
{
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * quadrature_points[q_point] (1)
- - 1.0,
- 2.0 * quadrature_points[q_point] (2));
-
- tmp (18) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (19) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (20) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
+ m = 0;
+
+ if (quadrature_points[q_point] (2) < 0.5)
+ l = 0;
+
+ else
+ {
+ quadrature_point (2) -= 1.0;
+ l = 1;
+ }
}
-
+
else
{
- const Point<dim> quadrature_point
- (2.0 * quadrature_points[q_point] (0)
- - 1.0,
- 2.0 * quadrature_points[q_point] (1)
- - 1.0,
- 2.0 * quadrature_points[q_point] (2)
- - 1.0);
-
- tmp (21) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 0);
- tmp (22) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 1);
- tmp (23) += 2.0 * this->shape_value_component
- (dof, quadrature_point, 2);
+ quadrature_point (1) -= 1.0;
+ m = 1;
+
+ if (quadrature_points[q_point] (2) < 0.5)
+ l = 0;
+
+ else
+ {
+ quadrature_point (2) -= 1.0;
+ l = 1;
+ }
}
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- for (unsigned int l = 0; l <= deg; ++l)
- {
- tmp (3 * (i + 2 * (j + 2 * k)))
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((4 * i + j + 2) * this->degree + l, dof)
- * this->shape_value_component
- ((4 * i + j + 2) * this->degree + l,
- quadrature_points[q_point], 0);
- tmp (3 * (i + 2 * (j + 2 * k)) + 1)
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((4 * i + k) * this->degree + l, dof)
- * this->shape_value_component
- ((4 * i + k) * this->degree + l,
- quadrature_points[q_point], 1);
- tmp (3 * (i + 2 * (j + 2 * k)) + 2)
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((2 * (j + 4) + k) * this->degree + l,
- dof)
- * this->shape_value_component
- ((2 * (j + 4) + k) * this->degree + l,
- quadrature_points[q_point], 2);
-
- for (unsigned int m = 0; m < deg; ++m)
- {
- tmp (3 * (i + 2 * (j + 2 * k)))
- -= this->restriction[index][2 * (2 * i + j) + k]
- (((2 * j + 5) * deg + m)
- * this->degree + l + n_edge_dofs,
- dof)
- * this->shape_value_component
- (((2 * j + 5) * deg + m)
- * this->degree + l + n_edge_dofs,
- quadrature_points[q_point], 0);
- tmp (3 * (i + 2 * (j + 2 * k)))
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((2 * (i + 4) * this->degree + l)
- * deg + m + n_edge_dofs, dof)
- * this->shape_value_component
- ((2 * (i + 4) * this->degree + l)
- * deg + m + n_edge_dofs,
- quadrature_points[q_point], 0);
- tmp (3 * (i + 2 * (j + 2 * k)) + 1)
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((2 * k * this->degree + l) * deg + m
- + n_edge_dofs,
- dof)
- * this->shape_value_component
- ((2 * k * this->degree + l) * deg + m
- + n_edge_dofs,
- quadrature_points[q_point], 1);
- tmp (3 * (i + 2 * (j + 2 * k)) + 1)
- -= this->restriction[index][2 * (2 * i + j) + k]
- (((2 * i + 9) * deg + m)
- * this->degree + l + n_edge_dofs,
- dof)
- * this->shape_value_component
- (((2 * i + 9) * deg + m)
- * this->degree + l + n_edge_dofs,
- quadrature_points[q_point], 1);
- tmp (3 * (i + 2 * (j + 2 * k)) + 2)
- -= this->restriction[index][2 * (2 * i + j) + k]
- (((2 * k + 1) * deg + m)
- * this->degree + l + n_edge_dofs,
- dof)
- * this->shape_value_component
- (((2 * k + 1) * deg + m)
- * this->degree + l + n_edge_dofs,
- quadrature_points[q_point], 2);
- tmp (3 * (i + 2 * (j + 2 * k)) + 2)
- -= this->restriction[index][2 * (2 * i + j) + k]
- ((2 * (j + 2) * this->degree + l)
- * deg + m + n_edge_dofs, dof)
- * this->shape_value_component
- ((2 * (j + 2) * this->degree + l)
- * deg + m + n_edge_dofs,
- quadrature_points[q_point], 2);
- }
- }
-
- tmp *= quadrature.weight (q_point);
-
- for (unsigned int i = 0; i <= deg; ++i)
+ }
+
+ const double weight = 2.0 * quadrature.weight (q_point);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
{
const double L_i_0
- = legendre_polynomials[i].value
- (quadrature_points[q_point] (0));
+ = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0));
const double L_i_1
- = legendre_polynomials[i].value
- (quadrature_points[q_point] (1));
+ = weight * legendre_polynomials[i].value (quadrature_points[q_point] (1));
const double L_i_2
- = legendre_polynomials[i].value
- (quadrature_points[q_point] (2));
-
+ = weight * legendre_polynomials[i].value (quadrature_points[q_point] (2));
+
for (unsigned int j = 0; j < deg; ++j)
{
- const double l_j_0
- = L_i_0 * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (1));
- const double l_j_1
- = L_i_1 * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (0));
- const double l_j_2
- = L_i_2 * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (0));
-
+ const double Le_j_0
+ = legendre_polynomials[j + 1].value (quadrature_points[q_point] (0));
+ const double Le_j_1
+ = legendre_polynomials[j + 1].value (quadrature_points[q_point] (1));
+ const double lo_j_0
+ = lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0));
+ const double lo_j_1
+ = lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1));
+
for (unsigned int k = 0; k < deg; ++k)
{
- const double l_k_0
- = l_j_0 * lobatto_polynomials[k + 2].value
- (quadrature_points[q_point] (2));
+ const double L_k
+ = legendre_polynomials[k + 1].value (quadrature_points[q_point] (2));
const double l_k_1
- = l_j_1 * lobatto_polynomials[k + 2].value
- (quadrature_points[q_point] (2));
+ = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (1));
const double l_k_2
- = l_j_2 * lobatto_polynomials[k + 2].value
- (quadrature_points[q_point] (1));
-
- for (unsigned int l = 0; l < 8; ++l)
- {
- system_rhs ((i * deg + j) * deg + k,
- 3 * l)
- += tmp (3 * l) * l_k_0;
- system_rhs ((i * deg + j) * deg + k,
- 3 * l + 1)
- += tmp (3 * l + 1) * l_k_1;
- system_rhs ((i * deg + j) * deg + k,
- 3 * l + 2)
- += tmp (3 * l + 2) * l_k_2;
- }
+ = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (2));
+
+ this->restriction[index][2 * (2 * l + m) + n]
+ ((i * deg + j) * deg + k + n_boundary_dofs, dof)
+ += L_i_0 * (this->shape_grad_component (dof, quadrature_point, 0)[2]
+ * lo_j_1 * L_k
+ + this->shape_grad_component (dof, quadrature_point, 0)[1]
+ * Le_j_1 * l_k_2
+ + this->shape_value_component (dof, quadrature_point, 0)
+ * lo_j_1 * l_k_2);
+ this->restriction[index][2 * (2 * l + m) + n]
+ ((i + (j + deg) * this->degree) * deg + k + n_boundary_dofs, dof)
+ += L_i_1 * (this->shape_grad_component (dof, quadrature_point, 1)[2]
+ * lo_j_0 * L_k
+ + this->shape_grad_component (dof, quadrature_point, 1)[0]
+ * Le_j_0 * l_k_2
+ + this->shape_value_component (dof, quadrature_point, 1)
+ * lo_j_0 * l_k_2);
+ this->restriction[index][2 * (2 * l + m) + n]
+ (i + (k + (j + 2 * deg) * deg) * this->degree + n_boundary_dofs, dof)
+ += L_i_2 * (this->shape_grad_component (dof, quadrature_point, 2)[1]
+ * lo_j_0
+ * legendre_polynomials[k + 1].value (quadrature_points[q_point] (1))
+ + this->shape_grad_component (dof, quadrature_point, 2)[0]
+ * Le_j_0 * l_k_1
+ + this->shape_value_component (dof, quadrature_point, 2)
+ * lo_j_0 * l_k_1);
}
}
}
}
-
- system_matrix_inv.mmult (solution, system_rhs);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- for (unsigned int l = 0; l <= deg; ++l)
- for (unsigned int m = 0; m < deg; ++m)
- for (unsigned int n = 0; n < deg; ++n)
- {
- if (std::abs (solution
- ((l * deg + m) * deg + n,
- 3 * (i + 2 * (j + 2 * k))))
- > 1e-14)
- this->restriction[index][2 * (2 * i + j) + k]
- ((l * deg + m) * deg + n + n_boundary_dofs,
- dof) = solution ((l * deg + m) * deg + n,
- 3 * (i + 2 * (j + 2 * k)));
-
- if (std::abs (solution
- ((l * deg + m) * deg + n,
- 3 * (i + 2 * (j + 2 * k)) + 1))
- > 1e-14)
- this->restriction[index][2 * (2 * i + j) + k]
- ((l + (m + deg) * this->degree) * deg + n
- + n_boundary_dofs,
- dof) = solution ((l * deg + m) * deg + n,
- 3 * (i + 2 * (j + 2 * k)) + 1);
-
- if (std::abs (solution
- ((l * deg + m) * deg + n,
- 3 * (i + 2 * (j + 2 * k)) + 2))
- > 1e-14)
- this->restriction[index][2 * (2 * i + j) + k]
- (l + ((m + 2 * deg) * deg + n) * this->degree
- + n_boundary_dofs, dof)
- = solution ((l * deg + m) * deg + n,
- 3 * (i + 2 * (j + 2 * k)) + 2);
- }
}
}
return dpo;
}
+
//---------------------------------------------------------------------------
// Data field initialization
//---------------------------------------------------------------------------
+ GeometryInfo<dim>::lines_per_face)
= 1.0;
interpolation_matrix
- (i + (j + GeometryInfo<dim>::lines_per_face) * source_fe.degree,
- i + (j + GeometryInfo<dim>::lines_per_face) * this->degree)
+ (i + (j + GeometryInfo<dim>::lines_per_face + source_fe.degree - 1)
+ * source_fe.degree,
+ i + (j + GeometryInfo<dim>::lines_per_face + this->degree-1) * this->degree)
= 1.0;
}
}
// In this function we compute the
// subface interpolation matrix.
- // This is done by a projection-
- // based interpolation. Therefore
- // we first interpolate the
- // shape functions of the higher
- // order element on the lowest
- // order edge shape functions.
- // Then the remaining part of
- // the interpolated shape
- // functions is projected on the
- // higher order edge shape
- // functions, the face shape
- // functions and the interior
- // shape functions (if they all
- // exist).
template <int dim>
void
FE_Nedelec<dim>::get_subface_interpolation_matrix(
const unsigned int subface,
FullMatrix<double>& interpolation_matrix) const
{
- // this is only implemented, if the
- // source FE is also a
- // Nedelec element
+ // this is only implemented, if the
+ // source FE is also a
+ // Nedelec element
typedef FE_Nedelec<dim> FEN;
typedef FiniteElement<dim> FEL;
ExcDimensionMismatch (interpolation_matrix.n (),
this->dofs_per_face));
- // ok, source is a Nedelec element, so
- // we will be able to do the work
+ // ok, source is a Nedelec element, so
+ // we will be able to do the work
const FE_Nedelec<dim> &source_fe
= dynamic_cast<const FE_Nedelec<dim>&> (source);
- // Make sure, that the element,
+ // Make sure, that the element,
// for which the DoFs should be
// constrained is the one with
// the higher polynomial degree.
// also if this assertion is not
// satisfied. But the matrices
// produced in that case might
- // lead to problems in the
+ // lead to problems in the
// hp procedures, which use this
- // method.
+ // method.
Assert (this->dofs_per_face <= source_fe.dofs_per_face,
typename FEL::ExcInterpolationNotImplemented ());
interpolation_matrix = 0.0;
- // Perform projection-based interpolation
- // as usual.
+ // Set the degrees of freedom
+ // by interpolation as usual.
const QGauss<1> edge_quadrature (source_fe.degree);
const std::vector<Point<1> >&
edge_quadrature_points = edge_quadrature.get_points ();
+ const std::vector<Polynomials::Polynomial<double> >&
+ legendre_polynomials
+ = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
const unsigned int& n_edge_quadrature_points = edge_quadrature.size ();
switch (dim)
{
case 2:
{
- for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
- for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
- ++q_point)
- {
- const Point<dim> quadrature_point (0.0,
- 0.5 * (edge_quadrature_points[q_point] (0)
- + subface));
-
- interpolation_matrix (0, dof) += 0.5
- * edge_quadrature.weight (q_point)
- * this->shape_value_component
- (dof, quadrature_point, 1);
- }
-
- if (source_fe.degree > 1)
+ for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
+ ++q_point)
{
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
- FullMatrix<double> system_matrix_inv (source_fe.degree - 1,
- source_fe.degree - 1);
-
- {
- FullMatrix<double> assembling_matrix (source_fe.degree - 1,
- n_edge_quadrature_points);
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
- {
- const double weight
- = std::sqrt (edge_quadrature.weight (q_point));
-
- for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
- assembling_matrix (i, q_point) = weight
- * legendre_polynomials[i + 1].value
- (edge_quadrature_points[q_point] (0));
- }
-
- FullMatrix<double> system_matrix (source_fe.degree - 1, source_fe.degree - 1);
-
- assembling_matrix.mTmult (system_matrix, assembling_matrix);
- system_matrix_inv.invert (system_matrix);
- }
-
- Vector<double> solution (source_fe.degree - 1);
- Vector<double> system_rhs (source_fe.degree - 1);
-
+ const Point<dim> quadrature_point (0.0,
+ 0.5 * (edge_quadrature_points[q_point] (0)
+ + subface));
+ const double weight = 0.5 * edge_quadrature.weight (q_point);
+
for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
{
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
- {
- const Point<dim> quadrature_point_0 (0.0,
- 0.5 * (edge_quadrature_points[q_point] (0)
- + subface));
- const Point<dim> quadrature_point_1 (0.0,
- edge_quadrature_points[q_point] (0));
- const double tmp = edge_quadrature.weight (q_point)
- * (0.5 * this->shape_value_component
- (dof, quadrature_point_0, 1)
- - interpolation_matrix (0,
- dof)
- * source_fe.shape_value_component
- (0, quadrature_point_1, 1));
-
- for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
- system_rhs (i) += tmp
- * legendre_polynomials[i + 1].value
- (edge_quadrature_points[q_point] (0));
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
- if (std::abs (solution (i)) > 1e-14)
- interpolation_matrix (i + 1, dof) = solution (i);
+ const double shape_value = weight * this->shape_value_component (dof, quadrature_point, 1);
+
+ for (unsigned int i = 0; i < source_fe.dofs_per_face; ++i)
+ interpolation_matrix (i, dof) += shape_value
+ * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
}
}
const double shifts[4][2] = { { 0.0, 0.0 }, { 1.0, 0.0 },
{ 0.0, 1.0 }, { 1.0, 1.0 } };
- for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
- for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
- ++q_point)
- {
- const double weight = 0.5 * edge_quadrature.weight (q_point);
-
- for (unsigned int i = 0; i < 2; ++i)
- {
- Point<dim>
- quadrature_point (0.5 * (i + shifts[subface][0]),
- 0.5 * (edge_quadrature_points[q_point] (0)
- + shifts[subface][1]),
- 0.0);
-
- interpolation_matrix (i * source_fe.degree, dof) += weight
- * this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- quadrature_point,
- 1);
- quadrature_point
- = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
- + shifts[subface][0]),
- 0.5 * (i + shifts[subface][1]), 0.0);
- interpolation_matrix ((i + 2) * source_fe.degree, dof)
- += weight * this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- quadrature_point, 0);
- }
- }
-
- if (source_fe.degree > 1)
+ for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
+ ++q_point)
{
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
- FullMatrix<double> system_matrix_inv (source_fe.degree - 1,
- source_fe.degree - 1);
-
- {
- FullMatrix<double> assembling_matrix (source_fe.degree - 1,
- n_edge_quadrature_points);
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
- {
- const double weight
- = std::sqrt (edge_quadrature.weight (q_point));
-
- for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
- assembling_matrix (i, q_point) = weight
- * legendre_polynomials[i + 1].value
- (edge_quadrature_points[q_point] (0));
- }
-
- FullMatrix<double> system_matrix (source_fe.degree - 1, source_fe.degree - 1);
+ const double weight = 0.5 * edge_quadrature.weight (q_point);
- assembling_matrix.mTmult (system_matrix, assembling_matrix);
- system_matrix_inv.invert (system_matrix);
- }
-
- FullMatrix<double> solution (source_fe.degree - 1,
- GeometryInfo<dim>::lines_per_face);
- FullMatrix<double> system_rhs (source_fe.degree - 1,
- GeometryInfo<dim>::lines_per_face);
- Vector<double> tmp (GeometryInfo<dim>::lines_per_face);
-
- for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
+ for (unsigned int i = 0; i < 2; ++i)
{
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_edge_quadrature_points; ++q_point)
+ Point<dim>
+ quadrature_point (0.5 * (i + shifts[subface][0]),
+ 0.5 * (edge_quadrature_points[q_point] (0)
+ + shifts[subface][1]),
+ 0.0);
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
{
- const double weight = edge_quadrature.weight (q_point);
-
- for (unsigned int i = 0; i < 2; ++i)
- {
- Point<dim>
- quadrature_point_0
- (0.5 * (i + shifts[subface][0]),
- 0.5 * (edge_quadrature_points[q_point] (0)
- + shifts[subface][1]), 0.0);
- Point<dim> quadrature_point_1 (i,
- edge_quadrature_points[q_point] (0),
- 0.0);
-
- tmp (i) = weight
- * (0.5 * this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- quadrature_point_0, 1)
- - interpolation_matrix
- (i * source_fe.degree, dof)
- * source_fe.shape_value_component
- (i * source_fe.degree,
- quadrature_point_1, 1));
- quadrature_point_0
- = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
- + shifts[subface][0]),
- 0.5 * (i + shifts[subface][1]),
- 0.0);
- quadrature_point_1
- = Point<dim> (edge_quadrature_points[q_point] (0),
- i, 0.0);
- tmp (i + 2) = weight
- * (0.5 * this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- quadrature_point_0, 0)
- - interpolation_matrix
- ((i + 2) * source_fe.degree,
- dof)
- * source_fe.shape_value_component
- ((i + 2) * source_fe.degree,
- quadrature_point_1, 0));
- }
-
- for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
- {
- const double L_i
- = legendre_polynomials[i + 1].value
- (edge_quadrature_points[q_point] (0));
-
- for (unsigned int j = 0;
- j < GeometryInfo<dim>::lines_per_face; ++j)
- system_rhs (i, j) += tmp (j) * L_i;
- }
+ const double shape_value
+ = weight * this->shape_value_component (this->face_to_cell_index (dof, 4),
+ quadrature_point, 1);
+
+ for (unsigned int j = 0; j < source_fe.degree; ++j)
+ interpolation_matrix (i * source_fe.degree + j, dof)
+ += shape_value
+ * legendre_polynomials[j].value (edge_quadrature_points[q_point] (0));
+ }
+
+ quadrature_point
+ = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
+ + shifts[subface][0]),
+ 0.5 * (i + shifts[subface][1]), 0.0);
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
+ {
+ const double shape_value
+ = weight * this->shape_value_component (this->face_to_cell_index (dof, 4),
+ quadrature_point, 0);
+
+ for (unsigned int j = 0; j < source_fe.degree; ++j)
+ interpolation_matrix ((i + 2) * source_fe.degree + j, dof)
+ += shape_value
+ * legendre_polynomials[j].value (edge_quadrature_points[q_point] (0));
}
-
- system_matrix_inv.mmult (solution, system_rhs);
-
- for (unsigned int i = 0;
- i < GeometryInfo<dim>::lines_per_face; ++i)
- for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
- if (std::abs (solution (j, i)) > 1e-14)
- interpolation_matrix (i * source_fe.degree + j + 1,
- dof) = solution (j, i);
}
+ }
+ if (source_fe.degree > 1)
+ {
const QGauss<2> quadrature (source_fe.degree);
const std::vector<Point<2> >&
quadrature_points = quadrature.get_points ();
lobatto_polynomials
= Polynomials::Lobatto::generate_complete_basis
(source_fe.degree);
- const unsigned int n_boundary_dofs
+ const unsigned int n_edge_dofs
= GeometryInfo<dim>::lines_per_face * source_fe.degree;
const unsigned int& n_quadrature_points = quadrature.size ();
- {
- FullMatrix<double>
- assembling_matrix (source_fe.degree * (source_fe.degree - 1),
- n_quadrature_points);
-
- for (unsigned int q_point = 0; q_point < n_quadrature_points;
- ++q_point)
- {
- const double weight = std::sqrt (quadrature.weight (q_point));
-
- for (unsigned int i = 0; i < source_fe.degree; ++i)
- {
- const double L_i = weight
- * legendre_polynomials[i].value
- (quadrature_points[q_point] (0));
-
- for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
- assembling_matrix (i * (source_fe.degree - 1) + j,
- q_point)
- = L_i * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (1));
- }
- }
-
- FullMatrix<double> system_matrix (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 (), 2);
- system_rhs.reinit (system_matrix_inv.m (), 2);
- tmp.reinit (2);
-
- for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
+ for (unsigned int q_point = 0; q_point < n_quadrature_points;
+ ++q_point)
{
- system_rhs = 0.0;
-
- for (unsigned int q_point = 0;
- q_point < n_quadrature_points; ++q_point)
+ const double weight = 0.5 * quadrature.weight (q_point);
+ const Point<dim>
+ quadrature_point
+ (0.5 * (quadrature_points[q_point] (0)
+ + shifts[subface][0]),
+ 0.5 * (quadrature_points[q_point] (1)
+ + shifts[subface][1]), 0.0);
+
+ for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
{
- Point<dim>
- quadrature_point
- (0.5 * (quadrature_points[q_point] (0)
- + shifts[subface][0]),
- 0.5 * (quadrature_points[q_point] (1)
- + shifts[subface][1]), 0.0);
- tmp (0) = 0.5 * this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- quadrature_point, 0);
- tmp (1) = 0.5 * this->shape_value_component
- (this->face_to_cell_index (dof, 4),
- quadrature_point, 1);
- quadrature_point
- = Point<dim> (quadrature_points[q_point] (0),
- quadrature_points[q_point] (1), 0.0);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < source_fe.degree; ++j)
- {
- tmp (0) -= interpolation_matrix
- ((i + 2) * source_fe.degree + j, dof)
- * source_fe.shape_value_component
- ((i + 2) * source_fe.degree + j,
- quadrature_point, 0);
- tmp (1) -= interpolation_matrix
- (i * source_fe.degree + j, dof)
- * source_fe.shape_value_component
- (i * source_fe.degree + j,
- quadrature_point, 1);
- }
-
- tmp *= quadrature.weight (q_point);
-
+ const double shape_grad_0
+ = this->shape_grad_component (this->face_to_cell_index (dof, 4),
+ quadrature_point, 0)[1];
+ const double shape_grad_1
+ = this->shape_grad_component (this->face_to_cell_index (dof, 4),
+ quadrature_point, 1)[0];
+ const double shape_value_0
+ = this->shape_value_component (this->face_to_cell_index (dof, 4),
+ quadrature_point, 0);
+ const double shape_value_1
+ = this->shape_value_component (this->face_to_cell_index (dof, 4),
+ quadrature_point, 1);
+
for (unsigned int i = 0; i < source_fe.degree; ++i)
{
- const double L_i_0 = legendre_polynomials[i].value
- (quadrature_points[q_point] (0));
- const double L_i_1 = legendre_polynomials[i].value
- (quadrature_points[q_point] (1));
-
+ const double L_i_0
+ = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0));
+ const double L_i_1
+ = weight * legendre_polynomials[i].value (quadrature_points[q_point] (1));
+
for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
{
- system_rhs (i * (source_fe.degree - 1) + j, 0)
- += tmp (0) * L_i_0
- * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (1));
- system_rhs (i * (source_fe.degree - 1) + j, 1)
- += tmp (1) * L_i_1
- * lobatto_polynomials[j + 2].value
- (quadrature_points[q_point] (0));
+ interpolation_matrix (i * (source_fe.degree - 1) + j + n_edge_dofs, dof)
+ += L_i_0 * (shape_grad_0
+ * legendre_polynomials[j + 1].value (quadrature_points[q_point] (1))
+ + shape_value_0
+ * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1)));
+ interpolation_matrix (i + (j + source_fe.degree - 1) * source_fe.degree
+ + n_edge_dofs, dof)
+ += L_i_1 * (shape_grad_1
+ * legendre_polynomials[j + 1].value (quadrature_points[q_point] (0))
+ + shape_value_1
+ * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0)));
}
}
}
-
- system_matrix_inv.mmult (solution, system_rhs);
-
- for (unsigned int i = 0; i < source_fe.degree; ++i)
- for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
- {
- if (std::abs (solution (i * (source_fe.degree - 1) + j, 0))
- > 1e-14)
- interpolation_matrix (i * (source_fe.degree - 1)
- + j + n_boundary_dofs, dof)
- = solution (i * (source_fe.degree - 1) + j, 0);
-
- if (std::abs (solution (i * (source_fe.degree - 1) + j, 1))
- > 1e-14)
- interpolation_matrix (i + (j + source_fe.degree - 1)
- * source_fe.degree
- + n_boundary_dofs, dof)
- = solution (i * (source_fe.degree - 1) + j, 1);
- }
}
}
// its values at the generalized support
// points in the finite element space on the
// reference cell.
- // This is done as usual by projection-based
- // interpolation.
+ // The interpolation on the edge degrees of
+ // freedom is done by direct calculation.
+ // For the interpolation on the remaining
+ // degrees of freedom we use a projection-
+ // based interpolation scheme. Therefore
+ // the remaining part of the interpolated
+ // function is projected on the face shape
+ // functions and the interior shape
+ // functions (if they exist).
template <int dim>
void
FE_Nedelec<dim>::interpolate (std::vector<double>& local_dofs,
const std::vector<Vector<double> >& values,
unsigned int offset) const
{
- const unsigned int deg = this->degree-1;
-
Assert (values.size () == this->generalized_support_points.size (),
ExcDimensionMismatch (values.size (),
this->generalized_support_points.size ()));
std::fill (local_dofs.begin (), local_dofs.end (), 0.);
if (offset < dim)
- switch (dim)
- {
- case 2:
- {
- const QGauss<1> reference_edge_quadrature (this->degree);
- const unsigned int& n_edge_points
- = reference_edge_quadrature.size ();
-
+ {
+ const unsigned int deg = this->degree-1;
+ const std::vector<Polynomials::Polynomial<double> >&
+ legendre_polynomials
+ = Polynomials::Legendre::generate_complete_basis (deg);
+ const QGauss<1> edge_quadrature (this->degree);
+ const std::vector<Point<1> >&
+ edge_quadrature_points = edge_quadrature.get_points ();
+ const unsigned int& n_edge_points = edge_quadrature.size ();
+
+ switch (dim)
+ {
+ case 2:
+ {
// Let us begin with the
- // interpolation part.
- for (unsigned int i = 0; i < 2; ++i)
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[i * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[q_point + i * n_edge_points] (1);
-
+ // edge degrees of freedom.
+ for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+ {
+ const double weight = edge_quadrature.weight (q_point);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ {
+ const double L_i
+ = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+
+ for (unsigned int j = 0; j < 2; ++j)
+ {
+ local_dofs[i + j * this->degree]
+ += L_i * values[q_point + j * n_edge_points] (1);
+
+ if (offset == 0)
+ local_dofs[i + (j + 2) * this->degree]
+ += L_i * values[q_point + (j + 2) * n_edge_points] (0);
+ }
+ }
+ }
+
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- if (std::abs (local_dofs[i * this->degree]) < 1e-14)
- local_dofs[i * this->degree] = 0.0;
- }
+ for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+ for (unsigned int i = 0; i < this->degree; ++i)
+ if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+ local_dofs[i + e * this->degree] = 0.0;
- if (offset == 0)
- for (unsigned int i = 0; i < 2; ++i)
+ // If the degree is greater
+ // than 0, then we have still
+ // some interpolations onto
+ // the interior shape
+ // functions left.
+ if (this->degree > 1)
{
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[(i + 2) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[q_point + (i + 2) * n_edge_points] (0);
+ const std::vector<Polynomials::Polynomial<double> >&
+ lobatto_polynomials
+ = Polynomials::Lobatto::generate_complete_basis
+ (this->degree);
+ std::vector<Polynomials::Polynomial<double> >
+ lobatto_polynomials_grad (this->degree);
+
+ for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
+ ++i)
+ lobatto_polynomials_grad[i]
+ = lobatto_polynomials[i + 1].derivative ();
+
+ // We set up the system
+ // matrix and use it for
+ // both, the horizontal
+ // and the vertical
+ // interior shape
+ // functions.
+ const QGauss<dim> reference_quadrature (this->degree);
+ const unsigned int& n_interior_points
+ = reference_quadrature.size ();
- if (std::abs (local_dofs[(i + 2) * this->degree]) < 1e-14)
- local_dofs[(i + 2) * this->degree] = 0.0;
- }
+ FullMatrix<double> system_matrix (deg * this->degree,
+ deg * this->degree);
- // If the degree is greater
- // than 0, then we have still
- // some higher order edge
- // shape functions to
- // consider.
- // Here the projection part
- // starts. The dof values
- // are obtained by solving
- // a linear system of
- // equations.
- if (this->degree > 1)
- {
- // We start with projection
- // on the higher order edge
- // shape function.
- const std::vector<Polynomials::Polynomial<double> >&
- lobatto_polynomials
- = Polynomials::Lobatto::generate_complete_basis
- (this->degree);
- const unsigned int
- line_coordinate[GeometryInfo<2>::lines_per_cell]
- = {1, 1, 0, 0};
- std::vector<Polynomials::Polynomial<double> >
- lobatto_polynomials_grad (this->degree);
-
- for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
- ++i)
- lobatto_polynomials_grad[i]
- = lobatto_polynomials[i + 1].derivative ();
-
- // Set up the system matrix.
- // This can be used for all
- // edges.
- FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
-
- for (unsigned int i = 0; i < system_matrix.m (); ++i)
- for (unsigned int j = 0; j < system_matrix.n (); ++j)
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- system_matrix (i, j)
- += boundary_weights (q_point, j)
- * lobatto_polynomials_grad[i + 1].value
- (this->generalized_face_support_points[q_point]
- (1));
-
- FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < this->degree; ++k)
+ for (unsigned int l = 0; l < deg; ++l)
+ for (unsigned int q_point = 0;
+ q_point < n_interior_points; ++q_point)
+ system_matrix (i * deg + j, k * deg + l)
+ += reference_quadrature.weight (q_point)
+ * legendre_polynomials[i].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points]
+ (0))
+ * lobatto_polynomials[j + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points]
+ (1))
+ * lobatto_polynomials_grad[k].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points]
+ (0))
+ * lobatto_polynomials[l + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points]
+ (1));
- system_matrix_inv.invert (system_matrix);
+ FullMatrix<double> system_matrix_inv (system_matrix.m (),
+ system_matrix.m ());
+
+ system_matrix_inv.invert (system_matrix);
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
+ Vector<double> solution (system_matrix_inv.m ());
+ Vector<double> system_rhs (system_matrix.m ());
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell; ++line)
- if ((line < 2) || (offset == 0))
+ if (offset == 0)
{
- // Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
+ // Set up the right hand side
+ // for the horizontal shape
+ // functions.
+ for (unsigned int q_point = 0;
+ q_point < n_interior_points; ++q_point)
{
- const double tmp
- = values[line * n_edge_points + q_point]
- (line_coordinate[line])
- - local_dofs[line * this->degree]
- * this->shape_value_component
- (line * this->degree,
- this->generalized_support_points[line
- * n_edge_points
- + q_point],
- line_coordinate[line]);
-
- for (unsigned int i = 0; i < system_rhs.size ();
- ++i)
- system_rhs (i) += boundary_weights (q_point, i)
- * tmp;
+ double tmp
+ = values[q_point + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points] (0);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(i + 2) * this->degree + j]
+ * this->shape_value_component
+ ((i + 2) * this->degree + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points],
+ 0);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < this->degree-1; ++j)
+ system_rhs (i * (this->degree-1) + j)
+ += reference_quadrature.weight (q_point) * tmp
+ * lobatto_polynomials_grad[i].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points]
+ (0))
+ * lobatto_polynomials[j + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points]
+ (1));
}
system_matrix_inv.vmult (solution, system_rhs);
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i < solution.size (); ++i)
- if (std::abs (solution (i)) > 1e-14)
- local_dofs[line * this->degree + i + 1]
- = solution (i);
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j)) > 1e-14)
+ local_dofs[(i + GeometryInfo<dim>::lines_per_cell)
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution (i * deg + j);
}
- // Then we go on to the
- // interior shape
- // functions. Again we
- // set up the system
- // matrix and use it
- // for both, the
- // horizontal and the
- // vertical, interior
- // shape functions.
- const QGauss<dim> reference_quadrature (this->degree);
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (this->degree-1);
- const unsigned int& n_interior_points
- = reference_quadrature.size ();
-
- system_matrix.reinit ((this->degree-1) * this->degree,
- (this->degree-1) * this->degree);
- system_matrix = 0;
+ // Set up the right hand side
+ // for the vertical shape
+ // functions.
+ system_rhs = 0;
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- for (unsigned int k = 0; k < this->degree; ++k)
- for (unsigned int l = 0; l < this->degree-1; ++l)
- for (unsigned int q_point = 0;
- q_point < n_interior_points; ++q_point)
- system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
- += reference_quadrature.weight (q_point)
- * legendre_polynomials[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (0))
- * lobatto_polynomials[j + 2].value
+ for (unsigned int q_point = 0; q_point < n_interior_points;
+ ++q_point)
+ {
+ double tmp
+ = values[q_point + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points] (1);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[i * this->degree + j]
+ * this->shape_value_component
+ (i * this->degree + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points],
+ 1);
+
+ for (unsigned i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += reference_quadrature.weight (q_point) * tmp
+ * lobatto_polynomials_grad[i].value
(this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points]
(1))
- * lobatto_polynomials_grad[k].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (0))
- * lobatto_polynomials[l + 2].value
+ * lobatto_polynomials[j + 2].value
(this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points]
- (1));
-
- 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.m ());
-
- if (offset == 0)
- {
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_interior_points; ++q_point)
- {
- double tmp
- = values[q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points] (0);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- tmp -= local_dofs[(i + 2) * this->degree + j]
- * this->shape_value_component
- ((i + 2) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- 0);
-
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- system_rhs (i * (this->degree-1) + j)
- += reference_quadrature.weight (q_point) * tmp
- * lobatto_polynomials_grad[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (0))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (1));
- }
+ (0));
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell)
- * (this->degree-1) + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * (this->degree-1) + j);
- }
-
- // Set up the right hand side
- // for the vertical shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points] (1);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- tmp -= local_dofs[i * this->degree + j]
- * this->shape_value_component
- (i * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- 1);
-
- for (unsigned i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- system_rhs (i * (this->degree-1) + j)
- += reference_quadrature.weight (q_point) * tmp
- * lobatto_polynomials_grad[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (1))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (0));
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j)) > 1e-14)
+ local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+ + deg) * this->degree]
+ = solution (i * deg + j);
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ break;
+ }
+ case 3:
+ {
+ // Let us begin with the
+ // edge degrees of freedom.
+ for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+ {
+ const double weight = edge_quadrature.weight (q_point);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ {
+ const double L_i
+ = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+
+ for (unsigned int j = 0; j < 4; ++j)
+ local_dofs[i + (j + 8) * this->degree]
+ += L_i * values[q_point + (j + 8) * n_edge_points] (2);
+
+ if (offset + 1 < dim)
+ {
+ for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int k = 0; k < 2; ++k)
+ local_dofs[i + (j + 4 * k) * this->degree]
+ += L_i * values[q_point + (j + 4 * k) * n_edge_points] (1);
+
+ if (offset == 0)
+ for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int k = 0; k < 2; ++k)
+ local_dofs[i + (j + 4 * k + 2) * this->degree]
+ += L_i * values[q_point + (j + 4 * k + 2) * n_edge_points] (0);
+ }
+ }
+ }
+
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
+ for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
- local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + this->degree-1) * this->degree]
- = solution (i * (this->degree-1) + j);
- }
-
- break;
- }
+ if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+ local_dofs[i + e * this->degree] = 0.0;
- case 3:
- {
- const QGauss<1>
- reference_edge_quadrature (this->degree);
- const unsigned int&
- n_edge_points = reference_edge_quadrature.size ();
+ // If the degree is greater
+ // than 0, then we have still
+ // some interpolation to the
+ // face and interior shape
+ // functions left.
+ if (this->degree > 1)
+ {
+ const std::vector<Polynomials::Polynomial<double> >&
+ lobatto_polynomials
+ = Polynomials::Lobatto::generate_complete_basis
+ (this->degree);
+ FullMatrix<double> system_matrix (deg * this->degree,
+ deg * this->degree);
+ std::vector<Polynomials::Polynomial<double> >
+ lobatto_polynomials_grad (this->degree);
+
+ for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
+ ++i)
+ lobatto_polynomials_grad[i]
+ = lobatto_polynomials[i + 1].derivative ();
+
+ // We set up the system
+ // matrix and use it for
+ // both, the horizontal
+ // and the vertical, shape
+ // functions.
+ const unsigned int
+ n_face_points = n_edge_points * n_edge_points;
- // Let us begin with the
- // interpolation part.
- for (unsigned int i = 0; i < 4; ++i)
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[(i + 8) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[q_point + (i + 8) * n_edge_points] (2);
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < this->degree; ++k)
+ for (unsigned int l = 0; l < deg; ++l)
+ for (unsigned int q_point = 0; q_point < n_face_points;
+ ++q_point)
+ system_matrix (i * deg + j, k * deg + l)
+ += boundary_weights (q_point + n_edge_points,
+ 2 * (k * deg + l))
+ * legendre_polynomials[i].value
+ (this->generalized_face_support_points[q_point
+ + 4
+ * n_edge_points]
+ (0))
+ * lobatto_polynomials[j + 2].value
+ (this->generalized_face_support_points[q_point
+ + 4
+ * n_edge_points]
+ (1));
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
- local_dofs[(i + 8) * this->degree] = 0.0;
- }
+ FullMatrix<double> system_matrix_inv (system_matrix.m (),
+ system_matrix.m ());
+
+ system_matrix_inv.invert (system_matrix);
+
+ Vector<double> solution (system_matrix.m ());
+ Vector<double> system_rhs (system_matrix.m ());
- if (offset + 1 < dim)
- {
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int face = 0;
+ face < GeometryInfo<dim>::faces_per_cell; ++face)
{
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[(i + 4 * j) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[q_point + (i + 4 * j) * n_edge_points]
- (1);
+ switch (face)
+ {
+ case 0:
+ {
+ if (offset + 1 < dim)
+ {
+ // Set up the right hand side
+ // for the horizontal shape
+ // functions.
+ system_rhs = 0;
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- if (std::abs (local_dofs[(i + 4 * j) * this->degree])
- < 1e-14)
- local_dofs[(i + 4 * j) * this->degree] = 0.0;
- }
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points] (1);
- if (offset == 0)
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- {
- for (unsigned int q_point = 0;
- q_point < n_edge_points; ++q_point)
- local_dofs[(i + 4 * j + 2) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[q_point + (i + 4 * j + 2)
- * n_edge_points] (0);
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp
+ -= local_dofs[4 * i * this->degree
+ + j]
+ * this->shape_value_component
+ (4 * i * this->degree + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points],
+ 1);
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- if (std::abs (local_dofs[(i + 4 * j + 2)
- * this->degree]) < 1e-14)
- local_dofs[(i + 4 * j + 2) * this->degree] = 0.0;
- }
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j)) * tmp;
+ }
- // If the degree is greater
- // than 0, then we have still
- // some higher order shape
- // functions to consider.
- // Here the projection part
- // starts. The dof values
- // are obtained by solving
- // a linear system of
- // equations.
- if (this->degree > 1)
- {
- // We start with projection
- // on the higher order edge
- // shape function.
- const std::vector<Polynomials::Polynomial<double> >&
- lobatto_polynomials
- = Polynomials::Lobatto::generate_complete_basis
- (this->degree);
- const unsigned int
- line_coordinate[GeometryInfo<3>::lines_per_cell]
- = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
- FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
- FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
- std::vector<Polynomials::Polynomial<double> >
- lobatto_polynomials_grad (this->degree);
-
- for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
- ++i)
- lobatto_polynomials_grad[i]
- = lobatto_polynomials[i + 1].derivative ();
-
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- // Set up the system matrix.
- // This can be used for all
- // edges.
- for (unsigned int i = 0; i < system_matrix.m (); ++i)
- for (unsigned int j = 0; j < system_matrix.n (); ++j)
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- system_matrix (i, j)
- += boundary_weights (q_point, j)
- * lobatto_polynomials_grad[i + 1].value
- (this->generalized_face_support_points[q_point]
- (1));
+ system_matrix_inv.vmult (solution, system_rhs);
- system_matrix_inv.invert (system_matrix);
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j))
+ > 1e-14)
+ local_dofs[(i
+ + GeometryInfo<dim>::lines_per_cell)
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution (i * deg + j);
+ }
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell; ++line)
- {
- // Set up the right hand side.
- system_rhs = 0;
+ // Set up the right hand side
+ // for the vertical shape
+ // functions.
+ system_rhs = 0;
- if ((((line == 0) || (line == 1) || (line == 4) ||
- (line == 5)) && (offset + 1 < dim)) ||
- (((line == 2) || (line == 3) || (line == 6) ||
- (line == 7)) && (offset == 0)) || (line > 7))
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- {
- double tmp
- = values[line * n_edge_points + q_point]
- (line_coordinate[line])
- - local_dofs[line * this->degree]
- * this->shape_value_component
- (line * this->degree,
- this->generalized_support_points[line
- * this->degree
- + q_point],
- line_coordinate[line]);
-
- for (unsigned int i = 0; i < system_rhs.size ();
- ++i)
- system_rhs (i)
- += boundary_weights (q_point, i) * tmp;
- }
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points] (2);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[2 * (i + 4)
+ * this->degree + j]
+ * this->shape_value_component
+ (2 * (i + 4) * this->degree + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points],
+ 2);
+
+ for (unsigned i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j) + 1)
+ * tmp;
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i < solution.size (); ++i)
- if (std::abs (solution (i)) > 1e-14)
- local_dofs[line * this->degree + i + 1]
- = solution (i);
- }
- }
-
- // Then we go on to the
- // face shape functions.
- // Again we set up the
- // system matrix and
- // use it for both, the
- // horizontal and the
- // vertical, shape
- // functions.
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (this->degree-1);
- const unsigned int
- n_face_points = n_edge_points * n_edge_points;
-
- system_matrix.reinit ((this->degree-1) * this->degree,
- (this->degree-1) * this->degree);
- system_matrix = 0;
-
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- for (unsigned int k = 0; k < this->degree; ++k)
- for (unsigned int l = 0; l < this->degree-1; ++l)
- for (unsigned int q_point = 0; q_point < n_face_points;
- ++q_point)
- system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
- += boundary_weights (q_point + n_edge_points,
- 2 * (k * (this->degree-1) + l))
- * legendre_polynomials[i].value
- (this->generalized_face_support_points[q_point
- + 4
- * n_edge_points]
- (0))
- * lobatto_polynomials[j + 2].value
- (this->generalized_face_support_points[q_point
- + 4
- * n_edge_points]
- (1));
-
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.n ());
- system_matrix_inv.invert (system_matrix);
- solution.reinit (system_matrix.m ());
- system_rhs.reinit (system_matrix.m ());
-
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell; ++face)
- {
- switch (face)
- {
- case 0:
- {
- if (offset + 1 < dim)
- {
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points] (1);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- tmp
- -= local_dofs[4 * i * this->degree
- + j]
- * this->shape_value_component
- (4 * i * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- 1);
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j)) > 1e-14)
+ local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+ + deg)
+ * this->degree]
+ = solution (i * deg + j);
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- system_rhs (i * (this->degree-1) + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * (this->degree-1) + j)) * tmp;
- }
+ break;
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ case 1:
+ {
+ if (offset + 1 < dim)
+ {
+ // Set up the right hand side
+ // for the horizontal shape
+ // functions.
+ system_rhs = 0;
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- if (std::abs (solution (i * (this->degree-1) + j))
- > 1e-14)
- local_dofs[(i
- + GeometryInfo<dim>::lines_per_cell)
- * (this->degree-1) + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * (this->degree-1) + j);
- }
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + n_face_points] (1);
- // Set up the right hand side
- // for the vertical shape
- // functions.
- system_rhs = 0;
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(4 * i + 1)
+ * this->degree + j]
+ * this->shape_value_component
+ ((4 * i + 1) * this->degree
+ + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + n_face_points],
+ 1);
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points] (2);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- tmp -= local_dofs[2 * (i + 4)
- * this->degree + j]
- * this->shape_value_component
- (2 * (i + 4) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- 2);
-
- for (unsigned i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- system_rhs (i * (this->degree-1) + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * (this->degree-1) + j) + 1)
- * tmp;
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j)) * tmp;
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
- local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + this->degree-1)
- * this->degree]
- = solution (i * (this->degree-1) + j);
-
- break;
- }
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j))
+ > 1e-14)
+ local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+ + 2 * this->degree) * deg + j
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution (i * deg + j);
+ }
- case 1:
- {
- if (offset + 1 < dim)
- {
// Set up the right hand side
- // for the horizontal shape
+ // for the vertical shape
// functions.
- system_rhs = 0;
+ system_rhs = 0;
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + n_face_points] (1);
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points + n_face_points]
+ (2);
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(4 * i + 1)
- * this->degree + j]
- * this->shape_value_component
- ((4 * i + 1) * this->degree
- + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + n_face_points],
- 1);
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(2 * (i + 4) + 1)
+ * this->degree + j]
+ * this->shape_value_component
+ ((2 * (i + 4) + 1) * this->degree
+ + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + n_face_points],
+ 2);
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j)) * tmp;
- }
+ for (unsigned i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j) + 1) * tmp;
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 2 * this->degree) * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
-
- // Set up the right hand side
- // for the vertical shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points + n_face_points]
- (2);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(2 * (i + 4) + 1)
- * this->degree + j]
- * this->shape_value_component
- ((2 * (i + 4) + 1) * this->degree
- + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + n_face_points],
- 2);
-
- for (unsigned i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j) + 1) * tmp;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j)) > 1e-14)
+ local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+ + 3 * deg)
+ * this->degree]
+ = solution (i * deg + j);
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j)) > 1e-14)
- local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + 3 * deg)
- * this->degree]
- = solution (i * deg + j);
-
- break;
- }
+ break;
+ }
- case 2:
- {
- if (offset == 0)
- {
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- system_rhs = 0;
+ case 2:
+ {
+ if (offset == 0)
+ {
+ // Set up the right hand side
+ // for the horizontal shape
+ // functions.
+ system_rhs = 0;
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points + 2 * n_face_points]
- (2);
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points + 2 * n_face_points]
+ (2);
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(i + 8) * this->degree
- + j]
- * this->shape_value_component
- ((i + 8) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 2
- * n_face_points],
- 2);
-
- for (unsigned i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j)) * tmp;
- }
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(i + 8) * this->degree
+ + j]
+ * this->shape_value_component
+ ((i + 8) * this->degree + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + 2
+ * n_face_points],
+ 2);
+
+ for (unsigned i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j)) * tmp;
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 4 * this->degree) * deg
- + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j))
+ > 1e-14)
+ local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+ + 4 * this->degree) * deg
+ + j
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution (i * deg + j);
+ }
// Set up the right hand side
// for the vertical shape
// functions.
- system_rhs = 0;
+ system_rhs = 0;
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 2 * n_face_points] (0);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(4 * i + 2)
- * this->degree + j]
- * this->shape_value_component
- ((4 * i + 2) * this->degree
- + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 2
- * n_face_points],
- 0);
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j) + 1) * tmp;
- }
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + 2 * n_face_points] (0);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(4 * i + 2)
+ * this->degree + j]
+ * this->shape_value_component
+ ((4 * i + 2) * this->degree
+ + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + 2
+ * n_face_points],
+ 0);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j) + 1) * tmp;
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j)) > 1e-14)
- local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + 5 * deg) * this->degree]
- = solution (i * deg + j);
-
- break;
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j)) > 1e-14)
+ local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+ + 5 * deg) * this->degree]
+ = solution (i * deg + j);
- case 3:
- {
- if (offset == 0)
- {
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- system_rhs = 0;
+ break;
+ }
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points + 3 * n_face_points]
- (2);
+ case 3:
+ {
+ if (offset == 0)
+ {
+ // Set up the right hand side
+ // for the horizontal shape
+ // functions.
+ system_rhs = 0;
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(i + 10) * this->degree
- + j]
- * this->shape_value_component
- ((i + 10) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 3
- * n_face_points],
- 2);
-
- for (unsigned i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j)) * tmp;
- }
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points + 3 * n_face_points]
+ (2);
- system_matrix_inv.vmult (solution, system_rhs);
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(i + 10) * this->degree
+ + j]
+ * this->shape_value_component
+ ((i + 10) * this->degree + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + 3
+ * n_face_points],
+ 2);
+
+ for (unsigned i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j)) * tmp;
+ }
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 6 * this->degree) * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
+ system_matrix_inv.vmult (solution, system_rhs);
+
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j))
+ > 1e-14)
+ local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+ + 6 * this->degree) * deg + j
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution (i * deg + j);
+ }
// Set up the right hand side
// for the vertical shape
// functions.
- system_rhs = 0;
+ system_rhs = 0;
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points + 3
- * n_face_points] (0);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(4 * i + 3)
- * this->degree + j]
- * this->shape_value_component
- ((4 * i + 3) * this->degree
- + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 3
- * n_face_points],
- 0);
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j) + 1) * tmp;
- }
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points + 3
+ * n_face_points] (0);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(4 * i + 3)
+ * this->degree + j]
+ * this->shape_value_component
+ ((4 * i + 3) * this->degree
+ + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + 3
+ * n_face_points],
+ 0);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j) + 1) * tmp;
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j)) > 1e-14)
- local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + 7 * deg) * this->degree]
- = solution (i * deg + j);
-
- break;
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j)) > 1e-14)
+ local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+ + 7 * deg) * this->degree]
+ = solution (i * deg + j);
- case 4:
- {
+ break;
+ }
+
+ case 4:
+ {
+ if (offset + 1 < dim)
+ {
+ // Set up the right hand side
+ // for the horizontal shape
+ // functions.
+ if (offset == 0)
+ {
+ system_rhs = 0;
+
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points + 4
+ * n_face_points] (0);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(i + 2)
+ * this->degree
+ + j]
+ * this->shape_value_component
+ ((i + 2) * this->degree
+ + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + 4
+ * n_face_points],
+ 0);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j)) * tmp;
+ }
+
+ system_matrix_inv.vmult
+ (solution, system_rhs);
+
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j))
+ > 1e-14)
+ local_dofs[(i + GeometryInfo<dim>::lines_per_cell
+ + 8 * this->degree) * deg
+ + j
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution (i * deg + j);
+ }
+
+ // Set up the right hand side
+ // for the vertical shape
+ // functions.
+ system_rhs = 0;
+
+ for (unsigned int q_point = 0;
+ q_point < n_face_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points + 4
+ * n_face_points] (1);
+
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[i * this->degree + j]
+ * this->shape_value_component
+ (i * this->degree + j,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + 4
+ * n_face_points],
+ 1);
+
+ for (unsigned i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ system_rhs (i * deg + j)
+ += boundary_weights
+ (q_point + n_edge_points,
+ 2 * (i * deg + j) + 1) * tmp;
+ }
+
+ system_matrix_inv.vmult (solution, system_rhs);
+
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ if (std::abs (solution (i * deg + j))
+ > 1e-14)
+ local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
+ + 9 * deg)
+ * this->degree]
+ = solution (i * deg + j);
+ }
+
+ break;
+ }
+
+ default:
if (offset + 1 < dim)
{
// Set up the right hand side
double tmp
= values[q_point
+ GeometryInfo<dim>::lines_per_cell
- * n_edge_points + 4
- * n_face_points] (0);
+ * n_edge_points
+ + 5 * n_face_points] (0);
for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(i + 2)
- * this->degree
- + j]
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(i + 6)
+ * this->degree + j]
* this->shape_value_component
- ((i + 2) * this->degree
- + j,
+ ((i + 6) * this->degree + j,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
- + 4
+ + 5
* n_face_points],
0);
- for (unsigned int i = 0; i <= deg; ++i)
+ for (unsigned int i = 0; i < this->degree; ++i)
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= boundary_weights
system_matrix_inv.vmult
(solution, system_rhs);
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
for (unsigned int j = 0; j < deg; ++j)
if (std::abs (solution (i * deg + j))
> 1e-14)
local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 8 * this->degree) * deg
- + j
- + GeometryInfo<dim>::lines_per_cell]
+ + 10 * this->degree)
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell]
= solution (i * deg + j);
}
double tmp
= values[q_point
+ GeometryInfo<dim>::lines_per_cell
- * n_edge_points + 4
+ * n_edge_points + 5
* n_face_points] (1);
for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[i * this->degree + j]
+ for (unsigned int j = 0; j < this->degree; ++j)
+ tmp -= local_dofs[(i + 4)
+ * this->degree + j]
* this->shape_value_component
- (i * this->degree + j,
+ ((i + 4) * this->degree + j,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
- + 4
+ + 5
* n_face_points],
1);
- for (unsigned i = 0; i <= deg; ++i)
+ for (unsigned i = 0; i < this->degree; ++i)
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= boundary_weights
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i <= deg; ++i)
+ for (unsigned int i = 0; i < this->degree; ++i)
for (unsigned int j = 0; j < deg; ++j)
if (std::abs (solution (i * deg + j))
> 1e-14)
local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + 9 * deg)
- * this->degree]
+ + 11 * deg) * this->degree]
= solution (i * deg + j);
}
+ }
+ }
- break;
- }
-
- default:
- if (offset + 1 < dim)
- {
- // Set up the right hand side
- // for the horizontal shape
+ // Finally we project
+ // the remaining parts
+ // of the function on
+ // the interior shape
// functions.
- if (offset == 0)
- {
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 5 * n_face_points] (0);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(i + 6)
- * this->degree + j]
- * this->shape_value_component
- ((i + 6) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 5
- * n_face_points],
- 0);
+ const QGauss<dim> reference_quadrature (this->degree);
+ const unsigned int&
+ n_interior_points = reference_quadrature.size ();
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j)) * tmp;
- }
+ // We create the
+ // system matrix.
+ system_matrix.reinit (this->degree * deg * deg,
+ this->degree * deg * deg);
+ system_matrix = 0;
- system_matrix_inv.vmult
- (solution, system_rhs);
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < deg; ++k)
+ for (unsigned int l = 0; l < this->degree; ++l)
+ for (unsigned int m = 0; m < deg; ++m)
+ for (unsigned int n = 0; n < deg; ++n)
+ for (unsigned int q_point = 0;
+ q_point < n_interior_points; ++q_point)
+ system_matrix ((i * deg + j) * deg + k,
+ (l * deg + m) * deg + n)
+ += reference_quadrature.weight (q_point)
+ * legendre_polynomials[i].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (0)) * lobatto_polynomials[j + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (1))
+ * lobatto_polynomials[k + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (2))
+ * lobatto_polynomials_grad[l].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (0))
+ * lobatto_polynomials[m + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (1))
+ * lobatto_polynomials[n + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (2));
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 10 * this->degree)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
+ system_matrix_inv.reinit (system_matrix.m (),
+ system_matrix.m ());
+ system_matrix_inv.invert (system_matrix);
+ system_rhs.reinit (system_matrix_inv.m ());
+ solution.reinit (system_matrix.m ());
- // Set up the right hand side
- // for the vertical shape
- // functions.
- system_rhs = 0;
+ if (offset + 1 < dim)
+ {
+ if (offset == 0)
+ {
+ // Set up the right hand side.
+ system_rhs = 0;
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
+ for (unsigned int q_point = 0;
+ q_point < n_interior_points; ++q_point)
+ {
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points] (0);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
{
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points + 5
- * n_face_points] (1);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(i + 4)
- * this->degree + j]
+ for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int k = 0; k < 2; ++k)
+ tmp -= local_dofs[i + (j + 4 * k + 2)
+ * this->degree]
* this->shape_value_component
- ((i + 4) * this->degree + j,
+ (i + (j + 4 * k + 2)
+ * this->degree,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
- + 5
+ + GeometryInfo<dim>::faces_per_cell
* n_face_points],
- 1);
+ 0);
- for (unsigned i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j) + 1) * tmp;
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < 4; ++k)
+ tmp -= local_dofs[(i + 2 * (k + 2)
+ * this->degree
+ + GeometryInfo<dim>::lines_per_cell)
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell]
+ * this->shape_value_component
+ ((i + 2 * (k + 2) * this->degree
+ + GeometryInfo<dim>::lines_per_cell)
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points],
+ 0);
}
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
+ for (unsigned int i = 0; i < this->degree; ++i)
for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + 11 * deg) * this->degree]
- = solution (i * deg + j);
- }
- }
- }
-
- // Finally we project
- // the remaining parts
- // of the function on
- // the interior shape
- // functions.
- const QGauss<dim> reference_quadrature (this->degree);
- const unsigned int&
- n_interior_points = reference_quadrature.size ();
-
- // We create the
- // system matrix.
- system_matrix.reinit (this->degree * deg * deg,
- this->degree * deg * deg);
- system_matrix = 0;
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- for (unsigned int l = 0; l <= deg; ++l)
- for (unsigned int m = 0; m < deg; ++m)
- for (unsigned int n = 0; n < deg; ++n)
- for (unsigned int q_point = 0;
- q_point < n_interior_points; ++q_point)
- system_matrix ((i * deg + j) * deg + k,
- (l * deg + m) * deg + n)
- += reference_quadrature.weight (q_point)
- * legendre_polynomials[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (0)) * lobatto_polynomials[j + 2].value
+ for (unsigned int k = 0; k < deg; ++k)
+ system_rhs ((i * deg + j) * deg + k)
+ += reference_quadrature.weight (q_point)
+ * tmp
+ * lobatto_polynomials_grad[i].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (0))
+ * lobatto_polynomials[j + 2].value
(this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
+ GeometryInfo<dim>::faces_per_cell
* n_face_points]
- (1))
+ (1))
* lobatto_polynomials[k + 2].value
(this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
+ GeometryInfo<dim>::faces_per_cell
* n_face_points]
- (2))
- * lobatto_polynomials_grad[l].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (0))
- * lobatto_polynomials[m + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (1))
- * lobatto_polynomials[n + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (2));
+ (2));
+ }
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
- system_matrix_inv.invert (system_matrix);
- system_rhs.reinit (system_matrix_inv.m ());
- solution.reinit (system_matrix.m ());
+ system_matrix_inv.vmult (solution, system_rhs);
+
+ // Add the computed values
+ // to the resulting vector
+ // only, if they are not
+ // too small.
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < deg; ++k)
+ if (std::abs (solution ((i * deg + j) * deg + k))
+ > 1e-14)
+ local_dofs[((i + 2
+ * GeometryInfo<dim>::faces_per_cell)
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell
+ + 2
+ * GeometryInfo<dim>::faces_per_cell)
+ * deg + k
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution ((i * deg + j) * deg + k);
+ }
- if (offset + 1 < dim)
- {
- if (offset == 0)
- {
// Set up the right hand side.
- system_rhs = 0;
+ system_rhs = 0;
- for (unsigned int q_point = 0;
- q_point < n_interior_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points] (0);
-
- for (unsigned int i = 0; i <= deg; ++i)
+ for (unsigned int q_point = 0; q_point < n_interior_points;
+ ++q_point)
+ {
+ double tmp
+ = values[q_point + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points] (1);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < 2; ++j)
{
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- tmp -= local_dofs[i + (j + 4 * k + 2)
+ for (unsigned int k = 0; k < 2; ++k)
+ tmp -= local_dofs[i + (4 * j + k)
* this->degree]
- * this->shape_value_component
- (i + (j + 4 * k + 2)
- * this->degree,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 0);
-
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < 4; ++k)
- tmp -= local_dofs[(i + 2 * (k + 2)
- * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- * this->shape_value_component
- ((i + 2 * (k + 2) * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 0);
- }
+ * this->shape_value_component
+ (i + (4 * j + k) * this->degree,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points],
+ 1);
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
for (unsigned int k = 0; k < deg; ++k)
- system_rhs ((i * deg + j) * deg + k)
- += reference_quadrature.weight (q_point)
- * tmp
- * lobatto_polynomials_grad[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (0))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
+ tmp -= local_dofs[(i + 2 * j * this->degree
+ + GeometryInfo<dim>::lines_per_cell)
+ * deg + k
+ + GeometryInfo<dim>::lines_per_cell]
+ * this->shape_value_component
+ ((i + 2 * j * this->degree
+ + GeometryInfo<dim>::lines_per_cell)
+ * deg + k
+ + GeometryInfo<dim>::lines_per_cell,
+ this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
+ GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (1))
- * lobatto_polynomials[k + 2].value
- (this->generalized_support_points[q_point
+ * n_face_points],
+ 1)
+ + local_dofs[i + ((2 * j + 9) * deg + k
+ + GeometryInfo<dim>::lines_per_cell)
+ * this->degree]
+ * this->shape_value_component
+ (i + ((2 * j + 9) * deg + k
+ + GeometryInfo<dim>::lines_per_cell)
+ * this->degree,
+ this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
+ GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (2));
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- if (std::abs (solution ((i * deg + j) * deg + k))
- > 1e-14)
- local_dofs[((i + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell
- + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell]
- = solution ((i * deg + j) * deg + k);
- }
-
- // Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points] (1);
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- {
- for (unsigned int k = 0; k < 2; ++k)
- tmp -= local_dofs[i + (4 * j + k)
- * this->degree]
- * this->shape_value_component
- (i + (4 * j + k) * this->degree,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 1);
+ * n_face_points],
+ 1);
+ }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
for (unsigned int k = 0; k < deg; ++k)
- tmp -= local_dofs[(i + 2 * j * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell]
- * this->shape_value_component
- ((i + 2 * j * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 1)
- + local_dofs[i + ((2 * j + 9) * deg + k
- + GeometryInfo<dim>::lines_per_cell)
- * this->degree]
- * this->shape_value_component
- (i + ((2 * j + 9) * deg + k
- + GeometryInfo<dim>::lines_per_cell)
- * this->degree,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 1);
- }
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- system_rhs ((i * deg + j) * deg + k)
- += reference_quadrature.weight (q_point) * tmp
- * lobatto_polynomials_grad[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (1))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (0))
- * lobatto_polynomials[k + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (2));
- }
+ system_rhs ((i * deg + j) * deg + k)
+ += reference_quadrature.weight (q_point) * tmp
+ * lobatto_polynomials_grad[i].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (1))
+ * lobatto_polynomials[j + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (0))
+ * lobatto_polynomials[k + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (2));
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- if (std::abs (solution ((i * deg + j) * deg + k))
- > 1e-14)
- local_dofs[((i + this->degree + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell]
- = solution ((i * deg + j) * deg + k);
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < deg; ++k)
+ if (std::abs (solution ((i * deg + j) * deg + k))
+ > 1e-14)
+ local_dofs[((i + this->degree + 2
+ * GeometryInfo<dim>::faces_per_cell)
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell + 2
+ * GeometryInfo<dim>::faces_per_cell)
+ * deg + k
+ + GeometryInfo<dim>::lines_per_cell]
+ = solution ((i * deg + j) * deg + k);
+ }
// Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points] (2);
+ system_rhs = 0;
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < 4; ++j)
- {
- tmp -= local_dofs[i + (j + 8) * this->degree]
- * this->shape_value_component
- (i + (j + 8) * this->degree,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 2);
+ for (unsigned int q_point = 0; q_point < n_interior_points;
+ ++q_point)
+ {
+ double tmp
+ = values[q_point + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points] (2);
- for (unsigned int k = 0; k < deg; ++k)
- tmp -= local_dofs[i + ((2 * j + 1) * deg + k
- + GeometryInfo<dim>::lines_per_cell)
- * this->degree]
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < 4; ++j)
+ {
+ tmp -= local_dofs[i + (j + 8) * this->degree]
* this->shape_value_component
- (i + ((2 * j + 1) * deg + k
- + GeometryInfo<dim>::lines_per_cell)
- * this->degree,
+ (i + (j + 8) * this->degree,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
+ GeometryInfo<dim>::faces_per_cell
* n_face_points],
2);
- }
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- system_rhs ((i * deg + j) * deg + k)
- += reference_quadrature.weight (q_point) * tmp
- * lobatto_polynomials_grad[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (2))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (0))
- * lobatto_polynomials[k + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (1));
- }
+ for (unsigned int k = 0; k < deg; ++k)
+ tmp -= local_dofs[i + ((2 * j + 1) * deg + k
+ + GeometryInfo<dim>::lines_per_cell)
+ * this->degree]
+ * this->shape_value_component
+ (i + ((2 * j + 1) * deg + k
+ + GeometryInfo<dim>::lines_per_cell)
+ * this->degree,
+ this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points],
+ 2);
+ }
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < deg; ++k)
+ system_rhs ((i * deg + j) * deg + k)
+ += reference_quadrature.weight (q_point) * tmp
+ * lobatto_polynomials_grad[i].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (2))
+ * lobatto_polynomials[j + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (0))
+ * lobatto_polynomials[k + 2].value
+ (this->generalized_support_points[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points
+ + GeometryInfo<dim>::faces_per_cell
+ * n_face_points]
+ (1));
+ }
- system_matrix_inv.vmult (solution, system_rhs);
+ system_matrix_inv.vmult (solution, system_rhs);
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- if (std::abs (solution ((i * deg + j) * deg + k))
- > 1e-14)
- local_dofs[i + ((j + 2
- * (deg + GeometryInfo<dim>::faces_per_cell))
- * deg + k
- + GeometryInfo<dim>::lines_per_cell)
- * this->degree]
- = solution ((i * deg + j) * deg + k);
- }
+ for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int j = 0; j < deg; ++j)
+ for (unsigned int k = 0; k < deg; ++k)
+ if (std::abs (solution ((i * deg + j) * deg + k))
+ > 1e-14)
+ local_dofs[i + ((j + 2
+ * (deg + GeometryInfo<dim>::faces_per_cell))
+ * deg + k
+ + GeometryInfo<dim>::lines_per_cell)
+ * this->degree]
+ = solution ((i * deg + j) * deg + k);
+ }
- break;
- }
+ break;
+ }
- default:
- Assert (false, ExcNotImplemented ());
- }
+ default:
+ Assert (false, ExcNotImplemented ());
+ }
+ }
}
const VectorSlice<const std::vector<std::vector<double> > >& values)
const
{
- const unsigned int deg = this->degree-1;
Assert (values.size () == this->n_components (),
ExcDimensionMismatch (values.size (), this->n_components ()));
Assert (values[0].size () == this->generalized_support_points.size (),
ExcDimensionMismatch (local_dofs.size (), this->dofs_per_cell));
std::fill (local_dofs.begin (), local_dofs.end (), 0.0);
+ const unsigned int deg = this->degree-1;
+ const std::vector<Polynomials::Polynomial<double> >&
+ legendre_polynomials
+ = Polynomials::Legendre::generate_complete_basis (deg);
+ const QGauss<1> edge_quadrature (this->degree);
+ const std::vector<Point<1> >&
+ edge_quadrature_points = edge_quadrature.get_points ();
+ const unsigned int& n_edge_points = edge_quadrature.size ();
+
switch (dim)
{
case 2:
{
// Let us begin with the
- // interpolation part.
- const QGauss<dim - 1> reference_edge_quadrature (this->degree);
- const unsigned int&
- n_edge_points = reference_edge_quadrature.size ();
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[(i + 2 * j) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[1 - j][q_point + (i + 2 * j) * n_edge_points];
-
+ // edge degrees of freedom.
+ for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+ {
+ const double weight = edge_quadrature.weight (q_point);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ {
+ const double L_i
+ = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+
+ for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int k = 0; k < 2; ++k)
+ local_dofs[i + (j + 2 * k) * this->degree]
+ += L_i * values[1 - k][q_point + (j + 2 * k) * n_edge_points];
+ }
+ }
+
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14)
- local_dofs[(i + 2 * j) * this->degree] = 0.0;
- }
+ for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+ for (unsigned int i = 0; i < this->degree; ++i)
+ if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+ local_dofs[i + e * this->degree] = 0.0;
// If the degree is greater
// than 0, then we have still
- // some higher order edge
- // shape functions to
- // consider.
- // Here the projection part
- // starts. The dof values
- // are obtained by solving
- // a linear system of
- // equations.
- if (this->degree-1 > 1)
+ // some interpolations onto
+ // the interior shape
+ // functions left.
+ if (this->degree > 1)
{
- // We start with projection
- // on the higher order edge
- // shape function.
const std::vector<Polynomials::Polynomial<double> >&
lobatto_polynomials
= Polynomials::Lobatto::generate_complete_basis
(this->degree);
- FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
std::vector<Polynomials::Polynomial<double> >
lobatto_polynomials_grad (this->degree);
lobatto_polynomials_grad[i]
= lobatto_polynomials[i + 1].derivative ();
- // Set up the system matrix.
- // This can be used for all
- // edges.
- for (unsigned int i = 0; i < system_matrix.m (); ++i)
- for (unsigned int j = 0; j < system_matrix.n (); ++j)
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- system_matrix (i, j)
- += boundary_weights (q_point, j)
- * lobatto_polynomials_grad[i + 1].value
- (this->generalized_face_support_points[q_point]
- (1));
-
- FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
-
- system_matrix_inv.invert (system_matrix);
-
- const unsigned int
- line_coordinate[GeometryInfo<2>::lines_per_cell]
- = {1, 1, 0, 0};
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell; ++line)
- {
- // Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- {
- const double tmp
- = values[line_coordinate[line]][line * n_edge_points
- + q_point]
- - local_dofs[line * this->degree]
- * this->shape_value_component
- (line * this->degree,
- this->generalized_support_points[line
- * n_edge_points
- + q_point],
- line_coordinate[line]);
-
- for (unsigned int i = 0; i < system_rhs.size (); ++i)
- system_rhs (i) += boundary_weights (q_point, i) * tmp;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < solution.size (); ++i)
- if (std::abs (solution (i)) > 1e-14)
- local_dofs[line * this->degree + i + 1] = solution (i);
- }
-
- // Then we go on to the
+ // We set up the system
+ // matrix and use it for
+ // both, the horizontal
+ // and the vertical
// interior shape
- // functions. Again we
- // set up the system
- // matrix and use it
- // for both, the
- // horizontal and the
- // vertical, interior
- // shape functions.
+ // functions.
const QGauss<dim> reference_quadrature (this->degree);
- const unsigned int&
- n_interior_points = reference_quadrature.size ();
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (this->degree-1);
+ const unsigned int& n_interior_points
+ = reference_quadrature.size ();
- system_matrix.reinit ((this->degree-1) * this->degree,
- (this->degree-1) * this->degree);
- system_matrix = 0;
+ FullMatrix<double> system_matrix (deg * this->degree,
+ deg * this->degree);
for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
+ for (unsigned int j = 0; j < deg; ++j)
for (unsigned int k = 0; k < this->degree; ++k)
- for (unsigned int l = 0; l < this->degree-1; ++l)
+ for (unsigned int l = 0; l < deg; ++l)
for (unsigned int q_point = 0;
q_point < n_interior_points; ++q_point)
- system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
+ system_matrix (i * deg + j, k * deg + l)
+= reference_quadrature.weight (q_point)
* legendre_polynomials[i].value
(this->generalized_support_points[q_point
* n_edge_points]
(1));
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
+ FullMatrix<double> system_matrix_inv (system_matrix.m (),
+ system_matrix.m ());
+
system_matrix_inv.invert (system_matrix);
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- system_rhs.reinit (system_matrix_inv.m ());
- system_rhs = 0;
+
+ Vector<double> solution (system_matrix_inv.m ());
+ Vector<double> system_rhs (system_matrix.m ());
for (unsigned int q_point = 0; q_point < n_interior_points;
++q_point)
case 3:
{
// Let us begin with the
- // interpolation part.
- const QGauss<1> reference_edge_quadrature (this->degree);
- const unsigned int&
- n_edge_points = reference_edge_quadrature.size ();
-
+ // edge degrees of freedom.
for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
{
- for (unsigned int i = 0; i < 4; ++i)
- local_dofs[(i + 8) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[2][q_point + (i + 8) * n_edge_points];
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- local_dofs[(i + 2 * (2 * j + k)) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[1 - k][q_point + (i + 2 * (2 * j + k))
- * n_edge_points];
+ const double weight = edge_quadrature.weight (q_point);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
+ {
+ const double L_i
+ = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
+
+ for (unsigned int j = 0; j < 4; ++j)
+ local_dofs[i + (j + 8) * this->degree]
+ += L_i * values[2][q_point + (j + 8) * n_edge_points];
+
+ for (unsigned int j = 0; j < 2; ++j)
+ for (unsigned int k = 0; k < 2; ++k)
+ for (unsigned int l = 0; l < 2; ++l)
+ local_dofs[i + (j + 4 * k + 2 * l) * this->degree]
+ += L_i * values[1 - l][q_point + (j + 4 * k + 2 * l) * n_edge_points];
+ }
}
-
+
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i < 4; ++i)
- if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
- local_dofs[(i + 8) * this->degree] = 0.0;
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- if (std::abs (local_dofs[(i + 2 * (2 * j + k)) * this->degree])
- < 1e-14)
- local_dofs[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
+ for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+ for (unsigned int i = 0; i < this->degree; ++i)
+ if (std::abs (local_dofs[i + e * this->degree]) < 1e-14)
+ local_dofs[i + e * this->degree] = 0.0;
// If the degree is greater
// than 0, then we have still
- // some higher order shape
- // functions to consider.
- // Here the projection part
- // starts. The dof values
- // are obtained by solving
- // a linear system of
- // equations.
+ // some interpolation to the
+ // face and interior shape
+ // functions left.
if (this->degree > 1)
{
- // We start with projection
- // on the higher order edge
- // shape function.
const std::vector<Polynomials::Polynomial<double> >&
lobatto_polynomials
= Polynomials::Lobatto::generate_complete_basis
(this->degree);
- FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
+ FullMatrix<double> system_matrix (deg * this->degree,
+ deg * this->degree);
std::vector<Polynomials::Polynomial<double> >
lobatto_polynomials_grad (this->degree);
lobatto_polynomials_grad[i]
= lobatto_polynomials[i + 1].derivative ();
- // Set up the system matrix.
- // This can be used for all
- // edges.
- for (unsigned int i = 0; i < system_matrix.m (); ++i)
- for (unsigned int j = 0; j < system_matrix.n (); ++j)
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- system_matrix (i, j)
- += boundary_weights (q_point, j)
- * lobatto_polynomials_grad[i + 1].value
- (this->generalized_face_support_points[q_point]
- (1));
-
- FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
-
- system_matrix_inv.invert (system_matrix);
-
- const unsigned int
- line_coordinate[GeometryInfo<3>::lines_per_cell]
- = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell; ++line)
- {
- // Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < this->degree; ++q_point)
- {
- const double tmp
- = values[line_coordinate[line]][line * this->degree
- + q_point]
- - local_dofs[line * this->degree]
- * this->shape_value_component
- (line * this->degree,
- this->generalized_support_points[line
- * this->degree
- + q_point],
- line_coordinate[line]);
-
- for (unsigned int i = 0; i < system_rhs.size (); ++i)
- system_rhs (i) += boundary_weights (q_point, i)
- * tmp;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < solution.size (); ++i)
- if (std::abs (solution (i)) > 1e-14)
- local_dofs[line * this->degree + i + 1] = solution (i);
- }
-
- // Then we go on to the
- // face shape functions.
- // Again we set up the
- // system matrix and
- // use it for both, the
- // horizontal and the
- // vertical, shape
+ // We set up the system
+ // matrix and use it for
+ // both, the horizontal
+ // and the vertical, shape
// functions.
- const std::vector<Polynomials::Polynomial<double> >&
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (this->degree-1);
- const unsigned int n_face_points = n_edge_points * n_edge_points;
-
- system_matrix.reinit ((this->degree-1) * this->degree,
- (this->degree-1) * this->degree);
- system_matrix = 0;
+ const unsigned int
+ n_face_points = n_edge_points * n_edge_points;
for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
+ for (unsigned int j = 0; j < deg; ++j)
for (unsigned int k = 0; k < this->degree; ++k)
- for (unsigned int l = 0; l < this->degree-1; ++l)
+ for (unsigned int l = 0; l < deg; ++l)
for (unsigned int q_point = 0; q_point < n_face_points;
++q_point)
- system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
+ system_matrix (i * deg + j, k * deg + l)
+= boundary_weights (q_point + n_edge_points,
- 2 * (k * (this->degree-1) + l))
+ 2 * (k * deg + l))
* legendre_polynomials[i].value
(this->generalized_face_support_points[q_point
+ 4
* n_edge_points]
(1));
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
+ FullMatrix<double> system_matrix_inv (system_matrix.m (),
+ system_matrix.m ());
+
system_matrix_inv.invert (system_matrix);
- solution.reinit (system_matrix.m ());
- system_rhs.reinit (system_matrix.m ());
-
+
const unsigned int
face_coordinates[GeometryInfo<3>::faces_per_cell][2]
= {{1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
edge_indices[GeometryInfo<3>::faces_per_cell][GeometryInfo<3>::lines_per_face]
= {{0, 4, 8, 10}, {1, 5, 9, 11}, {8, 9, 2, 6},
{10, 11, 3, 7}, {2, 3, 0, 1}, {6, 7, 4, 5}};
+ Vector<double> solution (system_matrix.m ());
+ Vector<double> system_rhs (system_matrix.m ());
for (unsigned int face = 0;
face < GeometryInfo<dim>::faces_per_cell; ++face)
projects / dealii-svn.git / commitdiff