00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
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;
00059 of << "generated using VtkExport::export_eigenfields" << endl;
00060 of << "ASCII" << endl << endl;
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
00072
00073 for (int i = 0; i < nofTets; i ++) {
00074 curTet = tmesh->get_tet(i);
00075 of << 4;
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
00084 of << "CELL_TYPES " << nofTets << endl;
00085 for (int i = 0; i < nofTets; i++) {
00086 of << 10 << endl;
00087 }
00088 of << endl;
00089
00090
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 }
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;
00123 of << "generated using VtkExport::export_eigenfields" << endl;
00124 of << "ASCII" << endl << endl;
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
00136
00137 for (int i = 0; i < nofTets; i ++) {
00138 curTet = tmesh->get_tet(i);
00139 of << 4;
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
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
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
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 }
00171
00172
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
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 }
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
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
00236
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
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
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 }
00286
00287
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
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 }
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 }
00351 of << endl;
00352 }
00353 }
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 }
00364 of << endl;
00365 }
00366 }
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
00378
00379
00380
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
00388 Epetra_Vector q0(TargetMap);
00389 q0.Export(q_dist, Epetra_Export(q_dist.Map(), TargetMap), Add);
00390
00391
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
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
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
00464 class UnitLine line;
00465 mesh::TetMesh* tmesh = nedelec_mesh_.getTetMesh();
00466 tmesh->find_boundary<UnitLine>(line, 0.0, 1.0);
00467 }
00468
00469 };