unstructured_mesh.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 <fstream>
22 #include <sstream>
23 #include <iomanip>
24 
25 // C includes
26 #include <unistd.h> // for unlink()
27 
28 // Local includes
29 #include "libmesh/boundary_info.h"
33 #include "libmesh/elem.h"
34 #include "libmesh/mesh_tools.h" // For n_levels
35 #include "libmesh/parallel.h"
36 #include "libmesh/remote_elem.h"
37 
38 #include "libmesh/diva_io.h"
39 #include "libmesh/exodusII_io.h"
40 #include "libmesh/gmv_io.h"
41 #include "libmesh/tecplot_io.h"
42 #include "libmesh/tetgen_io.h"
43 #include "libmesh/ucd_io.h"
44 #include "libmesh/unv_io.h"
45 #include "libmesh/matlab_io.h"
46 #include "libmesh/off_io.h"
47 #include "libmesh/medit_io.h"
48 #include "libmesh/nemesis_io.h"
49 #include "libmesh/gmsh_io.h"
50 #include "libmesh/fro_io.h"
51 #include "libmesh/xdr_io.h"
52 #include "libmesh/legacy_xdr_io.h"
53 #include "libmesh/vtk_io.h"
54 #include "libmesh/abaqus_io.h"
55 #include "libmesh/checkpoint_io.h"
56 
57 #include LIBMESH_INCLUDE_UNORDERED_MAP
58 
59 
60 
61 // ------------------------------------------------------------
62 // Anonymous namespace for implementation details
63 namespace {
64  bool is_parallel_file_format (const std::string &name)
65  {
66  // Certain mesh formats can support parallel I/O, including the
67  // "new" Xdr format and the Nemesis format.
68  return ((name.rfind(".xda") < name.size()) ||
69  (name.rfind(".xdr") < name.size()) ||
70  (name.rfind(".nem") < name.size()) ||
71  (name.rfind(".n") < name.size()) ||
72  (name.rfind(".cp") < name.size())
73  );
74  }
75 }
76 
77 
78 namespace libMesh
79 {
80 
81 
82 // ------------------------------------------------------------
83 // UnstructuredMesh class member functions
85  unsigned int d) :
86  MeshBase (comm,d)
87 {
89 }
90 
91 
92 
93 #ifndef LIBMESH_DISABLE_COMMWORLD
95  MeshBase (d)
96 {
98 }
99 #endif
100 
101 
102 
104  (const UnstructuredMesh& other_mesh)
105 {
106  // We're assuming our subclass data needs no copy
107  libmesh_assert_equal_to (_n_parts, other_mesh._n_parts);
108  libmesh_assert_equal_to (_dim, other_mesh._dim);
109  libmesh_assert_equal_to (_is_prepared, other_mesh._is_prepared);
110 
111  // We're assuming the other mesh has proper element number ordering,
112  // so that we add parents before their children.
113 #ifdef DEBUG
115 #endif
116 
117  //Copy in Nodes
118  {
119  //Preallocate Memory if necessary
120  this->reserve_nodes(other_mesh.n_nodes());
121 
122  const_node_iterator it = other_mesh.nodes_begin();
123  const_node_iterator end = other_mesh.nodes_end();
124 
125  for (; it != end; ++it)
126  {
127  const Node *oldn = *it;
128 
129  // Add new nodes in old node Point locations
130  /*Node *newn =*/ this->add_point(*oldn, oldn->id(), oldn->processor_id());
131 
132  // And start them off in the same subdomain
133 // newn->processor_id() = oldn->processor_id();
134  }
135  }
136 
137  //Copy in Elements
138  {
139  //Preallocate Memory if necessary
140  this->reserve_elem(other_mesh.n_elem());
141 
142  // Loop over the elements
145 
146  // FIXME: Where do we set element IDs??
147  for (; it != end; ++it)
148  {
149  //Look at the old element
150  const Elem *old = *it;
151  //Build a new element
152  Elem *newparent = old->parent() ?
153  this->elem(old->parent()->id()) : NULL;
154  AutoPtr<Elem> ap = Elem::build(old->type(), newparent);
155  Elem * el = ap.release();
156 
157  el->subdomain_id() = old->subdomain_id();
158 
159  for (unsigned int s=0; s != old->n_sides(); ++s)
160  if (old->neighbor(s) == remote_elem)
161  el->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
162 
163 #ifdef LIBMESH_ENABLE_AMR
164  if (old->has_children())
165  for (unsigned int c=0; c != old->n_children(); ++c)
166  if (old->child(c) == remote_elem)
167  el->add_child(const_cast<RemoteElem*>(remote_elem), c);
168 
169  //Create the parent's child pointers if necessary
170  if (newparent)
171  {
172  unsigned int oldc = old->parent()->which_child_am_i(old);
173  newparent->add_child(el, oldc);
174  }
175 
176  // Copy the refinement flags
177  el->set_refinement_flag(old->refinement_flag());
178  el->set_p_refinement_flag(old->p_refinement_flag());
179 #endif // #ifdef LIBMESH_ENABLE_AMR
180 
181  //Assign all the nodes
182  for(unsigned int i=0;i<el->n_nodes();i++)
183  el->set_node(i) = &this->node(old->node(i));
184 
185  // And start it off in the same subdomain
186  el->processor_id() = old->processor_id();
187 
188  // Give it the same id
189  el->set_id(old->id());
190 
191  //Hold onto it
192  this->add_elem(el);
193  }
194  }
195 
196  //Finally prepare the new Mesh for use. Keep the same numbering and
197  //partitioning but also the same renumbering and partitioning
198  //policies as our source mesh.
199  this->allow_renumbering(false);
200  this->skip_partitioning(true);
201  this->prepare_for_use();
202  this->allow_renumbering(other_mesh.allow_renumbering());
203  this->skip_partitioning(other_mesh.skip_partitioning());
204 }
205 
206 
207 
209 {
210 // this->clear (); // Nothing to clear at this level
211 
213 }
214 
215 
216 
217 
218 
219 void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
220  const bool reset_current_list)
221 {
222  // We might actually want to run this on an empty mesh
223  // (e.g. the boundary mesh for a nonexistant bcid!)
224  // libmesh_assert_not_equal_to (this->n_nodes(), 0);
225  // libmesh_assert_not_equal_to (this->n_elem(), 0);
226 
227  // This function must be run on all processors at once
228  parallel_object_only();
229 
230  START_LOG("find_neighbors()", "Mesh");
231 
232  const element_iterator el_end = this->elements_end();
233 
234  //TODO:[BSK] This should be removed later?!
235  if (reset_current_list)
236  for (element_iterator el = this->elements_begin(); el != el_end; ++el)
237  {
238  Elem* e = *el;
239  for (unsigned int s=0; s<e->n_neighbors(); s++)
240  if (e->neighbor(s) != remote_elem ||
241  reset_remote_elements)
242  e->set_neighbor(s,NULL);
243  }
244 
245  // Find neighboring elements by first finding elements
246  // with identical side keys and then check to see if they
247  // are neighbors
248  {
249  // data structures -- Use the hash_multimap if available
250  typedef unsigned int key_type;
251  typedef std::pair<Elem*, unsigned char> val_type;
252  typedef std::pair<key_type, val_type> key_val_pair;
253 
254  typedef LIBMESH_BEST_UNORDERED_MULTIMAP<key_type, val_type> map_type;
255 
256  // A map from side keys to corresponding elements & side numbers
257  map_type side_to_elem_map;
258 
259 
260 
261  for (element_iterator el = this->elements_begin(); el != el_end; ++el)
262  {
263  Elem* element = *el;
264 
265  for (unsigned int ms=0; ms<element->n_neighbors(); ms++)
266  {
267  next_side:
268  // If we haven't yet found a neighbor on this side, try.
269  // Even if we think our neighbor is remote, that
270  // information may be out of date.
271  if (element->neighbor(ms) == NULL ||
272  element->neighbor(ms) == remote_elem)
273  {
274  // Get the key for the side of this element
275  const unsigned int key = element->key(ms);
276 
277  // Look for elements that have an identical side key
278  std::pair <map_type::iterator, map_type::iterator>
279  bounds = side_to_elem_map.equal_range(key);
280 
281  // May be multiple keys, check all the possible
282  // elements which _might_ be neighbors.
283  if (bounds.first != bounds.second)
284  {
285  // Get the side for this element
286  const AutoPtr<Elem> my_side(element->side(ms));
287 
288  // Look at all the entries with an equivalent key
289  while (bounds.first != bounds.second)
290  {
291  // Get the potential element
292  Elem* neighbor = bounds.first->second.first;
293 
294  // Get the side for the neighboring element
295  const unsigned int ns = bounds.first->second.second;
296  const AutoPtr<Elem> their_side(neighbor->side(ns));
297  //libmesh_assert(my_side.get());
298  //libmesh_assert(their_side.get());
299 
300  // If found a match with my side
301  //
302  // We need special tests here for 1D:
303  // since parents and children have an equal
304  // side (i.e. a node), we need to check
305  // ns != ms, and we also check level() to
306  // avoid setting our neighbor pointer to
307  // any of our neighbor's descendants
308  if( (*my_side == *their_side) &&
309  (element->level() == neighbor->level()) &&
310  ((_dim != 1) || (ns != ms)) )
311  {
312  // So share a side. Is this a mixed pair
313  // of subactive and active/ancestor
314  // elements?
315  // If not, then we're neighbors.
316  // If so, then the subactive's neighbor is
317 
318  if (element->subactive() ==
319  neighbor->subactive())
320  {
321  // an element is only subactive if it has
322  // been coarsened but not deleted
323  element->set_neighbor (ms,neighbor);
324  neighbor->set_neighbor(ns,element);
325  }
326  else if (element->subactive())
327  {
328  element->set_neighbor(ms,neighbor);
329  }
330  else if (neighbor->subactive())
331  {
332  neighbor->set_neighbor(ns,element);
333  }
334  side_to_elem_map.erase (bounds.first);
335 
336  // get out of this nested crap
337  goto next_side;
338  }
339 
340  ++bounds.first;
341  }
342  }
343 
344  // didn't find a match...
345  // Build the map entry for this element
346  key_val_pair kvp;
347 
348  kvp.first = key;
349  kvp.second.first = element;
350  kvp.second.second = ms;
351 
352  // use the lower bound as a hint for
353  // where to put it.
354 #if defined(LIBMESH_HAVE_UNORDERED_MAP) || defined(LIBMESH_HAVE_TR1_UNORDERED_MAP) || defined(LIBMESH_HAVE_HASH_MAP) || defined(LIBMESH_HAVE_EXT_HASH_MAP)
355  side_to_elem_map.insert (kvp);
356 #else
357  side_to_elem_map.insert (bounds.first,kvp);
358 #endif
359  }
360  }
361  }
362  }
363 
364 #ifdef LIBMESH_ENABLE_AMR
365 
386  const unsigned int n_levels = MeshTools::n_levels(*this);
387  for (unsigned int level = 1; level < n_levels; ++level)
388  {
389  element_iterator end = this->level_elements_end(level);
390  for (element_iterator el = this->level_elements_begin(level);
391  el != end; ++el)
392  {
393  Elem* current_elem = *el;
394  libmesh_assert(current_elem);
395  Elem* parent = current_elem->parent();
396  libmesh_assert(parent);
397  const unsigned int my_child_num = parent->which_child_am_i(current_elem);
398 
399  for (unsigned int s=0; s < current_elem->n_neighbors(); s++)
400  {
401  if (current_elem->neighbor(s) == NULL ||
402  (current_elem->neighbor(s) == remote_elem &&
403  parent->is_child_on_side(my_child_num, s)))
404  {
405  Elem *neigh = parent->neighbor(s);
406 
407  // If neigh was refined and had non-subactive children
408  // made remote earlier, then a non-subactive elem should
409  // actually have one of those remote children as a
410  // neighbor
411  if (neigh && (neigh->ancestor()) && (!current_elem->subactive()))
412  {
413 #ifdef DEBUG
414  // Let's make sure that "had children made remote"
415  // situation is actually the case
416  libmesh_assert(neigh->has_children());
417  bool neigh_has_remote_children = false;
418  for (unsigned int c = 0; c != neigh->n_children(); ++c)
419  {
420  if (neigh->child(c) == remote_elem)
421  neigh_has_remote_children = true;
422  }
423  libmesh_assert(neigh_has_remote_children);
424 
425  // And let's double-check that we don't have
426  // a remote_elem neighboring a local element
427  libmesh_assert_not_equal_to (current_elem->processor_id(),
428  this->processor_id());
429 #endif // DEBUG
430  neigh = const_cast<RemoteElem*>(remote_elem);
431  }
432 
433  current_elem->set_neighbor(s, neigh);
434 #ifdef DEBUG
435  if (neigh != NULL && neigh != remote_elem)
436  // We ignore subactive elements here because
437  // we don't care about neighbors of subactive element.
438  if ((!neigh->active()) && (!current_elem->subactive()))
439  {
440  libMesh::err << "On processor " << this->processor_id()
441  << std::endl;
442  libMesh::err << "Bad element ID = " << current_elem->id()
443  << ", Side " << s << ", Bad neighbor ID = " << neigh->id() << std::endl;
444  libMesh::err << "Bad element proc_ID = " << current_elem->processor_id()
445  << ", Bad neighbor proc_ID = " << neigh->processor_id() << std::endl;
446  libMesh::err << "Bad element size = " << current_elem->hmin()
447  << ", Bad neighbor size = " << neigh->hmin() << std::endl;
448  libMesh::err << "Bad element center = " << current_elem->centroid()
449  << ", Bad neighbor center = " << neigh->centroid() << std::endl;
450  libMesh::err << "ERROR: "
451  << (current_elem->active()?"Active":"Ancestor")
452  << " Element at level "
453  << current_elem->level() << std::endl;
454  libMesh::err << "with "
455  << (parent->active()?"active":
456  (parent->subactive()?"subactive":"ancestor"))
457  << " parent share "
458  << (neigh->subactive()?"subactive":"ancestor")
459  << " neighbor at level " << neigh->level()
460  << std::endl;
461  GMVIO(*this).write ("bad_mesh.gmv");
462  libmesh_error();
463  }
464 #endif // DEBUG
465  }
466  }
467  }
468  }
469 
470 #endif // AMR
471 
472 
473 #ifdef DEBUG
475 #endif
476 
477  STOP_LOG("find_neighbors()", "Mesh");
478 }
479 
480 
481 
482 void UnstructuredMesh::read (const std::string& name,
483  MeshData* mesh_data,
484  bool skip_renumber_nodes_and_elements)
485 {
486  // See if the file exists. Perform this check on all processors
487  // so that the code is terminated properly in the case that the
488  // file does not exist.
489 
490  // For Nemesis files, the name we try to read will have suffixes
491  // identifying processor rank
492  if (name.rfind(".nem") + 4 == name.size() ||
493  name.rfind(".n") + 2 == name.size())
494  {
495  std::ostringstream full_name;
496 
497  // Find the length of a string which represents the highest processor ID
498  full_name << (this->n_processors());
499  unsigned field_width = full_name.str().size();
500 
501  // reset the string stream
502  full_name.str("");
503 
504  // And build up the full filename
505  full_name << name
506  << '.' << this->n_processors()
507  << '.' << std::setfill('0') << std::setw(field_width) << this->processor_id();
508 
509  std::ifstream in (full_name.str().c_str());
510 
511  if (!in.good())
512  {
513  libMesh::err << "ERROR: cannot locate specified file:\n\t"
514  << full_name.str()
515  << std::endl;
516  libmesh_error();
517  }
518  }
519  else if(name.rfind(".cp")) {} // Do error checking in the reader
520  else
521  {
522  std::ifstream in (name.c_str());
523 
524  if (!in.good())
525  {
526  libMesh::err << "ERROR: cannot locate specified file:\n\t"
527  << name
528  << std::endl;
529  libmesh_error();
530  }
531  }
532 
533  // Set the skip_renumber_nodes_and_elements flag on all processors.
534  // This ensures that renumber_nodes_and_elements is *not* called
535  // during prepare_for_use() for certain types of mesh files.
536  // This is required in cases where there is an associated solution
537  // file which expects a certain ordering of the nodes.
538  if(name.rfind(".gmv")+4==name.size())
539  {
540  skip_renumber_nodes_and_elements = true;
541  }
542 
543  // Look for parallel formats first
544  if (is_parallel_file_format(name))
545  {
546  // no need to handle bz2 files here -- the Xdr class does that.
547  if ((name.rfind(".xda") < name.size()) ||
548  (name.rfind(".xdr") < name.size()))
549  {
550  XdrIO xdr_io(*this);
551 
552  // .xda* ==> bzip2/gzip/ASCII flavors
553  if (name.rfind(".xda") < name.size())
554  {
555  xdr_io.binary() = false;
556  xdr_io.read (name);
557  }
558  else // .xdr* ==> true binary XDR file
559  {
560  xdr_io.binary() = true;
561  xdr_io.read (name);
562  }
563 
564  // The xdr_io object gets constructed with legacy() == false.
565  // if legacy() == true then it means that a legacy file was detected and
566  // thus processor 0 performed the read. We therefore need to broadcast the
567  // mesh. Further, for this flavor of mesh solution data ordering is tied
568  // to the node ordering, so we better not reorder the nodes!
569  if (xdr_io.legacy())
570  {
571  this->allow_renumbering(false);
572  MeshCommunication().broadcast(*this);
573  }
574 
575  // libHilbert-enabled libMesh builds should construct files
576  // with a canonical node ordering, which libHilbert-enabled
577  // builds will be able to read in again regardless of any
578  // renumbering. So in that case we're free to renumber.
579  // However, if either the writer or the reader of this file
580  // don't have libHilbert, then we'll have to skip
581  // renumbering because we need the numbering to remain
582  // consistent with any solution file we read in next.
583 #ifdef LIBMESH_HAVE_LIBHILBERT
584  // if (!xdr_io.libhilbert_ordering())
585  // skip_renumber_nodes_and_elements = true;
586 #else
587  this->allow_renumbering(false);
588 #endif
589  }
590  else if (name.rfind(".nem") < name.size() ||
591  name.rfind(".n") < name.size())
592  Nemesis_IO(*this).read (name);
593  else if (name.rfind(".cp") < name.size())
594  {
595  if(name.rfind(".cpa") < name.size())
596  CheckpointIO(*this, false).read(name);
597  else
598  CheckpointIO(*this, true).read(name);
599  }
600  }
601 
602  // Serial mesh formats
603  else
604  {
605  START_LOG("read()", "Mesh");
606 
607  // Read the file based on extension. Only processor 0
608  // needs to read the mesh. It will then broadcast it and
609  // the other processors will pick it up
610  if (this->processor_id() == 0)
611  {
612  std::ostringstream pid_suffix;
613  pid_suffix << '_' << getpid();
614  // Nasty hack for reading/writing zipped files
615  std::string new_name = name;
616  if (name.size() - name.rfind(".bz2") == 4)
617  {
618 #ifdef LIBMESH_HAVE_BZIP
619  new_name.erase(new_name.end() - 4, new_name.end());
620  new_name += pid_suffix.str();
621  std::string system_string = "bunzip2 -f -k -c ";
622  system_string += name + " > " + new_name;
623  START_LOG("system(bunzip2)", "Mesh");
624  if (std::system(system_string.c_str()))
625  libmesh_file_error(system_string);
626  STOP_LOG("system(bunzip2)", "Mesh");
627 #else
628  libMesh::err << "ERROR: need bzip2/bunzip2 to open .bz2 file "
629  << name << std::endl;
630  libmesh_error();
631 #endif
632  }
633  else if (name.size() - name.rfind(".xz") == 3)
634  {
635 #ifdef LIBMESH_HAVE_XZ
636  new_name.erase(new_name.end() - 3, new_name.end());
637  new_name += pid_suffix.str();
638  std::string system_string = "xz -f -d -k -c ";
639  system_string += name + " > " + new_name;
640  START_LOG("system(xz -d)", "XdrIO");
641  if (std::system(system_string.c_str()))
642  libmesh_file_error(system_string);
643  STOP_LOG("system(xz -d)", "XdrIO");
644 #else
645  libMesh::err << "ERROR: need xz to open .xz file "
646  << name << std::endl;
647  libmesh_error();
648 #endif
649  }
650 
651  if (new_name.rfind(".mat") < new_name.size())
652  MatlabIO(*this).read(new_name);
653 
654  else if (new_name.rfind(".ucd") < new_name.size())
655  UCDIO(*this).read (new_name);
656 
657  else if ((new_name.rfind(".off") < new_name.size()) ||
658  (new_name.rfind(".ogl") < new_name.size()) ||
659  (new_name.rfind(".oogl") < new_name.size()))
660  OFFIO(*this).read (new_name);
661 
662  else if (new_name.rfind(".mgf") < new_name.size())
663  LegacyXdrIO(*this,true).read_mgf (new_name);
664 
665  else if (new_name.rfind(".unv") < new_name.size())
666  {
667  if (mesh_data == NULL)
668  {
669  libMesh::err << "Error! You must pass a "
670  << "valid MeshData pointer to "
671  << "read UNV files!" << std::endl;
672  libmesh_error();
673  }
674  UNVIO(*this, *mesh_data).read (new_name);
675  }
676 
677  else if ((new_name.rfind(".node") < new_name.size()) ||
678  (new_name.rfind(".ele") < new_name.size()))
679  TetGenIO(*this,mesh_data).read (new_name);
680 
681  else if (new_name.rfind(".exd") < new_name.size() ||
682  new_name.rfind(".e") < new_name.size())
683  ExodusII_IO(*this).read (new_name);
684 
685  else if (new_name.rfind(".msh") < new_name.size())
686  GmshIO(*this).read (new_name);
687 
688  else if (new_name.rfind(".gmv") < new_name.size())
689  GMVIO(*this).read (new_name);
690 
691  else if (new_name.rfind(".vtu") < new_name.size())
692  VTKIO(*this).read(new_name);
693 
694  else if (new_name.rfind(".inp") < new_name.size())
695  AbaqusIO(*this).read(new_name);
696 
697  else
698  {
699  libMesh::err << " ERROR: Unrecognized file extension: " << name
700  << "\n I understand the following:\n\n"
701  << " *.e -- Sandia's ExodusII format\n"
702  << " *.exd -- Sandia's ExodusII format\n"
703  << " *.gmv -- LANL's General Mesh Viewer format\n"
704  << " *.mat -- Matlab triangular ASCII file\n"
705  << " *.n -- Sandia's Nemesis format\n"
706  << " *.nem -- Sandia's Nemesis format\n"
707  << " *.off -- OOGL OFF surface format\n"
708  << " *.ucd -- AVS's ASCII UCD format\n"
709  << " *.unv -- I-deas Universal format\n"
710  << " *.vtu -- Paraview VTK format\n"
711  << " *.inp -- Abaqus .inp format\n"
712  << " *.xda -- libMesh ASCII format\n"
713  << " *.xdr -- libMesh binary format\n"
714  << " *.gz -- any above format gzipped\n"
715  << " *.bz2 -- any above format bzip2'ed\n"
716  << " *.xz -- any above format xzipped\n"
717  << " *.cpa -- libMesh Checkpoint ASCII format\n"
718  << " *.cpr -- libMesh Checkpoint binary format\n"
719 
720  << std::endl;
721  libmesh_error();
722  }
723 
724  // If we temporarily decompressed a file, remove the
725  // uncompressed version
726  if (name.size() - name.rfind(".bz2") == 4)
727  std::remove(new_name.c_str());
728  if (name.size() - name.rfind(".xz") == 3)
729  std::remove(new_name.c_str());
730  }
731 
732 
733  STOP_LOG("read()", "Mesh");
734 
735  // Send the mesh & bcs (which are now only on processor 0) to the other
736  // processors
737  MeshCommunication().broadcast (*this);
738  }
739 
740  if (skip_renumber_nodes_and_elements)
741  {
742  // Use MeshBase::allow_renumbering() yourself instead.
743  libmesh_deprecated();
744  this->allow_renumbering(false);
745  }
746 
747  // Done reading the mesh. Now prepare it for use.
748  this->prepare_for_use();
749 }
750 
751 
752 
753 void UnstructuredMesh::write (const std::string& name,
754  MeshData* mesh_data)
755 {
756  // parallel formats are special -- they may choose to write
757  // separate files, let's not try to handle the zipping here.
758  if (is_parallel_file_format(name))
759  {
760  // no need to handle bz2 files here -- the Xdr class does that.
761  if (name.rfind(".xda") < name.size())
762  XdrIO(*this).write(name);
763 
764  else if (name.rfind(".xdr") < name.size())
765  XdrIO(*this,true).write(name);
766 
767  else if (name.rfind(".nem") < name.size() ||
768  name.rfind(".n") < name.size())
769  Nemesis_IO(*this).write(name);
770  }
771 
772  // serial file formats
773  else
774  {
775  START_LOG("write()", "Mesh");
776 
777  // Nasty hack for reading/writing zipped files
778  std::string new_name = name;
779  processor_id_type pid_0 = 0;
780  if (this->processor_id() == 0)
781  pid_0 = getpid();
782  this->comm().broadcast(pid_0);
783  std::ostringstream pid_suffix;
784  pid_suffix << '_' << pid_0;
785 
786  if (name.size() - name.rfind(".bz2") == 4)
787  {
788  new_name.erase(new_name.end() - 4, new_name.end());
789  new_name += pid_suffix.str();
790  }
791  else if (name.size() - name.rfind(".xz") == 3)
792  {
793  new_name.erase(new_name.end() - 3, new_name.end());
794  new_name += pid_suffix.str();
795  }
796 
797  // New scope so that io will close before we try to zip the file
798  {
799  // Write the file based on extension
800  if (new_name.rfind(".dat") < new_name.size())
801  TecplotIO(*this).write (new_name);
802 
803  else if (new_name.rfind(".plt") < new_name.size())
804  TecplotIO(*this,true).write (new_name);
805 
806  else if (new_name.rfind(".ucd") < new_name.size())
807  UCDIO (*this).write (new_name);
808 
809  else if (new_name.rfind(".gmv") < new_name.size())
810  if (this->n_partitions() > 1)
811  GMVIO(*this).write (new_name);
812  else
813  {
814  GMVIO io(*this);
815  io.partitioning() = false;
816  io.write (new_name);
817  }
818 
819  else if (new_name.rfind(".ugrid") < new_name.size())
820  DivaIO(*this).write(new_name);
821  else if (new_name.rfind(".exd") < new_name.size() ||
822  new_name.rfind(".e") < new_name.size())
823  ExodusII_IO(*this).write(new_name);
824  else if (new_name.rfind(".mgf") < new_name.size())
825  LegacyXdrIO(*this,true).write_mgf(new_name);
826 
827  else if (new_name.rfind(".unv") < new_name.size())
828  {
829  if (mesh_data == NULL)
830  {
831  libMesh::err << "Error! You must pass a "
832  << "valid MeshData pointer to "
833  << "write UNV files!" << std::endl;
834  libmesh_error();
835  }
836  UNVIO(*this, *mesh_data).write (new_name);
837  }
838 
839  else if (new_name.rfind(".mesh") < new_name.size())
840  MEDITIO(*this).write (new_name);
841 
842  else if (new_name.rfind(".poly") < new_name.size())
843  TetGenIO(*this).write (new_name);
844 
845  else if (new_name.rfind(".msh") < new_name.size())
846  GmshIO(*this).write (new_name);
847 
848  else if (new_name.rfind(".fro") < new_name.size())
849  FroIO(*this).write (new_name);
850 
851  else if (new_name.rfind(".vtu") < new_name.size())
852  VTKIO(*this).write (new_name);
853 
854  else
855  {
857  << " ERROR: Unrecognized file extension: " << name
858  << "\n I understand the following:\n\n"
859  << " *.dat -- Tecplot ASCII file\n"
860  << " *.e -- Sandia's ExodusII format\n"
861  << " *.exd -- Sandia's ExodusII format\n"
862  << " *.fro -- ACDL's surface triangulation file\n"
863  << " *.gmv -- LANL's GMV (General Mesh Viewer) format\n"
864  << " *.mesh -- MEdit mesh format\n"
865  << " *.mgf -- MGF binary mesh format\n"
866  << " *.msh -- GMSH ASCII file\n"
867  << " *.n -- Sandia's Nemesis format\n"
868  << " *.nem -- Sandia's Nemesis format\n"
869  << " *.plt -- Tecplot binary file\n"
870  << " *.poly -- TetGen ASCII file\n"
871  << " *.ucd -- AVS's ASCII UCD format\n"
872  << " *.ugrid -- Kelly's DIVA ASCII format\n"
873  << " *.unv -- I-deas Universal format\n"
874  << " *.vtu -- VTK (paraview-readable) format\n"
875  << " *.xda -- libMesh ASCII format\n"
876  << " *.xdr -- libMesh binary format,\n"
877  << std::endl
878  << "\n Exiting without writing output\n";
879  }
880  }
881 
882  // Nasty hack for reading/writing zipped files
883  if (name.size() - name.rfind(".bz2") == 4)
884  {
885  START_LOG("system(bzip2)", "Mesh");
886  if (this->processor_id() == 0)
887  {
888  std::string system_string = "bzip2 -f -c ";
889  system_string += new_name + " > " + name;
890  if (std::system(system_string.c_str()))
891  libmesh_file_error(system_string);
892  std::remove(new_name.c_str());
893  }
894  this->comm().barrier();
895  STOP_LOG("system(bzip2)", "Mesh");
896  }
897  if (name.size() - name.rfind(".xz") == 3)
898  {
899  START_LOG("system(xz)", "Mesh");
900  if (this->processor_id() == 0)
901  {
902  std::string system_string = "xz -f -c ";
903  system_string += new_name + " > " + name;
904  if (std::system(system_string.c_str()))
905  libmesh_file_error(system_string);
906  std::remove(new_name.c_str());
907  }
908  this->comm().barrier();
909  STOP_LOG("system(xz)", "Mesh");
910  }
911 
912  STOP_LOG("write()", "Mesh");
913  }
914 }
915 
916 
917 
918 void UnstructuredMesh::write (const std::string& name,
919  const std::vector<Number>& v,
920  const std::vector<std::string>& vn)
921 {
922  START_LOG("write()", "Mesh");
923 
924  // Write the file based on extension
925  if (name.rfind(".dat") < name.size())
926  TecplotIO(*this).write_nodal_data (name, v, vn);
927 
928  else if (name.rfind(".plt") < name.size())
929  TecplotIO(*this,true).write_nodal_data (name, v, vn);
930 
931  else if (name.rfind(".gmv") < name.size())
932  {
933  if (n_subdomains() > 1)
934  GMVIO(*this).write_nodal_data (name, v, vn);
935  else
936  {
937  GMVIO io(*this);
938  io.partitioning() = false;
939  io.write_nodal_data (name, v, vn);
940  }
941  }
942  else if (name.rfind(".pvtu") < name.size())
943  {
944  VTKIO(*this).write_nodal_data (name, v, vn);
945  }
946  else
947  {
949  << " ERROR: Unrecognized file extension: " << name
950  << "\n I understand the following:\n\n"
951  << " *.dat -- Tecplot ASCII file\n"
952  << " *.gmv -- LANL's GMV (General Mesh Viewer) format\n"
953  << " *.plt -- Tecplot binary file\n"
954  << " *.pvtu -- Paraview VTK file\n"
955  << "\n Exiting without writing output\n";
956  }
957 
958  STOP_LOG("write()", "Mesh");
959 }
960 
961 
962 
963 
964 
966  const processor_id_type pid) const
967 {
968 
969  // Issue a warning if the number the number of processors
970  // currently available is less that that requested for
971  // partitioning. This is not necessarily an error since
972  // you may run on one processor and still partition the
973  // mesh into several partitions.
974 #ifdef DEBUG
975  if (this->n_processors() < pid)
976  {
977  libMesh::out << "WARNING: You are creating a "
978  << "mesh for a processor id (="
979  << pid
980  << ") greater than "
981  << "the number of processors available for "
982  << "the calculation. (="
983  << this->n_processors()
984  << ")."
985  << std::endl;
986  }
987 #endif
988 
989  // Create iterators to loop over the list of elements
990 // const_active_pid_elem_iterator it(this->elements_begin(), pid);
991 // const const_active_pid_elem_iterator it_end(this->elements_end(), pid);
992 
994  const const_element_iterator it_end = this->active_pid_elements_end(pid);
995 
996  this->create_submesh (pid_mesh, it, it_end);
997 }
998 
999 
1000 
1001 
1002 
1003 
1004 
1007  const const_element_iterator& it_end) const
1008 {
1009  // Just in case the subdomain_mesh already has some information
1010  // in it, get rid of it.
1011  new_mesh.clear();
1012 
1013  // Fail if (*this == new_mesh), we cannot create a submesh inside ourself!
1014  // This may happen if the user accidently passes the original mesh into
1015  // this function! We will check this by making sure we did not just
1016  // clear ourself.
1017  libmesh_assert_not_equal_to (this->n_nodes(), 0);
1018  libmesh_assert_not_equal_to (this->n_elem(), 0);
1019 
1020  // How the nodes on this mesh will be renumbered to nodes
1021  // on the new_mesh.
1022  std::vector<dof_id_type> new_node_numbers (this->n_nodes());
1023 
1024  std::fill (new_node_numbers.begin(),
1025  new_node_numbers.end(),
1027 
1028 
1029 
1030  // the number of nodes on the new mesh, will be incremented
1031  dof_id_type n_new_nodes = 0;
1032  dof_id_type n_new_elem = 0;
1033 
1034  for (; it != it_end; ++it)
1035  {
1036  // increment the new element counter
1037  n_new_elem++;
1038 
1039  const Elem* old_elem = *it;
1040 
1041  // Add an equivalent element type to the new_mesh
1042  Elem* new_elem =
1043  new_mesh.add_elem (Elem::build(old_elem->type()).release());
1044 
1045  libmesh_assert(new_elem);
1046 
1047  // Loop over the nodes on this element.
1048  for (unsigned int n=0; n<old_elem->n_nodes(); n++)
1049  {
1050  libmesh_assert_less (old_elem->node(n), new_node_numbers.size());
1051 
1052  if (new_node_numbers[old_elem->node(n)] == DofObject::invalid_id)
1053  {
1054  new_node_numbers[old_elem->node(n)] = n_new_nodes;
1055 
1056  // Add this node to the new mesh
1057  new_mesh.add_point (old_elem->point(n));
1058 
1059  // Increment the new node counter
1060  n_new_nodes++;
1061  }
1062 
1063  // Define this element's connectivity on the new mesh
1064  libmesh_assert_less (new_node_numbers[old_elem->node(n)], new_mesh.n_nodes());
1065 
1066  new_elem->set_node(n) = new_mesh.node_ptr (new_node_numbers[old_elem->node(n)]);
1067  }
1068 
1069  // Copy ids for this element
1070  new_elem->subdomain_id() = old_elem->subdomain_id();
1071  new_elem->processor_id() = old_elem->processor_id();
1072 
1073  // Maybe add boundary conditions for this element
1074  for (unsigned int s=0; s<old_elem->n_sides(); s++)
1075 // We're supporting boundary ids on internal sides now
1076 // if (old_elem->neighbor(s) == NULL)
1077  {
1078  const std::vector<boundary_id_type>& bc_ids = this->boundary_info->boundary_ids(old_elem, s);
1079  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1080  {
1081  const boundary_id_type bc_id = *id_it;
1082  if (bc_id != this->boundary_info->invalid_id)
1083  new_mesh.boundary_info->add_side (new_elem,
1084  s,
1085  bc_id);
1086  }
1087  }
1088  } // end loop over elements
1089 
1090 
1091  // Prepare the new_mesh for use
1092  new_mesh.prepare_for_use(/*skip_renumber =*/false);
1093 
1094 }
1095 
1096 
1097 
1098 #ifdef LIBMESH_ENABLE_AMR
1100 {
1101  START_LOG ("contract()", "Mesh");
1102 
1103  // Flag indicating if this call actually changes the mesh
1104  bool mesh_changed = false;
1105 
1107  const element_iterator end = elements_end();
1108 
1109 #ifdef DEBUG
1110  for ( ; in != end; ++in)
1111  if (*in != NULL)
1112  {
1113  Elem* el = *in;
1114  libmesh_assert(el->active() || el->subactive() || el->ancestor());
1115  }
1116  in = elements_begin();
1117 #endif
1118 
1119  // Loop over the elements.
1120  for ( ; in != end; ++in)
1121  if (*in != NULL)
1122  {
1123  Elem* el = *in;
1124 
1125  // Delete all the subactive ones
1126  if (el->subactive())
1127  {
1128  // No level-0 element should be subactive.
1129  // Note that we CAN'T test elem->level(), as that
1130  // touches elem->parent()->dim(), and elem->parent()
1131  // might have already been deleted!
1132  libmesh_assert(el->parent());
1133 
1134  // Delete the element
1135  // This just sets a pointer to NULL, and doesn't
1136  // invalidate any iterators
1137  this->delete_elem(el);
1138 
1139  // the mesh has certainly changed
1140  mesh_changed = true;
1141  }
1142  else
1143  {
1144  // Compress all the active ones
1145  if (el->active())
1146  el->contract();
1147  else
1148  libmesh_assert (el->ancestor());
1149  }
1150  }
1151 
1152  // Strip any newly-created NULL voids out of the element array
1154 
1155  // FIXME: Need to understand why deleting subactive children
1156  // invalidates the point locator. For now we will clear it explicitly
1157  this->clear_point_locator();
1158 
1159  STOP_LOG ("contract()", "Mesh");
1160 
1161  return mesh_changed;
1162 }
1163 #endif // #ifdef LIBMESH_ENABLE_AMR
1164 
1165 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo