mesh_data_unv_support.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 // C++ includes
21 #include <cstdio> // for std::sprintf
22 #include <fstream>
23 
24 // Local includes
25 #include "libmesh/libmesh_config.h"
26 #include "libmesh/mesh_data.h"
27 #include "libmesh/auto_ptr.h"
28 
29 #ifdef LIBMESH_HAVE_GZSTREAM
30 # include "gzstream.h" // For reading/writing compressed streams
31 #endif
32 
33 
34 namespace libMesh
35 {
36 
37 //------------------------------------------------------
38 // MeshData UNV support functions
39 void MeshData::read_unv (const std::string& file_name)
40 {
41  /*
42  * we should better be active or in compatibility mode
43  */
45 
46  /*
47  * When reading data, make sure the id maps are ok
48  */
51 
52  /*
53  * clear the data, but keep the id maps
54  */
55  this->clear();
56 
57  /*
58  * We can read either ".unv", or ".unv.gz"
59  * files, provided zlib.h is there
60  */
61  if (file_name.rfind(".gz") < file_name.size())
62  {
63 #ifdef LIBMESH_HAVE_GZSTREAM
64  igzstream in_stream(file_name.c_str());
65  this->read_unv_implementation (in_stream);
66 #else
67  libMesh::err << "ERROR: You must have the zlib.h header "
68  << "files and libraries to read and write "
69  << "compressed streams."
70  << std::endl;
71  libmesh_error();
72 #endif
73  return;
74  }
75 
76  else
77  {
78  std::ifstream in_stream(file_name.c_str());
79  this->read_unv_implementation (in_stream);
80  return;
81  }
82 }
83 
84 
85 
86 
87 
88 
89 void MeshData::read_unv_implementation (std::istream& in_file)
90 {
91  /*
92  * This is the actual implementation of
93  * reading in UNV format. This enables
94  * to read either through the conventional
95  * C++ stream, or through a stream that
96  * allows to read .gz'ed files.
97  */
98  if ( !in_file.good() )
99  {
100  libMesh::err << "ERROR: Input file not good."
101  << std::endl;
102  libmesh_error();
103  }
104 
105  const std::string _label_dataset_mesh_data = "2414";
106 
107  /*
108  * locate the beginning of data set
109  * and read it.
110  */
111  {
112  std::string olds, news;
113 
114  while (true)
115  {
116  in_file >> olds >> news;
117 
118  /*
119  * Yes, really dirty:
120  *
121  * When we found a dataset, and the user does
122  * not want this dataset, we jump back here
123  */
124  go_and_find_the_next_dataset:
125 
126  /*
127  * a "-1" followed by a number means the beginning of a dataset
128  * stop combing at the end of the file
129  */
130  while( ((olds != "-1") || (news == "-1") ) && !in_file.eof() )
131  {
132  olds = news;
133  in_file >> news;
134  }
135 
136  if(in_file.eof())
137  break;
138 
139  /*
140  * if beginning of dataset
141  */
142  if (news == _label_dataset_mesh_data)
143  {
144 
145  /*
146  * Now read the data of interest.
147  * Start with the header. For
148  * explanation of the variable
149  * dataset_location, see below.
150  */
151  unsigned int dataset_location;
152 
153  /*
154  * the type of data (complex, real,
155  * float, double etc, see below)
156  */
157  unsigned int data_type;
158 
159  /*
160  * the number of floating-point values per entity
161  */
162  unsigned int NVALDC;
163 
164 
165  /*
166  * If there is no MeshDataUnvHeader object
167  * attached
168  */
169  if (_unv_header==NULL)
170  {
171  /*
172  * Ignore the first lines that stand for
173  * analysis dataset label and name.
174  */
175  for(unsigned int i=0; i<3; i++)
176  in_file.ignore(256,'\n');
177 
178  /*
179  * Read the dataset location, where
180  * 1: Data at nodes
181  * 2: Data on elements
182  * other sets are currently not supported.
183  */
184  in_file >> dataset_location;
185 
186  /*
187  * Ignore five ID lines.
188  */
189  for(unsigned int i=0; i<6; i++)
190  in_file.ignore(256,'\n');
191 
192  /*
193  * These data are all of no interest to us...
194  */
195  unsigned int model_type,
196  analysis_type,
197  data_characteristic,
198  result_type;
199 
200  /*
201  * Read record 9.
202  */
203  in_file >> model_type // not used here
204  >> analysis_type // not used here
205  >> data_characteristic // not used here
206  >> result_type // not used here
207  >> data_type
208  >> NVALDC;
209 
210 
211  /*
212  * Ignore record 10 and 11
213  * (Integer analysis type specific data).
214  */
215  for (unsigned int i=0; i<3; i++)
216  in_file.ignore(256,'\n');
217 
218  /*
219  * Ignore record 12 and record 13. Since there
220  * exist UNV files with 'D' instead of 'e' as
221  * 10th-power char, it is safer to use a string
222  * to read the dummy reals.
223  */
224  {
225  std::string dummy_Real;
226  for (unsigned int i=0; i<12; i++)
227  in_file >> dummy_Real;
228  }
229 
230  }
231  else
232  {
233 
234  /*
235  * the read() method returns false when
236  * the user wanted a special header, and
237  * when the current header is _not_ the correct
238  * header
239  */
240  if (_unv_header->read(in_file))
241  {
242  dataset_location = _unv_header->dataset_location;
243  NVALDC = _unv_header->nvaldc;
244  data_type = _unv_header->data_type;
245  }
246  else
247  {
248  /*
249  * This is not the correct header. Go
250  * and find the next. For this to
251  * work correctly, shift to the
252  * next line, so that the "-1"
253  * disappears from olds
254  */
255  olds = news;
256  in_file >> news;
257 
258  /*
259  * No good style, i know...
260  */
261  goto go_and_find_the_next_dataset;
262  }
263 
264  }
265 
266  /*
267  * Check the location of the dataset.
268  */
269  if (dataset_location != 1)
270  {
271  libMesh::err << "ERROR: Currently only Data at nodes is supported."
272  << std::endl;
273  libmesh_error();
274  }
275 
276 
277  /*
278  * Now get the foreign node id number and the respective nodal data.
279  */
280  int f_n_id;
281  std::vector<Number> values;
282 
283  while(true)
284  {
285  in_file >> f_n_id;
286 
287  /*
288  * if node_nr = -1 then we have reached the end of the dataset.
289  */
290  if (f_n_id==-1)
291  break;
292 
293  /*
294  * Resize the values vector (usually data in three
295  * principle directions, i.e. NVALDC = 3).
296  */
297  values.resize(NVALDC);
298 
299  /*
300  * Read the meshdata for the respective node.
301  */
302  for (unsigned int data_cnt=0; data_cnt<NVALDC; data_cnt++)
303  {
304  /*
305  * Check what data type we are reading.
306  * 2,4: Real
307  * 5,6: Complex
308  * other data types are not supported yet.
309  * As again, these floats may also be written
310  * using a 'D' instead of an 'e'.
311  */
312  if (data_type == 2 || data_type == 4)
313  {
314  std::string buf;
315  in_file >> buf;
317 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
318  values[data_cnt] = Complex(std::atof(buf.c_str()), 0.);
319 #else
320  values[data_cnt] = std::atof(buf.c_str());
321 #endif
322  }
323 
324  else if(data_type == 5 || data_type == 6)
325 
326  {
327 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
328  Real re_val, im_val;
329 
330  std::string buf;
331  in_file >> buf;
332 
334  {
335  re_val = std::atof(buf.c_str());
336  in_file >> buf;
338  im_val = std::atof(buf.c_str());
339  }
340  else
341  {
342  re_val = std::atof(buf.c_str());
343  in_file >> im_val;
344  }
345 
346  values[data_cnt] = Complex(re_val,im_val);
347 #else
348 
349  libMesh::err << "ERROR: Complex data only supported" << std::endl
350  << "when libMesh is configured with --enable-complex!"
351  << std::endl;
352  libmesh_error();
353 #endif
354  }
355 
356  else
357  {
358  libMesh::err << "ERROR: Data type not supported."
359  << std::endl;
360  libmesh_error();
361  }
362 
363  } // end loop data_cnt
364 
365  /*
366  * Add the values vector to the MeshData data structure.
367  */
368  const Node* node = foreign_id_to_node(f_n_id);
369  _node_data.insert (std::make_pair(node, values));
370 
371  } // while(true)
372  }
373 
374 
375  else
376  {
377  /*
378  * all other datasets are ignored
379  */
380  }
381 
382  }
383  }
384 
385 
386  /*
387  * finished reading. Ready for use, provided
388  * there was any data contained in the file.
389  */
390  libmesh_assert ((this->_node_data.size() != 0) || (this->_elem_data.size() != 0));
391 
392  this->_node_data_closed = true;
393  this->_elem_data_closed = true;
394 }
395 
396 
397 
398 
399 
400 
401 void MeshData::write_unv (const std::string& file_name)
402 {
403  /*
404  * we should better be active or in compatibility mode
405  */
406  libmesh_assert (this->_active || this->_compatibility_mode);
407 
408  /*
409  * make sure the id maps are ready
410  * and that we have data to write
411  */
414 
417 
418  if (file_name.rfind(".gz") < file_name.size())
419  {
420 #ifdef LIBMESH_HAVE_GZSTREAM
421  ogzstream out_stream(file_name.c_str());
422  this->write_unv_implementation (out_stream);
423 #else
424  libMesh::err << "ERROR: You must have the zlib.h header "
425  << "files and libraries to read and write "
426  << "compressed streams."
427  << std::endl;
428  libmesh_error();
429 #endif
430  return;
431 
432  }
433 
434  else
435  {
436  std::ofstream out_stream(file_name.c_str());
437  this->write_unv_implementation (out_stream);
438  return;
439  }
440 }
441 
442 
443 
444 
445 
446 
447 void MeshData::write_unv_implementation (std::ostream& out_file)
448 {
449  /*
450  * This is the actual implementation of writing
451  * unv files, either as .unv or as .unv.gz file
452  */
453  if ( !out_file.good() )
454  {
455  libMesh::err << "ERROR: Output file not good."
456  << std::endl;
457  libmesh_error();
458  }
459 
460 
461  /*
462  * the beginning marker of the dataset block for
463  * nodal/element-associated data (not to be confused
464  * with _desired_dataset_label!)
465  */
466  const std::string _label_dataset_mesh_data = "2414";
467 
468  /*
469  * Currently this function handles only nodal data.
470  */
471  libmesh_assert (!_node_data.empty());
472 
473  if (!_elem_data.empty())
474  libMesh::err << "WARNING: MeshData currently only supports nodal data for Universal files."
475  << std::endl
476  << " Will proceed writing only nodal data, ignoring element data."
477  << std::endl;
478 
479 
480  /*
481  * Write the beginning of the dataset.
482  */
483  out_file << " -1\n"
484  << " "
485  << _label_dataset_mesh_data
486  << "\n";
487 
488  /*
489  * Write the header
490  */
491  if (_unv_header==NULL)
492  {
493  /*
494  * create a header that holds at
495  * least sufficient data to specify
496  * what this data set currently holds.
497  *
498  * The empty constructor automatically
499  * takes care of \p dataset_location
500  * and \p data_type.
501  */
502  MeshDataUnvHeader my_header;
503 
504  /*
505  * It remains to set the correct nvaldc...
506  */
507  my_header.nvaldc = this->n_val_per_node();
508 
509  /*
510  * and the correct data type. By default
511  * only distinguish complex or real data.
512  */
513 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
514  my_header.data_type = 5;
515 #else
516  my_header.data_type = 2;
517 #endif
518 
519  /*
520  * write this default header, then let
521  * the AutoPtr go out of scope. This
522  * will take care of memory management.
523  */
524  my_header.write (out_file);
525  }
526 
527  else
528  {
529  /*
530  * make sure our nvaldc coincide.
531  */
532  if (this->n_val_per_node() != _unv_header->nvaldc)
533  {
534  libMesh::err << "WARNING: nvaldc=" << _unv_header->nvaldc
535  << " of attached MeshDataUnvHeader object not valid!" << std::endl
536  << " re-set nvaldc to " << this->n_val_per_node() << std::endl;
537  _unv_header->nvaldc = this->n_val_per_node();
538  }
539 
540 
541  /*
542  * only issue a warning when data_type does
543  * not coincide. Perhaps user provided some
544  * other header in order to convert complex
545  * to real...
546  */
547 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
548  const unsigned int my_data_type = 5;
549 #else
550  const unsigned int my_data_type = 2;
551 #endif
552  if (my_data_type != _unv_header->data_type)
553  {
554  libMesh::err << "WARNING: data_type=" << _unv_header->data_type
555  << " of attached MeshDataUnvHeader differs from" << std::endl
556  << " default value=" << my_data_type
557  << " Perhaps the user wanted this," << std::endl
558  << " so I use the value from the MeshDataUnvHeader."
559  << std::endl;
560  }
561  _unv_header->write (out_file);
562  }
563 
564 
565  /*
566  * Write the foreign node number and the respective data.
567  */
568  std::map<const Node*,
569  std::vector<Number> >::const_iterator nit = _node_data.begin();
570 
571  char buf[27];
572  for (; nit != _node_data.end(); ++nit)
573  {
574  const Node* node = (*nit).first;
575 
576  unsigned int f_n_id = node_to_foreign_id (node);
577  std::sprintf(buf, "%10i\n", f_n_id);
578  out_file << buf;
579 
580  /* since we are iterating over our own map, this libmesh_assert
581  * should never break...
582  */
583  libmesh_assert (this->has_data(node));
584 
585  // const reference to the nodal values
586  const std::vector<Number>& values = this->get_data(node);
587 
588  for (unsigned int v_cnt=0; v_cnt<values.size(); v_cnt++)
589  {
590 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
591  std::sprintf(buf, "%13.5E%13.5E", values[v_cnt].real(),
592  values[v_cnt].imag());
593  out_file << buf;
594 #else
595  std::sprintf(buf, "%13.5E",
596  static_cast<double>(values[v_cnt]));
597  out_file << buf;
598 #endif
599  }
600 
601  out_file << "\n";
602 
603 
604  }
605 
606  /*
607  * Write end of the dataset.
608  */
609  out_file << " -1\n";
610 }
611 
612 
613 
614 
615 
616 //------------------------------------------------------
617 // MeshDataUnvHeader functions
619  dataset_label (0),
620  dataset_name ("libMesh mesh data"),
621  dataset_location (1), // default to nodal data
622  model_type (0),
623  analysis_type (0),
624  data_characteristic (0),
625  result_type (0),
626 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
627  data_type (5), // default to single precision complex
628 #else
629  data_type (2), // default to single precision real
630 #endif
631  nvaldc (3), // default to 3 (principle directions)
632  _desired_dataset_label (libMesh::invalid_uint)
633 {
634  id_lines_1_to_5.resize(5);
635  std::fill (id_lines_1_to_5.begin(), id_lines_1_to_5.end(), std::string("libMesh default"));
636  /*
637  * resize analysis specific data.
638  */
639  record_10.resize(8);
640  record_11.resize(2);
641  record_12.resize(6);
642  record_13.resize(6);
643 }
644 
645 
646 
647 
648 
650 {
651  // empty
652 }
653 
654 
655 
656 
657 bool MeshDataUnvHeader::read (std::istream& in_file)
658 {
659  in_file >> this->dataset_label;
660 
661  /*
662  * currently, we compare only the
663  * dataset_label with the _desired_dataset_label,
664  * but it may be easy to also compare the
665  * dataset_name.
666  *
667  * When the user provided a dataset label, and
668  * the current label does _not_ match, then just
669  * return false.
670  *
671  * Otherwise: when the current label matches,
672  * or when there is no desired dataset label,
673  * simply proceed.
674  */
676  (this->dataset_label != this->_desired_dataset_label))
677  return false;
678 
679 
680  in_file.ignore(256,'\n');
681  std::getline(in_file, dataset_name, '\n');
682 
683  in_file >> this->dataset_location;
684  in_file.ignore(256,'\n');
685 
686 
687  for (unsigned int n=0; n<5; n++)
688  std::getline(in_file, this->id_lines_1_to_5[n], '\n');
689 
690 
691  in_file >> this->model_type
692  >> this->analysis_type
693  >> this->data_characteristic
694  >> this->result_type
695  >> this->data_type
696  >> this->nvaldc;
697 
698  for (unsigned int i=0; i<8; i++)
699  in_file >> this->record_10[i];
700 
701  for (unsigned int i=0; i<2; i++)
702  in_file >> this->record_11[i];
703 
704 
705  /*
706  * There are UNV-files where floats are
707  * written with 'D' as the 10th-power
708  * character. Replace this 'D' by 'e',
709  * so that std::atof() can work fine.
710  */
711  std::string buf;
712  in_file >> buf;
713 
714  if (need_D_to_e(buf))
715  {
716  // have to convert _all_ 'D' to 'e'
717  this->record_12[0] = std::atof(buf.c_str());
718 
719  for (unsigned int i=1; i<6; i++)
720  {
721  in_file >> buf;
722  need_D_to_e(buf);
723  this->record_12[i] = std::atof(buf.c_str());
724  }
725 
726  for (unsigned int i=0; i<6; i++)
727  {
728  in_file >> buf;
729  need_D_to_e(buf);
730  this->record_13[i] = std::atof(buf.c_str());
731  }
732  }
733  else
734  {
735  // no 'D', the stream will recognize the floats
736  this->record_12[0] = std::atof(buf.c_str());
737 
738  for (unsigned int i=1; i<6; i++)
739  in_file >> this->record_12[i];
740 
741  for (unsigned int i=0; i<6; i++)
742  in_file >> this->record_13[i];
743  }
744 
745  /*
746  * no matter whether the user provided a desired
747  * dataset label or not: return true, b/c the
748  * non-match was already caught before.
749  */
750  return true;
751 }
752 
753 
754 
755 
756 void MeshDataUnvHeader::write (std::ostream& out_file)
757 {
758 
759 
760  char buf[82];
761 
762  std::sprintf(buf, "%6i\n",this->dataset_label);
763 
764  out_file << buf;
765 
766  out_file << this->dataset_name << "\n";
767 
768  std::sprintf(buf, "%6i\n",this->dataset_location);
769 
770  out_file << buf;
771 
772  for (unsigned int n=0; n<5; n++)
773  out_file << this->id_lines_1_to_5[n] << "\n";
774 
775  std::sprintf(buf, "%10i%10i%10i%10i%10i%10i\n",
778 
779  out_file << buf;
780 
781  std::sprintf(buf, "%10i%10i%10i%10i%10i%10i%10i%10i\n",
782  record_10[0], record_10[1], record_10[2], record_10[3],
783  record_10[4], record_10[5], record_10[6], record_10[7]);
784 
785  out_file << buf;
786 
787  std::sprintf(buf, "%10i%10i\n", record_11[0], record_11[1]);
788  out_file << buf;
789 
790  std::sprintf(buf, "%13.5E%13.5E%13.5E%13.5E%13.5E%13.5E\n",
791  static_cast<double>(record_12[0]),
792  static_cast<double>(record_12[1]),
793  static_cast<double>(record_12[2]),
794  static_cast<double>(record_12[3]),
795  static_cast<double>(record_12[4]),
796  static_cast<double>(record_12[5]));
797 
798  out_file << buf;
799 
800  std::sprintf(buf, "%13.5E%13.5E%13.5E%13.5E%13.5E%13.5E\n",
801  static_cast<double>(record_13[0]),
802  static_cast<double>(record_13[1]),
803  static_cast<double>(record_13[2]),
804  static_cast<double>(record_13[3]),
805  static_cast<double>(record_13[4]),
806  static_cast<double>(record_13[5]));
807 
808  out_file << buf;
809 }
810 
811 
812 
813 
814 
815 bool MeshDataUnvHeader::need_D_to_e (std::string& number)
816 {
817  // find "D" in string, start looking at 6th element, to improve speed.
818  // We dont expect a "D" earlier
819 
820 // #ifdef __HP_aCC
821 // // Use an "int" instead of unsigned int,
822 // // otherwise HP aCC may crash!
823 // const int position = number.find("D",6);
824 // #else
825 // const unsigned int position = number.find("D",6);
826 // #endif
827  std::string::size_type position = number.find("D",6);
828 
829  if(position!=std::string::npos) // npos means no position
830  {
831  // replace "D" in string
832  number.replace(position,1,"e");
833  return true;
834  }
835  else
836  // we assume that if this one number is written correctly, all numbers are
837  return false;
838 }
839 
840 
841 
842 void MeshDataUnvHeader::which_dataset (const unsigned int ds_label)
843 {
844  this->_desired_dataset_label = ds_label;
845 }
846 
847 
848 
850 {
851  this->dataset_label = omduh.dataset_label;
852  this->dataset_name = omduh.dataset_name;
853  this->dataset_location = omduh.dataset_location;
854  this->id_lines_1_to_5 = omduh.id_lines_1_to_5;
855 
856  this->model_type = omduh.model_type;
857  this->analysis_type = omduh.analysis_type;
859  this->result_type = omduh.result_type;
860 
861 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
862  /*
863  * in complex mode allow only
864  * values 5 or 6 (complex) for data_type
865  */
866  if ((omduh.data_type == 5) ||
867  (omduh.data_type == 6))
868  this->data_type = omduh.data_type;
869  else
870  {
871 # ifdef DEBUG
872  libMesh::err << "WARNING: MeshDataUnvHeader::operator=(): Other object has data_type for" << std::endl
873  << " real values. Will use default data_type=5 during assignment." << std::endl
874  << std::endl;
875 # endif
876  this->data_type = 5;
877  }
878 
879 #else
880 
881  /*
882  * in real mode allow only
883  * values 2 or 4 (real) for data_type
884  */
885  if ((omduh.data_type == 2) ||
886  (omduh.data_type == 4))
887  this->data_type = omduh.data_type;
888  else
889  {
890 # ifdef DEBUG
891  libMesh::err << "WARNING: Other MeshDataUnvHeader has data_type for complex values." << std::endl
892  << " Data import will likely _not_ work and result in infinite loop," << std::endl
893  << " provided the user forgot to re-size nvaldc to 2*nvaldc_old!" << std::endl
894  << std::endl;
895 # endif
896  this->data_type = 2;
897  }
898 
899 #endif
900 
901  this->nvaldc = omduh.nvaldc;
902 
903  this->record_10 = omduh.record_10;
904  this->record_11 = omduh.record_11;
905  this->record_12 = omduh.record_12;
906  this->record_13 = omduh.record_13;
907 
909 }
910 
911 
912 
913 
915 {
916  return (this->dataset_label == omduh.dataset_label &&
917  this->dataset_name == omduh.dataset_name &&
918  this->dataset_location == omduh.dataset_location &&
919  this->id_lines_1_to_5 == omduh.id_lines_1_to_5 &&
920 
921  this->model_type == omduh.model_type &&
922  this->analysis_type == omduh.analysis_type &&
923  this->data_characteristic == omduh.data_characteristic &&
924  this->result_type == omduh.result_type &&
925 
926  this->data_type == omduh.data_type &&
927  this->nvaldc == omduh.nvaldc &&
928 
929  this->record_10 == omduh.record_10 &&
930  this->record_11 == omduh.record_11 &&
931  this->record_12 == omduh.record_12 &&
932  this->record_13 == omduh.record_13 &&
933 
934  this->_desired_dataset_label == omduh._desired_dataset_label);
935 }
936 
937 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo