examples/comet_stems.cpp

Go to the documentation of this file.

//
// C++ Implementation: comet_stems
//
// Description: Main program for eigenmode calculations for the COMET cyclotron.
//
// TODO: Voltage normalisation: 80kV
//
// Author: Roman Geus <roman.geus@psi.ch>, (C) 2004, 2005
//
// Copyright: See COPYING file that comes with this distribution
//
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <stdlib.h>

#include "rlog/rlog.h"
#include "rlog/rloglocation.h"
#include "rlog/Error.h"
#include "rlog/RLogChannel.h"
#include "rlog/StdioNode.h"
#include "rlog/RLogTime.h"
#include "myrlog.h"

#include <algorithm>
#include <iostream>
#include <fstream>
#include <cmath>

//#include "Epetra_ConfigDefs.h"
#include <ml_include.h>

#ifdef HAVE_MPI
#include <mpi.h>
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif

#include "optparse/TypedOptionParser.h"
#include "pbe.h"

#include <Trilinos_Util_CommandLineParser.h>
#include <Teuchos_ParameterList.hpp>
#include <ml_epetra_preconditioner.h>

#include "adaptsim.h"
#include "spline.h"
#include "romberg.h"
#include "CircleZ.h"
#include "EfieldIntegrand.h"
#include "LoggingFunction.h"
#include "clp_utils.h"
#include "vector.h"
#include "tetmesh.h"
#include "nedelecmesh.h"
#include "femaxmesh.h"
#include "vtkexport.h"
#include "femaxxdriver.h"

// Disable Romberg integration and use only adaptsim.
#define DISABLE_QROMB 1

using namespace rlog;
using namespace std;

void finalize() {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
}

static int parse_cmd_args(int argc, char *argv[], const Epetra_Comm& comm, Teuchos::ParameterList& params) {
    using optparse::STORE;
    using optparse::STORE_TRUE;
    using optparse::STORE_FALSE;
    using optparse::STRING;
    using optparse::BOOL;
    using optparse::INT;
    using optparse::DOUBLE;

    // Initialise parameter list with driver defaults.
    FemaxxDriver::set_defaults(params);
    
    // Add application-specific parameters
    params.set("stem_z_threshold", 0.25400);

    // Setup option parser using defaults from parameter list    
    string usage("comet_stems\n\nEigenmode solver application for the COMET cyclotron\n\nUsage: comet_stems [options]\n\noptions:");
    optparse::TypedOptionParser parser(usage);

    // Application-specific default values
    char* home_dir = getenv("HOME");        // Default mesh is expected to be in $HOME/meshes
    string mesh_path;
    if (home_dir != 0)
        mesh_path += string(home_dir) + "/meshes/";
    mesh_path += "P_1499_1641_1641_1612.ng";
    parser.add_option("-m", "--mesh", "mesh_file_name", "Mesh file in HDF5 format", STORE, STRING, mesh_path);
    parser.add_option("", "--sym-plane-config", "sym_plane_config", "Boundary conditions on symmetry planes", STORE, INT, "1");
    parser.add_option("-s", "--sigma", "sigma", "Shift for preconditioner of A - sigma*M", STORE, DOUBLE, "1.1");
    
    // Application-specific command line parameters
    parser.add_option("", "--stem_height_1", "stem_height_1", "Hight of stem 1 in meters", STORE, DOUBLE, "0.218854");
    parser.add_option("", "--stem_height_2", "stem_height_2", "Hight of stem 2 in meters", STORE, DOUBLE, "0.239586");
    parser.add_option("", "--stem_height_3", "stem_height_3", "Hight of stem 3 in meters", STORE, DOUBLE, "0.239586");
    parser.add_option("", "--stem_height_4", "stem_height_4", "Hight of stem 4 in meters", STORE, DOUBLE, "0.235352");
   
    // Options taken over from femaxx_driver
    parser.add_option("-h", "--help", "help", "Show this message", STORE_TRUE, BOOL, "0");
    //parser.add_option("-m", "--mesh", "mesh_file_name", "Mesh file in NETGEN format", STORE, STRING);
    //parser.add_option("-s", "--sigma", "sigma", "Shift for preconditioner of A - sigma*M", STORE, DOUBLE);
    parser.add_option(params, "-n", "--no-logfile", "disable_log", "Disable logfile", STORE_TRUE);
    parser.add_option(params, "-o", "--order", "element_order", "Finite element order", STORE, "1,2");
    parser.add_option(params, "-k", "--kmax", "kmax", "Number of eigensolutions to be computed", STORE);
    parser.add_option(params, "", "--tol", "eigtol", "Error tolerance for eigenpairs", STORE);
    parser.add_option(params, "", "--export-mtx", "export_mtx", "Export matrices in MatrixMarket format", STORE_TRUE);
    parser.add_option(params, "-v", "--vtk", "vtk_file_name", "Export eigenfields to VTK file", STORE);
    //parser.add_option(params, "", "--sym-plane-config", "sym_plane_config", "Boundary conditions on symmetry planes", STORE);
    parser.add_option(params, "-d", "--debug", "debug", "Enable debugging code", STORE_TRUE);
    parser.add_option(params, "", "--eigsolver", "eigsolver", "Eigensolver", STORE, "jdsym,lobpcg,knyazev");
    parser.add_option(params, "", "--aprec", "asigma_precon.type", "Preconditioner type for A - sigma*M", STORE, "neumann,ml,lu,2level,diag,if");
    parser.add_option(params, "", "--a11solver", "asigma_precon.solver11", "Solver for (1,1)-block of A - sigma*M", STORE, "lu,ml,if");
    parser.add_option(params, "", "--a22solver", "asigma_precon.solver22", "Solver for (2,2)-block of A - sigma*M", STORE, "diag,gmres_incomplete_block_jacobi,pcg_jacobi,blockdiag,none");
    parser.add_option(params, "", "--hsolver", "h_solver.type", "Solver type for H", STORE, "lu,pcg,prec");
    parser.add_option(params, "", "--hprec", "h_precon.type", "Preconditioner type for H", STORE, "ml,2level,diag,if");
    parser.add_option(params, "", "--h11solver", "h_precon.solver11", "Solver for (1,1)-block of H", STORE, "lu,ml,pcg_ml,if");
    parser.add_option(params, "", "--h22solver", "h_precon.solver22", "Solver for (2,2)-block of H", STORE, "diag,gmres_incomplete_block_jacobi,pcg_jacobi,blockdiag,none");
    parser.add_option(params, "", "--jtau", "jd_params.tau", "Target value", STORE);
    parser.add_option(params, "", "--jmin", "jd_params.jmin", "Min. size of search space", STORE);
    parser.add_option(params, "", "--jmax", "jd_params.jmax", "Max. size of search space", STORE);
    parser.add_option(params, "", "--jitmax", "jd_params.max_it", "Max. #JD iterations", STORE);
    parser.add_option(params, "", "--jlinitmax", "jd_params.max_inner_it", "Max. #of inner iterations per solve", STORE);
    parser.add_option(params, "", "--jtoldecay", "jd_params.tol_decay", "Factor by which inner error tol is decreasing", STORE);
    parser.add_option(params, "", "--eps_tr", "jd_params.eps_tr", "Threshold for shift in correction equation", STORE);
    try {
        parser.parse_args(argc, argv);
           
        // Test if all required options have been set
        if (parser.get_option("mesh_file_name") == "")
            throw optparse::OptionError("mesh", "Missing required option");
        if (parser.get_option("sigma") == "")
            throw optparse::OptionError("sigma", "Missing required option");
    
        if (parser.get_option("help") == "1") {
            if (comm.MyPID() == 0)
                parser.help(cerr);
            return 1;
        }
    }
    catch(optparse::OptionError e) {
        if (comm.MyPID() == 0) {
            cerr << e.what() << endl << endl;
            parser.help(cerr);
        }
        return 1;
    }
    
    if (parser.arguments.size() != 0) {
        if (comm.MyPID() == 0) {
            cerr << "Unexpected argument \"" << parser.arguments[0] << "\" found." << endl << endl;
            parser.help(cerr);
        }
        return 1;
    }
    
    // Update parameter list with values from parsed options.
    parser.set_parameter_list(params);
    
    return 0;
}

class CometGapFinder {
public:
    CometGapFinder(std::string mesh_type, const mesh::TetMesh& tmesh)
        : tmesh_(tmesh)
    {
        if (mesh_type == "P_1499_1641_1641_1612") {
            // Tabulated functions (used for spline interpolation
            double rad_entry_in[] = {0.10243483, 0.24759676, 0.37891511, 0.48161345, 0.58906747, 0.67218756, 0.74142964, 0.82499992};
            double phi_entry_in[] = {4.06171642, 4.37300720, 4.62760356, 4.77981775, 4.94923057, 5.10438120, 5.27123020, 5.32050234};
            double rad_entry_out[] = {0.08249764, 0.15229278, 0.39877397, 0.60499174, 0.81999979};
            double phi_entry_out[] = {4.59501326, 4.67390096, 4.91424909, 5.31862873, 5.59166061};
            double rad_exit_in[] = {0.13510732, 0.18577222, 0.24290023, 0.27314372, 0.33987641, 0.41934085, 0.46370478, 0.53079440, 0.61135473, 0.68823301, 0.73127536, 0.75081509, 0.76798017};
            double phi_exit_in[] = {3.74278291, 3.58438653, 3.62984776, 3.71049239, 3.76697089, 3.87803507, 3.95716874, 4.06445940, 4.22136555, 4.41000000, 4.52095706, 4.58974001, 4.66243071};
            double rad_exit_out[] = {0.06894548, 0.15229278, 0.27732126, 0.40163350, 0.60412626, 0.70856021, 0.81999968};
            double phi_exit_out[] = {2.96489808, 3.10310464, 3.22001175, 3.46130553, 3.77422794, 4.08845302, 4.24000150};

            // Copy arrays into vectors (insert rather than assign)
            std::copy(rad_entry_in, rad_entry_in + sizeof(rad_entry_in)/sizeof(rad_entry_in[0]),
                      std::inserter(rad_entry_in_, rad_entry_in_.begin()));
            std::copy(phi_entry_in, phi_entry_in + sizeof(phi_entry_in)/sizeof(phi_entry_in[0]),
                      std::inserter(phi_entry_in_, phi_entry_in_.begin()));
            std::copy(rad_entry_out, rad_entry_out + sizeof(rad_entry_out)/sizeof(rad_entry_out[0]),
                      insert_iterator< vector<double> > (rad_entry_out_, rad_entry_out_.begin()));
            std::copy(phi_entry_out, phi_entry_out + sizeof(phi_entry_out)/sizeof(phi_entry_out[0]),
                      insert_iterator< vector<double> > (phi_entry_out_, phi_entry_out_.begin()));
            std::copy(rad_exit_in, rad_exit_in + sizeof(rad_exit_in)/sizeof(rad_exit_in[0]),
                      insert_iterator< vector<double> > (rad_exit_in_, rad_exit_in_.begin()));
            std::copy(phi_exit_in, phi_exit_in + sizeof(phi_exit_in)/sizeof(phi_exit_in[0]),
                      insert_iterator< vector<double> > (phi_exit_in_, phi_exit_in_.begin()));
            std::copy(rad_exit_out, rad_exit_out + sizeof(rad_exit_out)/sizeof(rad_exit_out[0]),
                      insert_iterator< vector<double> > (rad_exit_out_, rad_exit_out_.begin()));
            std::copy(phi_exit_out, phi_exit_out + sizeof(phi_exit_out)/sizeof(phi_exit_out[0]),
                      insert_iterator< vector<double> > (phi_exit_out_, phi_exit_out_.begin()));

            // Spline preprocessing (compute coefficients)
            NR::spline(rad_entry_in_, phi_entry_in_, 1e40, 1e40, phi_entry_in2_);
            NR::spline(rad_entry_out_, phi_entry_out_, 1e40, 1e40, phi_entry_out2_);
            NR::spline(rad_exit_in_, phi_exit_in_, 1e40, 1e40, phi_exit_in2_);
            NR::spline(rad_exit_out_, phi_exit_out_, 1e40, 1e40, phi_exit_out2_);
        } else {
            // Mesh type not supported
            rExit("Mesh not supported by CometGapFinder.");
        }

        // Compute min and max radius
        min_radius_ = rad_entry_in_[0];
        min_radius_ = std::max<double>(min_radius_, rad_entry_out_[0]);
        min_radius_ = std::max<double>(min_radius_, rad_exit_in_[0]);
        min_radius_ = std::max<double>(min_radius_, rad_exit_out_[0]);

        max_radius_ = rad_entry_in_.end()[-1];
        max_radius_ = std::min<double>(max_radius_, rad_entry_out_.end()[-1]);
        max_radius_ = std::min<double>(max_radius_, rad_exit_in_.end()[-1]);
        max_radius_ = std::min<double>(max_radius_, rad_exit_out_.end()[-1]);
    }

