examples/postprocess_serial.cpp

Go to the documentation of this file.

#include "unistd.h"
#include <cmath>
#include <iostream> 
#include "Teuchos_ParameterList.hpp"
#include "rlog/rlog.h"
#include "rlog/StdioNode.h"
#include "rlog/RLogChannel.h"
#include "myrlog.h"
#include "pbe.h"
#include "h5_tools.h"
#include "tetmesh/tetmesh.h"
#include "tetmesh/vector3.h"
#include "tetmeshbuilder.h"
#include "netgenreader.h"
#include "femaxmesh.h"
#include "nedelecmesh.h"
#include "lagrangemesh.h"

using namespace rlog;
using namespace std;

//
// Global variables
//
mesh::TetMesh* tmesh(0);
FemaxMesh* femax_mesh(0);
int num_eigenpairs(0);
colarray::ColumnVector<double> lambda;
colarray::Matrix<double> Q;

void log_mem_footprint(string msg) {
    size_t local_size = get_mem_footprint();
    if (msg != "")
        msg = " (" + msg + ")";
    rDebug("Memory footprint%s: %s",
           msg.c_str(),
           rlog::bytes2str(local_size).c_str());
}

void load_mesh(const Teuchos::ParameterList& params_) {
    // 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());

    log_mem_footprint("before loading mesh");
    TetMeshBuilder* builder = new TetMeshBuilder();
    NetgenReader* reader = new NetgenReader(&istr, builder);
    reader->read();
    delete reader;
    tmesh = builder->get_mesh();
    delete builder;
    tmesh->log_mesh_info();
    log_mem_footprint("after loading tetmesh");

    // 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);
    femax_mesh = new FemaxMesh(*tmesh, params_.get<int>("element_order"), sym_plane_config);
    log_mem_footprint("after calc dof mappings");
}

double q_factor(int k) {
    colarray::ColumnVector<double> q(Q, k);
    return femax_mesh->get_nedelec_mesh().q_factor(lambda[k], q._v);
}

void eigenmode_summary() {
    // 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 < num_eigenpairs; ++ 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(12) << setprecision(3) << fixed << q_factor(k);
    }
    rInfo(buf.str().c_str());
}

void vtk_export_eigenfields(string file) {
    const int precision = 4;
    NedelecMesh* nedelec_mesh(&femax_mesh->get_nedelec_mesh());
    int nofPoints = tmesh->get_nof_points();
    int nofTets = tmesh->get_nof_tets();

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

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

    of.precision(precision);
    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_id(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;

    of << "POINT_DATA " << nofPoints << endl << endl;

    // E-field data is written into the vtk-file
    for (int k = 0; k < num_eigenpairs; k ++) {
        colarray::ColumnVector<double> q(Q, k);

        of << "VECTORS el_field_mode_" << k << " float" << endl;
        //get the value for each point
        for (int j =0; j < nofPoints; j++) {
            value = nedelec_mesh->eval(tmesh->get_point(j)->get_coord(), q);
            of << value.x << " " << value.y << " " << value.z << endl;
        }
        of << endl;
    } // for k

    // B-field data is written into the vtk-file
    for (int k = 0; k < num_eigenpairs; k ++) {
        colarray::ColumnVector<double> q(Q, k);

        of << "VECTORS mag_field_mode_" << k << " float" << endl;
        //get the value for each point
        for (int j = 0; j < nofPoints; j++) {
            value = nedelec_mesh->eval_curl(tmesh->get_point(j)->get_coord(), q);
            of << value.x << " " << value.y << " " << value.z << endl;
        }
        of << endl;
    } // for k

}

int main(int argc, char *argv[])
{
    Teuchos::ParameterList params;
    RLogInit(argc, argv);
    MyStdioNode stdLog;
    stdLog.subscribeTo(RLOG_CHANNEL("error"));
    stdLog.subscribeTo(RLOG_CHANNEL("warning"));
    stdLog.subscribeTo(RLOG_CHANNEL("info"));
    stdLog.subscribeTo(RLOG_CHANNEL("debug"));

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

    // Read parameter list
    h5_read_param_list("eigenmodes.h5", params);
    params.print(cout);
    
    // Read tetmesh
    pbe_start(52, "Mesh generation");
    load_mesh(params);
    pbe_stop(52);
    
    pbe_start(56, "Read eigenmodes");
    h5_read_eigenmodes("eigenmodes.h5", Q, lambda);
    pbe_stop(56);
    num_eigenpairs = lambda._n;
    eigenmode_summary();

    rInfo("Export modes mesh and eigenmodes to VTK file...");
    pbe_start(131, "Construction of octree");
    tmesh->construct_octree();
    pbe_stop(131);
    pbe_start(134, "VTK export");
    vtk_export_eigenfields("out.vtk");
    pbe_stop(134);
    
    delete femax_mesh;
    delete tmesh;

    pbe_stop(50);
    pbe_dump();
    pbe_finalize(0);

    return 0;
}

---