diff --git a/src/common/bart_tools.cc b/src/common/bart_tools.cc
index 73215a873..1cfeaa87f 100644
--- a/src/common/bart_tools.cc
+++ b/src/common/bart_tools.cc
@@ -7,6 +7,9 @@
 #include "../equations/even_parity.h"
 #include "../aqdata/aq_lsgc"
 
+/** \brief Function to build finite element for general dimensions specified by 
+ * user.
+ */
 template <int dim>
 FE_Poly<TensorProductPolynomials<dim>,dim,dim>* build_finite_element
 (ParameterHandler &prm)
@@ -20,6 +23,8 @@ FE_Poly<TensorProductPolynomials<dim>,dim,dim>* build_finite_element
     return new FE_Q<dim> (p_order);
 }
 
+/** \brief Function to build mesh in calculations for general dimensions
+ */
 template <int dim>
 std_cxx11::shared_ptr<MeshGenerator<dim> > build_mesh (ParameterHandler &prm)
 {
@@ -29,6 +34,8 @@ std_cxx11::shared_ptr<MeshGenerator<dim> > build_mesh (ParameterHandler &prm)
   return mesh_class;
 }
 
+/** \brief Function to build pointer to MaterialProperties class.
+ */
 std_cxx11::shared_ptr<MaterialProperties> build_material (ParameterHandler &prm)
 {
   std_cxx11::shared_ptr<MaterialProperties> material_class =
@@ -50,6 +57,13 @@ std_cxx11::shared_ptr<IterationBase<dim> > build_transport_iteration
   return iteration_class;
 }
 
+/** \brief Build specific transport model
+ * 
+ * \parameter msh_ptr shared_ptr for MeshGenerator<dim> instance
+ * \parameter aqd_ptr shared_ptr for AQBase<dim> instance
+ * \parameter mat_ptr shared_ptr for MaterialProperties instance
+ * \return a shared pointer to transport model
+ */
 template <int dim>
 std_cxx11::shared_ptr<TransportBase<dim> > build_transport_model
 (ParameterHandler &prm,
@@ -64,13 +78,15 @@ std_cxx11::shared_ptr<TransportBase<dim> > build_transport_model
   if (transport_model_name=="ep")
     transport_class =
     std_cxx11::shared_ptr<TransportBase<dim> >
-    (new EvenParity<dim> (prm,
-                          msh_ptr,
-                          aqd_ptr,
-                          mat_ptr));
+    (new EvenParity<dim> (prm, msh_ptr, aqd_ptr, mat_ptr));
   return transport_class;
 }
 
+/** \brief Function to build angular quadrature for general dimensions
+ *
+ * \parameter prm A processed ParameterHandler instance
+ * \return a shared pointer to specific angular quadrature
+ */
 template <int dim>
 std_cxx11::shared_ptr<AQBase<dim> >
 build_aq_model (ParameterHandler &prm)
diff --git a/src/transport/equation_base.cc b/src/transport/equation_base.cc
index 00b14a8c0..20eefb5fd 100644
--- a/src/transport/equation_base.cc
+++ b/src/transport/equation_base.cc
@@ -89,19 +89,19 @@ void EquationBase<dim>::process_input
 }
 
 template <int dim>
-void EquationBase<dim>::assemble_ho_system
+void EquationBase<dim>::assemble_system
 (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
  std::vector<bool> &is_cell_at_bd)
 {
   radio ("Assemble volumetric bilinear forms");
-  assemble_ho_volume_boundary (local_cells, is_cell_at_bd);
+  assemble_volume_boundary (local_cells, is_cell_at_bd);
 
   if (discretization=="dfem")
   {
     AssertThrow (transport_model_name=="ep",
                  ExcMessage("DFEM is only implemented for even parity"));
     radio ("Assemble cell interface bilinear forms for DFEM");
-    assemble_ho_interface (local_cells);
+    assemble_interface (local_cells);
   }
 }
 
@@ -140,7 +140,7 @@ void EquationBase<dim>::initialize_assembly_related_objects
 }
 
 template <int dim>
-void EquationBase<dim>::assemble_ho_volume_boundary
+void EquationBase<dim>::assemble_volume_boundary
 (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
  std::vector<bool> &is_cell_at_bd)
 {
@@ -178,10 +178,9 @@ void EquationBase<dim>::assemble_ho_volume_boundary
       integrate_cell_bilinear_form (fv,
                                     cell,
                                     local_mat,
-                                    i_dir,
-                                    g,
                                     streaming_at_qp,
-                                    collision_at_qp);
+                                    collision_at_qp,
+                                    g, i_dir);
 
       if (is_cell_at_bd[ic])
         for (unsigned int fn=0; fn<GeometryInfo<dim>::faces_per_cell; ++fn)
