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 face shape functions
- // of the parent cell.
+ // to the lowest order 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 <= 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);
- }
+ {
+ 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);
+ }
else
- 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);
- }
+ {
+ 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);
+ }
}
- // Then interpolate the shape functions
- // of the child cells to the interior
- // shape functions of the parent cell.
+ // Then project the shape functions
+ // of the child cells to the higher
+ // order 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 ();
- for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
- for (unsigned int q_point = 0; q_point < 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;
++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;
+ const double weight
+ = std::sqrt (quadrature.weight (q_point));
- if (quadrature_points[q_point] (0) < 0.5)
+ for (unsigned int i = 0; i < this->degree; ++i)
{
+ 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)
- k = 0;
-
+ {
+ 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);
+ }
+
else
{
- quadrature_point (1) -= 1.0;
- k = 1;
+ 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);
}
- }
-
- else
- if (quadrature_points[q_point] (1) < 0.5)
+
+ 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);
+ }
+
+ tmp *= quadrature.weight (q_point);
+
+ for (unsigned int i = 0; i < this->degree; ++i)
{
- quadrature_point (0) -= 1.0;
- k = 2;
+ 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;
+ }
+ }
}
-
- else
+ }
+
+ 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)
{
- quadrature_point (0) -= 1.0;
- quadrature_point (1) -= 1.0;
- k = 3;
+ 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);
}
-
- 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 edge shape functions
- // of the parent cell.
+ // to the lowest order 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 < 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 (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));
+ 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 j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
+ else
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < 2; ++j)
{
- Point<dim> quadrature_point (j,
+ Point<dim> quadrature_point (i,
2.0 * edge_quadrature_points[q_point] (0)
- - 1.0, k);
+ - 1.0, j);
- this->restriction[index][j + 4 * k + 2]
- (i + (j + 4 * k) * this->degree, dof)
- += L_i * this->shape_value_component (dof,
- quadrature_point,
- 1);
+ this->restriction[index][i + 4 * j + 2]
+ ((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) - 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,
+ = 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,
2.0 * edge_quadrature_points[q_point] (0)
- 1.0);
- 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);
+ this->restriction[index][i + 2 * (j + 2)]
+ ((i + 2 * (j + 4)) * this->degree, dof)
+ += weight * this->shape_value_component (dof,
+ quadrature_point,
+ 2);
}
- }
}
- // Then interpolate the shape functions
- // of the child cells to the face
- // and interior shape functions of
- // the parent cell.
+ // Then project the shape functions
+ // of the child cells to the higher
+ // order shape functions of the
+ // parent cell.
if (this->degree > 1)
{
- const QGauss<2> face_quadrature (2 * this->degree);
- const std::vector<Point<2> >& face_quadrature_points
- = face_quadrature.get_points ();
+ const unsigned int deg = this->degree-1;
const std::vector<Polynomials::Polynomial<double> >&
- 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 ();
- 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)
+ 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)
{
- const double weight = 2.0 * face_quadrature.weight (q_point);
-
- 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 (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)
+ 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)
{
- m = 0;
-
- if (face_quadrature_points[q_point] (1) < 0.5)
- l = 0;
-
- else
- {
- quadrature_point (2) -= 1.0;
- l = 1;
- }
+ 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
{
- 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;
- }
+ 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 j = 0; j < this->degree; ++j)
+
+ for (unsigned int k = 0; k < deg; ++k)
{
- 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);
- }
+ 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);
+ }
}
-
- // Then, the interpolation to
- // the interior shape functions.
- for (unsigned int q_point = 0; q_point < n_quadrature_points;
- ++q_point)
- {
- 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
- {
- 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
- {
- quadrature_point (2) -= 1.0;
- l = 1;
- }
- }
- }
-
- else
- {
- quadrature_point (0) -= 1.0;
- n = 1;
-
- if (quadrature_points[q_point] (1) < 0.5)
+
+ const QGauss<2> face_quadrature (2 * this->degree);
+ const std::vector<Point<2> >& face_quadrature_points
+ = face_quadrature.get_points ();
+ const std::vector<Polynomials::Polynomial<double> >&
+ lobatto_polynomials
+ = Polynomials::Lobatto::generate_complete_basis
+ (this->degree);
+ const unsigned int n_edge_dofs
+ = GeometryInfo<dim>::lines_per_cell * this->degree;
+ 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)
{
- m = 0;
-
- if (quadrature_points[q_point] (2) < 0.5)
- l = 0;
-
- else
- {
- quadrature_point (2) -= 1.0;
- l = 1;
- }
+ 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);
}
-
+
else
{
- quadrature_point (1) -= 1.0;
- m = 1;
-
+ 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);
+ }
+
+ 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)
+ {
+ 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)
+ {
+ tmp = 0.0;
+
+ if (quadrature_points[q_point] (0) < 0.5)
+ {
+ if (quadrature_points[q_point] (1) < 0.5)
+ {
if (quadrature_points[q_point] (2) < 0.5)
- l = 0;
-
+ {
+ 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);
+ }
+
else
{
- quadrature_point (2) -= 1.0;
- l = 1;
+ 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);
}
}
+
+ 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);
+ }
}
-
- const double weight = 2.0 * quadrature.weight (q_point);
-
- for (unsigned int i = 0; i < this->degree; ++i)
+
+ 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)
+ {
+ 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);
+ }
+
+ 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);
+ }
+
+ 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 L_i_0
- = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0));
+ = legendre_polynomials[i].value
+ (quadrature_points[q_point] (0));
const double L_i_1
- = weight * legendre_polynomials[i].value (quadrature_points[q_point] (1));
+ = legendre_polynomials[i].value
+ (quadrature_points[q_point] (1));
const double L_i_2
- = weight * legendre_polynomials[i].value (quadrature_points[q_point] (2));
-
+ = legendre_polynomials[i].value
+ (quadrature_points[q_point] (2));
+
for (unsigned int j = 0; j < deg; ++j)
{
- 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));
-
+ 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));
+
for (unsigned int k = 0; k < deg; ++k)
{
- const double L_k
- = legendre_polynomials[k + 1].value (quadrature_points[q_point] (2));
+ const double l_k_0
+ = l_j_0 * lobatto_polynomials[k + 2].value
+ (quadrature_points[q_point] (2));
const double l_k_1
- = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (1));
+ = l_j_1 * lobatto_polynomials[k + 2].value
+ (quadrature_points[q_point] (2));
const double 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);
+ = 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;
+ }
}
}
}
}
+
+ 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);
+ }
}
}
}
+
template <int dim>
std::vector<unsigned int>
FE_Nedelec<dim>::get_dpo_vector (const unsigned int degree, bool dg)
return dpo;
}
-
//---------------------------------------------------------------------------
// Data field initialization
//---------------------------------------------------------------------------
+ GeometryInfo<dim>::lines_per_face)
= 1.0;
interpolation_matrix
- (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)
+ (i + (j + GeometryInfo<dim>::lines_per_face) * source_fe.degree,
+ i + (j + GeometryInfo<dim>::lines_per_face) * 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;
- // Set the degrees of freedom
- // by interpolation as usual.
+ // Perform projection-based 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 q_point = 0; q_point < n_edge_quadrature_points;
- ++q_point)
+ 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)
{
- 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);
-
+ 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);
+
for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
{
- 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));
+ 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 shifts[4][2] = { { 0.0, 0.0 }, { 1.0, 0.0 },
{ 0.0, 1.0 }, { 1.0, 1.0 } };
- 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 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);
-
- 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, 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));
- }
- }
- }
+ 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)
{
- const QGauss<2> quadrature (source_fe.degree);
- const std::vector<Point<2> >&
- quadrature_points = quadrature.get_points ();
const std::vector<Polynomials::Polynomial<double> >&
- lobatto_polynomials
- = Polynomials::Lobatto::generate_complete_basis
- (source_fe.degree);
- const unsigned int n_boundary_dofs
- = GeometryInfo<dim>::lines_per_face * source_fe.degree;
- const unsigned int& n_quadrature_points = quadrature.size ();
- FullMatrix<double>
- system_matrix_inv (source_fe.degree * (source_fe.degree - 1),
- source_fe.degree * (source_fe.degree - 1));
+ 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 * (source_fe.degree - 1),
- n_quadrature_points);
+ FullMatrix<double> assembling_matrix (source_fe.degree - 1,
+ n_edge_quadrature_points);
- for (unsigned int q_point = 0; q_point < n_quadrature_points;
- ++q_point)
+ for (unsigned int q_point = 0;
+ q_point < n_edge_quadrature_points; ++q_point)
{
- const double weight = std::sqrt (quadrature.weight (q_point));
+ const double weight
+ = std::sqrt (edge_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 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);
+ }
+
+ 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)
+ {
+ 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);
+
+ 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;
+ }
+ }
+
+ 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);
+ }
+
+ const QGauss<2> quadrature (source_fe.degree);
+ const std::vector<Point<2> >&
+ quadrature_points = quadrature.get_points ();
+ const std::vector<Polynomials::Polynomial<double> >&
+ lobatto_polynomials
+ = Polynomials::Lobatto::generate_complete_basis
+ (source_fe.degree);
+ const unsigned int n_boundary_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,
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);
}
- FullMatrix<double> solution (system_matrix_inv.m (), 2);
- FullMatrix<double> system_rhs (system_matrix_inv.m (), 2);
- Vector<double> tmp (2);
+ 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)
{
// its values at the generalized support
// points in the finite element space on the
// reference cell.
- // 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).
+ // This is done as usual by projection-based
+ // interpolation.
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)
- {
- 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:
- {
+ switch (dim)
+ {
+ case 2:
+ {
+ const QGauss<1> reference_edge_quadrature (this->degree);
+ const unsigned int& n_edge_points
+ = reference_edge_quadrature.size ();
+
// Let us begin with the
- // edge degrees of freedom.
- for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+ // 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);
+
+ // 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;
+ }
+
+ if (offset == 0)
+ for (unsigned int i = 0; i < 2; ++i)
{
- const double weight = edge_quadrature.weight (q_point);
-
- for (unsigned int i = 0; i < this->degree; ++i)
+ 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);
+
+ if (std::abs (local_dofs[(i + 2) * this->degree]) < 1e-14)
+ local_dofs[(i + 2) * 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)
+ {
+ // 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);
+
+ system_matrix_inv.invert (system_matrix);
+
+ 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)
+ if ((line < 2) || (offset == 0))
{
- const double L_i
- = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0));
-
- for (unsigned int j = 0; j < 2; ++j)
+ // Set up the right hand side.
+ system_rhs = 0;
+
+ for (unsigned int q_point = 0; q_point < n_edge_points;
+ ++q_point)
{
- 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);
+ 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;
}
- }
- }
-
+
+ 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 < GeometryInfo<dim>::lines_per_cell * this->degree; ++i)
- if (std::abs (local_dofs[i]) < 1e-14)
- local_dofs[i] = 0.0;
+ 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);
+ }
- // If the degree is greater
- // than 0, then we have still
- // some interpolations onto
- // the interior shape
- // functions left.
- if (this->degree > 1)
- {
- 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
+ // 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;
+
+ 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
+ (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.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.
- const QGauss<dim> reference_quadrature (this->degree);
- const unsigned int& n_interior_points
- = reference_quadrature.size ();
+ system_rhs = 0;
- FullMatrix<double> system_matrix (deg * this->degree,
- deg * this->degree);
+ 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 < 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
+ 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))
- * 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));
+ }
- FullMatrix<double> system_matrix_inv (system_matrix.m (),
- system_matrix.m ());
-
- system_matrix_inv.invert (system_matrix);
-
- Vector<double> solution (system_matrix_inv.m ());
- Vector<double> system_rhs (system_matrix.m ());
-
- if (offset == 0)
- {
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- 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));
- }
-
- 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 < 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 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);
+ system_rhs = 0;
- 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 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 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[j + 2].value
- (this->generalized_support_points[q_point
+ 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]
- (0));
- }
+ * n_edge_points],
+ 1);
- system_matrix_inv.vmult (solution, system_rhs);
+ 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));
+ }
+
+ 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
- + 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)
+ 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;
- }
+ break;
+ }
+
+ case 3:
+ {
+ const QGauss<1>
+ reference_edge_quadrature (this->degree);
+ const unsigned int&
+ n_edge_points = reference_edge_quadrature.size ();
- 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)
+ // 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);
+
+ // 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;
+ }
+
+ if (offset + 1 < dim)
+ {
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j < 2; ++j)
{
- 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);
- }
- }
- }
-
+ 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);
+
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell * this->degree; ++i)
- if (std::abs (local_dofs[i]) < 1e-14)
- local_dofs[i] = 0.0;
+ if (std::abs (local_dofs[(i + 4 * j) * this->degree])
+ < 1e-14)
+ local_dofs[(i + 4 * j) * this->degree] = 0.0;
+ }
- // 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;
+ 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 < 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 + 4 * j + 2)
+ * this->degree]) < 1e-14)
+ local_dofs[(i + 4 * j + 2) * 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 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));
- 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;
+ system_matrix_inv.invert (system_matrix);
- 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 line = 0;
+ line < GeometryInfo<dim>::lines_per_cell; ++line)
+ {
+ // Set up the right hand side.
+ 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 * 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);
+ 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 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 < solution.size (); ++i)
+ if (std::abs (solution (i)) > 1e-14)
+ local_dofs[line * this->degree + i + 1]
+ = solution (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)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
+ // 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 vertical shape
+ // for the horizontal 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);
+ 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[2 * (i + 4)
- * this->degree + j]
+ 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
- (2 * (i + 4) * this->degree + j,
+ (4 * i * this->degree + j,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points],
- 2);
+ 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;
- }
+ 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;
+ }
- 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 < 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);
-
- break;
- }
-
- case 1:
- {
- 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
- + n_face_points] (1);
+ 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 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);
+ // 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 < 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] (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;
+ }
- 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 < 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);
- }
+ // 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;
+ }
+ case 1:
+ {
+ if (offset + 1 < dim)
+ {
// Set up the right hand side
- // for the vertical shape
+ // for the horizontal 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]
- (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
+ + 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 < 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 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;
- }
+ 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;
+ }
- 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 < 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;
- }
+ 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);
+ }
- case 2:
- {
- if (offset == 0)
- {
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- system_rhs = 0;
+ // 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 + 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 + 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 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;
- }
+ 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);
+ 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
- + 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 <= 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;
+ }
+ case 2:
+ {
+ if (offset == 0)
+ {
// Set up the right hand side
- // for the vertical shape
+ // for the horizontal 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 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 < 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;
- }
+ 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;
+ }
- 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 < 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;
- }
-
- 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 <= 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);
+ }
- 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);
+ // 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[(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;
- }
+ 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;
+ }
- 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 < 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);
- }
+ // 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;
+ }
+ case 3:
+ {
+ if (offset == 0)
+ {
// Set up the right hand side
- // for the vertical shape
+ // for the horizontal 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 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);
- 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;
- }
+ 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;
+ }
- 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 < 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;
- }
-
- 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 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);
+ }
- 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;
- }
+ // Set up the right hand side
+ // for the vertical shape
+ // functions.
+ system_rhs = 0;
- system_matrix_inv.vmult (solution, system_rhs);
+ 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;
+ }
- // 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);
- }
+ system_matrix_inv.vmult (solution, system_rhs);
- break;
- }
+ // 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;
+ }
- default:
+ case 4:
+ {
if (offset + 1 < dim)
{
// Set up the right hand side
double tmp
= values[q_point
+ GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 5 * n_face_points] (0);
+ * 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 + 6)
- * this->degree + j]
+ for (unsigned int j = 0; j <= deg; ++j)
+ tmp -= local_dofs[(i + 2)
+ * this->degree
+ + j]
* this->shape_value_component
- ((i + 6) * this->degree + j,
+ ((i + 2) * this->degree
+ + j,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
- + 5
+ + 4
* n_face_points],
0);
- for (unsigned int i = 0; i < this->degree; ++i)
+ for (unsigned int i = 0; i <= deg; ++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 < this->degree; ++i)
+ // 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]
+ + 8 * 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 + 5
+ * 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 + 4)
- * this->degree + j]
+ for (unsigned int j = 0; j <= deg; ++j)
+ tmp -= local_dofs[i * this->degree + j]
* this->shape_value_component
- ((i + 4) * this->degree + j,
+ (i * this->degree + j,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
- + 5
+ + 4
* n_face_points],
1);
- for (unsigned i = 0; i < this->degree; ++i)
+ for (unsigned i = 0; i <= deg; ++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 < this->degree; ++i)
+ 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
- + 11 * deg) * this->degree]
+ + 9 * deg)
+ * this->degree]
= solution (i * deg + j);
}
- }
- }
- // Finally we project
- // the remaining parts
- // of the function on
- // the interior shape
+ break;
+ }
+
+ default:
+ if (offset + 1 < dim)
+ {
+ // Set up the right hand side
+ // for the horizontal shape
// functions.
- const QGauss<dim> reference_quadrature (this->degree);
- const unsigned int&
- n_interior_points = reference_quadrature.size ();
+ if (offset == 0)
+ {
+ system_rhs = 0;
- // We create the
- // system matrix.
- system_matrix.reinit (this->degree * deg * deg,
- this->degree * deg * deg);
- system_matrix = 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 < 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));
+ 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);
+
+ 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;
+ }
- 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);
- if (offset + 1 < dim)
- {
- if (offset == 0)
- {
- // Set up the right hand side.
- system_rhs = 0;
+ // 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);
+ }
- 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)
+ // 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)
{
- 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 + (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);
+ double tmp
+ = values[q_point
+ + GeometryInfo<dim>::lines_per_cell
+ * n_edge_points + 5
+ * n_face_points] (1);
- 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]
+ for (unsigned int i = 0; i < 2; ++i)
+ for (unsigned int j = 0; j <= deg; ++j)
+ tmp -= local_dofs[(i + 4)
+ * this->degree + j]
* this->shape_value_component
- ((i + 2 * (k + 2) * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell,
+ ((i + 4) * this->degree + j,
this->generalized_support_points[q_point
+ GeometryInfo<dim>::lines_per_cell
* n_edge_points
- + GeometryInfo<dim>::faces_per_cell
+ + 5
* n_face_points],
- 0);
+ 1);
+
+ 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 i = 0; i < this->degree; ++i)
+ 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)
- 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
+ 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
(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));
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
+ (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 < 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);
- }
+ 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 ());
+ 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_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);
+ system_rhs = 0;
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < 2; ++j)
+ 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 k = 0; k < 2; ++k)
- tmp -= local_dofs[i + (4 * j + k)
+ 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 * 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);
+ * 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 k = 0; k < deg; ++k)
- tmp -= local_dofs[(i + 2 * j * this->degree
+ 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 + k
- + GeometryInfo<dim>::lines_per_cell]
- * this->shape_value_component
- ((i + 2 * j * this->degree
+ * deg + j
+ + GeometryInfo<dim>::lines_per_cell]
+ * this->shape_value_component
+ ((i + 2 * (k + 2) * this->degree
+ GeometryInfo<dim>::lines_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell,
- this->generalized_support_points[q_point
+ * 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);
+ }
+
+ 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
+ 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
+ * 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],
- 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)
- 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));
- }
+ * 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 < 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
+ 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]
+ * 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;
+ 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);
+ 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 < 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
+ 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);
+
+ 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],
- 2);
+ * 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));
+ }
- 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);
- }
+ 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)
- 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
+ // 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);
+ }
+
+ // 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);
+
+ 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]
- (1));
- }
+ * n_face_points],
+ 2);
- system_matrix_inv.vmult (solution, system_rhs);
+ 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 <= 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));
+ }
+
+ 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 + ((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 <= 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);
+ }
- 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
- // 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];
- }
- }
-
+ // 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];
+
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell * this->degree; ++i)
- if (std::abs (local_dofs[i]) < 1e-14)
- local_dofs[i] = 0.0;
+ if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14)
+ local_dofs[(i + 2 * j) * this->degree] = 0.0;
+ }
// If the degree is greater
// than 0, then we have still
- // some interpolations onto
- // the interior shape
- // functions left.
- if (this->degree > 1)
+ // 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)
{
+ // 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 ();
- // We set up the system
- // matrix and use it for
- // both, the horizontal
- // and the vertical
+ // 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
// interior shape
- // functions.
+ // 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 unsigned int& n_interior_points
- = reference_quadrature.size ();
+ 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);
- FullMatrix<double> system_matrix (deg * this->degree,
- deg * this->degree);
+ 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 < deg; ++j)
+ 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 < deg; ++l)
+ 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 * deg + j, k * deg + l)
+ 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
* n_edge_points]
(1));
- FullMatrix<double> system_matrix_inv (system_matrix.m (),
- system_matrix.m ());
-
+ system_matrix_inv.reinit (system_matrix.m (),
+ system_matrix.m ());
system_matrix_inv.invert (system_matrix);
-
- Vector<double> solution (system_matrix_inv.m ());
- Vector<double> system_rhs (system_matrix.m ());
+ // Set up the right hand side
+ // for the horizontal shape
+ // functions.
+ system_rhs.reinit (system_matrix_inv.m ());
+ system_rhs = 0;
for (unsigned int q_point = 0; q_point < n_interior_points;
++q_point)
case 3:
{
// Let us begin with the
- // edge degrees of freedom.
+ // interpolation part.
+ const QGauss<1> reference_edge_quadrature (this->degree);
+ const unsigned int&
+ n_edge_points = reference_edge_quadrature.size ();
+
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[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];
- }
+ 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];
}
-
+
// Add the computed values
// to the resulting vector
// only, if they are not
// too small.
- for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell * this->degree; ++i)
- if (std::abs (local_dofs[i]) < 1e-14)
- local_dofs[i] = 0.0;
+ 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;
// If the degree is greater
// than 0, then we have still
- // some interpolation to the
- // face and interior shape
- // functions left.
+ // 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);
- FullMatrix<double> system_matrix (deg * this->degree,
- deg * 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 ();
- // We set up the system
- // matrix and use it for
- // both, the horizontal
- // and the vertical, shape
- // functions.
+ // 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
- n_face_points = n_edge_points * n_edge_points;
+ 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
+ // 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 < deg; ++j)
+ 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 < deg; ++l)
+ 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 * deg + j, k * deg + l)
+ system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
+= boundary_weights (q_point + n_edge_points,
- 2 * (k * deg + l))
+ 2 * (k * (this->degree-1) + l))
* legendre_polynomials[i].value
(this->generalized_face_support_points[q_point
+ 4
* n_edge_points]
(1));
- FullMatrix<double> system_matrix_inv (system_matrix.m (),
- system_matrix.m ());
-
+ system_matrix_inv.reinit (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