    double get_min_radius() const {
        return min_radius_;
    }

    double get_max_radius() const {
        return max_radius_;
    }

    void get_gap(int gap_id, double radius, double& phi_entry, double& phi_exit) {
        const double pi = ::acos(-1.0);
        const double z_gap_find = 0.2;                 // Height at which the wing walls are detected
        const double phi_range = 0.6;                  // Gap width. Angle across which the voltage is computed
        double phase[4] = {0.0, 0.5*pi, pi, 1.5*pi};   // Phase of the four wings wrt to wing 0
        int wing_id = gap_id/2;
        
        assert(0 <= gap_id && gap_id < 8);
        assert(min_radius_ <= radius && radius <= max_radius_);

        // Find start values for bisection using spline functions
        double phi_entry_in, phi_entry_out, phi_exit_in, phi_exit_out;
        NR::splint(rad_entry_out_, phi_entry_out_, phi_entry_out2_, radius, phi_entry_out);
        NR::splint(rad_entry_in_, phi_entry_in_, phi_entry_in2_, radius, phi_entry_in);
        NR::splint(rad_exit_out_, phi_exit_out_, phi_exit_out2_, radius, phi_exit_out);
        NR::splint(rad_exit_in_, phi_exit_in_, phi_exit_in2_, radius, phi_exit_in);

        // Find mesh boundaries using bisection
        CircleZ circle_f(radius, z_gap_find);

        rDebug("radius=%lf, phi_entry_out=%lf [%s], phi_entry_in=%lf [%s], phi_exit_out=%lf [%s], phi_exit_in=%lf [%s]",
               radius,
               phi_entry_out, tmesh_.is_inside(circle_f(phi_entry_out)) ? "T" : "F",
               phi_entry_in, tmesh_.is_inside(circle_f(phi_entry_in)) ? "T" : "F",
               phi_exit_out, tmesh_.is_inside(circle_f(phi_exit_out)) ? "T" : "F",
               phi_exit_in, tmesh_.is_inside(circle_f(phi_exit_in)) ? "T" : "F");

        double phi_begin = tmesh_.find_boundary<CircleZ>(circle_f, phi_exit_in, phi_exit_out);
        double phi_end = tmesh_.find_boundary<CircleZ>(circle_f, phi_entry_in, phi_entry_out);
        rDebug("wing_id=%d, radius=%lf, phi_begin=%lf, phi_end=%lf", wing_id, radius, phi_begin, phi_end);

        phi_begin += phase[wing_id];
        phi_end += phase[wing_id];
        if (gap_id % 2 == 0) {
            // gap at beginning of the wing
            phi_entry = phi_begin - 0.5*phi_range;
            phi_exit = phi_begin + 0.5*phi_range;
        } else {
            // gap at end of wing
            phi_entry = phi_end - 0.5*phi_range;
            phi_exit = phi_end + 0.5*phi_range;
        }            
    }

private:
    const mesh::TetMesh& tmesh_;
    double min_radius_;
    double max_radius_;
    std::vector<double> rad_entry_in_, rad_entry_out_, phi_entry_in_, phi_entry_out_, phi_entry_in2_, phi_entry_out2_;
    std::vector<double> rad_exit_in_, rad_exit_out_, phi_exit_in_, phi_exit_out_, phi_exit_in2_, phi_exit_out2_;
};

