src/vtkexport.cpp

Go to the documentation of this file.

//
// C++ Implementation: vtkexport
//
// Description:
//
//
// Author: Roman Geus <roman.geus@psi.ch>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

// To get round from math.h
#define _ISOC99_SOURCE
       
#include <iostream>
#include <fstream>
#include <string>
#include <iomanip>
#include <vector>
#include <math.h>
#include "Epetra_Comm.h"
#include "Epetra_Export.h"
#include "Epetra_Map.h"
#include "Epetra_MultiVector.h"
#include "vector.h"
#include "romberg.h"
#include "EfieldIntegrand.h"
#include "vtkexport.h"

namespace postprocess {

VtkExport::VtkExport(const NedelecMesh::NedelecMesh& nedelec_mesh,
                     const Epetra_SerialDenseVector& lambda,
                     const Epetra_MultiVector& Q)
    : nedelec_mesh_(nedelec_mesh),
      lambda_(lambda),
      Q_(Q),
      comm_(Q_.Comm()),
      precision_(6)
{
}

void VtkExport::export_mesh(string file) {
    mesh::TetMesh* tmesh(nedelec_mesh_.getTetMesh());
    int nofPoints = tmesh->get_nof_points();
    int nofTets = tmesh->get_nof_tets();

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

    if (comm_.MyPID() == 0) {
        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;

        // 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;

    } // MyPID == 0
}

void VtkExport::export_eigenfields(string file) {
    double lambda;
    colarray::Vector<double>* q(0);

    mesh::TetMesh* tmesh(nedelec_mesh_.getTetMesh());
    int nofPoints = tmesh->get_nof_points();
    int nofTets = tmesh->get_nof_tets();

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

    if (comm_.MyPID() == 0) {
        q = new colarray::Vector<double>(get_nof_dofs());

        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 < get_nof_vectors(); k ++) {
        get_eigenpair(k, &lambda, q);

        if (comm_.MyPID() == 0) {
            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 < get_nof_vectors(); k ++) {
        get_eigenpair(k, &lambda, q);

        if (comm_.MyPID() == 0) {
            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

    if (comm_.MyPID() == 0)
        delete q;
}

void VtkExport:: export_eigenfields2(std::string file, int nof_cells) {
    double lambda;
    colarray::Vector<double>* q(0);

    int nofPoints;
    int dim[3];
    mesh::TetMesh* tmesh;
    mesh::Vector3 center;
    mesh::Vector3 box_len;
    mesh::Vector3 coord;
    mesh::Vector3 value;
    tmesh = nedelec_mesh_.getTetMesh();

    std::vector<double> xcoords, ycoords, zcoords;
    double curx, cury, curz;
    mesh::Vector3 point;
    mesh::Vector3 ref_point;
    std::vector<mesh::Tet *> tets;
    std::ofstream of;

    if (comm_.MyPID() == 0) {
        q = new colarray::Vector<double>(get_nof_dofs());

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

        of.precision(precision_);
        of << "# vtk DataFile Version 2.0" << endl;
        of << "generated using FemaXX" << endl;
        of << "ASCII" << endl << endl;
        of << "DATASET RECTILINEAR_GRID" << endl;

        // get the bounding box from tmesh and calculate the dimensions of the new grid
        tmesh->get_bounding_box(center, box_len);
        ref_point = center - box_len;
        box_len *= 2.0;

        getdimensions(nof_cells, box_len, dim);
        of << "DIMENSIONS " << dim[0] << " " << dim[1] << " " << dim[2] << endl;

        double xstep = box_len.x / dim[0];
        double ystep = box_len.y / dim[1];
        double zstep = box_len.z / dim[2];

        // create three vectors containing the x, y and z coordinates of the points
        // and do the insertion of these into the vtk file
        of << "X_COORDINATES " << dim[0] << " float" << endl;
        for (int i = 0; i <  dim[0]; i++) {
            curx = ref_point.x + (i * xstep);
            xcoords.push_back(curx);
            of << curx << " ";
        }
        of << endl;

        of << "Y_COORDINATES " << dim[1] << " float" << endl;
        for (int i = 0; i <  dim[1]; i++) {
            cury = ref_point.y + (i * ystep);
            ycoords.push_back(cury);
            of << cury << " ";
        }
        of << endl;

        of << "Z_COORDINATES " << dim[2] << " float" << endl;
        for (int i = 0; i <  dim[2]; i++) {
            curz = ref_point.z + (i * zstep);
            zcoords.push_back(curz);
            of << curz << " ";
        }
        of << endl;

        nofPoints = xcoords.size() *  ycoords.size() * zcoords.size() ;
        of << endl << "POINT_DATA " <<  nofPoints <<  endl << endl;
    }

    // writing the electric field values into the vtk file
    for (int k = 0; k < get_nof_vectors(); k++) {
        get_eigenpair(k, &lambda, q);

        if (comm_.MyPID() == 0) {
            of << "VECTORS el_field_mode_" << k << " float" << endl;
            //get the electric field  value for each point
            for (int i = 0; i < dim[2]; i ++) {
                point.z = zcoords.at(i);
                for (int j = 0; j < dim[1]; j ++) {
                    point.y = ycoords.at(j);
                    for (int h = 0; h < dim[0]; h ++) {
                        point.x = xcoords.at(h);
                        value = nedelec_mesh_.eval(point, *q);
                        of << value.x << " " << value.y << " " << value.z << endl;
                    }
                }
            }
            of << endl;
        }
    } // for k

    // writing the magnetic field values into the vtk file
    for (int k = 0; k < Q_.NumVectors(); k++) {
        get_eigenpair(k, &lambda, q);

        if (comm_.MyPID() == 0) {
            of << "VECTORS mag_field_mode_" << k << " float" << endl;
            //get the magnetic field  value for each point
            for (int i = 0; i < dim[2]; i ++) {
                point.z = zcoords.at(i);
                for (int j = 0; j < dim[1]; j ++) {
                    point.y = ycoords.at(j);
                    for (int h = 0; h < dim[0]; h ++) {
                        point.x = xcoords.at(h);
                        value = nedelec_mesh_.eval_curl(point, *q);
                        of << value.x << " " << value.y << " " << value.z << endl;
                    }
                }
            }
            of << endl;
        }
    } // for k

    if (comm_.MyPID() == 0)
        delete q;
}

void VtkExport::export_a_to_b(std::string file,
                              mesh::Vector3 a,
                              mesh::Vector3 b,
                              int nof_intervals)
{
    double lambda;
    colarray::Vector<double>* q(0);
    mesh::Vector3 d = b - a;
    mesh::Vector3 p;
    mesh::Vector3 value;
    double spacing = 1.0 / static_cast<double>(nof_intervals);
    ofstream of;

    if (comm_.MyPID() == 0) {
        q = new colarray::Vector<double>(get_nof_dofs());

        of.open(file.c_str());
        assert(of.is_open());
        of.precision(precision_);
        of << "# vtk DataFile Version 2.0" << endl;
        of << "Generated using FemaXX" << endl;
        of << "ASCII" << endl << endl;
        of << "DATASET STRUCTURED_POINTS" << endl;
        of << "DIMENSIONS " << nof_intervals + 1 << " " << 1 << " " << 1 << endl;
        of << "ORIGIN 0 0 0" << endl;
        of << "SPACING " << spacing << " 1 1" << endl;
        of << endl << "POINT_DATA " << nof_intervals + 1 <<  endl << endl;
    }

    for (int k = 0; k < get_nof_vectors(); k ++) {
        get_eigenpair(k, &lambda, q);
        if (comm_.MyPID() == 0) {
            of << "VECTORS el_field_mode_" << k << " float" << endl;
            for (int i = 0; i <= nof_intervals; i ++) {
                p = a + d * (static_cast<double>(i) * spacing);
                value = nedelec_mesh_.eval(p, *q);
                of << value.x << " " << value.y << " " << value.z << endl;
            } // for i
            of << endl;
        }
    } // for k

    for (int k = 0; k < get_nof_vectors(); k ++) {
        get_eigenpair(k, &lambda, q);
        if (comm_.MyPID() == 0) {
            of << "VECTORS mag_field_mode_" << k << " float" << endl;
            for (int i = 0; i <= nof_intervals; i ++) {
                p = a + d * (static_cast<double>(i) * spacing);
                value = nedelec_mesh_.eval_curl(p, *q);
                of << value.x << " " << value.y << " " << value.z << endl;
            } // for i
            of << endl;
        }
    } // for k

    if (comm_.MyPID() == 0)
        delete q;
}

void VtkExport::get_eigenpair(int k, double* lambda, colarray::Vector<double>* q) {
    assert(0 <= k && k < get_nof_vectors());
    assert(comm_.MyPID() > 0 || q->_n == get_nof_dofs());

    //
    // Bring k-th eigenvector to processor 0
    //

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

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

    // store eigenpair on processor 0
    if (comm_.MyPID() == 0) {
        *lambda = lambda_[k];
        for (int i = 0; i < get_nof_dofs(); i ++) {
            q->operator()(i) =  q0[i];
        }
    }
}

void VtkExport::getdimensions(int nof_cells, mesh::Vector3& box_len, int dim[3]) {

    cout << endl << "   Vtk Rectilinear grid: bounding box length l, gridsize, grid dimensions nx ny nz   "
         << endl
         <<         "-------------------------------------------------------------------------------------"
         << endl;
    cout << endl <<  "l.x = " << box_len.x << ", l.y = " << box_len.y << ", l.z = " << box_len.z;

    mesh::Vector3 d;
    d.x = ::pow(static_cast<double>(nof_cells) * box_len.x * box_len.x / box_len.y / box_len.z, 1.0/3.0);
    d.y = d.x * box_len.y / box_len.x;
    d.z = d.x * box_len.z / box_len.x ;

    // make sure that all components are at least 0.5
    for (int i = 0; i < 3; i ++)
        if (d[i] < 0.5)
            d *= 0.5/d[i];

    dim[0] = static_cast<int>(round(d.x));
    dim[1] = static_cast<int>(round(d.y));
    dim[2] = static_cast<int>(round(d.z));

    cout << endl << "vtk_gridsize = " << nof_cells << endl;
    cout << "nx = " << dim[0] << ", " << "ny = " << dim[1]
         << ", " << " nz = " << dim[2] << endl;
    cout << "nx * ny * nz = " << dim[0]*dim[1]*dim[2] << endl;
    cout << "lx/nx = " << box_len.x / dim[0] << ", ly/ny =" << box_len.y/ dim[1] <<
        ", lz/nz = " << box_len.z / dim[2] << endl;
}

class UnitLine {
public:
    mesh::Vector3 operator()(double x) const {
        return mesh::Vector3(x, 0, 0);
    }
};

void VtkExport::integrate_a_to_b(int mode_id,
                                 mesh::Vector3 a,
                                 mesh::Vector3 b,
                                 double& int_e)
{
    double lambda;
    colarray::Vector<double>* q(0);

    if (comm_.MyPID() == 0) {
        q = new colarray::Vector<double>(get_nof_dofs());
    }
    // bring eigenvector to proc 0
    get_eigenpair(mode_id, &lambda, q);

    if (comm_.MyPID() == 0) {
        mesh::Vector3 d = b - a;
        double len = d.length();
        d /= len;
        EfieldIntegrand romberg_integrand(a, d, nedelec_mesh_, *q);
        int info;
        int_e = NR::qromb<EfieldIntegrand>(romberg_integrand, 0, len, 1e-10, info);
    }

    if (comm_.MyPID() == 0)
        delete q;

    // debug code following...
    class UnitLine line;
    mesh::TetMesh* tmesh = nedelec_mesh_.getTetMesh();
    tmesh->find_boundary<UnitLine>(line, 0.0, 1.0);
}

}; // namespace

---