Branch data Line data Source code
1 : : /**
2 : : * \file mesh.cxx \brief Mesh Implementation
3 : : *
4 : : * This source files implements the different methods for the mesh class in order to be able to read
5 : : * the mesh data from the files
6 : : */
7 : :
8 : : #include "read_mesh.h"
9 : : #include "file_parser.h"
10 : :
11 : : // mesh Class
12 : :
13 : : //! Default Constructor
14 : 0 : mesh::mesh()
15 : : {
16 : 0 : _number_of_points = -1;
17 : 0 : _number_of_tetrahedra = -1;
18 : 0 : _number_of_materials = -1;
19 : :
20 : 0 : _index_of_refraction = -1;
21 : 0 : }
22 : :
23 : : /**
24 : : * Regular Constructor
25 : : * @param mesh_file [in]: Path to the mesh file
26 : : * @param fparser [in]: Pointer to the file parser to read the optical properties
27 : : * @param read_from_vtk [in]: Selects whether read from .mesh file or .vtk file
28 : : */
29 : 3248 : mesh::mesh(string mesh_file, Parser* fparser, bool read_from_vtk)
30 : : {
31 : : // Reading data
32 [ + - ]: 1624 : read_opt_props(fparser);
33 [ + + ]: 1624 : if (read_from_vtk) {
34 [ + - ][ + - ]: 728 : read_vtk(mesh_file);
35 : : }
36 : : else {
37 [ + - ][ + - ]: 896 : read_mesh(mesh_file); //opt_file_ss.str());
38 : : }
39 : 1624 : }
40 : :
41 : : //! compute the volume of a teterahedron
42 : 392870565 : double mesh::compute_tet_volume(tetrahedron* t)
43 : : {
44 : : SP_Point *p0, *p1, *p2, *p3;
45 : :
46 [ + - ]: 392870565 : p0 = _list_of_points[t->get_p1_id() - 1];
47 [ + - ]: 392870565 : p1 = _list_of_points[t->get_p2_id() - 1];
48 [ + - ]: 392870565 : p2 = _list_of_points[t->get_p3_id() - 1];
49 [ + - ]: 392870565 : p3 = _list_of_points[t->get_p4_id() - 1];
50 : :
51 : :
52 : : // volume = (1/6) * (AB x AC) dot AD
53 [ + - ][ + - ]: 392870565 : double ab[3] = { p1->get_xcoord() - p0->get_xcoord(),
54 [ + - ][ + - ]: 392870565 : p1->get_ycoord() - p0->get_ycoord(),
55 [ + - ][ + - ]: 392870565 : p1->get_zcoord() - p0->get_zcoord() };
56 [ + - ][ + - ]: 392870565 : double ac[3] = { p2->get_xcoord() - p0->get_xcoord(),
57 [ + - ][ + - ]: 392870565 : p2->get_ycoord() - p0->get_ycoord(),
58 [ + - ][ + - ]: 392870565 : p2->get_zcoord() - p0->get_zcoord() };
59 [ + - ][ + - ]: 392870565 : double ad[3] = { p3->get_xcoord() - p0->get_xcoord(),
60 [ + - ][ + - ]: 392870565 : p3->get_ycoord() - p0->get_ycoord(),
61 [ + - ][ + - ]: 392870565 : p3->get_zcoord() - p0->get_zcoord() };
62 : 1178611695 : double volume = abs((1.0 / 6.0) * (ad[0] * (ab[1] * ac[2] - ab[2] * ac[1]) +
63 : 785741130 : ad[1] * (ab[2] * ac[0] - ab[0] * ac[2]) +
64 : 392870565 : ad[2] * (ab[0] * ac[1] - ab[1] * ac[0])));
65 : :
66 : 392870565 : return volume;
67 : : }
68 : :
69 : : //! Destructor
70 : 1736 : mesh::~mesh()
71 : : {
72 [ + + ]: 60957904 : for (unsigned int i = 0; i < _list_of_points.size(); i++)
73 [ + - ]: 60957036 : delete _list_of_points[i];
74 : :
75 [ + + ]: 367491236 : for (unsigned int i = 0; i < _list_of_tetrahedra.size(); i++)
76 [ + - ]: 367490368 : delete _list_of_tetrahedra[i];
77 : :
78 [ + + ]: 5208 : for (unsigned int i = 0; i < _list_of_materials.size(); i++)
79 [ + - ]: 4340 : delete _list_of_materials[i];
80 : 868 : }
81 : :
82 : : // Getters
83 : :
84 : : //! This method returns a list of pointers to all points in the mesh
85 : 1288480496 : vector<SP_Point*> & mesh::get_list_of_points()
86 : : {
87 : 1288480496 : return _list_of_points;
88 : : }
89 : :
90 : : //! This method returns a list of pointers to all tetrahedra in the mesh
91 : 1104374852 : vector<tetrahedron*> & mesh::get_list_of_tetrahedra()
92 : : {
93 : 1104374852 : return _list_of_tetrahedra;
94 : : }
95 : :
96 : : //! This method returns a list of pointers to all materials in the mesh
97 : 15137224 : vector<material_opt*> & mesh::get_list_of_materials()
98 : : {
99 : 15137224 : return _list_of_materials;
100 : : }
101 : :
102 : : //! This method returns the number of points in the mesh
103 : 2128 : unsigned int mesh::get_number_of_points()
104 : : {
105 : 2128 : return _number_of_points;
106 : : }
107 : :
108 : : //! This method returns the mesh size in terms of number of tetrahedra in the mesh
109 : 3752 : unsigned int mesh::get_mesh_size()
110 : : {
111 : 3752 : return _number_of_tetrahedra;
112 : : }
113 : :
114 : : //! This method returns the number of materials in the mesh
115 : 10920 : unsigned int mesh::get_number_of_materials()
116 : : {
117 : 10920 : return _number_of_materials;
118 : : }
119 : :
120 : : //! This method returns the refraction index in the current mesh
121 : 0 : double mesh::get_index_of_refraction()
122 : : {
123 : 0 : return _index_of_refraction;
124 : : }
125 : :
126 : : //! Setters
127 : :
128 : : //!
129 : : //! This method adds a point to the list of points
130 : : //! @param p [in]: pointer to the point to add
131 : : //!
132 : 114048648 : void mesh::add_point(SP_Point *p)
133 : : {
134 : 114048648 : _list_of_points.push_back(p);
135 : 114048648 : }
136 : :
137 : : //!
138 : : //! This method adds a tetrahedron to the list of tetrahedra
139 : : //! @param t [in]: pointer to the tetrahedron to add
140 : : //!
141 : 687562624 : void mesh::add_tetrahedron(tetrahedron *t)
142 : : {
143 : 687562624 : _list_of_tetrahedra.push_back(t);
144 : 687562624 : }
145 : :
146 : : //!
147 : : //! This method adds a new material to the list of materials
148 : : //! @param m [in]: pointer to the new material to add
149 : : //!
150 : 8120 : void mesh::add_material(material_opt *m)
151 : : {
152 : 8120 : _list_of_materials.push_back(m);
153 : 8120 : }
154 : :
155 : : //!
156 : : //! This method sets the number of points in the mesh
157 : : //! @param np [in]: number of points
158 : : //!
159 : 0 : void mesh::set_number_of_points(unsigned int np)
160 : : {
161 : 0 : _number_of_points = np;
162 : 0 : }
163 : :
164 : : //!
165 : : //! This method sets the number of tetrahedra in the mesh
166 : : //! @param ms [in]: mesh size in terms of number of tetrahedra
167 : : //!
168 : 0 : void mesh::set_mesh_size(unsigned int ms)
169 : : {
170 : 0 : _number_of_tetrahedra = ms;
171 : 0 : }
172 : :
173 : : //!
174 : : //! This method sets the number of materials in the mesh
175 : : //! @param nm [in]: number of materials
176 : : //!
177 : 0 : void mesh::set_number_of_materials(unsigned int nm)
178 : : {
179 : 0 : _number_of_materials = nm;
180 : 0 : }
181 : :
182 : : //! Helper methods
183 : :
184 : : //!
185 : : //! This method converts the optical properties read by file parser to pdt-space optical properties
186 : : //! @param fparser [in]: Pointer to the file parser
187 : : //!
188 : 1624 : void mesh::read_opt_props(Parser* fparser)
189 : : {
190 : : // Reading the optical properties
191 [ + - ]: 1624 : fprintf(stderr, "\033[1;%dmSetting optical properties....\033[0m\n", 33);
192 : :
193 [ + - ]: 3248 : const vector<Parser_material> materials = fparser->get_materials();
194 : :
195 [ - + ]: 1624 : if (materials.size() == 0)
196 : : {
197 [ # # ]: 0 : fprintf(stderr, "\033[1;%dmNo materials read.\nExiting...\033[0m\n", 31);
198 : 0 : exit(-1);
199 : : }
200 : :
201 : 1624 : _number_of_materials = (unsigned) (materials.size()) - 1; // exterior does not count
202 : :
203 : : // Reading the optical properties for each material
204 : : material_opt * tmp_material;
205 [ + + ]: 9744 : for (unsigned int i = 1; i <= _number_of_materials; i++)
206 : : {
207 : 8120 : tmp_material = new material_opt(materials[i].id, materials[i].mu_a,
208 [ + - ][ + - ]: 8120 : materials[i].mu_s, materials[i].g, materials[i].n, this);
209 [ + - ]: 8120 : add_material(tmp_material);
210 : : }
211 : :
212 : : // PDT-SPACE assumes exterior to be air
213 : 1624 : _index_of_refraction = materials[0].n;
214 : 1624 : }
215 : :
216 : : //!
217 : : //! This method reads the data files and constructs the mesh
218 : : //! @param mesh_file [in]: file name of the mesh data (points and tetrahedra)
219 : : //! @param opt_file [in]: Metrials' properties data (optical properties)
220 : : //!
221 : 896 : void mesh::read_mesh(string mesh_file)
222 : : {
223 [ + - ]: 1792 : ifstream input;
224 : 896 : char str[255] = "";
225 : :
226 : : // Reading mesh data
227 [ + - ]: 896 : fprintf(stderr, "\033[1;%dmReading mesh file....\033[0m\n", 33);
228 [ + - ]: 896 : input.open(mesh_file);
229 [ + - ][ - + ]: 896 : if (!input)
230 : : {
231 [ # # ]: 0 : fprintf(stderr, "\033[1;%dmCould not open %s\nExiting...\033[0m\n", 31, mesh_file.c_str());
232 : 0 : exit(-1);
233 : : }
234 : :
235 [ + + ]: 2688 : for (int i = 0; i < 2; i++)
236 : : {
237 [ + - ]: 1792 : input.getline(str, 255);
238 [ + - ][ + - ]: 3584 : stringstream line(str);
239 : :
240 [ + + ]: 1792 : if (i == 0)
241 [ + - ]: 896 : line >> _number_of_points;
242 : : else
243 [ + - ]: 896 : line >> _number_of_tetrahedra;
244 : : }
245 : :
246 : : SP_Point * tmp_p;
247 : : float xcoord, ycoord, zcoord;
248 [ + + ]: 62924288 : for (unsigned int i = 0; i < _number_of_points; i++)
249 : : {
250 [ + - ]: 62923392 : input.getline(str, 255);
251 [ + - ][ + - ]: 125846784 : stringstream line(str);
252 : :
253 [ + - ][ + - ]: 62923392 : line >> xcoord >> ycoord >> zcoord;
[ + - ]
254 : :
255 [ + - ][ + - ]: 62923392 : tmp_p = new SP_Point(i+1, xcoord, ycoord, zcoord, this);
256 [ + - ]: 62923392 : add_point(tmp_p);
257 : : }
258 : :
259 : : tetrahedron* tmp_tet;
260 : : int id1, id2, id3, id4, mat_id;
261 [ + + ]: 379345792 : for(unsigned int i = 0; i < _number_of_tetrahedra; i++)
262 : : {
263 [ + - ]: 379344896 : input.getline(str, 255);
264 [ + - ][ + - ]: 758689792 : stringstream line(str);
265 : :
266 [ + - ][ + - ]: 379344896 : line >> id1 >> id2 >> id3 >> id4 >> mat_id;
[ + - ][ + - ]
[ + - ]
267 : :
268 [ + - ][ + - ]: 379344896 : tmp_tet = new tetrahedron(i+1, id1, id2, id3, id4, mat_id, this);
269 [ + - ]: 379344896 : add_tetrahedron(tmp_tet);
270 : : }
271 : :
272 [ + - ]: 896 : fprintf(stderr, "\033[1;%dmDone Reading files....\033[0m\n", 33);
273 : 896 : }
274 : :
275 : : //!
276 : : //! This method reads the vtk file and constructs the mesh
277 : : //! @param vtk_file [in]: file name of the mesh data (points and tetrahedra)
278 : : //!
279 : 728 : void mesh::read_vtk(string vtk_file)
280 : : {
281 : : // Reading from vtk file
282 [ + - ]: 728 : fprintf(stderr, "\033[1;%dmReading vtk file....\033[0m\n", 33);
283 : :
284 [ + - ]: 728 : vtkUnstructuredGridReader* reader = vtkUnstructuredGridReader::New();
285 [ + - ]: 728 : reader->vtkDataReader::SetFileName(vtk_file.c_str());
286 [ + - ]: 728 : reader->vtkDataReader::Update();
287 : :
288 [ + - ]: 728 : vtkUnstructuredGrid * vtk_input = reader->GetOutput();
289 : :
290 [ + - ]: 728 : _number_of_tetrahedra = (unsigned int) vtk_input->GetNumberOfCells();
291 [ + - ]: 728 : _number_of_points = (unsigned int) vtk_input->GetNumberOfPoints(); // vtkPoints
292 : :
293 : : // Converting VTK points into PDT-SPACE points
294 : : SP_Point * tmp_p;
295 : : double p_coord[3];
296 [ + + ]: 51125984 : for (unsigned int i = 0; i < _number_of_points; i++) {
297 [ + - ]: 51125256 : vtk_input->GetPoint(i, p_coord);
298 : :
299 [ + - ][ + - ]: 51125256 : tmp_p = new SP_Point(i+1, (float)(p_coord[0]), (float)(p_coord[1]), (float)(p_coord[2]), this);
300 [ + - ]: 51125256 : add_point(tmp_p);
301 : : }
302 : :
303 : : // Converting VTK cells into PDT-SPACE tetrahedra
304 : : tetrahedron * tmp_tet;
305 [ + - ]: 728 : vtkIdList* pt_ids = vtkIdList::New();
306 : :
307 : : // set region
308 : 728 : vtkCellData* cd = vtk_input->GetCellData();
309 : 728 : vtkDataArray* regions = nullptr;
310 : 728 : bool has_region_data = false;
311 : :
312 [ + - ][ + - ]: 728 : if (cd && (regions = cd->GetScalars("Region"))) {
[ + - ][ + - ]
313 : 728 : has_region_data = true;
314 [ # # ][ # # ]: 0 : } else if (cd && (regions = cd->GetScalars())) {
[ # # ][ # # ]
315 [ # # ]: 0 : const char *name = regions->GetName();
316 [ # # ]: 0 : fprintf(stderr, "\033[1;%dmRegion label not found, using data label %s instead.\033[0m\n", 33, name);
317 : 0 : has_region_data = true;
318 : : }
319 : : int id[4];
320 [ + + ]: 308218456 : for (unsigned int i = 0; i < _number_of_tetrahedra; i++) {
321 [ + - ]: 308217728 : vtk_input->GetCellPoints(i, pt_ids);
322 [ + + ]: 1541088640 : for (unsigned int j = 0; j < 4; j++) {
323 : 1232870912 : id[j] = (int) pt_ids->GetId(j);
324 : : }
325 : :
326 : 308217728 : int region = 0;
327 [ + - ]: 308217728 : if (has_region_data) {
328 [ + - ]: 308217728 : region = (int) regions->GetTuple1(i);
329 : : }
330 : :
331 [ + - ][ + - ]: 308217728 : tmp_tet = new tetrahedron(i+1, id[0]+1, id[1]+1, id[2]+1, id[3]+1, region, this);
332 [ + - ]: 308217728 : add_tetrahedron(tmp_tet);
333 : : }
334 : :
335 [ + - ]: 728 : fprintf(stderr, "\033[1;%dmDone Reading files....\033[0m\n", 33);
336 : :
337 : 728 : }
338 : :
339 : 0 : void mesh::print_vector_matlab(FILE* matlab,
340 : : vector<double> &v,
341 : : string v_name)
342 : : {
343 : 0 : fprintf(matlab, "%s = [", v_name.c_str());
344 [ # # ]: 0 : for (unsigned int i = 0; i < v.size() - 1; i++)
345 : : {
346 : 0 : fprintf(matlab, "%-15g;", v[i]);
347 [ # # ]: 0 : if ((i + 1) % 10 == 0)
348 : 0 : fprintf(matlab, "\n");
349 : : }
350 : 0 : fprintf(matlab, "%-15g", v[v.size()-1]);
351 : 0 : fprintf(matlab, "];\n\n");
352 : 0 : }
353 : :
354 : 0 : void mesh::print_vector_matlab(FILE* matlab,
355 : : vector<int> &v,
356 : : string v_name)
357 : : {
358 : 0 : fprintf(matlab, "%s = [", v_name.c_str());
359 [ # # ]: 0 : for (unsigned int i = 0; i < v.size() - 1; i++)
360 : : {
361 : 0 : fprintf(matlab, "%-15d;", v[i]);
362 [ # # ]: 0 : if ((i + 1) % 10 == 0)
363 : 0 : fprintf(matlab, "\n");
364 : : }
365 : 0 : fprintf(matlab, "%-15d", v[v.size()-1]);
366 : 0 : fprintf(matlab, "];\n\n");
367 [ + - ][ + - ]: 112 : }
|
