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 : // libMesh includes
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/ensight_io.h"
23 : #include "libmesh/equation_systems.h"
24 : #include "libmesh/fe_interface.h"
25 : #include "libmesh/libmesh.h"
26 : #include "libmesh/system.h"
27 : #include "libmesh/elem.h"
28 : #include "libmesh/enum_elem_type.h"
29 : #include "libmesh/int_range.h"
30 : #include "libmesh/utility.h" // map_find
31 :
32 : // C++ includes
33 : #include <sstream>
34 : #include <fstream>
35 : #include <string>
36 : #include <iomanip>
37 :
38 : namespace libMesh
39 : {
40 :
41 : // Initialize the static data members by calling the static build functions.
42 : std::map<ElemType, std::string> EnsightIO::_element_map = EnsightIO::build_element_map();
43 :
44 : // Static function used to build the _element_map.
45 16389 : std::map<ElemType, std::string> EnsightIO::build_element_map()
46 : {
47 462 : std::map<ElemType, std::string> ret;
48 16389 : ret[EDGE2] = "bar2";
49 16389 : ret[EDGE3] = "bar3";
50 16389 : ret[QUAD4] = "quad4";
51 16389 : ret[QUAD8] = "quad8";
52 : // ret[QUAD9] = "quad9"; // not supported
53 16389 : ret[TRI3] = "tria3";
54 16389 : ret[TRI6] = "tria6";
55 16389 : ret[TET4] = "tetra4";
56 16389 : ret[TET10] = "tetra10";
57 16389 : ret[HEX8] = "hexa8";
58 16389 : ret[HEX20] = "hexa20";
59 : // ret[HEX27] = "HEX27"; // not supported
60 16389 : ret[PYRAMID5] = "pyramid5";
61 16389 : return ret;
62 : }
63 :
64 :
65 0 : EnsightIO::EnsightIO (const std::string & filename,
66 0 : const EquationSystems & eq) :
67 : MeshOutput<MeshBase> (eq.get_mesh()),
68 0 : _equation_systems(eq)
69 : {
70 0 : if (_equation_systems.n_processors() == 1)
71 0 : _ensight_file_name = filename;
72 : else
73 : {
74 0 : std::ostringstream tmp_file;
75 0 : tmp_file << filename << "_rank" << _equation_systems.processor_id();
76 0 : _ensight_file_name = tmp_file.str();
77 0 : }
78 0 : }
79 :
80 :
81 :
82 0 : void EnsightIO::add_vector (const std::string & system_name,
83 : std::string_view vec_description,
84 : std::string u,
85 : std::string v)
86 : {
87 0 : libmesh_assert (_equation_systems.has_system(system_name));
88 0 : libmesh_assert (_equation_systems.get_system(system_name).has_variable(u));
89 0 : libmesh_assert (_equation_systems.get_system(system_name).has_variable(v));
90 :
91 0 : Vectors vec;
92 0 : vec.description = vec_description;
93 0 : vec.components.push_back(std::move(u));
94 0 : vec.components.push_back(std::move(v));
95 :
96 0 : _system_vars_map[system_name].EnsightVectors.push_back(std::move(vec));
97 0 : }
98 :
99 :
100 :
101 0 : void EnsightIO::add_vector (const std::string & system_name,
102 : std::string_view vec_name,
103 : std::string u,
104 : std::string v,
105 : std::string w)
106 : {
107 0 : libmesh_assert(_equation_systems.has_system(system_name));
108 0 : libmesh_assert(_equation_systems.get_system(system_name).has_variable(u));
109 0 : libmesh_assert(_equation_systems.get_system(system_name).has_variable(v));
110 0 : libmesh_assert(_equation_systems.get_system(system_name).has_variable(w));
111 :
112 0 : Vectors vec;
113 0 : vec.description = vec_name;
114 0 : vec.components.push_back(std::move(u));
115 0 : vec.components.push_back(std::move(v));
116 0 : vec.components.push_back(std::move(w));
117 0 : _system_vars_map[system_name].EnsightVectors.push_back(std::move(vec));
118 0 : }
119 :
120 :
121 :
122 0 : void EnsightIO::add_scalar(const std::string & system_name,
123 : std::string_view scl_description,
124 : std::string_view s)
125 : {
126 0 : libmesh_assert(_equation_systems.has_system(system_name));
127 0 : libmesh_assert(_equation_systems.get_system(system_name).has_variable(s));
128 :
129 0 : Scalars scl;
130 0 : scl.description = scl_description;
131 0 : scl.scalar_name = s;
132 :
133 0 : _system_vars_map[system_name].EnsightScalars.push_back(std::move(scl));
134 0 : }
135 :
136 :
137 :
138 : // This method must be implemented as it is pure virtual in
139 : // the MeshOutput base class.
140 0 : void EnsightIO::write (const std::string & name)
141 : {
142 : // We may need to gather a DistributedMesh to output it, making that
143 : // const qualifier in our constructor a dirty lie
144 0 : MeshSerializer serialize(const_cast<MeshBase &>(this->mesh()), !_is_parallel_format);
145 :
146 0 : _ensight_file_name = name;
147 0 : this->write();
148 0 : }
149 :
150 :
151 :
152 0 : void EnsightIO::write (Real time)
153 : {
154 0 : this->write_ascii(time);
155 0 : this->write_case();
156 0 : }
157 :
158 :
159 :
160 0 : void EnsightIO::write_ascii (Real time)
161 : {
162 0 : _time_steps.push_back(time);
163 :
164 0 : this->write_geometry_ascii();
165 0 : this->write_solution_ascii();
166 0 : }
167 :
168 :
169 :
170 0 : void EnsightIO::write_geometry_ascii()
171 : {
172 0 : std::ostringstream file;
173 0 : file << _ensight_file_name
174 : << ".geo"
175 : << std::setw(3)
176 : << std::setprecision(0)
177 0 : << std::setfill('0')
178 0 : << std::right
179 0 : << _time_steps.size()-1;
180 :
181 : // Open a stream to write the mesh
182 0 : std::ofstream mesh_stream(file.str().c_str());
183 :
184 0 : mesh_stream << "EnSight Gold Geometry File Format\n";
185 0 : mesh_stream << "Generated by \n";
186 0 : mesh_stream << "node id off\n";
187 0 : mesh_stream << "element id given\n";
188 0 : mesh_stream << "part\n";
189 0 : mesh_stream << std::setw(10) << 1 << "\n";
190 0 : mesh_stream << "uns-elements\n";
191 0 : mesh_stream << "coordinates\n";
192 :
193 : // mapping between nodal index and your coordinates
194 0 : std::map<int, Point> mesh_nodes_map;
195 :
196 : // Map for grouping elements of the same type
197 0 : std::map<ElemType, std::vector<const Elem *>> ensight_parts_map;
198 :
199 0 : const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
200 :
201 : // Construct the various required maps
202 0 : for (const auto & elem : the_mesh.active_local_element_ptr_range())
203 : {
204 0 : ensight_parts_map[elem->type()].push_back(elem);
205 :
206 0 : for (const Node & node : elem->node_ref_range())
207 0 : mesh_nodes_map[node.id()] = node;
208 0 : }
209 :
210 : // Write number of local points
211 0 : mesh_stream << std::setw(10) << mesh_nodes_map.size() << "\n";
212 :
213 : // write x, y, and z node positions, build mapping between
214 : // ensight and libmesh node numbers.
215 0 : std::map <int, int> ensight_node_index;
216 0 : for (unsigned direction=0; direction<3; ++direction)
217 : {
218 0 : int i = 1;
219 0 : for (const auto & [idx, pt] : mesh_nodes_map)
220 : {
221 : mesh_stream << std::setw(12)
222 0 : << std::setprecision(5)
223 0 : << std::scientific
224 0 : << pt(direction)
225 0 : << "\n";
226 0 : ensight_node_index[idx] = i++;
227 : }
228 : }
229 :
230 : // Write parts
231 0 : for (const auto & [elem_type, elem_ref] : ensight_parts_map)
232 : {
233 : // Look up this ElemType in the map, error if not present.
234 0 : std::string name = libmesh_map_find(_element_map, elem_type);
235 :
236 : // Write element type
237 0 : mesh_stream << "\n" << name << "\n";
238 :
239 : // Write number of element
240 0 : mesh_stream << std::setw(10) << elem_ref.size() << "\n";
241 :
242 : // Write element id
243 0 : for (const auto & elem : elem_ref)
244 0 : mesh_stream << std::setw(10) << elem->id() << "\n";
245 :
246 : // Write connectivity
247 0 : for (auto i : index_range(elem_ref))
248 : {
249 0 : for (const auto & node : elem_ref[i]->node_ref_range())
250 : {
251 : // tests!
252 0 : if (elem_type == QUAD9 && i==4)
253 0 : continue;
254 :
255 : // tests!
256 0 : if (elem_type == HEX27 &&
257 0 : (i==4 || i ==10 || i == 12 ||
258 0 : i == 13 || i ==14 || i == 16 || i == 22))
259 0 : continue;
260 :
261 0 : mesh_stream << std::setw(10) << ensight_node_index[node.id()];
262 : }
263 0 : mesh_stream << "\n";
264 : }
265 : }
266 0 : }
267 :
268 :
269 :
270 :
271 :
272 0 : void EnsightIO::write_case()
273 : {
274 0 : std::ostringstream case_file;
275 0 : case_file << _ensight_file_name << ".case";
276 :
277 : // Open a stream for writing the case file.
278 0 : std::ofstream case_stream(case_file.str().c_str());
279 :
280 0 : case_stream << "FORMAT\n";
281 0 : case_stream << "type: ensight gold\n\n";
282 0 : case_stream << "GEOMETRY\n";
283 0 : case_stream << "model: 1 " << _ensight_file_name << ".geo" << "***\n";
284 :
285 : // Write Variable per node section
286 0 : if (!_system_vars_map.empty())
287 0 : case_stream << "\n\nVARIABLE\n";
288 :
289 0 : for (const auto & pr : _system_vars_map)
290 : {
291 0 : for (const auto & scalar : pr.second.EnsightScalars)
292 : case_stream << "scalar per node: 1 "
293 0 : << scalar.description << " "
294 0 : << _ensight_file_name << "_" << scalar.scalar_name << ".scl***\n";
295 :
296 0 : for (const auto & vec : pr.second.EnsightVectors)
297 : case_stream << "vector per node: 1 "
298 0 : << vec.description << " "
299 0 : << _ensight_file_name << "_" << vec.description << ".vec***\n";
300 :
301 : // Write time step section
302 0 : if (_time_steps.size() != 0)
303 : {
304 0 : case_stream << "\n\nTIME\n";
305 0 : case_stream << "time set: 1\n";
306 0 : case_stream << "number of steps: " << std::setw(10) << _time_steps.size() << "\n";
307 0 : case_stream << "filename start number: " << std::setw(10) << 0 << "\n";
308 0 : case_stream << "filename increment: " << std::setw(10) << 1 << "\n";
309 0 : case_stream << "time values:\n";
310 0 : for (const auto & time : _time_steps)
311 0 : case_stream << std::setw(12) << std::setprecision(5) << std::scientific << time << "\n";
312 : }
313 : }
314 0 : }
315 :
316 :
317 : // Write scalar and vector solution
318 0 : void EnsightIO::write_solution_ascii()
319 : {
320 0 : for (const auto & [sys_name, sys_vars] : _system_vars_map)
321 : {
322 0 : for (const auto & scalar : sys_vars.EnsightScalars)
323 0 : this->write_scalar_ascii(sys_name,
324 : scalar.scalar_name);
325 :
326 0 : for (const auto & vec : sys_vars.EnsightVectors)
327 0 : this->write_vector_ascii(sys_name,
328 0 : vec.components,
329 : vec.description);
330 : }
331 0 : }
332 :
333 :
334 0 : void EnsightIO::write_scalar_ascii(std::string_view sys,
335 : std::string_view var_name)
336 : {
337 : // Construct scalar variable filename
338 0 : std::ostringstream scl_file;
339 0 : scl_file << _ensight_file_name
340 : << "_"
341 : << var_name
342 : << ".scl"
343 : << std::setw(3)
344 : << std::setprecision(0)
345 0 : << std::setfill('0')
346 0 : << std::right
347 0 : << _time_steps.size()-1;
348 :
349 : // Open a stream and start writing scalar variable info.
350 0 : std::ofstream scl_stream(scl_file.str().c_str());
351 0 : scl_stream << "Per node scalar value\n";
352 0 : scl_stream << "part\n";
353 0 : scl_stream << std::setw(10) << 1 << "\n";
354 0 : scl_stream << "coordinates\n";
355 :
356 0 : const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
357 0 : const unsigned int dim = the_mesh.mesh_dimension();
358 0 : const System & system = _equation_systems.get_system(sys);
359 0 : const DofMap & dof_map = system.get_dof_map();
360 0 : int var = system.variable_number(var_name);
361 :
362 0 : std::vector<dof_id_type> dof_indices_scl;
363 :
364 : // Map from node id -> solution value. We end up just writing this
365 : // map out in order, not sure what would happen if there were holes
366 : // in the numbering...
367 0 : std::map<int, Real> local_soln;
368 :
369 0 : std::vector<Number> elem_soln;
370 0 : std::vector<Number> nodal_soln;
371 :
372 : // Loop over active local elements, construct the nodal solution, and write it to file.
373 0 : for (const auto & elem : the_mesh.active_local_element_ptr_range())
374 : {
375 0 : const FEType & fe_type = system.variable_type(var);
376 :
377 0 : dof_map.dof_indices (elem, dof_indices_scl, var);
378 :
379 0 : elem_soln.resize(dof_indices_scl.size());
380 :
381 0 : for (auto i : index_range(dof_indices_scl))
382 0 : elem_soln[i] = system.current_solution(dof_indices_scl[i]);
383 :
384 0 : FEInterface::nodal_soln (dim, fe_type, elem, elem_soln, nodal_soln);
385 :
386 0 : libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
387 :
388 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
389 0 : libmesh_error_msg("Complex-valued Ensight output not yet supported");
390 : #endif
391 :
392 0 : for (auto n : elem->node_index_range())
393 0 : local_soln[elem->node_id(n)] = libmesh_real(nodal_soln[n]);
394 0 : }
395 :
396 0 : for (const auto & pr : local_soln)
397 : scl_stream << std::setw(12)
398 0 : << std::setprecision(5)
399 0 : << std::scientific
400 0 : << pr.second
401 0 : << "\n";
402 0 : }
403 :
404 :
405 0 : void EnsightIO::write_vector_ascii(std::string_view sys,
406 : const std::vector<std::string> & vec,
407 : std::string_view var_name)
408 : {
409 : // Construct vector variable filename
410 0 : std::ostringstream vec_file;
411 0 : vec_file << _ensight_file_name
412 : << "_"
413 : << var_name
414 : << ".vec"
415 : << std::setw(3)
416 : << std::setprecision(0)
417 0 : << std::setfill('0')
418 0 : << std::right
419 0 : << _time_steps.size()-1;
420 :
421 : // Open a stream and start writing vector variable info.
422 0 : std::ofstream vec_stream(vec_file.str().c_str());
423 0 : vec_stream << "Per vector per value\n";
424 0 : vec_stream << "part\n";
425 0 : vec_stream << std::setw(10) << 1 << "\n";
426 0 : vec_stream << "coordinates\n";
427 :
428 : // Get a constant reference to the mesh object.
429 0 : const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
430 :
431 : // The dimension that we are running
432 0 : const unsigned int dim = the_mesh.mesh_dimension();
433 :
434 0 : const System & system = _equation_systems.get_system(sys);
435 :
436 0 : const DofMap & dof_map = system.get_dof_map();
437 :
438 0 : const unsigned int u_var = system.variable_number(vec[0]);
439 0 : const unsigned int v_var = system.variable_number(vec[1]);
440 0 : const unsigned int w_var = (dim==3) ? system.variable_number(vec[2]) : 0;
441 :
442 0 : std::vector<dof_id_type> dof_indices_u;
443 0 : std::vector<dof_id_type> dof_indices_v;
444 0 : std::vector<dof_id_type> dof_indices_w;
445 :
446 : // Map from node id -> solution value. We end up just writing this
447 : // map out in order, not sure what would happen if there were holes
448 : // in the numbering...
449 0 : std::map<int,std::vector<Real>> local_soln;
450 :
451 : // Now we will loop over all the elements in the mesh.
452 0 : for (const auto & elem : the_mesh.active_local_element_ptr_range())
453 : {
454 0 : const FEType & fe_type = system.variable_type(u_var);
455 :
456 0 : dof_map.dof_indices (elem, dof_indices_u, u_var);
457 0 : dof_map.dof_indices (elem, dof_indices_v, v_var);
458 0 : if (dim==3)
459 0 : dof_map.dof_indices (elem, dof_indices_w, w_var);
460 :
461 0 : std::vector<Number> elem_soln_u;
462 0 : std::vector<Number> elem_soln_v;
463 0 : std::vector<Number> elem_soln_w;
464 :
465 0 : std::vector<Number> nodal_soln_u;
466 0 : std::vector<Number> nodal_soln_v;
467 0 : std::vector<Number> nodal_soln_w;
468 :
469 0 : elem_soln_u.resize(dof_indices_u.size());
470 0 : elem_soln_v.resize(dof_indices_v.size());
471 0 : if (dim == 3)
472 0 : elem_soln_w.resize(dof_indices_w.size());
473 :
474 0 : for (auto i : index_range(dof_indices_u))
475 : {
476 0 : elem_soln_u[i] = system.current_solution(dof_indices_u[i]);
477 0 : elem_soln_v[i] = system.current_solution(dof_indices_v[i]);
478 0 : if (dim==3)
479 0 : elem_soln_w[i] = system.current_solution(dof_indices_w[i]);
480 : }
481 :
482 0 : FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_u, nodal_soln_u);
483 0 : FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_v, nodal_soln_v);
484 0 : if (dim == 3)
485 0 : FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_w, nodal_soln_w);
486 :
487 0 : libmesh_assert_equal_to (nodal_soln_u.size(), elem->n_nodes());
488 0 : libmesh_assert_equal_to (nodal_soln_v.size(), elem->n_nodes());
489 :
490 : #ifdef LIBMESH_ENABLE_COMPLEX
491 : libmesh_error_msg("Complex-valued Ensight output not yet supported");
492 : #endif
493 :
494 0 : for (const auto & n : elem->node_index_range())
495 : {
496 0 : std::vector<Real> node_vec(3);
497 0 : node_vec[0] = libmesh_real(nodal_soln_u[n]);
498 0 : node_vec[1] = libmesh_real(nodal_soln_v[n]);
499 0 : node_vec[2] = 0.0;
500 0 : if (dim==3)
501 0 : node_vec[2] = libmesh_real(nodal_soln_w[n]);
502 0 : local_soln[elem->node_id(n)] = node_vec;
503 : }
504 0 : }
505 :
506 0 : for (unsigned dir=0; dir<3; ++dir)
507 : {
508 0 : for (const auto & pr : local_soln)
509 0 : vec_stream << std::setw(12)
510 0 : << std::scientific
511 0 : << std::setprecision(5)
512 0 : << pr.second[dir]
513 0 : << "\n";
514 : }
515 0 : }
516 :
517 : } // namespace libMesh
|
