From cf3ac8ac32d96d9ac686e9e6b2ec63844cbed6f0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 13 Mar 2017 17:41:06 -0600 Subject: [PATCH] Implement the new FiniteElement::convert_generalized_support_point_values_to_nodal_values interface. This is easily done by copying the implementation of the existing FiniteElement::interpolate() functions in the various derived classes, removing the 'offset' argument, and reordering+renaming the other two arguments. --- include/deal.II/fe/fe_abf.h | 6 + include/deal.II/fe/fe_bdm.h | 6 + include/deal.II/fe/fe_nedelec.h | 6 + include/deal.II/fe/fe_q_bubbles.h | 6 + include/deal.II/fe/fe_q_dg0.h | 6 + include/deal.II/fe/fe_rannacher_turek.h | 6 + include/deal.II/fe/fe_raviart_thomas.h | 11 + source/fe/fe_abf.cc | 65 ++ source/fe/fe_bdm.cc | 71 ++ source/fe/fe_nedelec.cc | 928 ++++++++++++++++++++++++ source/fe/fe_q_bubbles.cc | 28 + source/fe/fe_q_dg0.cc | 27 + source/fe/fe_rannacher_turek.cc | 20 + source/fe/fe_raviart_thomas.cc | 36 + source/fe/fe_raviart_thomas_nodal.cc | 49 ++ 15 files changed, 1271 insertions(+) diff --git a/include/deal.II/fe/fe_abf.h b/include/deal.II/fe/fe_abf.h index 24e3f0b8a2..8a6f0a7c16 100644 --- a/include/deal.II/fe/fe_abf.h +++ b/include/deal.II/fe/fe_abf.h @@ -123,6 +123,12 @@ public: virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; + virtual void interpolate(std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; diff --git a/include/deal.II/fe/fe_bdm.h b/include/deal.II/fe/fe_bdm.h index e64dd95172..b6a1d384dc 100644 --- a/include/deal.II/fe/fe_bdm.h +++ b/include/deal.II/fe/fe_bdm.h @@ -72,6 +72,12 @@ public: virtual FiniteElement *clone () const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; + virtual void interpolate(std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; diff --git a/include/deal.II/fe/fe_nedelec.h b/include/deal.II/fe/fe_nedelec.h index 52a3b0af9d..9ea468c122 100644 --- a/include/deal.II/fe/fe_nedelec.h +++ b/include/deal.II/fe/fe_nedelec.h @@ -255,6 +255,12 @@ public: get_prolongation_matrix (const unsigned int child, const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; + virtual void interpolate (std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; diff --git a/include/deal.II/fe/fe_q_bubbles.h b/include/deal.II/fe/fe_q_bubbles.h index 046d6fd676..27ca79a0b1 100644 --- a/include/deal.II/fe/fe_q_bubbles.h +++ b/include/deal.II/fe/fe_q_bubbles.h @@ -104,6 +104,12 @@ public: */ virtual std::string get_name () const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; + virtual void interpolate(std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; diff --git a/include/deal.II/fe/fe_q_dg0.h b/include/deal.II/fe/fe_q_dg0.h index fa3467301f..ddb27bdbc8 100644 --- a/include/deal.II/fe/fe_q_dg0.h +++ b/include/deal.II/fe/fe_q_dg0.h @@ -258,6 +258,12 @@ public: */ virtual std::string get_name () const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; + virtual void interpolate(std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; diff --git a/include/deal.II/fe/fe_rannacher_turek.h b/include/deal.II/fe/fe_rannacher_turek.h index 5c138583be..80c5b7a2ae 100644 --- a/include/deal.II/fe/fe_rannacher_turek.h +++ b/include/deal.II/fe/fe_rannacher_turek.h @@ -67,6 +67,12 @@ public: virtual std::string get_name() const; virtual FiniteElement *clone() const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; + virtual void interpolate(std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; diff --git a/include/deal.II/fe/fe_raviart_thomas.h b/include/deal.II/fe/fe_raviart_thomas.h index 1e76c5b4f0..bb7151844f 100644 --- a/include/deal.II/fe/fe_raviart_thomas.h +++ b/include/deal.II/fe/fe_raviart_thomas.h @@ -128,6 +128,12 @@ public: virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; + virtual void interpolate(std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; @@ -260,6 +266,11 @@ public: virtual FiniteElement *clone () const; + // documentation inherited from the base class + virtual + void + convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const; virtual void interpolate(std::vector &local_dofs, const std::vector &values) const DEAL_II_DEPRECATED; diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc index 3daec17566..86d0b0b8ac 100644 --- a/source/fe/fe_abf.cc +++ b/source/fe/fe_abf.cc @@ -508,6 +508,68 @@ FE_ABF::has_support_on_face (const unsigned int shape_index, +template +void +FE_ABF:: +convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const +{ + Assert (support_point_values.size() == this->generalized_support_points.size(), + ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size())); + Assert (support_point_values[0].size() == this->n_components(), + ExcDimensionMismatch(support_point_values[0].size(), this->n_components())); + Assert (nodal_values.size() == this->dofs_per_cell, + ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell)); + + std::fill(nodal_values.begin(), nodal_values.end(), 0.); + + const unsigned int n_face_points = boundary_weights.size(0); + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int k=0; kdofs_per_face] += boundary_weights(k,i) + * support_point_values[face*n_face_points+k][GeometryInfo::unit_normal_direction[face]]; + } + + const unsigned int start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; + const unsigned int start_cell_points = GeometryInfo::faces_per_cell*n_face_points; + + for (unsigned int k=0; k::faces_per_cell; ++face) + { + double n_orient = (double) GeometryInfo::unit_normal_orientation[face]; + for (unsigned int fp=0; fp < n_face_points; ++fp) + { + // TODO: Check what the face_orientation, face_flip and face_rotation have to be in 3D + unsigned int k = QProjector::DataSetDescriptor::face (face, false, false, false, n_face_points); + for (unsigned int i=0; i::unit_normal_direction[face]]; + } + } + + // TODO: Check if this "correction" can be removed. + for (unsigned int i=0; i void FE_ABF::interpolate( @@ -580,6 +642,8 @@ FE_ABF::interpolate( local_dofs[start_abf_dofs+i] = 0.0; } + + template void FE_ABF::interpolate( @@ -641,6 +705,7 @@ FE_ABF::interpolate( } + template std::size_t FE_ABF::memory_consumption () const diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index 120d2c8808..557c8fc7c3 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -122,6 +122,77 @@ FE_BDM::clone() const +template +void +FE_BDM:: +convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const +{ + Assert (support_point_values.size() == this->generalized_support_points.size(), + ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size())); + AssertDimension (support_point_values[0].size(), dim); + Assert (nodal_values.size() == this->dofs_per_cell, + ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell)); + + // First do interpolation on faces. There, the component evaluated + // depends on the face direction and orientation. + + // The index of the first dof on this face or the cell + unsigned int dbase = 0; + // The index of the first generalized support point on this face or the cell + unsigned int pbase = 0; + for (unsigned int f = 0; f::faces_per_cell; ++f) + { + // Old version with no moments in 2D. See comment below in + // initialize_support_points() + if (test_values_face.size() == 0) + { + for (unsigned int i=0; idofs_per_face; ++i) + nodal_values[dbase+i] = support_point_values[pbase+i][GeometryInfo::unit_normal_direction[f]]; + pbase += this->dofs_per_face; + } + else + { + for (unsigned int i=0; idofs_per_face; ++i) + { + double s = 0.; + for (unsigned int k=0; k::unit_normal_direction[f]] * test_values_face[k][i]; + nodal_values[dbase+i] = s; + } + pbase += test_values_face.size(); + } + dbase += this->dofs_per_face; + } + + AssertDimension (dbase, this->dofs_per_face * GeometryInfo::faces_per_cell); + AssertDimension (pbase, this->generalized_support_points.size() - test_values_cell.size()); + + // Done for BDM1 + if (dbase == this->dofs_per_cell) return; + + // What's missing are the interior + // degrees of freedom. In each + // point, we take all components of + // the solution. + Assert ((this->dofs_per_cell - dbase) % dim == 0, ExcInternalError()); + + for (unsigned int d=0; ddofs_per_cell, ExcInternalError()); +} + + + template void FE_BDM::interpolate( diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index 761cd84c2c..39fef89974 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -3086,6 +3086,934 @@ FE_Nedelec return this->restriction[refinement_case-1][child]; } + +// 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 +void +FE_Nedelec:: +convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, + std::vector &nodal_values) const +{ + const unsigned int deg = this->degree-1; + Assert (support_point_values.size () == this->generalized_support_points.size (), + ExcDimensionMismatch (support_point_values.size (), + this->generalized_support_points.size ())); + Assert (support_point_values[0].size () == this->n_components (), + ExcDimensionMismatch (support_point_values[0].size (), this->n_components ())); + Assert (nodal_values.size () == this->dofs_per_cell, + ExcDimensionMismatch (nodal_values.size (), this->dofs_per_cell)); + std::fill (nodal_values.begin (), nodal_values.end (), 0.0); + + switch (dim) + { + case 2: + { + // Let us begin with the + // interpolation part. + const QGauss 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) + nodal_values[(i + 2 * j) * this->degree] + += reference_edge_quadrature.weight (q_point) + * support_point_values[q_point + (i + 2 * j) * n_edge_points][1 - j]; + + // Add the computed support_point_values to the resulting vector + // only, if they are not + // too small. + if (std::abs (nodal_values[(i + 2 * j) * this->degree]) < 1e-14) + nodal_values[(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 support_point_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 > & + lobatto_polynomials + = Polynomials::Lobatto::generate_complete_basis + (this->degree); + FullMatrix system_matrix (this->degree-1, this->degree-1); + std::vector > + 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 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 system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int line = 0; + line < GeometryInfo::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 + = support_point_values[line * n_edge_points + q_point][line_coordinate[line]] + - nodal_values[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 support_point_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) + nodal_values[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 reference_quadrature (this->degree); + const unsigned int & + n_interior_points = reference_quadrature.size (); + const std::vector > & + 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::lines_per_cell + * n_edge_points] + (0)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (1)) + * lobatto_polynomials_grad[k].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (0)) + * lobatto_polynomials[l + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::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 + = support_point_values[q_point + GeometryInfo::lines_per_cell + * n_edge_points][0]; + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= nodal_values[(i + 2) * this->degree + j] + * this->shape_value_component + ((i + 2) * this->degree + j, + this->generalized_support_points[q_point + + GeometryInfo::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::lines_per_cell + * n_edge_points] + (0)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (1)); + } + + solution.reinit (system_matrix.m ()); + system_matrix_inv.vmult (solution, system_rhs); + + // Add the computed support_point_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) + nodal_values[(i + GeometryInfo::lines_per_cell) * deg + + j + GeometryInfo::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 + = support_point_values[q_point + GeometryInfo::lines_per_cell + * n_edge_points][1]; + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= nodal_values[i * this->degree + j] + * this->shape_value_component + (i * this->degree + j, + this->generalized_support_points[q_point + + GeometryInfo::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::lines_per_cell + * n_edge_points] + (1)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (0)); + } + + system_matrix_inv.vmult (solution, system_rhs); + + // Add the computed support_point_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) + nodal_values[i + (j + GeometryInfo::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) + nodal_values[(i + 8) * this->degree] + += reference_edge_quadrature.weight (q_point) + * support_point_values[q_point + (i + 8) * n_edge_points][2]; + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + nodal_values[(i + 2 * (2 * j + k)) * this->degree] + += reference_edge_quadrature.weight (q_point) + * support_point_values[q_point + (i + 2 * (2 * j + k)) + * n_edge_points][1 - k]; + } + + // Add the computed support_point_values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i < 4; ++i) + if (std::abs (nodal_values[(i + 8) * this->degree]) < 1e-14) + nodal_values[(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 (nodal_values[(i + 2 * (2 * j + k)) * this->degree]) + < 1e-14) + nodal_values[(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 support_point_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 > & + lobatto_polynomials + = Polynomials::Lobatto::generate_complete_basis + (this->degree); + FullMatrix system_matrix (this->degree-1, this->degree-1); + std::vector > + 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 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 system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int line = 0; + line < GeometryInfo::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 + = support_point_values[line * this->degree + + q_point][line_coordinate[line]] + - nodal_values[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) + nodal_values[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 > & + 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::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 + = support_point_values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points][face_coordinates[face][0]]; + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= nodal_values[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::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 support_point_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) + nodal_values[(2 * face * this->degree + i + + GeometryInfo::lines_per_cell) * deg + + j + GeometryInfo::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 + = support_point_values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points][face_coordinates[face][1]]; + + for (int i = 2; i < (int) GeometryInfo::lines_per_face; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= nodal_values[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::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 support_point_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) + nodal_values[((2 * face + 1) * deg + j + GeometryInfo::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 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::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (1)) + * lobatto_polynomials[k + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (2)) + * lobatto_polynomials_grad[l].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[m + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (1)) + * lobatto_polynomials[n + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::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 + = support_point_values[q_point + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::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 -= nodal_values[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::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 0); + + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < 4; ++k) + tmp -= nodal_values[(i + 2 * (k + 2) * this->degree + + GeometryInfo::lines_per_cell) + * deg + j + + GeometryInfo::lines_per_cell] + * this->shape_value_component + ((i + 2 * (k + 2) * this->degree + + GeometryInfo::lines_per_cell) + * deg + j + + GeometryInfo::lines_per_cell, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::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::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (1)) + * lobatto_polynomials[k + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::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) + nodal_values[((i + 2 * GeometryInfo::faces_per_cell) + * deg + j + GeometryInfo::lines_per_cell + + 2 * GeometryInfo::faces_per_cell) + * deg + k + GeometryInfo::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 + = support_point_values[q_point + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::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 -= nodal_values[i + (4 * j + k) * this->degree] + * this->shape_value_component + (i + (4 * j + k) * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 1); + + for (unsigned int k = 0; k < deg; ++k) + tmp -= nodal_values[(i + 2 * j * this->degree + + GeometryInfo::lines_per_cell) + * deg + k + + GeometryInfo::lines_per_cell] + * this->shape_value_component + ((i + 2 * j * this->degree + + GeometryInfo::lines_per_cell) + * deg + k + + GeometryInfo::lines_per_cell, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 1) + + nodal_values[i + ((2 * j + 9) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree] + * this->shape_value_component + (i + ((2 * j + 9) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::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::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (1)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[k + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (2)); + } + + system_matrix_inv.vmult (solution, system_rhs); + + // Add the computed support_point_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) + nodal_values[((i + this->degree + 2 + * GeometryInfo::faces_per_cell) * deg + + j + GeometryInfo::lines_per_cell + 2 + * GeometryInfo::faces_per_cell) * deg + + k + GeometryInfo::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 + = support_point_values[q_point + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points][2]; + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < 4; ++j) + { + tmp -= nodal_values[i + (j + 8) * this->degree] + * this->shape_value_component + (i + (j + 8) * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 2); + + for (unsigned int k = 0; k < deg; ++k) + tmp -= nodal_values[i + ((2 * j + 1) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree] + * this->shape_value_component + (i + ((2 * j + 1) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::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::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (2)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[k + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (1)); + } + + system_matrix_inv.vmult (solution, system_rhs); + + // Add the computed support_point_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) + nodal_values[i + ((j + 2 * (deg + + GeometryInfo::faces_per_cell)) + * deg + k + + GeometryInfo::lines_per_cell) + * this->degree] + = solution ((i * deg + j) * deg + k); + } + + break; + } + + default: + Assert (false, ExcNotImplemented ()); + } +} + + + + + // Since this is a vector valued element, // we cannot interpolate a scalar function. template diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc index 164864032e..e7314673b9 100644 --- a/source/fe/fe_q_bubbles.cc +++ b/source/fe/fe_q_bubbles.cc @@ -328,6 +328,34 @@ FE_Q_Bubbles::clone() const +template +void +FE_Q_Bubbles:: +convert_generalized_support_point_values_to_nodal_values(const std::vector > &support_point_values, + std::vector &nodal_values) const +{ + Assert (support_point_values.size() == this->unit_support_points.size(), + ExcDimensionMismatch(support_point_values.size(), + this->unit_support_points.size())); + Assert (nodal_values.size() == this->dofs_per_cell, + ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell)); + Assert (support_point_values[0].size() == this->n_components(), + ExcDimensionMismatch(support_point_values[0].size(),this->n_components())); + + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + const std::pair index + = this->system_to_component_index(i); + nodal_values[i] = support_point_values[i](index.first); + } + + // We don't use the bubble functions for local interpolation + for (unsigned int i = 0; i void FE_Q_Bubbles::interpolate(std::vector &local_dofs, diff --git a/source/fe/fe_q_dg0.cc b/source/fe/fe_q_dg0.cc index d708f6b09d..2866943612 100644 --- a/source/fe/fe_q_dg0.cc +++ b/source/fe/fe_q_dg0.cc @@ -174,6 +174,33 @@ FE_Q_DG0::clone() const +template +void +FE_Q_DG0:: +convert_generalized_support_point_values_to_nodal_values(const std::vector > &support_point_values, + std::vector &nodal_dofs) const +{ + Assert (support_point_values.size() == this->unit_support_points.size(), + ExcDimensionMismatch(support_point_values.size(), + this->unit_support_points.size())); + Assert (nodal_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(nodal_dofs.size(),this->dofs_per_cell)); + Assert (support_point_values[0].size() == this->n_components(), + ExcDimensionMismatch(support_point_values[0].size(),this->n_components())); + + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + const std::pair index + = this->system_to_component_index(i); + nodal_dofs[i] = support_point_values[i](index.first); + } + + // We don't need the discontinuous function for local interpolation + nodal_dofs[nodal_dofs.size()-1] = 0.; +} + + + template void FE_Q_DG0::interpolate(std::vector &local_dofs, diff --git a/source/fe/fe_rannacher_turek.cc b/source/fe/fe_rannacher_turek.cc index e158d8783e..4f7624ac4d 100644 --- a/source/fe/fe_rannacher_turek.cc +++ b/source/fe/fe_rannacher_turek.cc @@ -101,6 +101,26 @@ void FE_RannacherTurek::initialize_support_points() +template +void +FE_RannacherTurek:: +convert_generalized_support_point_values_to_nodal_values(const std::vector > &support_point_values, + std::vector &nodal_dofs) const +{ + AssertDimension(support_point_values.size(), this->generalized_support_points.size()); + AssertDimension(nodal_dofs.size(), this->dofs_per_cell); + + // extract component and call scalar version of this function + std::vector scalar_values(support_point_values.size()); + for (unsigned int q = 0; q < support_point_values.size(); ++q) + { + scalar_values[q] = support_point_values[q][0]; + } + this->interpolate(nodal_dofs, scalar_values); +} + + + template void FE_RannacherTurek::interpolate( std::vector &local_dofs, diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index 8c3c4a5188..772bfb5d01 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -460,6 +460,42 @@ FE_RaviartThomas::has_support_on_face ( } + +template +void +FE_RaviartThomas:: +convert_generalized_support_point_values_to_nodal_values(const std::vector > &support_point_values, + std::vector &nodal_values) const +{ + Assert (support_point_values.size() == this->generalized_support_points.size(), + ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size())); + Assert (nodal_values.size() == this->dofs_per_cell, + ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell)); + Assert (support_point_values[0].size() == this->n_components(), + ExcDimensionMismatch(support_point_values[0].size(),this->n_components())); + + std::fill(nodal_values.begin(), nodal_values.end(), 0.); + + const unsigned int n_face_points = boundary_weights.size(0); + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int k=0; kdofs_per_face] += boundary_weights(k,i) + * support_point_values[face*n_face_points+k](GeometryInfo::unit_normal_direction[face]); + } + + const unsigned int start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; + const unsigned int start_cell_points = GeometryInfo::faces_per_cell*n_face_points; + + for (unsigned int k=0; k diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index f208e6e623..56696f83bd 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -280,6 +280,55 @@ FE_RaviartThomasNodal::has_support_on_face ( } + +template +void +FE_RaviartThomasNodal:: +convert_generalized_support_point_values_to_nodal_values(const std::vector > &support_point_values, + std::vector &nodal_values) const +{ + Assert (support_point_values.size() == this->generalized_support_points.size(), + ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size())); + Assert (nodal_values.size() == this->dofs_per_cell, + ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell)); + Assert (support_point_values[0].size() == this->n_components(), + ExcDimensionMismatch(support_point_values[0].size(),this->n_components())); + + // First do interpolation on + // faces. There, the component + // evaluated depends on the face + // direction and orientation. + unsigned int fbase = 0; + unsigned int f=0; + for (; f::faces_per_cell; + ++f, fbase+=this->dofs_per_face) + { + for (unsigned int i=0; idofs_per_face; ++i) + { + nodal_values[fbase+i] = support_point_values[fbase+i](GeometryInfo::unit_normal_direction[f]); + } + } + + // The remaining points form dim + // chunks, one for each component. + const unsigned int istep = (this->dofs_per_cell - fbase) / dim; + Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError()); + + f = 0; + while (fbase < this->dofs_per_cell) + { + for (unsigned int i=0; idofs_per_cell, ExcInternalError()); +} + + + template void FE_RaviartThomasNodal::interpolate( -- 2.39.5