Added support for FE_DGQArbitraryNodes in codimension one.

git-svn-id: https://svn.dealii.org/trunk@25604 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index f52970946e14e7f3da8d8bd958a553f81ae7e0f5..62b37298a71a3b0b9f6b9f312542f12cc9e480a7 100644 (file)
--- 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
 <a name="specific"></a>
 <h3>Specific improvements</h3>
 
+<li> New: The finite element type FE_DGQArbitraryNodes is now
+working also in codimension one spaces. 
+<br>
+(Luca Heltai, Andrea Mola 2012/06/06)
+
 <ol>
 <li> 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 921a7fc88d98f5e2be83a2bdaef12b83a80c3fc6..0737b718a4fc925a462201663b17edb52a220062 100644 (file)
--- 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<TensorProductPolynomials<dim>, dim, spacedim>
  *
  * @author F. Prill 2006
  */
-template <int dim>
-class FE_DGQArbitraryNodes : public FE_DGQ<dim>
+template <int dim,int spacedim=dim>
+class FE_DGQArbitraryNodes : public FE_DGQ<dim,spacedim>
 {
   public:
                                      /**
@@ -448,7 +448,7 @@ class FE_DGQArbitraryNodes : public FE_DGQ<dim>
                                       * This function is needed by the
                                       * constructors of @p FESystem.
                                       */
-    virtual FiniteElement<dim> *clone() const;
+    virtual FiniteElement<dim,spacedim> *clone() const;
 };
 
 

diff --git a/deal.II/source/fe/fe_dgq.cc b/deal.II/source/fe/fe_dgq.cc
index 6bf205aa83c2b7319c3776a48b6e32afed758be6..48c79431a311ce6fc58acf49373a0491d969457e 100644 (file)
--- a/deal.II/source/fe/fe_dgq.cc
+++ b/deal.II/source/fe/fe_dgq.cc
@@ -227,9 +227,9 @@ FE_DGQ<dim, spacedim>::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<dim, double, spacedim> (*this, this->prolongation);
                                    // Fill restriction matrices with L2-projection
-  FETools::compute_projection_matrices (*this, this->restriction);
+  FETools::compute_projection_matrices<dim, double, spacedim> (*this, this->restriction);
 
                                    // Compute support points, which
                                    // are the tensor product of the
@@ -677,16 +677,16 @@ FE_DGQ<dim, spacedim>::memory_consumption () const
 
 
 
-template <int dim>
-FE_DGQArbitraryNodes<dim>::FE_DGQArbitraryNodes (const Quadrature<1>& points)
-                : FE_DGQ<dim>(points)
+template <int dim, int spacedim>
+FE_DGQArbitraryNodes<dim,spacedim>::FE_DGQArbitraryNodes (const Quadrature<1>& points)
+               : FE_DGQ<dim,spacedim>(points)
 {}
 
 
 
-template <int dim>
+template <int dim, int spacedim>
 std::string
-FE_DGQArbitraryNodes<dim>::get_name () const
+FE_DGQArbitraryNodes<dim,spacedim>::get_name () const
 {
                                    // note that the
                                    // FETools::get_fe_from_name
@@ -751,9 +751,9 @@ FE_DGQArbitraryNodes<dim>::get_name () const
 
 
 
-template <int dim>
-FiniteElement<dim> *
-FE_DGQArbitraryNodes<dim>::clone() const
+template <int dim, int spacedim>
+FiniteElement<dim,spacedim> *
+FE_DGQArbitraryNodes<dim,spacedim>::clone() const
 {
                                    // TODO[Prill] : There must be a better way
                                    // to extract 1D quadrature points from the
@@ -766,7 +766,7 @@ FE_DGQArbitraryNodes<dim>::clone() const
     qpoints[i] = Point<1>(this->unit_support_points[i][0]);
   Quadrature<1> pquadrature(qpoints);
 
-  return new FE_DGQArbitraryNodes<dim>(pquadrature);
+  return new FE_DGQArbitraryNodes<dim,spacedim>(pquadrature);
 }
 
 

diff --git a/deal.II/source/fe/fe_dgq.inst.in b/deal.II/source/fe/fe_dgq.inst.in
index 65b35cc07cd4ecf1bb5c6f1772d35da17ed216db..ebef67e7483257731912981c538fd8f3c6cbad6a 100644 (file)
--- 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<deal_II_dimension>;
+
+#if deal_II_dimension != 3
+    template class FE_DGQArbitraryNodes<deal_II_dimension, deal_II_dimension+1>;
+#endif
+
   }
 

