/**
* 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);
*
* 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.
*
* <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.
*
* @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
*
* @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.
* @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.
*/
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
}
+template <>
+unsigned int
+PolynomialSpace<0>::compute_n_pols (const unsigned int n)
+{
+ return 0;
+}
+
+
template <>
void
PolynomialSpace<1>::
// ---------------------------------------------------------------------
+#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>
// 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)
{
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;
// 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;
// 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());
}
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;
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.
}
+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);
}
#include "fe_bdm.inst"
DEAL_II_NAMESPACE_CLOSE
+
#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>
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"]
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 ();
// 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
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));
{
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();
}
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);
}
}
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;
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 ();
}
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());
// 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
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
\}
}
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
// ---------------------------------------------------------------------
//
-// 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.
//
-// 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"
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;
}
// ---------------------------------------------------------------------
-// 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.
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;
}
-// 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>
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;
}
for (unsigned int degree=1; degree<4; ++degree)
{
test<2>(degree);
-// test<3>(degree);
+ test<3>(degree);
}
return 0;
// ---------------------------------------------------------------------
-// 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"
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;
}
// ---------------------------------------------------------------------
-// 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"
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;
}
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);
projects / dealii.git / commitdiff