mesh_generation.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 <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::sqrt
23 
24 
25 // Local includes
28 // #include "libmesh/elem.h"
30 #include "libmesh/edge_edge2.h"
31 #include "libmesh/edge_edge3.h"
32 #include "libmesh/edge_edge4.h"
33 #include "libmesh/face_tri3.h"
34 #include "libmesh/face_tri6.h"
35 #include "libmesh/face_quad4.h"
36 #include "libmesh/face_quad8.h"
37 #include "libmesh/face_quad9.h"
38 #include "libmesh/cell_hex8.h"
39 #include "libmesh/cell_hex20.h"
40 #include "libmesh/cell_hex27.h"
41 #include "libmesh/cell_prism6.h"
42 #include "libmesh/cell_prism15.h"
43 #include "libmesh/cell_prism18.h"
44 #include "libmesh/cell_tet4.h"
45 #include "libmesh/cell_pyramid5.h"
47 #include "libmesh/boundary_info.h"
48 #include "libmesh/sphere.h"
51 #include "libmesh/node_elem.h"
52 #include "libmesh/vector_value.h"
53 
54 namespace libMesh
55 {
56 
57 namespace MeshTools {
58  namespace Generation {
59  namespace Private {
67  inline
68  unsigned int idx(const ElemType type,
69  const unsigned int nx,
70  const unsigned int i,
71  const unsigned int j)
72  {
73  switch(type)
74  {
75  case INVALID_ELEM:
76  case QUAD4:
77  case TRI3:
78  {
79  return i + j*(nx+1);
80  break;
81  }
82 
83  case QUAD8:
84  case QUAD9:
85  case TRI6:
86  {
87  return i + j*(2*nx+1);
88  break;
89  }
90 
91  default:
92  {
93  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
94  libmesh_error();
95  }
96  }
97 
98  return libMesh::invalid_uint;
99  }
100 
101 
102 
103  // Same as the function above, but for 3D elements
104  inline
105  unsigned int idx(const ElemType type,
106  const unsigned int nx,
107  const unsigned int ny,
108  const unsigned int i,
109  const unsigned int j,
110  const unsigned int k)
111  {
112  switch(type)
113  {
114  case INVALID_ELEM:
115  case HEX8:
116  case PRISM6:
117  {
118  return i + (nx+1)*(j + k*(ny+1));
119  break;
120  }
121 
122  case HEX20:
123  case HEX27:
124  case TET4: // TET4's are created from an initial HEX27 discretization
125  case TET10: // TET10's are created from an initial HEX27 discretization
126  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
127  case PYRAMID14:
128  case PRISM15:
129  case PRISM18:
130  {
131  return i + (2*nx+1)*(j + k*(2*ny+1));
132  break;
133  }
134 
135  default:
136  {
137  libMesh::err << "ERROR: Unrecognized element type." << std::endl;
138  libmesh_error();
139  }
140  }
141 
142  return libMesh::invalid_uint;
143  }
144  } // namespace Private
145  } // namespace Generation
146 } // namespace MeshTools
147 
148 // ------------------------------------------------------------
149 // MeshTools::Generation function for mesh generation
151  const unsigned int nx,
152  const unsigned int ny,
153  const unsigned int nz,
154  const Real xmin, const Real xmax,
155  const Real ymin, const Real ymax,
156  const Real zmin, const Real zmax,
157  const ElemType type,
158  const bool gauss_lobatto_grid)
159 {
160  START_LOG("build_cube()", "MeshTools::Generation");
161 
162  // Declare that we are using the indexing utility routine
163  // in the "Private" part of our current namespace. If this doesn't
164  // work in GCC 2.95.3 we can either remove it or stop supporting
165  // 2.95.3 altogether.
166  // Changing this to import the whole namespace... just importing idx
167  // causes an internal compiler error for Intel Compiler 11.0 on Linux
168  // in debug mode.
169  using namespace MeshTools::Generation::Private;
170 
171  // Clear the mesh and start from scratch
172  mesh.clear();
173 
174  if (nz != 0)
175  mesh.set_mesh_dimension(3);
176  else if (ny != 0)
177  mesh.set_mesh_dimension(2);
178  else if (nx != 0)
179  mesh.set_mesh_dimension(1);
180  else
181  mesh.set_mesh_dimension(0);
182 
183  switch (mesh.mesh_dimension())
184  {
185  //---------------------------------------------------------------------
186  // Build a 0D point
187  case 0:
188  {
189  libmesh_assert_equal_to (nx, 0);
190  libmesh_assert_equal_to (ny, 0);
191  libmesh_assert_equal_to (nz, 0);
192 
193  libmesh_assert (type == INVALID_ELEM || type == NODEELEM);
194 
195  // Build one nodal element for the mesh
196  mesh.add_point (Point(0, 0, 0), 0);
197  Elem* elem = mesh.add_elem (new NodeElem);
198  elem->set_node(0) = mesh.node_ptr(0);
199 
200  break;
201  }
202 
203 
204 
205  //---------------------------------------------------------------------
206  // Build a 1D line
207  case 1:
208  {
209  libmesh_assert_not_equal_to (nx, 0);
210  libmesh_assert_equal_to (ny, 0);
211  libmesh_assert_equal_to (nz, 0);
212  libmesh_assert_less (xmin, xmax);
213 
214  // Reserve elements
215  switch (type)
216  {
217  case INVALID_ELEM:
218  case EDGE2:
219  case EDGE3:
220  case EDGE4:
221  {
222  mesh.reserve_elem (nx);
223  break;
224  }
225 
226  default:
227  {
228  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
229  libmesh_error();
230  }
231  }
232 
233  // Reserve nodes
234  switch (type)
235  {
236  case INVALID_ELEM:
237  case EDGE2:
238  {
239  mesh.reserve_nodes(nx+1);
240  break;
241  }
242 
243  case EDGE3:
244  {
245  mesh.reserve_nodes(2*nx+1);
246  break;
247  }
248 
249  case EDGE4:
250  {
251  mesh.reserve_nodes(3*nx+1);
252  break;
253  }
254 
255  default:
256  {
257  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
258  libmesh_error();
259  }
260  }
261 
262 
263  // Build the nodes, depends on whether we're using linears,
264  // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
265  unsigned int node_id = 0;
266  switch(type)
267  {
268  case INVALID_ELEM:
269  case EDGE2:
270  {
271  for (unsigned int i=0; i<=nx; i++)
272  {
273  if (gauss_lobatto_grid)
274  mesh.add_point (Point(0.5*(std::cos(libMesh::pi*static_cast<Real>(nx-i)/static_cast<Real>(nx))+1.0),
275  0,
276  0), node_id++);
277  else
278  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
279  0,
280  0), node_id++);
281  }
282  break;
283  }
284 
285  case EDGE3:
286  {
287  for (unsigned int i=0; i<=2*nx; i++)
288  {
289  if (gauss_lobatto_grid)
290  {
291  // The x location of the point.
292  Real x=0.;
293 
294  // Shortcut quantities (do not depend on i)
295  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
296 
297  // If i is even, compute a normal Gauss-Lobatto point
298  if (i%2 == 0)
299  x = 0.5*(1.0 - c);
300 
301  // Otherwise, it is the average of the previous and next points
302  else
303  {
304  Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(2*nx) );
305  Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(2*nx) );
306 
307  Real gl_xmin = 0.5*(1.0 - cmin);
308  Real gl_xmax = 0.5*(1.0 - cmax);
309  x = 0.5*(gl_xmin + gl_xmax);
310  }
311 
312  mesh.add_point (Point(x,0.,0.), node_id++);
313  }
314  else
315  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
316  0,
317  0), node_id++);
318  }
319  break;
320  }
321 
322  case EDGE4:
323  {
324  for (unsigned int i=0; i<=3*nx; i++)
325  {
326  if (gauss_lobatto_grid)
327  {
328  // The x location of the point
329  Real x=0.;
330 
331  // Shortcut quantities
332  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(3*nx) );
333 
334  // If i is multiple of 3, compute a normal Gauss-Lobatto point
335  if (i%3 == 0)
336  x = 0.5*(1.0 - c);
337 
338  // Otherwise, distribute points evenly within the element
339  else
340  {
341  if(i%3 == 1)
342  {
343  Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(3*nx) );
344  Real cmax = std::cos( libMesh::pi*(i+2) / static_cast<Real>(3*nx) );
345 
346  Real gl_xmin = 0.5*(1.0 - cmin);
347  Real gl_xmax = 0.5*(1.0 - cmax);
348 
349  x = (2.*gl_xmin + gl_xmax)/3.;
350  }
351  else
352  if(i%3 == 2)
353  {
354  Real cmin = std::cos( libMesh::pi*(i-2) / static_cast<Real>(3*nx) );
355  Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(3*nx) );
356 
357  Real gl_xmin = 0.5*(1.0 - cmin);
358  Real gl_xmax = 0.5*(1.0 - cmax);
359 
360  x = (gl_xmin + 2.*gl_xmax)/3.;
361  }
362 
363  }
364 
365  mesh.add_point (Point(x,0.,0.), node_id++);
366  }
367  else
368  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(3*nx),
369  0,
370  0), node_id++);
371  }
372 
373 
374 
375  break;
376  }
377 
378  default:
379  {
380  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
381  libmesh_error();
382  }
383 
384  }
385 
386  // Build the elements of the mesh
387  switch(type)
388  {
389  case INVALID_ELEM:
390  case EDGE2:
391  {
392  for (unsigned int i=0; i<nx; i++)
393  {
394  Elem* elem = mesh.add_elem (new Edge2);
395  elem->set_node(0) = mesh.node_ptr(i);
396  elem->set_node(1) = mesh.node_ptr(i+1);
397 
398  if (i == 0)
399  mesh.boundary_info->add_side(elem, 0, 0);
400 
401  if (i == (nx-1))
402  mesh.boundary_info->add_side(elem, 1, 1);
403  }
404  break;
405  }
406 
407  case EDGE3:
408  {
409  for (unsigned int i=0; i<nx; i++)
410  {
411  Elem* elem = mesh.add_elem (new Edge3);
412  elem->set_node(0) = mesh.node_ptr(2*i);
413  elem->set_node(2) = mesh.node_ptr(2*i+1);
414  elem->set_node(1) = mesh.node_ptr(2*i+2);
415 
416  if (i == 0)
417  mesh.boundary_info->add_side(elem, 0, 0);
418 
419  if (i == (nx-1))
420  mesh.boundary_info->add_side(elem, 1, 1);
421  }
422  break;
423  }
424 
425  case EDGE4:
426  {
427  for (unsigned int i=0; i<nx; i++)
428  {
429  Elem* elem = mesh.add_elem (new Edge4);
430  elem->set_node(0) = mesh.node_ptr(3*i);
431  elem->set_node(2) = mesh.node_ptr(3*i+1);
432  elem->set_node(3) = mesh.node_ptr(3*i+2);
433  elem->set_node(1) = mesh.node_ptr(3*i+3);
434 
435  if (i == 0)
436  mesh.boundary_info->add_side(elem, 0, 0);
437 
438  if (i == (nx-1))
439  mesh.boundary_info->add_side(elem, 1, 1);
440  }
441  break;
442  }
443 
444  default:
445  {
446  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
447  libmesh_error();
448  }
449  }
450 
451  // Scale the nodal positions
452  for (unsigned int p=0; p<mesh.n_nodes(); p++)
453  mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
454 
455  // Add sideset names to boundary info
456  mesh.boundary_info->sideset_name(0) = "left";
457  mesh.boundary_info->sideset_name(1) = "right";
458 
459  // Add nodeset names to boundary info
460  mesh.boundary_info->nodeset_name(0) = "left";
461  mesh.boundary_info->nodeset_name(1) = "right";
462 
463  break;
464  }
465 
466 
467 
468 
469 
470 
471 
472 
473 
474 
475  //---------------------------------------------------------------------
476  // Build a 2D quadrilateral
477  case 2:
478  {
479  libmesh_assert_not_equal_to (nx, 0);
480  libmesh_assert_not_equal_to (ny, 0);
481  libmesh_assert_equal_to (nz, 0);
482  libmesh_assert_less (xmin, xmax);
483  libmesh_assert_less (ymin, ymax);
484 
485  // Reserve elements. The TRI3 and TRI6 meshes
486  // have twice as many elements...
487  switch (type)
488  {
489  case INVALID_ELEM:
490  case QUAD4:
491  case QUAD8:
492  case QUAD9:
493  {
494  mesh.reserve_elem (nx*ny);
495  break;
496  }
497 
498  case TRI3:
499  case TRI6:
500  {
501  mesh.reserve_elem (2*nx*ny);
502  break;
503  }
504 
505  default:
506  {
507  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
508  libmesh_error();
509  }
510  }
511 
512 
513 
514  // Reserve nodes. The quadratic element types
515  // need to reserve more nodes than the linear types.
516  switch (type)
517  {
518  case INVALID_ELEM:
519  case QUAD4:
520  case TRI3:
521  {
522  mesh.reserve_nodes( (nx+1)*(ny+1) );
523  break;
524  }
525 
526  case QUAD8:
527  case QUAD9:
528  case TRI6:
529  {
530  mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
531  break;
532  }
533 
534 
535  default:
536  {
537  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
538  libmesh_error();
539  }
540  }
541 
542 
543 
544  // Build the nodes. Depends on whether you are using a linear
545  // or quadratic element, and whether you are using a uniform
546  // grid or the Gauss-Lobatto grid points.
547  unsigned int node_id = 0;
548  switch (type)
549  {
550  case INVALID_ELEM:
551  case QUAD4:
552  case TRI3:
553  {
554  for (unsigned int j=0; j<=ny; j++)
555  for (unsigned int i=0; i<=nx; i++)
556  {
557  if (gauss_lobatto_grid)
558  {
559  mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
560  0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
561  0.), node_id++);
562  }
563 
564  else
565  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
566  static_cast<Real>(j)/static_cast<Real>(ny),
567  0.), node_id++);
568  }
569 
570  break;
571  }
572 
573  case QUAD8:
574  case QUAD9:
575  case TRI6:
576  {
577  for (unsigned int j=0; j<=(2*ny); j++)
578  for (unsigned int i=0; i<=(2*nx); i++)
579  {
580  if (gauss_lobatto_grid)
581  {
582  // The x,y locations of the point.
583  Real x=0., y=0.;
584 
585  // Shortcut quantities (do not depend on i,j)
586  const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
587  const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
588 
589  // Shortcut quantities (depend on i,j)
590  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
591  const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
592 
593  // If i is even, compute a normal Gauss-Lobatto point
594  if (i%2 == 0)
595  x = 0.5*(1.0 - c);
596 
597  // Otherwise, it is the average of the previous and next points
598  else
599  x = 0.5*(1.0 - a*c);
600 
601  // If j is even, compute a normal Gauss-Lobatto point
602  if (j%2 == 0)
603  y = 0.5*(1.0 - d);
604 
605  // Otherwise, it is the average of the previous and next points
606  else
607  y = 0.5*(1.0 - b*d);
608 
609 
610  mesh.add_point (Point(x,y,0.), node_id++);
611  }
612 
613 
614  else
615  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
616  static_cast<Real>(j)/static_cast<Real>(2*ny),
617  0), node_id++);
618  }
619 
620  break;
621  }
622 
623 
624  default:
625  {
626  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
627  libmesh_error();
628  }
629  }
630 
631 
632 
633 
634 
635 
636  // Build the elements. Each one is a bit different.
637  switch (type)
638  {
639 
640  case INVALID_ELEM:
641  case QUAD4:
642  {
643  for (unsigned int j=0; j<ny; j++)
644  for (unsigned int i=0; i<nx; i++)
645  {
646  Elem* elem = mesh.add_elem(new Quad4);
647 
648  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
649  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) );
650  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
651  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+1) );
652 
653  if (j == 0)
654  mesh.boundary_info->add_side(elem, 0, 0);
655 
656  if (j == (ny-1))
657  mesh.boundary_info->add_side(elem, 2, 2);
658 
659  if (i == 0)
660  mesh.boundary_info->add_side(elem, 3, 3);
661 
662  if (i == (nx-1))
663  mesh.boundary_info->add_side(elem, 1, 1);
664  }
665  break;
666  }
667 
668 
669  case TRI3:
670  {
671  for (unsigned int j=0; j<ny; j++)
672  for (unsigned int i=0; i<nx; i++)
673  {
674  Elem* elem = NULL;
675 
676  // Add first Tri3
677  elem = mesh.add_elem(new Tri3);
678 
679  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
680  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) );
681  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
682 
683  if (j == 0)
684  mesh.boundary_info->add_side(elem, 0, 0);
685 
686  if (i == (nx-1))
687  mesh.boundary_info->add_side(elem, 1, 1);
688 
689  // Add second Tri3
690  elem = mesh.add_elem(new Tri3);
691 
692  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
693  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j+1));
694  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+1) );
695 
696  if (j == (ny-1))
697  mesh.boundary_info->add_side(elem, 1, 2);
698 
699  if (i == 0)
700  mesh.boundary_info->add_side(elem, 2, 3);
701  }
702  break;
703  }
704 
705 
706 
707  case QUAD8:
708  case QUAD9:
709  {
710  for (unsigned int j=0; j<(2*ny); j += 2)
711  for (unsigned int i=0; i<(2*nx); i += 2)
712  {
713  Elem* elem = (type == QUAD8) ?
714  mesh.add_elem(new Quad8) :
715  mesh.add_elem(new Quad9);
716 
717 
718  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
719  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j) );
720  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
721  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+2) );
722  elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j) );
723  elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+2,j+1));
724  elem->set_node(6) = mesh.node_ptr(idx(type,nx,i+1,j+2));
725  elem->set_node(7) = mesh.node_ptr(idx(type,nx,i,j+1) );
726  if (type == QUAD9)
727  elem->set_node(8) = mesh.node_ptr(idx(type,nx,i+1,j+1));
728 
729 
730  if (j == 0)
731  mesh.boundary_info->add_side(elem, 0, 0);
732 
733  if (j == 2*(ny-1))
734  mesh.boundary_info->add_side(elem, 2, 2);
735 
736  if (i == 0)
737  mesh.boundary_info->add_side(elem, 3, 3);
738 
739  if (i == 2*(nx-1))
740  mesh.boundary_info->add_side(elem, 1, 1);
741  }
742  break;
743  }
744 
745 
746  case TRI6:
747  {
748  for (unsigned int j=0; j<(2*ny); j += 2)
749  for (unsigned int i=0; i<(2*nx); i += 2)
750  {
751  Elem* elem = NULL;
752 
753  // Add first Tri6
754  elem = mesh.add_elem(new Tri6);
755 
756  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
757  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j) );
758  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
759  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j) );
760  elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+2,j+1));
761  elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+1,j+1));
762 
763  if (j == 0)
764  mesh.boundary_info->add_side(elem, 0, 0);
765 
766  if (i == 2*(nx-1))
767  mesh.boundary_info->add_side(elem, 1, 1);
768 
769  // Add second Tri6
770  elem = mesh.add_elem(new Tri6);
771 
772  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
773  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j+2));
774  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+2) );
775  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j+1));
776  elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j+2));
777  elem->set_node(5) = mesh.node_ptr(idx(type,nx,i,j+1) );
778 
779  if (j == 2*(ny-1))
780  mesh.boundary_info->add_side(elem, 1, 2);
781 
782  if (i == 0)
783  mesh.boundary_info->add_side(elem, 2, 3);
784 
785  }
786  break;
787  };
788 
789 
790  default:
791  {
792  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
793  libmesh_error();
794  }
795  }
796 
797 
798 
799 
800  // Scale the nodal positions
801  for (unsigned int p=0; p<mesh.n_nodes(); p++)
802  {
803  mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
804  mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
805  }
806 
807  // Add sideset names to boundary info
808  mesh.boundary_info->sideset_name(0) = "bottom";
809  mesh.boundary_info->sideset_name(1) = "right";
810  mesh.boundary_info->sideset_name(2) = "top";
811  mesh.boundary_info->sideset_name(3) = "left";
812 
813  // Add nodeset names to boundary info
814  mesh.boundary_info->nodeset_name(0) = "bottom";
815  mesh.boundary_info->nodeset_name(1) = "right";
816  mesh.boundary_info->nodeset_name(2) = "top";
817  mesh.boundary_info->nodeset_name(3) = "left";
818 
819  break;
820  }
821 
822 
823 
824 
825 
826 
827 
828 
829 
830 
831 
832  //---------------------------------------------------------------------
833  // Build a 3D mesh using hexes, tets, prisms, or pyramids.
834  case 3:
835  {
836  libmesh_assert_not_equal_to (nx, 0);
837  libmesh_assert_not_equal_to (ny, 0);
838  libmesh_assert_not_equal_to (nz, 0);
839  libmesh_assert_less (xmin, xmax);
840  libmesh_assert_less (ymin, ymax);
841  libmesh_assert_less (zmin, zmax);
842 
843 
844  // Reserve elements. Meshes with prismatic elements require
845  // twice as many elements.
846  switch (type)
847  {
848  case INVALID_ELEM:
849  case HEX8:
850  case HEX20:
851  case HEX27:
852  case TET4: // TET4's are created from an initial HEX27 discretization
853  case TET10: // TET10's are created from an initial HEX27 discretization
854  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
855  case PYRAMID14:
856  {
857  mesh.reserve_elem(nx*ny*nz);
858  break;
859  }
860 
861  case PRISM6:
862  case PRISM15:
863  case PRISM18:
864  {
865  mesh.reserve_elem(2*nx*ny*nz);
866  break;
867  }
868 
869  default:
870  {
871  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
872  libmesh_error();
873  }
874  }
875 
876 
877 
878 
879 
880  // Reserve nodes. Quadratic elements need twice as many nodes as linear elements.
881  switch (type)
882  {
883  case INVALID_ELEM:
884  case HEX8:
885  case PRISM6:
886  {
887  mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
888  break;
889  }
890 
891  case HEX20:
892  case HEX27:
893  case TET4: // TET4's are created from an initial HEX27 discretization
894  case TET10: // TET10's are created from an initial HEX27 discretization
895  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
896  case PYRAMID14:
897  case PRISM15:
898  case PRISM18:
899  {
900  // FYI: The resulting TET4 mesh will have exactly
901  // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
902  // nodes once the additional mid-edge nodes for the HEX27 discretization
903  // have been deleted.
904  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
905  break;
906  }
907 
908  default:
909  {
910  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
911  libmesh_error();
912  }
913  }
914 
915 
916 
917 
918  // Build the nodes.
919  unsigned int node_id = 0;
920  switch (type)
921  {
922  case INVALID_ELEM:
923  case HEX8:
924  case PRISM6:
925  {
926  for (unsigned int k=0; k<=nz; k++)
927  for (unsigned int j=0; j<=ny; j++)
928  for (unsigned int i=0; i<=nx; i++)
929  {
930  if (gauss_lobatto_grid)
931  {
932  mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
933  0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
934  0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(k)/static_cast<Real>(nz)))), node_id++);
935  }
936 
937  else
938  mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(nx),
939  static_cast<Real>(j)/static_cast<Real>(ny),
940  static_cast<Real>(k)/static_cast<Real>(nz)), node_id++);
941  }
942  break;
943  }
944 
945  case HEX20:
946  case HEX27:
947  case TET4: // TET4's are created from an initial HEX27 discretization
948  case TET10: // TET10's are created from an initial HEX27 discretization
949  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
950  case PYRAMID14:
951  case PRISM15:
952  case PRISM18:
953  {
954  for (unsigned int k=0; k<=(2*nz); k++)
955  for (unsigned int j=0; j<=(2*ny); j++)
956  for (unsigned int i=0; i<=(2*nx); i++)
957  {
958  if (gauss_lobatto_grid)
959  {
960  // The x,y locations of the point.
961  Real x=0., y=0., z=0.;
962 
963  // Shortcut quantities (do not depend on i,j)
964  const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
965  const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
966 
967  // Shortcut quantities (depend on i,j)
968  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
969  const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
970 
971  // Additional shortcut quantities (for 3D)
972  const Real e = std::cos( libMesh::pi / static_cast<Real>(2*nz) );
973  const Real f = std::cos( libMesh::pi*k / static_cast<Real>(2*nz) );
974 
975  // If i is even, compute a normal Gauss-Lobatto point
976  if (i%2 == 0)
977  x = 0.5*(1.0 - c);
978 
979  // Otherwise, it is the average of the previous and next points
980  else
981  x = 0.5*(1.0 - a*c);
982 
983  // If j is even, compute a normal Gauss-Lobatto point
984  if (j%2 == 0)
985  y = 0.5*(1.0 - d);
986 
987  // Otherwise, it is the average of the previous and next points
988  else
989  y = 0.5*(1.0 - b*d);
990 
991  // If k is even, compute a normal Gauss-Lobatto point
992  if (k%2 == 0)
993  z = 0.5*(1.0 - f);
994 
995  // Otherwise, it is the average of the previous and next points
996  else
997  z = 0.5*(1.0 - e*f);
998 
999 
1000  mesh.add_point (Point(x,y,z), node_id++);
1001  }
1002 
1003  else
1004  mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
1005  static_cast<Real>(j)/static_cast<Real>(2*ny),
1006  static_cast<Real>(k)/static_cast<Real>(2*nz)), node_id++);
1007  }
1008  break;
1009  }
1010 
1011 
1012  default:
1013  {
1014  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
1015  libmesh_error();
1016  }
1017  }
1018 
1019 
1020 
1021 
1022  // Build the elements.
1023  switch (type)
1024  {
1025  case INVALID_ELEM:
1026  case HEX8:
1027  {
1028  for (unsigned int k=0; k<nz; k++)
1029  for (unsigned int j=0; j<ny; j++)
1030  for (unsigned int i=0; i<nx; i++)
1031  {
1032  Elem* elem = mesh.add_elem(new Hex8);
1033 
1034  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) );
1035  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
1036  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1037  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
1038  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) );
1039  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
1040  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1041  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
1042 
1043  if (k == 0)
1044  mesh.boundary_info->add_side(elem, 0, 0);
1045 
1046  if (k == (nz-1))
1047  mesh.boundary_info->add_side(elem, 5, 5);
1048 
1049  if (j == 0)
1050  mesh.boundary_info->add_side(elem, 1, 1);
1051 
1052  if (j == (ny-1))
1053  mesh.boundary_info->add_side(elem, 3, 3);
1054 
1055  if (i == 0)
1056  mesh.boundary_info->add_side(elem, 4, 4);
1057 
1058  if (i == (nx-1))
1059  mesh.boundary_info->add_side(elem, 2, 2);
1060  }
1061  break;
1062  }
1063 
1064 
1065 
1066 
1067  case PRISM6:
1068  {
1069  for (unsigned int k=0; k<nz; k++)
1070  for (unsigned int j=0; j<ny; j++)
1071  for (unsigned int i=0; i<nx; i++)
1072  {
1073  // First Prism
1074  Elem* elem = NULL;
1075  elem = mesh.add_elem(new Prism6);
1076 
1077  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) );
1078  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
1079  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
1080  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) );
1081  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
1082  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
1083 
1084  // Add sides for first prism to boundary info object
1085  if (i==0)
1086  mesh.boundary_info->add_side(elem, 3, 4);
1087 
1088  if (j==0)
1089  mesh.boundary_info->add_side(elem, 1, 1);
1090 
1091  if (k==0)
1092  mesh.boundary_info->add_side(elem, 0, 0);
1093 
1094  if (k == (nz-1))
1095  mesh.boundary_info->add_side(elem, 4, 5);
1096 
1097  // Second Prism
1098  elem = mesh.add_elem(new Prism6);
1099 
1100  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
1101  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1102  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
1103  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
1104  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1105  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
1106 
1107  // Add sides for second prism to boundary info object
1108  if (i == (nx-1))
1109  mesh.boundary_info->add_side(elem, 1, 2);
1110 
1111  if (j == (ny-1))
1112  mesh.boundary_info->add_side(elem, 2, 3);
1113 
1114  if (k==0)
1115  mesh.boundary_info->add_side(elem, 0, 0);
1116 
1117  if (k == (nz-1))
1118  mesh.boundary_info->add_side(elem, 4, 5);
1119  }
1120  break;
1121  }
1122 
1123 
1124 
1125 
1126 
1127 
1128  case HEX20:
1129  case HEX27:
1130  case TET4: // TET4's are created from an initial HEX27 discretization
1131  case TET10: // TET10's are created from an initial HEX27 discretization
1132  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
1133  case PYRAMID14:
1134  {
1135  for (unsigned int k=0; k<(2*nz); k += 2)
1136  for (unsigned int j=0; j<(2*ny); j += 2)
1137  for (unsigned int i=0; i<(2*nx); i += 2)
1138  {
1139  Elem* elem = (type == HEX20) ?
1140  mesh.add_elem(new Hex20) :
1141  mesh.add_elem(new Hex27);
1142 
1143  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) );
1144  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) );
1145  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) );
1146  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) );
1147  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2));
1148  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2));
1149  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2));
1150  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2));
1151  elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) );
1152  elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) );
1153  elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) );
1154  elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) );
1155  elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1));
1156  elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1));
1157  elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
1158  elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1));
1159  elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2));
1160  elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
1161  elem->set_node(18) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
1162  elem->set_node(19) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2));
1163  if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == PYRAMID5) || (type == PYRAMID14))
1164  {
1165  elem->set_node(20) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1166  elem->set_node(21) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1));
1167  elem->set_node(22) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
1168  elem->set_node(23) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
1169  elem->set_node(24) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1));
1170  elem->set_node(25) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
1171  elem->set_node(26) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1172  }
1173 
1174 
1175  if (k == 0)
1176  mesh.boundary_info->add_side(elem, 0, 0);
1177 
1178  if (k == 2*(nz-1))
1179  mesh.boundary_info->add_side(elem, 5, 5);
1180 
1181  if (j == 0)
1182  mesh.boundary_info->add_side(elem, 1, 1);
1183 
1184  if (j == 2*(ny-1))
1185  mesh.boundary_info->add_side(elem, 3, 3);
1186 
1187  if (i == 0)
1188  mesh.boundary_info->add_side(elem, 4, 4);
1189 
1190  if (i == 2*(nx-1))
1191  mesh.boundary_info->add_side(elem, 2, 2);
1192  }
1193  break;
1194  }
1195 
1196 
1197 
1198 
1199  case PRISM15:
1200  case PRISM18:
1201  {
1202  for (unsigned int k=0; k<(2*nz); k += 2)
1203  for (unsigned int j=0; j<(2*ny); j += 2)
1204  for (unsigned int i=0; i<(2*nx); i += 2)
1205  {
1206  // First Prism
1207  Elem* elem = NULL;
1208  elem = ((type == PRISM15) ?
1209  mesh.add_elem(new Prism15) :
1210  mesh.add_elem(new Prism18));
1211 
1212  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) );
1213  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) );
1214  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) );
1215  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2));
1216  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2));
1217  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2));
1218  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) );
1219  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1220  elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) );
1221  elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1));
1222  elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1));
1223  elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1));
1224  elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2));
1225  elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
1226  elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2));
1227  if (type == PRISM18)
1228  {
1229  elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1));
1230  elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1231  elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1));
1232  }
1233 
1234  // Add sides for first prism to boundary info object
1235  if (i==0)
1236  mesh.boundary_info->add_side(elem, 3, 4);
1237 
1238  if (j==0)
1239  mesh.boundary_info->add_side(elem, 1, 1);
1240 
1241  if (k==0)
1242  mesh.boundary_info->add_side(elem, 0, 0);
1243 
1244  if (k == 2*(nz-1))
1245  mesh.boundary_info->add_side(elem, 4, 5);
1246 
1247 
1248  // Second Prism
1249  elem = ((type == PRISM15) ?
1250  mesh.add_elem(new Prism15) :
1251  mesh.add_elem(new Prism18));
1252 
1253  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k) );
1254  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) );
1255  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k) );
1256  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) );
1257  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) );
1258  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) );
1259  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) );
1260  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) );
1261  elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1262  elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) );
1263  elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
1264  elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) );
1265  elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
1266  elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
1267  elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
1268  if (type == PRISM18)
1269  {
1270  elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
1271  elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
1272  elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1273  }
1274 
1275  // Add sides for second prism to boundary info object
1276  if (i == 2*(nx-1))
1277  mesh.boundary_info->add_side(elem, 1, 2);
1278 
1279  if (j == 2*(ny-1))
1280  mesh.boundary_info->add_side(elem, 2, 3);
1281 
1282  if (k==0)
1283  mesh.boundary_info->add_side(elem, 0, 0);
1284 
1285  if (k == 2*(nz-1))
1286  mesh.boundary_info->add_side(elem, 4, 5);
1287 
1288  }
1289  break;
1290  }
1291 
1292 
1293 
1294 
1295 
1296  default:
1297  {
1298  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
1299  libmesh_error();
1300  }
1301  }
1302 
1303 
1304 
1305 
1306  //.......................................
1307  // Scale the nodal positions
1308  for (unsigned int p=0; p<mesh.n_nodes(); p++)
1309  {
1310  mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
1311  mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
1312  mesh.node(p)(2) = (mesh.node(p)(2))*(zmax-zmin) + zmin;
1313  }
1314 
1315 
1316 
1317 
1318  // Additional work for tets and pyramids: we take the existing
1319  // HEX27 discretization and split each element into 24
1320  // sub-tets or 6 sub-pyramids.
1321  //
1322  // 24 isn't the minimum-possible number of tets, but it
1323  // obviates any concerns about the edge orientations between
1324  // the various elements.
1325  if ((type == TET4) ||
1326  (type == TET10) ||
1327  (type == PYRAMID5) ||
1328  (type == PYRAMID14))
1329  {
1330  // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
1331  std::vector<Elem*> new_elements;
1332 
1333  if ((type == TET4) || (type == TET10))
1334  new_elements.reserve(24*mesh.n_elem());
1335  else
1336  new_elements.reserve(6*mesh.n_elem());
1337 
1338  // Create tetrahedra or pyramids
1339  {
1341  const MeshBase::element_iterator end_el = mesh.elements_end();
1342 
1343  for ( ; el != end_el; ++el)
1344  {
1345  // Get a pointer to the HEX27 element.
1346  Elem* base_hex = *el;
1347 
1348  // Get a pointer to the node located at the HEX27 centroid
1349  Node* apex_node = base_hex->get_node(26);
1350 
1351  for (unsigned int s=0; s<base_hex->n_sides(); ++s)
1352  {
1353  // Get the boundary ID for this side
1354  boundary_id_type b_id = mesh.boundary_info->boundary_id(*el, s);
1355 
1356  // Need to build the full-ordered side!
1357  AutoPtr<Elem> side = base_hex->build_side(s);
1358 
1359  if ((type == TET4) || (type == TET10))
1360  {
1361  // Build 4 sub-tets per side
1362  for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
1363  {
1364  new_elements.push_back( new Tet4 );
1365  Elem* sub_elem = new_elements.back();
1366  sub_elem->set_node(0) = side->get_node(sub_tet);
1367  sub_elem->set_node(1) = side->get_node(8); // centroid of the face
1368  sub_elem->set_node(2) = side->get_node(sub_tet==3 ? 0 : sub_tet+1 ); // wrap-around
1369  sub_elem->set_node(3) = apex_node; // apex node always used!
1370 
1371  // If the original hex was a boundary hex, add the new sub_tet's side
1372  // 0 with the same b_id. Note: the tets are all aligned so that their
1373  // side 0 is on the boundary.
1374  if (b_id != BoundaryInfo::invalid_id)
1375  mesh.boundary_info->add_side(sub_elem, 0, b_id);
1376  }
1377  } // end if ((type == TET4) || (type == TET10))
1378 
1379  else // type==PYRAMID5 || type==PYRAMID14
1380  {
1381  // Build 1 sub-pyramid per side.
1382  new_elements.push_back(new Pyramid5);
1383  Elem* sub_elem = new_elements.back();
1384 
1385  // Set the base. Note that since the apex is *inside* the base_hex,
1386  // and the pyramid uses a counter-clockwise base numbering, we need to
1387  // reverse the [1] and [3] node indices.
1388  sub_elem->set_node(0) = side->get_node(0);
1389  sub_elem->set_node(1) = side->get_node(3);
1390  sub_elem->set_node(2) = side->get_node(2);
1391  sub_elem->set_node(3) = side->get_node(1);
1392 
1393  // Set the apex
1394  sub_elem->set_node(4) = apex_node;
1395 
1396  // If the original hex was a boundary hex, add the new sub_pyr's side
1397  // 4 (the square base) with the same b_id.
1398  if (b_id != BoundaryInfo::invalid_id)
1399  mesh.boundary_info->add_side(sub_elem, 4, b_id);
1400  } // end else type==PYRAMID5 || type==PYRAMID14
1401  }
1402  }
1403  }
1404 
1405 
1406  // Delete the original HEX27 elements from the mesh, and the boundary info structure.
1407  {
1409  const MeshBase::element_iterator end_el = mesh.elements_end();
1410 
1411  for ( ; el != end_el; ++el)
1412  {
1413  mesh.boundary_info->remove(*el); // Safe even if *el has no boundary info.
1414  mesh.delete_elem(*el);
1415  }
1416  }
1417 
1418  // Add the new elements
1419  for (unsigned int i=0; i<new_elements.size(); ++i)
1420  mesh.add_elem(new_elements[i]);
1421 
1422  } // end if (type == TET4,TET10,PYRAMID5,PYRAMID14
1423 
1424 
1425  // Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's
1426  if ((type == TET10) || (type == PYRAMID14))
1427  {
1428  mesh.all_second_order();
1429  }
1430 
1431  // Add sideset names to boundary info (Z axis out of the screen)
1432  mesh.boundary_info->sideset_name(0) = "back";
1433  mesh.boundary_info->sideset_name(1) = "bottom";
1434  mesh.boundary_info->sideset_name(2) = "right";
1435  mesh.boundary_info->sideset_name(3) = "top";
1436  mesh.boundary_info->sideset_name(4) = "left";
1437  mesh.boundary_info->sideset_name(5) = "front";
1438 
1439  // Add nodeset names to boundary info
1440  mesh.boundary_info->nodeset_name(0) = "back";
1441  mesh.boundary_info->nodeset_name(1) = "bottom";
1442  mesh.boundary_info->nodeset_name(2) = "right";
1443  mesh.boundary_info->nodeset_name(3) = "top";
1444  mesh.boundary_info->nodeset_name(4) = "left";
1445  mesh.boundary_info->nodeset_name(5) = "front";
1446 
1447  break;
1448  } // end case dim==3
1449 
1450  default:
1451  {
1452  libmesh_error();
1453  }
1454  }
1455 
1456  STOP_LOG("build_cube()", "MeshTools::Generation");
1457 
1458 
1459 
1460  // Done building the mesh. Now prepare it for use.
1461  mesh.prepare_for_use (/*skip_renumber =*/ false);
1462 }
1463 
1464 
1465 
1467  const ElemType type,
1468  const bool gauss_lobatto_grid)
1469 {
1470  // This method only makes sense in 0D!
1471  // But we now just turn a non-0D mesh into a 0D mesh
1472  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1473 
1474  build_cube(mesh,
1475  0, 0, 0,
1476  0., 0.,
1477  0., 0.,
1478  0., 0.,
1479  type,
1480  gauss_lobatto_grid);
1481 }
1482 
1483 
1485  const unsigned int nx,
1486  const Real xmin, const Real xmax,
1487  const ElemType type,
1488  const bool gauss_lobatto_grid)
1489 {
1490  // This method only makes sense in 1D!
1491  // But we now just turn a non-1D mesh into a 1D mesh
1492  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1493 
1494  build_cube(mesh,
1495  nx, 0, 0,
1496  xmin, xmax,
1497  0., 0.,
1498  0., 0.,
1499  type,
1500  gauss_lobatto_grid);
1501 }
1502 
1503 
1504 
1506  const unsigned int nx,
1507  const unsigned int ny,
1508  const Real xmin, const Real xmax,
1509  const Real ymin, const Real ymax,
1510  const ElemType type,
1511  const bool gauss_lobatto_grid)
1512 {
1513  // This method only makes sense in 2D!
1514  // But we now just turn a non-2D mesh into a 2D mesh
1515  //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1516 
1517  // Call the build_cube() member to actually do the work for us.
1518  build_cube (mesh,
1519  nx, ny, 0,
1520  xmin, xmax,
1521  ymin, ymax,
1522  0., 0.,
1523  type,
1524  gauss_lobatto_grid);
1525 }
1526 
1527 
1528 
1529 
1530 
1531 
1532 
1533 
1534 
1535 #ifndef LIBMESH_ENABLE_AMR
1537  const Real,
1538  const unsigned int,
1539  const ElemType,
1540  const unsigned int,
1541  const bool)
1542 {
1543  libMesh::out << "Building a circle/sphere only works with AMR." << std::endl;
1544  libmesh_error();
1545 }
1546 
1547 #else
1548 
1550  const Real rad,
1551  const unsigned int nr,
1552  const ElemType type,
1553  const unsigned int n_smooth,
1554  const bool flat)
1555 {
1556  libmesh_assert_greater (rad, 0.);
1557  //libmesh_assert_greater (nr, 0); // must refine at least once otherwise will end up with a square/cube
1558 
1559  START_LOG("build_sphere()", "MeshTools::Generation");
1560 
1561  // Clear the mesh and start from scratch
1562  mesh.clear();
1563 
1564  // Sphere is centered at origin by default
1565  const Point cent;
1566 
1567  const Sphere sphere (cent, rad);
1568 
1569  switch (mesh.mesh_dimension())
1570  {
1571  //-----------------------------------------------------------------
1572  // Build a line in one dimension
1573  case 1:
1574  {
1575  build_line (mesh, 3, -rad, rad, type);
1576  }
1577 
1578 
1579 
1580 
1581  //-----------------------------------------------------------------
1582  // Build a circle or hollow sphere in two dimensions
1583  case 2:
1584  {
1585  // For ParallelMesh, if we don't specify node IDs the Mesh
1586  // will try to pick an appropriate (unique) one for us. But
1587  // since we are adding these nodes on all processors, we want
1588  // to be sure they have consistent IDs across all processors.
1589  unsigned node_id = 0;
1590 
1591  if (flat)
1592  {
1593  const Real sqrt_2 = std::sqrt(2.);
1594  const Real rad_2 = .25*rad;
1595  const Real rad_sqrt_2 = rad/sqrt_2;
1596 
1597  // (Temporary) convenient storage for node pointers
1598  std::vector<Node*> nodes(8);
1599 
1600  // Point 0
1601  nodes[0] = mesh.add_point (Point(-rad_2,-rad_2, 0.), node_id++);
1602 
1603  // Point 1
1604  nodes[1] = mesh.add_point (Point( rad_2,-rad_2, 0.), node_id++);
1605 
1606  // Point 2
1607  nodes[2] = mesh.add_point (Point( rad_2, rad_2, 0.), node_id++);
1608 
1609  // Point 3
1610  nodes[3] = mesh.add_point (Point(-rad_2, rad_2, 0.), node_id++);
1611 
1612  // Point 4
1613  nodes[4] = mesh.add_point (Point(-rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
1614 
1615  // Point 5
1616  nodes[5] = mesh.add_point (Point( rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
1617 
1618  // Point 6
1619  nodes[6] = mesh.add_point (Point( rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
1620 
1621  // Point 7
1622  nodes[7] = mesh.add_point (Point(-rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
1623 
1624  // Build the elements & set node pointers
1625 
1626  // Element 0
1627  {
1628  Elem* elem0 = mesh.add_elem (new Quad4);
1629  elem0->set_node(0) = nodes[0];
1630  elem0->set_node(1) = nodes[1];
1631  elem0->set_node(2) = nodes[2];
1632  elem0->set_node(3) = nodes[3];
1633  }
1634 
1635  // Element 1
1636  {
1637  Elem* elem1 = mesh.add_elem (new Quad4);
1638  elem1->set_node(0) = nodes[4];
1639  elem1->set_node(1) = nodes[0];
1640  elem1->set_node(2) = nodes[3];
1641  elem1->set_node(3) = nodes[7];
1642  }
1643 
1644  // Element 2
1645  {
1646  Elem* elem2 = mesh.add_elem (new Quad4);
1647  elem2->set_node(0) = nodes[4];
1648  elem2->set_node(1) = nodes[5];
1649  elem2->set_node(2) = nodes[1];
1650  elem2->set_node(3) = nodes[0];
1651  }
1652 
1653  // Element 3
1654  {
1655  Elem* elem3 = mesh.add_elem (new Quad4);
1656  elem3->set_node(0) = nodes[1];
1657  elem3->set_node(1) = nodes[5];
1658  elem3->set_node(2) = nodes[6];
1659  elem3->set_node(3) = nodes[2];
1660  }
1661 
1662  // Element 4
1663  {
1664  Elem* elem4 = mesh.add_elem (new Quad4);
1665  elem4->set_node(0) = nodes[3];
1666  elem4->set_node(1) = nodes[2];
1667  elem4->set_node(2) = nodes[6];
1668  elem4->set_node(3) = nodes[7];
1669  }
1670 
1671  }
1672  else
1673  {
1674  // Create the 12 vertices of a regular unit icosahedron
1675  Real t = 0.5 * (1 + std::sqrt(5.0));
1676  Real s = rad / std::sqrt(1 + t*t);
1677  t *= s;
1678 
1679  mesh.add_point (Point(-s, t, 0), node_id++);
1680  mesh.add_point (Point( s, t, 0), node_id++);
1681  mesh.add_point (Point(-s, -t, 0), node_id++);
1682  mesh.add_point (Point( s, -t, 0), node_id++);
1683 
1684  mesh.add_point (Point( 0, -s, t), node_id++);
1685  mesh.add_point (Point( 0, s, t), node_id++);
1686  mesh.add_point (Point( 0, -s, -t), node_id++);
1687  mesh.add_point (Point( 0, s, -t), node_id++);
1688 
1689  mesh.add_point (Point( t, 0, -s), node_id++);
1690  mesh.add_point (Point( t, 0, s), node_id++);
1691  mesh.add_point (Point(-t, 0, -s), node_id++);
1692  mesh.add_point (Point(-t, 0, s), node_id++);
1693 
1694  // Create the 20 triangles of the icosahedron
1695  static const unsigned int idx1 [6] = {11, 5, 1, 7, 10, 11};
1696  static const unsigned int idx2 [6] = {9, 4, 2, 6, 8, 9};
1697  static const unsigned int idx3 [6] = {1, 5, 11, 10, 7, 1};
1698 
1699  for (unsigned int i = 0; i < 5; ++i)
1700  {
1701  // 5 elems around point 0
1702  Elem* new_elem = mesh.add_elem (new Tri3);
1703  new_elem->set_node(0) = mesh.node_ptr(0);
1704  new_elem->set_node(1) = mesh.node_ptr(idx1[i]);
1705  new_elem->set_node(2) = mesh.node_ptr(idx1[i+1]);
1706 
1707  // 5 adjacent elems
1708  new_elem = mesh.add_elem (new Tri3);
1709  new_elem->set_node(0) = mesh.node_ptr(idx3[i]);
1710  new_elem->set_node(1) = mesh.node_ptr(idx3[i+1]);
1711  new_elem->set_node(2) = mesh.node_ptr(idx2[i]);
1712 
1713  // 5 elems around point 3
1714  new_elem = mesh.add_elem (new Tri3);
1715  new_elem->set_node(0) = mesh.node_ptr(3);
1716  new_elem->set_node(1) = mesh.node_ptr(idx2[i]);
1717  new_elem->set_node(2) = mesh.node_ptr(idx2[i+1]);
1718 
1719  // 5 adjacent elems
1720  new_elem = mesh.add_elem (new Tri3);
1721  new_elem->set_node(0) = mesh.node_ptr(idx2[i+1]);
1722  new_elem->set_node(1) = mesh.node_ptr(idx2[i]);
1723  new_elem->set_node(2) = mesh.node_ptr(idx3[i+1]);
1724  }
1725  }
1726 
1727  break;
1728  } // end case 2
1729 
1730 
1731 
1732 
1733 
1734  //-----------------------------------------------------------------
1735  // Build a sphere in three dimensions
1736  case 3:
1737  {
1738  // (Currently) supported types
1739  if (!((type == HEX8) || (type == HEX27)))
1740  {
1741  // FIXME: We'd need an all_tet() routine (which could also be used by
1742  // build_square()) to do Tets, or Prisms for that matter.
1743  libMesh::err << "Error: Only HEX8/27 currently supported."
1744  << std::endl;
1745  libmesh_error();
1746  }
1747 
1748 
1749  // 3D analog of 2D initial grid:
1750  const Real
1751  r_small = 0.25*rad, // 0.25 *radius
1752  r_med = (0.125*std::sqrt(2.)+0.5)*rad; // .67677*radius
1753 
1754  // (Temporary) convenient storage for node pointers
1755  std::vector<Node*> nodes(16);
1756 
1757  // For ParallelMesh, if we don't specify node IDs the Mesh
1758  // will try to pick an appropriate (unique) one for us. But
1759  // since we are adding these nodes on all processors, we want
1760  // to be sure they have consistent IDs across all processors.
1761  unsigned node_id = 0;
1762 
1763  // Points 0-7 are the initial HEX8
1764  nodes[0] = mesh.add_point (Point(-r_small,-r_small, -r_small), node_id++);
1765  nodes[1] = mesh.add_point (Point( r_small,-r_small, -r_small), node_id++);
1766  nodes[2] = mesh.add_point (Point( r_small, r_small, -r_small), node_id++);
1767  nodes[3] = mesh.add_point (Point(-r_small, r_small, -r_small), node_id++);
1768  nodes[4] = mesh.add_point (Point(-r_small,-r_small, r_small), node_id++);
1769  nodes[5] = mesh.add_point (Point( r_small,-r_small, r_small), node_id++);
1770  nodes[6] = mesh.add_point (Point( r_small, r_small, r_small), node_id++);
1771  nodes[7] = mesh.add_point (Point(-r_small, r_small, r_small), node_id++);
1772 
1773  // Points 8-15 are for the outer hexes, we number them in the same way
1774  nodes[8] = mesh.add_point (Point(-r_med,-r_med, -r_med), node_id++);
1775  nodes[9] = mesh.add_point (Point( r_med,-r_med, -r_med), node_id++);
1776  nodes[10] = mesh.add_point (Point( r_med, r_med, -r_med), node_id++);
1777  nodes[11] = mesh.add_point (Point(-r_med, r_med, -r_med), node_id++);
1778  nodes[12] = mesh.add_point (Point(-r_med,-r_med, r_med), node_id++);
1779  nodes[13] = mesh.add_point (Point( r_med,-r_med, r_med), node_id++);
1780  nodes[14] = mesh.add_point (Point( r_med, r_med, r_med), node_id++);
1781  nodes[15] = mesh.add_point (Point(-r_med, r_med, r_med), node_id++);
1782 
1783  // Now create the elements and add them to the mesh
1784  // Element 0 - center element
1785  {
1786  Elem* elem0 = mesh.add_elem (new Hex8);
1787  elem0->set_node(0) = nodes[0];
1788  elem0->set_node(1) = nodes[1];
1789  elem0->set_node(2) = nodes[2];
1790  elem0->set_node(3) = nodes[3];
1791  elem0->set_node(4) = nodes[4];
1792  elem0->set_node(5) = nodes[5];
1793  elem0->set_node(6) = nodes[6];
1794  elem0->set_node(7) = nodes[7];
1795  }
1796 
1797  // Element 1 - "bottom"
1798  {
1799  Elem* elem1 = mesh.add_elem (new Hex8);
1800  elem1->set_node(0) = nodes[8];
1801  elem1->set_node(1) = nodes[9];
1802  elem1->set_node(2) = nodes[10];
1803  elem1->set_node(3) = nodes[11];
1804  elem1->set_node(4) = nodes[0];
1805  elem1->set_node(5) = nodes[1];
1806  elem1->set_node(6) = nodes[2];
1807  elem1->set_node(7) = nodes[3];
1808  }
1809 
1810  // Element 2 - "front"
1811  {
1812  Elem* elem2 = mesh.add_elem (new Hex8);
1813  elem2->set_node(0) = nodes[8];
1814  elem2->set_node(1) = nodes[9];
1815  elem2->set_node(2) = nodes[1];
1816  elem2->set_node(3) = nodes[0];
1817  elem2->set_node(4) = nodes[12];
1818  elem2->set_node(5) = nodes[13];
1819  elem2->set_node(6) = nodes[5];
1820  elem2->set_node(7) = nodes[4];
1821  }
1822 
1823  // Element 3 - "right"
1824  {
1825  Elem* elem3 = mesh.add_elem (new Hex8);
1826  elem3->set_node(0) = nodes[1];
1827  elem3->set_node(1) = nodes[9];
1828  elem3->set_node(2) = nodes[10];
1829  elem3->set_node(3) = nodes[2];
1830  elem3->set_node(4) = nodes[5];
1831  elem3->set_node(5) = nodes[13];
1832  elem3->set_node(6) = nodes[14];
1833  elem3->set_node(7) = nodes[6];
1834  }
1835 
1836  // Element 4 - "back"
1837  {
1838  Elem* elem4 = mesh.add_elem (new Hex8);
1839  elem4->set_node(0) = nodes[3];
1840  elem4->set_node(1) = nodes[2];
1841  elem4->set_node(2) = nodes[10];
1842  elem4->set_node(3) = nodes[11];
1843  elem4->set_node(4) = nodes[7];
1844  elem4->set_node(5) = nodes[6];
1845  elem4->set_node(6) = nodes[14];
1846  elem4->set_node(7) = nodes[15];
1847  }
1848 
1849  // Element 5 - "left"
1850  {
1851  Elem* elem5 = mesh.add_elem (new Hex8);
1852  elem5->set_node(0) = nodes[8];
1853  elem5->set_node(1) = nodes[0];
1854  elem5->set_node(2) = nodes[3];
1855  elem5->set_node(3) = nodes[11];
1856  elem5->set_node(4) = nodes[12];
1857  elem5->set_node(5) = nodes[4];
1858  elem5->set_node(6) = nodes[7];
1859  elem5->set_node(7) = nodes[15];
1860  }
1861 
1862  // Element 6 - "top"
1863  {
1864  Elem* elem6 = mesh.add_elem (new Hex8);
1865  elem6->set_node(0) = nodes[4];
1866  elem6->set_node(1) = nodes[5];
1867  elem6->set_node(2) = nodes[6];
1868  elem6->set_node(3) = nodes[7];
1869  elem6->set_node(4) = nodes[12];
1870  elem6->set_node(5) = nodes[13];
1871  elem6->set_node(6) = nodes[14];
1872  elem6->set_node(7) = nodes[15];
1873  }
1874 
1875  break;
1876  } // end case 3
1877 
1878  default:
1879  libmesh_error();
1880 
1881 
1882 
1883  } // end switch (dim)
1884 
1885  // Now we have the beginnings of a sphere.
1886  // Add some more elements by doing uniform refinements and
1887  // popping nodes to the boundary.
1888  MeshRefinement mesh_refinement (mesh);
1889 
1890  // Loop over the elements, refine, pop nodes to boundary.
1891  for (unsigned int r=0; r<nr; r++)
1892  {
1893  mesh_refinement.uniformly_refine(1);
1894 
1895  MeshBase::element_iterator it = mesh.active_elements_begin();
1896  const MeshBase::element_iterator end = mesh.active_elements_end();
1897 
1898  for (; it != end; ++it)
1899  {
1900  Elem* elem = *it;
1901 
1902  for (unsigned int s=0; s<elem->n_sides(); s++)
1903  if (elem->neighbor(s) == NULL || (mesh.mesh_dimension() == 2 && !flat))
1904  {
1905  AutoPtr<Elem> side(elem->build_side(s));
1906 
1907  // Pop each point to the sphere boundary
1908  for (unsigned int n=0; n<side->n_nodes(); n++)
1909  side->point(n) =
1910  sphere.closest_point(side->point(n));
1911  }
1912  }
1913  }
1914 
1915  // The mesh now contains a refinement hierarchy due to the refinements
1916  // used to generate the grid. In order to call other support functions
1917  // like all_tri() and all_second_order, you need a "flat" mesh file (with no
1918  // refinement trees) so
1920 
1921  // In 2D, convert all the quads to triangles if requested
1922  if (mesh.mesh_dimension()==2)
1923  {
1924  if ((type == TRI6) || (type == TRI3))
1925  {
1927  }
1928  }
1929 
1930 
1931  // Convert to second-order elements if the user requested it.
1933  {
1934  // type is already a second-order, determine if it is the
1935  // "full-ordered" second-order element, or the "serendipity"
1936  // second order element. Note also that all_second_order
1937  // can't be called once the mesh has been refined.
1938  bool full_ordered = !((type==QUAD8) || (type==HEX20));
1939  mesh.all_second_order(full_ordered);
1940 
1941  // And pop to the boundary again...
1942  MeshBase::element_iterator it = mesh.active_elements_begin();
1943  const MeshBase::element_iterator end = mesh.active_elements_end();
1944 
1945  for (; it != end; ++it)
1946  {
1947  Elem* elem = *it;
1948 
1949  for (unsigned int s=0; s<elem->n_sides(); s++)
1950  if (elem->neighbor(s) == NULL)
1951  {
1952  AutoPtr<Elem> side(elem->build_side(s));
1953 
1954  // Pop each point to the sphere boundary
1955  for (unsigned int n=0; n<side->n_nodes(); n++)
1956  side->point(n) =
1957  sphere.closest_point(side->point(n));
1958  }
1959  }
1960  }
1961 
1962 
1963  // The meshes could probably use some smoothing.
1964  LaplaceMeshSmoother smoother(mesh);
1965  smoother.smooth(n_smooth);
1966 
1967  // We'll give the whole sphere surface a boundary id of 0
1968  {
1969  MeshBase::element_iterator it = mesh.active_elements_begin();
1970  const MeshBase::element_iterator end = mesh.active_elements_end();
1971 
1972  for (; it != end; ++it)
1973  {
1974  Elem* elem = *it;
1975  for (unsigned int s=0; s != elem->n_sides(); ++s)
1976  if (!elem->neighbor(s))
1977  mesh.boundary_info->add_side(elem, s, 0);
1978  }
1979  }
1980 
1981  STOP_LOG("build_sphere()", "MeshTools::Generation");
1982 
1983 
1984  // Done building the mesh. Now prepare it for use.
1985  mesh.prepare_for_use(/*skip_renumber =*/ false);
1986 }
1987 
1988 #endif // #ifndef LIBMESH_ENABLE_AMR
1989 
1990 
1991 // Meshes the tensor product of a 1D and a 1D-or-2D domain.
1993  const MeshBase& cross_section,
1994  const unsigned int nz,
1995  RealVectorValue extrusion_vector)
1996 {
1997  if (!cross_section.n_elem())
1998  return;
1999 
2000  START_LOG("build_extrusion()", "MeshTools::Generation");
2001 
2002  dof_id_type orig_elem = cross_section.n_elem();
2003  dof_id_type orig_nodes = cross_section.n_nodes();
2004 
2005  unsigned int order = 1;
2006 
2007  // If cross_section is distributed, so is its extrusion
2008  if (!cross_section.is_serial())
2009  mesh.delete_remote_elements();
2010 
2011  // We know a priori how many elements we'll need
2012  mesh.reserve_elem(nz*orig_elem);
2013 
2014  // For straightforward meshes we need one or two additional layers per
2015  // element.
2016  if ((*cross_section.elements_begin())->default_order() == SECOND)
2017  order = 2;
2018 
2019  mesh.reserve_nodes((order*nz+1)*orig_nodes);
2020 
2021  MeshBase::const_node_iterator nd = cross_section.nodes_begin();
2022  const MeshBase::const_node_iterator nend = cross_section.nodes_end();
2023  for (; nd!=nend; ++nd)
2024  {
2025  const Node* node = *nd;
2026 
2027  for (unsigned int k=0; k != order*nz+1; ++k)
2028  {
2029  Node *new_node =
2030  mesh.add_point(*node +
2031  (extrusion_vector * k / nz / order),
2032  node->id() + (k * orig_nodes),
2033  node->processor_id());
2034 
2035  const std::vector<boundary_id_type> ids_to_copy =
2036  cross_section.boundary_info->boundary_ids(node);
2037 
2038  mesh.boundary_info->add_node(new_node, ids_to_copy);
2039  }
2040  }
2041 
2042  const std::set<boundary_id_type> &side_ids =
2043  cross_section.boundary_info->get_side_boundary_ids();
2044  const boundary_id_type next_side_id = side_ids.empty() ?
2045  0 : *side_ids.rbegin() + 1;
2046 
2047  MeshBase::const_element_iterator el = cross_section.elements_begin();
2048  const MeshBase::const_element_iterator end = cross_section.elements_end();
2049  for (; el!=end; ++el)
2050  {
2051  const Elem* elem = *el;
2052  const ElemType etype = elem->type();
2053 
2054  // build_extrusion currently only works on coarse meshes
2055  libmesh_assert (!elem->parent());
2056 
2057  // We need a map from low-D to high-D sides for boundary id
2058  // setting
2059  std::vector<unsigned char> sidemap(4);
2060 
2061  for (unsigned int k=0; k != nz; ++k)
2062  {
2063  Elem *new_elem;
2064  switch (etype)
2065  {
2066  case EDGE2:
2067  {
2068  new_elem = new Quad4;
2069  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
2070  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
2071  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
2072  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
2073  break;
2074  }
2075  case EDGE3:
2076  {
2077  new_elem = new Quad9;
2078  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
2079  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
2080  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
2081  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
2082  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
2083  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
2084  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
2085  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
2086  new_elem->set_node(8) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
2087  break;
2088  }
2089  case TRI3:
2090  {
2091  new_elem = new Prism6;
2092  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
2093  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
2094  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
2095  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
2096  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
2097  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
2098  break;
2099  }
2100  case TRI6:
2101  {
2102  new_elem = new Prism18;
2103  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
2104  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
2105  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
2106  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
2107  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
2108  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
2109  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
2110  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
2111  new_elem->set_node(8) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
2112  new_elem->set_node(9) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
2113  new_elem->set_node(10) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
2114  new_elem->set_node(11) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
2115  new_elem->set_node(12) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
2116  new_elem->set_node(13) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
2117  new_elem->set_node(14) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
2118  new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
2119  new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
2120  new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
2121  break;
2122  }
2123  case QUAD4:
2124  {
2125  new_elem = new Hex8;
2126  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
2127  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
2128  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
2129  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (k * orig_nodes));
2130  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
2131  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
2132  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
2133  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((k+1) * orig_nodes));
2134  break;
2135  }
2136  case QUAD9:
2137  {
2138  new_elem = new Hex27;
2139  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
2140  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
2141  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
2142  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
2143  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
2144  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
2145  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
2146  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
2147  new_elem->set_node(8) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
2148  new_elem->set_node(9) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
2149  new_elem->set_node(10) = mesh.node_ptr(elem->get_node(6)->id() + (2*k * orig_nodes));
2150  new_elem->set_node(11) = mesh.node_ptr(elem->get_node(7)->id() + (2*k * orig_nodes));
2151  new_elem->set_node(12) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
2152  new_elem->set_node(13) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
2153  new_elem->set_node(14) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
2154  new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
2155  new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
2156  new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
2157  new_elem->set_node(18) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+2) * orig_nodes));
2158  new_elem->set_node(19) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+2) * orig_nodes));
2159  new_elem->set_node(20) = mesh.node_ptr(elem->get_node(8)->id() + (2*k * orig_nodes));
2160  new_elem->set_node(21) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
2161  new_elem->set_node(22) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
2162  new_elem->set_node(23) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+1) * orig_nodes));
2163  new_elem->set_node(24) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+1) * orig_nodes));
2164  new_elem->set_node(25) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+2) * orig_nodes));
2165  new_elem->set_node(26) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+1) * orig_nodes));
2166  break;
2167  }
2168  default:
2169  {
2170  libmesh_not_implemented();
2171  break;
2172  }
2173  }
2174 
2175  new_elem->set_id(elem->id() + (k * orig_elem));
2176  new_elem->processor_id() = elem->processor_id();
2177 
2178  // maintain the subdomain_id
2179  new_elem->subdomain_id() = elem->subdomain_id();
2180 
2181  new_elem = mesh.add_elem(new_elem);
2182 
2183  // Copy any old boundary ids on all sides
2184  for (unsigned int s = 0; s != elem->n_sides(); ++s)
2185  {
2186  const std::vector<boundary_id_type> ids_to_copy =
2187  cross_section.boundary_info->boundary_ids(elem, s);
2188 
2189  mesh.boundary_info->add_side(new_elem, s+1, ids_to_copy);
2190  }
2191 
2192  // Give new boundary ids to bottom and top
2193  if (k == 0)
2194  mesh.boundary_info->add_side(new_elem, 0, next_side_id);
2195  if (k == nz-1)
2196  mesh.boundary_info->add_side(new_elem, elem->n_sides()+1, next_side_id+1);
2197  }
2198  }
2199 
2200  STOP_LOG("build_extrusion()", "MeshTools::Generation");
2201 
2202  // Done building the mesh. Now prepare it for use.
2203  mesh.prepare_for_use(/*skip_renumber =*/ false);
2204 }
2205 
2206 
2207 
2208 
2209 #ifdef LIBMESH_HAVE_TRIANGLE
2210 
2211 // Triangulates a 2D rectangular region with or without holes
2213  const unsigned int nx, // num. of elements in x-dir
2214  const unsigned int ny, // num. of elements in y-dir
2215  const Real xmin, const Real xmax,
2216  const Real ymin, const Real ymax,
2217  const ElemType type,
2218  const std::vector<TriangleInterface::Hole*>* holes)
2219 {
2220  // Check for reasonable size
2221  libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
2222  libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
2223  libmesh_assert_less (xmin, xmax);
2224  libmesh_assert_less (ymin, ymax);
2225 
2226  // Clear out any data which may have been in the Mesh
2227  mesh.clear();
2228 
2229  // Make sure the new Mesh will be 2D
2230  mesh.set_mesh_dimension(2);
2231 
2232  // The x and y spacing between boundary points
2233  const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
2234  const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
2235 
2236  // Bottom
2237  for (unsigned int p=0; p<=nx; ++p)
2238  mesh.add_point(Point(xmin + p*delta_x, ymin));
2239 
2240  // Right side
2241  for (unsigned int p=1; p<ny; ++p)
2242  mesh.add_point(Point(xmax, ymin + p*delta_y));
2243 
2244  // Top
2245  for (unsigned int p=0; p<=nx; ++p)
2246  mesh.add_point(Point(xmax - p*delta_x, ymax));
2247 
2248  // Left side
2249  for (unsigned int p=1; p<ny; ++p)
2250  mesh.add_point(Point(xmin, ymax - p*delta_y));
2251 
2252  // Be sure we added as many points as we thought we did
2253  libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
2254 
2255  // Construct the Triangle Interface object
2256  TriangleInterface t(mesh);
2257 
2258  // Set custom variables for the triangulation
2259  t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
2261  t.elem_type() = type;
2262 
2263  if (holes != NULL)
2264  t.attach_hole_list(holes);
2265 
2266  // Triangulate!
2267  t.triangulate();
2268 
2269  // The mesh is now generated, but we still need to mark the boundaries
2270  // to be consistent with the other build_square routines. Note that all
2271  // hole boundary elements get the same ID, 4.
2273  const MeshBase::element_iterator end_el = mesh.elements_end();
2274  for ( ; el != end_el; ++el)
2275  {
2276  const Elem* elem = *el;
2277 
2278  for (unsigned int s=0; s<elem->n_sides(); s++)
2279  if (elem->neighbor(s) == NULL)
2280  {
2281  AutoPtr<Elem> side (elem->build_side(s));
2282 
2283  // Check the location of the side's midpoint. Since
2284  // the square has straight sides, the midpoint is not
2285  // on the corner and thus it is uniquely on one of the
2286  // sides.
2287  Point side_midpoint= 0.5f*( (*side->get_node(0)) + (*side->get_node(1)) );
2288 
2289  // The boundary ids are set following the same convention as Quad4 sides
2290  // bottom = 0
2291  // right = 1
2292  // top = 2
2293  // left = 3
2294  // hole = 4
2296 
2297  // bottom
2298  if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
2299  bc_id=0;
2300 
2301  // right
2302  else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
2303  bc_id=1;
2304 
2305  // top
2306  else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
2307  bc_id=2;
2308 
2309  // left
2310  else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
2311  bc_id=3;
2312 
2313  // If the point is not on any of the external boundaries, it
2314  // is on one of the holes....
2315 
2316  // Finally, add this element's information to the boundary info object.
2317  mesh.boundary_info->add_side(elem->id(), s, bc_id);
2318  }
2319  }
2320 
2321 } // end build_delaunay_square
2322 
2323 #endif // LIBMESH_HAVE_TRIANGLE
2324 
2325 
2326 
2327 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo