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 : // Local includes
20 : #include "libmesh/libmesh_config.h"
21 : #include "libmesh/ucd_io.h"
22 : #include "libmesh/mesh_base.h"
23 : #include "libmesh/face_quad4.h"
24 : #include "libmesh/face_tri3.h"
25 : #include "libmesh/cell_tet4.h"
26 : #include "libmesh/cell_hex8.h"
27 : #include "libmesh/cell_prism6.h"
28 : #include "libmesh/enum_io_package.h"
29 : #include "libmesh/enum_elem_type.h"
30 : #include "libmesh/int_range.h"
31 : #include "libmesh/utility.h"
32 :
33 : #ifdef LIBMESH_HAVE_GZSTREAM
34 : # include "gzstream.h" // For reading/writing compressed streams
35 : #endif
36 :
37 : // C++ includes
38 : #include <array>
39 : #include <fstream>
40 :
41 :
42 : namespace libMesh
43 : {
44 :
45 : // Initialize the static data members by calling the static build functions.
46 : std::map<ElemType, std::string> UCDIO::_writing_element_map = UCDIO::build_writing_element_map();
47 : std::map<std::string, ElemType> UCDIO::_reading_element_map = UCDIO::build_reading_element_map();
48 :
49 :
50 :
51 : // Static function used to build the _writing_element_map.
52 16389 : std::map<ElemType, std::string> UCDIO::build_writing_element_map()
53 : {
54 462 : std::map<ElemType, std::string> ret;
55 16389 : ret[EDGE2] = "edge";
56 16389 : ret[TRI3] = "tri";
57 16389 : ret[QUAD4] = "quad";
58 16389 : ret[TET4] = "tet";
59 16389 : ret[HEX8] = "hex";
60 16389 : ret[PRISM6] = "prism";
61 16389 : ret[PYRAMID5] = "pyramid";
62 16389 : return ret;
63 : }
64 :
65 :
66 :
67 : // Static function used to build the _reading_element_map.
68 16389 : std::map<std::string, ElemType> UCDIO::build_reading_element_map()
69 : {
70 462 : std::map<std::string, ElemType> ret;
71 16389 : ret["edge"] = EDGE2;
72 16389 : ret["tri"] = TRI3;
73 16389 : ret["quad"] = QUAD4;
74 16389 : ret["tet"] = TET4;
75 16389 : ret["hex"] = HEX8;
76 16389 : ret["prism"] = PRISM6;
77 16389 : ret["pyramid"] = PYRAMID5;
78 16389 : return ret;
79 : }
80 :
81 :
82 24 : void UCDIO::read (const std::string & file_name)
83 : {
84 24 : if (Utility::contains(file_name, ".gz"))
85 : {
86 : #ifdef LIBMESH_HAVE_GZSTREAM
87 28 : igzstream in_stream (file_name.c_str());
88 24 : this->read_implementation (in_stream);
89 : #else
90 : libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
91 : #endif
92 20 : }
93 :
94 : else
95 : {
96 0 : std::ifstream in_stream (file_name.c_str());
97 0 : this->read_implementation (in_stream);
98 0 : }
99 24 : }
100 :
101 :
102 :
103 0 : void UCDIO::write (const std::string & file_name)
104 : {
105 0 : if (Utility::contains(file_name, ".gz"))
106 : {
107 : #ifdef LIBMESH_HAVE_GZSTREAM
108 0 : ogzstream out_stream (file_name.c_str());
109 0 : this->write_implementation (out_stream);
110 : #else
111 : libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
112 : #endif
113 0 : }
114 :
115 : else
116 : {
117 0 : std::ofstream out_stream (file_name.c_str());
118 0 : this->write_implementation (out_stream);
119 0 : }
120 0 : }
121 :
122 :
123 :
124 24 : void UCDIO::read_implementation (std::istream & in)
125 : {
126 : // This is a serial-only process for now;
127 : // the Mesh should be read on processor 0 and
128 : // broadcast later
129 2 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
130 :
131 : // Check input buffer
132 2 : libmesh_assert (in.good());
133 :
134 24 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
135 :
136 : // Keep track of what kinds of elements this file contains
137 24 : elems_of_dimension.clear();
138 24 : elems_of_dimension.resize(4, false);
139 :
140 24 : this->skip_comment_lines (in, '#');
141 :
142 24 : unsigned int nNodes=0, nElem=0, dummy=0;
143 :
144 2 : in >> nNodes // Read the number of nodes from the stream
145 2 : >> nElem // Read the number of elements from the stream
146 2 : >> dummy
147 2 : >> dummy
148 2 : >> dummy;
149 :
150 :
151 : // Read the nodal coordinates. Note that UCD format always
152 : // stores (x,y,z), and in 2D z=0. We don't need to store this,
153 : // however. So, we read in x,y,z for each node and make a point
154 : // in the proper way based on what dimension we're in
155 : {
156 73692 : for (unsigned int i=0; i<nNodes; i++)
157 : {
158 6139 : libmesh_assert (in.good());
159 :
160 : std::array<Real, 3> xyz;
161 :
162 6139 : in >> dummy // Point number
163 6139 : >> xyz[0] // x-coordinate value
164 6139 : >> xyz[1] // y-coordinate value
165 6139 : >> xyz[2]; // z-coordinate value
166 :
167 73668 : Point p(xyz[0]);
168 : #if LIBMESH_DIM > 1
169 73668 : p(1) = xyz[1];
170 : #else
171 : libmesh_assert_equal_to(xyz[1], 0);
172 : #endif
173 : #if LIBMESH_DIM > 2
174 73668 : p(2) = xyz[2];
175 : #else
176 : libmesh_assert_equal_to(xyz[2], 0);
177 : #endif
178 :
179 : // Build the node
180 73668 : mesh.add_point (p, i);
181 : }
182 : }
183 :
184 : // Read the elements from the stream. Notice that the UCD node-numbering
185 : // scheme is 1-based, and we just created a 0-based scheme above
186 : // (which is of course what we want). So, when we read in the nodal
187 : // connectivity for each element we need to take 1 off the value of
188 : // each node so that we get the right thing.
189 : {
190 24 : unsigned int material_id=0, node=0;
191 4 : std::string type;
192 :
193 105912 : for (unsigned int i=0; i<nElem; i++)
194 : {
195 8824 : libmesh_assert (in.good());
196 :
197 : // The cell type can be either tri, quad, tet, hex, or prism.
198 8824 : in >> dummy // Cell number, means nothing to us
199 8824 : >> material_id // We'll use this for the element subdomain id.
200 105888 : >> type; // string describing cell type
201 :
202 : // Convert the UCD type string to a libmesh ElemType
203 105888 : ElemType et = libmesh_map_find(_reading_element_map, type);
204 :
205 : // Build the required type and release it into our custody.
206 105888 : auto elem = Elem::build(et);
207 :
208 461184 : for (auto n : elem->node_index_range())
209 : {
210 29608 : libmesh_assert (in.good());
211 :
212 29608 : in >> node; // read the current node
213 355296 : node -= 1; // UCD is 1-based, so subtract
214 :
215 29608 : libmesh_assert_less (node, mesh.n_nodes());
216 :
217 : // assign the node
218 355296 : elem->set_node(n, mesh.node_ptr(node));
219 : }
220 :
221 105888 : elems_of_dimension[elem->dim()] = true;
222 :
223 : // Set the element's subdomain ID based on the material_id.
224 114712 : elem->subdomain_id() = cast_int<subdomain_id_type>(material_id);
225 :
226 : // Add the element to the mesh
227 88240 : elem->set_id(i);
228 123536 : mesh.add_elem (std::move(elem));
229 88240 : }
230 :
231 : // Set the mesh dimension to the largest encountered for an element
232 120 : for (unsigned char i=0; i!=4; ++i)
233 104 : if (elems_of_dimension[i])
234 24 : mesh.set_mesh_dimension(i);
235 :
236 : #if LIBMESH_DIM < 3
237 : libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
238 : "Cannot open dimension "
239 : << mesh.mesh_dimension()
240 : << " mesh file when configured without "
241 : << mesh.mesh_dimension()
242 : << "D support.");
243 : #endif
244 : }
245 24 : }
246 :
247 :
248 :
249 0 : void UCDIO::write_implementation (std::ostream & out_stream)
250 : {
251 0 : libmesh_assert (out_stream.good());
252 :
253 0 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
254 :
255 : // UCD doesn't work any dimension except 3?
256 0 : libmesh_error_msg_if(mesh.mesh_dimension() != 3,
257 : "Error: Can't write boundary elements for meshes of dimension less than 3. "
258 : "Mesh dimension = " << mesh.mesh_dimension());
259 :
260 : // Write header
261 0 : this->write_header(out_stream, mesh, mesh.n_elem(), 0);
262 :
263 : // Write the node coordinates
264 0 : this->write_nodes(out_stream, mesh);
265 :
266 : // Write the elements
267 0 : this->write_interior_elems(out_stream, mesh);
268 0 : }
269 :
270 :
271 :
272 12 : void UCDIO::write_header(std::ostream & out_stream,
273 : const MeshBase & mesh,
274 : dof_id_type n_elems,
275 : unsigned int n_vars)
276 : {
277 1 : libmesh_assert (out_stream.good());
278 : // TODO: We could print out the libmesh revision used to write this file here.
279 : out_stream << "# For a description of the UCD format see the AVS Developer's guide.\n"
280 12 : << "#\n";
281 :
282 : // Write the mesh info
283 22 : out_stream << mesh.n_nodes() << " "
284 11 : << n_elems << " "
285 1 : << n_vars << " "
286 12 : << " 0 0\n";
287 12 : }
288 :
289 :
290 :
291 12 : void UCDIO::write_nodes(std::ostream & out_stream,
292 : const MeshBase & mesh)
293 : {
294 : // 1-based node number for UCD
295 1 : unsigned int n=1;
296 :
297 : // Write the node coordinates
298 4775 : for (auto & node : mesh.node_ptr_range())
299 : {
300 216 : libmesh_assert (out_stream.good());
301 :
302 2592 : out_stream << n++ << "\t";
303 2592 : node->write_unformatted(out_stream);
304 10 : }
305 12 : }
306 :
307 :
308 :
309 12 : void UCDIO::write_interior_elems(std::ostream & out_stream,
310 : const MeshBase & mesh)
311 : {
312 : // 1-based element number for UCD
313 1 : unsigned int e=1;
314 :
315 : // Write element information
316 2773 : for (const auto & elem : mesh.element_ptr_range())
317 : {
318 125 : libmesh_assert (out_stream.good());
319 :
320 : // Look up the corresponding UCD element type in the static map.
321 1625 : std::string elem_string = libmesh_map_find(_writing_element_map, elem->type());
322 :
323 : // Write the element's subdomain ID as the UCD "material_id".
324 3000 : out_stream << e++ << " " << elem->subdomain_id() << " " << elem_string << "\t";
325 1500 : elem->write_connectivity(out_stream, UCD);
326 10 : }
327 12 : }
328 :
329 :
330 :
331 71 : void UCDIO::write_nodal_data(const std::string & fname,
332 : const std::vector<Number> & soln,
333 : const std::vector<std::string> & names)
334 : {
335 4 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
336 :
337 71 : const dof_id_type n_elem = mesh.n_elem();
338 :
339 : // Only processor 0 does the writing
340 73 : if (mesh.processor_id())
341 59 : return;
342 :
343 14 : std::ofstream out_stream(fname.c_str());
344 :
345 : // UCD doesn't work in 1D
346 1 : libmesh_assert (mesh.mesh_dimension() != 1);
347 :
348 : // Write header
349 13 : this->write_header(out_stream,mesh,n_elem,
350 : cast_int<unsigned int>(names.size()));
351 :
352 : // Write the node coordinates
353 12 : this->write_nodes(out_stream, mesh);
354 :
355 : // Write the elements
356 12 : this->write_interior_elems(out_stream, mesh);
357 :
358 : // Write the solution
359 12 : this->write_soln(out_stream, mesh, names, soln);
360 10 : }
361 :
362 :
363 :
364 12 : void UCDIO::write_soln(std::ostream & out_stream,
365 : const MeshBase & mesh,
366 : const std::vector<std::string> & names,
367 : const std::vector<Number> & soln)
368 : {
369 1 : libmesh_assert (out_stream.good());
370 :
371 : // First write out how many variables and how many components per variable
372 2 : out_stream << names.size();
373 48 : for (std::size_t i = 0, ns = names.size(); i < ns; i++)
374 : {
375 3 : libmesh_assert (out_stream.good());
376 : // Each named variable has only 1 component
377 36 : out_stream << " 1";
378 : }
379 1 : out_stream << std::endl;
380 :
381 : // Now write out variable names and units. Since we don't store units
382 : // We just write out dummy.
383 48 : for (const auto & name : names)
384 : {
385 3 : libmesh_assert (out_stream.good());
386 33 : out_stream << name << ", dummy" << std::endl;
387 : }
388 :
389 : // Now, for each node, write out the solution variables.
390 : // We use a 1-based node numbering for UCD.
391 2 : std::size_t nv = names.size();
392 2593 : for (auto n : IntRange<std::size_t>(1, mesh.n_nodes()))
393 : {
394 215 : libmesh_assert (out_stream.good());
395 215 : out_stream << n;
396 :
397 10320 : for (std::size_t var = 0; var != nv; var++)
398 : {
399 7740 : std::size_t idx = nv*(n-1) + var;
400 :
401 8385 : out_stream << " " << soln[idx];
402 : }
403 215 : out_stream << std::endl;
404 : }
405 12 : }
406 :
407 : } // namespace libMesh
|