diff --git a/deal.II/source/fe/fe_tools.cc b/deal.II/source/fe/fe_tools.cc
index facc4b8efd0b87470d8b60d9f1162a73c664fb4b..679aa7edf1aadeab05126a3860a2b19dbd53b71c 100644 (file)
--- 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<int dim, typename number, int spacedim>
@@ -688,7 +688,7 @@ namespace FETools
     compute_embedding_for_shape_function (
       const unsigned int i,
       const FiniteElement<dim, spacedim>& fe,
-      const FEValues<dim>& coarse,
+      const FEValues<dim, spacedim>& coarse,
       const Householder<double>& H,
       FullMatrix<number>& this_matrix)
     {
@@ -758,15 +758,15 @@ namespace FETools
       tria.begin_active()->set_refine_flag (RefinementCase<dim>(ref_case));
       tria.execute_coarsening_and_refinement ();
 
-      MappingCartesian<dim> mapping;
+      MappingQ1<dim,spacedim> mapping;
       const unsigned int degree = fe.degree;
       QGauss<dim> q_fine (degree+1);
       const unsigned int nq = q_fine.size();
 
-      FEValues<dim> fine (mapping, fe, q_fine,
-                          update_quadrature_points |
-                          update_JxW_values |
-                          update_values);
+      FEValues<dim,spacedim> 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<void> task_group;
 
-      for (typename Triangulation<dim>::active_cell_iterator
+      for (typename Triangulation<dim,spacedim>::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<dim> q_coarse (fine.get_quadrature_points (),
-                                          fine.get_JxW_values ());
-          FEValues<dim> 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<Point<spacedim> > &q_points_fine = fine.get_quadrature_points();
+          std::vector<Point<dim> > q_points_coarse(q_points_fine.size());
+          for (unsigned int i=0;i<q_points_fine.size();++i)
+              for (unsigned int j=0;j<dim;++j)
+                  q_points_coarse[i](j) = q_points_fine[i](j);
+         const Quadrature<dim> q_coarse (q_points_coarse,
+                                         fine.get_JxW_values ());
+         FEValues<dim,spacedim> 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 <int dim, typename number, int spacedim>
@@ -1112,7 +1117,7 @@ namespace FETools
 
                                      // prepare FEValues, quadrature etc on
                                      // coarse cell
-    MappingCartesian<dim> mapping;
+    MappingQ1<dim,spacedim> mapping;
     QGauss<dim> q_fine(degree+1);
     const unsigned int nq = q_fine.size();
 
@@ -1123,7 +1128,7 @@ namespace FETools
       Triangulation<dim,spacedim> tr;
       GridGenerator::hyper_cube (tr, 0, 1);
 
-      FEValues<dim> coarse (mapping, fe, q_fine,
+      FEValues<dim,spacedim> coarse (mapping, fe, q_fine,
                             update_JxW_values | update_values);
 
       typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
@@ -1181,9 +1186,9 @@ namespace FETools
         tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
         tr.execute_coarsening_and_refinement();
 
-        FEValues<dim> fine (mapping, fe, q_fine,
-                            update_quadrature_points | update_JxW_values |
-                            update_values);
+       FEValues<dim,spacedim> fine (mapping, fe, q_fine,
+                                    update_quadrature_points | update_JxW_values |
+                                    update_values);
 
         typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
           = tr.begin(0);
@@ -1195,15 +1200,20 @@ namespace FETools
           {
             FullMatrix<double> &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<dim> q_coarse (fine.get_quadrature_points(),
-                                      fine.get_JxW_values());
-            FEValues<dim> 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<Point<spacedim> > &q_points_fine = fine.get_quadrature_points();
+            std::vector<Point<dim> > q_points_coarse(q_points_fine.size());
+            for (unsigned int q=0;q<q_points_fine.size();++q)
+                for (unsigned int j=0;j<dim;++j)
+                    q_points_coarse[q](j) = q_points_fine[q](j);            
+           Quadrature<dim> q_coarse (q_points_coarse,
+                                     fine.get_JxW_values());
+           FEValues<dim,spacedim> coarse (mapping, fe, q_coarse, update_values);
+           coarse.reinit(coarse_cell);
 
                                              // Build RHS