From 37c2f810be21be57288af9038d3c6a26f60d334d Mon Sep 17 00:00:00 2001 From: Guido Kanschat 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] + *
In 3D: *
For any i=0,...,k the * curls of (0,0,xyi+1zk-i), * (xk-iyzi+1,0,0) and * (0,xi+1yk-iz,0) + * + * 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] + * * * * @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 * *

Degrees of freedom

* - * @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. + */ + std::vector > 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. */ - std::vector > test_values; + std::vector > 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 true 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 void compute_embedding_matrices(const FiniteElement &fe, std::vector > > &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 &fe, FullMatrix (&matrices)[GeometryInfo::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::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 #include +#include #include #include #include @@ -171,7 +173,7 @@ PolynomialsBDM::compute (const Point &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::FE_BDM (const unsigned int deg) std::vector(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::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 face_embeddings[GeometryInfo::max_children_per_face]; for (unsigned int i=0; i::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::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::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::faces_per_cell; ++f) { - for (unsigned int i=0; idofs_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; idofs_per_face; ++i) + local_dofs[dbase+i] = values[GeometryInfo::unit_normal_direction[f]][pbase+i]; + pbase += this->dofs_per_face; + } + else { - local_dofs[fbase+i] = values[GeometryInfo::unit_normal_direction[f]][fbase+i]; + 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 (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; ddofs_per_cell - dbase) % dim == 0, ExcInternalError()); + + for (unsigned int d=0; ddofs_per_cell, ExcInternalError()); + Assert (dbase == this->dofs_per_cell, ExcInternalError()); } @@ -215,20 +231,17 @@ template std::vector FE_BDM::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; d2) - { - 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::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::compute_n_pols(deg-1); std::vector dpo(dim+1); dpo[dim-1] = dofs_per_face; @@ -249,9 +262,12 @@ FE_BDM::get_ria_vector (const unsigned int deg) return std::vector(); } - Assert(dim==2, ExcNotImplemented()); const unsigned int dofs_per_cell = PolynomialsBDM::compute_n_pols(deg); - const unsigned int dofs_per_face = deg+1; + const unsigned int dofs_per_face = PolynomialSpace::compute_n_pols(deg); + + Assert(GeometryInfo::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::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 + void + initialize_test_values (std::vector > &test_values, + const Quadrature &quadrature, + const unsigned int deg) + { + PolynomialsP poly(deg); + std::vector > dummy1; + std::vector > dummy2; + + test_values.resize(quadrature.size()); + + for (unsigned int k=0; k + void + initialize_test_values (std::vector > &, + const Quadrature<0> &, + const unsigned int) + {} +} + + template void FE_BDM::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::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 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; kgeneralized_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 cell_points(deg-1); + + // Compute the size of the whole support point set + const unsigned int npoints + = cell_points.size() + GeometryInfo::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 face_points (deg+1); - Assert (face_points.size() == this->dofs_per_face, - ExcInternalError()); - for (unsigned int k=0; kdofs_per_face; ++k) - this->generalized_face_support_points[k] = face_points.point(k); - Quadrature faces = QProjector::project_to_all_faces(face_points); - for (unsigned int k=0; - kdofs_per_face*GeometryInfo::faces_per_cell; - ++k) - this->generalized_support_points[k] = faces.point(k+QProjector - ::DataSetDescriptor::face(0, - true, - false, - false, - this->dofs_per_face)); - - current = this->dofs_per_face*GeometryInfo::faces_per_cell; - } + + Quadrature faces = QProjector::project_to_all_faces(face_points); + for (unsigned int k=0; k < face_points.size()*GeometryInfo::faces_per_cell; ++k) + this->generalized_support_points[k] + = faces.point(k+QProjector + ::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 quadrature(deg); // Remember where interior points start - const unsigned int ibase=current; -// for (unsigned int k=0;kgeneralized_support_points[current] = quadrature.point(current-ibase); - ++current; - } - Assert(current == npoints, ExcInternalError()); - + const unsigned int ibase = face_points.size()*GeometryInfo::faces_per_cell; + for (unsigned int k=0; kgeneralized_support_points[ibase+k] = cell_points.point(k); + } // Finally, compute the values of // the test functions in the // interior quadrature points - PolynomialsP poly(deg-2); - test_values.resize(quadrature.size()); - std::vector > dummy1; - std::vector > dummy2; - - for (unsigned int k=0; k::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 #include #include +#include +#include #include #include #include @@ -123,7 +125,9 @@ namespace result["FE_Q_Hierarchical"] = FEFactoryPointer(new FETools::FEFactory >); result["FE_ABF"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_BDM"] + = FEFactoryPointer(new FETools::FEFactory >); result["FE_RaviartThomas"] = FEFactoryPointer(new FETools::FEFactory >); result["FE_RaviartThomasNodal"] @@ -708,7 +712,8 @@ namespace FETools const FiniteElement &fe, const FEValues &coarse, const Householder &H, - FullMatrix &this_matrix) + FullMatrix &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 &fe, std::vector > &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::n_children(RefinementCase(ref_case)); @@ -851,7 +859,7 @@ namespace FETools { task_group += Threads::new_task (&compute_embedding_for_shape_function, - 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 - (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 &fe, std::vector > > &matrices, - const bool isotropic_only) + const bool isotropic_only, + const double threshold) { Threads::TaskGroup task_group; @@ -893,7 +902,7 @@ namespace FETools for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case, - 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 &fe, FullMatrix (&matrices)[GeometryInfo::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 (const FiniteElement &, - std::vector > > &,bool); + std::vector > > &, 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 (const FiniteElement &, FullMatrix ( &)[GeometryInfo::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 @@ -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 *> &fe_datas) FE_DGQ (1), 1)); deallog << (*fe_datas.rbegin())->get_name() << std::endl; + FE_BDM *bdm1 = new FE_BDM(1); + fe_datas.push_back(bdm1); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + FE_BDM *bdm2 = new FE_BDM(2); + fe_datas.push_back(bdm2); + deallog << (*fe_datas.rbegin())->get_name() << std::endl; + // Hcurl elements FE_Nedelec *ned0 = new FE_Nedelec(0); fe_datas.push_back(ned0); -- 2.39.5