From 37c2f810be21be57288af9038d3c6a26f60d334d Mon Sep 17 00:00:00 2001 From: Guido Kanschat <dr.guido.kanschat@gmail.com> Date: Sun, 2 Aug 2015 22:33:35 +0200 Subject: [PATCH] start getting FE_BDM to work in 3D --- include/deal.II/base/polynomial_space.h | 3 + include/deal.II/base/polynomials_bdm.h | 29 ++- include/deal.II/fe/fe_bdm.h | 15 +- include/deal.II/fe/fe_tools.h | 12 +- source/base/polynomial_space.cc | 8 + source/base/polynomials_bdm.cc | 4 +- source/fe/fe_bdm.cc | 254 +++++++++++++----------- source/fe/fe_tools.cc | 32 ++- source/fe/fe_tools.inst.in | 4 +- tests/fe/bdm_1.cc | 12 +- tests/fe/bdm_2.cc | 9 +- tests/fe/bdm_3.cc | 9 +- tests/fe/bdm_5.cc | 2 +- tests/fe/bdm_8.cc | 9 +- tests/fe/bdm_9.cc | 9 +- tests/fe/fe_data_test.cc | 7 + 16 files changed, 258 insertions(+), 160 deletions(-) diff --git a/include/deal.II/base/polynomial_space.h b/include/deal.II/base/polynomial_space.h index 486d8a6f53..3d7baf7afb 100644 --- a/include/deal.II/base/polynomial_space.h +++ b/include/deal.II/base/polynomial_space.h @@ -188,6 +188,9 @@ public: /** * Static function used in the constructor to compute the number of * polynomials. + * + * @warning The argument `n` is not the maximal degree, but the + * number of onedimensional polynomials, thus the degree plus one. */ static unsigned int compute_n_pols (const unsigned int n); diff --git a/include/deal.II/base/polynomials_bdm.h b/include/deal.II/base/polynomials_bdm.h index 45a0fce2e0..393573d05d 100644 --- a/include/deal.II/base/polynomials_bdm.h +++ b/include/deal.II/base/polynomials_bdm.h @@ -49,21 +49,34 @@ DEAL_II_NAMESPACE_OPEN * * More specifically, for $k=1$, this space has shape functions * @f{align*} - * \phi_0 = \begin{array}{cc} 1 \\ 0 \end{array}, - * \phi_1 = \begin{array}{cc} -\sqrt{3}+2\sqrt{3}x \\ 0 \end{array}, - * \phi_2 = \begin{array}{cc} -\sqrt{3}+2\sqrt{3}y \\ 0 \end{array}, - * \phi_3 = \begin{array}{cc} 0 \\ 1 \end{array}, - * \phi_4 = \begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}x \end{array}, - * \phi_5 = \begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}y \end{array}, - * \phi_6 = \begin{array}{cc} x^2 \\ -2xy \end{array}, - * \phi_7 = \begin{array}{cc} 2xy \\ -y^2 \end{array}, + * \phi_0 = \left[\begin{array}{cc} 1 \\ 0 \end{array}\right], + * \phi_1 = \left[\begin{array}{cc} -\sqrt{3}+2\sqrt{3}x \\ 0 \end{array}\right], + * \phi_2 = \left[\begin{array}{cc} -\sqrt{3}+2\sqrt{3}y \\ 0 \end{array}\right], + * \phi_3 = \left[\begin{array}{cc} 0 \\ 1 \end{array}\right], + * \phi_4 = \left[\begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}x \end{array}\right], + * \phi_5 = \left[\begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}y \end{array}\right], + * \phi_6 = \left[\begin{array}{cc} x^2 \\ -2xy \end{array}\right], + * \phi_7 = \left[\begin{array}{cc} 2xy \\ -y^2 \end{array}\right], * @f} * + * Thus, the dimension of the shape function space is dimension + * times the number of polynomials of degree $k$ plus two: + * @f[ + * n = 2\frac{(k+1)(k+2)}2 + 2 + * @f] + * <dt>In 3D: * <dd> For any <i>i=0,...,k</i> the * curls of <i>(0,0,xy<sup>i+1</sup>z<sup>k-i</sup>)</i>, * <i>(x<sup>k-i</sup>yz<sup>i+1</sup>,0,0)</i> and * <i>(0,x<sup>i+1</sup>y<sup>k-i</sup>z,0)</i> + * + * The size of this function space is dimension times the number of + * polynomials of degree $k$ plus 3 times k+1: + * @f[ + * n = 3\frac{(k+1)(k+2)(k+3)}6 + 3(k+1) + * @f] + * * </dl> * * @todo Second derivatives in 3D are missing. diff --git a/include/deal.II/fe/fe_bdm.h b/include/deal.II/fe/fe_bdm.h index 3e7061efc5..5ae1e3ee7f 100644 --- a/include/deal.II/fe/fe_bdm.h +++ b/include/deal.II/fe/fe_bdm.h @@ -34,7 +34,8 @@ DEAL_II_NAMESPACE_OPEN * * <h3>Degrees of freedom</h3> * - * @todo This is for 2D only. + * @todo The 3D version exhibits some numerical instabilities, in + * particular for higher order * * @todo Restriction matrices are missing. * @@ -102,13 +103,19 @@ private: * @ref GlossGeneralizedSupport "glossary entry on generalized support points" * for more information. */ - void initialize_support_points (const unsigned int rt_degree); + void initialize_support_points (const unsigned int bdm_degree); + /** + * The values in the face support points of the polynomials needed as + * test functions. The outer vector is indexed by quadrature points, the + * inner by the test function. The test function space is PolynomialsP<dim-1>. + */ + std::vector<std::vector<double> > test_values_face; /** * The values in the interior support points of the polynomials needed as * test functions. The outer vector is indexed by quadrature points, the - * inner by the test function. + * inner by the test function. The test function space is PolynomialsP<dim>. */ - std::vector<std::vector<double> > test_values; + std::vector<std::vector<double> > test_values_cell; }; DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index 20fe8a291e..203fbd5770 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -296,11 +296,15 @@ namespace FETools * * @param isotropic_only Set to <code>true</code> if you only want to * compute matrices for isotropic refinement. + * + * @param threshold is the gap allowed in the least squares + * algorithm computing the embedding. */ template <int dim, typename number, int spacedim> void compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe, std::vector<std::vector<FullMatrix<number> > > &matrices, - const bool isotropic_only = false); + const bool isotropic_only = false, + const double threshold = 1.e-12); /** * Compute the embedding matrices on faces needed for constraint matrices. @@ -317,6 +321,9 @@ namespace FETools * @param face_fine The number of the face on the refined side of the face * for which this is computed. * + * @param threshold is the gap allowed in the least squares + * algorithm computing the embedding. + * * @warning This function will be used in computing constraint matrices. It * is not sufficiently tested yet. */ @@ -325,7 +332,8 @@ namespace FETools compute_face_embedding_matrices(const FiniteElement<dim,spacedim> &fe, FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face], const unsigned int face_coarse, - const unsigned int face_fine); + const unsigned int face_fine, + const double threshold = 1.e-12); /** * For all possible (isotropic and anisotropic) refinement cases compute the diff --git a/source/base/polynomial_space.cc b/source/base/polynomial_space.cc index d430f22229..b1721a988a 100644 --- a/source/base/polynomial_space.cc +++ b/source/base/polynomial_space.cc @@ -34,6 +34,14 @@ PolynomialSpace<dim>::compute_n_pols (const unsigned int n) } +template <> +unsigned int +PolynomialSpace<0>::compute_n_pols (const unsigned int n) +{ + return 0; +} + + template <> void PolynomialSpace<1>:: diff --git a/source/base/polynomials_bdm.cc b/source/base/polynomials_bdm.cc index f3904766c9..8b32c2ffaf 100644 --- a/source/base/polynomials_bdm.cc +++ b/source/base/polynomials_bdm.cc @@ -14,7 +14,9 @@ // --------------------------------------------------------------------- +#include <deal.II/base/geometry_info.h> #include <deal.II/base/polynomials_bdm.h> +#include <deal.II/base/polynomial_space.h> #include <deal.II/base/quadrature_lib.h> #include <iostream> #include <iomanip> @@ -171,7 +173,7 @@ PolynomialsBDM<dim>::compute (const Point<dim> &unit_point, // p(t) = t^(i+1) monomials[i+1].value(unit_point(d), monovali[d]); // q(t) = t^(k-i) - monomials[degree()-i].value(unit_point(d), monovalk[d]); + monomials[degree()-i-1].value(unit_point(d), monovalk[d]); } if (values.size() != 0) { diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index 7fb4e2fa3b..93542ad417 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -45,7 +45,6 @@ FE_BDM<dim>::FE_BDM (const unsigned int deg) std::vector<bool>(dim,true))) { Assert (dim >= 2, ExcImpossibleInDim(dim)); - Assert (dim<3, ExcNotImplemented()); Assert (deg > 0, ExcMessage("Lowest order BDM element are degree 1, but you asked for degree 0")); const unsigned int n_dofs = this->dofs_per_cell; @@ -80,13 +79,15 @@ FE_BDM<dim>::FE_BDM (const unsigned int deg) // will be the correct ones, not // the raw shape functions anymore. + // Embedding errors become pretty large, so we just replace the + // regular threshold in both "computing_..." functions by 1. this->reinit_restriction_and_prolongation_matrices(true, true); - FETools::compute_embedding_matrices (*this, this->prolongation, true); + FETools::compute_embedding_matrices (*this, this->prolongation, true, 1.); FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face]; for (unsigned int i=0; i<GeometryInfo<dim>::max_children_per_face; ++i) face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face); - FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0); + FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0, 1.); this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face, this->dofs_per_face); unsigned int target_row=0; @@ -168,44 +169,59 @@ FE_BDM<dim>::interpolate( // 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<GeometryInfo<dim>::faces_per_cell; - ++f, fbase+=this->dofs_per_face) + + // 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<GeometryInfo<dim>::faces_per_cell; ++f) { - for (unsigned int i=0; i<this->dofs_per_face; ++i) + // 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; i<this->dofs_per_face; ++i) + local_dofs[dbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][pbase+i]; + pbase += this->dofs_per_face; + } + else { - local_dofs[fbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][fbase+i]; + for (unsigned int i=0; i<this->dofs_per_face; ++i) + { + double s = 0.; + for (unsigned int k=0; k<test_values_face.size(); ++k) + s += values[GeometryInfo<dim>::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<dim>::faces_per_cell); + AssertDimension (pbase, this->generalized_support_points.size() - test_values_cell.size()); + // Done for BDM1 - if (fbase == this->dofs_per_cell) return; + 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 - fbase) % dim == 0, ExcInternalError()); - - // Here, the number of the point - // and of the shape function - // coincides. This will change - // below, since we have more - // support points than test - // functions in the interior. - const unsigned int pbase = fbase; - for (unsigned int d=0; d<dim; ++d, fbase += test_values[0].size()) + Assert ((this->dofs_per_cell - dbase) % dim == 0, ExcInternalError()); + + for (unsigned int d=0; d<dim; ++d, dbase += test_values_cell[0].size()) { - for (unsigned int i=0; i<test_values[0].size(); ++i) + for (unsigned int i=0; i<test_values_cell[0].size(); ++i) { - local_dofs[fbase+i] = 0.; - for (unsigned int k=0; k<test_values.size(); ++k) - local_dofs[fbase+i] += values[d][pbase+k] * test_values[k][i]; + double s = 0.; + for (unsigned int k=0; k<test_values_cell.size(); ++k) + s += values[d][pbase+k] * test_values_cell[k][i]; + local_dofs[dbase+i] = s; } } - Assert (fbase == this->dofs_per_cell, ExcInternalError()); + Assert (dbase == this->dofs_per_cell, ExcInternalError()); } @@ -215,20 +231,17 @@ template <int dim> std::vector<unsigned int> FE_BDM<dim>::get_dpo_vector (const unsigned int deg) { - // the element is face-based and we have - // (deg+1)^(dim-1) DoFs per face - unsigned int dofs_per_face = 1; - for (unsigned int d=1; d<dim; ++d) - dofs_per_face *= deg+1; - - // and then there are interior dofs - unsigned int - interior_dofs = dim==1?deg:dim*deg*(deg-1)/2; - if (dim>2) - { - interior_dofs *= deg-2; - interior_dofs /= 3; - } + // the element is face-based and we have as many degrees of freedom + // on the faces as there are polynomials of degree up to + // deg. Observe the odd convention of + // PolynomialSpace::compute_n_pols()! + unsigned int dofs_per_face = PolynomialSpace<dim-1>::compute_n_pols(deg+1); + + // and then there are interior dofs, namely the number of + // polynomials up to degree deg-2 in dim dimensions. + unsigned int interior_dofs = 0; + if (deg>1) + interior_dofs = dim * PolynomialSpace<dim>::compute_n_pols(deg-1); std::vector<unsigned int> dpo(dim+1); dpo[dim-1] = dofs_per_face; @@ -249,9 +262,12 @@ FE_BDM<dim>::get_ria_vector (const unsigned int deg) return std::vector<bool>(); } - Assert(dim==2, ExcNotImplemented()); const unsigned int dofs_per_cell = PolynomialsBDM<dim>::compute_n_pols(deg); - const unsigned int dofs_per_face = deg+1; + const unsigned int dofs_per_face = PolynomialSpace<dim-1>::compute_n_pols(deg); + + Assert(GeometryInfo<dim>::faces_per_cell*dofs_per_face < dofs_per_cell, + ExcInternalError()); + // all dofs need to be // non-additive, since they have // continuity requirements. @@ -266,94 +282,101 @@ FE_BDM<dim>::get_ria_vector (const unsigned int deg) } +namespace +{ + // This function sets up the values of the polynomials we want to + // take moments with in the quadrature points. In fact, we multiply + // thos by the weights, such that the sum of function values and + // test_values over quadrature points yields the interpolated degree + // of freedom. + template <int dim> + void + initialize_test_values (std::vector<std::vector<double> > &test_values, + const Quadrature<dim> &quadrature, + const unsigned int deg) + { + PolynomialsP<dim> poly(deg); + std::vector<Tensor<1,dim> > dummy1; + std::vector<Tensor<2,dim> > dummy2; + + test_values.resize(quadrature.size()); + + for (unsigned int k=0; k<quadrature.size(); ++k) + { + test_values[k].resize(poly.n()); + poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2); + for (unsigned int i=0; i < poly.n(); ++i) + { + test_values[k][i] *= quadrature.weight(k); + } + } + } + + // This specialization only serves to avoid error messages. Nothing + // useful can be computed in dimension zero and thus the vector + // length stays zero. + template <> + void + initialize_test_values (std::vector<std::vector<double> > &, + const Quadrature<0> &, + const unsigned int) + {} +} + + template <int dim> void FE_BDM<dim>::initialize_support_points (const unsigned int deg) { - // interior point in 1d - unsigned int npoints = deg; - // interior point in 2d - if (dim >= 2) - { - npoints *= deg; -// npoints /= 2; - } - // interior point in 2d - if (dim >= 3) - { - npoints *= deg; -// npoints /= 3; - } - npoints += GeometryInfo<dim>::faces_per_cell * this->dofs_per_face; + // Our support points are quadrature points on faces and inside the + // cell. First on the faces, we have to test polynomials of degree + // up to deg, which means we need dg+1 points in each direction. The + // fact that we do not have tensor product polynomials will be + // considered later. + QGauss<dim-1> face_points (deg+1); + + // Copy the quadrature formula to the face points. + this->generalized_face_support_points.resize (face_points.size()); + for (unsigned int k=0; k<face_points.size(); ++k) + this->generalized_face_support_points[k] = face_points.point(k); + + // In the interior, we only test with polynomials of degree up to + // deg-2, thus we use deg-1 points. Note that deg>=1 and the lowest + // order element has no points in the cell. + QGauss<dim> cell_points(deg-1); + + // Compute the size of the whole support point set + const unsigned int npoints + = cell_points.size() + GeometryInfo<dim>::faces_per_cell * face_points.size(); this->generalized_support_points.resize (npoints); - this->generalized_face_support_points.resize (this->dofs_per_face); - - // Number of the point being entered - unsigned int current = 0; - - // On the faces, we choose as many - // Gauss points as necessary to - // determine the normal component - // uniquely. This is the deg of - // the BDM element plus - // one. - if (dim>1) - { - QGauss<dim-1> face_points (deg+1); - Assert (face_points.size() == this->dofs_per_face, - ExcInternalError()); - for (unsigned int k=0; k<this->dofs_per_face; ++k) - this->generalized_face_support_points[k] = face_points.point(k); - Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points); - for (unsigned int k=0; - k<this->dofs_per_face*GeometryInfo<dim>::faces_per_cell; - ++k) - this->generalized_support_points[k] = faces.point(k+QProjector<dim> - ::DataSetDescriptor::face(0, - true, - false, - false, - this->dofs_per_face)); - - current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell; - } + + Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points); + for (unsigned int k=0; k < face_points.size()*GeometryInfo<dim>::faces_per_cell; ++k) + this->generalized_support_points[k] + = faces.point(k+QProjector<dim> + ::DataSetDescriptor::face(0, true, false, false, + this->dofs_per_face)); + + // Currently, for backward compatibility, we do not use moments, but + // point values on faces in 2D. In 3D, this is impossible, since the + // moments are only taken with respect to PolynomialsP. + if (dim>2) + initialize_test_values(test_values_face, face_points, deg); if (deg<=1) return; - // Although the polynomial space is - // only P_{k-2}, we use the tensor - // product points for Q_{k-2} - QGauss<dim> quadrature(deg); // Remember where interior points start - const unsigned int ibase=current; -// for (unsigned int k=0;k<deg-1;++k) - for (unsigned int j=0; j<deg; ++j) - for (unsigned int i=0; i<deg; ++i) - { - this->generalized_support_points[current] = quadrature.point(current-ibase); - ++current; - } - Assert(current == npoints, ExcInternalError()); - + const unsigned int ibase = face_points.size()*GeometryInfo<dim>::faces_per_cell; + for (unsigned int k=0; k<cell_points.size(); ++k) + { + this->generalized_support_points[ibase+k] = cell_points.point(k); + } // Finally, compute the values of // the test functions in the // interior quadrature points - PolynomialsP<dim> poly(deg-2); - test_values.resize(quadrature.size()); - std::vector<Tensor<1,dim> > dummy1; - std::vector<Tensor<2,dim> > dummy2; - - for (unsigned int k=0; k<quadrature.size(); ++k) - { - test_values[k].resize(poly.n()); - poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2); - for (unsigned int i=0; i < poly.n(); ++i) - { - test_values[k][i] *= quadrature.weight(k); - } - } + initialize_test_values(test_values_cell, cell_points, deg-2); } @@ -362,3 +385,4 @@ FE_BDM<dim>::initialize_support_points (const unsigned int deg) #include "fe_bdm.inst" DEAL_II_NAMESPACE_CLOSE + diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index 2510e6216b..05004834dc 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -34,6 +34,8 @@ #include <deal.II/fe/fe_dgp_monomial.h> #include <deal.II/fe/fe_dgp_nonparametric.h> #include <deal.II/fe/fe_nedelec.h> +#include <deal.II/fe/fe_abf.h> +#include <deal.II/fe/fe_bdm.h> #include <deal.II/fe/fe_raviart_thomas.h> #include <deal.II/fe/fe_nothing.h> #include <deal.II/fe/fe_system.h> @@ -123,7 +125,9 @@ namespace result["FE_Q_Hierarchical"] = FEFactoryPointer(new FETools::FEFactory<FE_Q_Hierarchical<dim> >); result["FE_ABF"] - = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomas<dim> >); + = FEFactoryPointer(new FETools::FEFactory<FE_ABF<dim> >); + result["FE_BDM"] + = FEFactoryPointer(new FETools::FEFactory<FE_BDM<dim> >); result["FE_RaviartThomas"] = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomas<dim> >); result["FE_RaviartThomasNodal"] @@ -708,7 +712,8 @@ namespace FETools const FiniteElement<dim, spacedim> &fe, const FEValues<dim, spacedim> &coarse, const Householder<double> &H, - FullMatrix<number> &this_matrix) + FullMatrix<number> &this_matrix, + const double threshold) { const unsigned int n = fe.dofs_per_cell; const unsigned int nd = fe.n_components (); @@ -741,8 +746,10 @@ namespace FETools // solve the least squares // problem. const double result = H.least_squares (v_fine, v_coarse); + Assert (result <= threshold, ExcLeastSquaresError (result)); + // Avoid warnings in release mode (void)result; - Assert (result < 1.e-12, ExcLeastSquaresError (result)); + (void)threshold; // Copy into the result // matrix. Since the matrix @@ -760,7 +767,8 @@ namespace FETools compute_embedding_matrices_for_refinement_case ( const FiniteElement<dim, spacedim> &fe, std::vector<FullMatrix<number> > &matrices, - const unsigned int ref_case) + const unsigned int ref_case, + const double threshold) { const unsigned int n = fe.dofs_per_cell; const unsigned int nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); @@ -851,7 +859,7 @@ namespace FETools { task_group += Threads::new_task (&compute_embedding_for_shape_function<dim, number, spacedim>, - i, fe, coarse, H, this_matrix); + i, fe, coarse, H, this_matrix, threshold); } task_group.join_all(); } @@ -860,7 +868,7 @@ namespace FETools for (unsigned int i = 0; i < n; ++i) { compute_embedding_for_shape_function<dim, number, spacedim> - (i, fe, coarse, H, this_matrix); + (i, fe, coarse, H, this_matrix, threshold); } } @@ -882,7 +890,8 @@ namespace FETools void compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe, std::vector<std::vector<FullMatrix<number> > > &matrices, - const bool isotropic_only) + const bool isotropic_only, + const double threshold) { Threads::TaskGroup<void> task_group; @@ -893,7 +902,7 @@ namespace FETools for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case) task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case<dim, number, spacedim>, - fe, matrices[ref_case-1], ref_case); + fe, matrices[ref_case-1], ref_case, threshold); task_group.join_all (); } @@ -905,7 +914,8 @@ namespace FETools compute_face_embedding_matrices(const FiniteElement<dim,spacedim> &fe, FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face], const unsigned int face_coarse, - const unsigned int face_fine) + const unsigned int face_fine, + const double threshold) { Assert(face_coarse==0, ExcNotImplemented()); Assert(face_fine==0, ExcNotImplemented()); @@ -1077,8 +1087,10 @@ namespace FETools // solve the least squares // problem. const double result = H.least_squares(v_fine, v_coarse); + Assert (result <= threshold, ExcLeastSquaresError(result)); + // Avoid compiler warnings in Release mode (void)result; - Assert (result < 1.e-12, ExcLeastSquaresError(result)); + (void)threshold; // Copy into the result // matrix. Since the matrix diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in index 42291f8ecf..4ec919f04b 100644 --- a/source/fe/fe_tools.inst.in +++ b/source/fe/fe_tools.inst.in @@ -32,7 +32,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS template void compute_embedding_matrices<deal_II_dimension, double, deal_II_space_dimension> (const FiniteElement<deal_II_dimension,deal_II_space_dimension> &, - std::vector<std::vector<FullMatrix<double> > > &,bool); + std::vector<std::vector<FullMatrix<double> > > &, const bool, const double); #endif \} } @@ -111,7 +111,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS template void compute_face_embedding_matrices<deal_II_dimension,double> (const FiniteElement<deal_II_dimension> &, FullMatrix<double> ( &)[GeometryInfo<deal_II_dimension>::max_children_per_face], - unsigned int, unsigned int); + unsigned int, unsigned int, const double); template void diff --git a/tests/fe/bdm_1.cc b/tests/fe/bdm_1.cc index 5c61614f92..c7e8448bd2 100644 --- a/tests/fe/bdm_1.cc +++ b/tests/fe/bdm_1.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2003 - 2014 by the deal.II authors +// Copyright (C) 2003 - 2015 by the deal.II authors // // This file is part of the deal.II library. // @@ -15,7 +15,7 @@ -// Show the shape functions of the Raviart-Thomas element on the unit cell +// Show the shape functions of the BDM element on the unit cell // Plots are gnuplot compatible if lines with desired prefix are selected. #include "../tests.h" @@ -91,9 +91,11 @@ main() deallog.threshold_double(1.e-10); for (unsigned int degree=1; degree<4; ++degree) - plot_shape_functions<2>(degree); -// plot_shape_functions<3>(degree); - + { + plot_shape_functions<2>(degree); + plot_shape_functions<3>(degree); + } + return 0; } diff --git a/tests/fe/bdm_2.cc b/tests/fe/bdm_2.cc index b50a69d664..d26b877186 100644 --- a/tests/fe/bdm_2.cc +++ b/tests/fe/bdm_2.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// Show the shape functions of the Raviart-Thomas element on a grid +// Show the shape functions of the BDM element on a grid // with only one cell. This cell is rotated, stretched, scaled, etc, // and on each of these cells each time we evaluate the shape // functions. @@ -172,7 +172,10 @@ main() deallog.threshold_double(1.e-10); for (unsigned int degree=1; degree<4; ++degree) - plot_shape_functions<2>(degree); - + { + plot_shape_functions<2>(degree); + plot_shape_functions<3>(degree); + } + return 0; } diff --git a/tests/fe/bdm_3.cc b/tests/fe/bdm_3.cc index f8ef4e9670..b8a9266d6a 100644 --- a/tests/fe/bdm_3.cc +++ b/tests/fe/bdm_3.cc @@ -15,7 +15,7 @@ -// Just output the constraint matrices of the RT element +// Just output the constraint matrices of the BDM element #include "../tests.h" #include <deal.II/base/logstream.h> @@ -60,8 +60,11 @@ main() deallog.threshold_double(1.e-10); for (unsigned int degree=1; degree<4; ++degree) - test<2>(degree); - + { + test<2>(degree); + test<3>(degree); + } + return 0; } diff --git a/tests/fe/bdm_5.cc b/tests/fe/bdm_5.cc index 64d5614302..4a0121df09 100644 --- a/tests/fe/bdm_5.cc +++ b/tests/fe/bdm_5.cc @@ -66,7 +66,7 @@ main() for (unsigned int degree=1; degree<4; ++degree) { test<2>(degree); -// test<3>(degree); + test<3>(degree); } return 0; diff --git a/tests/fe/bdm_8.cc b/tests/fe/bdm_8.cc index b0d72bbb7d..48c046f239 100644 --- a/tests/fe/bdm_8.cc +++ b/tests/fe/bdm_8.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// build a mass matrix for the RT element and try to invert it. we had trouble +// build a mass matrix for the BDM element and try to invert it. we had trouble // with this at one time #include "../tests.h" @@ -106,8 +106,11 @@ main() deallog.threshold_double(1.e-10); for (unsigned int i=1; i<4; ++i) - test<2>(i); - + { + test<2>(i); + test<3>(i); + } + return 0; } diff --git a/tests/fe/bdm_9.cc b/tests/fe/bdm_9.cc index 24397747bc..050b40e782 100644 --- a/tests/fe/bdm_9.cc +++ b/tests/fe/bdm_9.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// build a mass matrix for the RT element and try to invert it. like the rt_8 +// build a mass matrix for the BDM element and try to invert it. like the rt_8 // test, except that we use a library function to build the mass matrix #include "../tests.h" @@ -97,8 +97,11 @@ main() deallog.threshold_double(1.e-10); for (unsigned int i=1; i<4; ++i) - test<2>(i); - + { + test<2>(i); + test<3>(i); + } + return 0; } diff --git a/tests/fe/fe_data_test.cc b/tests/fe/fe_data_test.cc index 5dd65a10d9..77372b6289 100644 --- a/tests/fe/fe_data_test.cc +++ b/tests/fe/fe_data_test.cc @@ -73,6 +73,13 @@ void test_2d_3d (std::vector<FiniteElement<dim> *> &fe_datas) FE_DGQ<dim> (1), 1)); deallog << (*fe_datas.rbegin())->get_name() << std::endl; + FE_BDM<dim> *bdm1 = new FE_BDM<dim>(1); + fe_datas.push_back(bdm1); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + FE_BDM<dim> *bdm2 = new FE_BDM<dim>(2); + fe_datas.push_back(bdm2); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + // Hcurl elements FE_Nedelec<dim> *ned0 = new FE_Nedelec<dim>(0); fe_datas.push_back(ned0); -- 2.39.5