meshtool.cc

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

int main (int argc, char** argv)
{
  LibMeshInit init(argc, argv);

  PerfMon perfmon(argv[0]);
  
  unsigned int n_subdomains = 1;
  unsigned int n_rsteps = 0;
  unsigned int dim = static_cast<unsigned int>(-1); // invalid dimension
  double dist_fact = 0.;
  bool verbose = false;
  BoundaryMeshWriteMode write_bndry = BM_DISABLED;
  unsigned int convert_second_order = 0;
  bool addinfelems = false;
  bool triangulate = false;
  
#ifdef ENABLE_INFINITE_ELEMENTS
  InfElemBuilder::InfElemOriginValue origin_x(false, 0.);
  InfElemBuilder::InfElemOriginValue origin_y(false, 0.);
  InfElemBuilder::InfElemOriginValue origin_z(false, 0.);
#endif

  bool x_sym=false;
  bool y_sym=false;
  bool z_sym=false;
  
    
  std::vector<std::string> names;
  std::vector<std::string> var_names;
  std::vector<Number>      soln;

  process_cmd_line(argc, argv, names,
                   n_subdomains, n_rsteps,
                   dim, dist_fact, verbose, write_bndry, 
                   convert_second_order,

                   triangulate,

                   addinfelems,

#ifdef ENABLE_INFINITE_ELEMENTS
                   origin_x, origin_y, origin_z, 
#endif

                   x_sym, y_sym, z_sym);

  if (dim == static_cast<unsigned int>(-1))
    {
      std::cout << "ERROR:  you must specify the dimension on "
                << "the command line!\n\n"
                << argv[0] << " -d 3 ... for example\n\n";
      libmesh_error();
    }

  Mesh mesh(dim);
  MeshData mesh_data(mesh);
  
  /**
   * Read the input mesh
   */      
  if (!names.empty())
    {
      /*
       * activate the MeshData of the dim mesh,
       * so that it can be copied to the boundary
       */

      if (write_bndry == BM_WITH_MESHDATA)
        {
          mesh_data.activate();

          if (verbose)
            mesh_data.print_info();
        }


      mesh.read(names[0]);

      if (verbose)
        mesh.print_info();

    }
  
  else
    {
      std::cout << "No input specified." << std::endl;
      return 1;
    }
             


#ifdef ENABLE_INFINITE_ELEMENTS

  if(addinfelems)
    {
      if (names.size() == 3)
        {
          std::cout << "ERROR: Invalid combination: Building infinite elements " << std::endl
                    << "not compatible with solution import." << std::endl;
          exit(1);
        }
      
      if (write_bndry != BM_DISABLED)
        {
          std::cout << "ERROR: Invalid combination: Building infinite elements " << std::endl
                    << "not compatible with writing boundary conditions." << std::endl;
          exit(1);
        }

      /*
       * 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
        {
          std::cout << "ERROR: When x-symmetry is requested using -X, then" << std::endl
                    << "the option -x <coord> also has to be given." << std::endl
                    << "This holds obviously for y and z, too." << std::endl;
          exit(1);
        }


      // 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();

    }

  // sanity check
  else if((origin_x.first ||  origin_y.first || origin_z.first) ||
          (x_sym          ||  y_sym          || z_sym))
    {
      std::cout << "ERROR:  -x/-y/-z/-X/-Y/-Z is only to be used when" << std::endl
                << "the option -a is also specified!" << std::endl;
      exit(1);
    }
  
#endif


  /**
   * Possibly read the solution
   */
  if (names.size() == 3)
    LegacyXdrIO(mesh,true).read_mgf_soln(names[2],
                                         soln,
                                         var_names);



  /**
   * Maybe Triangulate
   */
  if (dim == 2 && triangulate)
    {
      if (verbose) std::cout << "...Converting to all triangles...\n";
                     
      MeshTools::Modification::all_tri(mesh);
    }
  
  /**
   * Compute Shape quality metrics
   */
  const bool do_quality = false;

  if (do_quality)
    {
      StatisticsVector<Real> sv;
      sv.resize(mesh.n_elem());
      
      const ElemQuality q = DIAGONAL;

      std::cout << "Quality type is: " << Quality::name(q) << std::endl;
      
      // What are the quality bounds for this element?
      std::pair<Real, Real> bounds = mesh.elem(0)->qual_bounds(q);
      std::cout << "Quality bounds for this element type are: (" << bounds.first
                << ", " << bounds.second << ") "
                << std::endl;

      for (unsigned int e=0; e<mesh.n_elem(); e++)
        {
          sv[e] = mesh.elem(e)->quality(q);
        }

      const unsigned int n_bins = 10;
      std::cout << "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<unsigned int> bad_elts  = sv.cut_below(0.8);
      std::cout << "Found " << bad_elts.size()
                << " of " << mesh.n_elem()
                << " elements below the cutoff." << std::endl;
      
      /*
      for (unsigned int i=0; i<bad_elts.size(); i++)
        std::cout << bad_elts[i] << " ";
      std::cout << std::endl;
      */
      
      // Compute the histogram for this distribution
      std::vector<unsigned int> histogram;
      sv.histogram(histogram, n_bins);
      
      /*
      for (unsigned int i=0; i<n_bins; i++)
        histogram[i] = histogram[i] / mesh.n_elem();
      */
      
      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(q) << "');" << std::endl;
          
        }
    }


  /**
   * 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();

      if (verbose)
        std::cout << message << std::endl;
      
      mesh.all_second_order(second_order_mode);
      
      if (verbose)
        mesh.print_info();
    }

  
#ifdef ENABLE_AMR
  
  /**
   * Possibly refine the mesh
   */
  if (n_rsteps > 0)
    {
      if (verbose)
        std::cout << "Refining the mesh "
                  << n_rsteps << " times"
                  << std::endl;

      MeshRefinement mesh_refinement (mesh);
      mesh_refinement.uniformly_refine(n_rsteps);
      
      if (verbose)
        mesh.print_info();
    }

  
  /**
   * Possibly distort the mesh
   */
  if (dist_fact > 0.)
    {
      std::cout << "Distoring the mesh by a factor of "
                << dist_fact
                << std::endl;
      
      MeshTools::Modification::distort(mesh,dist_fact);
    };


  /*
  char filechar[81];
  sprintf(filechar,"%s-%04d.plt", "out", 0);
  std::string oname(filechar);
  
  mesh.write(oname);
  
  for (unsigned int step=0; step<100; step++)
    {
//        const Real x = .5 + .25*cos((((Real) step)/100.)*3.1415927); 
//        const Real y = .5 + .25*sin((((Real) step)/100.)*3.1415927);
      const Real x = 2.5*cos((((Real) step)/100.)*3.1415927); 
      const Real y = 2.5*sin((((Real) step)/100.)*3.1415927);

      const Point p(x,y);
      
      for (unsigned int e=0; e<mesh.n_elem(); e++)
        if (mesh.elem(e)->active())
          mesh.elem(e)->set_refinement_flag() = -1;
      
      
      
      for (unsigned int e=0; e<mesh.n_elem(); e++)
        if (mesh.elem(e)->active())
          {
            const Point diff = mesh.elem(e)->centroid(mesh) - p;

            if (diff.size() < .5)
              {
                if (mesh.elem(e)->level() < 4)
                  mesh.elem(e)->set_refinement_flag() = 1;
                else if (mesh.elem(e)->level() == 4)
                  mesh.elem(e)->set_refinement_flag() = 0;
              }
          }
      

      mesh.mesh_refinement.refine_and_coarsen_elements();
      
      char filechar[81];
      sprintf(filechar,"%s-%04d.plt", "out", step+1);
      std::string oname(filechar);

      mesh.write(oname);
    }
  */

#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)
                std::cout << " Mesh got refined, will write only _active_ elements." << std::endl;

            Mesh new_mesh (dim);

            construct_mesh_of_active_elements(new_mesh, mesh);

            // 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();
          }
        else
          {
            if (names.size() == 2)
                mesh.write(names[1]);
            else if (names.size() == 3)
                mesh.write(names[1], soln, var_names);
            else
                libmesh_error();
          }



        /**
         * Possibly write the BCs
         */
        if (write_bndry != BM_DISABLED)
          {
            BoundaryMesh boundary_mesh (mesh.mesh_dimension()-1);
            MeshData boundary_mesh_data (boundary_mesh);
            
            std::string boundary_name = "bndry_";
            boundary_name += names[1];
            
            if (write_bndry == BM_MESH_ONLY)
              mesh.boundary_info->sync(boundary_mesh);
            
            else if  (write_bndry == BM_WITH_MESHDATA)
              mesh.boundary_info->sync(boundary_mesh, &boundary_mesh_data, &mesh_data);

            else
              libmesh_error();
            
            if (names.size() == 2)
              GMVIO(boundary_mesh).write(boundary_name);
            else if (names.size() == 3)
              GMVIO(boundary_mesh).write_nodal_data(boundary_name,
                                                    soln, var_names);
          }
      }
  };
        
  /*
  std::cout << "Infinite loop, look at memory footprint" << std::endl;
  for (;;)
    ;
  */

  return libMesh::close();
}