-
-
-// Since this is a vector valued element,
-// we cannot interpolate a scalar function.
-template <int dim>
-void FE_Nedelec<dim>::interpolate (std::vector<double> &, const std::vector<double> &) const
-{
- Assert(false, ExcNotImplemented ());
-}
-
-
-// Interpolate a function, which is given by
-// its values at the generalized support
-// points in the finite element space on the
-// reference cell.
-// This is done as usual by projection-based
-// interpolation.
-template <int dim>
-void
-FE_Nedelec<dim>::interpolate (std::vector<double> &local_dofs,
- const std::vector<Vector<double> > &values,
- const 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 ()));
- Assert (local_dofs.size () == this->dofs_per_cell,
- ExcDimensionMismatch (local_dofs.size (),this->dofs_per_cell));
- Assert (values[0].size () >= offset + this->n_components (),
- ExcDimensionMismatch (values[0].size (),
- offset + this->n_components ()));
- std::fill (local_dofs.begin (), local_dofs.end (), 0.);
-
- if (offset < dim)
- switch (dim)
- {
- case 2:
- {
- const QGauss<1> reference_edge_quadrature (this->degree);
- const unsigned int &n_edge_points
- = reference_edge_quadrature.size ();
-
- // Let us begin with the
- // interpolation part.
- for (unsigned int i = 0; i < 2; ++i)
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[i * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[q_point + i * n_edge_points] (1);
-
- // 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)
- {
- 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))
- {
- // 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 * 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 < 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. 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.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_interior_points; ++q_point)
- {
- double tmp
- = values[q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points] (0);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- tmp -= local_dofs[(i + 2) * this->degree + j]
- * this->shape_value_component
- ((i + 2) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- 0);
-
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- system_rhs (i * (this->degree-1) + j)
- += reference_quadrature.weight (q_point) * tmp
- * lobatto_polynomials_grad[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (0))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (1));
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell)
- * (this->degree-1) + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * (this->degree-1) + j);
- }
-
- // Set up the right hand side
- // for the vertical shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points] (1);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- tmp -= local_dofs[i * this->degree + j]
- * this->shape_value_component
- (i * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- 1);
-
- for (unsigned 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]
- (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 < 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 3:
- {
- const QGauss<1>
- reference_edge_quadrature (this->degree);
- const unsigned int &
- n_edge_points = reference_edge_quadrature.size ();
-
- // Let us begin with the
- // interpolation part.
- for (unsigned int i = 0; i < 4; ++i)
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[(i + 8) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[q_point + (i + 8) * n_edge_points] (2);
-
- // 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)
- {
- 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.
- if (std::abs (local_dofs[(i + 4 * j) * this->degree])
- < 1e-14)
- local_dofs[(i + 4 * j) * this->degree] = 0.0;
- }
-
- 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);
-
- // 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;
- }
- }
-
- // If the degree is greater
- // than 0, then we have still
- // some higher order shape
- // functions to consider.
- // Here the projection part
- // starts. The dof values
- // are obtained by solving
- // a linear system of
- // equations.
- if (this->degree > 1)
- {
- // We start with projection
- // on the higher order edge
- // shape function.
- const std::vector<Polynomials::Polynomial<double> > &
- lobatto_polynomials
- = Polynomials::Lobatto::generate_complete_basis
- (this->degree);
- const unsigned int
- line_coordinate[GeometryInfo<3>::lines_per_cell]
- = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
- FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
- FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
- std::vector<Polynomials::Polynomial<double> >
- lobatto_polynomials_grad (this->degree);
-
- for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
- ++i)
- lobatto_polynomials_grad[i]
- = lobatto_polynomials[i + 1].derivative ();
-
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- // Set up the system matrix.
- // This can be used for all
- // edges.
- for (unsigned int i = 0; i < system_matrix.m (); ++i)
- for (unsigned int j = 0; j < system_matrix.n (); ++j)
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- system_matrix (i, j)
- += boundary_weights (q_point, j)
- * lobatto_polynomials_grad[i + 1].value
- (this->generalized_face_support_points[q_point]
- (1));
-
- system_matrix_inv.invert (system_matrix);
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell; ++line)
- {
- // Set up the right hand side.
- system_rhs = 0;
-
- if ((((line == 0) || (line == 1) || (line == 4) ||
- (line == 5)) && (offset + 1 < dim)) ||
- (((line == 2) || (line == 3) || (line == 6) ||
- (line == 7)) && (offset == 0)) || (line > 7))
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- {
- double tmp
- = values[line * n_edge_points + q_point]
- (line_coordinate[line])
- - local_dofs[line * this->degree]
- * this->shape_value_component
- (line * this->degree,
- this->generalized_support_points[line
- * this->degree
- + q_point],
- line_coordinate[line]);
-
- for (unsigned int i = 0; i < system_rhs.size ();
- ++i)
- system_rhs (i)
- += boundary_weights (q_point, i) * tmp;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < solution.size (); ++i)
- if (std::abs (solution (i)) > 1e-14)
- local_dofs[line * this->degree + i + 1]
- = solution (i);
- }
- }
-
- // Then we go on to the
- // face shape functions.
- // Again we set up the
- // system matrix and
- // use it for both, the
- // horizontal and the
- // vertical, shape
- // functions.
- const std::vector<Polynomials::Polynomial<double> > &
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (this->degree-1);
- const unsigned int
- n_face_points = n_edge_points * n_edge_points;
-
- system_matrix.reinit ((this->degree-1) * this->degree,
- (this->degree-1) * this->degree);
- system_matrix = 0;
-
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- for (unsigned int k = 0; k < this->degree; ++k)
- for (unsigned int l = 0; l < this->degree-1; ++l)
- for (unsigned int q_point = 0; q_point < n_face_points;
- ++q_point)
- system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
- += boundary_weights (q_point + n_edge_points,
- 2 * (k * (this->degree-1) + l))
- * legendre_polynomials[i].value
- (this->generalized_face_support_points[q_point
- + 4
- * n_edge_points]
- (0))
- * lobatto_polynomials[j + 2].value
- (this->generalized_face_support_points[q_point
- + 4
- * n_edge_points]
- (1));
-
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.n ());
- system_matrix_inv.invert (system_matrix);
- solution.reinit (system_matrix.m ());
- system_rhs.reinit (system_matrix.m ());
-
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell; ++face)
- {
- switch (face)
- {
- case 0:
- {
- if (offset + 1 < dim)
- {
- // Set up the right hand side
- // for the horizontal shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points] (1);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < this->degree; ++j)
- tmp
- -= local_dofs[4 * i * this->degree
- + j]
- * this->shape_value_component
- (4 * i * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- 1);
-
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < 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);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < this->degree-1; ++j)
- if (std::abs (solution (i * (this->degree-1) + j))
- > 1e-14)
- local_dofs[(i
- + GeometryInfo<dim>::lines_per_cell)
- * (this->degree-1) + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * (this->degree-1) + j);
- }
-
- // Set up the right hand side
- // for the vertical shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_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 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) + 1)
- * tmp;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < this->degree; ++i)
- for (unsigned int j = 0; j < 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 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 < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(4 * i + 1)
- * this->degree + j]
- * this->shape_value_component
- ((4 * i + 1) * this->degree
- + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + n_face_points],
- 1);
-
- for (unsigned int i = 0; i <= 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);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 2 * this->degree) * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
-
- // Set up the right hand side
- // for the vertical shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points + n_face_points]
- (2);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(2 * (i + 4) + 1)
- * this->degree + j]
- * this->shape_value_component
- ((2 * (i + 4) + 1) * this->degree
- + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + n_face_points],
- 2);
-
- for (unsigned 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);
-
- // 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 horizontal shape
- // functions.
- system_rhs = 0;
-
- for (unsigned int q_point = 0;
- q_point < n_face_points; ++q_point)
- {
- double tmp
- = values[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points + 2 * n_face_points]
- (2);
-
- for (unsigned int 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 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);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 4 * this->degree) * deg
- + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
-
- // 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] (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);
-
- // 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 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 + 3 * n_face_points]
- (2);
-
- 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 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);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j))
- > 1e-14)
- local_dofs[(i + GeometryInfo<dim>::lines_per_cell
- + 6 * this->degree) * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
- }
-
- // 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 + 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;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j)) > 1e-14)
- local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
- + 7 * deg) * this->degree]
- = solution (i * deg + j);
-
- break;
- }
-
- 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 <= deg; ++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 <= 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);
-
- // 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
- + 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 <= deg; ++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)
- system_rhs (i * deg + j)
- += boundary_weights
- (q_point + n_edge_points,
- 2 * (i * deg + j) + 1) * tmp;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= 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
- + 9 * deg)
- * this->degree]
- = solution (i * deg + j);
- }
-
- break;
- }
-
- default:
- 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
- + 5 * n_face_points] (0);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(i + 6)
- * this->degree + j]
- * this->shape_value_component
- ((i + 6) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 5
- * n_face_points],
- 0);
-
- 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);
-
- // 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);
- }
-
- // 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 + 5
- * n_face_points] (1);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[(i + 4)
- * this->degree + j]
- * this->shape_value_component
- ((i + 4) * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + 5
- * n_face_points],
- 1);
-
- 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);
-
- // 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
- + 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))
- * 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));
-
- 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] (0);
-
- 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 + (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);
-
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < 4; ++k)
- tmp -= local_dofs[(i + 2 * (k + 2)
- * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- * this->shape_value_component
- ((i + 2 * (k + 2) * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 0);
- }
-
- 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))
- * 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);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- if (std::abs (solution ((i * deg + j) * deg + k))
- > 1e-14)
- local_dofs[((i + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell
- + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell]
- = solution ((i * deg + j) * deg + k);
- }
-
- // Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points] (1);
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- {
- for (unsigned int k = 0; k < 2; ++k)
- tmp -= local_dofs[i + (4 * j + k)
- * this->degree]
- * this->shape_value_component
- (i + (4 * j + k) * this->degree,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 1);
-
- for (unsigned int k = 0; k < deg; ++k)
- tmp -= local_dofs[(i + 2 * j * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell]
- * this->shape_value_component
- ((i + 2 * j * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 1)
- + local_dofs[i + ((2 * j + 9) * deg + k
- + GeometryInfo<dim>::lines_per_cell)
- * this->degree]
- * this->shape_value_component
- (i + ((2 * j + 9) * deg + k
- + GeometryInfo<dim>::lines_per_cell)
- * this->degree,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 1);
- }
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- system_rhs ((i * deg + j) * deg + k)
- += reference_quadrature.weight (q_point) * tmp
- * lobatto_polynomials_grad[i].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (1))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (0))
- * lobatto_polynomials[k + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (2));
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- if (std::abs (solution ((i * deg + j) * deg + k))
- > 1e-14)
- local_dofs[((i + this->degree + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell + 2
- * GeometryInfo<dim>::faces_per_cell)
- * deg + k
- + GeometryInfo<dim>::lines_per_cell]
- = solution ((i * deg + j) * deg + k);
- }
-
- // 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],
- 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);
- }
-
- 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 <= 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;
- }
-
- default:
- Assert (false, ExcNotImplemented ());
- }
-}
-
-
-// Interpolate a function, which is given by
-// its values at the generalized support
-// points in the finite element space on the
-// reference cell.
-// This is done as usual by projection-based
-// interpolation.
-template <int dim>
-void
-FE_Nedelec<dim>::interpolate (std::vector<double> &local_dofs,
- 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 (values[0].size (),
- this->generalized_support_points.size ()));
- Assert (local_dofs.size () == this->dofs_per_cell,
- ExcDimensionMismatch (local_dofs.size (), this->dofs_per_cell));
- std::fill (local_dofs.begin (), local_dofs.end (), 0.0);
-
- switch (dim)
- {
- case 2:
- {
- // Let us begin with the
- // interpolation part.
- const QGauss<dim - 1> reference_edge_quadrature (this->degree);
- const unsigned int &
- n_edge_points = reference_edge_quadrature.size ();
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- {
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- local_dofs[(i + 2 * j) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[1 - j][q_point + (i + 2 * j) * n_edge_points];
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14)
- local_dofs[(i + 2 * j) * this->degree] = 0.0;
- }
-
- // If the degree is greater
- // than 0, then we have still
- // some higher order edge
- // shape functions to
- // consider.
- // Here the projection part
- // starts. The dof values
- // are obtained by solving
- // a linear system of
- // equations.
- if (this->degree-1 > 1)
- {
- // 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);
-
- 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.
- 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. 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 std::vector<Polynomials::Polynomial<double> > &
- legendre_polynomials
- = Polynomials::Legendre::generate_complete_basis (this->degree-1);
-
- 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);
- // 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)
- {
- double tmp
- = values[0][q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points];
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++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 <= deg; ++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]
- (0))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points]
- (1));
- }
-
- solution.reinit (system_matrix.m ());
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= 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) * deg
- + j + GeometryInfo<dim>::lines_per_cell]
- = solution (i * deg + j);
-
- system_rhs = 0;
- // Set up the right hand side
- // for the vertical shape
- // functions.
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[1][q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points];
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++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 i = 0; i <= deg; ++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
- + 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 <= 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
- + deg) * this->degree]
- = solution (i * deg + j);
- }
-
- break;
- }
-
- case 3:
- {
- // Let us begin with the
- // interpolation part.
- const QGauss<1> reference_edge_quadrature (this->degree);
- const unsigned int &
- n_edge_points = reference_edge_quadrature.size ();
-
- for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
- {
- for (unsigned int i = 0; i < 4; ++i)
- local_dofs[(i + 8) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[2][q_point + (i + 8) * n_edge_points];
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- local_dofs[(i + 2 * (2 * j + k)) * this->degree]
- += reference_edge_quadrature.weight (q_point)
- * values[1 - k][q_point + (i + 2 * (2 * j + k))
- * n_edge_points];
- }
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < 4; ++i)
- if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
- local_dofs[(i + 8) * this->degree] = 0.0;
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j < 2; ++j)
- for (unsigned int k = 0; k < 2; ++k)
- if (std::abs (local_dofs[(i + 2 * (2 * j + k)) * this->degree])
- < 1e-14)
- local_dofs[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
-
- // 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);
- FullMatrix<double> system_matrix (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 ();
-
- // Set up the system matrix.
- // This can be used for all
- // edges.
- for (unsigned int i = 0; i < system_matrix.m (); ++i)
- for (unsigned int j = 0; j < system_matrix.n (); ++j)
- for (unsigned int q_point = 0; q_point < n_edge_points;
- ++q_point)
- system_matrix (i, j)
- += boundary_weights (q_point, j)
- * lobatto_polynomials_grad[i + 1].value
- (this->generalized_face_support_points[q_point]
- (1));
-
- FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
-
- system_matrix_inv.invert (system_matrix);
-
- const unsigned int
- line_coordinate[GeometryInfo<3>::lines_per_cell]
- = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
- Vector<double> system_rhs (system_matrix.m ());
- Vector<double> solution (system_rhs.size ());
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell; ++line)
- {
- // Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < this->degree; ++q_point)
- {
- const double tmp
- = values[line_coordinate[line]][line * this->degree
- + q_point]
- - local_dofs[line * this->degree]
- * this->shape_value_component
- (line * this->degree,
- this->generalized_support_points[line
- * this->degree
- + q_point],
- line_coordinate[line]);
-
- for (unsigned int i = 0; i < system_rhs.size (); ++i)
- system_rhs (i) += boundary_weights (q_point, i)
- * tmp;
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i < solution.size (); ++i)
- if (std::abs (solution (i)) > 1e-14)
- local_dofs[line * this->degree + i + 1] = solution (i);
- }
-
- // Then we go on to the
- // face shape functions.
- // Again we set up the
- // system matrix and
- // use it for both, the
- // horizontal and the
- // vertical, shape
- // 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.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}};
- const unsigned int
- 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}
- };
-
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell; ++face)
- {
- // 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[face_coordinates[face][0]][q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points];
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[edge_indices[face][i]
- * this->degree + j]
- * this->shape_value_component
- (edge_indices[face][i] * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- face_coordinates[face][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.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- if (std::abs (solution (i * deg + j)) > 1e-14)
- local_dofs[(2 * face * this->degree + i
- + GeometryInfo<dim>::lines_per_cell) * 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[face_coordinates[face][1]][q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points];
-
- for (int i = 2; i < (int) GeometryInfo<dim>::lines_per_face; ++i)
- for (unsigned int j = 0; j <= deg; ++j)
- tmp -= local_dofs[edge_indices[face][i]
- * this->degree + j]
- * this->shape_value_component
- (edge_indices[face][i] * this->degree + j,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points],
- face_coordinates[face][1]);
-
- 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);
-
- // 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[((2 * face + 1) * deg + j + GeometryInfo<dim>::lines_per_cell)
- * this->degree + i]
- = 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))
- * 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));
-
- system_matrix_inv.reinit (system_matrix.m (),
- system_matrix.m ());
- system_matrix_inv.invert (system_matrix);
- // Set up the right hand side.
- system_rhs.reinit (system_matrix.m ());
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[0][q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points];
-
- 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 + (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);
-
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < 4; ++k)
- tmp -= local_dofs[(i + 2 * (k + 2) * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell]
- * this->shape_value_component
- ((i + 2 * (k + 2) * this->degree
- + GeometryInfo<dim>::lines_per_cell)
- * deg + j
- + GeometryInfo<dim>::lines_per_cell,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 0);
- }
-
- 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))
- * 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));
- }
-
- solution.reinit (system_rhs.size ());
- system_matrix_inv.vmult (solution, system_rhs);
-
- // Add the computed values
- // to the resulting vector
- // only, if they are not
- // too small.
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < deg; ++j)
- for (unsigned int k = 0; k < deg; ++k)
- if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
- local_dofs[((i + 2 * GeometryInfo<dim>::faces_per_cell)
- * deg + j + GeometryInfo<dim>::lines_per_cell
- + 2 * GeometryInfo<dim>::faces_per_cell)
- * deg + k + GeometryInfo<dim>::lines_per_cell]
- = solution ((i * deg + j) * deg + k);
-
- // Set up the right hand side.
- system_rhs = 0;
-
- for (unsigned int q_point = 0; q_point < n_interior_points;
- ++q_point)
- {
- double tmp
- = values[1][q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points];
-
- 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]
- (1))
- * lobatto_polynomials[j + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (0))
- * lobatto_polynomials[k + 2].value
- (this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points]
- (2));
- }
-
- system_matrix_inv.vmult (solution, system_rhs);
-
- // 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[2][q_point + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points];
-
- for (unsigned int i = 0; i <= deg; ++i)
- for (unsigned int j = 0; j < 4; ++j)
- {
- tmp -= local_dofs[i + (j + 8) * this->degree]
- * this->shape_value_component
- (i + (j + 8) * this->degree,
- this->generalized_support_points[q_point
- + GeometryInfo<dim>::lines_per_cell
- * n_edge_points
- + GeometryInfo<dim>::faces_per_cell
- * n_face_points],
- 2);
-
- for (unsigned int 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 <= 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;
- }
-
- default:
- Assert (false, ExcNotImplemented ());
- }
-}
-
-
-
template <int dim>
std::pair<Table<2,bool>, std::vector<unsigned int> >
FE_Nedelec<dim>::get_constant_modes() const
projects / dealii.git / commitdiff