diff --git a/.gitignore b/.gitignore
index e87606706..4206447b1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,4 @@ doc/html/**
 *.gnuplot
 *.swp
 *.swo
+*.vtk
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b7735883b..6e9d2d687 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -30,6 +30,17 @@ IF(NOT ${deal.II_FOUND})
     )
 ENDIF()
 
+option(USE_LD_GOLD "Use GNU gold linker" ON)
+
+if(USE_LD_GOLD AND "${CMAKE_C_COMPILER_ID}" STREQUAL "GNU")
+  execute_process(COMMAND ${CMAKE_C_COMPILER} -fuse-ld=gold -Wl,--version OUTPUT_VARIABLE stdout ERROR_QUIET)
+  if("${stdout}" MATCHES "GNU gold")
+    set(CMAKE_CXX_FLAGS "${CMAKE_C_FLAGS} -fuse-ld=gold")
+    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fuse-ld=gold")
+  else()
+    message(WARNING "GNU gold linker isn't available, using the default system linker.")
+  endif()
+endif()
 
 DEAL_II_INITIALIZE_CACHED_VARIABLES()
 
diff --git a/TODO b/TODO
index a8423fec0..9f99d9458 100644
--- a/TODO
+++ b/TODO
@@ -3,6 +3,8 @@ https://stackoverflow.com/a/1776295/11047627
 
 Remove pointers to FEValues. Just create the object through initializer list.
 
+Alot of copied code between weak and strong form. Reduce copying somehow.
+
 
 Lost adjoint consistency by computing F* = F*(Uin, Ubc) instead of F*= F*(Ubc, Ubc). See comment in assemble_boundary_term_implicit.
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 333871fc2..1229cbd12 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -41,7 +41,7 @@ foreach(dim RANGE 1 3)
 #    # See one of the following links (translate chinese to english using Google)
 #    # https://zhuanlan.zhihu.com/p/21385662
 #    # https://www.getit01.com/p2018020821385662/
-#    target_link_libraries(${PHiLiPLib} tbb m mpi mpi_cxx)
+#    target_link_libraries(${PHiLiPLib} m mpi mpi_cxx)
 
 
     # #####################################################
@@ -62,7 +62,7 @@ foreach(dim RANGE 1 3)
     string(CONCAT TestsLib Tests_${dim}D)
     target_link_libraries(${MAIN_TARGET} ${ParameterLib})
     target_link_libraries(${MAIN_TARGET} ${TestsLib})
-    target_link_libraries(${MAIN_TARGET} tbb m mpi mpi_cxx)
+    target_link_libraries(${MAIN_TARGET} m mpi mpi_cxx)
     DEAL_II_SETUP_TARGET(${MAIN_TARGET})
     
     unset(MAIN_TARGET)
@@ -77,7 +77,7 @@ endforeach()
 #     temporary_test.cpp)
 # 
 # add_executable(${SMALL_TEST} ${SMALL_TEST_SRC})
-#     target_link_libraries(${SMALL_TEST} tbb m mpi mpi_cxx)
+#     target_link_libraries(${SMALL_TEST} m mpi mpi_cxx)
 # 
 # DEAL_II_SETUP_TARGET(${SMALL_TEST})
 # message("Adding executable " ${SMALL_TEST} " with files " ${SMALL_TEST_SRC} "\n")
diff --git a/src/dg/CMakeLists.txt b/src/dg/CMakeLists.txt
index 654ec78e1..1c7e24e55 100644
--- a/src/dg/CMakeLists.txt
+++ b/src/dg/CMakeLists.txt
@@ -1,6 +1,7 @@
 set(DG_SOURCE
     dg.cpp
-    assemble_implicit.cpp
+    weak_dg.cpp
+    strong_dg.cpp
     )
 
 foreach(dim RANGE 1 3)
diff --git a/src/dg/dg.cpp b/src/dg/dg.cpp
index eed2b252c..69d7a68e0 100644
--- a/src/dg/dg.cpp
+++ b/src/dg/dg.cpp
@@ -23,7 +23,7 @@
 #include <deal.II/meshworker/loop.h>
 
 
-// Finally, we take our exact solution from the library as well as quadrature
+// Finally, we take our exact solution from the library as well as volume_quadrature
 // and additional tools.
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/data_out_dof_data.h>
@@ -52,18 +52,34 @@ ::create_discontinuous_galerkin(
     //    return new DG<dim,1,real>(parameters_input, degree);
     //}
 
-    if (pde_type == PDE_enum::advection) {
-        return std::make_shared< DG<dim,1,real> >(parameters_input, degree);
-    } else if (pde_type == PDE_enum::advection_vector) {
-        return std::make_shared< DG<dim,2,real> >(parameters_input, degree);
-    } else if (pde_type == PDE_enum::diffusion) {
-        return std::make_shared< DG<dim,1,real> >(parameters_input, degree);
-    } else if (pde_type == PDE_enum::convection_diffusion) {
-        return std::make_shared< DG<dim,1,real> >(parameters_input, degree);
-    } else if (pde_type == PDE_enum::burgers_inviscid) {
-        return std::make_shared< DG<dim,dim,real> >(parameters_input, degree);
-    } else if (pde_type == PDE_enum::euler) {
-        return std::make_shared< DG<dim,dim+2,real> >(parameters_input, degree);
+    if (parameters_input->use_weak_form) {
+        if (pde_type == PDE_enum::advection) {
+            return std::make_shared< DGWeak<dim,1,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::advection_vector) {
+            return std::make_shared< DGWeak<dim,2,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::diffusion) {
+            return std::make_shared< DGWeak<dim,1,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::convection_diffusion) {
+            return std::make_shared< DGWeak<dim,1,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::burgers_inviscid) {
+            return std::make_shared< DGWeak<dim,dim,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::euler) {
+            return std::make_shared< DGWeak<dim,dim+2,real> >(parameters_input, degree);
+        }
+    } else {
+        if (pde_type == PDE_enum::advection) {
+            return std::make_shared< DGStrong<dim,1,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::advection_vector) {
+            return std::make_shared< DGStrong<dim,2,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::diffusion) {
+            return std::make_shared< DGStrong<dim,1,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::convection_diffusion) {
+            return std::make_shared< DGStrong<dim,1,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::burgers_inviscid) {
+            return std::make_shared< DGStrong<dim,dim,real> >(parameters_input, degree);
+        } else if (pde_type == PDE_enum::euler) {
+            return std::make_shared< DGStrong<dim,dim+2,real> >(parameters_input, degree);
+        }
     }
     std::cout << "Can't create DGBase in create_discontinuous_galerkin(). Invalid PDE type: " << pde_type << std::endl;
     return nullptr;
@@ -81,51 +97,231 @@ DGBase<dim,real>::DGBase(
     , fe_dg(degree)
     , fe_system(fe_dg, nstate)
     , all_parameters(parameters_input)
-    , quadrature (degree+1)
+    , oned_quadrature (degree+1)
+    , volume_quadrature (degree+1)
     , face_quadrature (degree+1)
-{
-
-    const dealii::UpdateFlags update_flags =
-        dealii::update_values | dealii::update_gradients
-        | dealii::update_quadrature_points | dealii::update_JxW_values;
-    const dealii::UpdateFlags face_update_flags =
-        dealii::update_values | dealii::update_gradients
-        | dealii::update_quadrature_points | dealii::update_JxW_values
-        | dealii::update_normal_vectors;
-    const dealii::UpdateFlags neighbor_face_update_flags =
-        dealii::update_values | dealii::update_gradients;
-
-    fe_values_cell          = new dealii::FEValues<dim,dim>
-        (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::quadrature, update_flags);
-    fe_values_face_int      = new dealii::FEFaceValues<dim,dim>
-        (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::face_quadrature, face_update_flags);
-    fe_values_subface_int   = new dealii::FESubfaceValues<dim,dim>
-        (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::face_quadrature, face_update_flags);
-    fe_values_face_ext      = new dealii::FEFaceValues<dim,dim>
-        (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::face_quadrature, neighbor_face_update_flags);
-}
+{ }
 
 // Destructor
 template <int dim, typename real>
-DGBase<dim,real>::~DGBase ()
-{ 
-    std::cout << "Destructing DGBase..." << std::endl;
-    delete_fe_values();
-}
+DGBase<dim,real>::~DGBase () { }
 
 template <int dim, typename real>
-void DGBase<dim,real>::delete_fe_values ()
+void DGBase<dim,real>::allocate_system ()
 {
-    if (fe_values_cell          != NULL) delete fe_values_cell;
-    if (fe_values_face_int      != NULL) delete fe_values_face_int;
-    if (fe_values_subface_int   != NULL) delete fe_values_subface_int;
-    if (fe_values_face_ext      != NULL) delete fe_values_face_ext;
-    fe_values_cell          = NULL; 
-    fe_values_face_int      = NULL;
-    fe_values_subface_int   = NULL;
-    fe_values_face_ext      = NULL;
+    std::cout << std::endl << "Allocating DGWeak system and initializing FEValues" << std::endl;
+    // This function allocates all the necessary memory to the 
+    // system matrices and vectors.
+
+    DGBase<dim,real>::dof_handler.initialize(*DGBase<dim,real>::triangulation, DGBase<dim,real>::fe_system);
+    //DGBase<dim,real>::dof_handler.initialize(*DGBase<dim,real>::triangulation, DGBase<dim,real>::fe_system);
+    // Allocates memory from triangulation and finite element space
+    // Use fe_system since it will have the (fe_system.n_dofs)*nstate
+    DGBase<dim,real>::dof_handler.distribute_dofs(DGBase<dim,real>::fe_system);
+
+    //std::vector<unsigned int> block_component(nstate,0);
+    //dealii::DoFRenumbering::component_wise(DGBase<dim,real>::dof_handler, block_component);
+
+    // Allocate matrix
+    unsigned int n_dofs = DGBase<dim,real>::dof_handler.n_dofs();
+    //DynamicSparsityPattern dsp(n_dofs, n_dofs);
+    DGBase<dim,real>::sparsity_pattern.reinit(n_dofs, n_dofs);
+
+    dealii::DoFTools::make_flux_sparsity_pattern(DGBase<dim,real>::dof_handler, DGBase<dim,real>::sparsity_pattern);
+
+    DGBase<dim,real>::system_matrix.reinit(DGBase<dim,real>::sparsity_pattern);
+
+    // Allocate vectors
+    DGBase<dim,real>::solution.reinit(n_dofs);
+    DGBase<dim,real>::right_hand_side.reinit(n_dofs);
+
 }
 
+
+template <int dim, typename real>
+void DGBase<dim,real>::assemble_residual_dRdW ()
+{
+    DGBase<dim,real>::system_matrix = 0;
+    DGBase<dim,real>::right_hand_side = 0;
+
+    // For now assume same polynomial degree across domain
+    const unsigned int dofs_per_cell = DGBase<dim,real>::dof_handler.get_fe().dofs_per_cell;
+    std::vector<dealii::types::global_dof_index> current_dofs_indices (dofs_per_cell);
+    std::vector<dealii::types::global_dof_index> neighbor_dofs_indices (dofs_per_cell);
+
+    // ACTIVE cells, therefore, no children
+    typename dealii::DoFHandler<dim>::active_cell_iterator
+        current_cell = DGBase<dim,real>::dof_handler.begin_active(),
+        endc = DGBase<dim,real>::dof_handler.end();
+
+    unsigned int n_cell_visited = 0;
+    unsigned int n_face_visited = 0;
+
+    dealii::FEValues<dim,dim>        fe_values_cell (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::volume_quadrature, this->update_flags);
+    dealii::FEFaceValues<dim,dim>    fe_values_face_int (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::face_quadrature, this->face_update_flags);
+    dealii::FESubfaceValues<dim,dim> fe_values_subface_int (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::face_quadrature, this->face_update_flags);
+    dealii::FEFaceValues<dim,dim>    fe_values_face_ext (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::face_quadrature, this->neighbor_face_update_flags);
+
+    for (; current_cell!=endc; ++current_cell) {
+        // std::cout << "Current cell index: " << current_cell->index() << std::endl;
+        n_cell_visited++;
+
+        // Local vector contribution from each cell
+        dealii::Vector<double> current_cell_rhs (dofs_per_cell); // Defaults to 0.0 initialization
+
+        fe_values_cell.reinit (current_cell);
+        current_cell->get_dof_indices (current_dofs_indices);
+
+        assemble_cell_terms_implicit (fe_values_cell, current_dofs_indices, current_cell_rhs);
+
+        for (unsigned int iface=0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface) {
+
+            typename dealii::DoFHandler<dim>::face_iterator current_face = current_cell->face(iface);
+            typename dealii::DoFHandler<dim>::cell_iterator neighbor_cell = current_cell->neighbor(iface);
+
+            // See tutorial step-30 for breakdown of 4 face cases
+
+            // Case 1:
+            // Face at boundary
+            if (current_face->at_boundary()) {
+
+                n_face_visited++;
+
+                fe_values_face_int.reinit (current_cell, iface);
+                const unsigned int degree_current = DGBase<dim,real>::fe_system.tensor_degree();
+                const unsigned int deg1sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
+                const unsigned int normal_direction = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
+                const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction);
+
+                real penalty = deg1sq / vol_div_facearea1;
+                //penalty = 1;//99;
+
+                // Need to somehow get boundary type from the mesh
+                assemble_boundary_term_implicit (fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs);
+
+            // Case 2:
+            // Neighbour is finer occurs if the face has children
+            // This is because we are looping over the current_cell's face, so 2, 4, and 6 faces.
+            } else if (current_face->has_children()) {
+                std::cout << "SHOULD NOT HAPPEN!!!!!!!!!!!! I haven't put in adaptatation yet" << std::endl;
+
+                dealii::Vector<double> neighbor_cell_rhs (dofs_per_cell); // Defaults to 0.0 initialization
+                Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError());
+
+                // Obtain cell neighbour
+                const unsigned int neighbor_face_no = current_cell->neighbor_face_no(iface);
+
+                for (unsigned int subface_no=0; subface_no < current_face->number_of_children(); ++subface_no) {
+
+                    n_face_visited++;
+
+                    typename dealii::DoFHandler<dim>::cell_iterator neighbor_child_cell = current_cell->neighbor_child_on_subface (iface, subface_no);
+
+                    Assert (!neighbor_child_cell->has_children(), dealii::ExcInternalError());
+
+                    neighbor_child_cell->get_dof_indices (neighbor_dofs_indices);
+
+                    fe_values_subface_int.reinit (current_cell, iface, subface_no);
+                    fe_values_face_ext.reinit (neighbor_child_cell, neighbor_face_no);
+
+                    const unsigned int normal_direction1 = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
+                    const unsigned int normal_direction2 = dealii::GeometryInfo<dim>::unit_normal_direction[neighbor_face_no];
+                    const unsigned int degree_current = DGBase<dim,real>::fe_system.tensor_degree();
+                    const unsigned int deg1sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
+                    const unsigned int deg2sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
+
+                    //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children();
+                    const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1);
+                    const real vol_div_facearea2 = neighbor_child_cell->extent_in_direction(normal_direction2);
+
+                    const real penalty1 = deg1sq / vol_div_facearea1;
+                    const real penalty2 = deg2sq / vol_div_facearea2;
+                    
+                    real penalty = 0.5 * ( penalty1 + penalty2 );
+                    //penalty = 1;
+
+                    assemble_face_term_implicit (
+                        fe_values_subface_int, fe_values_face_ext,
+                        penalty,
+                        current_dofs_indices, neighbor_dofs_indices,
+                        current_cell_rhs, neighbor_cell_rhs);
+
+                    // Add local contribution from neighbor cell to global vector
+                    for (unsigned int i=0; i<dofs_per_cell; ++i) {
+                        DGBase<dim,real>::right_hand_side(neighbor_dofs_indices[i]) += neighbor_cell_rhs(i);
+                    }
+                }
+
+            // Case 3:
+            // Neighbor cell is NOT coarser
+            // Therefore, they have the same coarseness, and we need to choose one of them to do the work
+            } else if (
+                !current_cell->neighbor_is_coarser(iface) &&
+                    // Cell with lower index does work
+                    (neighbor_cell->index() > current_cell->index() || 
+                    // If both cells have same index
+                    // See https://www.dealii.org/developer/doxygen/deal.II/classTriaAccessorBase.html#a695efcbe84fefef3e4c93ee7bdb446ad
+                    // then cell at the lower level does the work
+                        (neighbor_cell->index() == current_cell->index() && current_cell->level() < neighbor_cell->level())
+                    ) )
+            {
+                n_face_visited++;
+
+                dealii::Vector<double> neighbor_cell_rhs (dofs_per_cell); // Defaults to 0.0 initialization
+
+                Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError());
+                typename dealii::DoFHandler<dim>::cell_iterator neighbor_cell = current_cell->neighbor(iface);
+
+                neighbor_cell->get_dof_indices (neighbor_dofs_indices);
+
+                const unsigned int neighbor_face_no = current_cell->neighbor_of_neighbor(iface);
+
+                const unsigned int normal_direction1 = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
+                const unsigned int normal_direction2 = dealii::GeometryInfo<dim>::unit_normal_direction[neighbor_face_no];
+                const unsigned int degree_current = DGBase<dim,real>::fe_system.tensor_degree();
+                const unsigned int deg1sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
+                const unsigned int deg2sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
+
+                //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children();
+                const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1);
+                const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2);
+
+                const real penalty1 = deg1sq / vol_div_facearea1;
+                const real penalty2 = deg2sq / vol_div_facearea2;
+                
+                real penalty = 0.5 * ( penalty1 + penalty2 );
+                //penalty = 1;//99;
+
+                fe_values_face_int.reinit (current_cell, iface);
+                fe_values_face_ext.reinit (neighbor_cell, neighbor_face_no);
+                assemble_face_term_implicit (
+                        fe_values_face_int, fe_values_face_ext,
+                        penalty,
+                        current_dofs_indices, neighbor_dofs_indices,
+                        current_cell_rhs, neighbor_cell_rhs);
+
+                // Add local contribution from neighbor cell to global vector
+                for (unsigned int i=0; i<dofs_per_cell; ++i) {
+                    DGBase<dim,real>::right_hand_side(neighbor_dofs_indices[i]) += neighbor_cell_rhs(i);
+                }
+            } else {
+            // Do nothing
+            }
+            // Case 4: Neighbor is coarser
+            // Do nothing.
+            // The face contribution from the current cell will appear then the coarse neighbor checks for subfaces
+
+        } // end of face loop
+
+        for (unsigned int i=0; i<dofs_per_cell; ++i) {
+            DGBase<dim,real>::right_hand_side(current_dofs_indices[i]) += current_cell_rhs(i);
+        }
+
+    } // end of cell loop
+
+} // end of assemble_system_implicit ()
+
+
 template <int dim, typename real>
 double DGBase<dim,real>::get_residual_l2norm ()
 {
@@ -237,7 +433,7 @@ void DGBase<dim,real>::evaluate_inverse_mass_matrices ()
     // Invert and store mass matrix
     // Using Gauss-Jordan since it's deal.II's default invert function
     // Could store Cholesky decomposition for more efficient pre-processing
-    const int n_quad_pts      = quadrature.size();
+    const int n_quad_pts      = volume_quadrature.size();
     // Number of test functions, not total number of degrees of freedom
     const int n_tests_cell    = fe_system.dofs_per_cell;
 
@@ -248,19 +444,20 @@ void DGBase<dim,real>::evaluate_inverse_mass_matrices ()
     inv_mass_matrix.resize(triangulation->n_active_cells(),
                            dealii::FullMatrix<real>(n_tests_cell));
     dealii::FullMatrix<real> mass_matrix(n_tests_cell);
+    dealii::FEValues<dim,dim> fe_values_cell (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::volume_quadrature, this->update_flags);
     for (; cell!=endc; ++cell) {
 
         int cell_index = cell->index();
-        fe_values_cell->reinit(cell);
+        fe_values_cell.reinit(cell);
 
         for (int itest=0; itest<n_tests_cell; ++itest) {
             for (int itrial=itest; itrial<n_tests_cell; ++itrial) {
                 mass_matrix[itest][itrial] = 0.0;
                 for (int iquad=0; iquad<n_quad_pts; ++iquad) {
                     mass_matrix[itest][itrial] +=
-                        fe_values_cell->shape_value(itest,iquad)
-                        * fe_values_cell->shape_value(itrial,iquad)
-                        * fe_values_cell->JxW(iquad);
+                        fe_values_cell.shape_value(itest,iquad)
+                        * fe_values_cell.shape_value(itrial,iquad)
+                        * fe_values_cell.JxW(iquad);
                 }
                 mass_matrix[itrial][itest] = mass_matrix[itest][itrial];
             }
@@ -280,7 +477,7 @@ void DGBase<dim,real>::evaluate_mass_matrices ()
 
     global_mass_matrix.reinit(sp);
 
-    const int n_quad_pts      = quadrature.size();
+    const int n_quad_pts      = volume_quadrature.size();
 
     typename dealii::DoFHandler<dim>::active_cell_iterator
        cell = dof_handler.begin_active(),
@@ -288,11 +485,12 @@ void DGBase<dim,real>::evaluate_mass_matrices ()
 
     dealii::FullMatrix<real> local_mass_matrix(n_dofs_per_cell);
     std::vector<dealii::types::global_dof_index> dofs_indices (n_dofs_per_cell);
+    dealii::FEValues<dim,dim> fe_values_cell (DGBase<dim,real>::mapping, DGBase<dim,real>::fe_system, DGBase<dim,real>::volume_quadrature, this->update_flags);
     for (; cell!=endc; ++cell) {
 
         //int cell_index = cell->index();
         cell->get_dof_indices (dofs_indices);
-        fe_values_cell->reinit(cell);
+        fe_values_cell.reinit(cell);
 
 
         //const int n_dofs_per_state = fe_dg.dofs_per_cell;
@@ -302,9 +500,9 @@ void DGBase<dim,real>::evaluate_mass_matrices ()
         //            real value = 0.0;
         //            for (int iquad=0; iquad<n_quad_pts; ++iquad) {
         //                value +=
-        //                    fe_values_cell->shape_value_component(itest,iquad,istate)
-        //                    * fe_values_cell->shape_value_component(itrial,iquad,istate)
-        //                    * fe_values_cell->JxW(iquad);
+        //                    fe_values_cell.shape_value_component(itest,iquad,istate)
+        //                    * fe_values_cell.shape_value_component(itrial,iquad,istate)
+        //                    * fe_values_cell.JxW(iquad);
         //            }
         //            local_mass_matrix[istate*n_dofs_per_state+itrial][istate*n_dofs_per_state+itest] = value;
         //            local_mass_matrix[istate*n_dofs_per_state+itest][istate*n_dofs_per_state+itrial] = value;
@@ -313,15 +511,15 @@ void DGBase<dim,real>::evaluate_mass_matrices ()
         //}
 
         for (int itest=0; itest<n_dofs_per_cell; ++itest) {
-            const unsigned int istate_test = fe_values_cell->get_fe().system_to_component_index(itest).first;
+            const unsigned int istate_test = fe_values_cell.get_fe().system_to_component_index(itest).first;
             for (int itrial=itest; itrial<n_dofs_per_cell; ++itrial) {
-                const unsigned int istate_trial = fe_values_cell->get_fe().system_to_component_index(itrial).first;
+                const unsigned int istate_trial = fe_values_cell.get_fe().system_to_component_index(itrial).first;
                 real value = 0.0;
                 for (int iquad=0; iquad<n_quad_pts; ++iquad) {
                     value +=
-                        fe_values_cell->shape_value_component(itest,iquad,istate_test)
-                        * fe_values_cell->shape_value_component(itrial,iquad,istate_trial)
-                        * fe_values_cell->JxW(iquad);
+                        fe_values_cell.shape_value_component(itest,iquad,istate_test)
+                        * fe_values_cell.shape_value_component(itrial,iquad,istate_trial)
+                        * fe_values_cell.JxW(iquad);
                 }
                 local_mass_matrix[itrial][itest] = 0.0;
                 local_mass_matrix[itest][itrial] = 0.0;
@@ -366,249 +564,7 @@ std::vector<real> DGBase<dim,real>::evaluate_time_steps (const bool exact_time_s
     return time_steps;
 }
 
-
-// DG     ***************************************************************************
-// Constructor
-template <int dim, int nstate, typename real>
-DG<dim,nstate,real>::DG(
-    const Parameters::AllParameters *const parameters_input,
-    const unsigned int degree)
-    : DGBase<dim,real>::DGBase(nstate, parameters_input, degree) // Use DGBase constructor
-
-{
-    using ADtype = Sacado::Fad::DFad<real>;
-    pde_physics = Physics::PhysicsFactory<dim, nstate, ADtype >
-        ::create_Physics(parameters_input->pde_type);
-    conv_num_flux = NumericalFlux::NumericalFluxFactory<dim, nstate, ADtype>
-        ::create_convective_numerical_flux (parameters_input->conv_num_flux_type, pde_physics);
-    diss_num_flux = NumericalFlux::NumericalFluxFactory<dim, nstate, ADtype>
-        ::create_dissipative_numerical_flux (parameters_input->diss_num_flux_type, pde_physics);
-}
-// Destructor
-template <int dim, int nstate, typename real>
-DG<dim,nstate,real>::~DG ()
-{ 
-    std::cout << "Destructing DG..." << std::endl;
-    delete conv_num_flux;
-    delete diss_num_flux;
-    delete pde_physics;
-}
-
-template <int dim, int nstate, typename real>
-void DG<dim,nstate,real>::allocate_system ()
-{
-    std::cout << std::endl << "Allocating DG system and initializing FEValues" << std::endl;
-    // This function allocates all the necessary memory to the 
-    // system matrices and vectors.
-
-    DGBase<dim,real>::dof_handler.initialize(*DGBase<dim,real>::triangulation, DGBase<dim,real>::fe_system);
-    //DGBase<dim,real>::dof_handler.initialize(*DGBase<dim,real>::triangulation, DGBase<dim,real>::fe_system);
-    // Allocates memory from triangulation and finite element space
-    // Use fe_system since it will have the (fe_system.n_dofs)*nstate
-    DGBase<dim,real>::dof_handler.distribute_dofs(DGBase<dim,real>::fe_system);
-
-    //std::vector<unsigned int> block_component(nstate,0);
-    //dealii::DoFRenumbering::component_wise(DGBase<dim,real>::dof_handler, block_component);
-
-    // Allocate matrix
-    unsigned int n_dofs = DGBase<dim,real>::dof_handler.n_dofs();
-    //DynamicSparsityPattern dsp(n_dofs, n_dofs);
-    DGBase<dim,real>::sparsity_pattern.reinit(n_dofs, n_dofs);
-
-    dealii::DoFTools::make_flux_sparsity_pattern(DGBase<dim,real>::dof_handler, DGBase<dim,real>::sparsity_pattern);
-
-    DGBase<dim,real>::system_matrix.reinit(DGBase<dim,real>::sparsity_pattern);
-
-    // Allocate vectors
-    DGBase<dim,real>::solution.reinit(n_dofs);
-    DGBase<dim,real>::right_hand_side.reinit(n_dofs);
-
-}
-
-template <int dim, int nstate, typename real>
-void DG<dim,nstate,real>::assemble_system ()
-{
-    DGBase<dim,real>::system_matrix = 0;
-    DGBase<dim,real>::right_hand_side = 0;
-
-    // For now assume same polynomial degree across domain
-    const unsigned int dofs_per_cell = DGBase<dim,real>::dof_handler.get_fe().dofs_per_cell;
-    std::vector<dealii::types::global_dof_index> current_dofs_indices (dofs_per_cell);
-    std::vector<dealii::types::global_dof_index> neighbor_dofs_indices (dofs_per_cell);
-
-    // Local vector contribution from each cell
-    dealii::Vector<double> current_cell_rhs (dofs_per_cell);
-    dealii::Vector<double> neighbor_cell_rhs (dofs_per_cell);
-
-    // ACTIVE cells, therefore, no children
-    typename dealii::DoFHandler<dim>::active_cell_iterator
-        current_cell = DGBase<dim,real>::dof_handler.begin_active(),
-        endc = DGBase<dim,real>::dof_handler.end();
-
-    unsigned int n_cell_visited = 0;
-    unsigned int n_face_visited = 0;
-    for (; current_cell!=endc; ++current_cell) {
-        // std::cout << "Current cell index: " << current_cell->index() << std::endl;
-        n_cell_visited++;
-
-        current_cell_rhs = 0;
-        DGBase<dim,real>::fe_values_cell->reinit (current_cell);
-        current_cell->get_dof_indices (current_dofs_indices);
-
-        assemble_cell_terms_implicit (DGBase<dim,real>::fe_values_cell, current_dofs_indices, current_cell_rhs);
-
-        for (unsigned int iface=0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface) {
-
-            typename dealii::DoFHandler<dim>::face_iterator current_face = current_cell->face(iface);
-            typename dealii::DoFHandler<dim>::cell_iterator neighbor_cell = current_cell->neighbor(iface);
-
-            // See tutorial step-30 for breakdown of 4 face cases
-
-            // Case 1:
-            // Face at boundary
-            if (current_face->at_boundary()) {
-
-                n_face_visited++;
-
-                DGBase<dim,real>::fe_values_face_int->reinit (current_cell, iface);
-                const unsigned int degree_current = DGBase<dim,real>::fe_system.tensor_degree();
-                const unsigned int deg1sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
-                const unsigned int normal_direction = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
-                const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction);
-
-                real penalty = deg1sq / vol_div_facearea1;
-                //penalty = 1;//99;
-
-                // Need to somehow get boundary type from the mesh
-                assemble_boundary_term_implicit (DGBase<dim,real>::fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs);
-
-            // Case 2:
-            // Neighbour is finer occurs if the face has children
-            // This is because we are looping over the current_cell's face, so 2, 4, and 6 faces.
-            } else if (current_face->has_children()) {
-                std::cout << "SHOULD NOT HAPPEN!!!!!!!!!!!! I haven't put in adaptatation yet" << std::endl;
-                neighbor_cell_rhs = 0;
-                Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError());
-
-                // Obtain cell neighbour
-                const unsigned int neighbor_face_no = current_cell->neighbor_face_no(iface);
-
-                for (unsigned int subface_no=0; subface_no < current_face->number_of_children(); ++subface_no) {
-
-                    n_face_visited++;
-
-                    typename dealii::DoFHandler<dim>::cell_iterator neighbor_child_cell = current_cell->neighbor_child_on_subface (iface, subface_no);
-
-                    Assert (!neighbor_child_cell->has_children(), dealii::ExcInternalError());
-
-                    neighbor_child_cell->get_dof_indices (neighbor_dofs_indices);
-
-                    DGBase<dim,real>::fe_values_subface_int->reinit (current_cell, iface, subface_no);
-                    DGBase<dim,real>::fe_values_face_ext->reinit (neighbor_child_cell, neighbor_face_no);
-
-                    const unsigned int normal_direction1 = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
-                    const unsigned int normal_direction2 = dealii::GeometryInfo<dim>::unit_normal_direction[neighbor_face_no];
-                    const unsigned int degree_current = DGBase<dim,real>::fe_system.tensor_degree();
-                    const unsigned int deg1sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
-                    const unsigned int deg2sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
-
-                    //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children();
-                    const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1);
-                    const real vol_div_facearea2 = neighbor_child_cell->extent_in_direction(normal_direction2);
-
-                    const real penalty1 = deg1sq / vol_div_facearea1;
-                    const real penalty2 = deg2sq / vol_div_facearea2;
-                    
-                    real penalty = 0.5 * ( penalty1 + penalty2 );
-                    //penalty = 1;
-
-                    assemble_face_term_implicit (
-                        DGBase<dim,real>::fe_values_subface_int, DGBase<dim,real>::fe_values_face_ext,
-                        penalty,
-                        current_dofs_indices, neighbor_dofs_indices,
-                        current_cell_rhs, neighbor_cell_rhs);
-
-                    // Add local contribution from neighbor cell to global vector
-                    for (unsigned int i=0; i<dofs_per_cell; ++i) {
-                        DGBase<dim,real>::right_hand_side(neighbor_dofs_indices[i]) += neighbor_cell_rhs(i);
-                    }
-                }
-
-            // Case 3:
-            // Neighbor cell is NOT coarser
-            // Therefore, they have the same coarseness, and we need to choose one of them to do the work
-            } else if (
-                !current_cell->neighbor_is_coarser(iface) &&
-                    // Cell with lower index does work
-                    (neighbor_cell->index() > current_cell->index() || 
-                    // If both cells have same index
-                    // See https://www.dealii.org/developer/doxygen/deal.II/classTriaAccessorBase.html#a695efcbe84fefef3e4c93ee7bdb446ad
-                    // then cell at the lower level does the work
-                        (neighbor_cell->index() == current_cell->index() && current_cell->level() < neighbor_cell->level())
-                    ) )
-            {
-                n_face_visited++;
-
-                neighbor_cell_rhs = 0;
-                Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError());
-                typename dealii::DoFHandler<dim>::cell_iterator neighbor_cell = current_cell->neighbor(iface);
-
-                neighbor_cell->get_dof_indices (neighbor_dofs_indices);
-
-                const unsigned int neighbor_face_no = current_cell->neighbor_of_neighbor(iface);
-
-                const unsigned int normal_direction1 = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
-                const unsigned int normal_direction2 = dealii::GeometryInfo<dim>::unit_normal_direction[neighbor_face_no];
-                const unsigned int degree_current = DGBase<dim,real>::fe_system.tensor_degree();
-                const unsigned int deg1sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
-                const unsigned int deg2sq = (degree_current == 0) ? 1 : degree_current * (degree_current+1);
-
-                //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children();
-                const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1);
-                const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2);
-
-                const real penalty1 = deg1sq / vol_div_facearea1;
-                const real penalty2 = deg2sq / vol_div_facearea2;
-                
-                real penalty = 0.5 * ( penalty1 + penalty2 );
-                //penalty = 1;//99;
-
-                DGBase<dim,real>::fe_values_face_int->reinit (current_cell, iface);
-                DGBase<dim,real>::fe_values_face_ext->reinit (neighbor_cell, neighbor_face_no);
-                assemble_face_term_implicit (
-                        DGBase<dim,real>::fe_values_face_int, DGBase<dim,real>::fe_values_face_ext,
-                        penalty,
-                        current_dofs_indices, neighbor_dofs_indices,
-                        current_cell_rhs, neighbor_cell_rhs);
-
-                // Add local contribution from neighbor cell to global vector
-                for (unsigned int i=0; i<dofs_per_cell; ++i) {
-                    DGBase<dim,real>::right_hand_side(neighbor_dofs_indices[i]) += neighbor_cell_rhs(i);
-                }
-            } else {
-            // Do nothing
-            }
-            // Case 4: Neighbor is coarser
-            // Do nothing.
-            // The face contribution from the current cell will appear then the coarse neighbor checks for subfaces
-
-        } // end of face loop
-
-        for (unsigned int i=0; i<dofs_per_cell; ++i) {
-            DGBase<dim,real>::right_hand_side(current_dofs_indices[i]) += current_cell_rhs(i);
-        }
-
-    } // end of cell loop
-
-} // end of assemble_system_implicit ()
-
-
 template class DGBase <PHILIP_DIM, double>;
 template class DGFactory <PHILIP_DIM, double>;
-template class DG <PHILIP_DIM, 1, double>;
-template class DG <PHILIP_DIM, 2, double>;
-template class DG <PHILIP_DIM, 3, double>;
-template class DG <PHILIP_DIM, 4, double>;
-template class DG <PHILIP_DIM, 5, double>;
 
 } // PHiLiP namespace
diff --git a/src/dg/dg.h b/src/dg/dg.h
index 111cb830f..994fc6ddd 100644
--- a/src/dg/dg.h
+++ b/src/dg/dg.h
@@ -32,7 +32,7 @@ namespace PHiLiP {
 /**  This base class allows the use of arrays to efficiently allocate the data structures
   *  through std::array in the derived class DG.
   *  This class is the one being returned by the DGFactory and is the main
-  *  interface for a user to call its main functions such as "assemble_system".
+  *  interface for a user to call its main functions such as "assemble_residual".
   *
   *  Discretizes the problem
   *  \f[
@@ -42,9 +42,8 @@ namespace PHiLiP {
   *      + \mathbf{F}_{diss}(\mathbf{u},\boldsymbol\nabla\mathbf{u}) )
   *      = \mathbf{q}
   *  \f]
-  *
-  *  Note that it does not currently handle nonlinear or time-accurate problems
-  *  since mass matrices have not been used/verified yet.
+  *  
+  *  Also defines the main loop of the DGWeak class which is assemble_residual
   */
 template <int dim, typename real>
 class DGBase 
