test/test_field_interp.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 <stdlib.h>
00014 #include <iostream>
00015 #include <fstream>
00016 #include <sstream>
00017 #include <unistd.h>
00018 #include <tut.h>
00019 #include "myrlog.h"
00020 #include "tetmesh/point.h"
00021 #include "tetmesh/tet.h"
00022 #include "colarray/vector.h"
00023 #include "looseoctree.h"
00024 #include "tetmeshbuilder.h"
00025 #include "broadcastmeshbuilder.h"
00026 #include "netgenreader.h"
00027 #include "femaxmesh.h"
00028 #include "LinearFieldInterpolation.h"
00029 
00030 namespace tut {
00031     Epetra_Comm& get_comm_world();
00032 }
00033 
00034 namespace
00035 {
00036     using namespace mesh;
00037     using namespace postprocess;
00038 
00039     struct test_data
00040     {
00041         test_data()
00042         {
00043         }
00044         static mesh::TetMesh* tmesh_;
00045         static FemaxMesh* femax_mesh_;
00046     };
00047 
00048     mesh::TetMesh* test_data::tmesh_ = 0;
00049     FemaxMesh* test_data::femax_mesh_ = 0;
00050     
00051     typedef tut::test_group<test_data> testgroup;
00052     typedef testgroup::object testobject;
00053     testgroup group("field_interp");
00054 
00058     template<>
00059     template<>
00060     void testobject::test<1>()
00061     {
00062         // Read mesh data structure
00063         std::string mesh_file_name = "test/octree/cop40k.ng";
00064         rAssert(!access(mesh_file_name.c_str(), R_OK));
00065 
00066         std::ifstream istr(mesh_file_name.c_str());
00067 
00068         Epetra_Comm& comm(tut::get_comm_world());
00069         TetMeshBuilder* builder = new TetMeshBuilder();
00070 #ifdef HAVE_MPI
00071         BroadcastMeshBuilder* broadcast = new BroadcastMeshBuilder(&comm, builder);
00072         NetgenReader* reader(0);
00073         if (comm.MyPID() == 0) {
00074             reader = new NetgenReader(&istr, broadcast);
00075             reader->read();
00076         } else
00077             broadcast->receiver();
00078         delete broadcast;
00079         if (comm.MyPID() == 0)
00080             delete reader;
00081 #else
00082         NetgenReader* reader = new NetgenReader(&istr, builder);
00083         reader->read();
00084         delete reader;
00085 #endif
00086         tmesh_ = builder->get_mesh();
00087         delete builder;
00088         tmesh_->log_mesh_info();
00089 
00090         // Construct Octree
00091         tmesh_->construct_octree();
00092         
00093         // Construct FemaxMesh
00094         const int element_order = 1;
00095         const int sym_plane_config = 6;
00096         femax_mesh_ = new FemaxMesh(*tmesh_, element_order, sym_plane_config);
00097     }
00098 
00102     template<>
00103     template<>
00104     void testobject::test<10>()
00105     {
00106         const double diff_tol = 1e-14;
00107     
00108         // Create fake random eigenvector
00109         colarray::Vector<double> q(femax_mesh_->get_nedelec_mesh().get_num_global_mdofs(0));
00110         for (int i = 0; i < q._n; ++ i)
00111             q(i) = static_cast<double>(random()) / RAND_MAX;
00112         
00113         // Create LinearFieldInterpolation
00114         LinearFieldInterpolation eval_interp(femax_mesh_->get_nedelec_mesh(), q);
00115         
00116         // Evaluate and compare field at mesh points    
00117         for (TetMesh::point_iterator it = tmesh_->point_begin(); it != tmesh_->point_end(); ++ it) {
00118             Point& point = *it;
00119             Vector3 f1 = eval_interp.eval(point.get_coord());
00120             Vector3 f2 = femax_mesh_->get_nedelec_mesh().eval(point.get_coord(), q);
00121             Vector3 d = f1 - f2;
00122             double diff = d.length();
00123             ensure(diff < diff_tol);
00124         }
00125     }
00126 
00130     template<>
00131     template<>
00132     void testobject::test<11>()
00133     {
00134         const double diff_tol = 1e-12;
00135     
00136         // Create fake random eigenvector
00137         colarray::Vector<double> q(femax_mesh_->get_nedelec_mesh().get_num_global_mdofs(0));
00138         for (int i = 0; i < q._n; ++ i)
00139             q(i) = static_cast<double>(random()) / RAND_MAX;
00140         
00141         // Create LinearFieldInterpolation
00142         LinearFieldInterpolation eval_interp(femax_mesh_->get_nedelec_mesh(), q);
00143         
00144         // Evaluate and compare field at mesh points    
00145         for (TetMesh::tet_iterator it = tmesh_->tet_begin(); it != tmesh_->tet_end(); ++ it) {
00146             Tet* tet = &*it;
00147             // Interpolate field at mid-point of Tet
00148             Vector3 mid_coord(Vector3::ZERO);
00149             for (int i = 0; i < 4; ++ i)
00150                 mid_coord += tet->get_corner(i)->get_coord();
00151             mid_coord /= 4.0;
00152             Vector3 f1 = eval_interp.eval(mid_coord);
00153             // Build average of field evaluated at the four corners
00154             Vector3 f2(Vector3::ZERO);
00155             for (int i = 0; i < 4; ++ i)
00156                 f2 += eval_interp.eval(tet->get_corner(i)->get_coord());
00157             f2 /= 4.0;
00158             // Compare
00159             Vector3 d = f1 - f2;
00160             double diff = d.length();
00161             ensure(diff < diff_tol);
00162         }
00163     }
00164 
00168     template<>
00169     template<>
00170     void testobject::test<50>()
00171     {
00172         delete femax_mesh_;
00173         delete tmesh_;
00174     }
00175 
00176 }

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