src/femaxxdriver.cpp

Go to the documentation of this file.

//
// C++ Implementation: femaxxdriver
//
// Description:
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <string>
#include <cassert>
#include <unistd.h>

#include "myrlog.h"

#include <ml_include.h>
#include <Epetra_Operator.h>
#include <Epetra_MultiVector.h>
#include <Epetra_SerialDenseVector.h>
#include <Epetra_Vector.h>
#include <ml_epetra_preconditioner.h>

#include "OutputCapturer.h"
#include "pbedefs.h"

#include "looseoctree.h"
#include "debugmeshbuilder.h"
#include "tetmeshbuilder.h"
#include "broadcastmeshbuilder.h"
#include "netgenreader.h"
#include "femaxmesh.h"
#include "nedelecmesh.h"
#include "lagrangemesh.h"
#include "export_Epetra_CrsMatrix.h"
#include "operatortest.h"
#include "materials.h"

#include "AbstractEigsolvOperators.h"
#include "LinearEigsolvOperators.h"
#include "QuadraticEigsolvOperators.h"
#include "Projector.h"

#include "ModalAnalysisSolver.h"
#include "femaxxdriver.h"

#include "JDBSYM.h"
#include "KnyazevLOBPCG.h"
#include "LOBPCG.h"

#include "CheckingTools.h"
#include "vtkexport.h"
#include "h5_tools.h"

using namespace mesh;

FemaxxDriver::FemaxxDriver(const Epetra_Comm& comm, Teuchos::ParameterList& params_)
    : comm_(comm),
      params_(params_),
      Q_(0),
      lambda_(0),
      nof_converged_(0),
      tmesh_(0),
      femax_mesh_(0)
{
    // Init PBE
    //
    // Timers used by FemaxxDriver and sub-components:

    // Main program parts:
    // 50   Life of FemaxxDriver
    // 52   Generation of tetrahedral mesh
    // 53   DOF mappings (boundary conditions)
    // 54   Eigenvalue solver

    // Matrix construction:
    // 60   Matrix assembly A and M
    // 61   Matrix construction: Y
    // 62   Matrix multiplication: C and C^T
    // 63   Matrix multiplication: H

    // Preconditioner/Projector construction:
    // 70   Preconditioner construction of A-sigma*M
    // 71   Preconditioner construction of H
    // 72   Solver construction of H
    // 75   Nullspace projector construction

    // Operator application:
    // 80   Application of A
    // 81   Application of M
    // 82   Application of Preconditioner for A-sigma*M
    // 83   Application of Preconditioner for H
    // 84   Application of Nullspace Projector

    // Hierarchical precinditioner:
    // 110  Construction
    // 111  Operator
    // 112  Apply
    // 113  (1,1)-solve
    // 114  (2,2)-solve
    // 115  off-diagonal multiply

    // Utility routines:
    // 130  MatrixMarket export
    // 131  Octree construction
    // 132  Octree traversal
    // 133  Quality factors
    // 134  VTK export

    const int nof_pbe_timers = 150;
    pbe_init("femaxx", nof_pbe_timers, 0, NULL);
    pbe_start(50, "Whole program");
}

FemaxxDriver::~FemaxxDriver() {
    delete Q_;
    delete[] lambda_;
    delete femax_mesh_;
    delete tmesh_;

    // Finalize PBE
    pbe_stop(50);

    OutputCapturer cap;
    cap.begin_capture();
    pbe_dump();
    pbe_finalize(0);
    cap.end_capture();
    if (comm_.MyPID() == 0) {
        cap.log("PBE output", RLOG_CHANNEL("info"));
    }
}

void FemaxxDriver::set_defaults(Teuchos::ParameterList& params) {
    // Defaults for cmd line parser
    params.set("element_order", 1);
    params.set("kmax", 5);
    params.set("eigtol", 1e-8);
    params.set("tau", 0.0);
    params.set("sym_plane_config", 0);
    params.set("eigsolver", "jdsym");
    params.sublist("asigma_precon").set("type", "ml");
    params.sublist("asigma_precon").set("solver11", "lu");
    params.sublist("asigma_precon").set("solver22", "diag");
    params.sublist("h_solver").set("type", "pcg");
    params.sublist("h_precon").set("type", "ml");
    params.sublist("h_precon").set("solver11", "lu");
    params.sublist("h_precon").set("solver22", "diag");
    params.set("vtk_file_name", "");
    params.set("export_mtx", false);
    params.set("debug", false);
    params.set("disable_log", false);
        
    // default values for ML preconditioner of H for smoothed aggregation
    ML_Epetra::SetDefaults("SA", params.sublist("h_precon").sublist("ml_params"));
    // replace Gauss-Seidel smoother by MLS
    params.sublist("h_precon").sublist("ml_params").set("smoother: type", "MLS");

    // Amesos params used if h_solver.type == "lu"
    params.sublist("h_solver").sublist("amesos_params").set("solver_type", "Superludist");
    params.sublist("h_solver").sublist("amesos_params").set("MatrixType", "SPD");
    params.sublist("h_solver").sublist("amesos_params").set("OutputLevel", 2);

    // Amesos params used if asigma_precon.type == "lu"
    params.sublist("asigma_precon").sublist("amesos_params").set("solver_type", "Superludist");
    params.sublist("asigma_precon").sublist("amesos_params").set("MatrixType", "symmetric");
    params.sublist("asigma_precon").sublist("amesos_params").set("OutputLevel", 2);

    // ML params for (1,1)-block of A-sigma*M (entire matrix if element_order == 1)
    ML_Epetra::SetDefaults("maxwell", params.sublist("asigma_precon").sublist("ml_params"));

    // ML 4.0
    params.sublist("asigma_precon").sublist("ml_params").set("negative conductivity", true);
    params.sublist("asigma_precon").sublist("ml_params").set("subsmoother: type", "symmetric Gauss-Seidel");
#   if 0
    // FIXME: The repartitioning feature in ML crashes as of Trilinos 6.0.11 -> commented out.
    params.sublist("asigma_precon").sublist("ml_params").set("repartition: enable",true);
    //how strict the nodal load-balancing is
    params.sublist("asigma_precon").sublist("ml_params").set("repartition: node max min ratio",1.1);
    //minimum number of nodal unknowns per processor, you might have to change
    params.sublist("asigma_precon").sublist("ml_params").set("repartition: node min per proc",200);
    //how strict the edge load-balancing is
    params.sublist("asigma_precon").sublist("ml_params").set("repartition: edge max min ratio",1.1);
    //minimum number of edge unknowns per processor, you might have to change 
    params.sublist("asigma_precon").sublist("ml_params").set("repartition: edge min per proc",600);
    params.sublist("asigma_precon").sublist("ml_params").set("repartition: partitioner","ParMETIS");
#   endif
    //params.sublist("asigma_precon").sublist("ml_params").set("aggregation: type", "MIS");
    //params.sublist("asigma_precon").sublist("ml_params").set("coarse: max size", 30);
    //params.sublist("asigma_precon").sublist("ml_params").set("aggregation: threshold", 0.0);    

    // JDBSYM params
    JDBSYM::set_defaults(params.sublist("jd_params"), params.get<int>("kmax"));
    params.sublist("jd_params").set("clvl", 2);
}

void FemaxxDriver::load_mesh() {
    log_mem_footprint(comm_, "initial");
    rInfo("Read mesh file...");
    if (access(params_.get<string>("mesh_file_name").c_str(), R_OK)) {
        rError("Unable to access mesh file \"%s\".", params_.get<string>("mesh_file_name").c_str());
        exit(1);
    }

    //ifstream istr(params_.get<string>("mesh_file_name").c_str());
    const char* meshFileName = params_.get<string>("mesh_file_name").c_str();
    const char* materialFileName = params_.get<string>("material_file_name").c_str();

    //
    // Generate mesh data from file input on all processors
    //
    log_mem_footprint(comm_, "before loading mesh");
    pbe_start(52, "Mesh generation");
    TetMeshBuilder* builder = new TetMeshBuilder(dynamic_cast<const Epetra_MpiComm&>(comm_).GetMpiComm());
// #ifdef HAVE_MPI
//     BroadcastMeshBuilder* broadcast = new BroadcastMeshBuilder(&comm_, builder);
//     NetgenReader* reader(0);
//     if (comm_.MyPID() == 0) {
//         reader = new NetgenReader(&istr, broadcast);
//         reader->read();
//     } else
//         broadcast->receiver();
//     delete broadcast;
//     if (comm_.MyPID() == 0)
//         delete reader;
    HDF5ParallelReader* reader = new HDF5ParallelReader(meshFileName, builder, dynamic_cast<const Epetra_MpiComm&>(comm_).GetMpiComm());
    reader->read();
    delete reader;

// #else
//     NetgenReader* reader = new NetgenReader(&istr, builder);
//     reader->read();
//     delete reader;
// #endif

    tmesh_ = builder->get_mesh();
    delete builder;
    pbe_stop(52);
    if (comm_.MyPID() == 0)
        tmesh_->log_mesh_info();
    log_mem_footprint(comm_, "after loading tetmesh");

    // Load material properties (if any)
    tmesh_->load_materials(materialFileName);

    // create FE mesh data structures
    int sym_plane_config = params_.get<int>("sym_plane_config");
    rInfo("Calculate DOF mappings (sym_plane_config=%d)...", sym_plane_config);
    pbe_start(53, "Dof mappings");
    femax_mesh_ = new FemaxMesh(*tmesh_, params_.get<int>("element_order"), sym_plane_config);
    pbe_stop(53);
    log_mem_footprint(comm_, "after calc dof mappings");
}

void FemaxxDriver::calculate_eigenfields() {

    // create operators
    AbstractEigsolvOperators* lop;
    if (params_.get<int>("element_order") == 1)
        lop = new LinearEigsolvOperators(*femax_mesh_, params_, comm_);
    else
        lop = new QuadraticEigsolvOperators(*femax_mesh_, params_, comm_);
    log_mem_footprint(comm_);

    if (params_.get<bool>("debug")) {
        OperatorTest tester(lop->getA());
        tester.dump_stats(cout, "A original");
        tester.set_operator(lop->getM());
        tester.dump_stats(cout, "M original");
        tester.set_operator(lop->getH());
        tester.dump_stats(cout, "H original");
    }

    // Free mesh data structures
    // FIXME: these data structures could freed even earlier (after construction of A, M, Y),
    // but it would make the program ugly.
    log_mem_footprint(comm_, "before mesh data structure dealloc");
//    delete femax_mesh_; /** Do not delete femax_msh here because it will be used in the evaluation of the solution */
    
//    delete tmesh_;
    
//    femax_mesh_ = 0;
//    tmesh_ = 0;
    log_mem_footprint(comm_, "after mesh data structure dealloc");

    const Epetra_Operator* A = lop->getA();
    const Epetra_Operator* M = lop->getM();
    // Here the preconditioner is actually constructed
    if (params_.get<double>("sigma") == 0.0)
        rWarning("Shift sigma is zero: Calculating preconditioner for a singular matrix.");
    const Epetra_Operator* Prec = lop->getAsigmaPrec(params_.get<double>("sigma"));
    Epetra_Operator* H = lop->getH();
    Epetra_Operator* Y = lop->getY();
    Epetra_Operator* C = lop->getC();
    log_mem_footprint(comm_, "after constr of precon for A-sigma*M");

    if (params_.get<bool>("debug")) {
        OperatorTest tester(A);
        tester.dump_stats(cout, "A reordered");
        tester.set_operator(M);
        tester.dump_stats(cout, "M reordered");
        tester.set_operator(H);
        tester.dump_stats(cout, "H reordered");
    }

    //
    // Eigenvalue solver
    //
    int myPid = comm_.MyPID();
    int numProc = comm_.NumProc();

    // Initialize the random number generator
    srand((unsigned) time(0));

    // Initialize memory counters
    initMemCounters();

    double highMem = currentSize();
    highMem = (highMem < currentSize()) ? currentSize() : highMem;

    int globalSize = (A->OperatorRangeMap()).NumGlobalElements();

    // Get the first eigenvalue for the mass matrix
    int kmax = params_.get<int>("kmax");
    kmax = (kmax > globalSize) ? globalSize : kmax;

    // Define arrays for the eigensolver algorithm
    lambda_ = new (nothrow) double[kmax];
    assert(lambda_ != 0);

    Q_ = new Epetra_MultiVector(A->OperatorRangeMap(), kmax, false);
    Q_->Random();
    log_mem_footprint(comm_, "after alloc of Q and lambda");

    highMem = (highMem > currentSize()) ? highMem : currentSize();

    // Projector construction (including construction of solver for H)
    pbe_start(75, "Nullspace projector construction");
    rInfo("Calculate nullspace projector...");
    Projector YHC(Y, lop->getHSolver(), C);
    pbe_stop(75);
    log_mem_footprint(comm_, "after constr of projector");

    if (params_.get<bool>("debug")) {
        OperatorTest tester(&YHC);
        tester.dump_stats(cout, "Nullspace projector");
        // tester.set_operator(Prec);
        // tester.dump_stats(cout, "A-sigma*M operator");
        tester.set_operator(lop->getHSolver(), true);
        tester.dump_stats(cout, "H solver");
        tester.set_operator(Prec, true);
        tester.dump_stats(cout, "A-sigma*M precond");
    }

    // Define the GeneralEigenSolver object
    pbe_start(54, "Eigenvalue solver");
    double eigtol = params_.get<double>("eigtol");
    ModalAnalysisSolver *mySolver = 0;

    if (params_.get<string>("eigsolver") == "jdsym") {
        params_.sublist("jd_params").set("tol", eigtol);
        mySolver = new JDBSYM(comm_, A, M, Prec, &YHC, params_.sublist("jd_params"));
    } else if (params_.get<string>("eigsolver") == "knyazev") {
        const int max_iter = 300;
        const int clvl = 3;
        YHC.Apply(*Q_, *Q_); //applying projector step on blocked set of initial random vectors.
        mySolver = new KnyazevLOBPCG(comm_, A, M, Prec, &YHC, eigtol, max_iter, clvl); 
    } else if (params_.get<string>("eigsolver") == "lobpcg") {
        // FIXME: workspace mgmt of LOBPCG must be corrected to allow for any kmax
        const int max_iter = 300;
        const int clvl = 3;
        YHC.Apply(*Q_, *Q_);
        kmax = 2;
        mySolver = new LOBPCG(comm_, A, M, Prec, &YHC, kmax, eigtol, max_iter, clvl);
        //hack: 2 eigenvectors are computed. The three other Qs are used as working space
        //hack works only if computed eigenvectors are at least 5.
        //number of eigenvectors * 2 + 1 = needed size of parameter Q
        //if you want to compute x eigenvectors for this version of LOBPCG
        //change kmax=x above accordingly and run script
        //mpirun ... -kmax=2*x+1 ...
    } else {
        rError("Unknown eigensolver \"%s\".", params_.get<string>("eigsolver").c_str());
    }

    // Solve the eigenproblem
    nof_converged_ = mySolver->solve(kmax, *Q_, lambda_);
    pbe_stop(54);

    // Converged eigenpairs
    Epetra_MultiVector Q = get_eigenvectors();
    Epetra_SerialDenseVector lambda = get_eigenvalues();

    // Output information on simulation
    if (nof_converged_ < 0) {
        rError("Eigensolver exited with error %d!!!", nof_converged_);
    } else {
        rInfo("Number of computed eigenmodes: %d", nof_converged_);

        if (nof_converged_ > 0) {
            // Check the computed eigenmodes
            CheckingTools myVerify(comm_);
            // Check orthonormality of eigenvectors
            double tmp = myVerify.errorOrthonormality(&Q, M);
            if (myPid == 0)
                rInfo("Maximum coefficient in matrix Q^T M Q - I = %le", tmp);
            // Print out norm of residuals
            myVerify.errorEigenResiduals(Q, lambda_, A, M);
#if 0
            // Check C^T Q = 0
            Epetra_MultiVector CtQ(C->OperatorDomainMap(), nof_converged_);
            (const_cast<Epetra_Operator*>(C))->SetUseTranspose(true);
            C->Apply(Q, CtQ);

            // FIXME: Calling the Epetra_MultiVector::NormInf() leads to SIGSEGV for
            // small problems. This is probably a bug in Trilinos

            // Get the largest M-norm of Y from H
            double maxNormY = 0.0;
            Epetra_Vector Hdiag(dynamic_cast<Epetra_RowMatrix*>(H)->RowMatrixRowMap());
            dynamic_cast<Epetra_RowMatrix*>(H)->ExtractDiagonalCopy(Hdiag);
            Hdiag.NormInf(&maxNormY);
            double maxError = 0.0;
            for (int i = 0; i < nof_converged_; ++i) {
                double infNorm = 0.0;
                CtQ(i)->NormInf(&infNorm);
                maxError = (infNorm > maxError) ? infNorm : maxError;
            }
            if (myPid == 0) {
                cout << endl;
                cout << " Maximum M-norm of vectors Y = " << maxNormY << endl;
                cout << " Maximum coefficient in matrix C^T Q = " << maxError << endl;
                cout << endl;
            }
#endif
        }

        // Memory
        Teuchos::ParameterList& jd_params(params_.sublist("jd_params"));
        if (params_.get<string>("eigsolver") == "jdsym") {
            double maxHighMem = 0.0;
            comm_.MaxAll(&highMem, &maxHighMem, 1);
            rDebug("Memory: High water mark in set-up (EST)               : %.2lf MB", maxHighMem);
            rDebug("Memory requested per processor for eigenvectors  (EST): %.2lf MB",
                Q_->GlobalLength()*kmax*sizeof(double)/1024.0/1024.0/numProc);
            rDebug("Memory requested per processor for working space (EST): %.2lf MB",
                Q_->GlobalLength()*(2*kmax + jd_params.get<int>("jmax"))*sizeof(double)/1024.0/1024.0/numProc);
        }
        mySolver->memoryInfo();
        mySolver->operationInfo();
        mySolver->timeInfo();

    } // if (nof_converged_ < 0)

    // Release all objects
    log_mem_footprint(comm_, "before dealloc eigenvalue solver");
    delete mySolver; mySolver = 0;
    log_mem_footprint(comm_, "after dealloc eigenvalue solver");

    delete lop; lop = 0;
    log_mem_footprint(comm_, "after dealloc operators");

#if 0    
    // Export to hdf5
    if (myPid == 0) {
        h5_create_empty_file("eigenmodes.h5");
        h5_write_param_list("eigenmodes.h5", params_);
    }
    comm_.Barrier();
    h5_write_eigenmodes("eigenmodes.h5", Q, lambda);
#endif
}

