test/test_field_integrate.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: %{MODULE}
00003 //
00004 // Description:
00005 //
00006 //
00007 // Author: %{AUTHOR} <%{EMAIL}>, (C) %{YEAR}
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
00010 //
00011 //
00012 
00013 #include <fstream>
00014 #include <tut.h>
00015 #include "myrlog.h"
00016 #include "point.h"
00017 #include "tet.h"
00018 #include "looseoctree.h"
00019 #include "tetmeshbuilder.h"
00020 #include "femaxmesh.h"
00021 #include "LineAtoB.h"
00022 #include "CircleZ.h"
00023 #include "EfieldIntegrand.h"
00024 #include "LoggingFunction.h"
00025 #include "adaptsim.h"
00026 #include "romberg.h"
00027 #include "Epetra_Comm.h"
00028 
00029 namespace tut {
00030     Epetra_Comm& get_comm_world();
00031 }
00032 
00033 namespace
00034 {
00035     using namespace mesh;
00036 
00037     struct test_data
00038     {
00039         test_data()
00040         {
00041         }
00042     };
00043     
00044     typedef tut::test_group<test_data> testgroup;
00045     typedef testgroup::object testobject;
00046     testgroup group("field_integrate");
00047 
00049     mesh::TetMesh* get_mesh4() {
00050         // Generate Tetmesh
00051         TetMeshBuilder* builder = new TetMeshBuilder();
00052         builder->init_coord(6);
00053         builder->set_coord(0,  1,  0,  0);
00054         builder->set_coord(1,  0,  1,  0);
00055         builder->set_coord(2, -1,  0,  0);
00056         builder->set_coord(3,  0, -1,  0);
00057         builder->set_coord(4,  0,  0,  1);
00058         builder->set_coord(5,  0,  0, -1);
00059         builder->finalize_coord();
00060         builder->init_tet(4);
00061         builder->set_tet(0, 0, 1, 3, 4);
00062         builder->set_tet(1, 1, 2, 3, 4);
00063         builder->set_tet(2, 0, 1, 3, 5);
00064         builder->set_tet(3, 1, 2, 3, 5);
00065         builder->finalize_tet();
00066         builder->init_bc(8);
00067         builder->set_bc(0, 1, 4, 0);
00068         builder->set_bc(0, 1, 5, 0);
00069         builder->set_bc(0, 3, 4, 0);
00070         builder->set_bc(0, 3, 5, 0);
00071         builder->set_bc(1, 2, 4, 0);
00072         builder->set_bc(1, 2, 5, 0);
00073         builder->set_bc(2, 3, 4, 0);
00074         builder->set_bc(2, 3, 5, 0);
00075         builder->finalize_bc(0);
00076         TetMesh* mesh = builder->get_mesh();
00077         delete builder;
00078         return mesh;
00079     }
00080 
00081 void getdimensions(int nof_cells, mesh::Vector3& box_len, int dim[3]) {
00082     mesh::Vector3 d;
00083     d.x = ::pow(static_cast<double>(nof_cells) * box_len.x * box_len.x / box_len.y / box_len.z, 1.0/3.0);
00084     d.y = d.x * box_len.y / box_len.x;
00085     d.z = d.x * box_len.z / box_len.x ;
00086 
00087     // make sure that all components are at least 0.5
00088     for (int i = 0; i < 3; i ++)
00089         if (d[i] < 0.5)
00090             d *= 0.5/d[i];
00091 
00092     dim[0] = static_cast<int>(round(d.x));
00093     dim[1] = static_cast<int>(round(d.y));
00094     dim[2] = static_cast<int>(round(d.z));
00095 }
00096 
00097         
00098 void export_vtk(const FemaxMesh& femax_mesh, const colarray::Vector<double>& q, const std::string file, int nof_cells) {
00099     int nofPoints;
00100     int dim[3];
00101     mesh::TetMesh* tmesh;
00102     mesh::Vector3 center;
00103     mesh::Vector3 box_len;
00104     mesh::Vector3 value;
00105     tmesh = &(femax_mesh.get_tet_mesh());
00106 
00107     std::vector<double> xcoords, ycoords, zcoords;
00108     double curx, cury, curz;
00109     mesh::Vector3 point;
00110     mesh::Vector3 ref_point;
00111     std::vector<mesh::Tet *> tets;
00112     std::ofstream of;
00113     int k = 1;
00114 
00115         of.open(file.c_str());
00116         assert(of.is_open());
00117 
00118         of.precision(6);
00119         of << "# vtk DataFile Version 2.0" << endl;
00120         of << "generated using FemaXX" << endl;
00121         of << "ASCII" << endl << endl;
00122         of << "DATASET RECTILINEAR_GRID" << endl;
00123 
00124         // get the bounding box from tmesh and calculate the dimensions of the new grid
00125         tmesh->get_bounding_box(center, box_len);
00126         ref_point = center - box_len;
00127         box_len *= 2.0;
00128 
00129         getdimensions(nof_cells, box_len, dim);
00130         of << "DIMENSIONS " << dim[0] << " " << dim[1] << " " << dim[2] << endl;
00131 
00132         double xstep = box_len.x / dim[0];
00133         double ystep = box_len.y / dim[1];
00134         double zstep = box_len.z / dim[2];
00135 
00136         // create three vectors containing the x, y and z coordinates of the points
00137         // and do the insertion of these into the vtk file
00138         of << "X_COORDINATES " << dim[0] << " float" << endl;
00139         for (int i = 0; i <  dim[0]; i++) {
00140             curx = ref_point.x + (i * xstep);
00141             xcoords.push_back(curx);
00142             of << curx << " ";
00143         }
00144         of << endl;
00145 
00146         of << "Y_COORDINATES " << dim[1] << " float" << endl;
00147         for (int i = 0; i <  dim[1]; i++) {
00148             cury = ref_point.y + (i * ystep);
00149             ycoords.push_back(cury);
00150             of << cury << " ";
00151         }
00152         of << endl;
00153 
00154         of << "Z_COORDINATES " << dim[2] << " float" << endl;
00155         for (int i = 0; i <  dim[2]; i++) {
00156             curz = ref_point.z + (i * zstep);
00157             zcoords.push_back(curz);
00158             of << curz << " ";
00159         }
00160         of << endl;
00161 
00162         nofPoints = xcoords.size() *  ycoords.size() * zcoords.size() ;
00163         of << endl << "POINT_DATA " <<  nofPoints <<  endl << endl;
00164 
00165     // writing the electric field values into the vtk file
00166             of << "VECTORS el_field_mode_" << k << " float" << endl;
00167             //get the electric field  value for each point
00168             for (int i = 0; i < dim[2]; i ++) {
00169                 point.z = zcoords.at(i);
00170                 for (int j = 0; j < dim[1]; j ++) {
00171                     point.y = ycoords.at(j);
00172                     for (int h = 0; h < dim[0]; h ++) {
00173                         point.x = xcoords.at(h);
00174                         value = femax_mesh.get_nedelec_mesh().eval(point, q);
00175                         of << value.x << " " << value.y << " " << value.z << endl;
00176                     }
00177                 }
00178             }
00179             of << endl;
00180  
00181     // writing the magnetic field values into the vtk file
00182             of << "VECTORS mag_field_mode_" << k << " float" << endl;
00183             //get the magnetic field  value for each point
00184             for (int i = 0; i < dim[2]; i ++) {
00185                 point.z = zcoords.at(i);
00186                 for (int j = 0; j < dim[1]; j ++) {
00187                     point.y = ycoords.at(j);
00188                     for (int h = 0; h < dim[0]; h ++) {
00189                         point.x = xcoords.at(h);
00190                         value = femax_mesh.get_nedelec_mesh().eval_curl(point, q);
00191                         of << value.x << " " << value.y << " " << value.z << endl;
00192                     }
00193                 }
00194             }
00195             of << endl;
00196 }
00197 
00198     
00199     template<>
00200     template<>
00201     void testobject::test<1>()
00202     {
00203         const double tol = 1e-14;
00204         
00205         // Generate Tetmesh
00206         TetMesh* mesh = get_mesh4();
00207         mesh->log_mesh_info();
00208         mesh->construct_octree();
00209         
00210         // Create FemaxMesh
00211         FemaxMesh femax_mesh(*mesh, 1, 0x0);
00212         
00213         // Fake eigenvector
00214         colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
00215         ensure(q._n == 1);
00216         q(0) = 1.0;
00217         
00218         // Path with y==0 is orthogonal to E-field
00219         Vector3 A(0.2, 0.0, -0.5);
00220         Vector3 B(-0.2, 0.0,  0.5);
00221         LineAtoB line(A, B);
00222         EfieldDs<LineAtoB> integrand(line, femax_mesh.get_nedelec_mesh(), q);
00223 
00224         for (int i = 0; i <= 10; ++ i) {
00225             Vector3 x(line((double)i/10.0));
00226             Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
00227             rInfo("E(%f,%f,%f) = (%e,%e,%e)",
00228                   x.x, x.y, x.z,
00229                   E.x, E.y, E.z);
00230             ensure(fabs(E.x) < tol && fabs(E.z) < tol);
00231         }
00232       
00233         int info;
00234         double gap_voltage = NR::qromb(integrand, 0.0, 1.0, 1e-10, info);
00235         rInfo("qromb: voltage=%e", gap_voltage);
00236         ensure(fabs(gap_voltage) < tol);
00237 
00238         gap_voltage = NR::adaptsim(integrand, 0.0, 1.0, 1e-10, info);
00239         rInfo("adaptsim: voltage=%e, info=%d", gap_voltage, info);
00240         ensure(info == 0);
00241         ensure(fabs(gap_voltage) < tol);
00242 
00243         delete mesh;
00244     }
00245     
00246     template<>
00247     template<>
00248     void testobject::test<2>()
00249     {
00250         const double tol = 1e-14;
00251         
00252         // Generate Tetmesh
00253         TetMesh* mesh = get_mesh4();
00254         mesh->log_mesh_info();
00255         mesh->construct_octree();
00256         
00257         // Create FemaxMesh
00258         FemaxMesh femax_mesh(*mesh, 1, 0x0);
00259         
00260         // Fake eigenvector
00261         colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
00262         ensure(q._n == 1);
00263         q(0) = 1.0;
00264         
00265         // Path with y==0 is orthogonal to E-field
00266         Vector3 A(0.00001, 0.2, -0.5);
00267         Vector3 B(0.00001, -0.2,  0.5);
00268         LineAtoB line(A, B);
00269         EfieldDs<LineAtoB> integrand(line, femax_mesh.get_nedelec_mesh(), q);
00270 
00271         for (int i = 0; i <= 10; ++ i) {
00272             Vector3 x(line((double)i/10.0));
00273             Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
00274             rInfo("E(%f,%f,%f) = (%e,%e,%e)",
00275                   x.x, x.y, x.z,
00276                   E.x, E.y, E.z);
00277         }
00278       
00279         int info;
00280         double gap_voltage_1 = NR::qromb(integrand, 0.0, 1.0, 1e-10, info);
00281         rInfo("qromb: voltage=%e", gap_voltage_1);
00282 
00283         double gap_voltage_2 = NR::adaptsim(integrand, 0.0, 1.0, 1e-10, info);
00284         rInfo("adaptsim: voltage=%e, info=%d", gap_voltage_2, info);
00285         ensure(info == 0);
00286         ensure(fabs(gap_voltage_2 - gap_voltage_1) < tol);
00287 
00288         delete mesh;
00289     }
00290 
00291     template<>
00292     template<>
00293     void testobject::test<3>()
00294     {
00295         const double tol = 1e-14;
00296         
00297         // Generate Tetmesh
00298         TetMesh* mesh = get_mesh4();
00299         mesh->log_mesh_info();
00300         mesh->construct_octree();
00301         
00302         // Create FemaxMesh
00303         FemaxMesh femax_mesh(*mesh, 1, 0x0);
00304         
00305         // Fake eigenvector
00306         colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
00307         ensure(q._n == 1);
00308         q(0) = 100.0;
00309         
00310         // Path with y==0 is orthogonal to E-field
00311         Vector3 A(-0.0005, -0.3, 0.0001);
00312         Vector3 B(0.0005, 0.5, 0.0001);
00313         LineAtoB line(A, B);
00314         EfieldDs<LineAtoB> integrand(line, femax_mesh.get_nedelec_mesh(), q);
00315 
00316         for (int i = 0; i <= 10; ++ i) {
00317             double t = (double)i/10.0;
00318             Vector3 x(line(t));
00319             Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
00320             rInfo("E(%f,%f,%f), Eds = (%e,%e,%e), %e",
00321                   x.x, x.y, x.z,
00322                   E.x, E.y, E.z,
00323                   integrand(t));
00324         }
00325       
00326         int info;
00327         double gap_voltage_1 = NR::qromb(integrand, 0.0, 1.0, 1e-15, info);
00328         rInfo("qromb: voltage=%e", gap_voltage_1);
00329 
00330         double gap_voltage_2 = NR::adaptsim(integrand, 0.0, 1.0, 1e-15, info);
00331         rInfo("adaptsim: voltage=%e, info=%d", gap_voltage_2, info);
00332         ensure(info == 0);
00333         ensure(fabs(gap_voltage_2 - gap_voltage_1) < tol);
00334 
00335         delete mesh;
00336     }
00337 
00338     template<>
00339     template<>
00340     void testobject::test<10>()
00341     {
00342         const double int_tol = 2e-8;
00343         const double a = 0.0;
00344         //const double b = 4.0*acos(0.0);
00345         const double b = 6.0;
00346         
00347         // Generate Tetmesh
00348         TetMesh* mesh = get_mesh4();
00349         mesh->log_mesh_info();
00350         mesh->construct_octree();
00351         
00352         // Create FemaxMesh
00353         FemaxMesh femax_mesh(*mesh, 1, 0x0);
00354         
00355         // Fake eigenvector
00356         colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
00357         ensure(q._n == 1);
00358         q(0) = 1.0;
00359         
00360         // Circular path
00361         CircleZ curve(0.5, 0.01);
00362         EfieldDs<CircleZ> integrand(curve, femax_mesh.get_nedelec_mesh(), q);
00363         LoggingFunction<EfieldDs<CircleZ>, double> my_integrand(integrand);
00364 
00365         const int max_steps = 100;
00366         for (int i = 0; i <= max_steps; ++ i) {
00367             double t = (b-a) * (double)i/max_steps + a;
00368             Vector3 x(curve(t));
00369             Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
00370             rInfo("E(%f,%f,%f), Eds = (%e,%e,%e), %e",
00371                   x.x, x.y, x.z,
00372                   E.x, E.y, E.z,
00373                   integrand(t));
00374         }
00375       
00376         int info;
00377         double gap_voltage_1 = NR::qromb(my_integrand, a, b, int_tol, info);
00378         rInfo("qromb: voltage=%e", gap_voltage_1);
00379         rInfo("num_eval=%d, y_min=%e, y_max=%e", 
00380               my_integrand.get_num_eval(), 
00381               my_integrand.get_y_min(), 
00382               my_integrand.get_y_max());
00383         ensure(info == 0);
00384 
00385         my_integrand.reset();
00386         if (tut::get_comm_world().MyPID() == 0)
00387             my_integrand.enable_plot("adaptsim.tab");
00388         double gap_voltage_2 = NR::adaptsim(my_integrand, a, b, int_tol, info);
00389         rInfo("adaptsim: voltage=%e, info=%d", gap_voltage_2, info);
00390         rInfo("num_eval=%d, y_min=%e, y_max=%e", 
00391               my_integrand.get_num_eval(), 
00392               my_integrand.get_y_min(), 
00393               my_integrand.get_y_max());
00394         ensure(info == 0);
00395         
00396         my_integrand.reset();
00397         double gap_voltage_3 = NR::qtrap(my_integrand, a, b, int_tol, info);
00398         rInfo("qtrap: voltage=%e, info=%d", gap_voltage_3, info);
00399         rInfo("num_eval=%d, y_min=%e, y_max=%e", 
00400               my_integrand.get_num_eval(), 
00401               my_integrand.get_y_min(), 
00402               my_integrand.get_y_max());
00403         ensure(info == 0);
00404         
00405         if (!(fabs(gap_voltage_2 - gap_voltage_1) < 1e3*int_tol))
00406             rWarning("Results from qromb and adaptsim differ. qromb returns too early.");
00407         // ensure(fabs(gap_voltage_2 - gap_voltage_1) < 1e3*int_tol);
00408         ensure(fabs(gap_voltage_3 - gap_voltage_2) < 1e3*int_tol);
00409 
00410         delete mesh;
00411     }    
00412     
00414     template<>
00415     template<>
00416     void testobject::test<-1>()
00417     {
00418         // Generate Tetmesh
00419         TetMesh* mesh = get_mesh4();
00420         mesh->log_mesh_info();
00421         mesh->construct_octree();
00422         
00423         // Create FemaxMesh
00424         FemaxMesh femax_mesh(*mesh, 1, 0x0);
00425         
00426         // Fake eigenvector
00427         colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
00428         ensure(q._n == 1);
00429         q(0) = 1.0;
00430         
00431         // Export mesh/field to vtk file
00432         if (tut::get_comm_world().MyPID() == 0)
00433             export_vtk(femax_mesh, q, "test_fieldintegrate.vtk", 1000);
00434         
00435         delete mesh;
00436     }
00437     
00439     template<>
00440     template<>
00441     void testobject::test<-2>()
00442     {
00443         // Generate Tetmesh
00444         TetMesh* mesh = get_mesh4();
00445         mesh->log_mesh_info();
00446         mesh->construct_octree();
00447         
00448         // Create FemaxMesh
00449         FemaxMesh femax_mesh(*mesh, 1, 0x0);
00450         
00451         // Fake eigenvector
00452         colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
00453         ensure(q._n == 1);
00454         q(0) = 1.0;
00455         
00456         // Circular path
00457         CircleZ curve(0.5, 0.01);
00458         EfieldDs<CircleZ> integrand(curve, femax_mesh.get_nedelec_mesh(), q);
00459 
00460         for (double int_tol = 1.0; int_tol >= 1e-15; int_tol /= 10.0) {
00461             int info;
00462             double gap_voltage = NR::qromb(integrand, 0.0, 6.0, int_tol, info);
00463             rWarning("qromb: tol=%e, voltage=%.12e, info=%d", int_tol, gap_voltage, info);
00464             ensure(info == 0);
00465         }
00466         
00467         for (double int_tol = 1.0; int_tol >= 1e-15; int_tol /= 10.0) {
00468             int info;
00469             double gap_voltage = NR::qtrap(integrand, 0.0, 6.0, int_tol, info);
00470             rWarning("qtrap: tol=%e, voltage=%.12e, info=%d", int_tol, gap_voltage, info);
00471             ensure(info == 0);
00472         }
00473         
00474         for (double int_tol = 1.0; int_tol >= 1e-15; int_tol /= 10.0) {
00475             int info;
00476             double gap_voltage = NR::adaptsim(integrand, 0.0, 6.0, int_tol, info);
00477             rWarning("adaptsim: tol=%e, voltage=%.12e, info=%d", int_tol, gap_voltage, info);
00478             ensure(info == 0);
00479         }
00480         delete mesh;
00481     }
00482 
00483 }

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