// From now on, the shape functions
// will be the correct ones, not
// the raw shape functions anymore.
-
+
// Reinit the vectors of
// restriction and prolongation
// matrices to the right sizes.
// this function returns, so they
// have to be kept in synch
- std::ostringstream namebuf;
+ std::ostringstream namebuf;
namebuf << "FE_Nedelec<" << dim << ">(" << deg << ")";
return namebuf.str();
#if deal_II_dimension == 1
// Set the generalized support
- // points and precompute the
+ // points and precompute the
// parts of the projection-based
// interpolation, which does
// not depend on the interpolated
this->generalized_support_points[q_point + n_boundary_points]
= quadrature.point (q_point);
}
-
+
else
{
// In this case we only need support points
const Quadrature<dim>& edge_quadrature
= QProjector<dim>::project_to_all_faces
(reference_edge_quadrature);
-
+
this->generalized_support_points.resize (n_boundary_points);
for (unsigned int line = 0;
{
// Initialize quadratures to obtain
// quadrature points later on.
- const QGauss<dim - 2> reference_edge_quadrature (degree + 1);
+ const QGauss<1> reference_edge_quadrature (degree + 1);
const unsigned int& n_edge_points = reference_edge_quadrature.size ();
const Quadrature<dim - 1>& edge_quadrature
= QProjector<dim - 1>::project_to_all_faces
this->generalized_support_points[q_point + n_boundary_points]
= quadrature.point (q_point);
}
-
+
else
{
this->generalized_face_support_points.resize (4 * n_edge_points);
}
}
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}
#endif
- // Restriction operator
+ // Restriction operator
template <int dim>
void
FE_Nedelec<dim>::initialize_restriction ()
{
const unsigned int n_boundary_dofs
= GeometryInfo<dim>::lines_per_cell * this->degree;
-
+
for (unsigned int ref = RefinementCase<dim>::cut_x;
ref <= RefinementCase<dim>::isotropic_refinement; ++ref)
{
const unsigned int index = ref - 1;
-
+
switch (ref)
{
case RefinementCase<dim>::cut_x:
this->restriction[index][j] (i + j * this->degree,
i + j * this->degree)
= 2.0;
-
+
for (unsigned int j = 2;
j < GeometryInfo<dim>::lines_per_cell; ++j)
for (unsigned int k = 0; k < 2; ++k)
this->restriction[index][k]
(i + j * this->degree, i + j * this->degree)
= 1.0;
-
+
for (unsigned int j = 0; j < deg; ++j)
for (unsigned int k = 0; k < 2; ++k)
for (unsigned int child = 0;
(i + k * this->degree) * deg + j
+ n_boundary_dofs) = 1.0;
}
-
+
break;
}
-
+
case RefinementCase<dim>::cut_y:
{
for (unsigned int i = 0; i < this->degree; ++i)
this->restriction[index][k]
(i + j * this->degree, i + j * this->degree)
= 1.0;
-
+
this->restriction[index][j]
(i + (j + 2) * this->degree,
i + (j + 2) * this->degree) = 2.0;
}
-
+
for (unsigned int j = 0; j < deg; ++j)
for (unsigned int k = 0; k < 2; ++k)
for (unsigned int child = 0;
(i + k * this->degree) * deg + j
+ n_boundary_dofs) = 1.0;
}
-
+
break;
}
-
+
case RefinementCase<dim>::isotropic_refinement:
{
for (unsigned int i = 0; i < this->degree; ++i)
(i + 3 * this->degree, i + 3 * this->degree)
= 1.0;
}
-
+
for (unsigned int j = 0; j < deg; ++j)
for (unsigned int k = 0; k < 2; ++k)
for (unsigned int child = 0;
(i + k * this->degree) * deg + j
+ n_boundary_dofs) = 0.5;
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}
ref <= RefinementCase<dim>::isotropic_refinement; ++ref)
{
const unsigned int index = ref - 1;
-
+
switch (ref)
{
case RefinementCase<3>::cut_x:
(i + (j + 10) * this->degree,
i + (j + 10) * this->degree) = 2.0;
}
-
+
for (unsigned int i = 0; i < 2 * this->degree * deg; ++i)
for (unsigned int j = 0; j < 2; ++j)
{
this->restriction[index][j]
(i + j * this->degree * deg + n_edge_dofs,
i + j * this->degree + deg + n_edge_dofs) = 2.0;
-
+
for (unsigned int k = 0; k < 4; ++k)
this->restriction[index][j]
(i + (2 * k + 4) * this->degree * deg
i + (2 * k + 4) * this->degree * deg
+ n_edge_dofs) = 1.0;
}
-
+
break;
}
this->restriction[index][j] (i + 5 * this->degree,
i + 5 * this->degree)
= 1.0;
-
+
for (unsigned int k = 3; k < 6; ++k)
this->restriction[index][j]
(i + (j + 2 * k) * this->degree,
i + (j + 2 * k) * this->degree) = 2.0;
}
-
+
for (unsigned int i = 0; i < 2 * this->degree * deg; ++i)
for (unsigned int j = 0; j < 2; ++j)
{
(i + 10 * this->degree * deg + n_edge_dofs,
i + 10 * this->degree * deg + n_edge_dofs) = 1.0;
}
-
+
break;
}
(i + 7 * this->degree, i + 7 * this->degree)
= 1.0;
}
-
+
for (unsigned int j = 0; j < 4; ++j)
this->restriction[index][j]
(i + (j + 8) * this->degree,
i + (j + 8) * this->degree) = 2.0;
}
-
+
for (unsigned int i = 0; i < 2 * this->degree * deg; ++i)
for (unsigned int j = 0; j < 2; ++j)
{
this->restriction[index][j + 2]
(i + 6 * this->degree * deg + n_edge_dofs,
i + 6 * this->degree * deg + n_edge_dofs) = 1.0;
-
+
for (unsigned int k = 0; k < 4; ++k)
this->restriction[index][k]
(i + 2 * (j + 4) * this->degree * deg
i + 2 * (j + 4) * this->degree * deg
+ n_edge_dofs) = 0.5;
}
-
+
break;
}
(i + (j + 8) * this->degree,
i + (j + 8) * this->degree) = 1.0;
}
-
+
for (unsigned int i = 0; i < 2 * this->degree * deg; ++i)
for (unsigned int j = 0; j < 2; ++j)
{
(i + 2 * k * this->degree * deg + n_edge_dofs,
i + 2 * k * this->degree * deg + n_edge_dofs)
= 1.0;
-
+
this->restriction[index][j]
(i + 2 * (j + 4) * this->degree * deg
+ n_edge_dofs,
i + 2 * (j + 4) * this->degree * deg
+ n_edge_dofs) = 2.0;
}
-
+
break;
}
this->restriction[index][j + 2]
(i + (j + 4) * this->degree,
i + (j + 4) * this->degree) = 2.0;
-
+
for (unsigned int k = 0; k < 2; ++k)
{
this->restriction[index][j]
i + (2 * k + 9) * this->degree) = 1.0;
}
}
-
+
for (unsigned int i = 0; i < 2 * this->degree * deg; ++i)
for (unsigned int j = 0; j < 2; ++j)
{
this->restriction[index][2 * j + 1]
(i + 2 * this->degree * deg + n_edge_dofs,
i + 2 * this->degree * deg + n_edge_dofs) = 1.0;
-
+
for (unsigned int k = 0; k < 2; ++k)
{
this->restriction[index][j + 2 * k]
+ n_edge_dofs) = 1.0;
}
}
-
+
break;
}
(i + (k + 2 * (l + 4)) * this->degree,
i + (k + 2 * (l + 4)) * this->degree) = 1.0;
}
-
+
this->restriction[index][j + 2 * k]
(i + (j + 4 * k + 2) * this->degree,
i + (j + 4 * k + 2) * this->degree) = 2.0;
}
-
+
for (unsigned int i = 0; i < 2 * this->degree * deg; ++i)
for (unsigned int j = 0; j < 2; ++j)
{
(i + 2 * j * this->degree * deg + n_edge_dofs,
i + 2 * j * this->degree * deg + n_edge_dofs)
= 0.5;
-
+
for (unsigned int k = 0; k < 2; ++k)
{
this->restriction[index][j + 2 * k]
+ n_edge_dofs) = 1.0;
}
}
-
+
break;
}
(i + (k + 8) * this->degree,
i + (k + 8) * this->degree) = 1.0;
}
-
+
this->restriction[index][2 * (2 * j + 1)]
(i + 10 * this->degree, i + 10 * this->degree)
= 1.0;
(i + 11 * this->degree, i + 11 * this->degree)
= 1.0;
}
-
+
for (unsigned int i = 0; i < 2 * this->degree * deg; ++i)
{
for (unsigned int j = 0; j < 4; ++j)
i + 10 * this->degree * deg + n_edge_dofs)
= 0.5;
}
-
+
for (unsigned int j = 0; j < 2; ++j)
for (unsigned int k = 0; k < 2; ++k)
for (unsigned int l = 0; l < 2; ++l)
i + 2 * (l + 2) * this->degree * deg
+ n_edge_dofs) = 0.5;
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}
-
+
for (unsigned int i = 0; i < 3 * this->degree * deg * deg; ++i)
for (unsigned int child = 0;
child < GeometryInfo<dim>::n_children
= 2.0 / GeometryInfo<dim>::n_children
(RefinementCase<dim> (ref));
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}
// Chech wheter a given shape
// function has support on a
// given face.
-
+
// We just switch through the
// faces of the cell and return
// true, if the shape function
case 0:
if (!((shape_index > deg) && (shape_index < 2 * this->degree)))
return true;
-
+
else
return false;
(shape_index
< GeometryInfo<2>::lines_per_cell * this->degree))
return true;
-
+
else
return false;
case 2:
if (shape_index < 3 * this->degree)
return true;
-
+
else
return false;
if (!((shape_index >= 2 * this->degree) &&
(shape_index < 3 * this->degree)))
return true;
-
+
else
return false;
-
+
default:
{
Assert (false, ExcNotImplemented ());
< (GeometryInfo<3>::lines_per_cell + 11 * deg)
* this->degree)))
return false;
-
+
else
return true;
< (GeometryInfo<3>::lines_per_cell + 12 * deg)
* this->degree)))
return true;
-
+
else
return false;
< (GeometryInfo<3>::lines_per_cell + 11 * deg)
* this->degree)))
return true;
-
+
else
return false;
< (GeometryInfo<3>::lines_per_cell + 11 * deg)
* this->degree)))
return true;
-
+
else
return false;
< (GeometryInfo<3>::lines_per_cell + 10 * deg)
* this->degree)))
return true;
-
+
else
return false;
< (GeometryInfo<3>::lines_per_cell + 12 * deg)
* this->degree)))
return true;
-
+
else
return false;
-
+
default:
{
Assert (false, ExcNotImplemented ());
for (int j = 1; j < (int) GeometryInfo<dim>::lines_per_face; ++j)
interpolation_matrix (j * source_fe.degree + i,
j * this->degree + i) = 1.0;
-
+
for (unsigned int j = 0; j < deg; ++j)
{
interpolation_matrix
(dof, Point<dim> (0.0, quadrature_points[q_point] (0)),
1);
}
-
+
// Add the computed values
// to the interpolation
// matrix only, if they are
assembling_matrix (i, q_point)
= weight * lobatto_polynomials_grad[i + 1].value (tmp);
}
-
+
FullMatrix<double> system_matrix (deg, deg);
-
+
assembling_matrix.mTmult (system_matrix, assembling_matrix);
-
+
FullMatrix<double> system_matrix_inv (deg, deg);
-
+
system_matrix_inv.invert (system_matrix);
-
+
Vector<double> solution (deg);
Vector<double> system_rhs (deg);
-
+
for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
{
system_rhs = 0;
-
+
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
(dof, Point<dim> (0.0,
quadrature_points[q_point] (0)),
1) - interpolation_matrix (0, dof));
-
+
for (unsigned int i = 0; i < deg; ++i)
system_rhs (i)
+= weight
case 3:
{
- const QGauss<dim - 2> reference_edge_quadrature (this->degree);
+ const QGauss<1> reference_edge_quadrature (this->degree);
switch (subface)
{
case 0:
{
- const Quadrature<dim - 2>& edge_quadrature
- = QProjector<dim - 2>::project_to_child
+ const Quadrature<1>& edge_quadrature
+ = QProjector<1>::project_to_child
(reference_edge_quadrature, 0);
const unsigned int n_edge_points = edge_quadrature.size ();
- const std::vector<Point<dim - 2> >&
+ const std::vector<Point<1> >&
edge_quadrature_points = edge_quadrature.get_points ();
// Let us begin with the
0.5 * i, 0.0), 0);
}
}
-
+
// Add the computed values
// to the interpolation
// matrix only, if they are
(i * source_fe.degree, dof)) < 1e-14)
interpolation_matrix (i * source_fe.degree, dof)
= 0.0;
-
+
if (std::abs (interpolation_matrix
((i + 2) * source_fe.degree, dof))
< 1e-14)
FullMatrix<double> system_matrix_inv (deg, deg);
std::vector<Polynomials::Polynomial<double> >
lobatto_polynomials_grad (this->degree);
-
+
for (unsigned int i = 0; i <= deg; ++i)
lobatto_polynomials_grad[i]
= lobatto_polynomials[i + 1].derivative ();
//TODO:[Markus Buerg] We should not need those, since the projections
//on each face should just be copies of each other.
-
+
// Shifted and scaled
// quadrature points on
// the four edges of a
edge_quadrature_points_full_dim
(GeometryInfo<dim>::lines_per_face,
std::vector<Point<dim> > (n_edge_points));
-
+
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
= Point<dim> (edge_quadrature_points[q_point] (0),
0.5, 0.0);
}
-
+
// Set up the system matrix.
// This can be used for all
// edges.
= weight * lobatto_polynomials_grad[i + 1].value
(tmp);
}
-
+
assembling_matrix.mTmult (system_matrix,
assembling_matrix);
system_matrix_inv.invert (system_matrix);
-
+
Vector<double> solution (deg);
Vector<double> system_rhs (deg);
-
+
for (unsigned int dof = 0; dof < this->dofs_per_face;
++dof)
for (unsigned int line = 0;
{
// Set up the right hand side.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_edge_points; ++q_point)
{
(line * source_fe.degree, dof));
const double tmp
= 2.0 * edge_quadrature_points[q_point] (0);
-
+
for (unsigned int i = 0; i < deg; ++i)
system_rhs (i)
+= right_hand_side_value
assembling_matrix.reinit (deg * this->degree,
n_face_points);
-
+
for (unsigned int q_point = 0; q_point < n_face_points;
++q_point)
{
2.0 * face_quadrature_points[q_point] (1));
const double weight
= std::sqrt (face_quadrature.weight (q_point));
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= weight * legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
assembling_matrix (i * deg + j, q_point)
= tmp * lobatto_polynomials[j + 2].value
(quadrature_point (1));
}
}
-
+
system_matrix.reinit (assembling_matrix.m (),
assembling_matrix.m ());
assembling_matrix.mTmult (system_matrix,
system_matrix_inv.invert (system_matrix);
solution.reinit (system_matrix_inv.m ());
system_rhs.reinit (system_matrix_inv.m ());
-
+
for (unsigned int dof = 0; dof < this->dofs_per_face;
++dof)
{
// part on the face shape
// functions. First on the
// horizontal ones, then on
- // the vertical ones.
+ // the vertical ones.
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
* source_fe.shape_value_component
(i * source_fe.degree + j,
quadrature_point, 1);
-
+
right_hand_side_value
*= face_quadrature.weight (q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= right_hand_side_value
* legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= tmp
interpolation_matrix
((i + 4) * source_fe.degree + j - i, dof)
= solution (i * deg + j);
-
+
// Set up the right hand side
// for the vertical shape
// functions.
right_hand_side_value
*= face_quadrature.weight (q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= right_hand_side_value
* legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= tmp
case 1:
{
- const Quadrature<dim - 2>& edge_quadrature_x
- = QProjector<dim - 2>::project_to_child
+ const Quadrature<1>& edge_quadrature_x
+ = QProjector<1>::project_to_child
(reference_edge_quadrature, 1);
- const Quadrature<dim - 2>& edge_quadrature_y
- = QProjector<dim - 2>::project_to_child
+ const Quadrature<1>& edge_quadrature_y
+ = QProjector<1>::project_to_child
(reference_edge_quadrature, 0);
- const std::vector<Point<dim - 2> >&
+ const std::vector<Point<1> >&
edge_quadrature_x_points = edge_quadrature_x.get_points ();
- const std::vector<Point<dim - 2> >&
+ const std::vector<Point<1> >&
edge_quadrature_y_points = edge_quadrature_y.get_points ();
const unsigned int& n_edge_points
= edge_quadrature_x.size ();
0.5 * i, 0.0), 0);
}
}
-
+
// Add the computed values
// to the interpolation
// matrix only, if they are
(i * source_fe.degree, dof)) < 1e-14)
interpolation_matrix (i * source_fe.degree, dof)
= 0.0;
-
+
if (std::abs (interpolation_matrix
((i + 2) * source_fe.degree, dof))
< 1e-14)
i < lobatto_polynomials_grad.size (); ++i)
lobatto_polynomials_grad[i]
= lobatto_polynomials[i + 1].derivative ();
-
+
// Shifted and scaled
// quadrature points and
// weights on the four
edge_quadrature_points_full_dim
(GeometryInfo<dim>::lines_per_face,
std::vector<Point<dim> > (n_edge_points));
-
+
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
edge_quadrature_weights[3][q_point]
= edge_quadrature_weights[2][q_point];
}
-
+
// Set up the system matrix.
// This can be used for all
// edges.
= weight
* lobatto_polynomials_grad[i + 1].value (tmp);
}
-
+
assembling_matrix.mTmult (system_matrix,
assembling_matrix);
system_matrix_inv.invert (system_matrix);
-
+
Vector<double> system_rhs (system_matrix.m ());
Vector<double> solution (system_rhs.size ());
{
// Set up the right hand side.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_edge_points; ++q_point)
{
1) - interpolation_matrix
(line * source_fe.degree,
dof));
-
+
for (unsigned int i = 0; i < deg; ++i)
system_rhs (i)
+= right_hand_side_value
// part on the face shape
// functions. First on the
// horizontal ones, then on
- // the vertical ones.
+ // the vertical ones.
assembling_matrix.reinit (deg * this->degree,
n_face_points);
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
2.0 * face_quadrature_points[q_point] (1));
const double weight
= std::sqrt (face_quadrature.weight (q_point));
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= weight * legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
assembling_matrix (i * deg + j, q_point)
= tmp * lobatto_polynomials[j + 2].value
(quadrature_point (1));
}
}
-
+
system_matrix.reinit (assembling_matrix.m (),
assembling_matrix.m ());
assembling_matrix.mTmult (system_matrix,
system_matrix_inv.invert (system_matrix);
solution.reinit (system_matrix_inv.m ());
system_rhs.reinit (assembling_matrix.m ());
-
+
for (unsigned int dof = 0; dof < this->dofs_per_face;
++dof)
{
* source_fe.shape_value_component
(i * source_fe.degree + j,
quadrature_point, 1);
-
+
right_hand_side_value
*= face_quadrature.weight (q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= right_hand_side_value
* legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
= tmp
interpolation_matrix
((i + 4) * source_fe.degree + j - i, dof)
= solution (i * deg + j);
-
+
// Set up the right hand side
// for the vertical shape
// functions.
right_hand_side_value
*= face_quadrature.weight (q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= right_hand_side_value
* legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= tmp
case 2:
{
- const Quadrature<dim - 2>& edge_quadrature_x
- = QProjector<dim - 2>::project_to_child
+ const Quadrature<1>& edge_quadrature_x
+ = QProjector<1>::project_to_child
(reference_edge_quadrature, 0);
- const Quadrature<dim - 2>& edge_quadrature_y
- = QProjector<dim - 2>::project_to_child
+ const Quadrature<1>& edge_quadrature_y
+ = QProjector<1>::project_to_child
(reference_edge_quadrature, 1);
const unsigned int& n_edge_points
= edge_quadrature_x.size ();
- const std::vector<Point<dim - 2> >&
+ const std::vector<Point<1> >&
edge_quadrature_x_points = edge_quadrature_x.get_points ();
- const std::vector<Point<dim - 2> >&
+ const std::vector<Point<1> >&
edge_quadrature_y_points = edge_quadrature_y.get_points ();
// Let us begin with the
0.5 * (i + 1), 0.0), 0);
}
}
-
+
// Add the computed values
// to the interpolation
// matrix only, if they are
(i * source_fe.degree, dof)) < 1e-14)
interpolation_matrix (i * source_fe.degree, dof)
= 0.0;
-
+
if (std::abs (interpolation_matrix
((i + 2) * source_fe.degree, dof))
< 1e-14)
i < lobatto_polynomials_grad.size (); ++i)
lobatto_polynomials_grad[i]
= lobatto_polynomials[i + 1].derivative ();
-
+
// Shifted and scaled
// quadrature points and
// weights on the four
edge_quadrature_points_full_dim
(GeometryInfo<dim>::lines_per_face,
std::vector<Point<dim> > (n_edge_points));
-
+
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
edge_quadrature_weights[3][q_point]
= edge_quadrature_weights[2][q_point];
}
-
+
// Set up the system matrix.
// This can be used for all
// edges.
* lobatto_polynomials_grad[i + 1].value
(tmp);
}
-
+
assembling_matrix.mTmult (system_matrix,
assembling_matrix);
system_matrix_inv.invert (system_matrix);
-
+
Vector<double> system_rhs (system_matrix.m ());
Vector<double> solution (system_rhs.size ());
{
// Set up the right hand side.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_edge_points; ++q_point)
{
edge_quadrature_points_full_dim[line][q_point],
1) - interpolation_matrix
(line * source_fe.degree, dof));
-
+
for (unsigned int i = 0; i < deg; ++i)
system_rhs (i)
+= right_hand_side_value
assembling_matrix.reinit (deg * this->degree,
n_face_points);
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
const Point<dim - 1> quadrature_point
(2.0 * face_quadrature_points[q_point] (0),
2.0 * face_quadrature_points[q_point] (1) - 1.0);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= weight * legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
assembling_matrix (i * deg + j, q_point)
= tmp * lobatto_polynomials[j + 2].value
(quadrature_point (1));
}
}
-
+
system_matrix.reinit (assembling_matrix.m (),
assembling_matrix.m ());
assembling_matrix.mTmult (system_matrix,
// part on the face shape
// functions. First on the
// horizontal ones, then on
- // the vertical ones.
+ // the vertical ones.
for (unsigned int dof = 0; dof < this->dofs_per_face;
++dof)
{
* source_fe.shape_value_component
(i * source_fe.degree + j,
quadrature_point, 1);
-
+
right_hand_side_value
*= face_quadrature.weight (q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= right_hand_side_value
* legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= tmp
interpolation_matrix
((i + 4) * source_fe.degree + j - i, dof)
= solution (i * deg + j);
-
+
// Set up the right hand side
// for the vertical shape
// functions.
right_hand_side_value *= face_quadrature.weight
(q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= right_hand_side_value
* legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= tmp
case 3:
{
- const Quadrature<dim - 2>& edge_quadrature
- = QProjector<dim - 2>::project_to_child
+ const Quadrature<1>& edge_quadrature
+ = QProjector<1>::project_to_child
(reference_edge_quadrature, 1);
const unsigned int& n_edge_points = edge_quadrature.size ();
- const std::vector<Point<dim - 2> >&
+ const std::vector<Point<1> >&
edge_quadrature_points = edge_quadrature.get_points ();
// Let us begin with the
0.5 * (i + 1), 0.0), 0);
}
}
-
+
// Add the computed values
// to the interpolation
// matrix only, if they are
(i * source_fe.degree, dof)) < 1e-14)
interpolation_matrix (i * source_fe.degree, dof)
= 0.0;
-
+
if (std::abs (interpolation_matrix
((i + 2) * source_fe.degree, dof))
< 1e-14)
i < lobatto_polynomials_grad.size (); ++i)
lobatto_polynomials_grad[i]
= lobatto_polynomials[i + 1].derivative ();
-
+
// Shifted and scaled
// quadrature points on
// the four edges of a
edge_quadrature_points_full_dim
(GeometryInfo<dim>::lines_per_face,
std::vector<Point<dim> > (n_edge_points));
-
+
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
= Point<dim> (edge_quadrature_points[q_point] (0),
1.0, 0.0);
}
-
+
// Set up the system matrix.
// This can be used for all
// edges.
* lobatto_polynomials_grad[i + 1].value
(tmp);
}
-
+
assembling_matrix.mTmult (system_matrix,
assembling_matrix);
system_matrix_inv.invert (system_matrix);
-
+
Vector<double> system_rhs (system_matrix.m ());
Vector<double> solution (system_rhs.size ());
{
// 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
= 2.0 * edge_quadrature_points[q_point] (0)
- 1.0;
-
+
for (unsigned int i = 0; i < deg; ++i)
system_rhs (i)
+= right_hand_side_value
assembling_matrix.reinit (deg * this->degree,
n_face_points);
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
const Point<dim - 1> quadrature_point
(2.0 * face_quadrature_points[q_point] (0) - 1.0,
2.0 * face_quadrature_points[q_point] (1) - 1.0);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= weight * legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
assembling_matrix (i * deg + j, q_point)
= tmp * lobatto_polynomials[j + 2].value
// part on the face shape
// functions. First on the
// horizontal ones, then on
- // the vertical ones.
+ // the vertical ones.
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
* source_fe.shape_value_component
(i * source_fe.degree + j,
quadrature_point, 1);
-
+
right_hand_side_value
*= face_quadrature.weight (q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double tmp
= right_hand_side_value
* legendre_polynomials[i].value
(quadrature_point (0));
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= tmp
interpolation_matrix
((i + 4) * source_fe.degree + j - i, dof)
= solution (i * deg + j);
-
+
// Set up the right hand side
// for the vertical shape
// functions.
right_hand_side_value
*= face_quadrature.weight (q_point);
-
+
for (unsigned int i = 0; i <= deg; ++i)
{
const double L_i
(quadrature_point (0));
const double tmp
= right_hand_side_value * L_i;
-
+
for (unsigned int j = 0; j < deg; ++j)
system_rhs (i * deg + j)
+= tmp
= solution (i * deg + j);
}
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}
{
case 2:
{
- const QGauss<dim - 1> reference_edge_quadrature (this->degree);
+ const QGauss<1> reference_edge_quadrature (this->degree);
const unsigned int& n_edge_points
= reference_edge_quadrature.size ();
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
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;
}
++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_inv (deg, deg);
system_matrix_inv.invert (system_matrix);
-
+
Vector<double> system_rhs (system_matrix.m ());
Vector<double> solution (system_rhs.size ());
{
// Set up the right hand side.
system_rhs = 0;
-
+
for (unsigned int q_point = 0; q_point < n_edge_points;
++q_point)
{
}
system_matrix_inv.vmult (solution, system_rhs);
-
+
// Add the computed values
// to the resulting vector
// only, if they are not
// for the horizontal shape
// functions.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_interior_points; ++q_point)
{
case 3:
{
- const QGauss<dim - 2>
+ const QGauss<1>
reference_edge_quadrature (this->degree);
const unsigned int&
n_edge_points = reference_edge_quadrature.size ();
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
+= 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
+= 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
++i)
lobatto_polynomials_grad[i]
= lobatto_polynomials[i + 1].derivative ();
-
+
Vector<double> system_rhs (system_matrix.m ());
Vector<double> solution (system_rhs.size ());
}
system_matrix_inv.vmult (solution, system_rhs);
-
+
// Add the computed values
// to the resulting vector
// only, if they are not
// for the horizontal shape
// functions.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
+ GeometryInfo<dim>::lines_per_cell]
= solution (i * deg + j);
}
-
+
// Set up the right hand side
// for the vertical shape
// functions.
// for the horizontal shape
// functions.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
+ GeometryInfo<dim>::lines_per_cell]
= solution (i * deg + j);
}
-
+
// Set up the right hand side
// for the vertical shape
// functions.
// for the horizontal shape
// functions.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_face_points;
++q_point)
// for the vertical shape
// functions.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
// for the horizontal shape
// functions.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
+ GeometryInfo<dim>::lines_per_cell]
= solution (i * deg + j);
}
-
+
// Set up the right hand side
// for the vertical shape
// functions.
if (offset == 0)
{
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
+ GeometryInfo<dim>::lines_per_cell]
= solution (i * deg + j);
}
-
+
// Set up the right hand side
// for the vertical shape
// functions.
if (offset == 0)
{
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_face_points; ++q_point)
{
+ GeometryInfo<dim>::lines_per_cell]
= solution (i * deg + j);
}
-
+
// Set up the right hand side
// for the vertical shape
// functions.
{
// Set up the right hand side.
system_rhs = 0;
-
+
for (unsigned int q_point = 0;
q_point < n_interior_points; ++q_point)
{
+ GeometryInfo<dim>::lines_per_cell]
= solution ((i * deg + j) * deg + k);
}
-
+
// Set up the right hand side.
system_rhs = 0;
+ GeometryInfo<dim>::lines_per_cell]
= solution ((i * deg + j) * deg + k);
}
-
+
// Set up the right hand side.
system_rhs = 0;
* this->degree]
= solution ((i * deg + j) * deg + k);
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}
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
FullMatrix<double> system_matrix (deg, deg);
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]
(1));
FullMatrix<double> system_matrix_inv (deg, deg);
-
+
system_matrix_inv.invert (system_matrix);
-
+
const unsigned int
line_coordinate[GeometryInfo<2>::lines_per_cell]
= {1, 1, 0, 0};
}
system_matrix_inv.vmult (solution, system_rhs);
-
+
// Add the computed values
// to the resulting vector
// only, if they are not
const std::vector<Polynomials::Polynomial<double> >&
legendre_polynomials
= Polynomials::Legendre::generate_complete_basis (deg);
-
+
system_matrix.reinit (deg * this->degree, deg * this->degree);
system_matrix = 0;
{
// Let us begin with the
// interpolation part.
- const QGauss<dim - 2> reference_edge_quadrature (this->degree);
+ const QGauss<1> reference_edge_quadrature (this->degree);
const unsigned int&
n_edge_points = reference_edge_quadrature.size ();
* 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
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)
FullMatrix<double> system_matrix (deg, deg);
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.
(1));
FullMatrix<double> system_matrix_inv (deg, deg);
-
+
system_matrix_inv.invert (system_matrix);
const unsigned int
{
// Set up the right hand side.
system_rhs = 0;
-
+
for (unsigned int q_point = 0; q_point <= deg; ++q_point)
{
const double tmp
legendre_polynomials
= Polynomials::Legendre::generate_complete_basis (deg);
const unsigned int n_face_points = n_edge_points * n_edge_points;
-
+
system_matrix.reinit (deg * this->degree, deg * this->degree);
system_matrix = 0;
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}, {0, 2}, {0, 2}, {0, 1}, {0, 1}};
* this->degree]
= solution ((i * deg + j) * deg + k);
}
-
+
break;
}
-
+
default:
Assert (false, ExcNotImplemented ());
}