Implement interior penalty discontinuous Galerkin method for temperat…

…ure and composition fields. Includes 2 tests.

Run astyle

include/aspect/assembly.h

@@ -117,11 +117,16 @@ namespace aspect
                            const FiniteElement<dim> &advection_element,
                            const Mapping<dim>       &mapping,
                            const Quadrature<dim>    &quadrature,
+                           const Quadrature<dim-1>  &face_quadrature,
                            const unsigned int        n_compositional_fields);
           AdvectionSystem (const AdvectionSystem &data);
 
           FEValues<dim>               finite_element_values;
 
+          std::shared_ptr<FEFaceValues<dim> >           face_finite_element_values;
+          std::shared_ptr<FEFaceValues<dim> >           neighbor_face_finite_element_values;
+          std::shared_ptr<FESubfaceValues<dim> >        subface_finite_element_values;
+
           std::vector<types::global_dof_index>   local_dof_indices;
 
           /**
@@ -134,6 +139,10 @@ namespace aspect
            */
           std::vector<double>         phi_field;
           std::vector<Tensor<1,dim> > grad_phi_field;
+          std::vector<double>         face_phi_field;
+          std::vector<Tensor<1,dim> > face_grad_phi_field;
+          std::vector<double>         neighbor_face_phi_field;
+          std::vector<Tensor<1,dim> > neighbor_face_grad_phi_field;
 
           std::vector<Tensor<1,dim> > old_velocity_values;
           std::vector<Tensor<1,dim> > old_old_velocity_values;
@@ -161,7 +170,9 @@ namespace aspect
 
           std::vector<double>         current_temperature_values;
           std::vector<Tensor<1,dim> > current_velocity_values;
+          std::vector<Tensor<1,dim> > face_current_velocity_values;
           std::vector<Tensor<1,dim> > mesh_velocity_values;
+          std::vector<Tensor<1,dim> > face_mesh_velocity_values;
 
           std::vector<SymmetricTensor<2,dim> > current_strain_rates;
           std::vector<std::vector<double> > current_composition_values;
@@ -173,6 +184,12 @@ namespace aspect
           MaterialModel::MaterialModelInputs<dim> material_model_inputs;
           MaterialModel::MaterialModelOutputs<dim> material_model_outputs;
 
+          MaterialModel::MaterialModelInputs<dim> face_material_model_inputs;
+          MaterialModel::MaterialModelOutputs<dim> face_material_model_outputs;
+
+          MaterialModel::MaterialModelInputs<dim> neighbor_face_material_model_inputs;
+          MaterialModel::MaterialModelOutputs<dim> neighbor_face_material_model_outputs;
+
           /**
            * Material model inputs and outputs computed at a previous
            * time step's solution, or an extrapolation from previous
@@ -232,16 +249,37 @@ namespace aspect
            *    are trying to assemble a linear system. <b>Not</b> the global finite
            *    element.
            */
-          AdvectionSystem (const FiniteElement<dim> &finite_element);
+          AdvectionSystem (const FiniteElement<dim> &finite_element,
+                           const bool                field_is_discontinuous);
           AdvectionSystem (const AdvectionSystem &data);
 
           /**
            * Local contributions to the global matrix and right hand side
            * that correspond only to the variables listed in local_dof_indices
            */
           FullMatrix<double>          local_matrix;
+          /** Local contributions to the global matrix from the face terms in the
+           * discontinuous Galerkin method. The vectors are of length
+           * GeometryInfo<dim>::max_children_per_face * GeometryInfo<dim>::faces_per_cell
+           * so as to hold one matrix for each possible face or subface of the cell.
+           * The discontinuous Galerkin bilinear form contains terms arising from internal
+           * (to the cell) values and external (to the cell) values.
+           * _int_ext and ext_int hold the terms arising from the pairing between a cell
+           * and its neighbor, while _ext_ext is the pairing of the neighbor's dofs with
+           * themselves. In the continuous Galerkin case, these are unused, and set to size zero.
+           **/
+          std::vector<FullMatrix<double> >         local_matrices_int_ext;
+          std::vector<FullMatrix<double> >         local_matrices_ext_int;
+          std::vector<FullMatrix<double> >         local_matrices_ext_ext;
           Vector<double>              local_rhs;
 
