Implementat HighOrderGrid with all existing test

Had to slightly change the Gaussian bump test to ensure convergence by increasing the initial grid size.
As a result, we do one less grid level to keep it tractable on a personal computer.

Also fixed a bug regarding the DoFCellAccessor being used with the Lagrange FiniteElement.

src/dg/dg.h

@@ -25,6 +25,7 @@
 
 #include <Sacado.hpp>
 
+#include "dg/high_order_grid.h"
 #include "physics/physics.h"
 #include "numerical_flux/numerical_flux.h"
 #include "parameters/all_parameters.h"
@@ -100,9 +101,6 @@ class DGBase
     using MassiveCollectionTuple = std::tuple<
         //dealii::hp::MappingCollection<dim>, // Mapping
         dealii::hp::FECollection<dim>, // Solution FE
-        //dealii::hp::FECollection<dim>, // Grid FE
-        //dealii::FESystem<dim>, // Grid FE // Can't just return FESystem. For some reason the copy constructor is deleted
-        dealii::FE_Q<dim>, // Grid FE
         dealii::hp::QCollection<dim>,  // Volume quadrature
         dealii::hp::QCollection<dim-1>, // Face quadrature
         dealii::hp::QCollection<1>, // 1D quadrature for strong form
@@ -129,7 +127,7 @@ class DGBase
     /** Since we are interested in performing mesh movement for optimization purposes,
      *  this is not a constant member variables.
      */
-    dealii::hp::MappingCollection<dim> mapping_collection;
+    //dealii::hp::MappingCollection<dim> mapping_collection;
 
     void set_all_cells_fe_degree ( const unsigned int degree );
 
@@ -219,9 +217,6 @@ class DGBase
      */
     dealii::LinearAlgebra::distributed::Vector<double> solution;
 
-    /// Current nodal coefficients of the high-order grid.
-    dealii::LinearAlgebra::distributed::Vector<double> high_order_grid;
-
     void initialize_manufactured_solution (); ///< Virtual function defined in DG
 
     void output_results_vtk (const unsigned int ith_grid); ///< Output solution
@@ -271,13 +266,30 @@ class DGBase
      *  Therefore, every grid/cell will use the maximal polynomial mapping regardless of the solution order.
      */
     //const dealii::hp::FECollection<dim>    fe_collection_grid;
-    const dealii::FESystem<dim>    fe_grid;
+    //const dealii::FESystem<dim>    fe_grid;
 
     dealii::hp::QCollection<dim>     volume_quadrature_collection;
     dealii::hp::QCollection<dim-1>   face_quadrature_collection;
     dealii::hp::QCollection<1>       oned_quadrature_collection;
 
+protected:
+    /// Lagrange basis used in strong form
+    /** This is a collection of scalar Lagrange bases */
+    const dealii::hp::FECollection<dim>  fe_collection_lagrange;
 
+public:
+    /// Degrees of freedom handler
+    /*  Allows us to iterate over the finite elements' degrees of freedom.
+     *  Note that since we are not using FESystem, we need to multiply
+     *  the index by a factor of "nstate"
+     *
+     *  Must be defined after fe_dg since it is a subscriptor of fe_dg.
+     *  Destructor are called in reverse order in which they appear in class definition. 
+     */ 
+    dealii::hp::DoFHandler<dim> dof_handler;
+
+    /// High order grid that will provide the MappingFEField
+    HighOrderGrid<dim,real> high_order_grid;
 
 protected:
 
@@ -370,31 +382,8 @@ class DGBase
     const dealii::UpdateFlags neighbor_face_update_flags = dealii::update_values | dealii::update_gradients;
 
 
