EquationBase instance contains its own local cell iterators Slaybaugh…

…Lab#10

This is a major design change. It means we will initialize
cell iterators on current process when we instantiate an equation.

src/common/bart_driver.cc

@@ -82,18 +82,12 @@ template <int dim>
 void BartDriver<dim>::setup_system ()
 {
   radio ("setup system");
-  initialize_dealii_objects ();
-}
-
-template <int dim>
-void BartDriver<dim>::initialize_dealii_objects ()
-{
   dof_handler.distribute_dofs (*fe);
-
+  
   local_dofs = dof_handler.locally_owned_dofs ();
   DoFTools::extract_locally_relevant_dofs (dof_handler,
                                            relevant_dofs);
-
+  
   constraints.clear ();
   constraints.reinit (relevant_dofs);
   DoFTools::make_hanging_node_constraints (dof_handler,
@@ -120,6 +114,7 @@ void BartDriver<dim>::initialize_dealii_objects ()
    MPI_COMM_WORLD,
    relevant_dofs);
 
+  itr_ptr = build_iterations (prm, dof_handler, msh_ptr, aqd_ptr, mat_ptr);
   itr_ptr->initialize_system_matrices_vectors (dsp, local_dofs);
 }
 
@@ -171,10 +166,6 @@ void BartDriver<dim>::run ()
   msh_ptr->make_grid (triangulation);
   setup_system ();
   report_system ();
-  // get cell iterators on local processors
-  msh_ptr->get_relevant_cell_iterators (dof_handler,
-                                        local_cells,
-                                        is_cell_at_bd);
   // solve the problem using iterative methods specified in Iterations class
   itr_ptr->solve_problems (local_cells, is_cell_at_bd, sflx_proc);
   if (is_eigen_problem)

src/common/bart_tools.cc

@@ -32,7 +32,8 @@ std_cxx11::shared_ptr<PreconditionerSolver> build_linalg
 
 template <int dim>
 std_cxx11::shared_ptr<Iterations<dim> > build_iterations
