meshtool.cc

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

#include "libmesh_config.h"

// C++ includes
#include <iostream>
#include <vector>
#include <string>
#include <unistd.h>
#ifdef HAVE_GETOPT_H
// GCC 2.95.3 (and maybe others) do not include
// getopt.h in unistd.h...  Hower IBM xlC has no
// getopt.h!  This works around that.
#include <getopt.h>
#endif
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <fstream>
#include <algorithm>

// Local Includes
#include "libmesh.h"
#include "elem.h"
#include "mesh.h"
#include "mesh_modification.h"
#include "boundary_mesh.h"
#include "mesh_refinement.h"
#include "dof_map.h"
#include "perfmon.h"
#include "point.h"
#include "elem_quality.h"
#include "gmv_io.h"
#include "legacy_xdr_io.h"
#include "statistics.h"
#include "inf_elem_builder.h"
#include "mesh_data.h"
#include "boundary_info.h"



/*
 * convenient enum for the mode in which the boundary mesh
 * should be written
 */
enum BoundaryMeshWriteMode {BM_DISABLED=0, BM_MESH_ONLY, BM_WITH_MESHDATA};


void usage(char *progName)
{
  std::string baseName;
  static std::string helpList =
    "usage:\n"
    "        %s [options] ...\n"
    "\n"
    "options:\n"
    "    -i <string>                   Input file name\n"
    "    -o <string>                   Output file name\n"
    "    -s <string>                   Solution file name\n"
    "    -d <dim>                      <dim>-dimensional mesh\n"
  "\n    -D <factor>                   Randomly move interior nodes by D*hmin\n"
#ifdef ENABLE_AMR
    "    -r <count>                    Globally refine <count> times\n"
#endif
    "    -p <count>                    Partition into <count> subdomains\n"
    "    -t (-d 2 only)                Convert to triangles first\n"
    "    -b                            Write the boundary conditions\n"
    "    -B                            Like -b, but with activated MeshData\n"
    "                                  (allows to write .unv file of the\n"
    "                                  boundary with the correct node ids)\n"
    "    -2                            Converts a mesh of linear elements\n"
    "                                  to their second-order counterparts:\n"
    "                                  Quad4 -> Quad8, Tet4 -> Tet10 etc\n"
    "    -3                            Same, but to the highest possible:\n"
    "                                  Quad4 -> Quad9, Hex8 -> Hex27 etc\n"
#ifdef ENABLE_INFINITE_ELEMENTS
  "\n    -a                            Add infinite elements\n"
    "    -x <coord>                    Specify infinite element origin\n"
    "    -y <coord>                    coordinates. If none given, origin\n"
    "    -z <coord>                    is determined automatically.\n"
    "    -X                            When building infinite elements \n"
    "    -Y                            treat mesh as x/y/z-symmetric.\n"
    "    -Z                            When -X is given, -x <coord> also\n"
    "                                  has to be given.  Similar for y,z.\n"
#endif
/*  "\n    -l                            Build the L connectivity matrix \n"         */
/*    "    -L                            Build the script L connectivity matrix \n"  */
    "    -v                            Verbose\n"
    "    -h                            Print help menu\n"
    "\n"
    "\n"
    " This program is used to convert and partions from/to a variety of\n"
    " formats.  File types are inferred from file extensions.  For example,\n"
    " the command:\n"
    "\n"
    "  ./meshtool -d 2 -i in.exd -o out.plt\n"
    "\n"
    " will read a 2D mesh in the ExodusII format (from Cubit, for example)\n"
    " from the file in.exd.  It will then write the mesh in the Tecplot\n"
    " binary format to out.plt.\n"
    "\n"
    " and\n"
    "\n"
    "  ./meshtool -d 3 -i bench12.mesh.0000 -o out.gmv -s bench12.soln.0137\n"
    "\n"
    " will read a 3D MGF mesh from the file bench12.mesh.0000, read a\n"
    " solution from bench12.soln.0137, and write the output in GMV format\n"
    " to out.gmv\n"
#ifdef ENABLE_INFINITE_ELEMENTS
    "\n"
    " and\n"
    "\n"
    "  ./meshtool -d 3 -i dry.unv -o packed.gmv -a -x 30.5 -y -10.5 -X\n"
    "\n"
    " will read a 3D Universal file, determine the z-coordinate of the origin\n"
    " automatically, e.g. z_origin = 3., build infinite elements with the\n"
    " origin (30.5, -10.5, 3.) on top of volume elements, while preserving\n"
    " a symmetry plane through (30.5, -10.5, 3.) perpendicular to x.\n"
    " It is imperative that the origin lies _inside_ the given volume mesh.\n"
    " If not, infinite elements are not correctly built!\n"
#endif
    "\n"
    " Currently this program supports the following formats:\n"
    "\n"
    "INPUT:\n"
    "     .ucd -- AVS unstructured ASCII grid format\n"
    "     .exd -- Sandia's ExodusII binary grid format\n"
    "     .unv -- SDRC I-Deas Universal File ASCII format\n"
    "     .xdr -- Internal binary mesh format\n"
    "     .xda -- Same format, but ASCII and human-readable\n"
    "     .mgf -- MGF binary mesh format\n"
    "\n"
    "OUTPUT:\n"
    "     .plt   -- Tecplot binary format\n"
    "     .dat   -- Tecplot ASCII format\n"
    "     .gmv   -- LANL's General Mesh Viewer (~benkirk/work/GMV/linuxogl)\n"
    "     .ugrid -- Kelly's DIVA ASCII format (3D only)\n"
    "     .xdr   -- Internal binary mesh format\n"
    "     .xda   -- Same format, but ASCII and human-readable\n"
    "     .mgf   -- MGF binary mesh format\n"
    "     .fro   -- ACDL's .fro format\n"
    "\n"
    " Direct questions to:\n"
    " [e-mail]\n";

    
  if (progName == NULL)
    baseName = "UNKNOWN";
  else 
    baseName = progName;

  
  fprintf(stderr, helpList.c_str(), baseName.c_str());
  fflush(stderr);

  abort();
}