-    /// Lagrange basis used in strong form
-    /** This is a collection of scalar Lagrange bases */
-    const dealii::hp::FECollection<dim>  fe_collection_lagrange;
-    //dealii::hp::FEValues<dim,dim>        fe_values_collection_volume_lagrange;
-
-
-public:
-    /// Degrees of freedom handler
-    /*  Allows us to iterate over the finite elements' degrees of freedom.
-     *  Note that since we are not using FESystem, we need to multiply
-     *  the index by a factor of "nstate"
-     *
-     *  Must be defined after fe_dg since it is a subscriptor of fe_dg.
-     *  Destructor are called in reverse order in which they appear in class definition. 
-     */ 
-    dealii::hp::DoFHandler<dim> dof_handler;
-
-    /// Degrees of freedom handler for the high-order grid
-    dealii::DoFHandler<dim> dof_handler_grid;
 
 protected:
-    /// Main underlying mapping responsible for all mapping evaluations
-    //dealii::MappingFEField<dim,dim,dealii::LinearAlgebra::distributed::Vector<double>, dealii::hp::DoFHandler<dim> > high_order_grid_mapping;
-    dealii::MappingFEField<dim,dim,dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<dim> > high_order_grid_mapping;
-
     MPI_Comm mpi_communicator; ///< MPI communicator
     dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0
 private:

src/dg/high_order_grid.cpp

@@ -23,6 +23,10 @@ HighOrderGrid<dim,real,VectorType,DoFHandlerType>::HighOrderGrid(
     , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0)
 {
     allocate();
+    const dealii::ComponentMask mask(dim, true);
+    dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask);
+    nodes.update_ghost_values();
+    mapping_fe_field = std::make_shared< dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> > (dof_handler_grid,nodes,mask);
 }
 
 template <int dim, typename real, typename VectorType , typename DoFHandlerType>
@@ -32,41 +36,59 @@ HighOrderGrid<dim,real,VectorType,DoFHandlerType>::allocate()
     dof_handler_grid.initialize(*triangulation, fe_system);
     dof_handler_grid.distribute_dofs(fe_system);
 
+#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
+    nodes.reinit(dof_handler_grid.n_dofs());
+#else
     locally_owned_dofs_grid = dof_handler_grid.locally_owned_dofs();
     dealii::DoFTools::extract_locally_relevant_dofs(dof_handler_grid, ghost_dofs_grid);
     locally_relevant_dofs_grid = ghost_dofs_grid;
     ghost_dofs_grid.subtract_set(locally_owned_dofs_grid);
     nodes.reinit(locally_owned_dofs_grid, ghost_dofs_grid, mpi_communicator);
+#endif
 }
 
-template <int dim, typename real, typename VectorType , typename DoFHandlerType>
-dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> 
-HighOrderGrid<dim,real,VectorType,DoFHandlerType>::get_MappingFEField() {
-    const dealii::ComponentMask mask(dim, true);
-    dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask);
-
-    dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> mapping(dof_handler_grid,nodes,mask);
-
-    return mapping;
-}
+//template <int dim, typename real, typename VectorType , typename DoFHandlerType>
+//dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> 
+//HighOrderGrid<dim,real,VectorType,DoFHandlerType>::get_MappingFEField() {
+//    const dealii::ComponentMask mask(dim, true);
+//    dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask);
+//
+//    dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> mapping(dof_handler_grid,nodes,mask);
+//
+//    return mapping;
+//}
+//
 
 
 template <int dim, typename real, typename VectorType , typename DoFHandlerType>
 void HighOrderGrid<dim,real,VectorType,DoFHandlerType>::prepare_for_coarsening_and_refinement() {
 
     old_nodes = nodes;
     old_nodes.update_ghost_values();
+#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
+    solution_transfer.prepare_for_coarsening_and_refinement(old_nodes);
+#else
     solution_transfer.prepare_for_coarsening_and_refinement(old_nodes);
+#endif
 }
 
 template <int dim, typename real, typename VectorType , typename DoFHandlerType>
 void HighOrderGrid<dim,real,VectorType,DoFHandlerType>::execute_coarsening_and_refinement() {
     allocate();
+#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
+    solution_transfer.interpolate(old_nodes, nodes);
+#else
     solution_transfer.interpolate(nodes);
+#endif
     nodes.update_ghost_values();
+    mapping_fe_field = std::make_shared< dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> > (dof_handler_grid,nodes);
 }
 
 
 //template class HighOrderGrid<PHILIP_DIM, double>;