void eval_tangetial_efield(FemaxxDriver& driver, const colarray::Vector<double>& q, string gnuplot_filename) {
    const double pi = ::acos(-1.0);
    const int nof_radii = 10;
    const int nof_angles = 1440;
    ofstream of(gnuplot_filename.c_str());
    mesh::TetMesh& tmesh(driver.get_tet_mesh());
    NedelecMesh& nedelec_mesh = driver.get_femax_mesh().get_nedelec_mesh();
    
    rAssert(tmesh.get_octree() != 0);

    CometGapFinder gap_finder("P_1499_1641_1641_1612", tmesh);
    const double radius_min = gap_finder.get_min_radius();
    const double radius_max = gap_finder.get_max_radius();

    for (int j = 0; j < nof_angles; ++ j) {
        double phi = 2.0*pi*j/nof_angles;
        of << phi;
        
        for (int i = 0; i <= nof_radii; ++ i) {
            double radius = radius_min + i * (radius_max - radius_min) / nof_radii;
            CircleZ circle_f(radius, 0.0);
            EfieldDs<CircleZ> integrand(circle_f, nedelec_mesh, q);
                    
            double e_tangential = integrand(phi);
            of << " " << e_tangential;
        }
        of << endl;
    }
    of.close();
}


void eval_gap_voltage(FemaxxDriver& driver, const colarray::Vector<double>& q, string gnuplot_filename) {
    const double tol = 1e-10;
    int info;
    const int nof_radii = 20;
    mesh::TetMesh& tmesh(driver.get_tet_mesh());
    NedelecMesh& nedelec_mesh = driver.get_femax_mesh().get_nedelec_mesh();
    ofstream of(gnuplot_filename.c_str());
    
    rAssert(tmesh.get_octree() != 0);

    CometGapFinder gap_finder("P_1499_1641_1641_1612", tmesh);
    const double radius_min = gap_finder.get_min_radius();
    const double radius_max = gap_finder.get_max_radius();

    rDebug("radius_min=%lf, radius_max=%lf, steps=%d", radius_min, radius_max, nof_radii);

    for (int i = 0; i <= nof_radii; ++ i) {
        double radius = radius_min + i * (radius_max - radius_min) / nof_radii;
        double phi_entry, phi_exit;
        CircleZ circle_f(radius, 0.0);
        EfieldDs<CircleZ> integrand(circle_f, nedelec_mesh, q);
        LoggingFunction<EfieldDs<CircleZ>, double> log_integrand(integrand);
        
        of << radius;

        for (int gap_id = 0; gap_id < 8; ++ gap_id) {
            // Get entry and exit angle for gap gap_id
            gap_finder.get_gap(gap_id, radius, phi_entry, phi_exit);
            
            // Compute voltage using romberg
#ifndef DISABLE_QROMB
            double gap_voltage_qromb = NR::qromb(integrand, phi_entry, phi_exit, tol, info);
            rInfo("gap %d: radius=%f, phi=%f->%f, voltage=%e", gap_id, radius, phi_entry, phi_exit, gap_voltage_qromb);
#endif
            // Open new log file (for use with gnuplot)    
            if (driver.get_comm().MyPID() == 0) {
                stringstream os;
                os << "comet_stems-" << i << "-" << gap_id << ".tab";
                log_integrand.enable_plot(os.str());
            }
            
            // Compute voltage using adaptsim
            double gap_voltage_adaptsim = NR::adaptsim(log_integrand, phi_entry, phi_exit, tol, info);
            rInfo("gap %d: radius=%f, phi=%f->%f, voltage=%e", gap_id, radius, phi_entry, phi_exit, gap_voltage_adaptsim);
            rInfo("num_eval=%d, y_min=%e, y_max=%e", log_integrand.get_num_eval(), log_integrand.get_y_min(), log_integrand.get_y_max());

#ifndef DISABLE_QROMB
            if (fabs(gap_voltage_adaptsim - gap_voltage_qromb) > 0.001 * fabs(gap_voltage_adaptsim)) {
                rWarning("Voltages computed by qromb and adaptsim differ by more than 1/1000.");
            }
#endif
            
            // Reset statistics and close logfile
            log_integrand.reset();
            
            // write voltage to output file
            of << " " << gap_voltage_adaptsim;
        } // for gap_id
        of << endl;
    } // for i
    
    of.close();
}

void export_mesh(string file, const mesh::TetMesh* tmesh) {
    int nofPoints = tmesh->get_nof_points();
    int nofTets = tmesh->get_nof_tets();

    mesh::Vector3 coord;
    mesh::Tet* curTet;
    std::ofstream of;

    of.open(file.c_str());
    assert(of.is_open());

    of.precision(5);
    of << "# vtk DataFile Version 2.0" << endl;                      // first line of a vtk file
    of << "generated using VtkExport::export_eigenfields" << endl;   // header
    of << "ASCII" << endl << endl;                                   // file format
    of << "DATASET UNSTRUCTURED_GRID" << endl;
    of << "POINTS " << nofPoints << " float" << endl;
    of << scientific;
    for (int i = 0; i < nofPoints  ; i++) {
        coord   = tmesh->get_point(i)->get_coord();
        of << coord.x << " " << coord.y << " " << coord.z << endl;
    }
    of << endl;
    of << "CELLS " <<  nofTets << " " << 5*nofTets << endl;

    // for each cell and write IDs of its points (preceded by number of points
    // in that cell)
    for (int i = 0; i < nofTets; i ++) {
        curTet = tmesh->get_tet(i);
        of << 4;    //assume that it is always a tetrahedron, i.e. with 4 corners
        for (int j = 0; j < 4; j ++) {
            of << " " << curTet->get_corner(j);
        }
        of << endl;
    }
    of << endl;

    // cell-type always = 10 as it is a tetrahedron
    of << "CELL_TYPES " << nofTets << endl;
    for (int i = 0; i < nofTets; i++) {
        of << 10 << endl;
    }
    of << endl;

    // write dummy point data (zero vectors)
    of << "POINT_DATA " << nofPoints << endl << endl;

    of << "VECTORS dummy" << " float" << endl;
    for (int j =0; j < nofPoints; j++) {
        of << 0.0 << " " << 0.0 << " " << 0.0 << endl;
    }
    of << endl;
}

int main(int argc, char *argv[]) {

#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
    atexit(finalize);
    Epetra_MpiComm comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm comm;
#endif

#ifdef __GNUC__
    std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
#endif

    Teuchos::ParameterList params;
    // Parse command line using optparse and initialise parameter list
    if (parse_cmd_args(argc, argv, comm, params) != 0)
        // optparse discovered errors: exit program
        return 1;
    
    RLogInit(argc, argv);

    MyStdioNode stdLog;     // Logs to stderr
    StdioNode* dbgLog(0);   // Logs to file

    DEF_CHANNEL("warning/all", Log_Warning);
    DEF_CHANNEL("error/all", Log_Error);
    DEF_CHANNEL("info/all", Log_Info);
    DEF_CHANNEL("debug/all", Log_Debug);

    if (comm.MyPID() == 0) {
        stdLog.subscribeTo(RLOG_CHANNEL("warning"));
        stdLog.subscribeTo(RLOG_CHANNEL("error"));
        stdLog.subscribeTo(RLOG_CHANNEL("info"));
    } else {
        stdLog.subscribeTo(RLOG_CHANNEL("error"));
        stdLog.subscribeTo(RLOG_CHANNEL("warning/all"));
        stdLog.subscribeTo(RLOG_CHANNEL("info/all"));
    }
    string logname(get_log_file_name("comet_stems_"));
    if (!params.get<bool>("disable_log") && comm.MyPID() == 0) {
        int fid = ::open(logname.c_str(), O_WRONLY|O_CREAT|O_TRUNC, 0666);
        rAssert(fid != -1);
        dbgLog = new StdioNode(fid);
        dbgLog->subscribeTo(RLOG_CHANNEL("debug"));
        dbgLog->subscribeTo(RLOG_CHANNEL("info"));
        dbgLog->subscribeTo(RLOG_CHANNEL("warning"));
        dbgLog->subscribeTo(RLOG_CHANNEL("error"));
    }

    rInfo("Comet_stems running on %d processor(s)...", comm.NumProc());
    rLog(RLOG_CHANNEL("info/all"), "Hi from processor %d.", comm.MyPID());
    if (!params.get<bool>("disable_log"))
        rInfo("Logging debug output to %s", logname.c_str());

    // Log command-line options and parameter list
    log_command_line_args(argc, argv);
    if (comm.MyPID() == 0) {
        rDebug("Input parameters:");
        ostringstream buf;
        buf << endl;
        params.print(buf, 8);
        rDebug(buf.str().c_str());
    }

    FemaxxDriver driver(comm, params);
    driver.load_mesh();

    if (params.get<string>("vtk_file_name") != "") {
        string file_name = params.get<string>("vtk_file_name") + "_mesh_orig.vtk";
        if (comm.MyPID() == 0) {
            export_mesh(file_name, &driver.get_tet_mesh());
        }
    }

    rInfo("Rescale COMET stems...");
    // Extract mesh points belonging to stems
    const double z_threshold = params.get<double>("stem_z_threshold");

    int stem_count[4] = {0, 0, 0, 0};
    double stem_zmax[4] = {0.0, 0.0, 0.0, 0.0};
    double stem_factor[4] = {0.0, 0.0, 0.0, 0.0};
    double stem_newheight[4] = { params.get<double>("stem_height_1"),
                                 params.get<double>("stem_height_2"),
                                 params.get<double>("stem_height_3"),
                                 params.get<double>("stem_height_4") };

    pbe_start(10, "Stem node adjustment")
    mesh::TetMesh& tmesh(driver.get_tet_mesh());
    for (mesh::TetMesh::point_iterator it(tmesh.point_begin()); it < tmesh.point_end(); ++ it) {
        mesh::Vector3 coord((*it).get_coord());
        if (coord.z > z_threshold) {
            int stem_id = 0;
            if (coord.x > 0.0)
                stem_id += 1;
            if (coord.y > 0.0)
                stem_id += 2;
            stem_count[stem_id] ++;
            if (coord.z > stem_zmax[stem_id])
                stem_zmax[stem_id] = coord.z;
        }
    }

    if (comm.MyPID() == 0) {
        for (int i = 0; i < 4; i ++)
            rDebug("Stem %d: nof_points=%d, zmax=%lf", i, stem_count[i], stem_zmax[i]);
    }

    // reposition stem mesh points according to input parameters
    for (int i = 0; i < 4; i ++)
        stem_factor[i] = stem_newheight[i] / (stem_zmax[i] - z_threshold);

    for (mesh::TetMesh::point_iterator it(tmesh.point_begin()); it < tmesh.point_end(); ++ it) {
        mesh::Vector3& coord((*it).get_coord());
        if (coord.z > z_threshold) {
            int stem_id = 0;
            if (coord.x > 0.0)
                stem_id += 1;
            if (coord.y > 0.0)
                stem_id += 2;
            stem_count[stem_id] ++;
            if (coord.z > stem_zmax[stem_id])
                stem_zmax[stem_id] = coord.z;
            coord.z = (coord.z - z_threshold) * stem_factor[stem_id] + z_threshold;
        }
    }

    if (comm.MyPID() == 0) {
        for (int i = 0; i < 4; i ++)
            rInfo("Stem %d: old height=%lf, new height=%lf", i, stem_zmax[i] - z_threshold, stem_newheight[i]);
    }
    pbe_stop(10);

    if (params.get<string>("vtk_file_name") != "") {
        string file_name = params.get<string>("vtk_file_name") + "_mesh_modded.vtk";
        if (comm.MyPID() == 0) {
            export_mesh(file_name, &driver.get_tet_mesh());
        }
    }

    // calculate eigenfields
    driver.calculate_eigenfields();

    // reload mesh data structures
    driver.load_mesh();

    // print eigenfield info and export eigenfields
    driver.postprocess();

    // evaluate gap voltages
    rInfo("Evaluate gap voltages...");
    postprocess::VtkExport post(driver.get_femax_mesh().get_nedelec_mesh(),
                                driver.get_eigenvalues(),
                                driver.get_eigenvectors());
    colarray::Vector<double>* q(0);
    double lambda;
    if (comm.MyPID() == 0)
        q = new colarray::Vector<double>(post.get_nof_dofs());
    // get 1st eigenpair
    post.get_eigenpair(0, &lambda, q);
    if (comm.MyPID() == 0) {
        eval_tangetial_efield(driver, *q, "tangential.tab");
        eval_gap_voltage(driver, *q, "gap_voltages.tab");
        delete q;
    }

    return EXIT_SUCCESS;
}

---