void FemaxxDriver::postprocess() {
    rAssert(tmesh_);
    rAssert(femax_mesh_);

    Epetra_MultiVector Q = get_eigenvectors();
    Epetra_SerialDenseVector lambda = get_eigenvalues();

    pbe_start(131, "Construction of octree");
    if (comm_.MyPID() == 0) {
        tmesh_->construct_octree();
        // tmesh_->get_octree()->export_vtk();
    }

    pbe_stop(131);

    //
    // Q-factors
    //
    pbe_start(133, "Computation of quality factors");
    vector<double> q_factors;
    for (int k = 0; k < nof_converged_; ++ k) {
        //
        // Bring k-th eigenvector to processor 0
        //
        Epetra_MultiVector q_dist(View, Q, k, 1);

        // create a target map, for which all the elements are on proc 0
        int NumMyElements_target = 0;
        if (comm_.MyPID() == 0)
            NumMyElements_target = Q.GlobalLength();
        Epetra_Map TargetMap(-1, NumMyElements_target, 0, comm_);

        // redistribute vector to proc 0
        Epetra_Vector q(TargetMap);
        q.Export(q_dist, Epetra_Export(q_dist.Map(), TargetMap), Add);

        // Compute q-factor on processor 0
        if (comm_.MyPID() == 0) {
            q_factors.push_back(femax_mesh_->get_nedelec_mesh().q_factor(lambda[k], q.Values()));
        }
    }
    pbe_stop(133);

    //
    // Print eigenvalue summary
    //
    if (comm_.MyPID() == 0) {
        // pi
        const double pi = ::acos(-1.0);
        // speed of light [m / s]
        const double cLight = 2.99792458e8;

        ostringstream buf;
        buf << "Eigenvalue summary:" << endl;
        buf << "nr.   norm. frequency       freq. [MHz]            Q";
        for (int k = 0; k < nof_converged_; ++ k) {
            double omega = ::sqrt(lambda[k]) * cLight;
            double f = omega / 2.0 / pi;
            buf << endl << right
                << setw(3) << k+1 << " "
                << setw(17) << setprecision(8) << fixed << lambda[k] << " "
                << setw(17) << setprecision(8) << fixed << f / 1e6 << " "
                //                << setw(17) << setprecision(8) << scientific << 0.0 << " "
                << setw(12) << setprecision(3) << fixed << q_factors[k];
        }
        rInfo(buf.str().c_str());
    }

 
        // Export eigensolutions to VTK file
         //
        if (params_.get<string>("vtk_file_name") != "")
         {
                pbe_start(134, "VTK export");

                rInfo("Export eigensolutions to VTK file...");

                postprocess::VtkExport exporter(femax_mesh_->get_nedelec_mesh(),
                                        lambda,
                                        Q);

                exporter.export_eigenfields(params_.get<string>("vtk_file_name"));
#if 0
        int nof_cells = 1000;
        exporter.export_eigenfields2("test_rectangular.vtk", nof_cells);
#endif

#if 0
        mesh::Vector3 a(1.0, 1.0, 0.0);
        mesh::Vector3 b(1.0, 1.0, 0.77);
        exporter.export_a_to_b("test_a_b.vtk", a, b, 100);
#endif
        pbe_stop(134);
    }

//       h5_write_eigenmodes(params_.get<string>("vtk_file_name"), get_eigenvectors(), get_eigenvalues());

         // log unused parameters
    ostringstream buf;
    buf << "Parameter list:" << endl;
    buf << "        [unused] flag is not shown properly for ML paramter lists :-(" << endl;
    params_.print(buf, 8);
    rDebug(buf.str().c_str());
}