-(ParameterHandler &prm,
+(const ParameterHandler &prm,
+ const DoFHandler<dim> &dof_handler,
  const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
  const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
  const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)

src/iteration/in_group_base.h

@@ -10,7 +10,9 @@ class InGroupBase : public IterationBase<dim>
 
   // has to be provided
   virtual void solve_in_group
-  (std_cxx11::shared_ptr<EquationBase<dim> > equ_ptr,
+  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
+   std::vector<bool> &is_cell_at
+   std_cxx11::shared_ptr<EquationBase<dim> > equ_ptrs,
    unsigned int &g);
 
 protected:
@@ -27,7 +29,7 @@ class SourceIteration : public InGroupBase<dim>
   ~SourceIteration ();
 
   void solve_in_group
-  (std_cxx11::shared_ptr<EquationBase<dim> > equ_ptr,
+  (std_cxx11::shared_ptr<EquationBase<dim> > equ_ptrs,
    unsigned int &g);
 };
 

src/iteration/iterations.cc

@@ -16,7 +16,8 @@ using namespace dealii;
 
 template <int dim>
 Iterations<dim>::Iterations
-(ParameterHandler &prm,
+(const ParameterHandler &prm,
+ const DoFHandler<dim> &dof_handler
  const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
  const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
  const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)
@@ -28,9 +29,9 @@ do_nda(prm.get_bool("do NDA"))
   mat_ptr = build_material (prm);
   // vectors containing all the pointers to equations. Size is 2 if NDA is used.
   equ_ptrs.resize (do_nda?2:1);
-  equ_ptrs[0] = build_equation (transport_name, prm, msh_ptr, aqd_ptr, mat_ptr);
+  equ_ptrs[0] = build_equation (transport_name, prm, dof_handler, msh_ptr, aqd_ptr, mat_ptr);
   if (do_nda)
-    equ_ptrs[1] = build_equation ("nda", prm, msh_ptr, aqd_ptr, mat_ptr);
+    equ_ptrs[1] = build_equation ("nda", prm, dof_handler, msh_ptr, aqd_ptr, mat_ptr);
 }
 
 template <int dim>
@@ -42,29 +43,23 @@ template <int dim>
 void Iterations<dim>::initialize_system_matrices_vectors
 (SparsityPatternType &dsp, IndexSet &local_dofs)
 {
-  equ_ptrs[0]->initialize_system_matrices_vectors (dsp, local_dofs);
-  if (do_nda)
-    equ_ptrs[1]->initialize_system_matrices_vectors (dsp, local_dofs);
+  for (unsigned int i=0; i<equ_ptrs.size(); ++i)
+    equ_ptrs[i]->initialize_system_matrices_vectors (dsp, local_dofs);
 }
 
 template <int dim>
-void Iterations<dim>::solve_problems
-(std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
- std::vector<bool> &is_cell_at_bd,
- std::vector<Vector<double> > &sflx_proc)
+void Iterations<dim>::solve_problems (std::vector<Vector<double> > &sflx_proc)
 {
   if (is_eigen_problem)
   {
     std_cxx11::shared_ptr<EigenBase<dim> > pro_ptr = build_eigen_problem (prm);
     pro_ptr->do_iterations (local_cells, is_cell_at_bd, sflx_proc, equ_ptrs);
-    //pro_ptr->get_sflx_proc (sflx_proc);
     pro_ptr->get_keff (keff);
   }
   else
   {
     std_cxx11::shared_ptr<MGBase<dim> > pro_ptr = build_mg_problem (prm);
     pro_ptr->do_iterations (local_cells, is_cell_at_bd, sflx_proc, equ_ptrs);
-    //pro_ptr->get_sflx_proc (sflx_proc);
   }
 }
 

src/iteration/iterations.h

@@ -45,13 +45,15 @@ template <int dim>
 class Iterations
 {
 public:
-  Iterations (ParameterHandler &prm);
+  Iterations (std::string equation_name,
+              const ParameterHandler &prm,
+              const DoFHandler<dim> &dof_handler
+              const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
+              const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
+              const std_cxx11::shared_ptr<MaterialProperties> mat_ptr);
   virtual ~Iterations ();
 
-  void solve_problems
-  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
-   std::vector<bool> &is_cell_at_bd,
-   std::vector<Vector<double> > &sflx_proc);
+  void solve_problems (std::vector<Vector<double> > &sflx_proc);
 
   void initialize_system_matrices_vectors
   (SparsityPatternType &dsp, IndexSet &local_dofs);

src/transport/equation_base.cc

@@ -17,7 +17,8 @@ using namespace dealii;
 template <int dim>
 EquationBase<dim>::EquationBase
 (std::string equation_name,
- ParameterHandler &prm,
+ const ParameterHandler &prm,
+ const DoFHandler<dim> &dof_handler
  const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
  const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
  const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)
@@ -32,7 +33,10 @@ have_reflective_bc(prm.get_bool("have reflective BC")),
 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);
+  // process input for mesh, AQ and material related data
+  process_input (msh_ptr, aqd_ptr, mat_ptr);
+  // get cell iterators and booleans to tell if they are at boundary on this process
+  msh_ptr->get_relevant_cell_iterators (dof_handler, local_cells, is_cell_at_bd);
   alg_ptr = build_linalg (prm, equation_name, n_total_vars)
 }
 
@@ -93,6 +97,15 @@ void EquationBase<dim>::process_input
   }
 }
 
+template <int dim>
+void EquationBase<dim>::initialize_cell_iterators_this_proc
+(const DoFHandler<dim> &dof_handler)
+{
+  msh_ptr->get_relevant_cell_iterators (dof_handler,
+                                        local_cells,
+                                        is_cell_at_bd);
+}
+
 template <int dim>
 void EquationBase<dim>::initialize_system_matrices_vectors
 (SparsityPatternType &dsp, IndexSet &local_dofs)
@@ -440,9 +453,7 @@ void EquationBase<dim>::scale_fiss_transfer_matrices (double keff)
 // generate rhs for equation
 template <int dim>
 void EquationBase<dim>::assemble_linear_form
