Fix wrong mass matrix

attention: it is still present in old solver

mfds/active_set_operator.cpp

@@ -8,6 +8,9 @@
 #include "kronecker.hpp"
 
 #include "mfem.hpp"
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <unsupported/Eigen/KroneckerProduct>
 
 static std::vector<double> linspace(double const a, double const b, int const n) {
     std::vector<double> vec(n);
@@ -20,24 +23,33 @@ static std::vector<double> linspace(double const a, double const b, int const n)
 
 static Eigen::SparseMatrix<double> get_1d_mass_matrix_zero_initial_cond(std::vector<double> const &timepoints) {
     //assemble 1d mass matrix for piecewise linear functions
-    auto const n = timepoints.size();
-    Eigen::SparseMatrix<double> M(n - 1, n - 1);
-    std::vector<Eigen::Triplet<double>> triplets;
-    for (int j = 1; j < n - 1; ++j) {
-        int const i = j - 1;
-        auto const h = timepoints[j + 1] - timepoints[j];
-        triplets.emplace_back(i, i, h / 3.);
-        triplets.emplace_back(i, i + 1, h / 6.);
-        triplets.emplace_back(i + 1, i, h / 6.);
-        triplets.emplace_back(i + 1, i + 1, h / 3.);
-    }
-    M.setFromTriplets(triplets.begin(), triplets.end());
-    M.finalize();
-    return M;
+    auto const n = timepoints.size()-1;
+    mfem::Mesh mesh = mfem::Mesh::MakeCartesian1D(n, timepoints.back());
+    mfem::H1_FECollection fec(1, mesh.Dimension());
+    mfem::FiniteElementSpace fes(&mesh, &fec);
+    mfem::BilinearForm mass(&fes);
+    mass.AddDomainIntegrator(new mfem::MassIntegrator());
+    mass.Assemble();
+    mass.Finalize();
+    auto const M_full = vtint::mfem_to_eigen(mass.SpMat());
+    using triplet_t = Eigen::Triplet<double>;
+    auto const M_d = M_full.toDense();
+    std::vector<triplet_t> triplets;
+    for (int i = 1; i < n+1; ++i) {
+        for (int j = 1; j < n+1; ++j) {
+            if (std::abs(M_d.coeff(i, j)) > 1e-10) {
+                triplets.emplace_back(i-1, j-1, M_d.coeff(i, j));
+            }
+        }
+    }
+    Eigen::SparseMatrix<double> M0(n, n);
+    M0.setFromTriplets(triplets.begin(), triplets.end());
+    return M0;
 }
 
 namespace mfds {
-    heat_analytic_linear_operator_t::heat_analytic_linear_operator_t(mfem::FiniteElementSpace &X, const int nt, const double T, double rho_) :
+    heat_analytic_linear_operator_t::heat_analytic_linear_operator_t(mfem::FiniteElementSpace &X, const int nt,
+                                                                     const double T, double rho_) :
             rho(rho_) {
         t_points = linspace(0., T, nt + 1);
         std::shared_ptr<mfds::discrete_sine_transform_t> dst_ptr(new mfds::discrete_sine_transform_t(nt, true));
@@ -47,6 +59,8 @@ namespace mfds {
         //mass in time
         std::vector<double> const t_points = linspace(0., T, nt + 1);
         mass_t = get_1d_mass_matrix_zero_initial_cond(t_points);
+        std::cout << "mass_t" << std::endl;
+        std::cout << mass_t << std::endl;
         //get eigenvalues
         auto const values_and_vectors = get_analytic_eigenvalues_and_vectors(t_points);
         generalized_eigenvalues = std::get<0>(values_and_vectors);
@@ -67,13 +81,20 @@ namespace mfds {
     }
 
     void
-    heat_analytic_linear_operator_t::get_space_time_boundary_dofs_internal(mfem::FiniteElementSpace &X, const int nt, const double T) {
+    heat_analytic_linear_operator_t::get_space_time_boundary_dofs_internal(mfem::FiniteElementSpace &X, const int nt,
+                                                                           const double T) {
         mfem::Array<int> boundary_dofs;
         X.GetBoundaryTrueDofs(boundary_dofs);
+        std::cout << "boundary_dofs, size:" << boundary_dofs.Size() << std::endl;
+        for (int i = 0; i < boundary_dofs.Size(); ++i) {
+            std::cout << boundary_dofs[i] << std::endl;
+        }
+        std::cout << std::flush;
         mfem::Mesh mesh_t = mfem::Mesh::MakeCartesian1D(nt, T);
-        mfem::H1_FECollection T_fec(1, mesh_t.Dimension());
+        mfem::H1_FECollection T_fec(1, mesh_t.Dimension()); //we need here other boundary conditions
         mfem::FiniteElementSpace T_fes(&mesh_t, &T_fec);
-        mfem::Array<int> space_time_boundary_dofs_mfem = vtint::spatial_essential_dofs_to_spacetime_dofs(boundary_dofs, T_fes, X);
+        mfem::Array<int> space_time_boundary_dofs_mfem = vtint::spatial_essential_dofs_to_spacetime_dofs(boundary_dofs,
+                                                                                                         T_fes, X);
         for (int i = 0; i < space_time_boundary_dofs_mfem.Size(); ++i) {
             this->space_time_boundary_dofs.push_back(space_time_boundary_dofs_mfem[i]);
         }
@@ -95,16 +116,44 @@ namespace mfds {
     }
 
     vtint::block_vector_t heat_analytic_linear_operator_t::operator*(const vtint::block_vector_t &v) const {
+        if (v.get_n_blocks() != t_points.size() - 1) {
+            throw std::runtime_error("Block vector has wrong number of blocks");
+        }
+        if (generalized_eigenvalues.size() != v.get_n_blocks()) {
+            throw std::runtime_error("Generalized eigenvalues have wrong size");
+        }
         //first transform v to eigenbasis
         auto const v_tilda = apply_temporal_solver(v, dst);
+        //auto const v_tilda = apply_temporal_kronecker(v, dst);
+
         //now loop over all time-constant blocks and apply spatial operators
-        auto res = v;
+        auto res = v_tilda;
         for (size_t i = 0; i < v.get_n_blocks(); ++i) {
-            Eigen::VectorXd const v_i = v_tilda.get_block(i);
+            Eigen::VectorXd const &v_i = v_tilda.get_block(i);
             res.get_block(i) = (1. + rho * generalized_eigenvalues[i]) * mass_x * v_i + rho * stiffness_x * v_i;
         }
+
+
+        /*
+        Eigen::MatrixXd Lam(v.get_n_blocks(), v.get_n_blocks());
+        Lam.setZero();
+        for(size_t i=0; i<v.get_n_blocks(); ++i){
+            Lam(i, i) = generalized_eigenvalues[i];
+        }
+        auto const Id = Eigen::MatrixXd::Identity(v.get_n_blocks(), v.get_n_blocks());
+        //auto T1 = Eigen::KroneckerProduct(Lam, rho*mass_x);
+        //auto T2 = Eigen::KroneckerProduct(Id, rho*stiffness_x+mass_x);
+
+        auto const v1 = apply_kronecker_eigen_eigen(v_tilda, Lam, rho*mass_x);
+        auto const v2 = apply_kronecker_eigen_eigen(v_tilda, Id, rho*stiffness_x+mass_x);
+
+        auto res = v1+v2;
+
+        */
+
         //now transform back
         res = apply_temporal_kronecker(res, dst);
+        //res = apply_temporal_solver(res, dst);
         //and apply mass again
         res = apply_temporal_kronecker(res, mass_t);
 
@@ -122,11 +171,6 @@ namespace mfds {
         return res.to_vector();
     }
 
-    double heat_analytic_linear_operator_t::operator()(size_t i, size_t j) const {
-        //reverse kronecker index and get it via multiplication
-
-        return 1.;
-    }
 
     Eigen::VectorXd const &heat_analytic_linear_operator_t::get_diagonal_preconditioner() const {
         return mass_diagonal;

mfds/active_set_operator.hpp

@@ -142,8 +142,6 @@ namespace mfds {
 
         std::vector<int> get_space_time_boundary_dofs() const;
 
-        double operator()(size_t i, size_t j) const;
-
         Eigen::VectorXd const &get_diagonal_preconditioner() const;
 
 

mfds/heat_control.cpp

@@ -207,9 +207,6 @@ namespace mfds {
             return target_fun(p, t);
         };
         vtint::time_dependent_gridfunction_t u_bdr_gf(T_test, X, target_reverse_arguments);
-        vtint::block_vector_t u_bdr_transformed = !use_fft ? mfds::apply_temporal_kronecker(u_bdr_gf.get_blockvector(),
-                                                                                             S_inv)
-                : mfds::apply_temporal_solver(u_bdr_gf.get_blockvector(), dst);
         toc = std::chrono::high_resolution_clock::now();
         timings.add_timing("transform rhs to temporal eigenbasis (boundary data)",
                            std::chrono::duration_cast<std::chrono::nanoseconds>(toc - tic1));
@@ -225,7 +222,8 @@ namespace mfds {
         for (int i = 0; i < eigenvalues.size(); ++i) {
             double const lambda = eigenvalues[i];
             mfem::Vector const cur_rhs = vtint::eigen_to_mfem(rhs_transformed.get_block(i));
-            mfem::Vector const cur_u_brd_transformed = vtint::eigen_to_mfem(u_bdr_transformed.get_block(i));
+            mfem::Vector cur_u_brd_transformed = cur_rhs;
+            cur_u_brd_transformed *= 0.;
             auto const [sol, cg_iter] = solve_spatial_problem_heat_serial(X, cur_rhs, lambda, rho, cur_u_brd_transformed,
                                                                           false);
             state_transformed.get_block(i) = vtint::mfem_to_eigen(sol);

mfds/kronecker.hpp

@@ -59,6 +59,22 @@ namespace mfds {
         return res;
     }
 
+    template <class Eigen_OP_T, class Eigen_OP_T2>
+    vtint::block_vector_t apply_kronecker_eigen_eigen(vtint::block_vector_t const& v, Eigen_OP_T& B, Eigen_OP_T2& A) {
+        //calculate the kronecker product (B \otimes A) * v
+        Eigen::MatrixXd V = to_matrix(v);
+        Eigen::MatrixXd BV = B*V.transpose();
+        Eigen::MatrixXd VBT = BV.transpose();
+        vtint::block_vector_t res(VBT.cols());
+        for(Eigen::Index i=0; i<VBT.cols(); ++i) {
+            //get i-th column of VBT
+            Eigen::VectorXd v_i = VBT.col(i);
+            //write back to result
+            res.get_block(i) = A*v_i;
+        }
+        return res;
+    }
+
 
 } // mfds
 

mfds/tests/CMakeLists.txt

 project(mfds_test)
 
 add_executable(mfds_test_runner kronecker_test.cpp
-        analytical_eigenfuctions_test.cpp active_set_wrapper_test.cpp)
+        analytical_eigenfuctions_test.cpp boundary_condition_wrapper.cpp)
 target_link_libraries(mfds_test_runner PRIVATE Catch2::Catch2WithMain Eigen3::Eigen vtint mfds)
 target_include_directories(mfds_test_runner
         INTERFACE ${CMAKE_CURRENT_SOURCE_DIR})

mfds/tests/boundary_condition_wrapper.cpp

+//
+// Created by mreic on 12.08.2024.
+//
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include <Eigen/Dense>
+#include "mfem.hpp"
+
+#include "active_set_operator.hpp"
+
+class matrix_wrapper_t{
+public:
+    matrix_wrapper_t(mfem::SparseMatrix const matrix_){
+        matrix = matrix_;
+    }
+    Eigen::VectorXd operator*(Eigen::VectorXd const& v){
+        Eigen::VectorXd result(matrix.Height());
+        mfem::Vector mfem_v(v.size());
+        for(int i=0; i<v.size(); ++i){
+            mfem_v[i] = v[i];
+        }
+        mfem::Vector mfem_result(result.size());
+        matrix.Mult(mfem_v, mfem_result);
+        for(int i=0; i<result.size(); ++i){
+            result[i] = mfem_result[i];
+        }
+        return result;
+    }
+
+    Eigen::Index rows() const { return matrix.Height(); }
+    Eigen::Index cols() const { return matrix.Width(); }
+private:
+    mfem::SparseMatrix matrix;
+};
+
+template <class BW>
+bool check_matrix_are_same(BW& boundary_wrapper, mfem::SparseMatrix& matrix){
+    Eigen::VectorXd v = Eigen::VectorXd::Random(matrix.Width());
+    Eigen::VectorXd res1 = boundary_wrapper.mult_internal(v);
+    Eigen::VectorXd res2 = matrix_wrapper_t(matrix)*v;
+    return res1.isApprox(res2);
+}
+
+template <class Op>
+Eigen::MatrixXd to_matrix_brute_force(Op& op){
+    Eigen::MatrixXd mat = Eigen::MatrixXd::Zero(op.rows(), op.cols());
+    for(Eigen::Index i=0; i<op.rows(); ++i){
+        Eigen::VectorXd e = Eigen::VectorXd::Zero(op.rows());
+        e(i) = 1;
+        mat.col(i) = op*e;
+    }
+    return mat;
+}
+
+TEST_CASE("Boundary Condition Wrapper","[mfds][boundary_condition][laplace]"){
+    int const nx = 3;
+    mfem::Mesh mesh = mfem::Mesh::MakeCartesian2D(nx, nx, mfem::Element::TRIANGLE, true);
+    mfem::H1_FECollection fec(1, mesh.Dimension());
+    mfem::FiniteElementSpace fes(&mesh, &fec);
+    mfem::BilinearForm k(&fes);
+    k.AddDomainIntegrator(new mfem::DiffusionIntegrator());
+    k.Assemble();
+    mfem::Array<int> boundary_dofs_m;
+    fes.GetBoundaryTrueDofs(boundary_dofs_m);
+    std::vector<int> boundary_dofs(boundary_dofs_m.Size());
+    for(int i=0; i<boundary_dofs_m.Size(); ++i){
+        boundary_dofs[i] = boundary_dofs_m[i];
+    }
+    mfem::Vector rhs_m(fes.GetVSize()), rhs_mh(fes.GetVSize());
+    rhs_m = 1.0;
+    mfem::SparseMatrix Km(k.SpMat());
+    matrix_wrapper_t Kw(Km);
+    mfem::Vector u_m(rhs_m), u_mh(rhs_m);
+    u_m = 0.0;
+
+    //homogenize system
+    mfem::SparseMatrix Kmh;
+    k.FormLinearSystem(boundary_dofs_m, u_m, rhs_m, Kmh, u_mh, rhs_mh);
+
+    //Kmh.ToDenseMatrix()->Print();
+
+    //construct wrapper
+    mfds::boundary_condition_wrapper_t bc_wrapper(Kw, boundary_dofs);
+
+    //get the matrices as dense matrices
+    auto Kmhw = matrix_wrapper_t(Kmh);
+    Eigen::MatrixXd Kmh_dense = to_matrix_brute_force(Kmhw);
+    Eigen::MatrixXd Kh_dense = to_matrix_brute_force(bc_wrapper);
+
+    bool matrices_same = check_matrix_are_same(bc_wrapper, Kmh);
+    REQUIRE(!matrices_same); //matrices cannot be the same due to it keeping symmetry
+    Eigen::VectorXd rhs(rhs_m.Size());
+    rhs.setOnes();
+    Eigen::VectorXd u(rhs_m.Size());
+    u.setZero();
+    auto const rhs_h = bc_wrapper.modify_rhs(rhs, u);
+    REQUIRE(rhs_h.size() == rhs_mh.Size());
+    for(int i=0; i<rhs_h.size(); ++i){
+        REQUIRE_THAT(rhs_h(i), Catch::Matchers::WithinAbs(rhs_mh(i), 1e-10));
+    }
+
+    //check if multiplication yields the same
+    auto const v = bc_wrapper.mult_internal(rhs_h);
+    auto const v_m = Kmhw*rhs_h;
+
+    //check if the result is correct
+    REQUIRE(v.size() == v_m.size());
+    REQUIRE(v.isApprox(v_m));
+
+}
+

playground/control_analytic_performance_and_eoc_study.cpp

@@ -107,8 +107,8 @@ int main(int argc, char* argv[]) {
     int nx0 = 2;
     double const T = 1.;
     int max_refinements = 4;
-    int parabolic_scaling = 1;
-    target_t target_fun_type = target_t::moving_target;
+    int parabolic_scaling = 0;
+    target_t target_fun_type = target_t::H22;
     int target_fun_int = static_cast<int>(target_fun_type);
     const char * output_path = "./";
     int spatial_dim = 2;

playground/example_matrix_free_operator_2d.cpp

@@ -7,6 +7,7 @@
 
 #include "mfem.hpp"
 #include "active_set_operator.hpp"
+#include "heat_control.h"
 #include "target.hpp"
 #include "conjugate_gradient.hpp"
 #include "homogenization.hpp"
@@ -17,7 +18,11 @@ static double target_function(mfem::Vector const &p, double const t) {
     double const x = p(0);
     double const y = p(1);
     //return 16*t * t * x * (1 - x) * y * (1 - y);
-    return std::sin(M_PI * x) * std::sin(M_PI * y) * std::sin(M_PI * t);
+    double const offset = 0.;
+    if(t<offset){
+        return 0.;
+    }
+    return std::sin(M_PI * x) * std::sin(M_PI * y) * std::sin(M_PI * (t-offset));
 }
 
 static double lower_bound_fun(double const t, mfem::Vector const &p) {
@@ -29,7 +34,7 @@ static double upper_bound_fun(double const t, mfem::Vector const &p) {
 }
 
 int main(int argc, char *argv[]) {
-    int const refinement = 6;
+    int const refinement = 5;
     int const nx = 2 * std::pow(2, refinement);
     int const ny = 2 * std::pow(2, refinement);
     int const nt = 2 * std::pow(2, refinement);
@@ -39,7 +44,7 @@ int main(int argc, char *argv[]) {
     double const newton_relax_param = 1;
 
 
-    mfem::Mesh spatial_mesh = mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::QUADRILATERAL, true, 1.0, 1.0);
+    mfem::Mesh spatial_mesh = mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::TRIANGLE, true, 1.0, 1.0);
     int const dimX = spatial_mesh.Dimension();
     mfem::Mesh temporal_mesh = mfem::Mesh::MakeCartesian1D(nt, T);
 
@@ -50,13 +55,15 @@ int main(int argc, char *argv[]) {
     mfem::FiniteElementSpace T_test(&temporal_mesh, &T_fec);
 
     //construct matrix free operator
-    double const rho = std::pow(std::min(1. / nx, 1. / ny), 2);
+    double const rho = 0*std::pow(std::min(1. / nx, 1. / ny), 2);
     mfds::heat_analytic_linear_operator_t heat_operator(X, nt, T, rho);
 
+
     //get space time boundary
     auto const space_time_boundary_dofs = heat_operator.get_space_time_boundary_dofs();
+    std::cout << "n_stb: " << space_time_boundary_dofs.size() << std::endl;
     std::vector<bool> isBoundary(X.GetTrueVSize() * nt, false);
-    std::cout << "Boundary dofs:" << std::endl;
+    std::cout << "Space-Time Boundary dofs:" << std::endl;
     for (auto const dof: space_time_boundary_dofs) {
         std::cout << dof << std::endl;
         isBoundary[dof] = true;
@@ -72,6 +79,7 @@ int main(int argc, char *argv[]) {
     auto const rhs_bdr = boundary_condition_wrapper.modify_rhs(rhs_inhom.to_vector(), u_bdr);
 
 
+    auto test = heat_operator*rhs_bdr;
 
     Eigen::ConjugateGradient<mfds::boundary_condition_wrapper_t<mfds::heat_analytic_linear_operator_t>,
             Eigen::Lower | Eigen::Upper, Eigen::IdentityPreconditioner> cg;
@@ -80,6 +88,9 @@ int main(int argc, char *argv[]) {
     std::cout << "CG:       #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
 
 
+    auto const [sim_res, sim_timing] = mfds::solve_heat_analytic(X, T_test, target_function, rho, true);
+    auto const state_reference = sim_res.state_full;
+
     /*
 
     //wrap operator so that boundary conditions can be applied
@@ -114,6 +125,10 @@ int main(int argc, char *argv[]) {
     u0_vec = 0.;
     auto const state_full = vtint::add_initial_condition(state_block, u0_vec);
 
+    //calculate difference to reference:
+    double const diff_to_ref = (state_reference.to_vector()-state_full.to_vector()).lpNorm<Eigen::Infinity>();
+    std::cout << "max difference to reference in infinity norm: " << diff_to_ref << std::endl;
+
     //compute L2 error, to be sure everything is correct
     vtint::time_dependent_gridfunction_t state_gf(T_test, X, state_full);
     //state_gf.save_as_vtk("state", 0.01);
@@ -127,6 +142,9 @@ int main(int argc, char *argv[]) {
     //save function to paraview
     state_gf.save_as_vtk("state", 0.01);
 
+    vtint::time_dependent_gridfunction_t reference_state_gf(T_test, X, state_reference);
+    reference_state_gf.save_as_vtk("reference_state", 0.01);
+
     vtint::time_dependent_gridfunction_t target_gf(T_test, X, target_fun_reverse_arguments);
     target_gf.save_as_vtk("target", 0.01);
 

playground/mixed_iga_mfem.cpp

+//
+// Created by mreic on 21.05.2024.
+//
+
+#include <iostream>
+#include "mfem.hpp"
+
+constexpr int const n_refinements = 1;
+constexpr int const order_z = 0;
+
+double u_analytic(mfem::Vector const& p){
+    double const x = p(0);
+    double const y = p(1);
+    return std::sin(M_PI*x)*std::sin(M_PI*y);
+}
+
+
+//standard main with arguments
+int main(int argc, char **argv)
+{
+    std::string const mesh_path = MFDS_MESHES;
+    std::string const mesh_file = "square_c1.mesh";
+
+    //read mesh
+    mfem::Mesh mesh(mesh_path + mesh_file);
+
+    //uniform refinement
+    for (int i = 0; i < n_refinements; ++i) {
+        mesh.UniformRefinement();
+    }
+
+    //this below works as expected and provides a C^1 nurbs mesh
+    /*
+    //save mesh to file
+    std::string const mesh_out = "square_c1_refined.mesh";
+    mesh.Save(mesh_out.c_str());
+
+    //display mesh with GLVis
+    mfem::socketstream vis_stream("localhost", 19916);
+    if (!vis_stream.good()) {
+        std::cerr << "Cannot connect to GLVis server; exiting." << std::endl;
+        return 1;
+    }
+    vis_stream.precision(8);
+    vis_stream << "mesh\n" << mesh << "pause\n";
+     */
+
+    //setup nurbs fe space
+    if(!mesh.GetNodes()){
+        throw std::runtime_error("Mesh is no IGA mesh.");
+    }
+    mfem::FiniteElementCollection *fec = mesh.GetNodes()->OwnFEC();
+    mfem::FiniteElementSpace X(&mesh, fec);
+    mfem::L2_FECollection fecZ(order_z, mesh.Dimension());
+    //mfem::FiniteElementSpace Z(&mesh, &fecZ);
+    mfem::FiniteElementSpace Z(&mesh, fec);
+
+    //setup mass on Z
+    mfem::BilinearForm Mf(&Z);
+    Mf.AddDomainIntegrator(new mfem::MassIntegrator);
+    Mf.Assemble();
+    mfem::SparseMatrix M(Mf.SpMat());
+    M.Finalize();
+
+    //assemble laplacian matrix
+    //mfem::MixedBilinearForm Af(&Z, &X);
+    mfem::BilinearForm Af(&X);
+    mfem::ConstantCoefficient one(1.0);
+    Af.AddDomainIntegrator(new mfem::LaplaceIntegrator(one));
+    Af.Assemble();
+    mfem::SparseMatrix A(Af.SpMat());
+    A.Finalize();
+
+    std::cout << "A" << std::endl;
+    A.Print();
+    std::cout << "M" << std::endl;
+    M.Print();
+
+    //project u_analytic to X
+    mfem::GridFunction u_analytic_gf(&X);
+    mfem::FunctionCoefficient u_analytic_fc(u_analytic);
+    u_analytic_gf.ProjectCoefficient(u_analytic_fc);
+
+    //setup rhs
+    mfem::GridFunction rhs_gf(&Z);
+    mfem::Vector rhs(Z.GetVSize());
+    A.Mult(u_analytic_gf, rhs);
+
+    //solve Mz = rhs
+    mfem::CGSolver M_solver;
+    M_solver.SetOperator(M);
+    mfem::Vector z(rhs.Size());
+    M_solver.Mult(rhs, z);
+
+    //setup a grifunction for z
+    mfem::GridFunction z_gf(&Z);
+    z_gf = z;
+
+
+    //plot solution
+    mfem::socketstream vis_stream("localhost", 19916);
+    if (!vis_stream.good()) {
+        std::cerr << "Cannot connect to GLVis server; exiting." << std::endl;
+        return 1;
+    }
+    vis_stream.precision(8);
+    vis_stream << "solution\n" << mesh << z_gf << "pause\n";
+
+    return 0;
+}
\ No newline at end of file

vtint @ 5700b842

-Subproject commit a4e388b426af9c23a384d9d60635bd5c67645614
+Subproject commit 5700b842cd915f197e2227b2c683c7dbdb55546d