Added FE objects for NDA SlaybaughLab#5

src/transport/equation_base.cc

@@ -30,9 +30,11 @@ n_azi(prm.get_integer("angular quadrature order")),
 is_eigen_problem(prm.get_bool("do eigenvalue calculations")),
 do_nda(prm.get_bool("do NDA")),
 have_reflective_bc(prm.get_bool("have reflective BC")),
-p_order(prm.get_integer("finite element polynomial degree"))
+p_order(prm.get_integer("finite element polynomial degree")),
+nda_quadrature_order(p_order+3) //this is hard coded
 {
-  this->process_input (msh_ptr, aqd_ptr, mat_ptr);    
+  this->process_input (msh_ptr, aqd_ptr, mat_ptr);
+
 }
 
 template <int dim>
@@ -134,6 +136,25 @@ void EquationBase<dim>::initialize_assembly_related_objects
   dofs_per_cell = fe->dofs_per_cell;
   n_q = q_rule->size();
   n_qf = qf_rule->size();
+
+  if (do_nda)
+  {
+    qc_rule = std_cxx11::shared_ptr<QGauss<dim> > (new QGauss<dim> (nda_quadrature_order));
+    qfc_rule = std_cxx11::shared_ptr<QGauss<dim-1> > (new QGauss<dim-1> (nda_quadrature_order));
+    fvc = std_cxx11::shared_ptr<FEValues<dim> >
+    (new FEValues<dim> (*fe, *qc_rule,
+                        update_values | update_gradients |
+                        update_quadrature_points |
+                        update_JxW_values));
+
+    fvfc = std_cxx11::shared_ptr<FEFaceValues<dim> >
+    (new FEFaceValues<dim> (*fe, *qfc_rule,
+                            update_values | update_gradients |
+                            update_quadrature_points | update_normal_vectors |
+                            update_JxW_values));
+    n_qc = qc_rule->size ();
+    n_qfc = qfc_rule->size ();
+  }
 
   local_dof_indices.resize (dofs_per_cell);
   neigh_dof_indices.resize (dofs_per_cell);

src/transport/equation_base.h

@@ -164,11 +164,16 @@ class EquationBase
   unsigned int get_reflective_direction_index (unsigned int boundary_id,
                                                unsigned int incident_angle_index);
 
+  // "c" in the following quantities means "correction" for NDA use
   std_cxx11::shared_ptr<QGauss<dim> > q_rule;
   std_cxx11::shared_ptr<QGauss<dim-1> > qf_rule;
+  std_cxx11::shared_ptr<QGauss<dim> > qc_rule;
+  std_cxx11::shared_ptr<QGauss<dim-1> > qfc_rule;
   std_cxx11::shared_ptr<FEValues<dim> > fv;
   std_cxx11::shared_ptr<FEFaceValues<dim> > fvf;
   std_cxx11::shared_ptr<FEFaceValues<dim> > fvf_nei;
+  std_cxx11::shared_ptr<FEValues<dim> > fvc;
+  std_cxx11::shared_ptr<FEFaceValues<dim> > fvfc;
 
   double keff;
   double keff_prev_gen;
@@ -179,8 +184,11 @@ class EquationBase
   bool do_nda;
   bool have_reflective_bc;
 
+  const unsigned int nda_quadrature_order;
   unsigned int n_q;
   unsigned int n_qf;
+  unsigned int n_qc;
+  unsigned int n_qfc;
   unsigned int dofs_per_cell;
 
   unsigned int n_dir;