src/vtkexport.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: vtkexport
00003 //
00004 // Description:
00005 //
00006 //
00007 // Author: Roman Geus <roman.geus@psi.ch>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
00010 //
00011 //
00012 
00013 // To get round from math.h
00014 #define _ISOC99_SOURCE
00015        
00016 #include <iostream>
00017 #include <fstream>
00018 #include <string>
00019 #include <iomanip>
00020 #include <vector>
00021 #include <math.h>
00022 #include "Epetra_Comm.h"
00023 #include "Epetra_Export.h"
00024 #include "Epetra_Map.h"
00025 #include "Epetra_MultiVector.h"
00026 #include "vector.h"
00027 #include "romberg.h"
00028 #include "EfieldIntegrand.h"
00029 #include "vtkexport.h"
00030 
00031 namespace postprocess {
00032 
00033 VtkExport::VtkExport(const NedelecMesh::NedelecMesh& nedelec_mesh,
00034                      const Epetra_SerialDenseVector& lambda,
00035                      const Epetra_MultiVector& Q)
00036     : nedelec_mesh_(nedelec_mesh),
00037       lambda_(lambda),
00038       Q_(Q),
00039       comm_(Q_.Comm()),
00040       precision_(6)
00041 {
00042 }
00043 
00044 void VtkExport::export_mesh(string file) {
00045     mesh::TetMesh* tmesh(nedelec_mesh_.getTetMesh());
00046     int nofPoints = tmesh->get_nof_points();
00047     int nofTets = tmesh->get_nof_tets();
00048 
00049     mesh::Vector3 coord;
00050     mesh::Tet* curTet;
00051     std::ofstream of;
00052 
00053     if (comm_.MyPID() == 0) {
00054         of.open(file.c_str());
00055         assert(of.is_open());
00056 
00057         of.precision(precision_);
00058         of << "# vtk DataFile Version 2.0" << endl;                      // first line of a vtk file
00059         of << "generated using VtkExport::export_eigenfields" << endl;   // header
00060         of << "ASCII" << endl << endl;                                   // file format
00061         of << "DATASET UNSTRUCTURED_GRID" << endl;
00062         of << "POINTS " << nofPoints << " float" << endl;
00063         of << scientific;
00064         for (int i = 0; i < nofPoints  ; i++) {
00065             coord   = tmesh->get_point(i)->get_coord();
00066             of << coord.x << " " << coord.y << " " << coord.z << endl;
00067         }
00068         of << endl;
00069         of << "CELLS " <<  nofTets << " " << 5*nofTets << endl;
00070 
00071         // for each cell and write IDs of its points (preceded by number of points
00072         // in that cell)
00073         for (int i = 0; i < nofTets; i ++) {
00074             curTet = tmesh->get_tet(i);
00075             of << 4;    //assume that it is always a tetrahedron, i.e. with 4 corners
00076             for (int j = 0; j < 4; j ++) {
00077                 of << " " << curTet->get_corner_id(j);
00078             }
00079             of << endl;
00080         }
00081         of << endl;
00082 
00083         // cell-type always = 10 as it is a tetrahedron
00084         of << "CELL_TYPES " << nofTets << endl;
00085         for (int i = 0; i < nofTets; i++) {
00086             of << 10 << endl;
00087         }
00088         of << endl;
00089 
00090         // write dummy point data (zero vectors)
00091         of << "POINT_DATA " << nofPoints << endl << endl;
00092 
00093         of << "VECTORS dummy" << " float" << endl;
00094         for (int j =0; j < nofPoints; j++) {
00095             of << 0.0 << " " << 0.0 << " " << 0.0 << endl;
00096         }
00097         of << endl;
00098 
00099     } // MyPID == 0
00100 }
00101 
00102 void VtkExport::export_eigenfields(string file) {
00103     double lambda;
00104     colarray::Vector<double>* q(0);
00105 
00106     mesh::TetMesh* tmesh(nedelec_mesh_.getTetMesh());
00107     int nofPoints = tmesh->get_nof_points();
00108     int nofTets = tmesh->get_nof_tets();
00109 
00110     mesh::Vector3 coord;
00111     mesh::Vector3 value;
00112     mesh::Tet* curTet;
00113     std::ofstream of;
00114 
00115     if (comm_.MyPID() == 0) {
00116         q = new colarray::Vector<double>(get_nof_dofs());
00117 
00118         of.open(file.c_str());
00119         assert(of.is_open());
00120 
00121         of.precision(precision_);
00122         of << "# vtk DataFile Version 2.0" << endl;                      // first line of a vtk file
00123         of << "generated using VtkExport::export_eigenfields" << endl;   // header
00124         of << "ASCII" << endl << endl;                                   // file format
00125         of << "DATASET UNSTRUCTURED_GRID" << endl;
00126         of << "POINTS " << nofPoints << " float" << endl;
00127         of << scientific;
00128         for (int i = 0; i < nofPoints; ++ i) {
00129             coord   = tmesh->get_point(i)->get_coord();
00130             of << coord.x << " " << coord.y << " " << coord.z << endl;
00131         }
00132         of << endl;
00133         of << "CELLS " <<  nofTets << " " << 5*nofTets << endl;
00134 
00135         // for each cell and write IDs of its points (preceded by number of points
00136         // in that cell)
00137         for (int i = 0; i < nofTets; i ++) {
00138             curTet = tmesh->get_tet(i);
00139             of << 4;    //assume that it is always a tetrahedron, i.e. with 4 corners
00140             for (int j = 0; j < 4; j ++) {
00141                 of << " " << curTet->get_corner_id(j);
00142             }
00143             of << endl;
00144         }
00145         of << endl;
00146 
00147         // cell-type always = 10 as it is a tetrahedron
00148         of << "CELL_TYPES " << nofTets << endl;
00149         for (int i = 0; i < nofTets; i++) {
00150             of << 10 << endl;
00151         }
00152         of << endl;
00153 
00154         of << "POINT_DATA " << nofPoints << endl << endl;
00155     }
00156 
00157     // E-field data is written into the vtk-file
00158     for (int k = 0; k < get_nof_vectors(); k ++) {
00159         get_eigenpair(k, &lambda, q);
00160 
00161         if (comm_.MyPID() == 0) {
00162             of << "VECTORS el_field_mode_" << k << " float" << endl;
00163             //get the value for each point
00164             for (int j =0; j < nofPoints; j++) {
00165                 value = nedelec_mesh_.eval(tmesh->get_point(j)->get_coord(), *q);
00166                 of << value.x << " " << value.y << " " << value.z << endl;
00167             }
00168             of << endl;
00169         }
00170     } // for k
00171 
00172     // B-field data is written into the vtk-file
00173     for (int k = 0; k < get_nof_vectors(); k ++) {
00174         get_eigenpair(k, &lambda, q);
00175 
00176         if (comm_.MyPID() == 0) {
00177             of << "VECTORS mag_field_mode_" << k << " float" << endl;
00178             //get the value for each point
00179             for (int j = 0; j < nofPoints; j++) {
00180                 value = nedelec_mesh_.eval_curl(tmesh->get_point(j)->get_coord(), *q);
00181                 of << value.x << " " << value.y << " " << value.z << endl;
00182             }
00183             of << endl;
00184         }
00185     } // for k
00186 
00187     if (comm_.MyPID() == 0)
00188         delete q;
00189 }
00190 
00191 void VtkExport:: export_eigenfields2(std::string file, int nof_cells) {
00192     double lambda;
00193     colarray::Vector<double>* q(0);
00194 
00195     int nofPoints;
00196     int dim[3];
00197     mesh::TetMesh* tmesh;
00198     mesh::Vector3 center;
00199     mesh::Vector3 box_len;
00200     mesh::Vector3 coord;
00201     mesh::Vector3 value;
00202     tmesh = nedelec_mesh_.getTetMesh();
00203 
00204     std::vector<double> xcoords, ycoords, zcoords;
00205     double curx, cury, curz;
00206     mesh::Vector3 point;
00207     mesh::Vector3 ref_point;
00208     std::vector<mesh::Tet *> tets;
00209     std::ofstream of;
00210 
00211     if (comm_.MyPID() == 0) {
00212         q = new colarray::Vector<double>(get_nof_dofs());
00213 
00214         of.open(file.c_str());
00215         assert(of.is_open());
00216 
00217         of.precision(precision_);
00218         of << "# vtk DataFile Version 2.0" << endl;
00219         of << "generated using FemaXX" << endl;
00220         of << "ASCII" << endl << endl;
00221         of << "DATASET RECTILINEAR_GRID" << endl;
00222 
00223         // get the bounding box from tmesh and calculate the dimensions of the new grid
00224         tmesh->get_bounding_box(center, box_len);
00225         ref_point = center - box_len;
00226         box_len *= 2.0;
00227 
00228         getdimensions(nof_cells, box_len, dim);
00229         of << "DIMENSIONS " << dim[0] << " " << dim[1] << " " << dim[2] << endl;
00230 
00231         double xstep = box_len.x / dim[0];
00232         double ystep = box_len.y / dim[1];
00233         double zstep = box_len.z / dim[2];
00234 
00235         // create three vectors containing the x, y and z coordinates of the points
00236         // and do the insertion of these into the vtk file
00237         of << "X_COORDINATES " << dim[0] << " float" << endl;
00238         for (int i = 0; i <  dim[0]; i++) {
00239             curx = ref_point.x + (i * xstep);
00240             xcoords.push_back(curx);
00241             of << curx << " ";
00242         }
00243         of << endl;
00244 
00245         of << "Y_COORDINATES " << dim[1] << " float" << endl;
00246         for (int i = 0; i <  dim[1]; i++) {
00247             cury = ref_point.y + (i * ystep);
00248             ycoords.push_back(cury);
00249             of << cury << " ";
00250         }
00251         of << endl;
00252 
00253         of << "Z_COORDINATES " << dim[2] << " float" << endl;
00254         for (int i = 0; i <  dim[2]; i++) {
00255             curz = ref_point.z + (i * zstep);
00256             zcoords.push_back(curz);
00257             of << curz << " ";
00258         }
00259         of << endl;
00260 
00261         nofPoints = xcoords.size() *  ycoords.size() * zcoords.size() ;
00262         of << endl << "POINT_DATA " <<  nofPoints <<  endl << endl;
00263     }
00264 
00265     // writing the electric field values into the vtk file
00266     for (int k = 0; k < get_nof_vectors(); k++) {
00267         get_eigenpair(k, &lambda, q);
00268 
00269         if (comm_.MyPID() == 0) {
00270             of << "VECTORS el_field_mode_" << k << " float" << endl;
00271             //get the electric field  value for each point
00272             for (int i = 0; i < dim[2]; i ++) {
00273                 point.z = zcoords.at(i);
00274                 for (int j = 0; j < dim[1]; j ++) {
00275                     point.y = ycoords.at(j);
00276                     for (int h = 0; h < dim[0]; h ++) {
00277                         point.x = xcoords.at(h);
00278                         value = nedelec_mesh_.eval(point, *q);
00279                         of << value.x << " " << value.y << " " << value.z << endl;
00280                     }
00281                 }
00282             }
00283             of << endl;
00284         }
00285     } // for k
00286 
00287     // writing the magnetic field values into the vtk file
00288     for (int k = 0; k < Q_.NumVectors(); k++) {
00289         get_eigenpair(k, &lambda, q);
00290 
00291         if (comm_.MyPID() == 0) {
00292             of << "VECTORS mag_field_mode_" << k << " float" << endl;
00293             //get the magnetic field  value for each point
00294             for (int i = 0; i < dim[2]; i ++) {
00295                 point.z = zcoords.at(i);
00296                 for (int j = 0; j < dim[1]; j ++) {
00297                     point.y = ycoords.at(j);
00298                     for (int h = 0; h < dim[0]; h ++) {
00299                         point.x = xcoords.at(h);
00300                         value = nedelec_mesh_.eval_curl(point, *q);
00301                         of << value.x << " " << value.y << " " << value.z << endl;
00302                     }
00303                 }
00304             }
00305             of << endl;
00306         }
00307     } // for k
00308 
00309     if (comm_.MyPID() == 0)
00310         delete q;
00311 }
00312 
00313 void VtkExport::export_a_to_b(std::string file,
00314                               mesh::Vector3 a,
00315                               mesh::Vector3 b,
00316                               int nof_intervals)
00317 {
00318     double lambda;
00319     colarray::Vector<double>* q(0);
00320     mesh::Vector3 d = b - a;
00321     mesh::Vector3 p;
00322     mesh::Vector3 value;
00323     double spacing = 1.0 / static_cast<double>(nof_intervals);
00324     ofstream of;
00325 
00326     if (comm_.MyPID() == 0) {
00327         q = new colarray::Vector<double>(get_nof_dofs());
00328 
00329         of.open(file.c_str());
00330         assert(of.is_open());
00331         of.precision(precision_);
00332         of << "# vtk DataFile Version 2.0" << endl;
00333         of << "Generated using FemaXX" << endl;
00334         of << "ASCII" << endl << endl;
00335         of << "DATASET STRUCTURED_POINTS" << endl;
00336         of << "DIMENSIONS " << nof_intervals + 1 << " " << 1 << " " << 1 << endl;
00337         of << "ORIGIN 0 0 0" << endl;
00338         of << "SPACING " << spacing << " 1 1" << endl;
00339         of << endl << "POINT_DATA " << nof_intervals + 1 <<  endl << endl;
00340     }
00341 
00342     for (int k = 0; k < get_nof_vectors(); k ++) {
00343         get_eigenpair(k, &lambda, q);
00344         if (comm_.MyPID() == 0) {
00345             of << "VECTORS el_field_mode_" << k << " float" << endl;
00346             for (int i = 0; i <= nof_intervals; i ++) {
00347                 p = a + d * (static_cast<double>(i) * spacing);
00348                 value = nedelec_mesh_.eval(p, *q);
00349                 of << value.x << " " << value.y << " " << value.z << endl;
00350             } // for i
00351             of << endl;
00352         }
00353     } // for k
00354 
00355     for (int k = 0; k < get_nof_vectors(); k ++) {
00356         get_eigenpair(k, &lambda, q);
00357         if (comm_.MyPID() == 0) {
00358             of << "VECTORS mag_field_mode_" << k << " float" << endl;
00359             for (int i = 0; i <= nof_intervals; i ++) {
00360                 p = a + d * (static_cast<double>(i) * spacing);
00361                 value = nedelec_mesh_.eval_curl(p, *q);
00362                 of << value.x << " " << value.y << " " << value.z << endl;
00363             } // for i
00364             of << endl;
00365         }
00366     } // for k
00367 
00368     if (comm_.MyPID() == 0)
00369         delete q;
00370 }
00371 
00372 void VtkExport::get_eigenpair(int k, double* lambda, colarray::Vector<double>* q) {
00373     assert(0 <= k && k < get_nof_vectors());
00374     assert(comm_.MyPID() > 0 || q->_n == get_nof_dofs());
00375 
00376     //
00377     // Bring k-th eigenvector to processor 0
00378     //
00379 
00380     // create a target map, for which all the elements are on proc 0
00381     Epetra_MultiVector q_dist(View, Q_, k, 1);
00382     int NumMyElements_target = 0;
00383     if (comm_.MyPID() == 0)
00384         NumMyElements_target = Q_.GlobalLength();
00385     Epetra_Map TargetMap(-1, NumMyElements_target, 0, comm_);
00386 
00387     // redistribute vector to processor 0
00388     Epetra_Vector q0(TargetMap);
00389     q0.Export(q_dist, Epetra_Export(q_dist.Map(), TargetMap), Add);
00390 
00391     // store eigenpair on processor 0
00392     if (comm_.MyPID() == 0) {
00393         *lambda = lambda_[k];
00394         for (int i = 0; i < get_nof_dofs(); i ++) {
00395             q->operator()(i) =  q0[i];
00396         }
00397     }
00398 }
00399 
00400 void VtkExport::getdimensions(int nof_cells, mesh::Vector3& box_len, int dim[3]) {
00401 
00402     cout << endl << "   Vtk Rectilinear grid: bounding box length l, gridsize, grid dimensions nx ny nz   "
00403          << endl
00404          <<         "-------------------------------------------------------------------------------------"
00405          << endl;
00406     cout << endl <<  "l.x = " << box_len.x << ", l.y = " << box_len.y << ", l.z = " << box_len.z;
00407 
00408     mesh::Vector3 d;
00409     d.x = ::pow(static_cast<double>(nof_cells) * box_len.x * box_len.x / box_len.y / box_len.z, 1.0/3.0);
00410     d.y = d.x * box_len.y / box_len.x;
00411     d.z = d.x * box_len.z / box_len.x ;
00412 
00413     // make sure that all components are at least 0.5
00414     for (int i = 0; i < 3; i ++)
00415         if (d[i] < 0.5)
00416             d *= 0.5/d[i];
00417 
00418     dim[0] = static_cast<int>(round(d.x));
00419     dim[1] = static_cast<int>(round(d.y));
00420     dim[2] = static_cast<int>(round(d.z));
00421 
00422     cout << endl << "vtk_gridsize = " << nof_cells << endl;
00423     cout << "nx = " << dim[0] << ", " << "ny = " << dim[1]
00424          << ", " << " nz = " << dim[2] << endl;
00425     cout << "nx * ny * nz = " << dim[0]*dim[1]*dim[2] << endl;
00426     cout << "lx/nx = " << box_len.x / dim[0] << ", ly/ny =" << box_len.y/ dim[1] <<
00427         ", lz/nz = " << box_len.z / dim[2] << endl;
00428 }
00429 
00430 class UnitLine {
00431 public:
00432     mesh::Vector3 operator()(double x) const {
00433         return mesh::Vector3(x, 0, 0);
00434     }
00435 };
00436 
00437 void VtkExport::integrate_a_to_b(int mode_id,
00438                                  mesh::Vector3 a,
00439                                  mesh::Vector3 b,
00440                                  double& int_e)
00441 {
00442     double lambda;
00443     colarray::Vector<double>* q(0);
00444 
00445     if (comm_.MyPID() == 0) {
00446         q = new colarray::Vector<double>(get_nof_dofs());
00447     }
00448     // bring eigenvector to proc 0
00449     get_eigenpair(mode_id, &lambda, q);
00450 
00451     if (comm_.MyPID() == 0) {
00452         mesh::Vector3 d = b - a;
00453         double len = d.length();
00454         d /= len;
00455         EfieldIntegrand romberg_integrand(a, d, nedelec_mesh_, *q);
00456         int info;
00457         int_e = NR::qromb<EfieldIntegrand>(romberg_integrand, 0, len, 1e-10, info);
00458     }
00459 
00460     if (comm_.MyPID() == 0)
00461         delete q;
00462 
00463     // debug code following...
00464     class UnitLine line;
00465     mesh::TetMesh* tmesh = nedelec_mesh_.getTetMesh();
00466     tmesh->find_boundary<UnitLine>(line, 0.0, 1.0);
00467 }
00468 
00469 }; // namespace

Generated on Fri Oct 26 13:35:13 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7