From 2fe0c3c7a7f2edf1c723dc2a38ecdccc4b63cb01 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 5 Apr 2017 11:15:06 -0600 Subject: [PATCH] Remove the deprecated FiniteElement::interpolate() function. Also remove all of its implementations. --- .../incompatibilities/20170405Bangerth | 7 + include/deal.II/fe/fe.h | 44 - include/deal.II/fe/fe_abf.h | 10 - include/deal.II/fe/fe_bdm.h | 10 - include/deal.II/fe/fe_dg_vector.h | 10 - include/deal.II/fe/fe_dg_vector.templates.h | 29 - include/deal.II/fe/fe_nedelec.h | 10 - include/deal.II/fe/fe_q_bubbles.h | 10 - include/deal.II/fe/fe_q_hierarchical.h | 25 - include/deal.II/fe/fe_rannacher_turek.h | 9 - include/deal.II/fe/fe_raviart_thomas.h | 20 - source/fe/fe.cc | 69 - source/fe/fe_abf.cc | 136 - source/fe/fe_bdm.cc | 94 - source/fe/fe_nedelec.cc | 2499 ----------------- source/fe/fe_q_base.cc | 1 + source/fe/fe_q_bubbles.cc | 78 - source/fe/fe_q_hierarchical.cc | 35 +- source/fe/fe_rannacher_turek.cc | 69 +- source/fe/fe_raviart_thomas.cc | 82 - source/fe/fe_raviart_thomas_nodal.cc | 103 - 21 files changed, 15 insertions(+), 3335 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20170405Bangerth diff --git a/doc/news/changes/incompatibilities/20170405Bangerth b/doc/news/changes/incompatibilities/20170405Bangerth new file mode 100644 index 0000000000..2fb1492973 --- /dev/null +++ b/doc/news/changes/incompatibilities/20170405Bangerth @@ -0,0 +1,7 @@ +Removed: The FiniteElement::interpolate() function and all of its +implementations in derived classes has been removed. It was +previously already deprecated. Use +FiniteElement::convert_generalized_support_point_values_to_nodal_values() +instead. +
+(Wolfgang Bangerth, 2017/04/05) diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 1d7c58cbe0..dd0663a6a7 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -2008,50 +2008,6 @@ public: convert_generalized_support_point_values_to_nodal_values (const std::vector > &support_point_values, std::vector &nodal_values) const; - /** - * Interpolate a set of scalar values, computed in the generalized support - * points. - * - * @note This function is implemented in FiniteElement for the case that the - * element has support points. In this case, the resulting coefficients are - * just the values in the support points. All other elements must - * reimplement it. - * - * @deprecated Use convert_generalized_support_point_values_to_nodal_values() instead. - */ - virtual - void - interpolate(std::vector &local_dofs, - const std::vector &values) const DEAL_II_DEPRECATED; - - /** - * Interpolate a set of vector values, computed in the generalized support - * points. - * - * Since a finite element often only interpolates part of a vector, - * offset is used to determine the first component of the vector to - * be interpolated. Maybe consider changing your data structures to use the - * next function. - * - * @deprecated Use convert_generalized_support_point_values_to_nodal_values() instead. - */ - virtual - void - interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - /** - * Interpolate a set of vector values, computed in the generalized support - * points. - * - * @deprecated Use convert_generalized_support_point_values_to_nodal_values() instead. - */ - virtual - void - interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - //@} /** diff --git a/include/deal.II/fe/fe_abf.h b/include/deal.II/fe/fe_abf.h index 8a6f0a7c16..eef630f388 100644 --- a/include/deal.II/fe/fe_abf.h +++ b/include/deal.II/fe/fe_abf.h @@ -129,16 +129,6 @@ public: 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; - - virtual void interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - virtual std::size_t memory_consumption () const; virtual FiniteElement *clone() const; diff --git a/include/deal.II/fe/fe_bdm.h b/include/deal.II/fe/fe_bdm.h index b6a1d384dc..c7cc26e1ea 100644 --- a/include/deal.II/fe/fe_bdm.h +++ b/include/deal.II/fe/fe_bdm.h @@ -78,16 +78,6 @@ public: 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; - - virtual void interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - private: /** * Only for internal use. Its full name is @p get_dofs_per_object_vector diff --git a/include/deal.II/fe/fe_dg_vector.h b/include/deal.II/fe/fe_dg_vector.h index 516e14c03c..fdafa6a703 100644 --- a/include/deal.II/fe/fe_dg_vector.h +++ b/include/deal.II/fe/fe_dg_vector.h @@ -80,16 +80,6 @@ public: virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; - virtual void interpolate(std::vector &local_dofs, - const std::vector &values) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - virtual std::size_t memory_consumption () const; private: diff --git a/include/deal.II/fe/fe_dg_vector.templates.h b/include/deal.II/fe/fe_dg_vector.templates.h index ebaf40e9fb..4a67fc7b8f 100644 --- a/include/deal.II/fe/fe_dg_vector.templates.h +++ b/include/deal.II/fe/fe_dg_vector.templates.h @@ -90,35 +90,6 @@ FE_DGVector::has_support_on_face } -template -void -FE_DGVector::interpolate -(std::vector &, - const std::vector &) const -{ - Assert(false, ExcNotImplemented()); -} - - -template -void -FE_DGVector::interpolate -(std::vector & /*local_dofs*/, - const std::vector > & /*values*/, - const unsigned int /*offset*/) const -{ - Assert(false, ExcNotImplemented()); -} - -template -void -FE_DGVector::interpolate -(std::vector & /*local_dofs*/, - const VectorSlice > > & /*values*/) const -{ - Assert(false, ExcNotImplemented()); -} - template std::size_t diff --git a/include/deal.II/fe/fe_nedelec.h b/include/deal.II/fe/fe_nedelec.h index 9ea468c122..480b52e43e 100644 --- a/include/deal.II/fe/fe_nedelec.h +++ b/include/deal.II/fe/fe_nedelec.h @@ -261,16 +261,6 @@ public: 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; - - virtual void interpolate (std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate (std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - /** * Return a list of constant modes of the element. */ diff --git a/include/deal.II/fe/fe_q_bubbles.h b/include/deal.II/fe/fe_q_bubbles.h index 27ca79a0b1..20ac7c0720 100644 --- a/include/deal.II/fe/fe_q_bubbles.h +++ b/include/deal.II/fe/fe_q_bubbles.h @@ -110,16 +110,6 @@ public: 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; - - virtual void interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - /** * Return the matrix interpolating from the given finite element to the * present one. The size of the matrix is then @p dofs_per_cell times diff --git a/include/deal.II/fe/fe_q_hierarchical.h b/include/deal.II/fe/fe_q_hierarchical.h index b11550778d..6b956843e5 100644 --- a/include/deal.II/fe/fe_q_hierarchical.h +++ b/include/deal.II/fe/fe_q_hierarchical.h @@ -688,31 +688,6 @@ public: virtual std::pair, std::vector > get_constant_modes () const; - /** - * This function is not implemented and throws an exception if called. - */ - virtual - void interpolate(std::vector &local_dofs, - const std::vector &values) const DEAL_II_DEPRECATED; - - /** - * This function is not implemented and throws an exception if called. - */ - virtual - void - interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - /** - * This function is not implemented and throws an exception if called. - */ - virtual - void - interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - - protected: /** * @p clone function instead of a copy constructor. diff --git a/include/deal.II/fe/fe_rannacher_turek.h b/include/deal.II/fe/fe_rannacher_turek.h index 80c5b7a2ae..c7739b9952 100644 --- a/include/deal.II/fe/fe_rannacher_turek.h +++ b/include/deal.II/fe/fe_rannacher_turek.h @@ -73,15 +73,6 @@ public: 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; - - virtual void interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; private: /** * Order of this element. diff --git a/include/deal.II/fe/fe_raviart_thomas.h b/include/deal.II/fe/fe_raviart_thomas.h index bb7151844f..ebcab3ef3c 100644 --- a/include/deal.II/fe/fe_raviart_thomas.h +++ b/include/deal.II/fe/fe_raviart_thomas.h @@ -134,16 +134,6 @@ public: 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; - - virtual void interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - /** * Return a list of constant modes of the element. This method is currently * not correctly implemented because it returns ones for all components. @@ -271,16 +261,6 @@ public: 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; - - virtual void interpolate(std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset = 0) const DEAL_II_DEPRECATED; - - virtual void interpolate(std::vector &local_dofs, - const VectorSlice > > &values) const DEAL_II_DEPRECATED; - virtual void get_face_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; diff --git a/source/fe/fe.cc b/source/fe/fe.cc index 6b2bdd88a5..dbcd3c97ff 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -1113,75 +1113,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector -void -FiniteElement::interpolate( - std::vector &local_dofs, - const std::vector &values) const -{ - Assert (has_support_points(), ExcFEHasNoSupportPoints()); - Assert (values.size() == unit_support_points.size(), - ExcDimensionMismatch(values.size(), unit_support_points.size())); - Assert (local_dofs.size() == this->dofs_per_cell, - ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); - Assert (this->n_components() == 1, - ExcDimensionMismatch(this->n_components(), 1)); - - std::copy(values.begin(), values.end(), local_dofs.begin()); -} - - - - -template -void -FiniteElement::interpolate( - std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset) const -{ - Assert (has_support_points(), ExcFEHasNoSupportPoints()); - Assert (values.size() == unit_support_points.size(), - ExcDimensionMismatch(values.size(), unit_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())); - - for (unsigned int i=0; idofs_per_cell; ++i) - { - const std::pair index - = this->system_to_component_index(i); - local_dofs[i] = values[i](offset+index.first); - } -} - - - - -template -void -FiniteElement::interpolate( - std::vector &local_dofs, - const VectorSlice > > &values) const -{ - Assert (has_support_points(), ExcFEHasNoSupportPoints()); - Assert (values[0].size() == unit_support_points.size(), - ExcDimensionMismatch(values.size(), unit_support_points.size())); - Assert (local_dofs.size() == this->dofs_per_cell, - ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); - Assert (values.size() == this->n_components(), - ExcDimensionMismatch(values.size(), this->n_components())); - - for (unsigned int i=0; idofs_per_cell; ++i) - { - const std::pair index - = this->system_to_component_index(i); - local_dofs[i] = values[index.first][i]; - } -} - - template std::size_t diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc index 4ad9967d18..03fa7b2620 100644 --- a/source/fe/fe_abf.cc +++ b/source/fe/fe_abf.cc @@ -570,142 +570,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector -void -FE_ABF::interpolate( - std::vector &, - const std::vector &) const -{ - Assert(false, ExcNotImplemented()); -} - - - -template -void -FE_ABF::interpolate( - std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset) const -{ - 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.); - - 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) - * values[face*n_face_points+k](GeometryInfo::unit_normal_direction[face]+offset); - } - - 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]+offset); - } - } - - // TODO: Check if this "correction" can be removed. - for (unsigned int i=0; i -void -FE_ABF::interpolate( - std::vector &local_dofs, - const VectorSlice > > &values) const -{ - 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.); - - 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) - * values[GeometryInfo::unit_normal_direction[face]][face*n_face_points+k]; - } - - 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]][face*n_face_points+fp]; - } - } - - // TODO: Check if this "correction" can be removed. - for (unsigned int i=0; i std::size_t FE_ABF::memory_consumption () const diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index 557c8fc7c3..3b1909aadf 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -193,100 +193,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector -void -FE_BDM::interpolate( - std::vector &, - const std::vector &) const -{ - Assert(false, ExcNotImplemented()); -} - - -template -void -FE_BDM::interpolate( - std::vector &, - const std::vector > &, - unsigned int) const -{ - Assert(false, ExcNotImplemented()); -} - - - -template -void -FE_BDM::interpolate( - std::vector &local_dofs, - const VectorSlice > > &values) const -{ - AssertDimension (values.size(), dim); - Assert (values[0].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)); - - // 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) - local_dofs[dbase+i] = values[GeometryInfo::unit_normal_direction[f]][pbase+i]; - 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]][pbase+k] * test_values_face[k][i]; - local_dofs[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 std::vector FE_BDM::get_dpo_vector (const unsigned int deg) diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index 39fef89974..e4aedd7233 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -4012,2505 +4012,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector -void FE_Nedelec::interpolate (std::vector &, const std::vector &) 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 -void -FE_Nedelec::interpolate (std::vector &local_dofs, - const std::vector > &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 > & - 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 > - 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 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 system_matrix_inv (this->degree-1, this->degree-1); - - system_matrix_inv.invert (system_matrix); - - Vector system_rhs (system_matrix.m ()); - Vector solution (system_rhs.size ()); - - for (unsigned int line = 0; - line < GeometryInfo::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 reference_quadrature (this->degree); - const std::vector > & - 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::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); - 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::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::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::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)); - } - - 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::lines_per_cell) - * (this->degree-1) + j - + GeometryInfo::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::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::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::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 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::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 > & - 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 system_matrix (this->degree-1, this->degree-1); - FullMatrix system_matrix_inv (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 (); - - Vector system_rhs (system_matrix.m ()); - Vector 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::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 > & - 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::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::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::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::lines_per_cell) - * (this->degree-1) + j - + GeometryInfo::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::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::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::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::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::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::lines_per_cell - + 2 * this->degree) * 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 - = values[q_point - + GeometryInfo::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::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::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::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::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::lines_per_cell - + 4 * this->degree) * 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 - = values[q_point - + GeometryInfo::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::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::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::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::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::lines_per_cell - + 6 * this->degree) * 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 - = values[q_point - + GeometryInfo::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::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::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::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::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::lines_per_cell - + 8 * this->degree) * 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 - = values[q_point - + GeometryInfo::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::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::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::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::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::lines_per_cell - + 10 * this->degree) - * 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 - = values[q_point - + GeometryInfo::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::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::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 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); - 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::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 -= 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::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 -= local_dofs[(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)); - } - - 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::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 - = 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 -= 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::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points], - 1); - - for (unsigned int k = 0; k < deg; ++k) - tmp -= local_dofs[(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) - + local_dofs[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 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::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 - = 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 -= local_dofs[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 -= local_dofs[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 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::faces_per_cell)) - * deg + k - + GeometryInfo::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 -void -FE_Nedelec::interpolate (std::vector &local_dofs, - const VectorSlice > > &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 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 > & - 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 - = 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 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 - = values[0][q_point + GeometryInfo::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::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 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::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 - = values[1][q_point + GeometryInfo::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::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 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::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 > & - 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 - = 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 > & - 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 - = values[face_coordinates[face][0]][q_point - + GeometryInfo::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::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::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 - = values[face_coordinates[face][1]][q_point - + GeometryInfo::lines_per_cell - * n_edge_points]; - - for (int i = 2; i < (int) GeometryInfo::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::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::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 - = values[0][q_point + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::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::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 -= local_dofs[(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) - local_dofs[((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 - = values[1][q_point + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::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::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points], - 1); - - for (unsigned int k = 0; k < deg; ++k) - tmp -= local_dofs[(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) - + local_dofs[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 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::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 - = values[2][q_point + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::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::lines_per_cell - * n_edge_points - + GeometryInfo::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::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 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::faces_per_cell)) - * deg + k - + GeometryInfo::lines_per_cell) - * this->degree] - = solution ((i * deg + j) * deg + k); - } - - break; - } - - default: - Assert (false, ExcNotImplemented ()); - } -} - - - template std::pair, std::vector > FE_Nedelec::get_constant_modes() const diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index cb0b1d848b..097a6e5c23 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -1601,6 +1601,7 @@ FE_Q_Base::get_constant_modes () const } + // explicit instantiations #include "fe_q_base.inst" diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc index da85d736f2..eccf22f2ff 100644 --- a/source/fe/fe_q_bubbles.cc +++ b/source/fe/fe_q_bubbles.cc @@ -356,84 +356,6 @@ convert_generalized_support_point_values_to_nodal_values(const std::vector -void -FE_Q_Bubbles::interpolate(std::vector &local_dofs, - const std::vector &values) const -{ - Assert (values.size() == this->unit_support_points.size(), - ExcDimensionMismatch(values.size(), - this->unit_support_points.size())); - Assert (local_dofs.size() == this->dofs_per_cell, - ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); - Assert (this->n_components() == 1, - ExcDimensionMismatch(this->n_components(), 1)); - - std::copy(values.begin(), values.end(), local_dofs.begin()); - - // 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, - const std::vector > &values, - const unsigned int offset) const -{ - Assert (values.size() == this->unit_support_points.size(), - ExcDimensionMismatch(values.size(), - this->unit_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())); - - for (unsigned int i=0; idofs_per_cell-1; ++i) - { - const std::pair index - = this->system_to_component_index(i); - local_dofs[i] = values[i](offset+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, - const VectorSlice > > &values) const -{ - Assert (values[0].size() == this->unit_support_points.size(), - ExcDimensionMismatch(values.size(), - this->unit_support_points.size())); - Assert (local_dofs.size() == this->dofs_per_cell, - ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); - Assert (values.size() == this->n_components(), - ExcDimensionMismatch(values.size(), this->n_components())); - - for (unsigned int i=0; idofs_per_cell-1; ++i) - { - const std::pair index - = this->system_to_component_index(i); - local_dofs[i] = values[index.first][i]; - } - - // We don't use the bubble functions for local interpolation - for (unsigned int i = 0; i void FE_Q_Bubbles:: diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index 519b600994..0e4e9b9088 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -123,39 +123,6 @@ FE_Q_Hierarchical::get_name () const } -template -void -FE_Q_Hierarchical::interpolate( - std::vector &, - const std::vector &) const -{ - // The default implementation assumes that the FE has a delta property, - // i.e., values at the support points equal the corresponding DoFs. This - // is obviously not the case here. - Assert (false, ExcNotImplemented()); -} - - -template -void -FE_Q_Hierarchical::interpolate( - std::vector &, - const std::vector > &, - unsigned int) const -{ - Assert (false, ExcNotImplemented()); -} - - -template -void -FE_Q_Hierarchical::interpolate( - std::vector &, - const VectorSlice > > &) const -{ - Assert (false, ExcNotImplemented()); -} - template FiniteElement * @@ -164,6 +131,8 @@ FE_Q_Hierarchical::clone() const return new FE_Q_Hierarchical(*this); } + + template void FE_Q_Hierarchical::get_interpolation_matrix(const FiniteElement< dim> &source, diff --git a/source/fe/fe_rannacher_turek.cc b/source/fe/fe_rannacher_turek.cc index 4f7624ac4d..8b1b8e7bfc 100644 --- a/source/fe/fe_rannacher_turek.cc +++ b/source/fe/fe_rannacher_turek.cc @@ -105,34 +105,15 @@ template void FE_RannacherTurek:: convert_generalized_support_point_values_to_nodal_values(const std::vector > &support_point_values, - std::vector &nodal_dofs) const + std::vector &nodal_values) 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, - const std::vector &values) const -{ - AssertDimension(values.size(), this->generalized_support_points.size()); - AssertDimension(local_dofs.size(), this->dofs_per_cell); + AssertDimension(nodal_values.size(), this->dofs_per_cell); const unsigned int q_points_per_face = this->weights.size(); - std::fill(local_dofs.begin(), local_dofs.end(), 0.0); + std::fill(nodal_values.begin(), nodal_values.end(), 0.0); - std::vector::const_iterator value = values.begin(); + std::vector >::const_iterator value = support_point_values.begin(); for (unsigned int face = 0; face < dealii::GeometryInfo::faces_per_cell; ++face) @@ -141,7 +122,7 @@ void FE_RannacherTurek::interpolate( q < q_points_per_face; ++q) { - local_dofs[face] += (*value) * this->weights[q]; + nodal_values[face] += (*value)[0] * this->weights[q]; ++value; } } @@ -149,46 +130,6 @@ void FE_RannacherTurek::interpolate( -template -void FE_RannacherTurek::interpolate( - std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset) const -{ - AssertDimension(values.size(), this->generalized_support_points.size()); - AssertDimension(local_dofs.size(), this->dofs_per_cell); - - // extract component at offset and call scalar version of this function - std::vector scalar_values(values.size()); - for (unsigned int q = 0; q < values.size(); ++q) - { - scalar_values[q] = values[q][offset]; - } - this->interpolate(local_dofs, scalar_values); -} - - - -template -void FE_RannacherTurek::interpolate( - std::vector &local_dofs, - const VectorSlice > > &values) const -{ - AssertDimension(values.size(), 1); - AssertDimension(values[0].size(), this->generalized_support_points.size()); - AssertDimension(local_dofs.size(), this->dofs_per_cell); - - // convert data structure to use scalar version of this function - std::vector scalar_values(values[0].size()); - for (unsigned int q = 0; q < values[0].size(); ++q) - { - scalar_values[q] = values[0][q]; - } - this->interpolate(local_dofs, scalar_values); -} - - - // explicit instantiations #include "fe_rannacher_turek.inst" diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index 772bfb5d01..7964b5a6ad 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -496,88 +496,6 @@ convert_generalized_support_point_values_to_nodal_values(const std::vector -void -FE_RaviartThomas::interpolate( - std::vector &, - const std::vector &) const -{ - Assert(false, ExcNotImplemented()); -} - - -template -void -FE_RaviartThomas::interpolate( - std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset) const -{ - 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.); - - 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) - * values[face*n_face_points+k](GeometryInfo::unit_normal_direction[face]+offset); - } - - 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 -void -FE_RaviartThomas::interpolate( - std::vector &local_dofs, - const VectorSlice > > &values) const -{ - 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.); - - 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) - * values[GeometryInfo::unit_normal_direction[face]][face*n_face_points+k]; - } - - 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 std::size_t FE_RaviartThomas::memory_consumption () const diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 56696f83bd..b5d7f269e8 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -329,109 +329,6 @@ convert_generalized_support_point_values_to_nodal_values(const std::vector -void -FE_RaviartThomasNodal::interpolate( - std::vector &, - const std::vector &) const -{ - Assert(false, ExcNotImplemented()); -} - - -template -void -FE_RaviartThomasNodal::interpolate( - std::vector &local_dofs, - const std::vector > &values, - const unsigned int offset) const -{ - 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())); - - // 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) - { - local_dofs[fbase+i] = values[fbase+i](offset+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( - std::vector &local_dofs, - const VectorSlice > > &values) const -{ - Assert (values.size() == this->n_components(), - ExcDimensionMismatch(values.size(), this->n_components())); - Assert (values[0].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)); - // 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) - { - local_dofs[fbase+i] = values[GeometryInfo::unit_normal_direction[f]][fbase+i]; - } - } - // 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()); -} - - //TODO: There are tests that check that the following few functions don't produce assertion failures, but none that actually check whether they do the right thing. one example for such a test would be to project a function onto an hp space and make sure that the convergence order is correct with regard to the lowest used polynomial degree template -- 2.39.5