@@ -59,7 +58,7 @@ class DGBase
     /// Constructor. Deleted the default constructor since it should not be used
     DGBase () = delete;
     /// Principal constructor.
-    /** Will initialize mapping, fe_dg, all_parameters, quadrature, and face_quadrature
+    /** Will initialize mapping, fe_dg, all_parameters, volume_quadrature, and face_quadrature
      *  from DGBase. The it will new some FEValues that will be used to retrieve the
      *  finite element values at physical locations.
      */
@@ -68,15 +67,15 @@ class DGBase
         const Parameters::AllParameters *const parameters_input, 
         const unsigned int degree);
 
-    virtual ~DGBase(); ///< Destructor. Also calls delete_fe_values().
-    void delete_fe_values (); ///< Delete the new-ed pointers to FEValues
+    virtual ~DGBase(); ///< Destructor.
 
     /// Sets the triangulation. Should be done before allocate system
     void set_triangulation(dealii::Triangulation<dim> *triangulation_input)
     { triangulation = triangulation_input; } ;
 
-    virtual void allocate_system () = 0; ///< Virtual function defined in DG
-    virtual void assemble_system () = 0; ///< Virtual function defined in DG
+    /// Allocates the system.
+    /** Must be done after setting the mesh and before assembling the system. */
+    virtual void allocate_system ();
 
     /// Allocates and evaluates the inverse mass matrices for the entire grid
     /*  Although straightforward, this has not been tested yet.
@@ -196,19 +195,79 @@ class DGBase
      */ 
     dealii::DoFHandler<dim> dof_handler;
 
+    /// Main loop of the DG class.
+    /** Evaluates the right-hand-side \f$ \mathbf{R(\mathbf{u}}) \f$ of the system
+     *
+     *  \f[
+     *      \frac{\partial \mathbf{u}}{\partial t} = \mathbf{R(\mathbf{u}}) = 
+     *      - \boldsymbol\nabla \cdot
+     *      ( \mathbf{F}_{conv}(\mathbf{u})
+     *      + \mathbf{F}_{diss}(\mathbf{u},\boldsymbol\nabla\mathbf{u}) )
+     *      + \mathbf{q}
+     *  \f]
+     *
+     *  As well as sets the
+     *  \f[
+     *  \mathbf{\text{system_matrix}} = \frac{\partial \mathbf{R}}{\partial \mathbf{u}}
+     *  \f]
+     *
+     * It loops over all the cells, evaluates the volume contributions,
+     * then loops over the faces of the current cell. Four scenarios may happen
+     *
+     * 1. Boundary condition.
+     *
+     * 2. Current face has children. Therefore, neighbor is finer. In that case,
+     * loop over neighbor faces to compute its face contributions.
+     *
+     * 3. Neighbor has same coarseness. Cell with lower global index will be used
+     * to compute the face contribution.
+     *
+     * 4. Neighbor is coarser. Therefore, the current cell is the finer one.
+     * Do nothing since this cell will be taken care of by scenario 2.
+     *    
+     */
+    void assemble_residual_dRdW ();
 
 protected:
 
+    /// Evaluate the integral over the cell volume
+    virtual void assemble_cell_terms_implicit(
+        const dealii::FEValues<dim,dim> &fe_values_cell,
+        const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
+        dealii::Vector<real> &current_cell_rhs) = 0;
+    /// Evaluate the integral over the cell edges that are on domain boundaries
+    virtual void assemble_boundary_term_implicit(
+        const dealii::FEFaceValues<dim,dim> &fe_values_face_int,
+        const real penalty,
+        const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
+        dealii::Vector<real> &current_cell_rhs) = 0;
+    /// Evaluate the integral over the internal cell edges
+    virtual void assemble_face_term_implicit(
+        const dealii::FEValuesBase<dim,dim>     &fe_values_face_int,
+        const dealii::FEFaceValues<dim,dim>     &fe_values_face_ext,
+        const real penalty,
+        const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
+        const std::vector<dealii::types::global_dof_index> &neighbor_dofs_indices,
+        dealii::Vector<real>          &current_cell_rhs,
+        dealii::Vector<real>          &neighbor_cell_rhs) = 0;
+
     // QGauss is Gauss-Legendre quadrature nodes
-    const dealii::QGauss<dim>   quadrature;
+    const dealii::QGauss<1>     oned_quadrature; // For the strong form
+    const dealii::QGauss<dim>   volume_quadrature;
     const dealii::QGauss<dim-1> face_quadrature;
-    // const dealii::QGaussLobatto<dim>   quadrature;
+    // const dealii::QGaussLobatto<dim>   volume_quadrature;
     // const dealii::QGaussLobatto<dim-1> face_quadrature;
 
-    dealii::FEValues<dim,dim>         *fe_values_cell;
-    dealii::FEFaceValues<dim,dim>     *fe_values_face_int;
-    dealii::FESubfaceValues<dim,dim>  *fe_values_subface_int;
-    dealii::FEFaceValues<dim,dim>     *fe_values_face_ext;
+    const dealii::UpdateFlags update_flags =
+        dealii::update_values | dealii::update_gradients
+        | dealii::update_quadrature_points | dealii::update_JxW_values;
+    const dealii::UpdateFlags face_update_flags =
+        dealii::update_values | dealii::update_gradients
+        | dealii::update_quadrature_points | dealii::update_JxW_values
+        | dealii::update_normal_vectors;
+    const dealii::UpdateFlags neighbor_face_update_flags =
+        dealii::update_values | dealii::update_gradients;
+
     /// Main loop of the DGBase class.
     /** It loops over all the cells, evaluates the volume contributions,
      * then loops over the faces of the current cell. Four scenarios may happen
@@ -228,96 +287,100 @@ class DGBase
     //virtual void allocate_system_implicit () = 0;
     //virtual void assemble_system_implicit () = 0;
 
-}; // end of DG class
+}; // end of DGBase class
 
-/// DG class templated on the number of state variables
+/// DGWeak class templated on the number of state variables
 /*  Contains the functions that need to be templated on the number of state variables.
- *  
- *  Also defines the main loop of the DG class which is assemble_system
  */
 template <int dim, int nstate, typename real>