void FemaxxDriver::run() 
{
    load_mesh();

  rInfo("starting calculation of eigenmodal solutions");
  calculate_eigenfields();

  rInfo("preparing storage of computed eigenmodal solution");
        h5_create_empty_file(params_.get<string>("eigensolutiondata_filename").c_str());
        
  rInfo("writing eigenmesh into HDF5 output file = %s", params_.get<string>("eigensolutiondata_filename").c_str());
        h5_write_eigenmesh(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_);
                 
  rInfo("writing mesh nodal coordinates into HDF5 output file");
        h5_write_eigenpoint(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_);
                 

  rInfo("writing eigenvalues into HDF5 file [master process only]");
  h5_write_eigenvalue(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_, femax_mesh_->get_nedelec_mesh(),*Q_,lambda_);
  rInfo("wrote eigenvalues into HDF5 file [master process only]");

  rInfo("writing eigenmodal fields sampled at centroids into HDF5 output file");
  h5_write_eigenfield(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_, femax_mesh_->get_nedelec_mesh(),*Q_,lambda_);


  rInfo("sampling computed modal solution on a cartesian grid");
  rInfo("this is contribution of processor %d", comm_.MyPID());

  double xmin=params_.get<double>("xaxis_cartesian_sampling_start_location");
  double xmax=params_.get<double>("xaxis_cartesian_sampling_stop_location");
  int nx=params_.get<int>("number_samples_x_axis");
  rInfo("sampling interval on x-axis = [%12.10e, %12.10e] & number of samples on x-axis = %d",xmin, xmax, nx);

  double ymin=params_.get<double>("yaxis_cartesian_sampling_start_location");
  double ymax=params_.get<double>("yaxis_cartesian_sampling_stop_location");
  int ny=params_.get<int>("number_samples_y_axis");
  rInfo("sampling interval on y-axis = [%12.10e, %12.10e] & number of samples on y-axis = %d",ymin, ymax, ny);

  double zmin=params_.get<double>("zaxis_cartesian_sampling_start_location");
  double zmax=params_.get<double>("zaxis_cartesian_sampling_stop_location");
  int nz=params_.get<int>("number_samples_z_axis");
  rInfo("sampling interval on z-axis = [%12.10e, %12.10e] & number of samples on z-axis = %d",zmin, zmax, nz);

  rInfo("total number of samples[nx*ny*nz] = %d",nx*ny*nz);

  h5_cartesian_sampling(params_.get<string>("eigensolutiondata_filename").c_str(), comm_,
                        tmesh_, femax_mesh_->get_nedelec_mesh(),
                        *Q_, lambda_,
                        xmin,  xmax,  nx,
                        ymin,  ymax,  ny,
                        zmin,  zmax,  nz);

}

const Epetra_MultiVector FemaxxDriver::get_eigenvectors() {
    return Epetra_MultiVector(View, *Q_, 0, nof_converged_);
}

Epetra_SerialDenseVector FemaxxDriver::get_eigenvalues() {
    return Epetra_SerialDenseVector(View, lambda_, nof_converged_);
}

---