+          /** Denotes which face matrices have actually been assembled in the DG field
+           * assembly. Entries for matrices not used (for example, those corresponding
+           * to non-existent subfaces; or faces being assembled by the neighboring cell)
+           * are set to false.
+           **/
+          std::vector<bool>               assembled_matrices;
+
           /**
            * Indices of those degrees of freedom that actually correspond
            * to the temperature or compositional field. since this structure
@@ -252,6 +290,13 @@ namespace aspect
            * any other variable outside the block we are currently considering)
            */
           std::vector<types::global_dof_index>   local_dof_indices;
+          /** Indices of the degrees of freedom corresponding to the temperature
+           * or composition field on all possible neighboring cells. This is used
+           * in the discontinuous Galerkin method. The outer std::vector has
+           *  length GeometryInfo<dim>::max_children_per_face * GeometryInfo<dim>::faces_per_cell,
+           * and has size zero if in the continuous Galerkin case.
+           **/
+          std::vector<std::vector<types::global_dof_index> >   neighbor_dof_indices;
         };
       }
 

include/aspect/introspection.h

@@ -77,6 +77,18 @@ namespace aspect
        */
       const unsigned int n_components;
 
+      /**
+       * A variable that holds whether the temperature field should use a
+       * discontinuous discretization.
+       */
+      const bool use_discontinuous_temperature_discretization;
+
+      /**
+       * A variable that holds whether the composition field(s) should use a
+       * discontinuous discretization.
+       */
+      const bool use_discontinuous_composition_discretization;
+
       /**
        * A structure that enumerates the vector components of the finite
        * element that correspond to each of the variables in this problem.

include/aspect/parameters.h

@@ -248,6 +248,7 @@ namespace aspect
     unsigned int                   stabilization_alpha;
     double                         stabilization_c_R;
     double                         stabilization_beta;
+    double                         discontinuous_penalty;
     /**
      * @}
      */
@@ -268,6 +269,8 @@ namespace aspect
      */
     unsigned int                   stokes_velocity_degree;
     bool                           use_locally_conservative_discretization;
+    bool                           use_discontinuous_temperature_discretization;
+    bool                           use_discontinuous_composition_discretization;
     unsigned int                   temperature_degree;
     unsigned int                   composition_degree;
     std::string                    pressure_normalization;

include/aspect/simulator.h

@@ -232,6 +232,12 @@ namespace aspect
         bool
         is_temperature () const;
 
+        /**
+         * Return whether this object refers to a field discretized by discontinuous finite elements.
+         */
+        bool
+        is_discontinuous (const Introspection<dim> &introspection) const;
+
         /**
          * Look up the component index for this temperature or compositional
          * field. See Introspection::component_indices for more information.
@@ -643,6 +649,18 @@ namespace aspect
       void
       copy_local_to_global_stokes_system (const internal::Assembly::CopyData::StokesSystem<dim> &data);
 
+      /**
+       * Compute the integrals for one advection matrix and right hand side on
+       * the faces of a single cell.
+       *
+       * This function is implemented in
+       * <code>source/simulator/assembly.cc</code>.
+       */
+      void
+      local_assemble_advection_face_terms(const AdvectionField &advection_field,
+                                          const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                          internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
+                                          internal::Assembly::CopyData::AdvectionSystem<dim> &data);
       /**
        * Compute the integrals for one advection matrix and right hand side on
        * a single cell.
@@ -665,7 +683,8 @@ namespace aspect
        * <code>source/simulator/assembly.cc</code>.
        */
       void
-      copy_local_to_global_advection_system (const internal::Assembly::CopyData::AdvectionSystem<dim> &data);
+      copy_local_to_global_advection_system (const AdvectionField &advection_field,
+                                             const internal::Assembly::CopyData::AdvectionSystem<dim> &data);
 
       /**
        * @}