From b63094d3cfed6720e5f4fc58891cdaaaa5a701fc Mon Sep 17 00:00:00 2001 From: heltai Date: Wed, 6 Jun 2012 12:36:34 +0000 Subject: [PATCH] Added support for FE_DGQArbitraryNodes in codimension one. git-svn-id: https://svn.dealii.org/trunk@25604 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 5 ++ deal.II/include/deal.II/fe/fe_dgq.h | 6 +-- deal.II/source/fe/fe_dgq.cc | 22 ++++----- deal.II/source/fe/fe_dgq.inst.in | 5 ++ deal.II/source/fe/fe_tools.cc | 74 ++++++++++++++++------------- 5 files changed, 66 insertions(+), 46 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index f52970946e..62b37298a7 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -207,6 +207,11 @@ enabled due to a missing include file in file

Specific improvements

+
  • New: The finite element type FE_DGQArbitraryNodes is now +working also in codimension one spaces. +
    +(Luca Heltai, Andrea Mola 2012/06/06) +
    1. Fixed: Computing the $W^{1,\infty}$ norm and seminorm in VectorTools::integrate_difference was not implemented. This is now diff --git a/deal.II/include/deal.II/fe/fe_dgq.h b/deal.II/include/deal.II/fe/fe_dgq.h index 921a7fc88d..0737b718a4 100644 --- a/deal.II/include/deal.II/fe/fe_dgq.h +++ b/deal.II/include/deal.II/fe/fe_dgq.h @@ -413,8 +413,8 @@ class FE_DGQ : public FE_Poly, dim, spacedim> * * @author F. Prill 2006 */ -template -class FE_DGQArbitraryNodes : public FE_DGQ +template +class FE_DGQArbitraryNodes : public FE_DGQ { public: /** @@ -448,7 +448,7 @@ class FE_DGQArbitraryNodes : public FE_DGQ * This function is needed by the * constructors of @p FESystem. */ - virtual FiniteElement *clone() const; + virtual FiniteElement *clone() const; }; diff --git a/deal.II/source/fe/fe_dgq.cc b/deal.II/source/fe/fe_dgq.cc index 6bf205aa83..48c79431a3 100644 --- a/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/source/fe/fe_dgq.cc @@ -227,9 +227,9 @@ FE_DGQ::FE_DGQ (const Quadrature<1>& points) // matrices to the right sizes this->reinit_restriction_and_prolongation_matrices(); // Fill prolongation matrices with embedding operators - FETools::compute_embedding_matrices (*this, this->prolongation); + FETools::compute_embedding_matrices (*this, this->prolongation); // Fill restriction matrices with L2-projection - FETools::compute_projection_matrices (*this, this->restriction); + FETools::compute_projection_matrices (*this, this->restriction); // Compute support points, which // are the tensor product of the @@ -677,16 +677,16 @@ FE_DGQ::memory_consumption () const -template -FE_DGQArbitraryNodes::FE_DGQArbitraryNodes (const Quadrature<1>& points) - : FE_DGQ(points) +template +FE_DGQArbitraryNodes::FE_DGQArbitraryNodes (const Quadrature<1>& points) + : FE_DGQ(points) {} -template +template std::string -FE_DGQArbitraryNodes::get_name () const +FE_DGQArbitraryNodes::get_name () const { // note that the // FETools::get_fe_from_name @@ -751,9 +751,9 @@ FE_DGQArbitraryNodes::get_name () const -template -FiniteElement * -FE_DGQArbitraryNodes::clone() const +template +FiniteElement * +FE_DGQArbitraryNodes::clone() const { // TODO[Prill] : There must be a better way // to extract 1D quadrature points from the @@ -766,7 +766,7 @@ FE_DGQArbitraryNodes::clone() const qpoints[i] = Point<1>(this->unit_support_points[i][0]); Quadrature<1> pquadrature(qpoints); - return new FE_DGQArbitraryNodes(pquadrature); + return new FE_DGQArbitraryNodes(pquadrature); } diff --git a/deal.II/source/fe/fe_dgq.inst.in b/deal.II/source/fe/fe_dgq.inst.in index 65b35cc07c..ebef67e748 100644 --- a/deal.II/source/fe/fe_dgq.inst.in +++ b/deal.II/source/fe/fe_dgq.inst.in @@ -23,5 +23,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS for (deal_II_dimension : DIMENSIONS) { template class FE_DGQArbitraryNodes; + +#if deal_II_dimension != 3 + template class FE_DGQArbitraryNodes; +#endif + } diff --git a/deal.II/source/fe/fe_tools.cc b/deal.II/source/fe/fe_tools.cc index facc4b8efd..679aa7edf1 100644 --- a/deal.II/source/fe/fe_tools.cc +++ b/deal.II/source/fe/fe_tools.cc @@ -649,7 +649,7 @@ namespace FETools } - +/* template<> void compute_embedding_matrices(const FiniteElement<1,2> &, @@ -680,7 +680,7 @@ namespace FETools Assert(false, ExcNotImplemented()); } - +*/ namespace { template @@ -688,7 +688,7 @@ namespace FETools compute_embedding_for_shape_function ( const unsigned int i, const FiniteElement& fe, - const FEValues& coarse, + const FEValues& coarse, const Householder& H, FullMatrix& this_matrix) { @@ -758,15 +758,15 @@ namespace FETools tria.begin_active()->set_refine_flag (RefinementCase(ref_case)); tria.execute_coarsening_and_refinement (); - MappingCartesian mapping; + MappingQ1 mapping; const unsigned int degree = fe.degree; QGauss q_fine (degree+1); const unsigned int nq = q_fine.size(); - FEValues fine (mapping, fe, q_fine, - update_quadrature_points | - update_JxW_values | - update_values); + FEValues fine (mapping, fe, q_fine, + update_quadrature_points | + update_JxW_values | + update_values); // We search for the polynomial on // the small cell, being equal to @@ -795,19 +795,24 @@ namespace FETools Threads::TaskGroup task_group; - for (typename Triangulation::active_cell_iterator + for (typename Triangulation::active_cell_iterator fine_cell = tria.begin_active (); fine_cell != tria.end (); ++fine_cell, ++cell_number) { fine.reinit (fine_cell); - // evaluate on the coarse cell (which - // is the first -- inactive -- cell on - // the lowest level of the - // triangulation we have created) - const Quadrature q_coarse (fine.get_quadrature_points (), - fine.get_JxW_values ()); - FEValues coarse (mapping, fe, q_coarse, update_values); + // evaluate on the coarse cell (which + // is the first -- inactive -- cell on + // the lowest level of the + // triangulation we have created) + const std::vector > &q_points_fine = fine.get_quadrature_points(); + std::vector > q_points_coarse(q_points_fine.size()); + for (unsigned int i=0;i q_coarse (q_points_coarse, + fine.get_JxW_values ()); + FEValues coarse (mapping, fe, q_coarse, update_values); coarse.reinit (tria.begin (0)); @@ -1072,7 +1077,7 @@ namespace FETools } - +/* template <> void compute_projection_matrices(const FiniteElement<1,2>&, @@ -1097,7 +1102,7 @@ namespace FETools { Assert(false, ExcNotImplemented()); } - +*/ template @@ -1112,7 +1117,7 @@ namespace FETools // prepare FEValues, quadrature etc on // coarse cell - MappingCartesian mapping; + MappingQ1 mapping; QGauss q_fine(degree+1); const unsigned int nq = q_fine.size(); @@ -1123,7 +1128,7 @@ namespace FETools Triangulation tr; GridGenerator::hyper_cube (tr, 0, 1); - FEValues coarse (mapping, fe, q_fine, + FEValues coarse (mapping, fe, q_fine, update_JxW_values | update_values); typename Triangulation::cell_iterator coarse_cell @@ -1181,9 +1186,9 @@ namespace FETools tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); tr.execute_coarsening_and_refinement(); - FEValues fine (mapping, fe, q_fine, - update_quadrature_points | update_JxW_values | - update_values); + FEValues fine (mapping, fe, q_fine, + update_quadrature_points | update_JxW_values | + update_values); typename Triangulation::cell_iterator coarse_cell = tr.begin(0); @@ -1195,15 +1200,20 @@ namespace FETools { FullMatrix &this_matrix = matrices[ref_case-1][cell_number]; - // Compute right hand side, - // which is a fine level basis - // function tested with the - // coarse level functions. - fine.reinit(coarse_cell->child(cell_number)); - Quadrature q_coarse (fine.get_quadrature_points(), - fine.get_JxW_values()); - FEValues coarse (mapping, fe, q_coarse, update_values); - coarse.reinit(coarse_cell); + // Compute right hand side, + // which is a fine level basis + // function tested with the + // coarse level functions. + fine.reinit(coarse_cell->child(cell_number)); + const std::vector > &q_points_fine = fine.get_quadrature_points(); + std::vector > q_points_coarse(q_points_fine.size()); + for (unsigned int q=0;q q_coarse (q_points_coarse, + fine.get_JxW_values()); + FEValues coarse (mapping, fe, q_coarse, update_values); + coarse.reinit(coarse_cell); // Build RHS -- 2.39.5