-(std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
- std::vector<bool> &is_cell_at_bd,
- std::vector<Vector<double> > &sflx_this_proc,
+(std::vector<Vector<double> > &sflx_this_proc,
  unsigned int &g)
 {
   for (unsigned int k=0; k<this->n_total_vars; ++k)
@@ -487,9 +498,7 @@ void EquationBase<dim>::integrate_scattering_linear_form
 
 template <int dim>
 void EquationBase<dim>::assemble_fixed_linear_form
-(std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
- std::vector<bool> &is_cell_at_bd,
- std::vector<Vector<double> > &sflx_prev)
+(std::vector<Vector<double> > &sflx_prev)
 {
   for (unsigned int k=0; k<n_total_vars; ++k)
   {
@@ -531,8 +540,7 @@ void EquationBase<dim>::initialize_preconditioners
 
 template <int dim>
 double EquationBase<dim>::estimate_fiss_source
-(std::vector<Vector<double> > &phis_this_process,
- std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells)
+(std::vector<Vector<double> > &phis_this_process)
 {
   // first, estimate local fission source
   double fiss_source = 0.0;

src/transport/equation_base.h

@@ -46,31 +46,40 @@ class EquationBase
 {
 public:
   // For SN and NDA (NDA still needs quadrature to calculate corrections)
-  EquationBase (ParameterHandler &prm,
+  EquationBase (std::string equation_name,
+                const ParameterHandler &prm,
+                const DoFHandler<dim> &dof_handler
                 const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
                 const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
-                const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)
-
-  // For diffusion (future work for someone)
-  EquationBase (ParameterHandler &prm,
+                const std_cxx11::shared_ptr<MaterialProperties> mat_ptr);
+
+  // TODO: for diffusion (future work for someone)
+  // instead of previous constructor, diffusion overrides a constructor without
+  // asking for AQ instance.
+  EquationBase (std::string equation_name,
+                const ParameterHandler &prm,
+                const DoFHandler<dim> &dof_handler
                 const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
                 const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)
 
   virtual ~EquationBase ();
 
   void run ();
 
-  virtual void assemble_bilinear_form
-  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
-   std::vector<bool> &is_cell_at_bd);
+  virtual void assemble_bilinear_form ();
+
+  virtual void assemble_volume_boundary_bilinear_form ();
+
+  virtual void assemble_interface_bilinear_form ();
+
+  virtual void assemble_linear_form
+  (std::vector<Vector<double> > &sflx_this_proc,
+   unsigned int &g);
+
+  virtual void assemble_fixed_linear_form
+  (std::vector<Vector<double> > &sflx_prev);
 
-  virtual void assemble_volume_boundary_bilinear_form
-  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
-   std::vector<bool> &is_cell_at_bd);
 
-  virtual void assemble_interface_bilinear_form
-  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
-   std::vector<bool> &is_cell_at_bd);
 
   virtual void pre_assemble_cell_matrices
   (typename DoFHandler<dim>::active_cell_iterator &cell,
@@ -190,12 +199,11 @@ class EquationBase
   unsigned int p_order;
   unsigned int global_refinements;
 
-  std::vector<unsigned int> linear_iters;
-
+  std::vector<typename DoFHandler<dim>::active_cell_iterator> local_cells;
+  std::vector<bool> is_cell_at_bd;
   std::vector<types::global_dof_index> local_dof_indices;
   std::vector<types::global_dof_index> neigh_dof_indices;
 
-
   std::vector<Tensor<1, dim> > omega_i;
   std::vector<double> wi;
   std::vector<double> tensor_norms;

src/transport/even_parity.cc

@@ -3,7 +3,8 @@
 template <int dim>
 EvenParity<dim>::EvenParity
 (std::string equation_name,
- ParameterHandler &prm,
+ const ParameterHandler &prm,
+ const DoFHandler<dim> &dof_handler
  const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
  const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
  const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)

src/transport/even_parity.h

@@ -8,7 +8,8 @@ class EvenParity : public EquationBase<dim>
 {
 public:
   EvenParity (std::string equation_name,
-              ParameterHandler &prm,
+              const ParameterHandler &prm,
+              const DoFHandler<dim> &dof_handler
               const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
               const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
               const std_cxx11::shared_ptr<MaterialProperties> mat_ptr);