@@ -192,8 +191,7 @@ void EquationBase<dim>::assemble_ho_volume_boundary
                                               cell,
                                               fn,
                                               local_mat,
-                                              i_dir,
-                                              g);
+                                              g, i_dir);
           }
       
       if (k==0)
@@ -228,45 +226,60 @@ void EquationBase<dim>::integrate_cell_bilinear_form
 (const std_cxx11::shared_ptr<FEValues<dim> > fv,
  typename DoFHandler<dim>::active_cell_iterator &cell,
  FullMatrix<double> &cell_matrix,
- unsigned int &i_dir,
- unsigned int &g,
  std::vector<std::vector<FullMatrix<double> > > &streaming_at_qp,
- std::vector<FullMatrix<double> > &collision_at_qp)
+ std::vector<FullMatrix<double> > &collision_at_qp,
+ const unsigned int &g,
+ const unsigned int &i_dir=0)
 {
 }
 
-// The following is a virtual function for integraing boundary bilinear form;
-// It must be overriden
+/** \brief Integrator for boundary weak form per boundary face per angular/group
+ *
+ * The function is a virtual function. For diffusion-like system, i_dir is set 
+ * to 0 by default.
+ */
 template <int dim>
 void EquationBase<dim>::integrate_boundary_bilinear_form
 (const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
  typename DoFHandler<dim>::active_cell_iterator &cell,
  unsigned int &fn,/*face number*/
  FullMatrix<double> &cell_matrix,
- unsigned int &i_dir,
- unsigned int &g)
+ const unsigned int &g,
+ const unsigned int &i_dir=0)
 {// this is a virtual function
 }
 
+/** \brief Right hand side integrator specifically for reflective boundary.
+ *
+ */
 template <int dim>
 void EquationBase<dim>::integrate_reflective_boundary_linear_form
 (const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
  typename DoFHandler<dim>::active_cell_iterator &cell,
  unsigned int &fn,/*face number*/
  std::vector<Vector<double> > &cell_rhses,
- unsigned int &i_dir,
- unsigned int &g)
+ const unsigned int &g,
+ const unsigned int &i_dir)
 {// this is a virtual function
 }
 
+/** \brief Interface weak form assembly driver.
+ * Member function used to assemble interface weak forms. The main functionality
+ * is to go through all non-boundary interfaces of the cells owned on current 
+ * processor and assemble the weak form using interface assembler. 
+ *
+ * There is no need to override this function for SN calculations. Yet, for PN,
+ * diffusion etc., this function must be overriden to correctly take care of the
+ * angular component.
+ */
 template <int dim>
-void EquationBase<dim>::assemble_ho_interface
+void EquationBase<dim>::assemble_interface
 (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells)
 {
-  FullMatrix<double> vp_up (dofs_per_cell, dofs_per_cell);
-  FullMatrix<double> vp_un (dofs_per_cell, dofs_per_cell);
-  FullMatrix<double> vn_up (dofs_per_cell, dofs_per_cell);
-  FullMatrix<double> vn_un (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> vi_ui (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> vi_ue (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> ve_ui (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> ve_ue (dofs_per_cell, dofs_per_cell);
 
   for (unsigned int k=0; k<n_total_ho_vars; ++k)
   {
@@ -288,37 +301,44 @@ void EquationBase<dim>::assemble_ho_interface
           neigh->get_dof_indices (neigh_dof_indices);
           fvf_nei->reinit (neigh, cell->neighbor_face_no(fn));
 
-          vp_up = 0;
-          vp_un = 0;
-          vn_up = 0;
-          vn_un = 0;
+          vi_ui = 0;
+          vi_ue = 0;
+          ve_ui = 0;
+          ve_ue = 0;
 
           integrate_interface_bilinear_form (fvf, fvf_nei,/*FEFaceValues objects*/
                                              cell, neigh,/*cell iterators*/
                                              fn,
-                                             i_dir, g,/*specific component*/
-                                             vp_up, vp_un, vn_up, vn_un);
+                                             vi_ui, vi_ue, ve_ui, ve_ue,
+                                             g, i_dir/*specific component*/);
           vec_ho_sys[k]->add (local_dof_indices,
                               local_dof_indices,
-                              vp_up);
+                              vi_ui);
 
           vec_ho_sys[k]->add (local_dof_indices,
                               neigh_dof_indices,
-                              vp_un);
+                              vi_ue);
 
           vec_ho_sys[k]->add (neigh_dof_indices,
                               local_dof_indices,
-                              vn_up);
+                              ve_ui);
 
           vec_ho_sys[k]->add (neigh_dof_indices,
                               neigh_dof_indices,
-                              vn_un);
+                              ve_ue);
         }// target faces
     }
     vec_ho_sys[k]->compress(VectorOperation::add);
   }// component
 }
 
