mesh_generation.c

一个用来实现偏微分方程中网格的计算库

// $Id: mesh_generation.C 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// C++ includes
#include <cmath> // for std::sqrt


// Local includes
#include "mesh_generation.h"
#include "unstructured_mesh.h"
// #include "elem.h"
#include "mesh_refinement.h"
#include "edge_edge2.h"
#include "edge_edge3.h"
#include "edge_edge4.h"
#include "face_tri3.h"
#include "face_tri6.h"
#include "face_quad4.h"
#include "face_quad8.h"
#include "face_quad9.h"
#include "cell_hex8.h"
#include "cell_hex20.h"
#include "cell_hex27.h"
#include "cell_prism6.h"
#include "cell_prism15.h"
#include "cell_prism18.h"
#include "cell_tet4.h"
#include "libmesh_logging.h"
#include "boundary_info.h"
#include "sphere.h"
#include "mesh_modification.h"
#include "mesh_smoother_laplace.h"

// ------------------------------------------------------------
// MeshTools::Generation function for mesh generation
void MeshTools::Generation::build_cube(UnstructuredMesh& mesh,
				       const unsigned int nx,
				       const unsigned int ny,
				       const unsigned int nz,
				       const Real xmin, const Real xmax,
				       const Real ymin, const Real ymax,
				       const Real zmin, const Real zmax,
				       const ElemType type,
				       const bool gauss_lobatto_grid)
{
  START_LOG("build_cube()", "MeshTools::Generation");

  // Declare that we are using the indexing utility routine
  // in the "Private" part of our current namespace.  If this doesn't
  // work in GCC 2.95.3 we can either remove it or stop supporting
  // 2.95.3 altogether.
  using MeshTools::Generation::Private::idx;
  
  // Clear the mesh and start from scratch
  mesh.clear();
  
  switch (mesh.mesh_dimension())
    {

      
      //---------------------------------------------------------------------
      // Build a 1D line
    case 1:      
      {
	libmesh_assert (nx != 0);

        // Reserve elements
        switch (type)
          {
          case INVALID_ELEM:
          case EDGE2:
          case EDGE3:
          case EDGE4:
            {
	      mesh.reserve_elem (nx);
              break;
            }

          default:
            {
	      std::cerr << "ERROR: Unrecognized 1D element type." << std::endl;
	      libmesh_error();
	    }
	  }

        // Reserve nodes
        switch (type)
          {
          case INVALID_ELEM:
          case EDGE2:
            {
              mesh.reserve_nodes(nx+1);
              break;
            }

          case EDGE3:
            {
              mesh.reserve_nodes(2*nx+1);
              break;
            }

          case EDGE4:
            {
              mesh.reserve_nodes(3*nx+1);
              break;
            }

          default:
            {
              std::cerr << "ERROR: Unrecognized 1D element type." << std::endl;
              libmesh_error();
            }
          }


        // Build the nodes, depends on whether we're using linears, 
        // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
        unsigned int node_id = 0;
        switch(type)
        {
          case INVALID_ELEM:
          case EDGE2:
            {
              for (unsigned int i=0; i<=nx; i++)
              {
                if (gauss_lobatto_grid)
                  mesh.add_point (Point(0.5*(std::cos(libMesh::pi*static_cast<Real>(nx-i)/static_cast<Real>(nx))+1.0),
                        0, 
                        0), node_id++);
                else
                  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx), 
                        0, 
                        0), node_id++);
              }
              break;
            }

          case EDGE3:
            {
              for (unsigned int i=0; i<=2*nx; i++)
              {
                if (gauss_lobatto_grid)
                {
                  // The x location of the point.
                  Real x=0.;

                  // Shortcut variable
                  const Real pi = libMesh::pi;

                  // Shortcut quantities (do not depend on i)
                  const Real c = std::cos( pi*i / static_cast<Real>(2*nx) );

                  // If i is even, compute a normal Gauss-Lobatto point
                  if (i%2 == 0)
                    x = 0.5*(1.0 - c);

                  // Otherwise, it is the average of the previous and next points
                  else
                  {
                    Real cmin = std::cos( pi*(i-1) / static_cast<Real>(2*nx) );
                    Real cmax = std::cos( pi*(i+1) / static_cast<Real>(2*nx) );
                    
                    Real xmin = 0.5*(1.0 - cmin);
                    Real xmax = 0.5*(1.0 - cmax);
                    x = 0.5*(xmin + xmax);
                  }

                  mesh.add_point (Point(x,0.,0.), node_id++);
                }
                else
                  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
                        0,
                        0), node_id++);
              }
              break;
            }

          case EDGE4:
            {
              for (unsigned int i=0; i<=3*nx; i++)
              {
                if (gauss_lobatto_grid)
                {
                  // The x location of the point
                  Real x=0.;
                  
                  const Real pi = libMesh::pi;

                  // Shortcut quantities
                  const Real c = std::cos( pi*i / static_cast<Real>(3*nx) );

                  // If i is multiple of 3, compute a normal Gauss-Lobatto point
                  if (i%3 == 0)
                    x = 0.5*(1.0 - c);

                  // Otherwise, distribute points evenly within the element
                  else
                  {
                    if(i%3 == 1)
                    {
                      Real cmin = std::cos( pi*(i-1) / static_cast<Real>(3*nx) );
                      Real cmax = std::cos( pi*(i+2) / static_cast<Real>(3*nx) );

                      Real xmin = 0.5*(1.0 - cmin);
                      Real xmax = 0.5*(1.0 - cmax);

                      x = (2.*xmin + xmax)/3.;
                    }
                    else
                    if(i%3 == 2)
                    {
                      Real cmin = std::cos( pi*(i-2) / static_cast<Real>(3*nx) );
                      Real cmax = std::cos( pi*(i+1) / static_cast<Real>(3*nx) );

                      Real xmin = 0.5*(1.0 - cmin);
                      Real xmax = 0.5*(1.0 - cmax);

                      x = (xmin + 2.*xmax)/3.;
                    }

                  }

                  mesh.add_point (Point(x,0.,0.), node_id++);
                }
                else
                mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(3*nx),
                        0,
                        0), node_id++);
              }


            
              break;
            }
            
          default:
            {
              std::cerr << "ERROR: Unrecognized 1D element type." << std::endl;
              libmesh_error();
            }
              
        }

        // Build the elements of the mesh
        switch(type)
          {
            case INVALID_ELEM:
            case EDGE2:
              {
                for (unsigned int i=0; i<nx; i++)
                {
                  Elem* elem = mesh.add_elem (new Edge2);
                  elem->set_node(0) = mesh.node_ptr(i);
                  elem->set_node(1) = mesh.node_ptr(i+1);
                }
              break;
              }

            case EDGE3:
              {
                for (unsigned int i=0; i<nx; i++)
                {
                  Elem* elem = mesh.add_elem (new Edge3);
                  elem->set_node(0) = mesh.node_ptr(2*i);
                  elem->set_node(2) = mesh.node_ptr(2*i+1);
                  elem->set_node(1) = mesh.node_ptr(2*i+2);
                }
              break;
              }

            case EDGE4:
              {
                for (unsigned int i=0; i<nx; i++)
                {
                  Elem* elem = mesh.add_elem (new Edge4);
                  elem->set_node(0) = mesh.node_ptr(3*i);
                  elem->set_node(2) = mesh.node_ptr(3*i+1);
                  elem->set_node(3) = mesh.node_ptr(3*i+2);
                  elem->set_node(1) = mesh.node_ptr(3*i+3);
                }
              break;
              }

            default:
              {
                std::cerr << "ERROR: Unrecognized 1D element type." << std::endl;
                libmesh_error();                
              }
          }


	// Scale the nodal positions
	for (unsigned int p=0; p<mesh.n_nodes(); p++)
	  mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
	
	break;
      }







      
      
      
      //---------------------------------------------------------------------
      // Build a 2D quadrilateral
    case 2:
      {
	libmesh_assert (nx != 0);
	libmesh_assert (ny != 0);

	// Reserve elements.  The TRI3 and TRI6 meshes
	// have twice as many elements...
	switch (type)
	  {
	  case INVALID_ELEM:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    {
	      mesh.reserve_elem (nx*ny);
	      break;
	    }

	  case TRI3:
	  case TRI6:
	    {
	      mesh.reserve_elem (2*nx*ny);
	      break;
	    }

	  default:
	    {
	      std::cerr << "ERROR: Unrecognized 2D element type." << std::endl;
	      libmesh_error();
	    }
	  }


	
	// Reserve nodes.  The quadratic element types
	// need to reserve more nodes than the linear types.
	switch (type)
	  {
	  case INVALID_ELEM:
	  case QUAD4:
	  case TRI3:
	    {
	      mesh.reserve_nodes( (nx+1)*(ny+1) );
	      break;
	    }
	    
	  case QUAD8:
	  case QUAD9:
	  case TRI6:
	    {
	      mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
	      break;
	    }

	    
	  default:
	    {
	      std::cerr << "ERROR: Unrecognized 2D element type." << std::endl;
	      libmesh_error();
	    }
	  }



	// Build the nodes. Depends on whether you are using a linear
	// or quadratic element, and whether you are using a uniform
	// grid or the Gauss-Lobatto grid points.
        unsigned int node_id = 0;
	switch (type)
	  {
	  case INVALID_ELEM:
	  case QUAD4:
	  case TRI3:
	    {
	      for (unsigned int j=0; j<=ny; j++)
		for (unsigned int i=0; i<=nx; i++)
		  {
		    if (gauss_lobatto_grid)
		      {
			// Shortcut variable
			const Real pi = libMesh::pi;

			mesh.add_point (Point(0.5*(1.0 - std::cos(pi*static_cast<Real>(i)/static_cast<Real>(nx))),
					      0.5*(1.0 - std::cos(pi*static_cast<Real>(j)/static_cast<Real>(ny))),
					      0.), node_id++);