meshtool.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// 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
 
#include "libmesh/libmesh_config.h"
 
// C++ includes
#include <algorithm>
#include <fstream>
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <string>
#include <unistd.h>
#include <vector>
 
// Local Includes
#include "libmesh/boundary_info.h"
#include "libmesh/boundary_mesh.h"
#include "libmesh/dof_map.h"
#include "libmesh/elem.h"
#include "libmesh/elem_quality.h"
#include "libmesh/gmv_io.h"
#include "libmesh/inf_elem_builder.h"
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/statistics.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/enum_elem_quality.h"
#include "libmesh/getpot.h"
 
using namespace libMesh;
 
// convenient enum for the mode in which the boundary mesh
// should be written
enum BoundaryMeshWriteMode {BM_DISABLED=0, BM_MESH_ONLY};
 
// Print a usage message before exiting the program.
void usage(const std::string & prog_name);
 
int main (int argc, char ** argv)
{
  LibMeshInit init(argc, argv);
 
  unsigned int n_subdomains = 1;
  unsigned int n_rsteps = 0;
  double dist_fact = 0.;
  bool verbose = false;
  BoundaryMeshWriteMode write_bndry = BM_DISABLED;
  bool convert_first_order = false;
  unsigned int convert_second_order = 0;
  bool triangulate = false;
  bool do_quality = false;
  ElemQuality quality_type = DIAGONAL;
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
  bool addinfelems = false;
  InfElemBuilder::InfElemOriginValue origin_x(false, 0.);
  InfElemBuilder::InfElemOriginValue origin_y(false, 0.);
  InfElemBuilder::InfElemOriginValue origin_z(false, 0.);
  bool x_sym=false;
  bool y_sym=false;
  bool z_sym=false;
#endif
 
  std::vector<std::string> names;
  std::vector<std::string> var_names;
  std::vector<Number>      soln;
 
  // Check for minimum number of command line args
  if (argc < 3)
    usage(std::string(argv[0]));
 
  // Create a GetPot object to parse the command line
  GetPot command_line (argc, argv);
 
  // Print usage and exit immediately if user asked for help.
  if (command_line.search(2, "-h", "-?"))
    usage(argv[0]);
 
  // Input file name
  if (command_line.search(1, "-i"))
    {
      std::string tmp;
      tmp = command_line.next(tmp);
      if (names.empty())
        names.push_back(tmp);
      else
        libmesh_error_msg("ERROR: Input name must precede output name!");
    }
 
  // Output file name
  if (command_line.search(1, "-o"))
    {
      std::string tmp;
      tmp = command_line.next(tmp);
      if (!names.empty())
        names.push_back(tmp);
      else
        libmesh_error_msg("ERROR: Input name must precede output name!");
    }
 
  // Get the mesh distortion factor
  if (command_line.search(1, "-D"))
    dist_fact = command_line.next(dist_fact);
 
  // Number of refinements
  if (command_line.search(1, "-r"))
    {
      int tmp;
      tmp = command_line.next(tmp);
      n_rsteps = cast_int<unsigned int>(tmp);
    }
 
  // Number of subdomains for partitioning
  if (command_line.search(1, "-p"))
    {
      int tmp;
      tmp = command_line.next(tmp);
      n_subdomains = cast_int<unsigned int>(tmp);
    }
 
  // Should we call all_tri()?
  if (command_line.search(1, "-t"))
    triangulate = true;
 
  // Should we calculate element quality?
  if (command_line.search(1, "-q"))
    {
      do_quality = true;
      std::string tmp;
      tmp = command_line.next(tmp);
      if (tmp != "")
        quality_type = Utility::string_to_enum<ElemQuality>(tmp);
    }
 
  // Should we be verbose?
  if (command_line.search(1, "-v"))
    verbose = true;
 
  // Should we write the boundary?
  if (command_line.search(1, "-b"))
    write_bndry = BM_MESH_ONLY;
 
  // Should we convert all elements to 1st order?
  if (command_line.search(1, "-1"))
    convert_first_order = true;
 
  // Should we convert all elements to 2nd order?
  if (command_line.search(1, "-2"))
    convert_second_order = 2;
 
  // Should we convert all elements to "full" 2nd order?
  if (command_line.search(1, "-3"))
    convert_second_order = 22;
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  // Add infinte elements
  if (command_line.search(1, "-a"))
    addinfelems = true;
 
  // Specify origin x coordinate
  if (command_line.search(1, "-x"))
    {
      origin_x.first  = true;
      origin_x.second = command_line.next(origin_x.second);
    }
 
  // Specify origin y coordinate
  if (command_line.search(1, "-y"))
    {
      origin_y.first  = true;
      origin_y.second = command_line.next(origin_y.second);
    }
 
  // Specify origin z coordinate
  if (command_line.search(1, "-z"))
    {
      origin_z.first  = true;
      origin_z.second = command_line.next(origin_z.second);
    }
 
  // Symmetries
  if (command_line.search(1, "-X"))
    x_sym = true;
  if (command_line.search(1, "-Y"))
    y_sym = true;
  if (command_line.search(1, "-Z"))
    z_sym = true;
 
#endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  // Read the input mesh
  Mesh mesh(init.comm());
  if (!names.empty())
    {
      mesh.read(names[0]);
 
      if (verbose)
        {
          mesh.print_info();
          mesh.get_boundary_info().print_summary();
        }
    }
 
  else
    {
      libMesh::out << "No input specified." << std::endl;
      return 1;
    }
 
 
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  if (addinfelems)
    {
      if (names.size() == 3)
        libmesh_error_msg("ERROR: Invalid combination: Building infinite elements\n"
                          << "not compatible with solution import.");
 
      if (write_bndry != BM_DISABLED)
        libmesh_error_msg("ERROR: Invalid combination: Building infinite elements\n"
                          << "not compatible with writing boundary conditions.");
 
      // Sanity checks: -X/Y/Z can only be used, when the
      // corresponding coordinate is also given (using -x/y/z)
      if ((x_sym && !origin_x.first) ||     // claim x-symmetry, but x-coordinate of origin not given!
          (y_sym && !origin_y.first) ||     // the same for y
          (z_sym && !origin_z.first))       // the same for z
        libmesh_error_msg("ERROR: When x-symmetry is requested using -X, then\n"
                          << "the option -x <coord> also has to be given.\n"
                          << "This holds obviously for y and z, too.");
 
      // build infinite elements
      InfElemBuilder(mesh).build_inf_elem(origin_x, origin_y, origin_z,
                                          x_sym, y_sym, z_sym,
                                          verbose);
 
      if (verbose)
        {
          mesh.print_info();
          mesh.get_boundary_info().print_summary();
        }
 
    }
 
  // sanity check
  else if ((origin_x.first ||  origin_y.first || origin_z.first) ||
           (x_sym          ||  y_sym          || z_sym))
    libmesh_error_msg("ERROR:  -x/-y/-z/-X/-Y/-Z is only to be used when\n"
                      << "the option -a is also specified!");
 
#endif
 
 
  // Maybe Triangulate
  if (triangulate)
    {
      if (verbose)
        libMesh::out << "...Converting to all simplices...\n";
 
      MeshTools::Modification::all_tri(mesh);
    }
 
  // Compute Shape quality metrics
  if (do_quality)
    {
      StatisticsVector<Real> sv;
      sv.reserve(mesh.n_elem());
 
      libMesh::out << "Quality type is: " << Quality::name(quality_type) << std::endl;
 
      // What are the quality bounds for this element?
      std::pair<Real, Real> bounds = mesh.elem_ref(0).qual_bounds(quality_type);
      libMesh::out << "Quality bounds for this element type are: ("
                   << bounds.first
                   << ", "
                   << bounds.second
                   << ") "
                   << std::endl;
 
      for (const auto & elem : mesh.active_element_ptr_range())
        sv.push_back(elem->quality(quality_type));
 
      const unsigned int n_bins = 10;
      libMesh::out << "Avg. shape quality: " << sv.mean() << std::endl;
 
      // Find element indices below the specified cutoff.
      // These might be considered "bad" elements which need refinement.
      std::vector<dof_id_type> bad_elts  = sv.cut_below(0.8);
      libMesh::out << "Found " << bad_elts.size()
                   << " of " << mesh.n_elem()
                   << " elements below the cutoff." << std::endl;
 
      // Compute the histogram for this distribution
      std::vector<dof_id_type> histogram;
      sv.histogram(histogram, n_bins);
 
      const bool do_matlab = true;
 
      if (do_matlab)
        {
          std::ofstream out ("histo.m");
 
          out << "% This is a sample histogram plot for Matlab." << std::endl;
          out << "bin_members = [" << std::endl;
          for (unsigned int i=0; i<n_bins; i++)
            out << static_cast<Real>(histogram[i]) / static_cast<Real>(mesh.n_elem())
                << std::endl;
          out << "];" << std::endl;
 
          std::vector<Real> bin_coords(n_bins);
          const Real max   = *(std::max_element(sv.begin(), sv.end()));
          const Real min   = *(std::min_element(sv.begin(), sv.end()));
          const Real delta = (max - min) / static_cast<Real>(n_bins);
          for (unsigned int i=0; i<n_bins; i++)
            bin_coords[i] = min + (i * delta) + delta / 2.0 ;
 
          out << "bin_coords = [" << std::endl;
          for (unsigned int i=0; i<n_bins; i++)
            out << bin_coords[i] << std::endl;
          out << "];" << std::endl;
 
          out << "bar(bin_coords, bin_members, 1);" << std::endl;
          out << "hold on" << std::endl;
          out << "plot (bin_coords, 0, 'kx');" << std::endl;
          out << "xlabel('Quality (0=Worst, 1=Best)');" << std::endl;
          out << "ylabel('Percentage of elements in each bin');" << std::endl;
          out << "axis([" << min << "," << max << ",0, max(bin_members)]);" << std::endl;
 
          out << "title('" << Quality::name(quality_type) << "');" << std::endl;
 
        }
    }
 
   // Possibly convert all linear elements to second-order
   // counterparts
  if (convert_first_order)
    {
      if (verbose)
        libMesh::out << "Converting elements to first order counterparts\n";
 
      mesh.all_first_order();
 
      if (verbose)
        {
          mesh.print_info();
          mesh.get_boundary_info().print_summary();
        }
    }
 
  // Possibly convert all linear elements to second-order counterparts
  if (convert_second_order > 0)
    {
      bool second_order_mode = true;
      std:: string message = "Converting elements to second order counterparts";
      if (convert_second_order == 2)
        {
          second_order_mode = false;
          message += ", lower version: Quad4 -> Quad8, not Quad9";
        }
 
      else if (convert_second_order == 22)
        {
          second_order_mode = true;
          message += ", highest version: Quad4 -> Quad9";
        }
 
      else
        libmesh_error_msg("Invalid value, convert_second_order = " << convert_second_order);
 
      if (verbose)
        libMesh::out << message << std::endl;
 
      mesh.all_second_order(second_order_mode);
 
      if (verbose)
        {
          mesh.print_info();
          mesh.get_boundary_info().print_summary();
        }
    }
 
 
#ifdef LIBMESH_ENABLE_AMR
 
  // Possibly refine the mesh
  if (n_rsteps > 0)
    {
      if (verbose)
        libMesh::out << "Refining the mesh "
                     << n_rsteps << " times"
                     << std::endl;
 
      MeshRefinement mesh_refinement (mesh);
      mesh_refinement.uniformly_refine(n_rsteps);
 
      if (verbose)
        {
          mesh.print_info();
          mesh.get_boundary_info().print_summary();
        }
    }
 
 
  // Possibly distort the mesh
  if (dist_fact > 0.)
    {
      libMesh::out << "Distorting the mesh by a factor of "
                   << dist_fact
                   << std::endl;
 
      MeshTools::Modification::distort(mesh,dist_fact);
    }
 
#endif
 
 
 
  // Possibly partition the mesh
  if (n_subdomains > 1)
    mesh.partition(n_subdomains);
 
 
  // Possibly write the mesh
  if (names.size() >= 2)
    {
      // When the mesh got refined, it is likely that
      // the user does _not_ want to write also
      // the coarse elements, but only the active ones.
      // Use Mesh::create_submesh() to create a mesh
      // of only active elements, and then write _this_
      // new mesh.
      if (n_rsteps > 0)
        {
          if (verbose)
            libMesh::out << " Mesh got refined, will write only _active_ elements." << std::endl;
 
          Mesh new_mesh (init.comm(), cast_int<unsigned char>(mesh.mesh_dimension()));
 
          mesh.create_submesh
            (new_mesh,
             mesh.active_elements_begin(),
             mesh.active_elements_end());
 
          // now write the new_mesh
          if (names.size() == 2)
            new_mesh.write(names[1]);
          else if (names.size() == 3)
            new_mesh.write(names[1], soln, var_names);
          else
            libmesh_error_msg("Invalid names.size() = " << names.size());
        }
      else
        {
          if (names.size() == 2)
            mesh.write(names[1]);
          else if (names.size() == 3)
            mesh.write(names[1], soln, var_names);
          else
            libmesh_error_msg("Invalid names.size() = " << names.size());
        }
 
 
 
      // Possibly write the BCs
      if (write_bndry != BM_DISABLED)
        {
          BoundaryMesh boundary_mesh
            (mesh.comm(), cast_int<unsigned char>(mesh.mesh_dimension()-1));
 
          std::string boundary_name = "bndry_";
          boundary_name += names[1];
 
          if (write_bndry == BM_MESH_ONLY)
            mesh.get_boundary_info().sync(boundary_mesh);
 
          else
            libmesh_error_msg("Invalid value write_bndry = " << write_bndry);
 
          if (names.size() == 2)
            boundary_mesh.write(boundary_name);
          else if (names.size() == 3)
            boundary_mesh.write(boundary_name, soln, var_names);
        }
    }
 
  return 0;
}
 
 
void usage(const std::string & prog_name)
{
  std::ostringstream helpList;
  helpList << "usage:\n"
           << "        "
           << prog_name
           << " [options] ...\n"
           << "\n"
           << "options:\n"
           << "    -d <dim>                      <dim>-dimensional mesh\n"
           << "    -i <string>                   Input file name\n"
           << "    -o <string>                   Output file name\n"
           << "\n    -b                            Write the boundary conditions\n"
           << "    -D <factor>                   Randomly move interior nodes by D*hmin\n"
           << "    -h                            Print help menu\n"
           << "    -p <count>                    Partition into <count> subdomains\n"
#ifdef LIBMESH_ENABLE_AMR
           << "    -r <count>                    Globally refine <count> times\n"
#endif
           << "    -t (-d 2 only)                Convert to triangles first\n"
           << "                                  (allows to write .unv file of the\n"
           << "                                  boundary with the correct node ids)\n"
           << "    -v                            Verbose\n"
           << "    -q <metric>                   Evaluates the named element quality metric\n"
           << "    -1                            Converts a mesh of higher order elements\n"
           << "                                  to their first-order counterparts:\n"
           << "                                  Quad8 -> Quad4, Tet10 -> Tet4 etc\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 LIBMESH_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"
           << "\n"
           << " This program is used to convert and partition 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.e -o out.plt\n"
           << "\n"
           << " will read a 2D mesh in the ExodusII format (from Cubit, for example)\n"
           << " from the file in.e.  It will then write the mesh in the Tecplot\n"
           << " binary format to out.plt.\n"
           << "\n"
           << " and\n"
           << "\n"
           << "  ./meshtool -d 3 -i cylinder.msh -o out.e -2\n"
           << "\n"
           << " will read the Gmsh formatted mesh file cylinder.msh, convert\n"
           << " the elements to second-order, and write the output in Exodus format\n"
           << " to out.e.\n"
#ifdef LIBMESH_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"
           << "     .e   -- Sandia's ExodusII binary grid format\n"
           << "     .inp -- Abaqus mesh file\n"
           << "     .msh -- GMSH ASCII file\n"
           << "     .ucd -- AVS unstructured ASCII grid format\n"
           << "     .unv -- SDRC I-Deas Universal File ASCII format\n"
           << "     .xda -- libMesh human-readable ASCII format\n"
           << "     .xdr -- libMesh binary format\n"
           << "\n"
           << "OUTPUT:\n"
           << "     .dat   -- Tecplot ASCII format\n"
           << "     .e     -- Sandia's ExodusII format\n"
           << "     .exd   -- Sandia's ExodusII format\n"
           << "     .fro   -- ACDL's .fro format\n"
           << "     .gmv   -- LANL's General Mesh Viewer format\n"
           << "     .mesh  -- MEdit mesh format\n"
           << "     .msh   -- GMSH ASCII file\n"
           << "     .plt   -- Tecplot binary format\n"
           << "     .poly  -- TetGen ASCII file\n"
           << "     .pvtu  -- Paraview VTK format\n"
           << "     .ucd   -- AVS's ASCII UCD format\n"
           << "     .unv   -- I-deas Universal format\n"
           << "     .xda   -- libMesh ASCII format\n"
           << "     .xdr   -- libMesh binary format\n"
           << "     .gz    -- any above format gzipped\n"
           << "     .bz2   -- any above format bzip2'ed\n"
           << "\n";
 
  libMesh::out << helpList.str();
  exit(0);
}