00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
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
00137
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
00166 of << "VECTORS el_field_mode_" << k << " float" << endl;
00167
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
00182 of << "VECTORS mag_field_mode_" << k << " float" << endl;
00183
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
00206 TetMesh* mesh = get_mesh4();
00207 mesh->log_mesh_info();
00208 mesh->construct_octree();
00209
00210
00211 FemaxMesh femax_mesh(*mesh, 1, 0x0);
00212
00213
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
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
00253 TetMesh* mesh = get_mesh4();
00254 mesh->log_mesh_info();
00255 mesh->construct_octree();
00256
00257
00258 FemaxMesh femax_mesh(*mesh, 1, 0x0);
00259
00260
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
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
00298 TetMesh* mesh = get_mesh4();
00299 mesh->log_mesh_info();
00300 mesh->construct_octree();
00301
00302
00303 FemaxMesh femax_mesh(*mesh, 1, 0x0);
00304
00305
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
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
00345 const double b = 6.0;
00346
00347
00348 TetMesh* mesh = get_mesh4();
00349 mesh->log_mesh_info();
00350 mesh->construct_octree();
00351
00352
00353 FemaxMesh femax_mesh(*mesh, 1, 0x0);
00354
00355
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
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
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
00419 TetMesh* mesh = get_mesh4();
00420 mesh->log_mesh_info();
00421 mesh->construct_octree();
00422
00423
00424 FemaxMesh femax_mesh(*mesh, 1, 0x0);
00425
00426
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
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
00444 TetMesh* mesh = get_mesh4();
00445 mesh->log_mesh_info();
00446 mesh->construct_octree();
00447
00448
00449 FemaxMesh femax_mesh(*mesh, 1, 0x0);
00450
00451
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
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 }