meshbcid.C File Reference

Go to the source code of this file.

Functions
void	usage_error (const char *progname)
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 48 of file meshbcid.C.

References libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::FEGenericBase< OutputType >::build(), libMesh::CONSTANT, libMesh::BoundingBox::contains_point(), dim, libMesh::err, libMesh::FIRST, libMesh::MeshBase::get_boundary_info(), libMesh::TriangleWrapper::init(), libMesh::BoundaryInfo::invalid_id, libMesh::LAGRANGE, libMesh::libmesh_assert(), libMesh::BoundingBox::max(), mesh, libMesh::BoundingBox::min(), libMesh::out, libMesh::MeshBase::read(), libMesh::BoundaryInfo::regenerate_id_sets(), libMesh::BoundaryInfo::remove_side(), usage_error(), and libMesh::MeshBase::write().

{
  LibMeshInit init(argc, argv);

  GetPot cl(argc, argv);

  unsigned char dim = -1;
  if (!cl.search("--dim"))
    {
      libMesh::err << "No --dim argument found!" << std::endl;
      usage_error(argv[0]);
    }
  dim = cl.next(dim);

  Mesh mesh(init.comm(), dim);

  if (!cl.search("--input"))
    {
      libMesh::err << "No --input argument found!" << std::endl;
      usage_error(argv[0]);
    }
  const char * meshname = cl.next("mesh.xda");

  mesh.read(meshname);
  libMesh::out << "Loaded mesh " << meshname << std::endl;

  if (!cl.search("--newbcid"))
    {
      libMesh::err << "No --bcid argument found!" << std::endl;
      usage_error(argv[0]);
    }
  boundary_id_type bcid = 0;
  bcid = cl.next(bcid);

  Point minpt(-std::numeric_limits<Real>::max());
#if LIBMESH_DIM > 1
  minpt(1) = -std::numeric_limits<Real>::max();
#endif
#if LIBMESH_DIM > 2
  minpt(2) = -std::numeric_limits<Real>::max();
#endif
  Point maxpt = -minpt;

  BoundingBox normals(minpt, maxpt),
              points(minpt, maxpt);

  if (cl.search("--minnormalx"))
    normals.min()(0) = cl.next(normals.min()(0));
  if (cl.search("--maxnormalx"))
    normals.max()(0) = cl.next(normals.max()(0));

  if (cl.search("--minpointx"))
    points.min()(0) = cl.next(points.min()(0));
  if (cl.search("--maxpointx"))
    points.max()(0) = cl.next(points.max()(0));

#if LIBMESH_DIM > 1
  if (cl.search("--minnormaly"))
    normals.min()(1) = cl.next(normals.min()(1));
  if (cl.search("--maxnormaly"))
    normals.max()(1) = cl.next(normals.max()(1));

  if (cl.search("--minpointy"))
    points.min()(1) = cl.next(points.min()(1));
  if (cl.search("--maxpointy"))
    points.max()(1) = cl.next(points.max()(1));
#endif

#if LIBMESH_DIM > 2
  if (cl.search("--minnormalz"))
    normals.min()(2) = cl.next(normals.min()(2));
  if (cl.search("--maxnormalz"))
    normals.max()(2) = cl.next(normals.max()(2));

  if (cl.search("--minpointz"))
    points.min()(2) = cl.next(points.min()(2));
  if (cl.search("--maxpointz"))
    points.max()(2) = cl.next(points.max()(2));
#endif

  libMesh::out << "min point = " << points.min() << std::endl;
  libMesh::out << "max point = " << points.max() << std::endl;
  libMesh::out << "min normal = " << normals.min() << std::endl;
  libMesh::out << "max normal = " << normals.max() << std::endl;

  bool matcholdbcid = false;
  boundary_id_type oldbcid = 0;
  if (cl.search("--oldbcid"))
    {
      matcholdbcid = true;
      oldbcid = cl.next(oldbcid);
      if (oldbcid < 0)
        oldbcid = BoundaryInfo::invalid_id;
    }

  std::unique_ptr<FEBase> fe = FEBase::build(dim, FEType(FIRST,LAGRANGE));
  QGauss qface(dim-1, CONSTANT);
  fe->attach_quadrature_rule(&qface);
  const std::vector<Point> & face_points = fe->get_xyz();
  const std::vector<Point> & face_normals = fe->get_normals();

  for (auto & elem : mesh.element_ptr_range())
    {
      unsigned int n_sides = elem->n_sides();

      // Container to catch ids handed back from BoundaryInfo
      std::vector<boundary_id_type> ids;

      for (unsigned short s=0; s != n_sides; ++s)
        {
          if (elem->neighbor_ptr(s))
            continue;

          fe->reinit(elem,s);
          const Point & p = face_points[0];
          const Point & n = face_normals[0];

          //libMesh::out << "elem = " << elem->id() << std::endl;
          //libMesh::out << "vertex average = " << elem->vertex_average() << std::endl;
          //libMesh::out << "p = " << p << std::endl;
          //libMesh::out << "n = " << n << std::endl;

          if (points.contains_point(p) &&
              normals.contains_point(n))
            {
              // Get the list of boundary ids for this side
              mesh.get_boundary_info().boundary_ids(elem, s, ids);

              // There should be at most one value present, otherwise the
              // logic here won't work.
              libmesh_assert(ids.size() <= 1);

              // A convenient name for the side's ID.
              boundary_id_type b_id = ids.empty() ? BoundaryInfo::invalid_id : ids[0];

              if (matcholdbcid && b_id != oldbcid)
                continue;

              mesh.get_boundary_info().remove_side(elem, s);
              mesh.get_boundary_info().add_side(elem, s, bcid);
              //libMesh::out << "Set element " << elem->id() << " side " << s <<
              //                " to boundary " << bcid << std::endl;
            }
        }
    }

  // We might have removed *every* instance of a given id, and if that
  // happened then we should make sure that file formats which write
  // out id sets do not write out the removed id.
  mesh.get_boundary_info().regenerate_id_sets();

  std::string outputname;
  if (cl.search("--output"))
    {
      outputname = cl.next("mesh.xda");
    }
  else
    {
      outputname = "new.";
      outputname += meshname;
    }


  mesh.write(outputname.c_str());
  libMesh::out << "Wrote mesh " << outputname << std::endl;

  return 0;
}

usage_error()

void usage_error

(

const char *

progname

)

Definition at line 39 of file meshbcid.C.

References libMesh::out.

Referenced by main().

{
  libMesh::out << "Usage: " << progname
               << " --dim d --input inputmesh --output outputmesh --newbcid idnum --tests --moretests"
               << std::endl;

  exit(1);
}