+/** \brief Virtual function for interface integrator.
+ * When DFEM is used, this function can be overridden as interface weak form
+ * assembler per face per angular and group component.
+ *
+ * When overridden for diffusion calculations, direction component is set to be
+ * zero by default
+ */
 // The following is a virtual function for integrating DG interface for HO system
 // it must be overriden
 template <int dim>
@@ -328,12 +348,12 @@ void EquationBase<dim>::integrate_interface_bilinear_form
  typename DoFHandler<dim>::active_cell_iterator &cell,
  typename DoFHandler<dim>::cell_iterator &neigh,/*cell iterator for cell*/
  unsigned int &fn,/*concerning face number in local cell*/
- unsigned int &i_dir,
- unsigned int &g,
- FullMatrix<double> &vp_up,
- FullMatrix<double> &vp_un,
- FullMatrix<double> &vn_up,
- FullMatrix<double> &vn_un)
+ FullMatrix<double> &vi_ui,
+ FullMatrix<double> &vi_ue,
+ FullMatrix<double> &ve_ui,
+ FullMatrix<double> &ve_ue,
+ const unsigned int &g,
+ const unsigned int &i_dir=0)
 {
 }
 
diff --git a/src/transport/equation_base.h b/src/transport/equation_base.h
index f3b2a4f7e..54ca3d7d4 100644
--- a/src/transport/equation_base.h
+++ b/src/transport/equation_base.h
@@ -54,6 +54,15 @@ class EquationBase
   
   void run ();
   
