mesh_generation.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 00020 // C++ includes 00021 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC 00022 #include <cmath> // for std::sqrt 00023 00024 00025 // Local includes 00026 #include "libmesh/mesh_generation.h" 00027 #include "libmesh/unstructured_mesh.h" 00028 // #include "libmesh/elem.h" 00029 #include "libmesh/mesh_refinement.h" 00030 #include "libmesh/edge_edge2.h" 00031 #include "libmesh/edge_edge3.h" 00032 #include "libmesh/edge_edge4.h" 00033 #include "libmesh/face_tri3.h" 00034 #include "libmesh/face_tri6.h" 00035 #include "libmesh/face_quad4.h" 00036 #include "libmesh/face_quad8.h" 00037 #include "libmesh/face_quad9.h" 00038 #include "libmesh/cell_hex8.h" 00039 #include "libmesh/cell_hex20.h" 00040 #include "libmesh/cell_hex27.h" 00041 #include "libmesh/cell_prism6.h" 00042 #include "libmesh/cell_prism15.h" 00043 #include "libmesh/cell_prism18.h" 00044 #include "libmesh/cell_tet4.h" 00045 #include "libmesh/cell_pyramid5.h" 00046 #include "libmesh/libmesh_logging.h" 00047 #include "libmesh/boundary_info.h" 00048 #include "libmesh/sphere.h" 00049 #include "libmesh/mesh_modification.h" 00050 #include "libmesh/mesh_smoother_laplace.h" 00051 #include "libmesh/node_elem.h" 00052 #include "libmesh/vector_value.h" 00053 00054 namespace libMesh 00055 { 00056 00057 namespace MeshTools { 00058 namespace Generation { 00059 namespace Private { 00067 inline 00068 unsigned int idx(const ElemType type, 00069 const unsigned int nx, 00070 const unsigned int i, 00071 const unsigned int j) 00072 { 00073 switch(type) 00074 { 00075 case INVALID_ELEM: 00076 case QUAD4: 00077 case TRI3: 00078 { 00079 return i + j*(nx+1); 00080 } 00081 00082 case QUAD8: 00083 case QUAD9: 00084 case TRI6: 00085 { 00086 return i + j*(2*nx+1); 00087 } 00088 00089 default: 00090 libmesh_error_msg("ERROR: Unrecognized 2D element type."); 00091 } 00092 00093 return libMesh::invalid_uint; 00094 } 00095 00096 00097 00098 // Same as the function above, but for 3D elements 00099 inline 00100 unsigned int idx(const ElemType type, 00101 const unsigned int nx, 00102 const unsigned int ny, 00103 const unsigned int i, 00104 const unsigned int j, 00105 const unsigned int k) 00106 { 00107 switch(type) 00108 { 00109 case INVALID_ELEM: 00110 case HEX8: 00111 case PRISM6: 00112 { 00113 return i + (nx+1)*(j + k*(ny+1)); 00114 } 00115 00116 case HEX20: 00117 case HEX27: 00118 case TET4: // TET4's are created from an initial HEX27 discretization 00119 case TET10: // TET10's are created from an initial HEX27 discretization 00120 case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization 00121 case PYRAMID13: 00122 case PYRAMID14: 00123 case PRISM15: 00124 case PRISM18: 00125 { 00126 return i + (2*nx+1)*(j + k*(2*ny+1)); 00127 } 00128 00129 default: 00130 libmesh_error_msg("ERROR: Unrecognized element type."); 00131 } 00132 00133 return libMesh::invalid_uint; 00134 } 00135 } // namespace Private 00136 } // namespace Generation 00137 } // namespace MeshTools 00138 00139 // ------------------------------------------------------------ 00140 // MeshTools::Generation function for mesh generation 00141 void MeshTools::Generation::build_cube(UnstructuredMesh& mesh, 00142 const unsigned int nx, 00143 const unsigned int ny, 00144 const unsigned int nz, 00145 const Real xmin, const Real xmax, 00146 const Real ymin, const Real ymax, 00147 const Real zmin, const Real zmax, 00148 const ElemType type, 00149 const bool gauss_lobatto_grid) 00150 { 00151 START_LOG("build_cube()", "MeshTools::Generation"); 00152 00153 // Declare that we are using the indexing utility routine 00154 // in the "Private" part of our current namespace. If this doesn't 00155 // work in GCC 2.95.3 we can either remove it or stop supporting 00156 // 2.95.3 altogether. 00157 // Changing this to import the whole namespace... just importing idx 00158 // causes an internal compiler error for Intel Compiler 11.0 on Linux 00159 // in debug mode. 00160 using namespace MeshTools::Generation::Private; 00161 00162 // Clear the mesh and start from scratch 00163 mesh.clear(); 00164 00165 BoundaryInfo& boundary_info = mesh.get_boundary_info(); 00166 00167 if (nz != 0) 00168 mesh.set_mesh_dimension(3); 00169 else if (ny != 0) 00170 mesh.set_mesh_dimension(2); 00171 else if (nx != 0) 00172 mesh.set_mesh_dimension(1); 00173 else 00174 mesh.set_mesh_dimension(0); 00175 00176 switch (mesh.mesh_dimension()) 00177 { 00178 //--------------------------------------------------------------------- 00179 // Build a 0D point 00180 case 0: 00181 { 00182 libmesh_assert_equal_to (nx, 0); 00183 libmesh_assert_equal_to (ny, 0); 00184 libmesh_assert_equal_to (nz, 0); 00185 00186 libmesh_assert (type == INVALID_ELEM || type == NODEELEM); 00187 00188 // Build one nodal element for the mesh 00189 mesh.add_point (Point(0, 0, 0), 0); 00190 Elem* elem = mesh.add_elem (new NodeElem); 00191 elem->set_node(0) = mesh.node_ptr(0); 00192 00193 break; 00194 } 00195 00196 00197 00198 //--------------------------------------------------------------------- 00199 // Build a 1D line 00200 case 1: 00201 { 00202 libmesh_assert_not_equal_to (nx, 0); 00203 libmesh_assert_equal_to (ny, 0); 00204 libmesh_assert_equal_to (nz, 0); 00205 libmesh_assert_less (xmin, xmax); 00206 00207 // Reserve elements 00208 switch (type) 00209 { 00210 case INVALID_ELEM: 00211 case EDGE2: 00212 case EDGE3: 00213 case EDGE4: 00214 { 00215 mesh.reserve_elem (nx); 00216 break; 00217 } 00218 00219 default: 00220 libmesh_error_msg("ERROR: Unrecognized 1D element type."); 00221 } 00222 00223 // Reserve nodes 00224 switch (type) 00225 { 00226 case INVALID_ELEM: 00227 case EDGE2: 00228 { 00229 mesh.reserve_nodes(nx+1); 00230 break; 00231 } 00232 00233 case EDGE3: 00234 { 00235 mesh.reserve_nodes(2*nx+1); 00236 break; 00237 } 00238 00239 case EDGE4: 00240 { 00241 mesh.reserve_nodes(3*nx+1); 00242 break; 00243 } 00244 00245 default: 00246 libmesh_error_msg("ERROR: Unrecognized 1D element type."); 00247 } 00248 00249 00250 // Build the nodes, depends on whether we're using linears, 00251 // quadratics or cubics and whether using uniform grid or Gauss-Lobatto 00252 unsigned int node_id = 0; 00253 switch(type) 00254 { 00255 case INVALID_ELEM: 00256 case EDGE2: 00257 { 00258 for (unsigned int i=0; i<=nx; i++) 00259 { 00260 if (gauss_lobatto_grid) 00261 mesh.add_point (Point(0.5*(std::cos(libMesh::pi*static_cast<Real>(nx-i)/static_cast<Real>(nx))+1.0), 00262 0, 00263 0), node_id++); 00264 else 00265 mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx), 00266 0, 00267 0), node_id++); 00268 } 00269 break; 00270 } 00271 00272 case EDGE3: 00273 { 00274 for (unsigned int i=0; i<=2*nx; i++) 00275 { 00276 if (gauss_lobatto_grid) 00277 { 00278 // The x location of the point. 00279 Real x=0.; 00280 00281 // Shortcut quantities (do not depend on i) 00282 const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) ); 00283 00284 // If i is even, compute a normal Gauss-Lobatto point 00285 if (i%2 == 0) 00286 x = 0.5*(1.0 - c); 00287 00288 // Otherwise, it is the average of the previous and next points 00289 else 00290 { 00291 Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(2*nx) ); 00292 Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(2*nx) ); 00293 00294 Real gl_xmin = 0.5*(1.0 - cmin); 00295 Real gl_xmax = 0.5*(1.0 - cmax); 00296 x = 0.5*(gl_xmin + gl_xmax); 00297 } 00298 00299 mesh.add_point (Point(x,0.,0.), node_id++); 00300 } 00301 else 00302 mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx), 00303 0, 00304 0), node_id++); 00305 } 00306 break; 00307 } 00308 00309 case EDGE4: 00310 { 00311 for (unsigned int i=0; i<=3*nx; i++) 00312 { 00313 if (gauss_lobatto_grid) 00314 { 00315 // The x location of the point 00316 Real x=0.; 00317 00318 // Shortcut quantities 00319 const Real c = std::cos( libMesh::pi*i / static_cast<Real>(3*nx) ); 00320 00321 // If i is multiple of 3, compute a normal Gauss-Lobatto point 00322 if (i%3 == 0) 00323 x = 0.5*(1.0 - c); 00324 00325 // Otherwise, distribute points evenly within the element 00326 else 00327 { 00328 if(i%3 == 1) 00329 { 00330 Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(3*nx) ); 00331 Real cmax = std::cos( libMesh::pi*(i+2) / static_cast<Real>(3*nx) ); 00332 00333 Real gl_xmin = 0.5*(1.0 - cmin); 00334 Real gl_xmax = 0.5*(1.0 - cmax); 00335 00336 x = (2.*gl_xmin + gl_xmax)/3.; 00337 } 00338 else 00339 if(i%3 == 2) 00340 { 00341 Real cmin = std::cos( libMesh::pi*(i-2) / static_cast<Real>(3*nx) ); 00342 Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(3*nx) ); 00343 00344 Real gl_xmin = 0.5*(1.0 - cmin); 00345 Real gl_xmax = 0.5*(1.0 - cmax); 00346 00347 x = (gl_xmin + 2.*gl_xmax)/3.; 00348 } 00349 00350 } 00351 00352 mesh.add_point (Point(x,0.,0.), node_id++); 00353 } 00354 else 00355 mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(3*nx), 00356 0, 00357 0), node_id++); 00358 } 00359 00360 00361 00362 break; 00363 } 00364 00365 default: 00366 libmesh_error_msg("ERROR: Unrecognized 1D element type."); 00367 00368 } 00369 00370 // Build the elements of the mesh 00371 switch(type) 00372 { 00373 case INVALID_ELEM: 00374 case EDGE2: 00375 { 00376 for (unsigned int i=0; i<nx; i++) 00377 { 00378 Elem* elem = mesh.add_elem (new Edge2); 00379 elem->set_node(0) = mesh.node_ptr(i); 00380 elem->set_node(1) = mesh.node_ptr(i+1); 00381 00382 if (i == 0) 00383 boundary_info.add_side(elem, 0, 0); 00384 00385 if (i == (nx-1)) 00386 boundary_info.add_side(elem, 1, 1); 00387 } 00388 break; 00389 } 00390 00391 case EDGE3: 00392 { 00393 for (unsigned int i=0; i<nx; i++) 00394 { 00395 Elem* elem = mesh.add_elem (new Edge3); 00396 elem->set_node(0) = mesh.node_ptr(2*i); 00397 elem->set_node(2) = mesh.node_ptr(2*i+1); 00398 elem->set_node(1) = mesh.node_ptr(2*i+2); 00399 00400 if (i == 0) 00401 boundary_info.add_side(elem, 0, 0); 00402 00403 if (i == (nx-1)) 00404 boundary_info.add_side(elem, 1, 1); 00405 } 00406 break; 00407 } 00408 00409 case EDGE4: 00410 { 00411 for (unsigned int i=0; i<nx; i++) 00412 { 00413 Elem* elem = mesh.add_elem (new Edge4); 00414 elem->set_node(0) = mesh.node_ptr(3*i); 00415 elem->set_node(2) = mesh.node_ptr(3*i+1); 00416 elem->set_node(3) = mesh.node_ptr(3*i+2); 00417 elem->set_node(1) = mesh.node_ptr(3*i+3); 00418 00419 if (i == 0) 00420 boundary_info.add_side(elem, 0, 0); 00421 00422 if (i == (nx-1)) 00423 boundary_info.add_side(elem, 1, 1); 00424 } 00425 break; 00426 } 00427 00428 default: 00429 libmesh_error_msg("ERROR: Unrecognized 1D element type."); 00430 } 00431 00432 // Scale the nodal positions 00433 for (unsigned int p=0; p<mesh.n_nodes(); p++) 00434 mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin; 00435 00436 // Add sideset names to boundary info 00437 boundary_info.sideset_name(0) = "left"; 00438 boundary_info.sideset_name(1) = "right"; 00439 00440 // Add nodeset names to boundary info 00441 boundary_info.nodeset_name(0) = "left"; 00442 boundary_info.nodeset_name(1) = "right"; 00443 00444 break; 00445 } 00446 00447 00448 00449 00450 00451 00452 00453 00454 00455 00456 //--------------------------------------------------------------------- 00457 // Build a 2D quadrilateral 00458 case 2: 00459 { 00460 libmesh_assert_not_equal_to (nx, 0); 00461 libmesh_assert_not_equal_to (ny, 0); 00462 libmesh_assert_equal_to (nz, 0); 00463 libmesh_assert_less (xmin, xmax); 00464 libmesh_assert_less (ymin, ymax); 00465 00466 // Reserve elements. The TRI3 and TRI6 meshes 00467 // have twice as many elements... 00468 switch (type) 00469 { 00470 case INVALID_ELEM: 00471 case QUAD4: 00472 case QUAD8: 00473 case QUAD9: 00474 { 00475 mesh.reserve_elem (nx*ny); 00476 break; 00477 } 00478 00479 case TRI3: 00480 case TRI6: 00481 { 00482 mesh.reserve_elem (2*nx*ny); 00483 break; 00484 } 00485 00486 default: 00487 libmesh_error_msg("ERROR: Unrecognized 2D element type."); 00488 } 00489 00490 00491 00492 // Reserve nodes. The quadratic element types 00493 // need to reserve more nodes than the linear types. 00494 switch (type) 00495 { 00496 case INVALID_ELEM: 00497 case QUAD4: 00498 case TRI3: 00499 { 00500 mesh.reserve_nodes( (nx+1)*(ny+1) ); 00501 break; 00502 } 00503 00504 case QUAD8: 00505 case QUAD9: 00506 case TRI6: 00507 { 00508 mesh.reserve_nodes( (2*nx+1)*(2*ny+1) ); 00509 break; 00510 } 00511 00512 00513 default: 00514 libmesh_error_msg("ERROR: Unrecognized 2D element type."); 00515 } 00516 00517 00518 00519 // Build the nodes. Depends on whether you are using a linear 00520 // or quadratic element, and whether you are using a uniform 00521 // grid or the Gauss-Lobatto grid points. 00522 unsigned int node_id = 0; 00523 switch (type) 00524 { 00525 case INVALID_ELEM: 00526 case QUAD4: 00527 case TRI3: 00528 { 00529 for (unsigned int j=0; j<=ny; j++) 00530 for (unsigned int i=0; i<=nx; i++) 00531 { 00532 if (gauss_lobatto_grid) 00533 { 00534 mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))), 00535 0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))), 00536 0.), node_id++); 00537 } 00538 00539 else 00540 mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx), 00541 static_cast<Real>(j)/static_cast<Real>(ny), 00542 0.), node_id++); 00543 } 00544 00545 break; 00546 } 00547 00548 case QUAD8: 00549 case QUAD9: 00550 case TRI6: 00551 { 00552 for (unsigned int j=0; j<=(2*ny); j++) 00553 for (unsigned int i=0; i<=(2*nx); i++) 00554 { 00555 if (gauss_lobatto_grid) 00556 { 00557 // The x,y locations of the point. 00558 Real x=0., y=0.; 00559 00560 // Shortcut quantities (do not depend on i,j) 00561 const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) ); 00562 const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) ); 00563 00564 // Shortcut quantities (depend on i,j) 00565 const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) ); 00566 const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) ); 00567 00568 // If i is even, compute a normal Gauss-Lobatto point 00569 if (i%2 == 0) 00570 x = 0.5*(1.0 - c); 00571 00572 // Otherwise, it is the average of the previous and next points 00573 else 00574 x = 0.5*(1.0 - a*c); 00575 00576 // If j is even, compute a normal Gauss-Lobatto point 00577 if (j%2 == 0) 00578 y = 0.5*(1.0 - d); 00579 00580 // Otherwise, it is the average of the previous and next points 00581 else 00582 y = 0.5*(1.0 - b*d); 00583 00584 00585 mesh.add_point (Point(x,y,0.), node_id++); 00586 } 00587 00588 00589 else 00590 mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx), 00591 static_cast<Real>(j)/static_cast<Real>(2*ny), 00592 0), node_id++); 00593 } 00594 00595 break; 00596 } 00597 00598 00599 default: 00600 libmesh_error_msg("ERROR: Unrecognized 2D element type."); 00601 } 00602 00603 00604 00605 00606 00607 00608 // Build the elements. Each one is a bit different. 00609 switch (type) 00610 { 00611 00612 case INVALID_ELEM: 00613 case QUAD4: 00614 { 00615 for (unsigned int j=0; j<ny; j++) 00616 for (unsigned int i=0; i<nx; i++) 00617 { 00618 Elem* elem = mesh.add_elem(new Quad4); 00619 00620 elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) ); 00621 elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) ); 00622 elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1)); 00623 elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+1) ); 00624 00625 if (j == 0) 00626 boundary_info.add_side(elem, 0, 0); 00627 00628 if (j == (ny-1)) 00629 boundary_info.add_side(elem, 2, 2); 00630 00631 if (i == 0) 00632 boundary_info.add_side(elem, 3, 3); 00633 00634 if (i == (nx-1)) 00635 boundary_info.add_side(elem, 1, 1); 00636 } 00637 break; 00638 } 00639 00640 00641 case TRI3: 00642 { 00643 for (unsigned int j=0; j<ny; j++) 00644 for (unsigned int i=0; i<nx; i++) 00645 { 00646 Elem* elem = NULL; 00647 00648 // Add first Tri3 00649 elem = mesh.add_elem(new Tri3); 00650 00651 elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) ); 00652 elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) ); 00653 elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1)); 00654 00655 if (j == 0) 00656 boundary_info.add_side(elem, 0, 0); 00657 00658 if (i == (nx-1)) 00659 boundary_info.add_side(elem, 1, 1); 00660 00661 // Add second Tri3 00662 elem = mesh.add_elem(new Tri3); 00663 00664 elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) ); 00665 elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j+1)); 00666 elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+1) ); 00667 00668 if (j == (ny-1)) 00669 boundary_info.add_side(elem, 1, 2); 00670 00671 if (i == 0) 00672 boundary_info.add_side(elem, 2, 3); 00673 } 00674 break; 00675 } 00676 00677 00678 00679 case QUAD8: 00680 case QUAD9: 00681 { 00682 for (unsigned int j=0; j<(2*ny); j += 2) 00683 for (unsigned int i=0; i<(2*nx); i += 2) 00684 { 00685 Elem* elem = (type == QUAD8) ? 00686 mesh.add_elem(new Quad8) : 00687 mesh.add_elem(new Quad9); 00688 00689 00690 elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) ); 00691 elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j) ); 00692 elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2)); 00693 elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+2) ); 00694 elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j) ); 00695 elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+2,j+1)); 00696 elem->set_node(6) = mesh.node_ptr(idx(type,nx,i+1,j+2)); 00697 elem->set_node(7) = mesh.node_ptr(idx(type,nx,i,j+1) ); 00698 if (type == QUAD9) 00699 elem->set_node(8) = mesh.node_ptr(idx(type,nx,i+1,j+1)); 00700 00701 00702 if (j == 0) 00703 boundary_info.add_side(elem, 0, 0); 00704 00705 if (j == 2*(ny-1)) 00706 boundary_info.add_side(elem, 2, 2); 00707 00708 if (i == 0) 00709 boundary_info.add_side(elem, 3, 3); 00710 00711 if (i == 2*(nx-1)) 00712 boundary_info.add_side(elem, 1, 1); 00713 } 00714 break; 00715 } 00716 00717 00718 case TRI6: 00719 { 00720 for (unsigned int j=0; j<(2*ny); j += 2) 00721 for (unsigned int i=0; i<(2*nx); i += 2) 00722 { 00723 Elem* elem = NULL; 00724 00725 // Add first Tri6 00726 elem = mesh.add_elem(new Tri6); 00727 00728 elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) ); 00729 elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j) ); 00730 elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2)); 00731 elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j) ); 00732 elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+2,j+1)); 00733 elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+1,j+1)); 00734 00735 if (j == 0) 00736 boundary_info.add_side(elem, 0, 0); 00737 00738 if (i == 2*(nx-1)) 00739 boundary_info.add_side(elem, 1, 1); 00740 00741 // Add second Tri6 00742 elem = mesh.add_elem(new Tri6); 00743 00744 elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) ); 00745 elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j+2)); 00746 elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+2) ); 00747 elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j+1)); 00748 elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j+2)); 00749 elem->set_node(5) = mesh.node_ptr(idx(type,nx,i,j+1) ); 00750 00751 if (j == 2*(ny-1)) 00752 boundary_info.add_side(elem, 1, 2); 00753 00754 if (i == 0) 00755 boundary_info.add_side(elem, 2, 3); 00756 00757 } 00758 break; 00759 }; 00760 00761 00762 default: 00763 libmesh_error_msg("ERROR: Unrecognized 2D element type."); 00764 } 00765 00766 00767 00768 00769 // Scale the nodal positions 00770 for (unsigned int p=0; p<mesh.n_nodes(); p++) 00771 { 00772 mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin; 00773 mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin; 00774 } 00775 00776 // Add sideset names to boundary info 00777 boundary_info.sideset_name(0) = "bottom"; 00778 boundary_info.sideset_name(1) = "right"; 00779 boundary_info.sideset_name(2) = "top"; 00780 boundary_info.sideset_name(3) = "left"; 00781 00782 // Add nodeset names to boundary info 00783 boundary_info.nodeset_name(0) = "bottom"; 00784 boundary_info.nodeset_name(1) = "right"; 00785 boundary_info.nodeset_name(2) = "top"; 00786 boundary_info.nodeset_name(3) = "left"; 00787 00788 break; 00789 } 00790 00791 00792 00793 00794 00795 00796 00797 00798 00799 00800 00801 //--------------------------------------------------------------------- 00802 // Build a 3D mesh using hexes, tets, prisms, or pyramids. 00803 case 3: 00804 { 00805 libmesh_assert_not_equal_to (nx, 0); 00806 libmesh_assert_not_equal_to (ny, 0); 00807 libmesh_assert_not_equal_to (nz, 0); 00808 libmesh_assert_less (xmin, xmax); 00809 libmesh_assert_less (ymin, ymax); 00810 libmesh_assert_less (zmin, zmax); 00811 00812 00813 // Reserve elements. Meshes with prismatic elements require 00814 // twice as many elements. 00815 switch (type) 00816 { 00817 case INVALID_ELEM: 00818 case HEX8: 00819 case HEX20: 00820 case HEX27: 00821 case TET4: // TET4's are created from an initial HEX27 discretization 00822 case TET10: // TET10's are created from an initial HEX27 discretization 00823 case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization 00824 case PYRAMID13: 00825 case PYRAMID14: 00826 { 00827 mesh.reserve_elem(nx*ny*nz); 00828 break; 00829 } 00830 00831 case PRISM6: 00832 case PRISM15: 00833 case PRISM18: 00834 { 00835 mesh.reserve_elem(2*nx*ny*nz); 00836 break; 00837 } 00838 00839 default: 00840 libmesh_error_msg("ERROR: Unrecognized 3D element type."); 00841 } 00842 00843 00844 00845 00846 00847 // Reserve nodes. Quadratic elements need twice as many nodes as linear elements. 00848 switch (type) 00849 { 00850 case INVALID_ELEM: 00851 case HEX8: 00852 case PRISM6: 00853 { 00854 mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) ); 00855 break; 00856 } 00857 00858 case HEX20: 00859 case HEX27: 00860 case TET4: // TET4's are created from an initial HEX27 discretization 00861 case TET10: // TET10's are created from an initial HEX27 discretization 00862 case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization 00863 case PYRAMID13: 00864 case PYRAMID14: 00865 case PRISM15: 00866 case PRISM18: 00867 { 00868 // FYI: The resulting TET4 mesh will have exactly 00869 // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1 00870 // nodes once the additional mid-edge nodes for the HEX27 discretization 00871 // have been deleted. 00872 mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) ); 00873 break; 00874 } 00875 00876 default: 00877 libmesh_error_msg("ERROR: Unrecognized 3D element type."); 00878 } 00879 00880 00881 00882 00883 // Build the nodes. 00884 unsigned int node_id = 0; 00885 switch (type) 00886 { 00887 case INVALID_ELEM: 00888 case HEX8: 00889 case PRISM6: 00890 { 00891 for (unsigned int k=0; k<=nz; k++) 00892 for (unsigned int j=0; j<=ny; j++) 00893 for (unsigned int i=0; i<=nx; i++) 00894 { 00895 if (gauss_lobatto_grid) 00896 { 00897 mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))), 00898 0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))), 00899 0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(k)/static_cast<Real>(nz)))), node_id++); 00900 } 00901 00902 else 00903 mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(nx), 00904 static_cast<Real>(j)/static_cast<Real>(ny), 00905 static_cast<Real>(k)/static_cast<Real>(nz)), node_id++); 00906 } 00907 break; 00908 } 00909 00910 case HEX20: 00911 case HEX27: 00912 case TET4: // TET4's are created from an initial HEX27 discretization 00913 case TET10: // TET10's are created from an initial HEX27 discretization 00914 case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization 00915 case PYRAMID13: 00916 case PYRAMID14: 00917 case PRISM15: 00918 case PRISM18: 00919 { 00920 for (unsigned int k=0; k<=(2*nz); k++) 00921 for (unsigned int j=0; j<=(2*ny); j++) 00922 for (unsigned int i=0; i<=(2*nx); i++) 00923 { 00924 if (gauss_lobatto_grid) 00925 { 00926 // The x,y locations of the point. 00927 Real x=0., y=0., z=0.; 00928 00929 // Shortcut quantities (do not depend on i,j) 00930 const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) ); 00931 const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) ); 00932 00933 // Shortcut quantities (depend on i,j) 00934 const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) ); 00935 const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) ); 00936 00937 // Additional shortcut quantities (for 3D) 00938 const Real e = std::cos( libMesh::pi / static_cast<Real>(2*nz) ); 00939 const Real f = std::cos( libMesh::pi*k / static_cast<Real>(2*nz) ); 00940 00941 // If i is even, compute a normal Gauss-Lobatto point 00942 if (i%2 == 0) 00943 x = 0.5*(1.0 - c); 00944 00945 // Otherwise, it is the average of the previous and next points 00946 else 00947 x = 0.5*(1.0 - a*c); 00948 00949 // If j is even, compute a normal Gauss-Lobatto point 00950 if (j%2 == 0) 00951 y = 0.5*(1.0 - d); 00952 00953 // Otherwise, it is the average of the previous and next points 00954 else 00955 y = 0.5*(1.0 - b*d); 00956 00957 // If k is even, compute a normal Gauss-Lobatto point 00958 if (k%2 == 0) 00959 z = 0.5*(1.0 - f); 00960 00961 // Otherwise, it is the average of the previous and next points 00962 else 00963 z = 0.5*(1.0 - e*f); 00964 00965 00966 mesh.add_point (Point(x,y,z), node_id++); 00967 } 00968 00969 else 00970 mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(2*nx), 00971 static_cast<Real>(j)/static_cast<Real>(2*ny), 00972 static_cast<Real>(k)/static_cast<Real>(2*nz)), node_id++); 00973 } 00974 break; 00975 } 00976 00977 00978 default: 00979 libmesh_error_msg("ERROR: Unrecognized 3D element type."); 00980 } 00981 00982 00983 00984 00985 // Build the elements. 00986 switch (type) 00987 { 00988 case INVALID_ELEM: 00989 case HEX8: 00990 { 00991 for (unsigned int k=0; k<nz; k++) 00992 for (unsigned int j=0; j<ny; j++) 00993 for (unsigned int i=0; i<nx; i++) 00994 { 00995 Elem* elem = mesh.add_elem(new Hex8); 00996 00997 elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) ); 00998 elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ); 00999 elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ); 01000 elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ); 01001 elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ); 01002 elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ); 01003 elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)); 01004 elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ); 01005 01006 if (k == 0) 01007 boundary_info.add_side(elem, 0, 0); 01008 01009 if (k == (nz-1)) 01010 boundary_info.add_side(elem, 5, 5); 01011 01012 if (j == 0) 01013 boundary_info.add_side(elem, 1, 1); 01014 01015 if (j == (ny-1)) 01016 boundary_info.add_side(elem, 3, 3); 01017 01018 if (i == 0) 01019 boundary_info.add_side(elem, 4, 4); 01020 01021 if (i == (nx-1)) 01022 boundary_info.add_side(elem, 2, 2); 01023 } 01024 break; 01025 } 01026 01027 01028 01029 01030 case PRISM6: 01031 { 01032 for (unsigned int k=0; k<nz; k++) 01033 for (unsigned int j=0; j<ny; j++) 01034 for (unsigned int i=0; i<nx; i++) 01035 { 01036 // First Prism 01037 Elem* elem = NULL; 01038 elem = mesh.add_elem(new Prism6); 01039 01040 elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) ); 01041 elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ); 01042 elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ); 01043 elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ); 01044 elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ); 01045 elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ); 01046 01047 // Add sides for first prism to boundary info object 01048 if (i==0) 01049 boundary_info.add_side(elem, 3, 4); 01050 01051 if (j==0) 01052 boundary_info.add_side(elem, 1, 1); 01053 01054 if (k==0) 01055 boundary_info.add_side(elem, 0, 0); 01056 01057 if (k == (nz-1)) 01058 boundary_info.add_side(elem, 4, 5); 01059 01060 // Second Prism 01061 elem = mesh.add_elem(new Prism6); 01062 01063 elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ); 01064 elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ); 01065 elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ); 01066 elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ); 01067 elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)); 01068 elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ); 01069 01070 // Add sides for second prism to boundary info object 01071 if (i == (nx-1)) 01072 boundary_info.add_side(elem, 1, 2); 01073 01074 if (j == (ny-1)) 01075 boundary_info.add_side(elem, 2, 3); 01076 01077 if (k==0) 01078 boundary_info.add_side(elem, 0, 0); 01079 01080 if (k == (nz-1)) 01081 boundary_info.add_side(elem, 4, 5); 01082 } 01083 break; 01084 } 01085 01086 01087 01088 01089 01090 01091 case HEX20: 01092 case HEX27: 01093 case TET4: // TET4's are created from an initial HEX27 discretization 01094 case TET10: // TET10's are created from an initial HEX27 discretization 01095 case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization 01096 case PYRAMID13: 01097 case PYRAMID14: 01098 { 01099 for (unsigned int k=0; k<(2*nz); k += 2) 01100 for (unsigned int j=0; j<(2*ny); j += 2) 01101 for (unsigned int i=0; i<(2*nx); i += 2) 01102 { 01103 Elem* elem = (type == HEX20) ? 01104 mesh.add_elem(new Hex20) : 01105 mesh.add_elem(new Hex27); 01106 01107 elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) ); 01108 elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ); 01109 elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ); 01110 elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ); 01111 elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2)); 01112 elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)); 01113 elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2)); 01114 elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)); 01115 elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ); 01116 elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ); 01117 elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ); 01118 elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ); 01119 elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1)); 01120 elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)); 01121 elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)); 01122 elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)); 01123 elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)); 01124 elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)); 01125 elem->set_node(18) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)); 01126 elem->set_node(19) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)); 01127 if ((type == HEX27) || (type == TET4) || (type == TET10) || 01128 (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14)) 01129 { 01130 elem->set_node(20) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ); 01131 elem->set_node(21) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)); 01132 elem->set_node(22) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)); 01133 elem->set_node(23) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)); 01134 elem->set_node(24) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)); 01135 elem->set_node(25) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)); 01136 elem->set_node(26) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)); 01137 } 01138 01139 01140 if (k == 0) 01141 boundary_info.add_side(elem, 0, 0); 01142 01143 if (k == 2*(nz-1)) 01144 boundary_info.add_side(elem, 5, 5); 01145 01146 if (j == 0) 01147 boundary_info.add_side(elem, 1, 1); 01148 01149 if (j == 2*(ny-1)) 01150 boundary_info.add_side(elem, 3, 3); 01151 01152 if (i == 0) 01153 boundary_info.add_side(elem, 4, 4); 01154 01155 if (i == 2*(nx-1)) 01156 boundary_info.add_side(elem, 2, 2); 01157 } 01158 break; 01159 } 01160 01161 01162 01163 01164 case PRISM15: 01165 case PRISM18: 01166 { 01167 for (unsigned int k=0; k<(2*nz); k += 2) 01168 for (unsigned int j=0; j<(2*ny); j += 2) 01169 for (unsigned int i=0; i<(2*nx); i += 2) 01170 { 01171 // First Prism 01172 Elem* elem = NULL; 01173 elem = ((type == PRISM15) ? 01174 mesh.add_elem(new Prism15) : 01175 mesh.add_elem(new Prism18)); 01176 01177 elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) ); 01178 elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ); 01179 elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ); 01180 elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2)); 01181 elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)); 01182 elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)); 01183 elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ); 01184 elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ); 01185 elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ); 01186 elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1)); 01187 elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)); 01188 elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)); 01189 elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)); 01190 elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)); 01191 elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)); 01192 if (type == PRISM18) 01193 { 01194 elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)); 01195 elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)); 01196 elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)); 01197 } 01198 01199 // Add sides for first prism to boundary info object 01200 if (i==0) 01201 boundary_info.add_side(elem, 3, 4); 01202 01203 if (j==0) 01204 boundary_info.add_side(elem, 1, 1); 01205 01206 if (k==0) 01207 boundary_info.add_side(elem, 0, 0); 01208 01209 if (k == 2*(nz-1)) 01210 boundary_info.add_side(elem, 4, 5); 01211 01212 01213 // Second Prism 01214 elem = ((type == PRISM15) ? 01215 mesh.add_elem(new Prism15) : 01216 mesh.add_elem(new Prism18)); 01217 01218 elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k) ); 01219 elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ); 01220 elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k) ); 01221 elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) ); 01222 elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) ); 01223 elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) ); 01224 elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ); 01225 elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ); 01226 elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ); 01227 elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) ); 01228 elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)); 01229 elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) ); 01230 elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)); 01231 elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)); 01232 elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)); 01233 if (type == PRISM18) 01234 { 01235 elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)); 01236 elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)); 01237 elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)); 01238 } 01239 01240 // Add sides for second prism to boundary info object 01241 if (i == 2*(nx-1)) 01242 boundary_info.add_side(elem, 1, 2); 01243 01244 if (j == 2*(ny-1)) 01245 boundary_info.add_side(elem, 2, 3); 01246 01247 if (k==0) 01248 boundary_info.add_side(elem, 0, 0); 01249 01250 if (k == 2*(nz-1)) 01251 boundary_info.add_side(elem, 4, 5); 01252 01253 } 01254 break; 01255 } 01256 01257 01258 01259 01260 01261 default: 01262 libmesh_error_msg("ERROR: Unrecognized 3D element type."); 01263 } 01264 01265 01266 01267 01268 //....................................... 01269 // Scale the nodal positions 01270 for (unsigned int p=0; p<mesh.n_nodes(); p++) 01271 { 01272 mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin; 01273 mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin; 01274 mesh.node(p)(2) = (mesh.node(p)(2))*(zmax-zmin) + zmin; 01275 } 01276 01277 01278 01279 01280 // Additional work for tets and pyramids: we take the existing 01281 // HEX27 discretization and split each element into 24 01282 // sub-tets or 6 sub-pyramids. 01283 // 01284 // 24 isn't the minimum-possible number of tets, but it 01285 // obviates any concerns about the edge orientations between 01286 // the various elements. 01287 if ((type == TET4) || 01288 (type == TET10) || 01289 (type == PYRAMID5) || 01290 (type == PYRAMID13) || 01291 (type == PYRAMID14)) 01292 { 01293 // Temporary storage for new elements. (24 tets per hex, 6 pyramids) 01294 std::vector<Elem*> new_elements; 01295 01296 if ((type == TET4) || (type == TET10)) 01297 new_elements.reserve(24*mesh.n_elem()); 01298 else 01299 new_elements.reserve(6*mesh.n_elem()); 01300 01301 // Create tetrahedra or pyramids 01302 { 01303 MeshBase::element_iterator el = mesh.elements_begin(); 01304 const MeshBase::element_iterator end_el = mesh.elements_end(); 01305 01306 for ( ; el != end_el; ++el) 01307 { 01308 // Get a pointer to the HEX27 element. 01309 Elem* base_hex = *el; 01310 01311 // Get a pointer to the node located at the HEX27 centroid 01312 Node* apex_node = base_hex->get_node(26); 01313 01314 for (unsigned short s=0; s<base_hex->n_sides(); ++s) 01315 { 01316 // Get the boundary ID for this side 01317 boundary_id_type b_id = boundary_info.boundary_id(*el, s); 01318 01319 // Need to build the full-ordered side! 01320 UniquePtr<Elem> side = base_hex->build_side(s); 01321 01322 if ((type == TET4) || (type == TET10)) 01323 { 01324 // Build 4 sub-tets per side 01325 for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet) 01326 { 01327 new_elements.push_back( new Tet4 ); 01328 Elem* sub_elem = new_elements.back(); 01329 sub_elem->set_node(0) = side->get_node(sub_tet); 01330 sub_elem->set_node(1) = side->get_node(8); // centroid of the face 01331 sub_elem->set_node(2) = side->get_node(sub_tet==3 ? 0 : sub_tet+1 ); // wrap-around 01332 sub_elem->set_node(3) = apex_node; // apex node always used! 01333 01334 // If the original hex was a boundary hex, add the new sub_tet's side 01335 // 0 with the same b_id. Note: the tets are all aligned so that their 01336 // side 0 is on the boundary. 01337 if (b_id != BoundaryInfo::invalid_id) 01338 boundary_info.add_side(sub_elem, 0, b_id); 01339 } 01340 } // end if ((type == TET4) || (type == TET10)) 01341 01342 else // type==PYRAMID5 || type==PYRAMID13 || type==PYRAMID14 01343 { 01344 // Build 1 sub-pyramid per side. 01345 new_elements.push_back(new Pyramid5); 01346 Elem* sub_elem = new_elements.back(); 01347 01348 // Set the base. Note that since the apex is *inside* the base_hex, 01349 // and the pyramid uses a counter-clockwise base numbering, we need to 01350 // reverse the [1] and [3] node indices. 01351 sub_elem->set_node(0) = side->get_node(0); 01352 sub_elem->set_node(1) = side->get_node(3); 01353 sub_elem->set_node(2) = side->get_node(2); 01354 sub_elem->set_node(3) = side->get_node(1); 01355 01356 // Set the apex 01357 sub_elem->set_node(4) = apex_node; 01358 01359 // If the original hex was a boundary hex, add the new sub_pyr's side 01360 // 4 (the square base) with the same b_id. 01361 if (b_id != BoundaryInfo::invalid_id) 01362 boundary_info.add_side(sub_elem, 4, b_id); 01363 } // end else type==PYRAMID5 || type==PYRAMID13 || type==PYRAMID14 01364 } 01365 } 01366 } 01367 01368 01369 // Delete the original HEX27 elements from the mesh, and the boundary info structure. 01370 { 01371 MeshBase::element_iterator el = mesh.elements_begin(); 01372 const MeshBase::element_iterator end_el = mesh.elements_end(); 01373 01374 for ( ; el != end_el; ++el) 01375 { 01376 boundary_info.remove(*el); // Safe even if *el has no boundary info. 01377 mesh.delete_elem(*el); 01378 } 01379 } 01380 01381 // Add the new elements 01382 for (unsigned int i=0; i<new_elements.size(); ++i) 01383 mesh.add_elem(new_elements[i]); 01384 01385 } // end if (type == TET4,TET10,PYRAMID5,PYRAMID13,PYRAMID14 01386 01387 01388 // Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's 01389 if ((type == TET10) || (type == PYRAMID14)) 01390 mesh.all_second_order(); 01391 01392 else if (type == PYRAMID13) 01393 mesh.all_second_order(/*full_ordered=*/false); 01394 01395 01396 // Add sideset names to boundary info (Z axis out of the screen) 01397 boundary_info.sideset_name(0) = "back"; 01398 boundary_info.sideset_name(1) = "bottom"; 01399 boundary_info.sideset_name(2) = "right"; 01400 boundary_info.sideset_name(3) = "top"; 01401 boundary_info.sideset_name(4) = "left"; 01402 boundary_info.sideset_name(5) = "front"; 01403 01404 // Add nodeset names to boundary info 01405 boundary_info.nodeset_name(0) = "back"; 01406 boundary_info.nodeset_name(1) = "bottom"; 01407 boundary_info.nodeset_name(2) = "right"; 01408 boundary_info.nodeset_name(3) = "top"; 01409 boundary_info.nodeset_name(4) = "left"; 01410 boundary_info.nodeset_name(5) = "front"; 01411 01412 break; 01413 } // end case dim==3 01414 01415 default: 01416 libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension()); 01417 } 01418 01419 STOP_LOG("build_cube()", "MeshTools::Generation"); 01420 01421 01422 01423 // Done building the mesh. Now prepare it for use. 01424 mesh.prepare_for_use (/*skip_renumber =*/ false); 01425 } 01426 01427 01428 01429 void MeshTools::Generation::build_point (UnstructuredMesh& mesh, 01430 const ElemType type, 01431 const bool gauss_lobatto_grid) 01432 { 01433 // This method only makes sense in 0D! 01434 // But we now just turn a non-0D mesh into a 0D mesh 01435 //libmesh_assert_equal_to (mesh.mesh_dimension(), 1); 01436 01437 build_cube(mesh, 01438 0, 0, 0, 01439 0., 0., 01440 0., 0., 01441 0., 0., 01442 type, 01443 gauss_lobatto_grid); 01444 } 01445 01446 01447 void MeshTools::Generation::build_line (UnstructuredMesh& mesh, 01448 const unsigned int nx, 01449 const Real xmin, const Real xmax, 01450 const ElemType type, 01451 const bool gauss_lobatto_grid) 01452 { 01453 // This method only makes sense in 1D! 01454 // But we now just turn a non-1D mesh into a 1D mesh 01455 //libmesh_assert_equal_to (mesh.mesh_dimension(), 1); 01456 01457 build_cube(mesh, 01458 nx, 0, 0, 01459 xmin, xmax, 01460 0., 0., 01461 0., 0., 01462 type, 01463 gauss_lobatto_grid); 01464 } 01465 01466 01467 01468 void MeshTools::Generation::build_square (UnstructuredMesh& mesh, 01469 const unsigned int nx, 01470 const unsigned int ny, 01471 const Real xmin, const Real xmax, 01472 const Real ymin, const Real ymax, 01473 const ElemType type, 01474 const bool gauss_lobatto_grid) 01475 { 01476 // This method only makes sense in 2D! 01477 // But we now just turn a non-2D mesh into a 2D mesh 01478 //libmesh_assert_equal_to (mesh.mesh_dimension(), 2); 01479 01480 // Call the build_cube() member to actually do the work for us. 01481 build_cube (mesh, 01482 nx, ny, 0, 01483 xmin, xmax, 01484 ymin, ymax, 01485 0., 0., 01486 type, 01487 gauss_lobatto_grid); 01488 } 01489 01490 01491 01492 01493 01494 01495 01496 01497 01498 #ifndef LIBMESH_ENABLE_AMR 01499 void MeshTools::Generation::build_sphere (UnstructuredMesh&, 01500 const Real, 01501 const unsigned int, 01502 const ElemType, 01503 const unsigned int, 01504 const bool) 01505 { 01506 libmesh_error_msg("Building a circle/sphere only works with AMR."); 01507 } 01508 01509 #else 01510 01511 void MeshTools::Generation::build_sphere (UnstructuredMesh& mesh, 01512 const Real rad, 01513 const unsigned int nr, 01514 const ElemType type, 01515 const unsigned int n_smooth, 01516 const bool flat) 01517 { 01518 libmesh_assert_greater (rad, 0.); 01519 //libmesh_assert_greater (nr, 0); // must refine at least once otherwise will end up with a square/cube 01520 01521 START_LOG("build_sphere()", "MeshTools::Generation"); 01522 01523 // Clear the mesh and start from scratch 01524 mesh.clear(); 01525 01526 BoundaryInfo& boundary_info = mesh.get_boundary_info(); 01527 01528 // Sphere is centered at origin by default 01529 const Point cent; 01530 01531 const Sphere sphere (cent, rad); 01532 01533 switch (mesh.mesh_dimension()) 01534 { 01535 //----------------------------------------------------------------- 01536 // Build a line in one dimension 01537 case 1: 01538 { 01539 build_line (mesh, 3, -rad, rad, type); 01540 } 01541 01542 01543 01544 01545 //----------------------------------------------------------------- 01546 // Build a circle or hollow sphere in two dimensions 01547 case 2: 01548 { 01549 // For ParallelMesh, if we don't specify node IDs the Mesh 01550 // will try to pick an appropriate (unique) one for us. But 01551 // since we are adding these nodes on all processors, we want 01552 // to be sure they have consistent IDs across all processors. 01553 unsigned node_id = 0; 01554 01555 if (flat) 01556 { 01557 const Real sqrt_2 = std::sqrt(2.); 01558 const Real rad_2 = .25*rad; 01559 const Real rad_sqrt_2 = rad/sqrt_2; 01560 01561 // (Temporary) convenient storage for node pointers 01562 std::vector<Node*> nodes(8); 01563 01564 // Point 0 01565 nodes[0] = mesh.add_point (Point(-rad_2,-rad_2, 0.), node_id++); 01566 01567 // Point 1 01568 nodes[1] = mesh.add_point (Point( rad_2,-rad_2, 0.), node_id++); 01569 01570 // Point 2 01571 nodes[2] = mesh.add_point (Point( rad_2, rad_2, 0.), node_id++); 01572 01573 // Point 3 01574 nodes[3] = mesh.add_point (Point(-rad_2, rad_2, 0.), node_id++); 01575 01576 // Point 4 01577 nodes[4] = mesh.add_point (Point(-rad_sqrt_2,-rad_sqrt_2, 0.), node_id++); 01578 01579 // Point 5 01580 nodes[5] = mesh.add_point (Point( rad_sqrt_2,-rad_sqrt_2, 0.), node_id++); 01581 01582 // Point 6 01583 nodes[6] = mesh.add_point (Point( rad_sqrt_2, rad_sqrt_2, 0.), node_id++); 01584 01585 // Point 7 01586 nodes[7] = mesh.add_point (Point(-rad_sqrt_2, rad_sqrt_2, 0.), node_id++); 01587 01588 // Build the elements & set node pointers 01589 01590 // Element 0 01591 { 01592 Elem* elem0 = mesh.add_elem (new Quad4); 01593 elem0->set_node(0) = nodes[0]; 01594 elem0->set_node(1) = nodes[1]; 01595 elem0->set_node(2) = nodes[2]; 01596 elem0->set_node(3) = nodes[3]; 01597 } 01598 01599 // Element 1 01600 { 01601 Elem* elem1 = mesh.add_elem (new Quad4); 01602 elem1->set_node(0) = nodes[4]; 01603 elem1->set_node(1) = nodes[0]; 01604 elem1->set_node(2) = nodes[3]; 01605 elem1->set_node(3) = nodes[7]; 01606 } 01607 01608 // Element 2 01609 { 01610 Elem* elem2 = mesh.add_elem (new Quad4); 01611 elem2->set_node(0) = nodes[4]; 01612 elem2->set_node(1) = nodes[5]; 01613 elem2->set_node(2) = nodes[1]; 01614 elem2->set_node(3) = nodes[0]; 01615 } 01616 01617 // Element 3 01618 { 01619 Elem* elem3 = mesh.add_elem (new Quad4); 01620 elem3->set_node(0) = nodes[1]; 01621 elem3->set_node(1) = nodes[5]; 01622 elem3->set_node(2) = nodes[6]; 01623 elem3->set_node(3) = nodes[2]; 01624 } 01625 01626 // Element 4 01627 { 01628 Elem* elem4 = mesh.add_elem (new Quad4); 01629 elem4->set_node(0) = nodes[3]; 01630 elem4->set_node(1) = nodes[2]; 01631 elem4->set_node(2) = nodes[6]; 01632 elem4->set_node(3) = nodes[7]; 01633 } 01634 01635 } 01636 else 01637 { 01638 // Create the 12 vertices of a regular unit icosahedron 01639 Real t = 0.5 * (1 + std::sqrt(5.0)); 01640 Real s = rad / std::sqrt(1 + t*t); 01641 t *= s; 01642 01643 mesh.add_point (Point(-s, t, 0), node_id++); 01644 mesh.add_point (Point( s, t, 0), node_id++); 01645 mesh.add_point (Point(-s, -t, 0), node_id++); 01646 mesh.add_point (Point( s, -t, 0), node_id++); 01647 01648 mesh.add_point (Point( 0, -s, t), node_id++); 01649 mesh.add_point (Point( 0, s, t), node_id++); 01650 mesh.add_point (Point( 0, -s, -t), node_id++); 01651 mesh.add_point (Point( 0, s, -t), node_id++); 01652 01653 mesh.add_point (Point( t, 0, -s), node_id++); 01654 mesh.add_point (Point( t, 0, s), node_id++); 01655 mesh.add_point (Point(-t, 0, -s), node_id++); 01656 mesh.add_point (Point(-t, 0, s), node_id++); 01657 01658 // Create the 20 triangles of the icosahedron 01659 static const unsigned int idx1 [6] = {11, 5, 1, 7, 10, 11}; 01660 static const unsigned int idx2 [6] = {9, 4, 2, 6, 8, 9}; 01661 static const unsigned int idx3 [6] = {1, 5, 11, 10, 7, 1}; 01662 01663 for (unsigned int i = 0; i < 5; ++i) 01664 { 01665 // 5 elems around point 0 01666 Elem* new_elem = mesh.add_elem (new Tri3); 01667 new_elem->set_node(0) = mesh.node_ptr(0); 01668 new_elem->set_node(1) = mesh.node_ptr(idx1[i]); 01669 new_elem->set_node(2) = mesh.node_ptr(idx1[i+1]); 01670 01671 // 5 adjacent elems 01672 new_elem = mesh.add_elem (new Tri3); 01673 new_elem->set_node(0) = mesh.node_ptr(idx3[i]); 01674 new_elem->set_node(1) = mesh.node_ptr(idx3[i+1]); 01675 new_elem->set_node(2) = mesh.node_ptr(idx2[i]); 01676 01677 // 5 elems around point 3 01678 new_elem = mesh.add_elem (new Tri3); 01679 new_elem->set_node(0) = mesh.node_ptr(3); 01680 new_elem->set_node(1) = mesh.node_ptr(idx2[i]); 01681 new_elem->set_node(2) = mesh.node_ptr(idx2[i+1]); 01682 01683 // 5 adjacent elems 01684 new_elem = mesh.add_elem (new Tri3); 01685 new_elem->set_node(0) = mesh.node_ptr(idx2[i+1]); 01686 new_elem->set_node(1) = mesh.node_ptr(idx2[i]); 01687 new_elem->set_node(2) = mesh.node_ptr(idx3[i+1]); 01688 } 01689 } 01690 01691 break; 01692 } // end case 2 01693 01694 01695 01696 01697 01698 //----------------------------------------------------------------- 01699 // Build a sphere in three dimensions 01700 case 3: 01701 { 01702 // (Currently) supported types 01703 if (!((type == HEX8) || (type == HEX27))) 01704 { 01705 // FIXME: We'd need an all_tet() routine (which could also be used by 01706 // build_square()) to do Tets, or Prisms for that matter. 01707 libmesh_error_msg("Error: Only HEX8/27 currently supported."); 01708 } 01709 01710 01711 // 3D analog of 2D initial grid: 01712 const Real 01713 r_small = 0.25*rad, // 0.25 *radius 01714 r_med = (0.125*std::sqrt(2.)+0.5)*rad; // .67677*radius 01715 01716 // (Temporary) convenient storage for node pointers 01717 std::vector<Node*> nodes(16); 01718 01719 // For ParallelMesh, if we don't specify node IDs the Mesh 01720 // will try to pick an appropriate (unique) one for us. But 01721 // since we are adding these nodes on all processors, we want 01722 // to be sure they have consistent IDs across all processors. 01723 unsigned node_id = 0; 01724 01725 // Points 0-7 are the initial HEX8 01726 nodes[0] = mesh.add_point (Point(-r_small,-r_small, -r_small), node_id++); 01727 nodes[1] = mesh.add_point (Point( r_small,-r_small, -r_small), node_id++); 01728 nodes[2] = mesh.add_point (Point( r_small, r_small, -r_small), node_id++); 01729 nodes[3] = mesh.add_point (Point(-r_small, r_small, -r_small), node_id++); 01730 nodes[4] = mesh.add_point (Point(-r_small,-r_small, r_small), node_id++); 01731 nodes[5] = mesh.add_point (Point( r_small,-r_small, r_small), node_id++); 01732 nodes[6] = mesh.add_point (Point( r_small, r_small, r_small), node_id++); 01733 nodes[7] = mesh.add_point (Point(-r_small, r_small, r_small), node_id++); 01734 01735 // Points 8-15 are for the outer hexes, we number them in the same way 01736 nodes[8] = mesh.add_point (Point(-r_med,-r_med, -r_med), node_id++); 01737 nodes[9] = mesh.add_point (Point( r_med,-r_med, -r_med), node_id++); 01738 nodes[10] = mesh.add_point (Point( r_med, r_med, -r_med), node_id++); 01739 nodes[11] = mesh.add_point (Point(-r_med, r_med, -r_med), node_id++); 01740 nodes[12] = mesh.add_point (Point(-r_med,-r_med, r_med), node_id++); 01741 nodes[13] = mesh.add_point (Point( r_med,-r_med, r_med), node_id++); 01742 nodes[14] = mesh.add_point (Point( r_med, r_med, r_med), node_id++); 01743 nodes[15] = mesh.add_point (Point(-r_med, r_med, r_med), node_id++); 01744 01745 // Now create the elements and add them to the mesh 01746 // Element 0 - center element 01747 { 01748 Elem* elem0 = mesh.add_elem (new Hex8); 01749 elem0->set_node(0) = nodes[0]; 01750 elem0->set_node(1) = nodes[1]; 01751 elem0->set_node(2) = nodes[2]; 01752 elem0->set_node(3) = nodes[3]; 01753 elem0->set_node(4) = nodes[4]; 01754 elem0->set_node(5) = nodes[5]; 01755 elem0->set_node(6) = nodes[6]; 01756 elem0->set_node(7) = nodes[7]; 01757 } 01758 01759 // Element 1 - "bottom" 01760 { 01761 Elem* elem1 = mesh.add_elem (new Hex8); 01762 elem1->set_node(0) = nodes[8]; 01763 elem1->set_node(1) = nodes[9]; 01764 elem1->set_node(2) = nodes[10]; 01765 elem1->set_node(3) = nodes[11]; 01766 elem1->set_node(4) = nodes[0]; 01767 elem1->set_node(5) = nodes[1]; 01768 elem1->set_node(6) = nodes[2]; 01769 elem1->set_node(7) = nodes[3]; 01770 } 01771 01772 // Element 2 - "front" 01773 { 01774 Elem* elem2 = mesh.add_elem (new Hex8); 01775 elem2->set_node(0) = nodes[8]; 01776 elem2->set_node(1) = nodes[9]; 01777 elem2->set_node(2) = nodes[1]; 01778 elem2->set_node(3) = nodes[0]; 01779 elem2->set_node(4) = nodes[12]; 01780 elem2->set_node(5) = nodes[13]; 01781 elem2->set_node(6) = nodes[5]; 01782 elem2->set_node(7) = nodes[4]; 01783 } 01784 01785 // Element 3 - "right" 01786 { 01787 Elem* elem3 = mesh.add_elem (new Hex8); 01788 elem3->set_node(0) = nodes[1]; 01789 elem3->set_node(1) = nodes[9]; 01790 elem3->set_node(2) = nodes[10]; 01791 elem3->set_node(3) = nodes[2]; 01792 elem3->set_node(4) = nodes[5]; 01793 elem3->set_node(5) = nodes[13]; 01794 elem3->set_node(6) = nodes[14]; 01795 elem3->set_node(7) = nodes[6]; 01796 } 01797 01798 // Element 4 - "back" 01799 { 01800 Elem* elem4 = mesh.add_elem (new Hex8); 01801 elem4->set_node(0) = nodes[3]; 01802 elem4->set_node(1) = nodes[2]; 01803 elem4->set_node(2) = nodes[10]; 01804 elem4->set_node(3) = nodes[11]; 01805 elem4->set_node(4) = nodes[7]; 01806 elem4->set_node(5) = nodes[6]; 01807 elem4->set_node(6) = nodes[14]; 01808 elem4->set_node(7) = nodes[15]; 01809 } 01810 01811 // Element 5 - "left" 01812 { 01813 Elem* elem5 = mesh.add_elem (new Hex8); 01814 elem5->set_node(0) = nodes[8]; 01815 elem5->set_node(1) = nodes[0]; 01816 elem5->set_node(2) = nodes[3]; 01817 elem5->set_node(3) = nodes[11]; 01818 elem5->set_node(4) = nodes[12]; 01819 elem5->set_node(5) = nodes[4]; 01820 elem5->set_node(6) = nodes[7]; 01821 elem5->set_node(7) = nodes[15]; 01822 } 01823 01824 // Element 6 - "top" 01825 { 01826 Elem* elem6 = mesh.add_elem (new Hex8); 01827 elem6->set_node(0) = nodes[4]; 01828 elem6->set_node(1) = nodes[5]; 01829 elem6->set_node(2) = nodes[6]; 01830 elem6->set_node(3) = nodes[7]; 01831 elem6->set_node(4) = nodes[12]; 01832 elem6->set_node(5) = nodes[13]; 01833 elem6->set_node(6) = nodes[14]; 01834 elem6->set_node(7) = nodes[15]; 01835 } 01836 01837 break; 01838 } // end case 3 01839 01840 default: 01841 libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension()); 01842 01843 01844 01845 } // end switch (dim) 01846 01847 // Now we have the beginnings of a sphere. 01848 // Add some more elements by doing uniform refinements and 01849 // popping nodes to the boundary. 01850 MeshRefinement mesh_refinement (mesh); 01851 01852 // Loop over the elements, refine, pop nodes to boundary. 01853 for (unsigned int r=0; r<nr; r++) 01854 { 01855 mesh_refinement.uniformly_refine(1); 01856 01857 MeshBase::element_iterator it = mesh.active_elements_begin(); 01858 const MeshBase::element_iterator end = mesh.active_elements_end(); 01859 01860 for (; it != end; ++it) 01861 { 01862 Elem* elem = *it; 01863 01864 for (unsigned int s=0; s<elem->n_sides(); s++) 01865 if (elem->neighbor(s) == NULL || (mesh.mesh_dimension() == 2 && !flat)) 01866 { 01867 UniquePtr<Elem> side(elem->build_side(s)); 01868 01869 // Pop each point to the sphere boundary 01870 for (unsigned int n=0; n<side->n_nodes(); n++) 01871 side->point(n) = 01872 sphere.closest_point(side->point(n)); 01873 } 01874 } 01875 } 01876 01877 // The mesh now contains a refinement hierarchy due to the refinements 01878 // used to generate the grid. In order to call other support functions 01879 // like all_tri() and all_second_order, you need a "flat" mesh file (with no 01880 // refinement trees) so 01881 MeshTools::Modification::flatten(mesh); 01882 01883 // In 2D, convert all the quads to triangles if requested 01884 if (mesh.mesh_dimension()==2) 01885 { 01886 if ((type == TRI6) || (type == TRI3)) 01887 { 01888 MeshTools::Modification::all_tri(mesh); 01889 } 01890 } 01891 01892 01893 // Convert to second-order elements if the user requested it. 01894 if (Elem::second_order_equivalent_type(type) == INVALID_ELEM) 01895 { 01896 // type is already a second-order, determine if it is the 01897 // "full-ordered" second-order element, or the "serendipity" 01898 // second order element. Note also that all_second_order 01899 // can't be called once the mesh has been refined. 01900 bool full_ordered = !((type==QUAD8) || (type==HEX20)); 01901 mesh.all_second_order(full_ordered); 01902 01903 // And pop to the boundary again... 01904 MeshBase::element_iterator it = mesh.active_elements_begin(); 01905 const MeshBase::element_iterator end = mesh.active_elements_end(); 01906 01907 for (; it != end; ++it) 01908 { 01909 Elem* elem = *it; 01910 01911 for (unsigned int s=0; s<elem->n_sides(); s++) 01912 if (elem->neighbor(s) == NULL) 01913 { 01914 UniquePtr<Elem> side(elem->build_side(s)); 01915 01916 // Pop each point to the sphere boundary 01917 for (unsigned int n=0; n<side->n_nodes(); n++) 01918 side->point(n) = 01919 sphere.closest_point(side->point(n)); 01920 } 01921 } 01922 } 01923 01924 01925 // The meshes could probably use some smoothing. 01926 LaplaceMeshSmoother smoother(mesh); 01927 smoother.smooth(n_smooth); 01928 01929 // We'll give the whole sphere surface a boundary id of 0 01930 { 01931 MeshBase::element_iterator it = mesh.active_elements_begin(); 01932 const MeshBase::element_iterator end = mesh.active_elements_end(); 01933 01934 for (; it != end; ++it) 01935 { 01936 Elem* elem = *it; 01937 for (unsigned short s=0; s != elem->n_sides(); ++s) 01938 if (!elem->neighbor(s)) 01939 boundary_info.add_side(elem, s, 0); 01940 } 01941 } 01942 01943 STOP_LOG("build_sphere()", "MeshTools::Generation"); 01944 01945 01946 // Done building the mesh. Now prepare it for use. 01947 mesh.prepare_for_use(/*skip_renumber =*/ false); 01948 } 01949 01950 #endif // #ifndef LIBMESH_ENABLE_AMR 01951 01952 01953 // Meshes the tensor product of a 1D and a 1D-or-2D domain. 01954 void MeshTools::Generation::build_extrusion (UnstructuredMesh& mesh, 01955 const MeshBase& cross_section, 01956 const unsigned int nz, 01957 RealVectorValue extrusion_vector) 01958 { 01959 if (!cross_section.n_elem()) 01960 return; 01961 01962 START_LOG("build_extrusion()", "MeshTools::Generation"); 01963 01964 dof_id_type orig_elem = cross_section.n_elem(); 01965 dof_id_type orig_nodes = cross_section.n_nodes(); 01966 01967 unsigned int order = 1; 01968 01969 BoundaryInfo& boundary_info = mesh.get_boundary_info(); 01970 const BoundaryInfo& cross_section_boundary_info = cross_section.get_boundary_info(); 01971 01972 // If cross_section is distributed, so is its extrusion 01973 if (!cross_section.is_serial()) 01974 mesh.delete_remote_elements(); 01975 01976 // We know a priori how many elements we'll need 01977 mesh.reserve_elem(nz*orig_elem); 01978 01979 // For straightforward meshes we need one or two additional layers per 01980 // element. 01981 if ((*cross_section.elements_begin())->default_order() == SECOND) 01982 order = 2; 01983 01984 mesh.reserve_nodes((order*nz+1)*orig_nodes); 01985 01986 MeshBase::const_node_iterator nd = cross_section.nodes_begin(); 01987 const MeshBase::const_node_iterator nend = cross_section.nodes_end(); 01988 for (; nd!=nend; ++nd) 01989 { 01990 const Node* node = *nd; 01991 01992 for (unsigned int k=0; k != order*nz+1; ++k) 01993 { 01994 Node *new_node = 01995 mesh.add_point(*node + 01996 (extrusion_vector * k / nz / order), 01997 node->id() + (k * orig_nodes), 01998 node->processor_id()); 01999 02000 const std::vector<boundary_id_type> ids_to_copy = 02001 cross_section_boundary_info.boundary_ids(node); 02002 02003 boundary_info.add_node(new_node, ids_to_copy); 02004 } 02005 } 02006 02007 const std::set<boundary_id_type> &side_ids = 02008 cross_section_boundary_info.get_side_boundary_ids(); 02009 const boundary_id_type next_side_id = side_ids.empty() ? 02010 0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1); 02011 02012 MeshBase::const_element_iterator el = cross_section.elements_begin(); 02013 const MeshBase::const_element_iterator end = cross_section.elements_end(); 02014 for (; el!=end; ++el) 02015 { 02016 const Elem* elem = *el; 02017 const ElemType etype = elem->type(); 02018 02019 // build_extrusion currently only works on coarse meshes 02020 libmesh_assert (!elem->parent()); 02021 02022 for (unsigned int k=0; k != nz; ++k) 02023 { 02024 Elem *new_elem; 02025 switch (etype) 02026 { 02027 case EDGE2: 02028 { 02029 new_elem = new Quad4; 02030 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes)); 02031 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes)); 02032 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes)); 02033 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes)); 02034 break; 02035 } 02036 case EDGE3: 02037 { 02038 new_elem = new Quad9; 02039 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes)); 02040 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes)); 02041 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes)); 02042 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes)); 02043 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes)); 02044 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes)); 02045 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes)); 02046 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes)); 02047 new_elem->set_node(8) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes)); 02048 break; 02049 } 02050 case TRI3: 02051 { 02052 new_elem = new Prism6; 02053 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes)); 02054 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes)); 02055 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes)); 02056 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes)); 02057 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes)); 02058 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes)); 02059 break; 02060 } 02061 case TRI6: 02062 { 02063 new_elem = new Prism18; 02064 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes)); 02065 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes)); 02066 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes)); 02067 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes)); 02068 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes)); 02069 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes)); 02070 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes)); 02071 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes)); 02072 new_elem->set_node(8) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes)); 02073 new_elem->set_node(9) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes)); 02074 new_elem->set_node(10) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes)); 02075 new_elem->set_node(11) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes)); 02076 new_elem->set_node(12) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes)); 02077 new_elem->set_node(13) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes)); 02078 new_elem->set_node(14) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes)); 02079 new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes)); 02080 new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes)); 02081 new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes)); 02082 break; 02083 } 02084 case QUAD4: 02085 { 02086 new_elem = new Hex8; 02087 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes)); 02088 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes)); 02089 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes)); 02090 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (k * orig_nodes)); 02091 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes)); 02092 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes)); 02093 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes)); 02094 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((k+1) * orig_nodes)); 02095 break; 02096 } 02097 case QUAD9: 02098 { 02099 new_elem = new Hex27; 02100 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes)); 02101 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes)); 02102 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes)); 02103 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes)); 02104 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes)); 02105 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes)); 02106 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes)); 02107 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes)); 02108 new_elem->set_node(8) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes)); 02109 new_elem->set_node(9) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes)); 02110 new_elem->set_node(10) = mesh.node_ptr(elem->get_node(6)->id() + (2*k * orig_nodes)); 02111 new_elem->set_node(11) = mesh.node_ptr(elem->get_node(7)->id() + (2*k * orig_nodes)); 02112 new_elem->set_node(12) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes)); 02113 new_elem->set_node(13) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes)); 02114 new_elem->set_node(14) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes)); 02115 new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes)); 02116 new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes)); 02117 new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes)); 02118 new_elem->set_node(18) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+2) * orig_nodes)); 02119 new_elem->set_node(19) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+2) * orig_nodes)); 02120 new_elem->set_node(20) = mesh.node_ptr(elem->get_node(8)->id() + (2*k * orig_nodes)); 02121 new_elem->set_node(21) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes)); 02122 new_elem->set_node(22) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes)); 02123 new_elem->set_node(23) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+1) * orig_nodes)); 02124 new_elem->set_node(24) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+1) * orig_nodes)); 02125 new_elem->set_node(25) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+2) * orig_nodes)); 02126 new_elem->set_node(26) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+1) * orig_nodes)); 02127 break; 02128 } 02129 default: 02130 { 02131 libmesh_not_implemented(); 02132 break; 02133 } 02134 } 02135 02136 new_elem->set_id(elem->id() + (k * orig_elem)); 02137 new_elem->processor_id() = elem->processor_id(); 02138 02139 // maintain the subdomain_id 02140 new_elem->subdomain_id() = elem->subdomain_id(); 02141 02142 new_elem = mesh.add_elem(new_elem); 02143 02144 // Copy any old boundary ids on all sides 02145 for (unsigned short s = 0; s != elem->n_sides(); ++s) 02146 { 02147 const std::vector<boundary_id_type> ids_to_copy = 02148 cross_section_boundary_info.boundary_ids(elem, s); 02149 02150 if (new_elem->dim() == 3) 02151 { 02152 // For 2D->3D extrusion, we give the boundary IDs 02153 // for side s on the old element to side s+1 on the 02154 // new element. This is just a happy coincidence as 02155 // far as I can tell... 02156 boundary_info.add_side 02157 (new_elem, cast_int<unsigned short>(s+1), 02158 ids_to_copy); 02159 } 02160 else 02161 { 02162 // For 1D->2D extrusion, the boundary IDs map as: 02163 // Old elem -> New elem 02164 // 0 -> 3 02165 // 1 -> 1 02166 libmesh_assert_less(s, 2); 02167 const unsigned short sidemap[2] = {3, 1}; 02168 boundary_info.add_side(new_elem, sidemap[s], ids_to_copy); 02169 } 02170 } 02171 02172 // Give new boundary ids to bottom and top 02173 if (k == 0) 02174 boundary_info.add_side(new_elem, 0, next_side_id); 02175 if (k == nz-1) 02176 { 02177 // For 2D->3D extrusion, the "top" ID is 1+the original 02178 // element's number of sides. For 1D->2D extrusion, the 02179 // "top" ID is side 2. 02180 const unsigned short top_id = new_elem->dim() == 3 ? 02181 cast_int<unsigned short>(elem->n_sides()+1) : 2; 02182 boundary_info.add_side 02183 (new_elem, top_id, 02184 cast_int<boundary_id_type>(next_side_id+1)); 02185 } 02186 } 02187 } 02188 02189 STOP_LOG("build_extrusion()", "MeshTools::Generation"); 02190 02191 // Done building the mesh. Now prepare it for use. 02192 mesh.prepare_for_use(/*skip_renumber =*/ false); 02193 } 02194 02195 02196 02197 02198 #ifdef LIBMESH_HAVE_TRIANGLE 02199 02200 // Triangulates a 2D rectangular region with or without holes 02201 void MeshTools::Generation::build_delaunay_square(UnstructuredMesh& mesh, 02202 const unsigned int nx, // num. of elements in x-dir 02203 const unsigned int ny, // num. of elements in y-dir 02204 const Real xmin, const Real xmax, 02205 const Real ymin, const Real ymax, 02206 const ElemType type, 02207 const std::vector<TriangleInterface::Hole*>* holes) 02208 { 02209 // Check for reasonable size 02210 libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction 02211 libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction 02212 libmesh_assert_less (xmin, xmax); 02213 libmesh_assert_less (ymin, ymax); 02214 02215 // Clear out any data which may have been in the Mesh 02216 mesh.clear(); 02217 02218 BoundaryInfo& boundary_info = mesh.get_boundary_info(); 02219 02220 // Make sure the new Mesh will be 2D 02221 mesh.set_mesh_dimension(2); 02222 02223 // The x and y spacing between boundary points 02224 const Real delta_x = (xmax-xmin) / static_cast<Real>(nx); 02225 const Real delta_y = (ymax-ymin) / static_cast<Real>(ny); 02226 02227 // Bottom 02228 for (unsigned int p=0; p<=nx; ++p) 02229 mesh.add_point(Point(xmin + p*delta_x, ymin)); 02230 02231 // Right side 02232 for (unsigned int p=1; p<ny; ++p) 02233 mesh.add_point(Point(xmax, ymin + p*delta_y)); 02234 02235 // Top 02236 for (unsigned int p=0; p<=nx; ++p) 02237 mesh.add_point(Point(xmax - p*delta_x, ymax)); 02238 02239 // Left side 02240 for (unsigned int p=1; p<ny; ++p) 02241 mesh.add_point(Point(xmin, ymax - p*delta_y)); 02242 02243 // Be sure we added as many points as we thought we did 02244 libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny)); 02245 02246 // Construct the Triangle Interface object 02247 TriangleInterface t(mesh); 02248 02249 // Set custom variables for the triangulation 02250 t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny); 02251 t.triangulation_type() = TriangleInterface::PSLG; 02252 t.elem_type() = type; 02253 02254 if (holes != NULL) 02255 t.attach_hole_list(holes); 02256 02257 // Triangulate! 02258 t.triangulate(); 02259 02260 // The mesh is now generated, but we still need to mark the boundaries 02261 // to be consistent with the other build_square routines. Note that all 02262 // hole boundary elements get the same ID, 4. 02263 MeshBase::element_iterator el = mesh.elements_begin(); 02264 const MeshBase::element_iterator end_el = mesh.elements_end(); 02265 for ( ; el != end_el; ++el) 02266 { 02267 const Elem* elem = *el; 02268 02269 for (unsigned int s=0; s<elem->n_sides(); s++) 02270 if (elem->neighbor(s) == NULL) 02271 { 02272 UniquePtr<Elem> side (elem->build_side(s)); 02273 02274 // Check the location of the side's midpoint. Since 02275 // the square has straight sides, the midpoint is not 02276 // on the corner and thus it is uniquely on one of the 02277 // sides. 02278 Point side_midpoint= 0.5f*( (*side->get_node(0)) + (*side->get_node(1)) ); 02279 02280 // The boundary ids are set following the same convention as Quad4 sides 02281 // bottom = 0 02282 // right = 1 02283 // top = 2 02284 // left = 3 02285 // hole = 4 02286 boundary_id_type bc_id=4; 02287 02288 // bottom 02289 if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE) 02290 bc_id=0; 02291 02292 // right 02293 else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE) 02294 bc_id=1; 02295 02296 // top 02297 else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE) 02298 bc_id=2; 02299 02300 // left 02301 else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE) 02302 bc_id=3; 02303 02304 // If the point is not on any of the external boundaries, it 02305 // is on one of the holes.... 02306 02307 // Finally, add this element's information to the boundary info object. 02308 boundary_info.add_side(elem->id(), s, bc_id); 02309 } 02310 } 02311 02312 } // end build_delaunay_square 02313 02314 #endif // LIBMESH_HAVE_TRIANGLE 02315 02316 02317 02318 } // namespace libMesh
generated by
