00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
00091 tmesh_->construct_octree();
00092
00093
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
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
00114 LinearFieldInterpolation eval_interp(femax_mesh_->get_nedelec_mesh(), q);
00115
00116
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
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
00142 LinearFieldInterpolation eval_interp(femax_mesh_->get_nedelec_mesh(), q);
00143
00144
00145 for (TetMesh::tet_iterator it = tmesh_->tet_begin(); it != tmesh_->tet_end(); ++ it) {
00146 Tet* tet = &*it;
00147
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
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
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 }