void process_cmd_line(int argc, char **argv,
                      std::vector<std::string>& names,
                      unsigned int& n_subdomains,
                      unsigned int& n_rsteps,
                      unsigned int& dim,
                      double& dist_fact,
                      bool& verbose,
                      BoundaryMeshWriteMode& write_bndry,
                      unsigned int& convert_second_order,
                      bool& triangulate,
                      bool& addinfelems,

#ifdef ENABLE_INFINITE_ELEMENTS
                      InfElemBuilder::InfElemOriginValue& origin_x,
                      InfElemBuilder::InfElemOriginValue& origin_y,
                      InfElemBuilder::InfElemOriginValue& origin_z,
#endif

                      bool& x_sym,
                      bool& y_sym,
                      bool& z_sym
                      )
{

#ifndef ENABLE_INFINITE_ELEMENTS

  /*
   * initialize these to some values,
   * so that the compiler does not complain
   */
  triangulate     = false;
  addinfelems     = false;
  x_sym           = y_sym           = z_sym           = false;

  char optionStr[] =
    "i:o:s:d:D:r:p:tbB23vlLm?h";

#else

  char optionStr[] =
    "i:o:s:d:D:r:p:tbB23a::x:y:z:XYZvlLm?h";

#endif

  bool b_mesh_b_given = false;
  bool b_mesh_B_given = false;

  int opt;  
  
  if (argc < 3)
    usage(argv[0]);


  
  while ((opt = getopt(argc, argv, optionStr)) != -1)
    {
      switch (opt)
        {
          
          /**
           * Get input file name
           */
        case 'i':
          {
            if (names.empty())
              names.push_back(optarg);
            else
              {
                std::cout << "ERROR: Input name must preceed output name!"
                          << std::endl;
                exit(1);
              }
            break;
          }
          
          /**
           * Get output file name
           */
        case 'o':
          {
            if (!names.empty())
              names.push_back(optarg);
            else
              {
                std::cout << "ERROR: Input name must preceed output name!"
                          << std::endl;
                exit(1);
              }
            break;
          }

          /**
           * Get solution file name
           */
        case 's':
          {
            if (names.size() == 2)
              names.push_back(optarg);
            else
              {
                std::cout << "ERROR: Input and output names must preceed solution name!"
                          << std::endl;
                exit(1);
              }
            break;
          }

          /**
           * Get the mesh dimension
           */
        case 'd':
          {
            dim = atoi(optarg);
            break;
          }

          /**
           * Get the mesh distortion factor
           */
        case 'D':
          {
            dist_fact = atof(optarg);
            break;
          }

          /**
           * Get the number of refinements to do
           */
        case 'r':
          {
            n_rsteps = atoi(optarg);
            break;
          }

          /**
           * Get the number of subdomains for partitioning
           */
        case 'p':
          {
            n_subdomains = atoi(optarg);
            break;
          }
          
          /**
           * Triangulate in 2D
           */
        case 't':
          {
            triangulate = true;
            break;
          }
          
          /**
           * Be verbose
           */
        case 'v':
          {
            verbose = true;
            break;
          }
          
          /**
           * Try to write the boundary
           */
        case 'b':
          {
            if (b_mesh_B_given)
              {
                std::cout << "ERROR: Do not use -b and -B concurrently!"
                          << std::endl;
                exit(1);
              }
                  
            b_mesh_b_given = true;
            write_bndry = BM_MESH_ONLY;
            break;
          }
          
          /**
           * Try to write the boundary,
           * but copy also the MeshData
           */
        case 'B':
          {
            if (b_mesh_b_given)
              {
                std::cout << "ERROR: Do not use -b and -B concurrently!"
                          << std::endl;
                exit(1);
              }
                  
            b_mesh_B_given = true;
            write_bndry = BM_WITH_MESHDATA;
            break;
          }
          

          /**
           * Convert elements to second-order
           * counterparts
           */
        case '2':
          {
            convert_second_order = 2;
            break;
          }

          /**
           * Convert elements to second-order
           * counterparts, highest order possible
           */
        case '3':
          {
            convert_second_order = 22;
            break;
          }


#ifdef ENABLE_INFINITE_ELEMENTS

          /**
           * Add infinite elements
           */
        case 'a':
          {
            addinfelems = true;
            break;
          }
          
          /**
           * Specify origin coordinates
           */          
        case 'x':
          {
            origin_x.first  = true;
            origin_x.second = atof(optarg);
            break;
          }
          
        case 'y':
          {
            origin_y.first  = true;
            origin_y.second = atof(optarg);
            break;
          }
          
        case 'z':
          {
            origin_z.first  = true;
            origin_z.second = atof(optarg);
            break;
          } 
          
          /**
           *  Symmetries
           */                    
        case 'X':
          {
            x_sym = true;
            break;
          }
          case 'Y':
          {
            y_sym = true;
            break;
          }
          case 'Z':
          {
            z_sym = true;
            break;
          }

#endif //ifdef ENABLE_INFINITE_ELEMENTS



        case 'h':
        case '?':
          usage(argv[0]);
          
        default:
          return;
        };
    };

}



// A helper function for creating a submesh of active elements.  Why do we need this?
// Well, the create submesh function is expecting const_element_iterators,
// and this is just one way to get them...
void construct_mesh_of_active_elements(Mesh& new_mesh, const Mesh& mesh)
{
  MeshBase::const_element_iterator       it     = mesh.active_elements_begin();
  const MeshBase::const_element_iterator it_end = mesh.active_elements_end();
  mesh.create_submesh(new_mesh, it, it_end);
}