}
-namespace
+namespace internal
{
- // 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)
+ namespace FE_BDM
{
- PolynomialsP<dim> poly(deg);
- std::vector<Tensor<1,dim> > dummy1;
- std::vector<Tensor<2,dim> > dummy2;
- std::vector<Tensor<3,dim> > dummy3;
- std::vector<Tensor<4,dim> > dummy4;
+ 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;
+ std::vector<Tensor<3,dim> > dummy3;
+ std::vector<Tensor<4,dim> > dummy4;
- test_values.resize(quadrature.size());
+ 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,
- dummy3, dummy4);
- for (unsigned int i=0; i < poly.n(); ++i)
+ for (unsigned int k=0; k<quadrature.size(); ++k)
{
- test_values[k][i] *= quadrature.weight(k);
+ test_values[k].resize(poly.n());
+ poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2,
+ dummy3, dummy4);
+ 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)
- {}
+ // 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)
+ {}
+ }
+ }
}
// 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);
+ internal::FE_BDM::initialize_test_values(test_values_face, face_points, deg);
if (deg<=1) return;
// the test functions in the
// interior quadrature points
- initialize_test_values(test_values_cell, cell_points, deg-2);
+ internal::FE_BDM::initialize_test_values(test_values_cell, cell_points, deg-2);
}
#include "fe_bdm.inst"
DEAL_II_NAMESPACE_CLOSE
-
DEAL_II_NAMESPACE_OPEN
-// namespace for some functions that are used in this file.
-namespace
+namespace internal
{
- // storage of hand-chosen support
- // points
- //
- // For dim=2, dofs_per_cell of
- // FE_DGPMonomial(k) is given by
- // 0.5(k+1)(k+2), i.e.
- //
- // k 0 1 2 3 4 5 6 7
- // dofs 1 3 6 10 15 21 28 36
- //
- // indirect access of unit points:
- // the points for degree k are
- // located at
- //
- // points[start_index[k]..start_index[k+1]-1]
- const unsigned int start_index2d[6]= {0,1,4,10,20,35};
- const double points2d[35][2]=
+ namespace FE_DGPMonomial
{
- {0,0},
- {0,0},{1,0},{0,1},
- {0,0},{1,0},{0,1},{1,1},{0.5,0},{0,0.5},
- {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5},
- {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5}
- };
-
- // For dim=3, dofs_per_cell of
- // FE_DGPMonomial(k) is given by
- // 1./6.(k+1)(k+2)(k+3), i.e.
- //
- // k 0 1 2 3 4 5 6 7
- // dofs 1 4 10 20 35 56 84 120
- const unsigned int start_index3d[6]= {0,1,5,15/*,35*/};
- const double points3d[35][3]=
- {
- {0,0,0},
- {0,0,0},{1,0,0},{0,1,0},{0,0,1},
- {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1}
- };
+ namespace
+ {
+ // storage of hand-chosen support
+ // points
+ //
+ // For dim=2, dofs_per_cell of
+ // FE_DGPMonomial(k) is given by
+ // 0.5(k+1)(k+2), i.e.
+ //
+ // k 0 1 2 3 4 5 6 7
+ // dofs 1 3 6 10 15 21 28 36
+ //
+ // indirect access of unit points:
+ // the points for degree k are
+ // located at
+ //
+ // points[start_index[k]..start_index[k+1]-1]
+ const unsigned int start_index2d[6]= {0,1,4,10,20,35};
+ const double points2d[35][2]=
+ {
+ {0,0},
+ {0,0},{1,0},{0,1},
+ {0,0},{1,0},{0,1},{1,1},{0.5,0},{0,0.5},
+ {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5},
+ {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5}
+ };
+
+ // For dim=3, dofs_per_cell of
+ // FE_DGPMonomial(k) is given by
+ // 1./6.(k+1)(k+2)(k+3), i.e.
+ //
+ // k 0 1 2 3 4 5 6 7
+ // dofs 1 4 10 20 35 56 84 120
+ const unsigned int start_index3d[6]= {0,1,5,15/*,35*/};
+ const double points3d[35][3]=
+ {
+ {0,0,0},
+ {0,0,0},{1,0,0},{0,1,0},{0,0,1},
+ {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1}
+ };
- template <int dim>
- void generate_unit_points (const unsigned int,
- std::vector<Point<dim> > &);
+ template <int dim>
+ void generate_unit_points (const unsigned int,
+
std::vector<Point<dim> > &);
- template <>
- void generate_unit_points (const unsigned int k,
- std::vector<Point<1> > &p)
- {
- Assert(p.size()==k+1, ExcDimensionMismatch(p.size(), k+1));
- const double h = 1./k;
- for (unsigned int i=0; i<p.size(); ++i)
- p[i](0)=i*h;
- }
+ template <>
+ void generate_unit_points (const unsigned int k,
+
std::vector<Point<1> > &p)
+ {
+ Assert(p.size()==k+1, ExcDimensionMismatch(p.size(), k+1));
+ const double h = 1./k;
+
for (unsigned int i=0; i<p.size(); ++i)
+ p[i](0)=i*h;
+ }
- template <>
- void generate_unit_points (const unsigned int k,
- std::vector<Point<2> > &p)
- {
- Assert(k<=4, ExcNotImplemented());
- Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError());
- for (unsigned int i=0; i<p.size(); ++i)
+ template <>
+ void generate_unit_points (const unsigned int k,
+ std::vector<Point<2> > &p)
{
- p[i](0)=points2d[start_index2d[k]+i][0];
- p[i](1)=points2d[start_index2d[k]+i][1];
+ Assert(k<=4, ExcNotImplemented());
+ Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError());
+ for (unsigned int i=0; i<p.size(); ++i)
+ {
+ p[i](0)=points2d[start_index2d[k]+i][0];
+ p[i](1)=points2d[start_index2d[k]+i][1];
+ }
}
- }
- template <>
- void generate_unit_points (const unsigned int k,
- std::vector<Point<3> > &p)
- {
- Assert(k<=2, ExcNotImplemented());
- Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError());
- for (unsigned int i=0; i<p.size(); ++i)
+ template <>
+ void generate_unit_points (const unsigned int k,
+ std::vector<Point<3> > &p)
{
- p[i](0)=points3d[start_index3d[k]+i][0];
- p[i](1)=points3d[start_index3d[k]+i][1];
- p[i](2)=points3d[start_index3d[k]+i][2];
+ Assert(k<=2, ExcNotImplemented());
+ Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError());
+ for (unsigned int i=0; i<p.size(); ++i)
+ {
+ p[i](0)=points3d[start_index3d[k]+i][0];
+ p[i](1)=points3d[start_index3d[k]+i][1];
+ p[i](2)=points3d[start_index3d[k]+i][2];
+ }
}
+ }
}
}
else
{
std::vector<Point<dim> > unit_points(this->dofs_per_cell);
- generate_unit_points(this->degree, unit_points);
+ internal::FE_DGPMonomial::generate_unit_points(this->degree, unit_points);
FullMatrix<double> source_fe_matrix(unit_points.size(), source_fe.dofs_per_cell);
for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
- std::vector<Point<1> >
- get_QGaussLobatto_points (const unsigned int degree)
+ namespace FE_DGQ
{
- if (degree > 0)
- return QGaussLobatto<1>(degree+1).get_points();
- else
- return std::vector<Point<1> >(1, Point<1>(0.5));
+ namespace
+ {
+ std::vector<Point<1> >
+ get_QGaussLobatto_points (const unsigned int degree)
+ {
+ if (degree > 0)
+ return QGaussLobatto<1>(degree+1).get_points();
+ else
+ return std::vector<Point<1> >(1, Point<1>(0.5));
+ }
+ }
}
}
FE_DGQ<dim, spacedim>::FE_DGQ (const unsigned int degree)
:
FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>
- (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
+ (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis
+ (internal::FE_DGQ::get_QGaussLobatto_points(degree))),
FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell, true),
std::vector<ComponentMask>(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(1,true)))
{
// Compute support points, which are the tensor product of the Lagrange
// interpolation points in the constructor.
- Quadrature<dim> support_quadrature(get_QGaussLobatto_points(degree));
+ Quadrature<dim> support_quadrature(internal::FE_DGQ::get_QGaussLobatto_points(degree));
Assert (support_quadrature.get_points().size() > 0,
(typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints ()));
this->unit_support_points = support_quadrature.get_points();
template <int dim, int spacedim>
FE_DGQHermite<dim,spacedim>::FE_DGQHermite (const unsigned int degree)
: FE_DGQ<dim,spacedim>(degree < 3 ?
- Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))
+ Polynomials::generate_complete_Lagrange_basis
+ (internal::FE_DGQ::get_QGaussLobatto_points(degree))
:
Polynomials::HermiteInterpolation::generate_complete_basis(degree))
{}
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
- /**
- * Auxiliary function to create multiplicity vector from input enrichment functions.
- */
- template <typename T>
- std::vector<unsigned int>
- build_multiplicities(const std::vector<std::vector<T > > &functions )
+ namespace FE_Enriched
{
- std::vector<unsigned int> multiplicities;
- multiplicities.push_back(1); // the first one is non-enriched FE
- for (unsigned int i = 0; i < functions.size(); i++)
- multiplicities.push_back(functions[i].size());
-
- return multiplicities;
- }
-
-
- /**
- * Auxiliary function to build FiniteElement's vector
- */
- template <int dim, int spacedim>
- std::vector< const FiniteElement< dim, spacedim > * >
- build_fes(const FiniteElement<dim,spacedim> *fe_base,
- const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched)
- {
- std::vector< const FiniteElement< dim, spacedim > * > fes;
- fes.push_back(fe_base);
- for (unsigned int i = 0; i < fe_enriched.size(); i++)
- fes.push_back(fe_enriched[i]);
-
- return fes;
- }
+ namespace
+ {
+ /**
+ * Auxiliary function to create multiplicity vector from input enrichment functions.
+ */
+ template <typename T>
+ std::vector<unsigned int>
+ build_multiplicities(const std::vector<std::vector<T > > &functions )
+ {
+ std::vector<unsigned int> multiplicities;
+ multiplicities.push_back(1); // the first one is non-enriched FE
+ for (unsigned int i = 0; i < functions.size(); i++)
+ multiplicities.push_back(functions[i].size());
+ return multiplicities;
+ }
- /**
- * Auxiliary function which check consistency of the input parameters.
- * Returns true if everything is ok.
- */
- template <int dim, int spacedim>
- bool
- consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
- const std::vector< unsigned int > &multiplicities,
- const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
- {
- AssertThrow(fes.size() > 0,
- ExcMessage("FEs size should be >=1"));
- AssertThrow(fes.size() == multiplicities.size(),
- ExcMessage("FEs and multiplicities should have the same size"));
- AssertThrow (functions.size() == fes.size() - 1,
- ExcDimensionMismatch(functions.size(), fes.size()-1));
+ /**
+ * Auxiliary function to build FiniteElement's vector
+ */
+ template <int dim, int spacedim>
+ std::vector< const FiniteElement< dim, spacedim > * >
+ build_fes(const FiniteElement<dim,spacedim> *fe_base,
+ const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched)
+ {
+ std::vector< const FiniteElement< dim, spacedim > * > fes;
+ fes.push_back(fe_base);
+ for (unsigned int i = 0; i < fe_enriched.size(); i++)
+ fes.push_back(fe_enriched[i]);
- AssertThrow(multiplicities[0] == 1,
- ExcMessage("First multiplicity should be 1"));
+ return fes;
+ }
- const unsigned int n_comp_base = fes[0]->n_components();
- // start from fe=1 as 0th is always non-enriched FE.
- for (unsigned int fe=1; fe < fes.size(); fe++)
+ /**
+ * Auxiliary function which check consistency of the input parameters.
+ * Returns true if everything is ok.
+ */
+ template <int dim, int spacedim>
+ bool
+ consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes,
+ const std::vector< unsigned int > &multiplicities,
+ const std::vector<std::vector<std::function<const Function<spacedim> *(const typename dealii::Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
{
- const FE_Nothing<dim> *fe_nothing = dynamic_cast<const FE_Nothing<dim>*>(fes[fe]);
- if (fe_nothing)
- AssertThrow (fe_nothing->is_dominating(),
- ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched"));
-
- AssertThrow (fes[fe]->n_components() == n_comp_base,
- ExcMessage("All elements must have the same number of components"));
+ AssertThrow(fes.size() > 0,
+ ExcMessage("FEs size should be >=1"));
+ AssertThrow(fes.size() == multiplicities.size(),
+ ExcMessage("FEs and multiplicities should have the same size"));
+
+ AssertThrow (functions.size() == fes.size() - 1,
+ ExcDimensionMismatch(functions.size(), fes.size()-1));
+
+ AssertThrow(multiplicities[0] == 1,
+ ExcMessage("First multiplicity should be 1"));
+
+ const unsigned int n_comp_base = fes[0]->n_components();
+
+ // start from fe=1 as 0th is always non-enriched FE.
+ for (unsigned int fe=1; fe < fes.size(); fe++)
+ {
+ const FE_Nothing<dim> *fe_nothing = dynamic_cast<const FE_Nothing<dim>*>(fes[fe]);
+ if (fe_nothing)
+ AssertThrow (fe_nothing->is_dominating(),
+ ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched"));
+
+ AssertThrow (fes[fe]->n_components() == n_comp_base,
+ ExcMessage("All elements must have the same number of components"));
+ }
+ return true;
}
- return true;
- }
- /**
- * Auxiliary function which determines whether the FiniteElement will be enriched.
- */
- template <int dim, int spacedim>
- bool
- check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes)
- {
- // start from fe=1 as 0th is always non-enriched FE.
- for (unsigned int fe=1; fe < fes.size(); fe++)
- if (dynamic_cast<const FE_Nothing<dim>*>(fes[fe]) == nullptr)
- // this is not FE_Nothing => there will be enrichment
- return true;
+ /**
+ * Auxiliary function which determines whether the FiniteElement will be enriched.
+ */
+ template <int dim, int spacedim>
+ bool
+ check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes)
+ {
+ // start from fe=1 as 0th is always non-enriched FE.
+ for (unsigned int fe=1; fe < fes.size(); fe++)
+ if (dynamic_cast<const FE_Nothing<dim>*>(fes[fe]) == nullptr)
+ // this is not FE_Nothing => there will be enrichment
+ return true;
- return false;
+ return false;
+ }
+ }
}
}
const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched,
const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)
:
- FE_Enriched<dim,spacedim> (build_fes(fe_base,fe_enriched),
- build_multiplicities(functions),
+ FE_Enriched<dim,spacedim> (internal::FE_Enriched::build_fes(fe_base,fe_enriched),
+ internal::FE_Enriched::build_multiplicities(functions),
functions)
{}
FETools::Compositing::compute_restriction_is_additive_flags(fes,multiplicities),
FETools::Compositing::compute_nonzero_components(fes,multiplicities,false)),
enrichments(functions),
- is_enriched(check_if_enriched(fes)),
+ is_enriched(internal::FE_Enriched::check_if_enriched(fes)),
fe_system(fes,multiplicities)
{
// descriptive error are thrown within the function.
- Assert(consistency_check(fes,multiplicities,functions),
+ Assert(internal::FE_Enriched::consistency_check(fes,multiplicities,functions),
ExcInternalError());
initialize(fes, multiplicities);
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
- std::vector<Point<1> >
- get_QGaussLobatto_points (const unsigned int degree)
+ namespace FE_FaceQ
{
- if (degree > 0)
- return QGaussLobatto<1>(degree+1).get_points();
- else
- return std::vector<Point<1> >(1, Point<1>(0.5));
+ namespace
+ {
+ std::vector<Point<1> >
+ get_QGaussLobatto_points (const unsigned int degree)
+ {
+ if (degree > 0)
+ return QGaussLobatto<1>(degree+1).get_points();
+ else
+ return std::vector<Point<1> >(1, Point<1>(0.5));
+ }
+ }
}
}
template <int dim, int spacedim>
FE_FaceQ<dim,spacedim>::FE_FaceQ (const unsigned int degree)
:
- FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim> (
- TensorProductPolynomials<dim-1>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
- FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
- std::vector<bool>(1,true))
+ FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim>
+ (TensorProductPolynomials<dim-1>(Polynomials::generate_complete_Lagrange_basis(internal::FE_FaceQ::get_QGaussLobatto_points(degree))),
+ FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
+ std::vector<bool>(1,true))
{
// initialize unit face support points
const unsigned int codim = dim-1;
this->unit_face_support_points[0][d] = 0.5;
else
{
- std::vector<Point<1> > points = get_QGaussLobatto_points(degree);
+ std::vector<Point<1> > points = internal::FE_FaceQ::get_QGaussLobatto_points(degree);
unsigned int k=0;
for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
namespace internal
{
- namespace
+ namespace FE_Nedelec
{
- double
- get_embedding_computation_tolerance(const unsigned int p)
+ namespace
{
- // This heuristic was computed by monitoring the worst residual
- // resulting from the least squares computation when computing
- // the face embedding matrices in the FE_Nedelec constructor.
- // The residual growth is exponential, but is bounded by this
- // function up to degree 12.
- return 1.e-15*std::exp(std::pow(p,1.075));
+ double
+ get_embedding_computation_tolerance(const unsigned int p)
+ {
+ // This heuristic was computed by monitoring the worst residual
+ // resulting from the least squares computation when computing
+ // the face embedding matrices in the FE_Nedelec constructor.
+ // The residual growth is exponential, but is bounded by this
+ // function up to degree 12.
+ return 1.e-15*std::exp(std::pow(p,1.075));
+ }
}
}
}
FETools::compute_face_embedding_matrices<dim,double>
(*this, face_embeddings, 0, 0,
- internal::get_embedding_computation_tolerance(order));
+ internal::FE_Nedelec::get_embedding_computation_tolerance(order));
switch (dim)
{
this_nonconst.reinit_restriction_and_prolongation_matrices ();
// Fill prolongation matrices with embedding operators
FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
- internal::get_embedding_computation_tolerance(this->degree));
+ internal::FE_Nedelec::get_embedding_computation_tolerance(this->degree));
#ifdef DEBUG_NEDELEC
deallog << "Restriction" << std::endl;
#endif
this_nonconst.reinit_restriction_and_prolongation_matrices ();
// Fill prolongation matrices with embedding operators
FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
- internal::get_embedding_computation_tolerance(this->degree));
+ internal::FE_Nedelec::get_embedding_computation_tolerance(this->degree));
#ifdef DEBUG_NEDELEC
deallog << "Restriction" << std::endl;
#endif
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
- const char *
- zero_dof_message = "This element has no shape functions.";
+ namespace FE_Nothing
+ {
+ namespace
+ {
+ const char *
+ zero_dof_message = "This element has no shape functions.";
+ }
+ }
}
FE_Nothing<dim,spacedim>::shape_value (const unsigned int /*i*/,
const Point<dim> & /*p*/) const
{
- (void)zero_dof_message;
- Assert(false,ExcMessage(zero_dof_message));
+ (void)internal::FE_Nothing::zero_dof_message;
+ Assert(false,ExcMessage(internal::FE_Nothing::zero_dof_message));
return 0;
}
DEAL_II_NAMESPACE_CLOSE
-
#include <deal.II/fe/fe_poly_tensor.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_cartesian.h>
+#include <deal.II/grid/tria.h>
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
-//---------------------------------------------------------------------------
-// Utility method, which is used to determine the change of sign for
-// the DoFs on the faces of the given cell.
-//---------------------------------------------------------------------------
-
- /**
- * On noncartesian grids, the sign of the DoFs associated with the faces of
- * the elements has to be changed in some cases. This procedure implements an
- * algorithm, which determines the DoFs, which need this sign change for a
- * given cell.
- */
- void
- get_face_sign_change_rt (const Triangulation<1>::cell_iterator &,
- const unsigned int,
- std::vector<double> &face_sign)
+ namespace FE_PolyTensor
{
- // nothing to do in 1d
- std::fill (face_sign.begin (), face_sign.end (), 1.0);
- }
-
-
+ namespace
+ {
+ //---------------------------------------------------------------------------
+ // Utility method, which is used to determine the change of sign for
+ // the DoFs on the faces of the given cell.
+ //---------------------------------------------------------------------------
+
+ /**
+ * On noncartesian grids, the sign of the DoFs associated with the faces of
+ * the elements has to be changed in some cases. This procedure implements an
+ * algorithm, which determines the DoFs, which need this sign change for a
+ * given cell.
+ */
+ void
+ get_face_sign_change_rt (const dealii::Triangulation<1>::cell_iterator &,
+ const unsigned int,
+ std::vector<double> &face_sign)
+ {
+ // nothing to do in 1d
+ std::fill (face_sign.begin (), face_sign.end (), 1.0);
+ }
- void
- get_face_sign_change_rt (const Triangulation<2>::cell_iterator &cell,
- const unsigned int dofs_per_face,
- std::vector<double> &face_sign)
- {
- const unsigned int dim = 2;
- const unsigned int spacedim = 2;
- // Default is no sign
- // change. I.e. multiply by one.
- std::fill (face_sign.begin (), face_sign.end (), 1.0);
- for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2;
- f < GeometryInfo<dim>::faces_per_cell; ++f)
+ void
+ get_face_sign_change_rt (const dealii::Triangulation<2>::cell_iterator &cell,
+ const unsigned int dofs_per_face,
+ std::vector<double> &face_sign)
{
- Triangulation<dim,spacedim>::face_iterator face = cell->face (f);
- if (!face->at_boundary ())
- {
- const unsigned int nn = cell->neighbor_face_no(f);
+ const unsigned int dim = 2;
+ const unsigned int spacedim = 2;
- if (nn < GeometryInfo<dim>::faces_per_cell / 2)
- for (unsigned int j = 0; j < dofs_per_face; ++j)
- {
- Assert (f * dofs_per_face + j < face_sign.size(),
- ExcInternalError());
+ // Default is no sign
+ // change. I.e. multiply by one.
+ std::fill (face_sign.begin (), face_sign.end (), 1.0);
-//TODO: This is probably only going to work for those elements for which all dofs are face dofs
- face_sign[f * dofs_per_face + j] = -1.0;
- }
+ for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2;
+ f < GeometryInfo<dim>::faces_per_cell; ++f)
+ {
+ dealii::Triangulation<dim,spacedim>::face_iterator face = cell->face (f);
+ if (!face->at_boundary ())
+ {
+ const unsigned int nn = cell->neighbor_face_no(f);
+
+ if (nn < GeometryInfo<dim>::faces_per_cell / 2)
+ for (unsigned int j = 0; j < dofs_per_face; ++j)
+ {
+ Assert (f * dofs_per_face + j < face_sign.size(),
+ ExcInternalError());
+
+ //TODO: This is probably only going to work for those elements for which all dofs are face dofs
+ face_sign[f * dofs_per_face + j] = -1.0;
+ }
+ }
}
}
- }
- void
- get_face_sign_change_rt (const Triangulation<3>::cell_iterator &/*cell*/,
- const unsigned int /*dofs_per_face*/,
- std::vector<double> &face_sign)
- {
- std::fill (face_sign.begin (), face_sign.end (), 1.0);
-//TODO: think about what it would take here
- }
+ void
+ get_face_sign_change_rt (const dealii::Triangulation<3>::cell_iterator &/*cell*/,
+ const unsigned int /*dofs_per_face*/,
+ std::vector<double> &face_sign)
+ {
+ std::fill (face_sign.begin (), face_sign.end (), 1.0);
+ //TODO: think about what it would take here
+ }
- void
- get_face_sign_change_nedelec (const Triangulation<1>::cell_iterator &/*cell*/,
- const unsigned int /*dofs_per_face*/,
- std::vector<double> &face_sign)
- {
- // nothing to do in 1d
- std::fill (face_sign.begin (), face_sign.end (), 1.0);
- }
+ void
+ get_face_sign_change_nedelec (const dealii::Triangulation<1>::cell_iterator &/*cell*/,
+ const unsigned int /*dofs_per_face*/,
+ std::vector<double> &face_sign)
+ {
+ // nothing to do in 1d
+ std::fill (face_sign.begin (), face_sign.end (), 1.0);
+ }
- void
- get_face_sign_change_nedelec (const Triangulation<2>::cell_iterator &/*cell*/,
- const unsigned int /*dofs_per_face*/,
- std::vector<double> &face_sign)
- {
- std::fill (face_sign.begin (), face_sign.end (), 1.0);
-//TODO: think about what it would take here
- }
+ void
+ get_face_sign_change_nedelec (const dealii::Triangulation<2>::cell_iterator &/*cell*/,
+ const unsigned int /*dofs_per_face*/,
+ std::vector<double> &face_sign)
+ {
+ std::fill (face_sign.begin (), face_sign.end (), 1.0);
+ //TODO: think about what it would take here
+ }
- void
- get_face_sign_change_nedelec (const Triangulation<3>::cell_iterator &cell,
- const unsigned int /*dofs_per_face*/,
- std::vector<double> &face_sign)
- {
- const unsigned int dim = 3;
- std::fill (face_sign.begin (), face_sign.end (), 1.0);
-//TODO: This is probably only going to work for those elements for which all dofs are face dofs
- for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
- if (!(cell->line_orientation (l)))
- face_sign[l] = -1.0;
+ void
+ get_face_sign_change_nedelec (const dealii::Triangulation<3>::cell_iterator &cell,
+ const unsigned int /*dofs_per_face*/,
+ std::vector<double> &face_sign)
+ {
+ const unsigned int dim = 3;
+ std::fill (face_sign.begin (), face_sign.end (), 1.0);
+ //TODO: This is probably only going to work for those elements for which all dofs are face dofs
+ for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+ if (!(cell->line_orientation (l)))
+ face_sign[l] = -1.0;
+ }
+ }
}
}
std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 );
if (mapping_type == mapping_raviart_thomas)
- get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
+ internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
else if (mapping_type == mapping_nedelec)
- get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
+ internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 );
if (mapping_type == mapping_raviart_thomas)
- get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
+ internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
else if (mapping_type == mapping_nedelec)
- get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
+ internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
{
std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 );
if (mapping_type == mapping_raviart_thomas)
- get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
+ internal::FE_PolyTensor::get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change);
else if (mapping_type == mapping_nedelec)
- get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
+ internal::FE_PolyTensor::get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change);
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
{
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
- std::vector<Point<1> >
- get_QGaussLobatto_points (const unsigned int degree)
+ namespace FE_Q
{
- if (degree > 0)
- return QGaussLobatto<1>(degree+1).get_points();
- else
+ namespace
+ {
+ std::vector<Point<1> >
+ get_QGaussLobatto_points (const unsigned int degree)
{
- typedef FE_Q_Base<TensorProductPolynomials<1>, 1, 1> FEQ;
- AssertThrow(false, FEQ::ExcFEQCannotHaveDegree0());
+ if (degree > 0)
+ return QGaussLobatto<1>(degree+1).get_points();
+ else
+ {
+ typedef FE_Q_Base<TensorProductPolynomials<1>, 1, 1> FEQ;
+ AssertThrow(false, FEQ::ExcFEQCannotHaveDegree0());
+ }
+ return std::vector<Point<1> >();
}
- return std::vector<Point<1> >();
+ }
}
}
FE_Q<dim,spacedim>::FE_Q (const unsigned int degree)
:
FE_Q_Base<TensorProductPolynomials<dim>, dim, spacedim>
- (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
+ (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(internal::FE_Q::get_QGaussLobatto_points(degree))),
FiniteElementData<dim>(this->get_dpo_vector(degree),
1, degree,
FiniteElementData<dim>::H1),
std::vector<bool> (1, false))
{
- this->initialize(get_QGaussLobatto_points(degree));
+ this->initialize(internal::FE_Q::get_QGaussLobatto_points(degree));
}
DEAL_II_NAMESPACE_OPEN
-namespace FE_Q_Helper
+namespace internal
{
- namespace
+ namespace FE_Q_Base
{
- // get the renumbering for faces
- template <int dim>
- inline
- std::vector<unsigned int>
- face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+ namespace
{
- std::vector<unsigned int> dpo(dim, 1U);
- for (unsigned int i=1; i<dpo.size(); ++i)
- dpo[i]=dpo[i-1]*(degree-1);
- const dealii::FiniteElementData<dim-1> face_data(dpo,1,degree);
- std::vector<unsigned int> face_renumber (face_data.dofs_per_cell);
- FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber);
- return face_renumber;
- }
+ // get the renumbering for faces
+ template <int dim>
+ inline
+ std::vector<unsigned int>
+ face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+ {
+ std::vector<unsigned int> dpo(dim, 1U);
+ for (unsigned int i=1; i<dpo.size(); ++i)
+ dpo[i]=dpo[i-1]*(degree-1);
+ const dealii::FiniteElementData<dim-1> face_data(dpo,1,degree);
+ std::vector<unsigned int> face_renumber (face_data.dofs_per_cell);
+ FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber);
+ return face_renumber;
+ }
- // dummy specialization for dim == 1 to avoid linker errors
- template <>
- inline
- std::vector<unsigned int>
- face_lexicographic_to_hierarchic_numbering<1> (const unsigned int)
- {
- return std::vector<unsigned int>();
- }
+ // dummy specialization for dim == 1 to avoid linker errors
+ template <>
+ inline
+ std::vector<unsigned int>
+ face_lexicographic_to_hierarchic_numbering<1> (const unsigned int)
+ {
+ return std::vector<unsigned int>();
+ }
- // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
- // tensorization on inner loops for performance reasons. this clears a
- // dim-array
- template <int dim>
- inline
- void
- zero_indices (unsigned int (&indices)[dim])
- {
- for (unsigned int d=0; d<dim; ++d)
- indices[d] = 0;
- }
+ // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
+ // tensorization on inner loops for performance reasons. this clears a
+ // dim-array
+ template <int dim>
+ inline
+ void
+ zero_indices (unsigned int (&indices)[dim])
+ {
+ for (unsigned int d=0; d<dim; ++d)
+ indices[d] = 0;
+ }
- // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
- // tensorization on inner loops for performance reasons. this increments
- // tensor product indices
- template <int dim>
- inline
- void
- increment_indices (unsigned int (&indices)[dim],
- const unsigned int dofs1d)
- {
- ++indices[0];
- for (int d=0; d<dim-1; ++d)
- if (indices[d]==dofs1d)
- {
- indices[d] = 0;
- indices[d+1]++;
- }
+ // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
+ // tensorization on inner loops for performance reasons. this increments
+ // tensor product indices
+ template <int dim>
+ inline
+ void
+ increment_indices (unsigned int (&indices)[dim],
+ const unsigned int dofs1d)
+ {
+ ++indices[0];
+ for (int d=0; d<dim-1; ++d)
+ if (indices[d]==dofs1d)
+ {
+ indices[d] = 0;
+ indices[d+1]++;
+ }
+ }
}
}
}
const std::vector<unsigned int> &index_map_inverse =
fe.poly_space.get_numbering_inverse();
const std::vector<unsigned int> face_index_map =
- FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
+ internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point<dim>())
- 1.) < 1e-14,
ExcInternalError());
const std::vector<unsigned int> &index_map_inverse =
fe.poly_space.get_numbering_inverse();
const std::vector<unsigned int> face_index_map =
- FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
+ internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point<dim>())
- 1.) < 1e-14,
ExcInternalError());
// find renumbering of faces and assign from values of quadrature
std::vector<unsigned int> face_index_map =
- FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_degree);
+ internal::FE_Q_Base::face_lexicographic_to_hierarchic_numbering<dim>(q_degree);
Quadrature<1> support_1d(points);
Quadrature<codim> support_quadrature(support_1d);
this->unit_face_support_points.resize(support_quadrature.size());
// now expand from 1D info. block innermost dimension (x_0) in order to
// avoid difficult checks at innermost loop
unsigned int j_indices[dim];
- FE_Q_Helper::zero_indices<dim> (j_indices);
+ internal::FE_Q_Base::zero_indices<dim> (j_indices);
for (unsigned int j=0; j<q_dofs_per_cell; j+=dofs1d)
{
unsigned int i_indices[dim];
- FE_Q_Helper::zero_indices<dim> (i_indices);
+ internal::FE_Q_Base::zero_indices<dim> (i_indices);
for (unsigned int i=0; i<q_dofs_per_cell; i+=dofs1d)
{
double val_extra_dim = 1.;
// update indices that denote the tensor product position. a bit
// fuzzy and therefore not done for innermost x_0 direction
- FE_Q_Helper::increment_indices<dim> (i_indices, dofs1d);
+ internal::FE_Q_Base::increment_indices<dim> (i_indices, dofs1d);
}
Assert (i_indices[dim-1] == 1, ExcInternalError());
- FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
+ internal::FE_Q_Base::increment_indices<dim> (j_indices, dofs1d);
}
// the discontinuous node is simply mapped on the discontinuous node on
this->poly_space.compute_value(index_map_inverse[j], point);
}
unsigned int j_indices[dim];
- FE_Q_Helper::zero_indices<dim> (j_indices);
+ internal::FE_Q_Base::zero_indices<dim> (j_indices);
double sum_check = 0;
for (unsigned int j = 0; j<q_dofs_per_cell; j += dofs1d)
{
my_restriction(mother_dof,child_dof)=val;
sum_check += val;
}
- FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
+ internal::FE_Q_Base::increment_indices<dim> (j_indices, dofs1d);
}
Assert (std::fabs(sum_check-1) <
std::max(eps, 5e-16*std::sqrt(this->dofs_per_cell)),
DEAL_II_NAMESPACE_OPEN
-namespace FE_Q_Bubbles_Helper
+namespace internal
{
- namespace
+ namespace FE_Q_Bubbles
{
- template <int dim, int spacedim>
- inline
- void
- compute_embedding_matrices(const FE_Q_Bubbles<dim, spacedim> &fe,
- std::vector<std::vector<FullMatrix<double> > > &matrices,
- const bool isotropic_only)
+ namespace
{
- const unsigned int dpc = fe.dofs_per_cell;
- const unsigned int degree = fe.degree;
+ template <int dim, int spacedim>
+ inline
+ void
+ compute_embedding_matrices(const dealii::FE_Q_Bubbles<dim, spacedim> &fe,
+ std::vector<std::vector<FullMatrix<double> > > &matrices,
+ const bool isotropic_only)
+ {
+ const unsigned int dpc = fe.dofs_per_cell;
+ const unsigned int degree = fe.degree;
- // Initialize quadrature formula on fine cells
- std::unique_ptr<Quadrature<dim> > q_fine;
- Quadrature<1> q_dummy(std::vector<Point<1> >(1), std::vector<double> (1,1.));
- switch (dim)
- {
- case 1:
- if (spacedim==1)
- q_fine.reset(new QGauss<dim> (degree+1));
- else if (spacedim==2)
- q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy));
- else
- q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy, q_dummy));
- break;
- case 2:
- if (spacedim==2)
+ // Initialize quadrature formula on fine cells
+ std::unique_ptr<Quadrature<dim> > q_fine;
+ Quadrature<1> q_dummy(std::vector<Point<1> >(1), std::vector<double> (1,1.));
+ switch (dim)
+ {
+ case 1:
+ if (spacedim==1)
+ q_fine.reset(new QGauss<dim> (degree+1));
+ else if (spacedim==2)
+ q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy));
+ else
+ q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy, q_dummy));
+ break;
+ case 2:
+ if (spacedim==2)
+ q_fine.reset(new QGauss<dim> (degree+1));
+ else
+ q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy));
+ break;
+ case 3:
q_fine.reset(new QGauss<dim> (degree+1));
- else
- q_fine.reset(new QAnisotropic<dim>(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy));
- break;
- case 3:
- q_fine.reset(new QGauss<dim> (degree+1));
- break;
- default:
- Assert(false, ExcInternalError());
- }
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
- Assert(q_fine.get() != nullptr, ExcInternalError());
- const unsigned int nq = q_fine->size();
+ Assert(q_fine.get() != nullptr, ExcInternalError());
+ const unsigned int nq = q_fine->size();
- // loop over all possible refinement cases
- unsigned int ref_case = (isotropic_only)
- ? RefinementCase<dim>::isotropic_refinement
- : RefinementCase<dim>::cut_x;
- for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
- {
- const unsigned int nc
- = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
-
- for (unsigned int i=0; i<nc; ++i)
- {
- Assert(matrices[ref_case-1][i].n() == dpc,
- ExcDimensionMismatch(matrices[ref_case-1][i].n(),dpc));
- Assert(matrices[ref_case-1][i].m() == dpc,
- ExcDimensionMismatch(matrices[ref_case-1][i].m(),dpc));
- }
-
- // create a respective refinement on the triangulation
- Triangulation<dim, spacedim> tr;
- GridGenerator::hyper_cube (tr, 0, 1);
- tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
- tr.execute_coarsening_and_refinement();
-
- DoFHandler<dim, spacedim> dh(tr);
- dh.distribute_dofs(fe);
-
- FEValues<dim, spacedim> fine (StaticMappingQ1<dim,spacedim>::mapping, fe, *q_fine,
- update_quadrature_points
- | update_JxW_values | update_values);
-
- const unsigned int n_dofs = dh.n_dofs();
-
- FullMatrix<double> fine_mass(n_dofs);
- FullMatrix<double> coarse_rhs_matrix(n_dofs, dpc);
-
- std::vector<std::vector<types::global_dof_index> > child_ldi
- (nc, std::vector<types::global_dof_index>(fe.dofs_per_cell));
-
- //now create the mass matrix and all the right_hand sides
- unsigned int child_no = 0;
- typename dealii::DoFHandler<dim>::active_cell_iterator cell
- = dh.begin_active();
- for (; cell!=dh.end(); ++cell, ++child_no)
- {
- fine.reinit(cell);
- cell->get_dof_indices(child_ldi[child_no]);
-
- for (unsigned int q=0; q<nq; ++q)
- for (unsigned int i=0; i<dpc; ++i)
- for (unsigned int j=0; j<dpc; ++j)
- {
- const unsigned int gdi=child_ldi[child_no][i];
- const unsigned int gdj=child_ldi[child_no][j];
- fine_mass(gdi, gdj)+=fine.shape_value(i,q)
- *fine.shape_value(j,q)
- *fine.JxW(q);
- Point<dim> quad_tmp;
- for (unsigned int k=0; k<dim; ++k)
- quad_tmp(k) = fine.quadrature_point(q)(k);
- coarse_rhs_matrix(gdi, j)
- +=fine.shape_value(i,q)
- *fe.shape_value(j, quad_tmp)
- *fine.JxW(q);
- }
- }
-
- //now solve for all right-hand sides simultaneously
- dealii::FullMatrix<double> solution (n_dofs, dpc);
- fine_mass.gauss_jordan();
- fine_mass.mmult(solution, coarse_rhs_matrix);
-
- //and distribute to the fine cell matrices
- for (unsigned int child_no=0; child_no<nc; ++child_no)
- for (unsigned int i=0; i<dpc; ++i)
- for (unsigned int j=0; j<dpc; ++j)
- {
- const unsigned int gdi=child_ldi[child_no][i];
- //remove small entries
- if (std::fabs(solution(gdi, j)) > 1.e-12)
- matrices[ref_case-1][child_no](i,j)=solution(gdi, j);
- }
- }
+ // loop over all possible refinement cases
+ unsigned int ref_case = (isotropic_only)
+ ? RefinementCase<dim>::isotropic_refinement
+ : RefinementCase<dim>::cut_x;
+ for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
+ {
+ const unsigned int nc
+ = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
+
+ for (unsigned int i=0; i<nc; ++i)
+ {
+ Assert(matrices[ref_case-1][i].n() == dpc,
+ ExcDimensionMismatch(matrices[ref_case-1][i].n(),dpc));
+ Assert(matrices[ref_case-1][i].m() == dpc,
+ ExcDimensionMismatch(matrices[ref_case-1][i].m(),dpc));
+ }
+
+ // create a respective refinement on the triangulation
+ dealii::Triangulation<dim, spacedim> tr;
+ GridGenerator::hyper_cube (tr, 0, 1);
+ tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
+ tr.execute_coarsening_and_refinement();
+
+ dealii::DoFHandler<dim, spacedim> dh(tr);
+ dh.distribute_dofs(fe);
+
+ dealii::FEValues<dim, spacedim> fine
+ (StaticMappingQ1<dim,spacedim>::mapping, fe, *q_fine,
+ update_quadrature_points
+ | update_JxW_values | update_values);
+
+ const unsigned int n_dofs = dh.n_dofs();
+
+ FullMatrix<double> fine_mass(n_dofs);
+ FullMatrix<double> coarse_rhs_matrix(n_dofs, dpc);
+
+ std::vector<std::vector<types::global_dof_index> > child_ldi
+ (nc, std::vector<types::global_dof_index>(fe.dofs_per_cell));
+
+ //now create the mass matrix and all the right_hand sides
+ unsigned int child_no = 0;
+ typename dealii::DoFHandler<dim>::active_cell_iterator cell
+ = dh.begin_active();
+ for (; cell!=dh.end(); ++cell, ++child_no)
+ {
+ fine.reinit(cell);
+ cell->get_dof_indices(child_ldi[child_no]);
+
+ for (unsigned int q=0; q<nq; ++q)
+ for (unsigned int i=0; i<dpc; ++i)
+ for (unsigned int j=0; j<dpc; ++j)
+ {
+ const unsigned int gdi=child_ldi[child_no][i];
+ const unsigned int gdj=child_ldi[child_no][j];
+ fine_mass(gdi, gdj)+=fine.shape_value(i,q)
+ *fine.shape_value(j,q)
+ *fine.JxW(q);
+ Point<dim> quad_tmp;
+ for (unsigned int k=0; k<dim; ++k)
+ quad_tmp(k) = fine.quadrature_point(q)(k);
+ coarse_rhs_matrix(gdi, j)
+ +=fine.shape_value(i,q)
+ *fe.shape_value(j, quad_tmp)
+ *fine.JxW(q);
+ }
+ }
+
+ //now solve for all right-hand sides simultaneously
+ dealii::FullMatrix<double> solution (n_dofs, dpc);
+ fine_mass.gauss_jordan();
+ fine_mass.mmult(solution, coarse_rhs_matrix);
+
+ //and distribute to the fine cell matrices
+ for (unsigned int child_no=0; child_no<nc; ++child_no)
+ for (unsigned int i=0; i<dpc; ++i)
+ for (unsigned int j=0; j<dpc; ++j)
+ {
+ const unsigned int gdi=child_ldi[child_no][i];
+ //remove small entries
+ if (std::fabs(solution(gdi, j)) > 1.e-12)
+ matrices[ref_case-1][child_no](i,j)=solution(gdi, j);
+ }
+ }
+ }
}
}
}
this->reinit_restriction_and_prolongation_matrices();
if (dim == spacedim)
{
- FE_Q_Bubbles_Helper::compute_embedding_matrices
+ internal::FE_Q_Bubbles::compute_embedding_matrices
(*this, this->prolongation, false);
// Fill restriction matrices with L2-projection
FETools::compute_projection_matrices (*this, this->restriction);
this->reinit_restriction_and_prolongation_matrices();
if (dim == spacedim)
{
- FE_Q_Bubbles_Helper::compute_embedding_matrices
+ internal::FE_Q_Bubbles::compute_embedding_matrices
(*this, this->prolongation, false);
// Fill restriction matrices with L2-projection
FETools::compute_projection_matrices (*this, this->restriction);
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
- /**
- * A function which maps in[i] to i,i.e. output[in[i]] = i;
- */
- inline
- std::vector<unsigned int>
- invert_numbering (const std::vector<unsigned int> &in)
+ namespace FE_Q_Hierarchical
{
- std::vector<unsigned int> out (in.size());
- for (unsigned int i=0; i<in.size(); ++i)
+ namespace
+ {
+ /**
+ * A function which maps in[i] to i,i.e. output[in[i]] = i;
+ */
+ inline
+ std::vector<unsigned int>
+ invert_numbering (const std::vector<unsigned int> &in)
{
- Assert (in[i] < out.size(),
- dealii::ExcIndexRange(in[i],0,out.size()));
- out[in[i]]=i;
+ std::vector<unsigned int> out (in.size());
+ for (unsigned int i=0; i<in.size(); ++i)
+ {
+ Assert (in[i] < out.size(),
+ dealii::ExcIndexRange(in[i],0,out.size()));
+ out[in[i]]=i;
+ }
+ return out;
}
- return out;
+ }
}
}
face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree)
{
FiniteElementData<dim-1> fe_data(FE_Q_Hierarchical<dim-1>::get_dpo_vector(degree),1,degree);
- return invert_numbering(FE_Q_Hierarchical<dim-1>::
- hierarchic_to_fe_q_hierarchical_numbering (fe_data));
+ return internal::FE_Q_Hierarchical::invert_numbering
+ (FE_Q_Hierarchical<dim-1>::hierarchic_to_fe_q_hierarchical_numbering (fe_data));
}
DEAL_II_NAMESPACE_OPEN
-namespace
+namespace internal
{
- bool IsNonZero (unsigned int i)
+ namespace FESystem
{
- return i>0;
- }
-
- unsigned int count_nonzeros(const std::vector<unsigned int> &vec)
- {
- return std::count_if(vec.begin(), vec.end(), IsNonZero);
+ namespace
+ {
+ unsigned int count_nonzeros(const std::vector<unsigned int> &vec)
+ {
+ return std::count_if(vec.begin(), vec.end(), [](const unsigned int i)
+ {
+ return i > 0;
+ });
+ }
+ }
}
-
}
/* ----------------------- FESystem::InternalData ------------------- */
FiniteElement<dim,spacedim> (FETools::Compositing::multiply_dof_numbers(fes, multiplicities),
FETools::Compositing::compute_restriction_is_additive_flags (fes, multiplicities),
FETools::Compositing::compute_nonzero_components(fes, multiplicities)),
- base_elements(count_nonzeros(multiplicities))
+ base_elements(internal::FESystem::count_nonzeros(multiplicities))
{
initialize(fes, multiplicities);
}
ExcDimensionMismatch (fes.size(), multiplicities.size()) );
Assert (fes.size() > 0,
ExcMessage ("Need to pass at least one finite element."));
- Assert (count_nonzeros(multiplicities) > 0,
+ Assert (internal::FESystem::count_nonzeros(multiplicities) > 0,
ExcMessage("You only passed FiniteElements with multiplicity 0."));
// Note that we need to skip every fe with multiplicity 0 in the following block of code
projects / dealii.git / commitdiff