+#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
+template class HighOrderGrid<PHILIP_DIM, double, dealii::Vector<double>, dealii::DoFHandler<PHILIP_DIM>>;
+#else
 template class HighOrderGrid<PHILIP_DIM, double, dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<PHILIP_DIM>>;
+#endif
 } // namespace PHiLiP

src/dg/high_order_grid.h

@@ -8,6 +8,7 @@
 
 #include <deal.II/dofs/dof_handler.h>
 
+#include <deal.II/numerics/solution_transfer.h>
 #include <deal.II/distributed/solution_transfer.h>
 
 #include <deal.II/lac/vector.h>
@@ -22,10 +23,24 @@ namespace PHiLiP {
  *  nodes vector. The dof_handler_grid is used to access and loop through those nodes.
  *  This will especially be useful when performing shape optimization, where the surface and volume 
  *  nodes will need to be displaced.
+ *  Note that there are a lot of pre-processor statements, and that is because the SolutionTransfer class
+ *  and the Vector class act quite differently between serial and parallel implementation. Hopefully,
+ *  deal.II will change this one day such that we have one interface for both.
  */
+#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
+template <int dim = PHILIP_DIM, typename real = double, typename VectorType = dealii::Vector<double>, typename DoFHandlerType = dealii::DoFHandler<PHILIP_DIM>>
+#else
 template <int dim = PHILIP_DIM, typename real = double, typename VectorType = dealii::LinearAlgebra::distributed::Vector<double>, typename DoFHandlerType = dealii::DoFHandler<PHILIP_DIM>>
+#endif
 class HighOrderGrid
 {
+#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
+    using Vector = dealii::Vector<double>;
+    using SolutionTransfer = dealii::SolutionTransfer<dim, dealii::Vector<double>, dealii::DoFHandler<dim>>;
+#else
+    using Vector = dealii::LinearAlgebra::distributed::Vector<double>;
+    using SolutionTransfer = dealii::parallel::distributed::SolutionTransfer<dim, dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<dim>>;
+#endif
 public:
     /// Principal constructor that will call delegated constructor.
     HighOrderGrid(const Parameters::AllParameters *const parameters_input, const unsigned int max_degree, dealii::Triangulation<dim> *const triangulation_input);
@@ -36,6 +51,9 @@ class HighOrderGrid
     /// Return a MappingFEField that corresponds to the current node locations
     dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType> get_MappingFEField();
 
+    /// Return a MappingFEField that corresponds to the current node locations
+    void update_MappingFEField();
+
     const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters
 
     /// Maximum degree of the geometry polynomial representing the grid.
@@ -47,13 +65,13 @@ class HighOrderGrid
     dealii::DoFHandler<dim> dof_handler_grid;
 
     /// Current nodal coefficients of the high-order grid.
-    dealii::LinearAlgebra::distributed::Vector<double> nodes;
+    Vector nodes;
 
     /// List of surface points
-    dealii::LinearAlgebra::distributed::Vector<double> surface_nodes;
+    Vector surface_nodes;
 
     /// RBF mesh deformation  -  To be done
-    void deform_mesh(dealii::LinearAlgebra::distributed::Vector<double> surface_displacements);
+    void deform_mesh(Vector surface_displacements);
 
     /// Evaluate cell metric Jacobian
     /** The metric Jacobian is given by the gradient of the physical location
@@ -84,18 +102,26 @@ class HighOrderGrid
     /// Using system of polynomials to represent the x, y, and z directions.
     const dealii::FESystem<dim> fe_system;
 
+
+    /** MappingFEField that will provide the polynomial-based grid.
+     *  It is a shared smart pointer because the constructor requires the dof_handler_grid to be properly initialized.
+     *  See discussion in the following thread:
+     *  https://stackoverflow.com/questions/7557153/defining-an-object-without-calling-its-constructor-in-c
+     */
+    std::shared_ptr<dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType>> mapping_fe_field;
+
     dealii::IndexSet locally_owned_dofs_grid; ///< Locally own degrees of freedom for the grid
     dealii::IndexSet ghost_dofs_grid; ///< Locally relevant ghost degrees of freedom for the grid
     dealii::IndexSet locally_relevant_dofs_grid; ///< Union of locally owned degrees of freedom and relevant ghost degrees of freedom for the grid
 protected:
 
     /// Used for the SolutionTransfer when performing grid adaptation.
-    dealii::LinearAlgebra::distributed::Vector<double> old_nodes;
+    Vector old_nodes;
 
     /** Transfers the coarse curved curve onto the fine curved grid.
      *  Used in prepare_for_coarsening_and_refinement() and execute_coarsening_and_refinement()
      */
-    dealii::parallel::distributed::SolutionTransfer<dim, dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<dim>> solution_transfer;
+    SolutionTransfer solution_transfer;
 
     MPI_Comm mpi_communicator; ///< MPI communicator
     dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0

src/dg/strong_dg.cpp

@@ -581,8 +581,6 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
 
             // Convective
             // Now minus such 2 integrations by parts
-            assert(JxW[iquad] - fe_values_lagrange.JxW(iquad) < 1e-14);
-
             rhs = rhs - fe_values_vol.shape_value_component(itest,iquad,istate) * flux_divergence[iquad][istate] * JxW[iquad];
 
             //// Diffusive

src/testing/euler_cylinder.cpp

@@ -206,7 +206,9 @@ ::run_test () const
                 dealii::LinearAlgebra::distributed::Vector<double> old_solution(dg->solution);
                 dealii::parallel::distributed::SolutionTransfer<dim, dealii::LinearAlgebra::distributed::Vector<double>, dealii::hp::DoFHandler<dim>> solution_transfer(dg->dof_handler);
                 solution_transfer.prepare_for_coarsening_and_refinement(old_solution);
+                dg->high_order_grid.prepare_for_coarsening_and_refinement();
                 grid.refine_global (1);
+                dg->high_order_grid.execute_coarsening_and_refinement();
                 dg->allocate_system ();
                 dg->solution.zero_out_ghosts();
                 solution_transfer.interpolate(dg->solution);
@@ -236,7 +238,7 @@ ::run_test () const
             // Overintegrate the error to make sure there is not integration error in the error estimate
             int overintegrate = 10;
             dealii::QGauss<dim> quad_extra(dg->max_degree+1+overintegrate);
-            dealii::FEValues<dim,dim> fe_values_extra(dg->mapping_collection[poly_degree], dg->fe_collection[poly_degree], quad_extra, 
+            dealii::FEValues<dim,dim> fe_values_extra(*(dg->high_order_grid.mapping_fe_field), dg->fe_collection[poly_degree], quad_extra, 
                     dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
             const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points;
             std::array<double,nstate> soln_at_q;

src/testing/euler_gaussian_bump.cpp

@@ -191,7 +191,7 @@ ::run_test () const
         //n_subdivisions[1] = n_1d_cells[0]; // y-direction
         //n_subdivisions[0] = 4*n_subdivisions[1]; // x-direction
         n_subdivisions[1] = n_1d_cells[0]; // y-direction
-        n_subdivisions[0] = 3*n_subdivisions[1]; // x-direction
+        n_subdivisions[0] = 9*n_subdivisions[1]; // x-direction
         dealii::Point<2> p1(-1.5,0.0), p2(1.5,y_height);
         const bool colorize = true;
         dealii::parallel::distributed::Triangulation<dim> grid(this->mpi_communicator,
@@ -224,7 +224,6 @@ ::run_test () const
         // Create DG object
         std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, &grid);
 
-
         // Initialize coarse grid solution with free-stream
         dg->allocate_system ();
         dealii::VectorTools::interpolate(dg->dof_handler, initial_conditions, dg->solution);
@@ -240,7 +239,9 @@ ::run_test () const
                 dealii::LinearAlgebra::distributed::Vector<double> old_solution(dg->solution);
                 dealii::parallel::distributed::SolutionTransfer<dim, dealii::LinearAlgebra::distributed::Vector<double>, dealii::hp::DoFHandler<dim>> solution_transfer(dg->dof_handler);
                 solution_transfer.prepare_for_coarsening_and_refinement(old_solution);
+                dg->high_order_grid.prepare_for_coarsening_and_refinement();
                 grid.refine_global (1);
+                dg->high_order_grid.execute_coarsening_and_refinement();
                 dg->allocate_system ();
                 dg->solution.zero_out_ghosts();
                 solution_transfer.interpolate(dg->solution);
@@ -346,12 +347,13 @@ ::run_test () const
 
         const double last_slope = log(entropy_error[n_grids-1]/entropy_error[n_grids-2])
                                   / log(grid_size[n_grids-1]/grid_size[n_grids-2]);
-        double before_last_slope = last_slope;
-        if ( n_grids > 2 ) {
-        before_last_slope = log(entropy_error[n_grids-2]/entropy_error[n_grids-3])
-                            / log(grid_size[n_grids-2]/grid_size[n_grids-3]);
-        }
-        const double slope_avg = 0.5*(before_last_slope+last_slope);
+        //double before_last_slope = last_slope;
+        //if ( n_grids > 2 ) {
+        //    before_last_slope = log(entropy_error[n_grids-2]/entropy_error[n_grids-3])
+        //                        / log(grid_size[n_grids-2]/grid_size[n_grids-3]);
+        //}
+        //const double slope_avg = 0.5*(before_last_slope+last_slope);
+        const double slope_avg = last_slope;
         const double slope_diff = slope_avg-expected_slope;
 
         double slope_deficit_tolerance = -std::abs(manu_grid_conv_param.slope_deficit_tolerance);

src/testing/grid_study.cpp

@@ -63,7 +63,7 @@ ::integrate_solution_over_domain(DGBase<dim,double> &dg) const
     int overintegrate = 10;
     dealii::QGauss<dim> quad_extra(dg.max_degree+1+overintegrate);
     //dealii::MappingQ<dim,dim> mappingq_temp(dg.max_degree+1);
-    dealii::FEValues<dim,dim> fe_values_extra(dg.mapping_collection[dg.max_degree], dg.fe_collection[dg.max_degree], quad_extra, 
+    dealii::FEValues<dim,dim> fe_values_extra(*(dg.high_order_grid.mapping_fe_field), dg.fe_collection[dg.max_degree], quad_extra, 
             dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
     const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points;
     std::array<double,nstate> soln_at_q;
@@ -260,7 +260,7 @@ ::run_test () const
             int overintegrate = 10;
             dealii::QGauss<dim> quad_extra(dg->max_degree+1+overintegrate);
             //dealii::MappingQ<dim,dim> mappingq(dg->max_degree+1);
-            dealii::FEValues<dim,dim> fe_values_extra(dg->mapping_collection[poly_degree], dg->fe_collection[poly_degree], quad_extra, 
+            dealii::FEValues<dim,dim> fe_values_extra(*(dg->high_order_grid.mapping_fe_field), dg->fe_collection[poly_degree], quad_extra, 
                     dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
             const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points;
             std::array<double,nstate> soln_at_q;

tests/integration_tests_control_files/euler_integration/2d_euler_gaussian_bump.prm

@@ -60,9 +60,9 @@ subsection manufactured solution convergence study
   set grid_progression_add  = 0
 
   # Initial grid of size (initial_grid_size)^dim
-  set initial_grid_size = 3
+  set initial_grid_size = 7
 
   # Number of grids in grid study
-  set number_of_grids   = 4
+  set number_of_grids   = 3
 end