-class DG : public DGBase<dim, real>
+class DGWeak : public DGBase<dim, real>
 {
 public:
     /// Constructor
-    DG(
+    DGWeak(
         const Parameters::AllParameters *const parameters_input, 
         const unsigned int degree);
 
     /// Destructor
-    ~DG();
+    ~DGWeak();
 
 private:
     /// Contains the physics of the PDE
-    Physics::PhysicsBase<dim, nstate, Sacado::Fad::DFad<real> >  *pde_physics;
+    std::shared_ptr < Physics::PhysicsBase<dim, nstate, Sacado::Fad::DFad<real> > > pde_physics;
     /// Convective numerical flux
     NumericalFlux::NumericalFluxConvective<dim, nstate, Sacado::Fad::DFad<real> > *conv_num_flux;
     /// Dissipative numerical flux
     NumericalFlux::NumericalFluxDissipative<dim, nstate, Sacado::Fad::DFad<real> > *diss_num_flux;
 
+    /// Evaluate the integral over the cell volume
+    void assemble_cell_terms_implicit(
+        const dealii::FEValues<dim,dim> &fe_values_cell,
+        const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
+        dealii::Vector<real> &current_cell_rhs);
+    /// Evaluate the integral over the cell edges that are on domain boundaries
+    void assemble_boundary_term_implicit(
+        const dealii::FEFaceValues<dim,dim> &fe_values_face_int,
+        const real penalty,
+        const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
+        dealii::Vector<real> &current_cell_rhs);
+    /// Evaluate the integral over the internal cell edges
+    void assemble_face_term_implicit(
+        const dealii::FEValuesBase<dim,dim>     &fe_values_face_int,
+        const dealii::FEFaceValues<dim,dim>     &fe_values_face_ext,
+        const real penalty,
+        const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
+        const std::vector<dealii::types::global_dof_index> &neighbor_dofs_indices,
+        dealii::Vector<real>          &current_cell_rhs,
+        dealii::Vector<real>          &neighbor_cell_rhs);
+
+}; // end of DGWeak class
 
+/// DGStrong class templated on the number of state variables
+/*  Contains the functions that need to be templated on the number of state variables.
+ */
+template <int dim, int nstate, typename real>
+class DGStrong : public DGBase<dim, real>
+{
+public:
+    /// Constructor
+    DGStrong(
+        const Parameters::AllParameters *const parameters_input, 
+        const unsigned int degree);
 
-    /// Allocates the system.
-    /** Must be done after setting the mesh and before assembling the system.
-     */
-    void allocate_system ();
+    /// Destructor
+    ~DGStrong();
 
-    /// Main loop of the DG class.
-    /** Evaluates the right-hand-side \f$ \mathbf{R(\mathbf{u}}) \f$ of the system
-     *
-     *  \f[
-     *      \frac{\partial \mathbf{u}}{\partial t} = \mathbf{R(\mathbf{u}}) = 
-     *      - \boldsymbol\nabla \cdot
-     *      ( \mathbf{F}_{conv}(\mathbf{u})
-     *      + \mathbf{F}_{diss}(\mathbf{u},\boldsymbol\nabla\mathbf{u}) )
-     *      + \mathbf{q}
-     *  \f]
-     *
-     *  As well as sets the
-     *  \f[
-     *  \mathbf{\text{system_matrix}} = \frac{\partial \mathbf{R}}{\partial \mathbf{u}}
-     *  \f]
-     *
-     * It loops over all the cells, evaluates the volume contributions,
-     * then loops over the faces of the current cell. Four scenarios may happen
-     *
-     * 1. Boundary condition.
-     *
-     * 2. Current face has children. Therefore, neighbor is finer. In that case,
-     * loop over neighbor faces to compute its face contributions.
-     *
-     * 3. Neighbor has same coarseness. Cell with lower global index will be used
-     * to compute the face contribution.
-     *
-     * 4. Neighbor is coarser. Therefore, the current cell is the finer one.
-     * Do nothing since this cell will be taken care of by scenario 2.
-     *    
-     */
-    void assemble_system ();
+private:
+    /// Contains the physics of the PDE
+    std::shared_ptr < Physics::PhysicsBase<dim, nstate, Sacado::Fad::DFad<real> > > pde_physics;
+    /// Convective numerical flux
+    NumericalFlux::NumericalFluxConvective<dim, nstate, Sacado::Fad::DFad<real> > *conv_num_flux;
+    /// Dissipative numerical flux
+    NumericalFlux::NumericalFluxDissipative<dim, nstate, Sacado::Fad::DFad<real> > *diss_num_flux;
 
 
     /// Evaluate the integral over the cell volume
     void assemble_cell_terms_implicit(
-        const dealii::FEValues<dim,dim> *fe_values_cell,
+        const dealii::FEValues<dim,dim> &fe_values_cell,
         const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
         dealii::Vector<real> &current_cell_rhs);
     /// Evaluate the integral over the cell edges that are on domain boundaries
     void assemble_boundary_term_implicit(
-        const dealii::FEFaceValues<dim,dim> *fe_values_face_int,
+        const dealii::FEFaceValues<dim,dim> &fe_values_face_int,
         const real penalty,
         const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
         dealii::Vector<real> &current_cell_rhs);
     /// Evaluate the integral over the internal cell edges
     void assemble_face_term_implicit(
-        const dealii::FEValuesBase<dim,dim>     *fe_values_face_int,
-        const dealii::FEFaceValues<dim,dim>     *fe_values_face_ext,
+        const dealii::FEValuesBase<dim,dim>     &fe_values_face_int,
+        const dealii::FEFaceValues<dim,dim>     &fe_values_face_ext,
         const real penalty,
         const std::vector<dealii::types::global_dof_index> &current_dofs_indices,
         const std::vector<dealii::types::global_dof_index> &neighbor_dofs_indices,
         dealii::Vector<real>          &current_cell_rhs,
         dealii::Vector<real>          &neighbor_cell_rhs);
 
-}; // end of DG class
+}; // end of DGStrong class
 
 /// This class creates a new DGBase object
 /** This allows the DGBase to not be templated on the number of state variables
diff --git a/src/dg/strong_dg.cpp b/src/dg/strong_dg.cpp
new file mode 100644
index 000000000..8d452f5c6
--- /dev/null
+++ b/src/dg/strong_dg.cpp
@@ -0,0 +1,460 @@
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <Sacado.hpp>
+//#include <deal.II/differentiation/ad/sacado_math.h>
+//#include <deal.II/differentiation/ad/sacado_number_types.h>
+#include <deal.II/differentiation/ad/sacado_product_types.h>
+
+#include <deal.II/fe/fe_dgq.h> // Used for flux interpolation
+
+#include "dg.h"
+
+namespace PHiLiP {
+
+
+template <int dim, int nstate, typename real>
+DGStrong<dim,nstate,real>::DGStrong(
+    const Parameters::AllParameters *const parameters_input,
+    const unsigned int degree)
+    : DGBase<dim,real>::DGBase(nstate, parameters_input, degree) // Use DGBase constructor
+{
+    using ADtype = Sacado::Fad::DFad<real>;
+    pde_physics = Physics::PhysicsFactory<dim,nstate,ADtype> 
+        ::create_Physics(parameters_input->pde_type);
+    conv_num_flux = NumericalFlux::NumericalFluxFactory<dim, nstate, ADtype>
+        ::create_convective_numerical_flux (parameters_input->conv_num_flux_type, pde_physics);
+    diss_num_flux = NumericalFlux::NumericalFluxFactory<dim, nstate, ADtype>
+        ::create_dissipative_numerical_flux (parameters_input->diss_num_flux_type, pde_physics);
+}
+
+template <int dim, int nstate, typename real>
+DGStrong<dim,nstate,real>::~DGStrong ()
+{ 
+    std::cout << "Destructing DGStrong..." << std::endl;
+    delete conv_num_flux;
+    delete diss_num_flux;
+}
+
+template <int dim, int nstate, typename real>
+void DGStrong<dim,nstate,real>::assemble_cell_terms_implicit(
+    const dealii::FEValues<dim,dim> &fe_values_vol,
+    const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
+    dealii::Vector<real> &local_rhs_int_cell)
+{
+    using ADtype = Sacado::Fad::DFad<real>;
+    using ADArray = std::array<ADtype,nstate>;
+    using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
+
+    const unsigned int n_quad_pts      = fe_values_vol.n_quadrature_points;
+    const unsigned int n_dofs_cell     = fe_values_vol.dofs_per_cell;
+
+    AssertDimension (n_dofs_cell, cell_dofs_indices.size());
+
+    const std::vector<real> &JxW = fe_values_vol.get_JxW_values ();
+
+
+    std::vector<real> residual_derivatives(n_dofs_cell);
+
+    std::vector< ADArray > soln_at_q(n_quad_pts);
+    std::vector< ADArrayTensor1 > soln_grad_at_q(n_quad_pts); // Tensor initialize with zeros
+
+    std::vector< ADArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
+    std::vector< ADArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
+    std::vector< ADArray > source_at_q(n_quad_pts);
+
+
+    // AD variable
+    std::vector< ADtype > soln_coeff(n_dofs_cell);
+    for (unsigned int idof = 0; idof < n_dofs_cell; ++idof) {
+        soln_coeff[idof] = DGBase<dim,real>::solution(cell_dofs_indices[idof]);
+        soln_coeff[idof].diff(idof, n_dofs_cell);
+    }
+    for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
+        for (int istate=0; istate<nstate; istate++) { 
+            // Interpolate solution to the volume quadrature points
+            soln_at_q[iquad][istate]      = 0;
+            soln_grad_at_q[iquad][istate] = 0;
+        }
+    }
+    // Interpolate solution to face
+    for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
+        for (unsigned int idof=0; idof<n_dofs_cell; ++idof) {
+              const unsigned int istate = fe_values_vol.get_fe().system_to_component_index(idof).first;
+              soln_at_q[iquad][istate]      += soln_coeff[idof] * fe_values_vol.shape_value_component(idof, iquad, istate);
+              soln_grad_at_q[iquad][istate] += soln_coeff[idof] * fe_values_vol.shape_grad_component(idof, iquad, istate);
+        }
+        //std::cout << "Density " << soln_at_q[iquad][0] << std::endl;
+        //if(nstate>1) std::cout << "Momentum " << soln_at_q[iquad][1] << std::endl;
+        //std::cout << "Energy " << soln_at_q[iquad][nstate-1] << std::endl;
+        // Evaluate physical convective flux and source term
+        conv_phys_flux_at_q[iquad] = pde_physics->convective_flux (soln_at_q[iquad]);
+        diss_phys_flux_at_q[iquad] = pde_physics->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad]);
+        source_at_q[iquad] = pde_physics->source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad]);
+    }
+
+
+    // Evaluate flux divergence by interpolating the flux
+    // Since we have nodal values of the flux, we use the Lagrange polynomials to obtain the gradients at the quadrature points.
+    dealii::FE_DGQArbitraryNodes<dim,dim> lagrange_poly(this->oned_quadrature);
+    dealii::FEValues<dim,dim> fe_values_lagrange (this->mapping, lagrange_poly, this->volume_quadrature, this->update_flags);
+    fe_values_lagrange.reinit(fe_values_vol.get_cell());
+    std::vector<ADArray> flux_divergence(n_quad_pts);
+    for (int istate = 0; istate<nstate; ++istate) {
+        for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
+            flux_divergence[iquad][istate] = 0.0;
+            for ( unsigned int flux_basis = 0; flux_basis < n_quad_pts; ++flux_basis ) {
+                flux_divergence[iquad][istate] += conv_phys_flux_at_q[flux_basis][istate] * fe_values_lagrange.shape_grad(flux_basis,iquad);
+            }
+        }
+    }
+
+    // Strong form
+    // The right-hand side sends all the term to the side of the source term
+    // Therefore, 
+    // \divergence ( Fconv + Fdiss ) = source 
+    // has the right-hand side
+    // rhs = - \divergence( Fconv + Fdiss ) + source 
+    // Since we have done an integration by parts, the volume term resulting from the divergence of Fconv and Fdiss
+    // is negative. Therefore, negative of negative means we add that volume term to the right-hand-side
+    for (unsigned int itest=0; itest<n_dofs_cell; ++itest) {
+
+        ADtype rhs = 0;
+
+        const unsigned int istate = fe_values_vol.get_fe().system_to_component_index(itest).first;
+
+        for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
+
+            // 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
+            //// Note that for diffusion, the negative is defined in the physics
+            rhs = rhs + fe_values_vol.shape_grad_component(itest,iquad,istate) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
+            // Source
+            rhs = rhs + fe_values_vol.shape_value_component(itest,iquad,istate) * source_at_q[iquad][istate] * JxW[iquad];
+        }
+
+        local_rhs_int_cell(itest) += rhs.val();
+
+        for (unsigned int idof = 0; idof < n_dofs_cell; ++idof) {
+            //residual_derivatives[idof] = rhs.fastAccessDx(idof);
+            residual_derivatives[idof] = rhs.dx(idof);
+        }
+        this->system_matrix.add(cell_dofs_indices[itest], cell_dofs_indices, residual_derivatives);
+    }
+}
+
+
+template <int dim, int nstate, typename real>
+void DGStrong<dim,nstate,real>::assemble_boundary_term_implicit(
+    const dealii::FEFaceValues<dim,dim> &fe_values_boundary,
+    const real penalty,
+    const std::vector<dealii::types::global_dof_index> &dof_indices_int,
+    dealii::Vector<real> &local_rhs_int_cell)
+{
+    using ADtype = Sacado::Fad::DFad<real>;
+    using ADArray = std::array<ADtype,nstate>;
+    using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
+
+    const unsigned int n_dofs_cell = fe_values_boundary.dofs_per_cell;
+    const unsigned int n_face_quad_pts = fe_values_boundary.n_quadrature_points;
+
+    AssertDimension (n_dofs_cell, dof_indices_int.size());
+
+    const std::vector<real> &JxW = fe_values_boundary.get_JxW_values ();
+    const std::vector<dealii::Tensor<1,dim>> &normals = fe_values_boundary.get_normal_vectors ();
+
+    std::vector<real> residual_derivatives(n_dofs_cell);
+
+    std::vector<ADArray> soln_int(n_face_quad_pts);
+    std::vector<ADArray> soln_ext(n_face_quad_pts);
+
+    std::vector<ADArrayTensor1> soln_grad_int(n_face_quad_pts);
+    std::vector<ADArrayTensor1> soln_grad_ext(n_face_quad_pts);
+
+    std::vector<ADArray> conv_num_flux_dot_n(n_face_quad_pts);
+    std::vector<ADArray> diss_soln_num_flux(n_face_quad_pts); // u*
+    std::vector<ADArrayTensor1> diss_flux_jump_int(n_face_quad_pts); // u*-u_int
+    std::vector<ADArray> diss_auxi_num_flux_dot_n(n_face_quad_pts); // sigma*
+
+    std::vector<ADArrayTensor1> conv_phys_flux(n_face_quad_pts);
+
+    // AD variable
+    std::vector< ADtype > soln_coeff_int(n_dofs_cell);
+    const unsigned int n_total_indep = n_dofs_cell;
+    for (unsigned int idof = 0; idof < n_dofs_cell; ++idof) {
+        soln_coeff_int[idof] = DGBase<dim,real>::solution(dof_indices_int[idof]);
+        soln_coeff_int[idof].diff(idof, n_total_indep);
+    }
+
+    for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+        for (int istate=0; istate<nstate; istate++) { 
+            // Interpolate solution to the face quadrature points
+            soln_int[iquad][istate]      = 0;
+            soln_grad_int[iquad][istate] = 0;
+        }
+    }
+    // Interpolate solution to face
+    const std::vector< dealii::Point<dim,real> > quad_pts = fe_values_boundary.get_quadrature_points();
+    for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+
+        const dealii::Tensor<1,dim,ADtype> normal_int = normals[iquad];
+        const dealii::Tensor<1,dim,ADtype> normal_ext = -normal_int;
+
+        for (unsigned int idof=0; idof<n_dofs_cell; ++idof) {
+            const int istate = fe_values_boundary.get_fe().system_to_component_index(idof).first;
+            soln_int[iquad][istate]      += soln_coeff_int[idof] * fe_values_boundary.shape_value_component(idof, iquad, istate);
+            soln_grad_int[iquad][istate] += soln_coeff_int[idof] * fe_values_boundary.shape_grad_component(idof, iquad, istate);
+        }
+
+        const unsigned int boundary_id = 1; // dummy value for now
+        const dealii::Point<dim, real> x_quad = quad_pts[iquad];
+        pde_physics->boundary_face_values (boundary_id, x_quad, normal_int, soln_int[iquad], soln_grad_int[iquad], soln_ext[iquad], soln_grad_ext[iquad]);
+
+        //
+        // Evaluate physical convective flux, physical dissipative flux
+        // Following the the boundary treatment given by 
+        //      Hartmann, R., Numerical Analysis of Higher Order Discontinuous Galerkin Finite Element Methods,
+        //      Institute of Aerodynamics and Flow Technology, DLR (German Aerospace Center), 2008.
+        //      Details given on page 93
+        //conv_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_ext[iquad], soln_ext[iquad], normal_int);
+
+        // So, I wasn't able to get Euler manufactured solutions to converge when F* = F*(Ubc, Ubc)
+        // Changing it back to the standdard F* = F*(Uin, Ubc)
+        // This is known not be adjoint consistent as per the paper above. Page 85, second to last paragraph.
+        // Losing 2p+1 OOA on functionals for all PDEs.
+        conv_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+
+        // Used for strong form
+        // Which physical convective flux to use?
+        conv_phys_flux[iquad] = pde_physics->convective_flux (soln_int[iquad]);
+
+        // Notice that the flux uses the solution given by the Dirichlet or Neumann boundary condition
+        diss_soln_num_flux[iquad] = diss_num_flux->evaluate_solution_flux(soln_ext[iquad], soln_ext[iquad], normal_int);
+
+        ADArrayTensor1 diss_soln_jump_int;
+        for (int s=0; s<nstate; s++) {
+            diss_soln_jump_int[s] = (diss_soln_num_flux[iquad][s] - soln_int[iquad][s]) * normal_int;
+        }
+        diss_flux_jump_int[iquad] = pde_physics->dissipative_flux (soln_int[iquad], diss_soln_jump_int);
+
+        diss_auxi_num_flux_dot_n[iquad] = diss_num_flux->evaluate_auxiliary_flux(
+            soln_int[iquad], soln_ext[iquad],
+            soln_grad_int[iquad], soln_grad_ext[iquad],
+            normal_int, penalty, true);
+    }
+
+    // Boundary integral
+    for (unsigned int itest=0; itest<n_dofs_cell; ++itest) {
+
+        ADtype rhs = 0.0;
+
+        const unsigned int istate = fe_values_boundary.get_fe().system_to_component_index(itest).first;
+
+        for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+
+            // Convection
+            const ADtype flux_diff = conv_num_flux_dot_n[iquad][istate] - conv_phys_flux[iquad][istate]*normals[iquad];
+            rhs = rhs - fe_values_boundary.shape_value_component(itest,iquad,istate) * flux_diff * JxW[iquad];
+            // Diffusive
+            rhs = rhs - fe_values_boundary.shape_value_component(itest,iquad,istate) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW[iquad];
+            rhs = rhs + fe_values_boundary.shape_grad_component(itest,iquad,istate) * diss_flux_jump_int[iquad][istate] * JxW[iquad];
+        }
+        // *******************
+
+        local_rhs_int_cell(itest) += rhs.val();
+
+        for (unsigned int idof = 0; idof < n_dofs_cell; ++idof) {
+            //residual_derivatives[idof] = rhs.fastAccessDx(idof);
+            residual_derivatives[idof] = rhs.dx(idof);
+        }
+        this->system_matrix.add(dof_indices_int[itest], dof_indices_int, residual_derivatives);
+    }
+}
+
+template <int dim, int nstate, typename real>
+void DGStrong<dim,nstate,real>::assemble_face_term_implicit(
+    const dealii::FEValuesBase<dim,dim>     &fe_values_int,
+    const dealii::FEFaceValues<dim,dim>     &fe_values_ext,
+    const real penalty,
+    const std::vector<dealii::types::global_dof_index> &dof_indices_int,
+    const std::vector<dealii::types::global_dof_index> &dof_indices_ext,
+    dealii::Vector<real>          &local_rhs_int_cell,
+    dealii::Vector<real>          &local_rhs_ext_cell)
+{
+    using ADtype = Sacado::Fad::DFad<real>;
+    using ADArray = std::array<ADtype,nstate>;
+    using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
+
+    // Use quadrature points of neighbor cell
+    // Might want to use the maximum n_quad_pts1 and n_quad_pts2
+    const unsigned int n_face_quad_pts = fe_values_ext.n_quadrature_points;
+
+    const unsigned int n_dofs_int = fe_values_int.dofs_per_cell;
+    const unsigned int n_dofs_ext = fe_values_ext.dofs_per_cell;
+
+    AssertDimension (n_dofs_int, dof_indices_int.size());
+    AssertDimension (n_dofs_ext, dof_indices_ext.size());
+
+    // Jacobian and normal should always be consistent between two elements
+    // even for non-conforming meshes?
+    const std::vector<real> &JxW_int = fe_values_int.get_JxW_values ();
+    const std::vector<dealii::Tensor<1,dim> > &normals_int = fe_values_int.get_normal_vectors ();
+
+    // AD variable
+    std::vector<ADtype> soln_coeff_int_ad(n_dofs_int);
+    std::vector<ADtype> soln_coeff_ext_ad(n_dofs_ext);
+
+
+    // Jacobian blocks
+    std::vector<real> dR1_dW1(n_dofs_int);
+    std::vector<real> dR1_dW2(n_dofs_ext);
+    std::vector<real> dR2_dW1(n_dofs_int);
+    std::vector<real> dR2_dW2(n_dofs_ext);
+
+    std::vector<ADArray> conv_num_flux_dot_n(n_face_quad_pts);
+    std::vector<ADArrayTensor1> conv_phys_flux_int(n_face_quad_pts);
+    std::vector<ADArrayTensor1> conv_phys_flux_ext(n_face_quad_pts);
+
+    // Interpolate solution to the face quadrature points
+    std::vector< ADArray > soln_int(n_face_quad_pts);
+    std::vector< ADArray > soln_ext(n_face_quad_pts);
+
+    std::vector< ADArrayTensor1 > soln_grad_int(n_face_quad_pts); // Tensor initialize with zeros
+    std::vector< ADArrayTensor1 > soln_grad_ext(n_face_quad_pts); // Tensor initialize with zeros
+
+    std::vector<ADArray> diss_soln_num_flux(n_face_quad_pts); // u*
+    std::vector<ADArray> diss_auxi_num_flux_dot_n(n_face_quad_pts); // sigma*
+
+    std::vector<ADArrayTensor1> diss_flux_jump_int(n_face_quad_pts); // u*-u_int
+    std::vector<ADArrayTensor1> diss_flux_jump_ext(n_face_quad_pts); // u*-u_ext
+    // AD variable
+    const unsigned int n_total_indep = n_dofs_int + n_dofs_ext;
+    for (unsigned int idof = 0; idof < n_dofs_int; ++idof) {
+        soln_coeff_int_ad[idof] = DGBase<dim,real>::solution(dof_indices_int[idof]);
+        soln_coeff_int_ad[idof].diff(idof, n_total_indep);
+    }
+    for (unsigned int idof = 0; idof < n_dofs_ext; ++idof) {
+        soln_coeff_ext_ad[idof] = DGBase<dim,real>::solution(dof_indices_ext[idof]);
+        soln_coeff_ext_ad[idof].diff(idof+n_dofs_int, n_total_indep);
+    }
+    for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+        for (int istate=0; istate<nstate; istate++) { 
+            soln_int[iquad][istate]      = 0;
+            soln_grad_int[iquad][istate] = 0;
+            soln_ext[iquad][istate]      = 0;
+            soln_grad_ext[iquad][istate] = 0;
+        }
+    }
+    for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+
+        const dealii::Tensor<1,dim,ADtype> normal_int = normals_int[iquad];
+        const dealii::Tensor<1,dim,ADtype> normal_ext = -normal_int;
+
+        // Interpolate solution to face
+        for (unsigned int idof=0; idof<n_dofs_int; ++idof) {
+            const unsigned int istate = fe_values_int.get_fe().system_to_component_index(idof).first;
+            soln_int[iquad][istate]      += soln_coeff_int_ad[idof] * fe_values_int.shape_value_component(idof, iquad, istate);
+            soln_grad_int[iquad][istate] += soln_coeff_int_ad[idof] * fe_values_int.shape_grad_component(idof, iquad, istate);
+        }
+        for (unsigned int idof=0; idof<n_dofs_ext; ++idof) {
+            const unsigned int istate = fe_values_ext.get_fe().system_to_component_index(idof).first;
+            soln_ext[iquad][istate]      += soln_coeff_ext_ad[idof] * fe_values_ext.shape_value_component(idof, iquad, istate);
+            soln_grad_ext[iquad][istate] += soln_coeff_ext_ad[idof] * fe_values_ext.shape_grad_component(idof, iquad, istate);
+        }
+        //std::cout << "Density int" << soln_int[iquad][0] << std::endl;
+        //if(nstate>1) std::cout << "Momentum int" << soln_int[iquad][1] << std::endl;
+        //std::cout << "Energy int" << soln_int[iquad][nstate-1] << std::endl;
+        //std::cout << "Density ext" << soln_ext[iquad][0] << std::endl;
+        //if(nstate>1) std::cout << "Momentum ext" << soln_ext[iquad][1] << std::endl;
+        //std::cout << "Energy ext" << soln_ext[iquad][nstate-1] << std::endl;
+
+        // Evaluate physical convective flux, physical dissipative flux, and source term
+        conv_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+
+        conv_phys_flux_int[iquad] = pde_physics->convective_flux (soln_int[iquad]);
+        conv_phys_flux_ext[iquad] = pde_physics->convective_flux (soln_ext[iquad]);
+
+        diss_soln_num_flux[iquad] = diss_num_flux->evaluate_solution_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+
+        ADArrayTensor1 diss_soln_jump_int, diss_soln_jump_ext;
+        for (int s=0; s<nstate; s++) {
+            diss_soln_jump_int[s] = (diss_soln_num_flux[iquad][s] - soln_int[iquad][s]) * normal_int;
+            diss_soln_jump_ext[s] = (diss_soln_num_flux[iquad][s] - soln_ext[iquad][s]) * normal_ext;
+        }
+        diss_flux_jump_int[iquad] = pde_physics->dissipative_flux (soln_int[iquad], diss_soln_jump_int);
+        diss_flux_jump_ext[iquad] = pde_physics->dissipative_flux (soln_ext[iquad], diss_soln_jump_ext);
+
+        diss_auxi_num_flux_dot_n[iquad] = diss_num_flux->evaluate_auxiliary_flux(
+            soln_int[iquad], soln_ext[iquad],
+            soln_grad_int[iquad], soln_grad_ext[iquad],
+            normal_int, penalty);
+    }
+
+    // From test functions associated with interior cell point of view
+    for (unsigned int itest_int=0; itest_int<n_dofs_int; ++itest_int) {
+        ADtype rhs = 0.0;
+        const unsigned int istate = fe_values_int.get_fe().system_to_component_index(itest_int).first;
+
+        for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+            // Convection
+            const ADtype flux_diff = conv_num_flux_dot_n[iquad][istate] - conv_phys_flux_int[iquad][istate]*normals_int[iquad];
+            rhs = rhs - fe_values_int.shape_value_component(itest_int,iquad,istate) * flux_diff * JxW_int[iquad];
+            // Diffusive
+            rhs = rhs - fe_values_int.shape_value_component(itest_int,iquad,istate) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW_int[iquad];
+            rhs = rhs + fe_values_int.shape_grad_component(itest_int,iquad,istate) * diss_flux_jump_int[iquad][istate] * JxW_int[iquad];
+        }
+
+        local_rhs_int_cell(itest_int) += rhs.val();
+        for (unsigned int idof = 0; idof < n_dofs_int; ++idof) {
+            dR1_dW1[idof] = rhs.dx(idof);
+        }
+        for (unsigned int idof = 0; idof < n_dofs_ext; ++idof) {
+            dR1_dW2[idof] = rhs.dx(n_dofs_int+idof);
+        }
+        this->system_matrix.add(dof_indices_int[itest_int], dof_indices_int, dR1_dW1);
+        this->system_matrix.add(dof_indices_int[itest_int], dof_indices_ext, dR1_dW2);
+    }
+
+    // From test functions associated with neighbour cell point of view
+    for (unsigned int itest_ext=0; itest_ext<n_dofs_ext; ++itest_ext) {
+        ADtype rhs = 0.0;
+        const unsigned int istate = fe_values_int.get_fe().system_to_component_index(itest_ext).first;
+
+        for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+            // Convection
+            const ADtype flux_diff = (-conv_num_flux_dot_n[iquad][istate]) - conv_phys_flux_ext[iquad][istate]*(-normals_int[iquad]);
+            rhs = rhs - fe_values_ext.shape_value_component(itest_ext,iquad,istate) * flux_diff * JxW_int[iquad];
+            // Diffusive
+            rhs = rhs - fe_values_ext.shape_value_component(itest_ext,iquad,istate) * (-diss_auxi_num_flux_dot_n[iquad][istate]) * JxW_int[iquad];
+            rhs = rhs + fe_values_ext.shape_grad_component(itest_ext,iquad,istate) * diss_flux_jump_ext[iquad][istate] * JxW_int[iquad];
+        }
+
+        local_rhs_ext_cell(itest_ext) += rhs.val();
+        for (unsigned int idof = 0; idof < n_dofs_int; ++idof) {
+            dR2_dW1[idof] = rhs.dx(idof);
+        }
+        for (unsigned int idof = 0; idof < n_dofs_ext; ++idof) {
+            dR2_dW2[idof] = rhs.dx(n_dofs_int+idof);
+        }
+        this->system_matrix.add(dof_indices_ext[itest_ext], dof_indices_int, dR2_dW1);
+        this->system_matrix.add(dof_indices_ext[itest_ext], dof_indices_ext, dR2_dW2);
+    }
+}
+template class DGStrong <PHILIP_DIM, 1, double>;
+template class DGStrong <PHILIP_DIM, 2, double>;
+template class DGStrong <PHILIP_DIM, 3, double>;
+template class DGStrong <PHILIP_DIM, 4, double>;
+template class DGStrong <PHILIP_DIM, 5, double>;
+
+} // PHiLiP namespace
+
diff --git a/src/dg/assemble_implicit.cpp b/src/dg/weak_dg.cpp
similarity index 73%
rename from src/dg/assemble_implicit.cpp
rename to src/dg/weak_dg.cpp
index b039d8d9a..10bd76806 100644
--- a/src/dg/assemble_implicit.cpp
+++ b/src/dg/weak_dg.cpp
@@ -17,8 +17,31 @@
 namespace PHiLiP {
 
 template <int dim, int nstate, typename real>
-void DG<dim,nstate,real>::assemble_cell_terms_implicit(
-    const dealii::FEValues<dim,dim> *fe_values_vol,
+DGWeak<dim,nstate,real>::DGWeak(
+    const Parameters::AllParameters *const parameters_input,
+    const unsigned int degree)
+    : DGBase<dim,real>::DGBase(nstate, parameters_input, degree) // Use DGBase constructor
+{
+    using ADtype = Sacado::Fad::DFad<real>;
+    pde_physics = Physics::PhysicsFactory<dim,nstate,ADtype> 
+        ::create_Physics(parameters_input->pde_type);
+    conv_num_flux = NumericalFlux::NumericalFluxFactory<dim, nstate, ADtype>
+        ::create_convective_numerical_flux (parameters_input->conv_num_flux_type, pde_physics);
+    diss_num_flux = NumericalFlux::NumericalFluxFactory<dim, nstate, ADtype>
+        ::create_dissipative_numerical_flux (parameters_input->diss_num_flux_type, pde_physics);
+}
+// Destructor
+template <int dim, int nstate, typename real>
+DGWeak<dim,nstate,real>::~DGWeak ()
+{ 
+    std::cout << "Destructing DGWeak..." << std::endl;
+    delete conv_num_flux;
+    delete diss_num_flux;
+}
+
+template <int dim, int nstate, typename real>
+void DGWeak<dim,nstate,real>::assemble_cell_terms_implicit(
+    const dealii::FEValues<dim,dim> &fe_values_vol,
     const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
     dealii::Vector<real> &local_rhs_int_cell)
 {
@@ -26,12 +49,12 @@ void DG<dim,nstate,real>::assemble_cell_terms_implicit(
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
 
-    const unsigned int n_quad_pts      = fe_values_vol->n_quadrature_points;
-    const unsigned int n_dofs_cell     = fe_values_vol->dofs_per_cell;
+    const unsigned int n_quad_pts      = fe_values_vol.n_quadrature_points;
+    const unsigned int n_dofs_cell     = fe_values_vol.dofs_per_cell;
 
     AssertDimension (n_dofs_cell, cell_dofs_indices.size());
 
-    const std::vector<real> &JxW = fe_values_vol->get_JxW_values ();
+    const std::vector<real> &JxW = fe_values_vol.get_JxW_values ();
 
 
     std::vector<real> residual_derivatives(n_dofs_cell);
@@ -60,9 +83,9 @@ void DG<dim,nstate,real>::assemble_cell_terms_implicit(
     // Interpolate solution to face
     for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
         for (unsigned int idof=0; idof<n_dofs_cell; ++idof) {
-              const unsigned int istate = fe_values_vol->get_fe().system_to_component_index(idof).first;
-              soln_at_q[iquad][istate]      += soln_coeff[idof] * fe_values_vol->shape_value_component(idof, iquad, istate);
-              soln_grad_at_q[iquad][istate] += soln_coeff[idof] * fe_values_vol->shape_grad_component(idof, iquad, istate);
+              const unsigned int istate = fe_values_vol.get_fe().system_to_component_index(idof).first;
+              soln_at_q[iquad][istate]      += soln_coeff[idof] * fe_values_vol.shape_value_component(idof, iquad, istate);
+              soln_grad_at_q[iquad][istate] += soln_coeff[idof] * fe_values_vol.shape_grad_component(idof, iquad, istate);
         }
         //std::cout << "Density " << soln_at_q[iquad][0] << std::endl;
         //if(nstate>1) std::cout << "Momentum " << soln_at_q[iquad][1] << std::endl;
@@ -70,7 +93,7 @@ void DG<dim,nstate,real>::assemble_cell_terms_implicit(
         // Evaluate physical convective flux and source term
         conv_phys_flux_at_q[iquad] = pde_physics->convective_flux (soln_at_q[iquad]);
         diss_phys_flux_at_q[iquad] = pde_physics->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad]);
-        source_at_q[iquad] = pde_physics->source_term (fe_values_vol->quadrature_point(iquad), soln_at_q[iquad]);
+        source_at_q[iquad] = pde_physics->source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad]);
     }
 
     // Weak form
@@ -85,17 +108,17 @@ void DG<dim,nstate,real>::assemble_cell_terms_implicit(
 
         ADtype rhs = 0;
 
-        const unsigned int istate = fe_values_vol->get_fe().system_to_component_index(itest).first;
+        const unsigned int istate = fe_values_vol.get_fe().system_to_component_index(itest).first;
 
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
 
             // Convective
-            rhs = rhs + fe_values_vol->shape_grad_component(itest,iquad,istate) * conv_phys_flux_at_q[iquad][istate] * JxW[iquad];
+            rhs = rhs + fe_values_vol.shape_grad_component(itest,iquad,istate) * conv_phys_flux_at_q[iquad][istate] * JxW[iquad];
             //// Diffusive
             //// Note that for diffusion, the negative is defined in the physics
-            rhs = rhs + fe_values_vol->shape_grad_component(itest,iquad,istate) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
+            rhs = rhs + fe_values_vol.shape_grad_component(itest,iquad,istate) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
             // Source
-            rhs = rhs + fe_values_vol->shape_value_component(itest,iquad,istate) * source_at_q[iquad][istate] * JxW[iquad];
+            rhs = rhs + fe_values_vol.shape_value_component(itest,iquad,istate) * source_at_q[iquad][istate] * JxW[iquad];
         }
 
         local_rhs_int_cell(itest) += rhs.val();
@@ -110,8 +133,8 @@ void DG<dim,nstate,real>::assemble_cell_terms_implicit(
 
 
 template <int dim, int nstate, typename real>
-void DG<dim,nstate,real>::assemble_boundary_term_implicit(
-    const dealii::FEFaceValues<dim,dim> *fe_values_boundary,
+void DGWeak<dim,nstate,real>::assemble_boundary_term_implicit(
+    const dealii::FEFaceValues<dim,dim> &fe_values_boundary,
     const real penalty,
     const std::vector<dealii::types::global_dof_index> &dof_indices_int,
     dealii::Vector<real> &local_rhs_int_cell)
@@ -120,13 +143,13 @@ void DG<dim,nstate,real>::assemble_boundary_term_implicit(
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
 
-    const unsigned int n_dofs_cell = fe_values_boundary->dofs_per_cell;
-    const unsigned int n_face_quad_pts = fe_values_boundary->n_quadrature_points;
+    const unsigned int n_dofs_cell = fe_values_boundary.dofs_per_cell;
+    const unsigned int n_face_quad_pts = fe_values_boundary.n_quadrature_points;
 
     AssertDimension (n_dofs_cell, dof_indices_int.size());
 
-    const std::vector<real> &JxW = fe_values_boundary->get_JxW_values ();
-    const std::vector<dealii::Tensor<1,dim>> &normals = fe_values_boundary->get_normal_vectors ();
+    const std::vector<real> &JxW = fe_values_boundary.get_JxW_values ();
+    const std::vector<dealii::Tensor<1,dim>> &normals = fe_values_boundary.get_normal_vectors ();
 
     std::vector<real> residual_derivatives(n_dofs_cell);
 
@@ -157,16 +180,16 @@ void DG<dim,nstate,real>::assemble_boundary_term_implicit(
         }
     }
     // Interpolate solution to face
-    const std::vector< dealii::Point<dim,real> > quad_pts = fe_values_boundary->get_quadrature_points();
+    const std::vector< dealii::Point<dim,real> > quad_pts = fe_values_boundary.get_quadrature_points();
     for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
 
         const dealii::Tensor<1,dim,ADtype> normal_int = normals[iquad];
         const dealii::Tensor<1,dim,ADtype> normal_ext = -normal_int;
 
         for (unsigned int idof=0; idof<n_dofs_cell; ++idof) {
-            const int istate = fe_values_boundary->get_fe().system_to_component_index(idof).first;
-            soln_int[iquad][istate]      += soln_coeff_int[idof] * fe_values_boundary->shape_value_component(idof, iquad, istate);
-            soln_grad_int[iquad][istate] += soln_coeff_int[idof] * fe_values_boundary->shape_grad_component(idof, iquad, istate);
+            const int istate = fe_values_boundary.get_fe().system_to_component_index(idof).first;
+            soln_int[iquad][istate]      += soln_coeff_int[idof] * fe_values_boundary.shape_value_component(idof, iquad, istate);
+            soln_grad_int[iquad][istate] += soln_coeff_int[idof] * fe_values_boundary.shape_grad_component(idof, iquad, istate);
         }
 
         const unsigned int boundary_id = 1; // dummy value for now
@@ -206,15 +229,15 @@ void DG<dim,nstate,real>::assemble_boundary_term_implicit(
 
         ADtype rhs = 0.0;
 
-        const unsigned int istate = fe_values_boundary->get_fe().system_to_component_index(itest).first;
+        const unsigned int istate = fe_values_boundary.get_fe().system_to_component_index(itest).first;
 
         for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
 
             // Convection
-            rhs = rhs - fe_values_boundary->shape_value_component(itest,iquad,istate) * conv_num_flux_dot_n[iquad][istate] * JxW[iquad];
+            rhs = rhs - fe_values_boundary.shape_value_component(itest,iquad,istate) * conv_num_flux_dot_n[iquad][istate] * JxW[iquad];
             // Diffusive
-            rhs = rhs - fe_values_boundary->shape_value_component(itest,iquad,istate) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW[iquad];
-            rhs = rhs + fe_values_boundary->shape_grad_component(itest,iquad,istate) * diss_flux_jump_int[iquad][istate] * JxW[iquad];
+            rhs = rhs - fe_values_boundary.shape_value_component(itest,iquad,istate) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW[iquad];
+            rhs = rhs + fe_values_boundary.shape_grad_component(itest,iquad,istate) * diss_flux_jump_int[iquad][istate] * JxW[iquad];
         }
         // *******************
 
@@ -229,9 +252,9 @@ void DG<dim,nstate,real>::assemble_boundary_term_implicit(
 }
 
 template <int dim, int nstate, typename real>
-void DG<dim,nstate,real>::assemble_face_term_implicit(
-    const dealii::FEValuesBase<dim,dim>     *fe_values_int,
-    const dealii::FEFaceValues<dim,dim>     *fe_values_ext,
+void DGWeak<dim,nstate,real>::assemble_face_term_implicit(
+    const dealii::FEValuesBase<dim,dim>     &fe_values_int,
+    const dealii::FEFaceValues<dim,dim>     &fe_values_ext,
     const real penalty,
     const std::vector<dealii::types::global_dof_index> &dof_indices_int,
     const std::vector<dealii::types::global_dof_index> &dof_indices_ext,
@@ -244,18 +267,18 @@ void DG<dim,nstate,real>::assemble_face_term_implicit(
 
     // Use quadrature points of neighbor cell
     // Might want to use the maximum n_quad_pts1 and n_quad_pts2
-    const unsigned int n_face_quad_pts = fe_values_ext->n_quadrature_points;
+    const unsigned int n_face_quad_pts = fe_values_ext.n_quadrature_points;
 
-    const unsigned int n_dofs_int = fe_values_int->dofs_per_cell;
-    const unsigned int n_dofs_ext = fe_values_ext->dofs_per_cell;
+    const unsigned int n_dofs_int = fe_values_int.dofs_per_cell;
+    const unsigned int n_dofs_ext = fe_values_ext.dofs_per_cell;
 
     AssertDimension (n_dofs_int, dof_indices_int.size());
     AssertDimension (n_dofs_ext, dof_indices_ext.size());
 
     // Jacobian and normal should always be consistent between two elements
     // even for non-conforming meshes?
-    const std::vector<real> &JxW_int = fe_values_int->get_JxW_values ();
-    const std::vector<dealii::Tensor<1,dim> > &normals_int = fe_values_int->get_normal_vectors ();
+    const std::vector<real> &JxW_int = fe_values_int.get_JxW_values ();
+    const std::vector<dealii::Tensor<1,dim> > &normals_int = fe_values_int.get_normal_vectors ();
 
     // AD variable
     std::vector<ADtype> soln_coeff_int_ad(n_dofs_int);
@@ -307,14 +330,14 @@ void DG<dim,nstate,real>::assemble_face_term_implicit(
 
         // Interpolate solution to face
         for (unsigned int idof=0; idof<n_dofs_int; ++idof) {
-            const unsigned int istate = fe_values_int->get_fe().system_to_component_index(idof).first;
-            soln_int[iquad][istate]      += soln_coeff_int_ad[idof] * fe_values_int->shape_value_component(idof, iquad, istate);
-            soln_grad_int[iquad][istate] += soln_coeff_int_ad[idof] * fe_values_int->shape_grad_component(idof, iquad, istate);
+            const unsigned int istate = fe_values_int.get_fe().system_to_component_index(idof).first;
+            soln_int[iquad][istate]      += soln_coeff_int_ad[idof] * fe_values_int.shape_value_component(idof, iquad, istate);
+            soln_grad_int[iquad][istate] += soln_coeff_int_ad[idof] * fe_values_int.shape_grad_component(idof, iquad, istate);
         }
         for (unsigned int idof=0; idof<n_dofs_ext; ++idof) {
-            const unsigned int istate = fe_values_ext->get_fe().system_to_component_index(idof).first;
-            soln_ext[iquad][istate]      += soln_coeff_ext_ad[idof] * fe_values_ext->shape_value_component(idof, iquad, istate);
-            soln_grad_ext[iquad][istate] += soln_coeff_ext_ad[idof] * fe_values_ext->shape_grad_component(idof, iquad, istate);
+            const unsigned int istate = fe_values_ext.get_fe().system_to_component_index(idof).first;
+            soln_ext[iquad][istate]      += soln_coeff_ext_ad[idof] * fe_values_ext.shape_value_component(idof, iquad, istate);
+            soln_grad_ext[iquad][istate] += soln_coeff_ext_ad[idof] * fe_values_ext.shape_grad_component(idof, iquad, istate);
         }
         //std::cout << "Density int" << soln_int[iquad][0] << std::endl;
         //if(nstate>1) std::cout << "Momentum int" << soln_int[iquad][1] << std::endl;
@@ -344,14 +367,14 @@ void DG<dim,nstate,real>::assemble_face_term_implicit(
     // From test functions associated with interior cell point of view
     for (unsigned int itest_int=0; itest_int<n_dofs_int; ++itest_int) {
         ADtype rhs = 0.0;
-        const unsigned int istate = fe_values_int->get_fe().system_to_component_index(itest_int).first;
+        const unsigned int istate = fe_values_int.get_fe().system_to_component_index(itest_int).first;
 
         for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
             // Convection
-            rhs = rhs - fe_values_int->shape_value_component(itest_int,iquad,istate) * conv_num_flux_dot_n[iquad][istate] * JxW_int[iquad];
+            rhs = rhs - fe_values_int.shape_value_component(itest_int,iquad,istate) * conv_num_flux_dot_n[iquad][istate] * JxW_int[iquad];
             // Diffusive
-            rhs = rhs - fe_values_int->shape_value_component(itest_int,iquad,istate) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW_int[iquad];
-            rhs = rhs + fe_values_int->shape_grad_component(itest_int,iquad,istate) * diss_flux_jump_int[iquad][istate] * JxW_int[iquad];
+            rhs = rhs - fe_values_int.shape_value_component(itest_int,iquad,istate) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW_int[iquad];
+            rhs = rhs + fe_values_int.shape_grad_component(itest_int,iquad,istate) * diss_flux_jump_int[iquad][istate] * JxW_int[iquad];
         }
 
         local_rhs_int_cell(itest_int) += rhs.val();
@@ -368,14 +391,14 @@ void DG<dim,nstate,real>::assemble_face_term_implicit(
     // From test functions associated with neighbour cell point of view
     for (unsigned int itest_ext=0; itest_ext<n_dofs_ext; ++itest_ext) {
         ADtype rhs = 0.0;
-        const unsigned int istate = fe_values_int->get_fe().system_to_component_index(itest_ext).first;
+        const unsigned int istate = fe_values_int.get_fe().system_to_component_index(itest_ext).first;
 
         for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
             // Convection
-            rhs = rhs - fe_values_ext->shape_value_component(itest_ext,iquad,istate) * (-conv_num_flux_dot_n[iquad][istate]) * JxW_int[iquad];
+            rhs = rhs - fe_values_ext.shape_value_component(itest_ext,iquad,istate) * (-conv_num_flux_dot_n[iquad][istate]) * JxW_int[iquad];
             // Diffusive
-            rhs = rhs - fe_values_ext->shape_value_component(itest_ext,iquad,istate) * (-diss_auxi_num_flux_dot_n[iquad][istate]) * JxW_int[iquad];
-            rhs = rhs + fe_values_ext->shape_grad_component(itest_ext,iquad,istate) * diss_flux_jump_ext[iquad][istate] * JxW_int[iquad];
+            rhs = rhs - fe_values_ext.shape_value_component(itest_ext,iquad,istate) * (-diss_auxi_num_flux_dot_n[iquad][istate]) * JxW_int[iquad];
+            rhs = rhs + fe_values_ext.shape_grad_component(itest_ext,iquad,istate) * diss_flux_jump_ext[iquad][istate] * JxW_int[iquad];
         }
 
         local_rhs_ext_cell(itest_ext) += rhs.val();
@@ -389,11 +412,11 @@ void DG<dim,nstate,real>::assemble_face_term_implicit(
         this->system_matrix.add(dof_indices_ext[itest_ext], dof_indices_ext, dR2_dW2);
     }
 }
-template class DG <PHILIP_DIM, 1, double>;
-template class DG <PHILIP_DIM, 2, double>;
-template class DG <PHILIP_DIM, 3, double>;
-template class DG <PHILIP_DIM, 4, double>;
-template class DG <PHILIP_DIM, 5, double>;
+template class DGWeak <PHILIP_DIM, 1, double>;
+template class DGWeak <PHILIP_DIM, 2, double>;
+template class DGWeak <PHILIP_DIM, 3, double>;
+template class DGWeak <PHILIP_DIM, 4, double>;
+template class DGWeak <PHILIP_DIM, 5, double>;
 
 } // PHiLiP namespace
 
diff --git a/src/numerical_flux/numerical_flux.cpp b/src/numerical_flux/numerical_flux.cpp
index 7e7ae1bc3..0b8fd9e65 100644
--- a/src/numerical_flux/numerical_flux.cpp
+++ b/src/numerical_flux/numerical_flux.cpp
@@ -27,7 +27,7 @@ NumericalFluxConvective<dim,nstate,real>*
 NumericalFluxFactory<dim, nstate, real>
 ::create_convective_numerical_flux(
     AllParam::ConvectiveNumericalFlux conv_num_flux_type,
-    Physics::PhysicsBase<dim, nstate, real> *physics_input)
+    std::shared_ptr<Physics::PhysicsBase<dim, nstate, real>> physics_input)
 {
     if(conv_num_flux_type == AllParam::lax_friedrichs) {
         return new LaxFriedrichs<dim, nstate, real>(physics_input);
@@ -40,7 +40,7 @@ NumericalFluxDissipative<dim,nstate,real>*
 NumericalFluxFactory<dim, nstate, real>
 ::create_dissipative_numerical_flux(
     AllParam::DissipativeNumericalFlux diss_num_flux_type,
-    Physics::PhysicsBase<dim, nstate, real> *physics_input)
+    std::shared_ptr <Physics::PhysicsBase<dim, nstate, real>> physics_input)
 {
     if(diss_num_flux_type == AllParam::symm_internal_penalty) {
         return new SymmetricInternalPenalty<dim, nstate, real>(physics_input);
diff --git a/src/numerical_flux/numerical_flux.h b/src/numerical_flux/numerical_flux.h
index 63cc7091b..8edabb2d9 100644
--- a/src/numerical_flux/numerical_flux.h
+++ b/src/numerical_flux/numerical_flux.h
@@ -33,7 +33,7 @@ class LaxFriedrichs: public NumericalFluxConvective<dim, nstate, real>
 public:
 
 /// Constructor
-LaxFriedrichs(Physics::PhysicsBase<dim, nstate, real> *physics_input)
+LaxFriedrichs(std::shared_ptr <Physics::PhysicsBase<dim, nstate, real>> physics_input)
 :
 pde_physics(physics_input)
 {};
@@ -48,7 +48,7 @@ std::array<real, nstate> evaluate_flux (
 
 protected:
 /// Numerical flux requires physics to evaluate convective eigenvalues.
-const Physics::PhysicsBase<dim, nstate, real> *pde_physics;
+const std::shared_ptr < Physics::PhysicsBase<dim, nstate, real> > pde_physics;
 
 };
 
@@ -62,12 +62,12 @@ class NumericalFluxFactory
     static NumericalFluxConvective<dim,nstate,real>*
         create_convective_numerical_flux
             (AllParam::ConvectiveNumericalFlux conv_num_flux_type,
-            Physics::PhysicsBase<dim, nstate, real> *physics_input);
+            std::shared_ptr<Physics::PhysicsBase<dim, nstate, real>> physics_input);
     /// Creates dissipative numerical flux based on input.
     static NumericalFluxDissipative<dim,nstate,real>*
         create_dissipative_numerical_flux
             (AllParam::DissipativeNumericalFlux diss_num_flux_type,
-            Physics::PhysicsBase<dim, nstate, real> *physics_input);
+            std::shared_ptr<Physics::PhysicsBase<dim, nstate, real>> physics_input);
 };
 
 } // NumericalFlux namespace
diff --git a/src/numerical_flux/viscous_numerical_flux.h b/src/numerical_flux/viscous_numerical_flux.h
index 0a5cc8aa6..717c48cdd 100644
--- a/src/numerical_flux/viscous_numerical_flux.h
+++ b/src/numerical_flux/viscous_numerical_flux.h
@@ -43,7 +43,7 @@ class SymmetricInternalPenalty: public NumericalFluxDissipative<dim, nstate, rea
 {
 public:
 /// Constructor
-SymmetricInternalPenalty(Physics::PhysicsBase<dim, nstate, real> *physics_input)
+SymmetricInternalPenalty(std::shared_ptr<Physics::PhysicsBase<dim, nstate, real>> physics_input)
 :
 pde_physics(physics_input)
 {};
@@ -74,7 +74,7 @@ std::array<real, nstate> evaluate_auxiliary_flux (
     const bool on_boundary = false) const;
     
 protected:
-const Physics::PhysicsBase<dim, nstate, real> *pde_physics;
+const std::shared_ptr < Physics::PhysicsBase<dim, nstate, real> > pde_physics;
 
 };
 
diff --git a/src/ode_solver/ode_solver.cpp b/src/ode_solver/ode_solver.cpp
index cd6ffca31..6eef51866 100644
--- a/src/ode_solver/ode_solver.cpp
+++ b/src/ode_solver/ode_solver.cpp
@@ -32,7 +32,7 @@ int Implicit_ODESolver<dim,real>::steady_state ()
         if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
             (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
         std::cout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
-        ODESolver<dim,real>::dg->assemble_system ();
+        ODESolver<dim,real>::dg->assemble_residual_dRdW ();
         // (M/dt - dRdW) dw = R
         ODESolver<dim,real>::dg->system_matrix *= -1.0;
 
@@ -63,7 +63,7 @@ int Explicit_ODESolver<dim,real>::steady_state ()
     Parameters::ODESolverParam ode_param = ODESolver<dim,real>::all_parameters->ode_solver_param;
     allocate_ode_system ();
 
-    ODESolver<dim,real>::dg->assemble_system ();
+    ODESolver<dim,real>::dg->assemble_residual_dRdW ();
 
     this->residual_norm = 1.0;
     this->current_iteration = 0;
@@ -84,7 +84,7 @@ int Explicit_ODESolver<dim,real>::steady_state ()
         evaluate_solution_update ();
         ODESolver<dim,real>::dg->solution += this->solution_update;
 
-        ODESolver<dim,real>::dg->assemble_system ();
+        ODESolver<dim,real>::dg->assemble_residual_dRdW ();
         this->residual_norm = ODESolver<dim,real>::dg->get_residual_l2norm();
     }
     return 1;
diff --git a/src/parameters/all_parameters.cpp b/src/parameters/all_parameters.cpp
index 2a75bf127..e2fed1b64 100644
--- a/src/parameters/all_parameters.cpp
+++ b/src/parameters/all_parameters.cpp
@@ -15,13 +15,17 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm)
     prm.declare_entry("dimension", "1",
                       dealii::Patterns::Integer(),
                       "Number of dimensions");
+
+    prm.declare_entry("use_weak_form", "true",
+                      dealii::Patterns::Bool(),
+                      "Use weak form by default. If false, use strong form.");
     prm.declare_entry("test_type", "run_control",
                       dealii::Patterns::Selection(
                       " run_control | "
                       " numerical_flux_convervation | "
                       " jacobian_regression "),
                       "The type of test we want to solve. "
-                      "Choices are " 
+                      "Choices are (only run control has been coded up for now)" 
                       " <run_control | " 
                       "  numerical_flux_convervation | "
                       "  jacobian_regression>.");
@@ -89,6 +93,7 @@ void AllParameters::parse_parameters (dealii::ParameterHandler &prm)
         pde_type = euler;
         nstate = dimension+2;
     }
+    use_weak_form = prm.get_bool("use_weak_form");
 
     const std::string conv_num_flux_string = prm.get("conv_num_flux");
     if (conv_num_flux_string == "lax_friedrichs") conv_num_flux_type = lax_friedrichs;
diff --git a/src/parameters/all_parameters.h b/src/parameters/all_parameters.h
index f7dad49e5..0518c4b0d 100644
--- a/src/parameters/all_parameters.h
+++ b/src/parameters/all_parameters.h
@@ -28,6 +28,8 @@ class AllParameters
     /// Number of dimensions. Note that it has to match the executable PHiLiP_xD
     unsigned int dimension;
 
+    /// Flag to use weak or strong form of DG
+    bool use_weak_form;
 
     /// Number of state variables. Will depend on PDE
     int nstate;
diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt
index 95504be94..64165948e 100644
--- a/src/physics/CMakeLists.txt
+++ b/src/physics/CMakeLists.txt
@@ -15,7 +15,7 @@ foreach(dim RANGE 1 3)
     DEAL_II_SETUP_TARGET(${PhysicsLib})
 
     target_compile_definitions(${PhysicsLib} PRIVATE PHILIP_DIM=${dim})
-    target_link_libraries(${PhysicsLib} tbb m mpi mpi_cxx)
+    target_link_libraries(${PhysicsLib} m mpi mpi_cxx)
 
     # Setup target with deal.II
     unset(PhysicsLib)
diff --git a/src/physics/physics.cpp b/src/physics/physics.cpp
index ebcd7cd71..dc6b25e3e 100644
--- a/src/physics/physics.cpp
+++ b/src/physics/physics.cpp
@@ -13,22 +13,22 @@ namespace PHiLiP {
 namespace Physics {
 
 template <int dim, int nstate, typename real>
-PhysicsBase<dim,nstate,real>* // returns points to base class PhysicsBase
+std::shared_ptr < PhysicsBase<dim,nstate,real> >
 PhysicsFactory<dim,nstate,real>
 ::create_Physics(Parameters::AllParameters::PartialDifferentialEquation pde_type)
 {
     using PDE_enum = Parameters::AllParameters::PartialDifferentialEquation;
 
     if (pde_type == PDE_enum::advection || pde_type == PDE_enum::advection_vector) {
-        if constexpr (nstate<=2) return new ConvectionDiffusion<dim,nstate,real>(true,false);
+        if constexpr (nstate<=2) return std::make_shared < ConvectionDiffusion<dim,nstate,real> >(true,false);
     } else if (pde_type == PDE_enum::diffusion) {
-        if constexpr (nstate==1) return new ConvectionDiffusion<dim,nstate,real>(false,true);
+        if constexpr (nstate==1) return std::make_shared < ConvectionDiffusion<dim,nstate,real> >(false,true);
     } else if (pde_type == PDE_enum::convection_diffusion) {
-        if constexpr (nstate==1) return new ConvectionDiffusion<dim,nstate,real>(true,true);
+        if constexpr (nstate==1) return std::make_shared < ConvectionDiffusion<dim,nstate,real> >(true,true);
     } else if (pde_type == PDE_enum::burgers_inviscid) {
-        if constexpr (nstate==dim) return new Burgers<dim,nstate,real>(true,false);
+        if constexpr (nstate==dim) return std::make_shared < Burgers<dim,nstate,real> >(true,false);
     } else if (pde_type == PDE_enum::euler) {
-        if constexpr (nstate==dim+2) return new Euler<dim,nstate,real>;
+        if constexpr (nstate==dim+2) return std::make_shared < Euler<dim,nstate,real> >();
     }
     std::cout << "Can't create PhysicsBase, invalid PDE type: " << pde_type << std::endl;
     assert(0==1 && "Can't create PhysicsBase, invalid PDE type");
diff --git a/src/physics/physics.h b/src/physics/physics.h
index 10f4babb1..2e99ea191 100644
--- a/src/physics/physics.h
+++ b/src/physics/physics.h
@@ -115,7 +115,7 @@ template <int dim, int nstate, typename real>
 class PhysicsFactory
 {
 public:
-    static PhysicsBase<dim,nstate,real>*
+    static std::shared_ptr< PhysicsBase<dim,nstate,real> >
         create_Physics(Parameters::AllParameters::PartialDifferentialEquation pde_type);
 };
 
diff --git a/src/tests/grid_study.cpp b/src/tests/grid_study.cpp
index 1fc5953b8..74c38d24d 100644
--- a/src/tests/grid_study.cpp
+++ b/src/tests/grid_study.cpp
@@ -110,7 +110,7 @@ ::run_test () const
     const double       grid_progression    = manu_grid_conv_param.grid_progression;
 
     std::cout<<"Test Physics nstate" << nstate << std::endl;
-    Physics::PhysicsBase<dim,nstate,double> *physics_double = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(param.pde_type);
+    std::shared_ptr <Physics::PhysicsBase<dim,nstate,double>> physics_double = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(param.pde_type);
 
     // Evaluate solution integral on really fine mesh
     double exact_solution_integral;
diff --git a/tests/advection_implicit/1d_advection_implicit_strong.prm b/tests/advection_implicit/1d_advection_implicit_strong.prm
new file mode 100644
index 000000000..3a6989f34
--- /dev/null
+++ b/tests/advection_implicit/1d_advection_implicit_strong.prm
@@ -0,0 +1,44 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 1
+
+# The PDE we want to solve. Choices are
+# <advection|diffusion|convection_diffusion>.
+set pde_type  = advection
+
+set use_weak_form = false
+
+set conv_num_flux = lax_friedrichs
+
+subsection ODE solver
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 500000
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-13
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set ode_solver_type                         = implicit
+end
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end        = 3
+
+  # Starting degree for convergence study
+  set degree_start      = 0
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 2
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 2
+
+  # Number of grids in grid study
+  set number_of_grids   = 4
+end
diff --git a/tests/advection_implicit/2d_advection_implicit_strong.prm b/tests/advection_implicit/2d_advection_implicit_strong.prm
new file mode 100644
index 000000000..dd98672b2
--- /dev/null
+++ b/tests/advection_implicit/2d_advection_implicit_strong.prm
@@ -0,0 +1,45 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 2
+
+# The PDE we want to solve. Choices are
+# <advection|diffusion|convection_diffusion>.
+set pde_type  = advection
+
+set use_weak_form = false
+
+set conv_num_flux = lax_friedrichs
+
+subsection ODE solver
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 500000
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-13
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set ode_solver_type                         = implicit
+end
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end        = 4
+
+  # Starting degree for convergence study
+  set degree_start      = 0
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 1.5
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 2
+
+  # Number of grids in grid study
+  set number_of_grids   = 6
+end
+
diff --git a/tests/advection_implicit/3d_advection_implicit_strong.prm b/tests/advection_implicit/3d_advection_implicit_strong.prm
new file mode 100644
index 000000000..910522f7e
--- /dev/null
+++ b/tests/advection_implicit/3d_advection_implicit_strong.prm
@@ -0,0 +1,44 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 3
+
+# The PDE we want to solve. Choices are
+# <advection|diffusion|convection_diffusion>.
+set pde_type  = advection
+
+set use_weak_form = false
+
+set conv_num_flux = lax_friedrichs
+
+subsection ODE solver
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 500000
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-13
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  #set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set ode_solver_type                         = implicit
+end
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end        = 2
+
+  # Starting degree for convergence study
+  set degree_start      = 0
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 1.5
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 2
+
+  # Number of grids in grid study
+  set number_of_grids   = 5
+end
diff --git a/tests/advection_implicit/CMakeLists.txt b/tests/advection_implicit/CMakeLists.txt
index 2e836ea67..ce7839c07 100644
--- a/tests/advection_implicit/CMakeLists.txt
+++ b/tests/advection_implicit/CMakeLists.txt
@@ -20,8 +20,25 @@ add_test(
   COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_3D -i ${CMAKE_CURRENT_BINARY_DIR}/3d_advection_implicit.prm
   WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
 )
-#add_test(
-#  NAME 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION
-#  COMMAND PHiLiP -i 3d_advection_implicit.prm
-#  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
-#)
+
+# Tests for strong form
+configure_file(1d_advection_implicit_strong.prm 1d_advection_implicit_strong.prm COPYONLY)
+add_test(
+  NAME 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG
+  COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_1D -i ${CMAKE_CURRENT_BINARY_DIR}/1d_advection_implicit_strong.prm
+  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+)
+
+configure_file(2d_advection_implicit_strong.prm 2d_advection_implicit_strong.prm COPYONLY)
+add_test(
+  NAME 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG
+  COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_2D -i ${CMAKE_CURRENT_BINARY_DIR}/2d_advection_implicit_strong.prm
+  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+)
+
+configure_file(3d_advection_implicit_strong.prm 3d_advection_implicit_strong.prm COPYONLY)
+add_test(
+  NAME 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG
+  COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_3D -i ${CMAKE_CURRENT_BINARY_DIR}/3d_advection_implicit_strong.prm
+  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+)
diff --git a/tests/numerical_flux/CMakeLists.txt b/tests/numerical_flux/CMakeLists.txt
index 0a1d07ee2..4f6590391 100644
--- a/tests/numerical_flux/CMakeLists.txt
+++ b/tests/numerical_flux/CMakeLists.txt
@@ -31,7 +31,7 @@ foreach(dim RANGE 1 3)
     # See one of the following links (translate chinese to english using Google)
     # https://zhuanlan.zhihu.com/p/21385662
     # https://www.getit01.com/p2018020821385662/
-    target_link_libraries(${TEST_TARGET} tbb m mpi mpi_cxx)
+    target_link_libraries(${TEST_TARGET} m mpi mpi_cxx)
 
 
     add_test(
diff --git a/tests/numerical_flux/numerical_flux_conservation.cpp b/tests/numerical_flux/numerical_flux_conservation.cpp
index a7230607b..c7ba5d074 100644
--- a/tests/numerical_flux/numerical_flux_conservation.cpp
+++ b/tests/numerical_flux/numerical_flux_conservation.cpp
@@ -63,7 +63,7 @@ int test_dissipative_numerical_flux_conservation (const PDEType pde_type, const
 
 
     using namespace PHiLiP;
-    Physics::PhysicsBase<dim, nstate, double> *pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
+    std::shared_ptr < Physics::PhysicsBase<dim, nstate, double> > pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
 
     dealii::Tensor<1,dim,double> normal_int;
     std::array<double, nstate> soln_int, soln_ext;
@@ -106,7 +106,6 @@ int test_dissipative_numerical_flux_conservation (const PDEType pde_type, const
     compare_array<dim,nstate> (diss_auxi_num_flux_dot_n_1, diss_auxi_num_flux_dot_n_2, -1.0);
 
     delete diss_num_flux;
-    delete pde_physics;
 
     return 0;
 }
@@ -115,7 +114,7 @@ template<int dim, int nstate>
 int test_dissipative_numerical_flux_consistency (const PDEType pde_type, const DissType diss_type)
 {
     using namespace PHiLiP;
-    Physics::PhysicsBase<dim, nstate, double> *pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
+    std::shared_ptr <Physics::PhysicsBase<dim, nstate, double>> pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
 
     NumericalFlux::NumericalFluxDissipative<dim, nstate, double> *diss_num_flux = 
         NumericalFlux::NumericalFluxFactory<dim, nstate, double>
@@ -166,7 +165,6 @@ int test_dissipative_numerical_flux_consistency (const PDEType pde_type, const D
     compare_array<dim,nstate> (diss_auxi_num_flux_dot_n, diss_phys_flux_int, 1.0);
 
     delete diss_num_flux;
-    delete pde_physics;
 
     return 0;
 }
@@ -175,7 +173,7 @@ template<int dim, int nstate>
 int test_convective_numerical_flux_conservation (const PDEType pde_type, const ConvType conv_type)
 {
     using namespace PHiLiP;
-    Physics::PhysicsBase<dim, nstate, double> *pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
+    std::shared_ptr <Physics::PhysicsBase<dim, nstate, double>> pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
 
     NumericalFlux::NumericalFluxConvective<dim, nstate, double> *conv_num_flux = 
         NumericalFlux::NumericalFluxFactory<dim, nstate, double>
@@ -205,7 +203,6 @@ int test_convective_numerical_flux_conservation (const PDEType pde_type, const C
     compare_array<dim,nstate> (conv_num_flux_dot_n_1, conv_num_flux_dot_n_2, -1.0);
 
     delete conv_num_flux;
-    delete pde_physics;
 
     return 0;
 }
@@ -214,7 +211,7 @@ template<int dim, int nstate>
 int test_convective_numerical_flux_consistency (const PDEType pde_type, const ConvType conv_type)
 {
     using namespace PHiLiP;
-    Physics::PhysicsBase<dim, nstate, double> *pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
+    std::shared_ptr <Physics::PhysicsBase<dim, nstate, double>> pde_physics = Physics::PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
 
     NumericalFlux::NumericalFluxConvective<dim, nstate, double> *conv_num_flux = 
         NumericalFlux::NumericalFluxFactory<dim, nstate, double>
@@ -254,7 +251,6 @@ int test_convective_numerical_flux_consistency (const PDEType pde_type, const Co
     compare_array<dim,nstate> (conv_num_flux_dot_n, conv_phys_flux_dot_n, 1.0);
 
     delete conv_num_flux;
-    delete pde_physics;
 
     return 0;
 }
diff --git a/tests/regression/CMakeLists.txt b/tests/regression/CMakeLists.txt
index cae8957d6..82f8e5892 100644
--- a/tests/regression/CMakeLists.txt
+++ b/tests/regression/CMakeLists.txt
@@ -22,7 +22,7 @@ foreach(dim RANGE 1 3)
     string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D)
     target_link_libraries(${TEST_TARGET} ${ParametersLib})
     target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib})
-    target_link_libraries(${TEST_TARGET} tbb m mpi mpi_cxx)
+    target_link_libraries(${TEST_TARGET} m mpi mpi_cxx)
     # Setup target with deal.II
     DEAL_II_SETUP_TARGET(${TEST_TARGET})
 
diff --git a/tests/regression/jacobian_matrix_regression.cpp b/tests/regression/jacobian_matrix_regression.cpp
index 9024b0182..06928dcfa 100644
--- a/tests/regression/jacobian_matrix_regression.cpp
+++ b/tests/regression/jacobian_matrix_regression.cpp
@@ -54,7 +54,7 @@ int main (int argc, char * argv[])
                 std::shared_ptr < DGBase<PHILIP_DIM, double> > dg = DGFactory<PHILIP_DIM,double>::create_discontinuous_galerkin(&all_parameters, poly_degree);
                 dg->set_triangulation(&grid);
                 dg->allocate_system ();
-                dg->assemble_system ();
+                dg->assemble_residual_dRdW ();
 
                 const int nrows = dg->system_matrix.m();
 
diff --git a/tests/regression/matrix_data/1d_advection_poly_1_gridsize_2.mat b/tests/regression/matrix_data/1d_advection_poly_1_gridsize_2.mat
index 5c09f6ddb..39d2a047d 100644
Binary files a/tests/regression/matrix_data/1d_advection_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/1d_advection_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_poly_1_gridsize_3.mat b/tests/regression/matrix_data/1d_advection_poly_1_gridsize_3.mat
index b2553e93b..73610627a 100644
Binary files a/tests/regression/matrix_data/1d_advection_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/1d_advection_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_poly_1_gridsize_4.mat b/tests/regression/matrix_data/1d_advection_poly_1_gridsize_4.mat
index 0d0418a8d..b20d26c74 100644
Binary files a/tests/regression/matrix_data/1d_advection_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/1d_advection_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_poly_2_gridsize_2.mat b/tests/regression/matrix_data/1d_advection_poly_2_gridsize_2.mat
index 5410ded8d..f35053ec6 100644
Binary files a/tests/regression/matrix_data/1d_advection_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/1d_advection_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_poly_2_gridsize_3.mat b/tests/regression/matrix_data/1d_advection_poly_2_gridsize_3.mat
index f04fc3e35..986a5d95c 100644
Binary files a/tests/regression/matrix_data/1d_advection_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/1d_advection_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_poly_2_gridsize_4.mat b/tests/regression/matrix_data/1d_advection_poly_2_gridsize_4.mat
index fd6089ca6..3393a46c9 100644
Binary files a/tests/regression/matrix_data/1d_advection_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/1d_advection_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_2.mat b/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_2.mat
index 60403af6e..29e635fbd 100644
Binary files a/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_3.mat b/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_3.mat
index fb0bb961d..552cf4fe3 100644
Binary files a/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_4.mat b/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_4.mat
index b14bef3e0..a6dbdd5ed 100644
Binary files a/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/1d_advection_vector_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_2.mat b/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_2.mat
index 14c56a351..6a5df32d5 100644
Binary files a/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_3.mat b/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_3.mat
index ecbe38909..9cdfb57c2 100644
Binary files a/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_4.mat b/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_4.mat
index 725a1527d..18020aa4c 100644
Binary files a/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/1d_advection_vector_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_2.mat b/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_2.mat
index 0c764e08c..62f4c48bb 100644
Binary files a/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_3.mat b/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_3.mat
index 4614bb1cb..6a32c0a0f 100644
Binary files a/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_4.mat b/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_4.mat
index 5bd47a807..f59db23ad 100644
Binary files a/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/1d_convection_diffusion_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_2.mat b/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_2.mat
index 27f78dba2..e8998aaf7 100644
Binary files a/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_3.mat b/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_3.mat
index ba32aff1e..f0e4b8e4e 100644
Binary files a/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_4.mat b/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_4.mat
index 5e29249b7..171437686 100644
Binary files a/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/1d_convection_diffusion_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_2.mat b/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_2.mat
index 7c7bf46af..c7f69d4ac 100644
Binary files a/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_3.mat b/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_3.mat
index 0d9b61411..5fb0ab819 100644
Binary files a/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_4.mat b/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_4.mat
index 935ea1a02..d3748b0e3 100644
Binary files a/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/1d_diffusion_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_2.mat b/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_2.mat
index 8a8d8ba82..23ef30257 100644
Binary files a/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_3.mat b/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_3.mat
index ac3f3286c..c88821a6c 100644
Binary files a/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_4.mat b/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_4.mat
index d3931b23a..d6a75467c 100644
Binary files a/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/1d_diffusion_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_poly_1_gridsize_2.mat b/tests/regression/matrix_data/2d_advection_poly_1_gridsize_2.mat
index d5bb829f8..898a67431 100644
Binary files a/tests/regression/matrix_data/2d_advection_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/2d_advection_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_poly_1_gridsize_3.mat b/tests/regression/matrix_data/2d_advection_poly_1_gridsize_3.mat
index f53c24e25..247d75267 100644
Binary files a/tests/regression/matrix_data/2d_advection_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/2d_advection_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_poly_1_gridsize_4.mat b/tests/regression/matrix_data/2d_advection_poly_1_gridsize_4.mat
index 3e6138aa6..44d43255a 100644
Binary files a/tests/regression/matrix_data/2d_advection_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/2d_advection_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_poly_2_gridsize_2.mat b/tests/regression/matrix_data/2d_advection_poly_2_gridsize_2.mat
index 6e15f94d9..6a036a8ca 100644
Binary files a/tests/regression/matrix_data/2d_advection_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/2d_advection_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_poly_2_gridsize_3.mat b/tests/regression/matrix_data/2d_advection_poly_2_gridsize_3.mat
index ddbf3afd4..ea5f5801a 100644
Binary files a/tests/regression/matrix_data/2d_advection_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/2d_advection_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_poly_2_gridsize_4.mat b/tests/regression/matrix_data/2d_advection_poly_2_gridsize_4.mat
index a7d6645d8..a7f774bd9 100644
Binary files a/tests/regression/matrix_data/2d_advection_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/2d_advection_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_2.mat b/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_2.mat
index 8ef057539..551e54aa0 100644
Binary files a/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_3.mat b/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_3.mat
index 4569efd89..2f29805ec 100644
Binary files a/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_4.mat b/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_4.mat
index de8c11667..85fafc06c 100644
Binary files a/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/2d_advection_vector_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_2.mat b/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_2.mat
index 952e1aaf0..2ef12df66 100644
Binary files a/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_3.mat b/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_3.mat
index 43e1f8e14..1e61c9b40 100644
Binary files a/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_4.mat b/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_4.mat
index 5ad5a21d5..14c05f879 100644
Binary files a/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/2d_advection_vector_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_2.mat b/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_2.mat
index 051e18d5c..e3ab4ac3e 100644
Binary files a/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_3.mat b/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_3.mat
index 3b46268f9..e03c60d5a 100644
Binary files a/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_4.mat b/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_4.mat
index ea955eac1..c2f776cd0 100644
Binary files a/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/2d_convection_diffusion_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_2.mat b/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_2.mat
index 2e536f4cc..7f457af1b 100644
Binary files a/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_3.mat b/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_3.mat
index 344c478cb..2138fb9ce 100644
Binary files a/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_4.mat b/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_4.mat
index f4c0f7bf2..2917e9d9d 100644
Binary files a/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/2d_convection_diffusion_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_2.mat b/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_2.mat
index 268be3366..5b304c94d 100644
Binary files a/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_3.mat b/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_3.mat
index 4ed23dced..078322dea 100644
Binary files a/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_4.mat b/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_4.mat
index b231bfefc..9eae02a3b 100644
Binary files a/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/2d_diffusion_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_2.mat b/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_2.mat
index 53c3fc086..023c53d94 100644
Binary files a/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_3.mat b/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_3.mat
index eb0ef3055..77e3f3cbf 100644
Binary files a/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_4.mat b/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_4.mat
index 55db77e09..fe7e76934 100644
Binary files a/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/2d_diffusion_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_poly_1_gridsize_2.mat b/tests/regression/matrix_data/3d_advection_poly_1_gridsize_2.mat
index 2eddb3ce5..347530890 100644
Binary files a/tests/regression/matrix_data/3d_advection_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/3d_advection_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_poly_1_gridsize_3.mat b/tests/regression/matrix_data/3d_advection_poly_1_gridsize_3.mat
index b58ca8fe2..6ae88d143 100644
Binary files a/tests/regression/matrix_data/3d_advection_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/3d_advection_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_poly_1_gridsize_4.mat b/tests/regression/matrix_data/3d_advection_poly_1_gridsize_4.mat
index 4fb4fc7d2..277508103 100644
Binary files a/tests/regression/matrix_data/3d_advection_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/3d_advection_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_poly_2_gridsize_2.mat b/tests/regression/matrix_data/3d_advection_poly_2_gridsize_2.mat
index 69308bedf..c4b7b82f6 100644
Binary files a/tests/regression/matrix_data/3d_advection_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/3d_advection_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_poly_2_gridsize_3.mat b/tests/regression/matrix_data/3d_advection_poly_2_gridsize_3.mat
index 03a73baf5..ed8f677b1 100644
Binary files a/tests/regression/matrix_data/3d_advection_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/3d_advection_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_poly_2_gridsize_4.mat b/tests/regression/matrix_data/3d_advection_poly_2_gridsize_4.mat
index 5f6cc2d80..0cc42ca2a 100644
Binary files a/tests/regression/matrix_data/3d_advection_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/3d_advection_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_2.mat b/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_2.mat
index f496eb9c9..a4bd19203 100644
Binary files a/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_3.mat b/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_3.mat
index fa8ba653e..b09460508 100644
Binary files a/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_4.mat b/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_4.mat
index 5948c399f..b3dca60c4 100644
Binary files a/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/3d_advection_vector_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_2.mat b/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_2.mat
index 691d593b6..389ad6f4f 100644
Binary files a/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_3.mat b/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_3.mat
index 48b44f13f..97e119913 100644
Binary files a/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_4.mat b/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_4.mat
index e7f939322..8067db548 100644
Binary files a/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/3d_advection_vector_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_2.mat b/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_2.mat
index 8e4aa527d..236eefe46 100644
Binary files a/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_3.mat b/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_3.mat
index 148cad28d..92166ed63 100644
Binary files a/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_4.mat b/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_4.mat
index 23184c80a..52948f63c 100644
Binary files a/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/3d_convection_diffusion_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_2.mat b/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_2.mat
index 9a6d0ef27..f8b1804d0 100644
Binary files a/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_3.mat b/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_3.mat
index f3c831039..64962a765 100644
Binary files a/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_4.mat b/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_4.mat
index 5ac5f0de4..6850e5b14 100644
Binary files a/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/3d_convection_diffusion_poly_2_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_2.mat b/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_2.mat
index 118532f88..1405dc17f 100644
Binary files a/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_2.mat and b/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_3.mat b/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_3.mat
index 7da545f25..0aef35586 100644
Binary files a/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_3.mat and b/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_4.mat b/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_4.mat
index d78ed93be..15df0dda1 100644
Binary files a/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_4.mat and b/tests/regression/matrix_data/3d_diffusion_poly_1_gridsize_4.mat differ
diff --git a/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_2.mat b/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_2.mat
index c379d9d05..c75daabf8 100644
Binary files a/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_2.mat and b/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_2.mat differ
diff --git a/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_3.mat b/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_3.mat
index 46cbee827..229980263 100644
Binary files a/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_3.mat and b/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_3.mat differ
diff --git a/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_4.mat b/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_4.mat
index 1488ef30f..cdaca19e3 100644
Binary files a/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_4.mat and b/tests/regression/matrix_data/3d_diffusion_poly_2_gridsize_4.mat differ