+  virtual void assemble_volume_boundary
+  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
+   std::vector<bool> &is_cell_at_bd);
+  virtual void assemble_interface
+  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells);
+  virtual void assemble_system
+  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
+   std::vector<bool> &is_cell_at_bd);
+  
   virtual void pre_assemble_cell_matrices
   (const std_cxx11::shared_ptr<FEValues<dim> > fv,
    typename DoFHandler<dim>::active_cell_iterator &cell,
@@ -64,26 +73,26 @@ class EquationBase
   (const std_cxx11::shared_ptr<FEValues<dim> > fv,
    typename DoFHandler<dim>::active_cell_iterator &cell,
    FullMatrix<double> &cell_matrix,
-   unsigned int &i_dir,
-   unsigned int &g,
    std::vector<std::vector<FullMatrix<double> > > &streaming_at_qp,
-   std::vector<FullMatrix<double> > &collision_at_qp);
+   std::vector<FullMatrix<double> > &collision_at_qp,
+   const unsigned int &g,
+   const unsigned int &i_dir=0);
   
   virtual void integrate_boundary_bilinear_form
   (const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
    typename DoFHandler<dim>::active_cell_iterator &cell,
    unsigned int &fn,/*face number*/
    FullMatrix<double> &cell_matrix,
-   unsigned int &i_dir,
-   unsigned int &g);
+   const unsigned int &g,
+   const unsigned int &i_dir=0);
   
   virtual void integrate_reflective_boundary_linear_form
   (const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
    typename DoFHandler<dim>::active_cell_iterator &cell,
    unsigned int &fn,/*face number*/
    std::vector<Vector<double> > &cell_rhses,
-   unsigned int &i_dir,
-   unsigned int &g);
+   const unsigned int &g,
+   const unsigned int &i_dir=0);
   
   virtual void integrate_interface_bilinear_form
   (const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
@@ -118,9 +127,7 @@ class EquationBase
   void setup_boundary_ids ();
   void get_cell_mfps (unsigned int &material_id, double &cell_dimension,
                       std::vector<double> &local_mfps);
-  void assemble_ho_volume_boundary ();
-  void assemble_ho_interface ();
-  void assemble_ho_system ();
+  
   void process_input (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);
diff --git a/src/transport/even_parity.cc b/src/transport/even_parity.cc
index 5ba1b8a2f..f3ba6700b 100644
--- a/src/transport/even_parity.cc
+++ b/src/transport/even_parity.cc
@@ -52,10 +52,10 @@ void EvenParity<dim>::integrate_cell_bilinear_form
 (const std_cxx11::shared_ptr<FEValues<dim> > fv,
  typename DoFHandler<dim>::active_cell_iterator &cell,
  FullMatrix<double> &cell_matrix,
- unsigned int &i_dir,
- unsigned int &g,
  std::vector<std::vector<FullMatrix<double> > > &streaming_at_qp,
- std::vector<FullMatrix<double> > &collision_at_qp)
+ std::vector<FullMatrix<double> > &collision_at_qp,
+ const unsigned int &g,
+ const unsigned int &i_dir)
 {
   unsigned int mid = cell->material_id ();
   for (unsigned int qi=0; qi<this->n_q; ++qi)
@@ -74,8 +74,8 @@ void EvenParity<dim>::integrate_boundary_bilinear_form
  typename DoFHandler<dim>::active_cell_iterator &cell,
  unsigned int &fn,/*face number*/
  FullMatrix<double> &cell_matrix,
- unsigned int &i_dir,
- unsigned int &g)
+ const unsigned int &g,
+ const unsigned int &i_dir)
 {
   unsigned int bd_id = cell->face(fn)->boundary_id ();
   const Tensor<1,dim> vec_n = fvf->normal_vector(0);
@@ -118,12 +118,12 @@ void EvenParity<dim>::integrate_interface_bilinear_form
  typename DoFHandler<dim>::active_cell_iterator &cell,
  typename DoFHandler<dim>::cell_iterator &neigh,/*cell iterator for cell*/
  unsigned int &fn,/*concerning face number in local cell*/
- unsigned int &i_dir,
- unsigned int &g,
  FullMatrix<double> &vp_up,
  FullMatrix<double> &vp_un,
  FullMatrix<double> &vn_up,
- FullMatrix<double> &vn_un)
+ FullMatrix<double> &vn_un,
+ const unsigned int &g,
+ const unsigned int &i_dir)
 {
   const Tensor<1,dim> vec_n = fvf->normal_vector (0);
   unsigned int mid = cell->material_id ();
diff --git a/src/transport/even_parity.h b/src/transport/even_parity.h
index d86f3a6c2..881c96d6e 100644
--- a/src/transport/even_parity.h
+++ b/src/transport/even_parity.h
@@ -23,18 +23,18 @@ class EvenParity : public TransportBase<dim>
   (const std_cxx11::shared_ptr<FEValues<dim> > fv,
    typename DoFHandler<dim>::active_cell_iterator &cell,
    FullMatrix<double> &cell_matrix,
-   unsigned int &i_dir,
-   unsigned int &g,
    std::vector<std::vector<FullMatrix<double> > > &streaming_at_qp,
-   std::vector<FullMatrix<double> > &collision_at_qp);
+   std::vector<FullMatrix<double> > &collision_at_qp,
+   const unsigned int &g,
+   const unsigned int &i_dir);
   
   void integrate_boundary_bilinear_form
   (const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
    typename DoFHandler<dim>::active_cell_iterator &cell,
    unsigned int &fn,/*face number*/
    FullMatrix<double> &cell_matrix,
-   unsigned int &i_dir,
-   unsigned int &g);
+   const unsigned int &g,
+   const unsigned int &i_dir);
   
   void integrate_interface_bilinear_form
   (const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
@@ -42,21 +42,22 @@ class EvenParity : public TransportBase<dim>
    typename DoFHandler<dim>::active_cell_iterator &cell,
    typename DoFHandler<dim>::cell_iterator &neigh,/*cell iterator for cell*/
    unsigned int &fn,/*concerning face number in local cell*/
-   unsigned int &i_dir,
-   unsigned int &g,
    FullMatrix<double> &vp_up,
    FullMatrix<double> &vp_un,
    FullMatrix<double> &vn_up,
-   FullMatrix<double> &vn_un);
+   FullMatrix<double> &vn_un,
+   const unsigned int &g,
+   const unsigned int &i_dir);
   
   void generate_ho_fixed_source
   (std::vector<PETScWrappers::MPI::Vector*> &vec_ho_fixed_rhs,
    std::vector<Vector<double> > &sflx_this_proc);
+  
   void generate_ho_rhs
   (std::vector<PETScWrappers::MPI::Vector*> &vec_ho_rhs,
    std::vector<PETScWrappers::MPI::Vector*> &vec_ho_fixed_rhs,
    std::vector<Vector> &sflx_this_proc);
-
+  
 private:
   double c_penalty;
   std::vector<double> tensor_norms;