Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/plt_loader.h"
22 :
23 : #include "libmesh/int_range.h"
24 :
25 : // C++ includes
26 : #include <fstream>
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 :
33 : // extern "C" {
34 :
35 : // void open_ (const int *, const char *, int *);
36 : // void idata_ (const int *, int *, int *);
37 : // void fdata_ (const int *, float*, int *);
38 : // void ddata_ (const int *, double *, int *);
39 : // void close_ (const int *);
40 : // }
41 :
42 :
43 :
44 : //-----------------------------------------------------------------------------
45 : // PltLoader write members
46 : // void PltLoader::write_plot3d (const std::string & basename,
47 : // const bool reverse,
48 : // const bool gridonly) const
49 : // {
50 : // const std::string gname = basename + ".g";
51 : // const std::string qname = basename + ".q";
52 :
53 : // // FORTRAN file unit numbers
54 : // const int gunit = 25;
55 : // const int qunit = 26;
56 :
57 : // // Tell the user what files we are creating
58 : // if (this->verbose())
59 : // {
60 : // libMesh::out << "Plot3D output will be written to " << gname;
61 :
62 : // if (!gridonly)
63 : // libMesh::out << " and " << qname;
64 :
65 : // libMesh::out << std::endl;
66 : // }
67 :
68 :
69 : // // Open the FORTRAN unformatted file
70 : // {
71 : // libmesh_assert_equal_to (gname.size(), qname.size());
72 : // int len = gname.size();
73 :
74 : // open_ (&gunit, gname.c_str(), &len);
75 :
76 : // if (!gridonly)
77 : // open_ (&qunit, qname.c_str(), &len);
78 : // }
79 :
80 : // // Write the headers
81 : // {
82 : // std::vector<int> ints;
83 : // ints.reserve (3*this->n_zones());
84 :
85 : // for (unsigned int zn=0; zn<this->n_zones(); zn++)
86 : // {
87 : // ints.push_back(this->imax(zn));
88 : // ints.push_back(this->jmax(zn));
89 : // ints.push_back(this->kmax(zn));
90 : // }
91 :
92 : // int nb = this->n_zones();
93 : // int one = 1;
94 : // int len = ints.size();
95 :
96 : // libmesh_assert_equal_to (static_cast<unsigned int>(len), 3*this->n_zones());
97 :
98 : // idata_ (&gunit, &nb, &one);
99 : // idata_ (&gunit, ints.data(), &len);
100 :
101 : // if (!gridonly)
102 : // {
103 : // idata_ (&qunit, &nb, &one);
104 : // idata_ (&qunit, ints.data(), &len);
105 : // }
106 : // }
107 :
108 :
109 : // // Variables to write to plot3D file
110 : // std::vector<unsigned int> write_vars;
111 : // write_vars.reserve (5);
112 :
113 : // std::fill (write_vars.begin(), write_vars.end(), 0);
114 :
115 :
116 : // //------------------------------------------------------------------------
117 : // // Ask the user which variables to write
118 : // if (!gridonly)
119 : // {
120 : // libMesh::out << "Variables:" << std::endl;
121 :
122 : // for (unsigned int v=0; v<this->n_vars(); v++)
123 : // libMesh::out << " " << v << ") \"" << this->var_name(v) << "\""
124 : // << std::endl;
125 : // libMesh::out << std::endl;
126 :
127 : // int n_write_vars = 0;
128 :
129 : // while (true)
130 : // {
131 : // libMesh::out << "How many variables to write to the Plot3D file? 1<=n<=" << this->n_vars()
132 : // << " "
133 : // << std::endl
134 : // << "(-1 writes them all): ";
135 :
136 : // std::cin >> n_write_vars;
137 :
138 : // if (n_write_vars == -1)
139 : // break;
140 :
141 : // if ((n_write_vars >= 1) &&
142 : // (n_write_vars <= static_cast<int>(this->n_vars())))
143 : // break;
144 : // };
145 :
146 :
147 : // // The user wants all the variables
148 : // if ((n_write_vars == -1) ||
149 : // (n_write_vars == static_cast<int>(this->n_vars())))
150 : // {
151 : // for (unsigned int wv=0; wv<this->n_vars(); wv++)
152 : // write_vars.push_back (wv);
153 : // }
154 :
155 : // // The user wants a subset of the variables
156 : // else
157 : // {
158 : // libmesh_assert_greater_equal (n_write_vars, 1);
159 : // libmesh_assert_less (n_write_vars, static_cast<int>(this->n_vars()));
160 :
161 : // libMesh::out << "Select the " << n_write_vars << " variables to write to the Plot3D file: "
162 : // << std::endl;
163 :
164 : // for (int wv=0; wv<n_write_vars; wv++)
165 : // {
166 : // int num=0;
167 :
168 : // std::cin >> num;
169 :
170 : // libmesh_assert_less (num, static_cast<int>(this->n_vars()));
171 :
172 : // write_vars.push_back (num);
173 : // }
174 : // }
175 :
176 : // libMesh::out << std::endl;
177 : // } // if (!gridonly)
178 :
179 :
180 :
181 : // //------------------------------------------------------------------------
182 : // // Write the coordinates & data for each block
183 : // for (unsigned int zn=0; zn<this->n_zones(); zn++)
184 : // {
185 : // // Write the coordinates
186 : // {
187 : // std::vector<float> coords; // the nodal coordinates
188 : // coords.reserve (3*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
189 :
190 : // for (unsigned int v=0; v<3; v++)
191 : // {
192 : // unsigned int l=0;
193 :
194 : // for (unsigned int k=0; k<this->kmax(zn); k++)
195 : // for (unsigned int j=0; j<this->jmax(zn); j++)
196 : // for (unsigned int i=0; i<this->imax(zn); i++)
197 : // {
198 : // libmesh_assert_less (l, _data[zn][v].size());
199 : // coords.push_back (_data[zn][v][l++]);
200 : // }
201 : // }
202 :
203 : // // Write to the grid file
204 : // {
205 : // int len = coords.size();
206 : // libmesh_assert_equal_to (static_cast<unsigned int>(len),
207 : // 3*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
208 :
209 : // fdata_ (&gunit, coords.data(), &len);
210 : // }
211 : // }
212 :
213 :
214 : // //------------------------------------------------------------------------
215 : // // Write the data
216 : // if (!gridonly)
217 : // {
218 : // std::vector<float> data; // arbitrary data
219 : // std::vector<float> conds(4); // plot3D conditions [FSMACH, ALPHA, RE, TIME]
220 : // data.reserve (write_vars.size()*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
221 : // std::fill (conds.begin(), conds.end(), 0.);
222 :
223 : // if (zn == 0)
224 : // libMesh::out << " Writing ";
225 :
226 : // for (std::size_t i=0; i<write_vars.size(); i++)
227 : // {
228 : // // Number of the variable to write
229 : // const unsigned int v = write_vars[i];
230 :
231 : // libmesh_assert_less (v, this->n_vars());
232 :
233 : // // Tell the user what variable we are writing, but only
234 : // // once per file.
235 : // if (zn == 0)
236 : // libMesh::out << "\"" << this->var_name(v) << "\" ";
237 :
238 : // unsigned int l=0;
239 :
240 : // for (unsigned int k=0; k<this->kmax(zn); k++)
241 : // for (unsigned int j=0; j<this->jmax(zn); j++)
242 : // for (unsigned int i=0; i<this->imax(zn); i++)
243 : // {
244 : // libmesh_assert_less (l, _data[zn][v].size());
245 : // data.push_back ((v < this->n_vars()) ?
246 : // _data[zn][v][l++] : 0.);
247 : // }
248 : // }
249 :
250 : // if (zn == 0)
251 : // libMesh::out << "to " << qname << std::endl;
252 :
253 : // // Write to the solution file
254 : // {
255 : // int len = conds.size();
256 :
257 : // fdata_ (&qunit, conds.data(), &len);
258 : // }
259 :
260 : // // Write to the solution file
261 : // {
262 : // int len = data.size();
263 : // libmesh_assert_equal_to (static_cast<unsigned int>(len),
264 : // write_vars.size()*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
265 :
266 : // fdata_ (&qunit, data.data(), &len);
267 : // }
268 : // }
269 : // }
270 :
271 : // // Close the FORTRAN files
272 : // close_ (&gunit);
273 :
274 : // if (!gridonly)
275 : // close_ (&qunit);
276 :
277 : // // Possibly reverse the orders
278 : // if (reverse)
279 : // {
280 : // if (this->verbose())
281 : // libMesh::out << "Reversing byte-ordering for output files."
282 : // << std::endl;
283 :
284 : // Utility::reverse_endian (gname);
285 :
286 : // if (!gridonly)
287 : // Utility::reverse_endian (qname);
288 : // }
289 : // }
290 :
291 :
292 :
293 : // void PltLoader::write_tri (const std::string & name,
294 : // const bool reverse,
295 : // const bool gridonly) const
296 : // {
297 : // // Check out
298 : // // http://people.nas.nasa.gov/~aftosmis/cart3d/cart3dTriangulations.html
299 : // // for the .tri, .triq format
300 :
301 : // // FORTRAN file unit numbers
302 : // const int gunit = 25;
303 :
304 : // if (this->verbose())
305 : // {
306 : // libMesh::out << "Writing unformatted .tri file " << name
307 : // << std::endl;
308 :
309 : // if (gridonly)
310 : // libMesh::out << "Only writing the grid to " << name
311 : // << std::endl;
312 : // }
313 :
314 :
315 : // // Open the FORTRAN unformatted file
316 : // {
317 : // int len = name.size();
318 :
319 : // open_ (&gunit, name.c_str(), &len);
320 : // }
321 :
322 :
323 : // // Write the header
324 : // unsigned int n_nodes =0;
325 : // unsigned int n_tri =0;
326 : // unsigned int n_scalar=this->n_vars()-3;
327 :
328 : // {
329 : // std::vector<int> ints;
330 :
331 : // for (unsigned int zone=0; zone<this->n_zones(); zone++)
332 : // {
333 : // libmesh_assert_equal_to (this->elem_type(zone), TRI);
334 : // n_nodes += this->n_nodes(zone);
335 : // n_tri += this->n_elem(zone);
336 : // }
337 :
338 : // ints.push_back (n_nodes);
339 : // ints.push_back (n_tri);
340 :
341 : // if (!gridonly)
342 : // if (this->n_vars() > 3)
343 : // ints.push_back(n_scalar);
344 :
345 : // int len = ints.size();
346 : // idata_ (&gunit, ints.data(), &len);
347 : // }
348 :
349 : // // Write the nodal values.
350 : // {
351 : // std::vector<float> coords;
352 : // coords.reserve (3*n_nodes);
353 :
354 : // for (unsigned int zone=0; zone<this->n_zones(); zone++)
355 : // for (unsigned int n=0; n<this->n_nodes(zone); n++)
356 : // {
357 : // coords.push_back (_data[zone][0][n]);
358 : // coords.push_back (_data[zone][1][n]);
359 : // coords.push_back (_data[zone][2][n]);
360 : // }
361 : // // Got all the nodes for all the zones
362 :
363 :
364 : // int len = coords.size();
365 : // fdata_ (&gunit, coords.data(), &len);
366 : // }
367 :
368 : // // Write the connectivity
369 : // {
370 : // std::vector<int> conn;
371 : // conn.reserve (3*n_tri);
372 :
373 : // for (unsigned int zone=0; zone<this->n_zones(); zone++)
374 : // {
375 : // // The connectivity for this zone
376 : // const std::vector<int> & zconn = _conn[zone];
377 :
378 : // libmesh_assert (!zconn.empty());
379 :
380 : // // Append the connectivity for this zone to the connectivity
381 : // // array
382 : // conn.insert (conn.end(), zconn.begin(), zconn.end());
383 : // }
384 :
385 : // int len = conn.size();
386 : // libmesh_assert_equal_to (static_cast<unsigned int>(len), 3*n_tri);
387 : // idata_ (&gunit, conn.data(), &len);
388 : // }
389 :
390 :
391 : // // Write the component index for each triangle
392 : // {
393 : // std::vector<int> comp;
394 : // comp.reserve (n_tri);
395 :
396 : // for (unsigned int zone=0; zone<this->n_zones(); zone++)
397 : // comp.insert (comp.end(), this->n_elem(zone), zone+1);
398 :
399 : // int len = comp.size();
400 : // libmesh_assert_equal_to (static_cast<unsigned int>(len), n_tri);
401 : // idata_ (&gunit, comp.data(), &len);
402 : // }
403 :
404 :
405 : // // Possibly write additional values for each node
406 : // if (!gridonly)
407 : // if (this->n_vars() > 3)
408 : // {
409 : // if (this->verbose())
410 : // {
411 : // libMesh::out << "Writing variables ";
412 :
413 : // for (unsigned int v=3; v<this->n_vars(); v++)
414 : // libMesh::out << "\"" << this->var_name(v) << "\" ";
415 :
416 : // libMesh::out << "to the output file " << name
417 : // << std::endl;
418 : // }
419 :
420 : // std::vector<float> data;
421 :
422 : // data.reserve (n_nodes*(this->n_vars()-3));
423 :
424 : // for (unsigned int zone=0; zone<this->n_zones(); zone++)
425 : // for (unsigned int n=0; n<this->n_nodes(zone); n++)
426 : // for (unsigned int v=3; v<this->n_vars(); v++)
427 : // data.push_back (_data[zone][v][n]);
428 :
429 : // int len = data.size();
430 : // libmesh_assert_equal_to (static_cast<unsigned int>(len),
431 : // n_nodes*(this->n_vars()-3));
432 : // fdata_ (&gunit, data.data(), &len);
433 : // }
434 :
435 :
436 : // // Close the FORTRAN file
437 : // close_ (&gunit);
438 :
439 :
440 : // // Possibly reverse the orders
441 : // if (reverse)
442 : // {
443 : // if (this->verbose())
444 : // libMesh::out << "Reversing byte-ordering for output files."
445 : // << std::endl;
446 :
447 : // Utility::reverse_endian (name);
448 : // }
449 : // }
450 :
451 :
452 :
453 0 : void PltLoader::write_dat (const std::string & name,
454 : const unsigned int version_in) const
455 : {
456 0 : std::ofstream out_stream (name.c_str());
457 :
458 : out_stream << "TITLE=\""
459 0 : << this->title()
460 : << "\""
461 0 : << '\n';
462 :
463 0 : out_stream << "VARIABLES = ";
464 :
465 0 : for (auto v : make_range(this->n_vars()))
466 0 : out_stream << "\"" << this->var_name(v) << "\"\n";
467 :
468 0 : for (auto z : make_range(this->n_zones()))
469 : {
470 0 : out_stream << "ZONE T=\"" << this->zone_name(z) << "\"\n";
471 0 : out_stream << " I=" << this->imax(z)
472 0 : << ", J=" << this->jmax(z)
473 0 : << ", K=" << this->kmax(z);
474 :
475 : // Write BLOCK data for this zone
476 0 : if (this->zone_type(z) == BLOCK)
477 : {
478 0 : if (version_in < 10)
479 : {
480 0 : out_stream << ", F=BLOCK\n";
481 : }
482 : else
483 : {
484 : out_stream << ", ZONETYPE=Ordered\n"
485 0 : << "DATAPACKING=BLOCK\n";
486 : }
487 :
488 0 : out_stream << "DT=(";
489 0 : for (unsigned int v=0; v<this->n_vars(); v++)
490 0 : out_stream << "SINGLE ";
491 0 : out_stream << ")\n";
492 :
493 0 : out_stream.precision(9);
494 :
495 0 : for (auto v : make_range(this->n_vars()))
496 : {
497 0 : unsigned int l=0;
498 :
499 0 : for (unsigned int k=0; k<this->kmax(z); k++)
500 0 : for (unsigned int j=0; j<this->jmax(z); j++)
501 0 : for (unsigned int i=0; i<this->imax(z); i++)
502 : {
503 0 : out_stream << std::scientific
504 0 : << _data[z][v][l++] << " ";
505 :
506 : // Throw in a newline every 5 entries to
507 : // avoid really long lines.
508 0 : if (l%5 == 0)
509 0 : out_stream << '\n';
510 : }
511 :
512 0 : if (l%5 != 0)
513 0 : out_stream << '\n';
514 : }
515 : } // end if (this->zone_type(z) == BLOCK)
516 :
517 : // Write POINT data for this zone
518 0 : else if (this->zone_type(z) == POINT)
519 : {
520 0 : if (version_in < 10)
521 : {
522 0 : out_stream << ", F=POINT\n";
523 : }
524 : else
525 : {
526 : out_stream << ", ZONETYPE=Ordered\n"
527 0 : << "DATAPACKING=POINT\n";
528 : }
529 :
530 0 : out_stream << "DT=(";
531 0 : for (unsigned int v=0; v<this->n_vars(); v++)
532 0 : out_stream << "SINGLE ";
533 0 : out_stream << ")\n";
534 :
535 0 : out_stream.precision(9);
536 :
537 : {
538 0 : unsigned int l=0;
539 :
540 0 : for (unsigned int k=0; k<this->kmax(z); k++)
541 0 : for (unsigned int j=0; j<this->jmax(z); j++)
542 0 : for (unsigned int i=0; i<this->imax(z); i++)
543 : {
544 0 : for (auto v : make_range(this->n_vars()))
545 0 : out_stream << std::scientific
546 0 : << _data[z][v][l] << " ";
547 :
548 0 : out_stream << '\n';
549 :
550 0 : l++;
551 : }
552 : }
553 : } // end else if (this->zone_type(z) == POINT)
554 :
555 : // Otherwise, unrecognized zone type
556 : else
557 0 : libmesh_error_msg("Unrecognized zone type: this->zone_type(z)==" << this->zone_type(z));
558 : }
559 0 : }
560 :
561 : } // namespace libMesh
|