C++ meshtools::BoundingBox类代码示例

void
PointSamplerBase::execute()
{
  MeshTools::BoundingBox bbox = _mesh.getInflatedProcessorBoundingBox();

  for (unsigned int i=0; i<_points.size(); i++)
  {
    Point & p = _points[i];

    // Do a bounding box check so we're not doing unnecessary PointLocator lookups
    if (bbox.contains_point(p))
    {
      std::vector<Real> & values = _values[i];

      if (values.empty())
        values.resize(_coupled_moose_vars.size());

      // First find the element the hit lands in
      const Elem * elem = getLocalElemContainingPoint(p, i);

      if (elem)
      {
        // We have to pass a vector of points into reinitElemPhys
        _point_vec[0] = p;

        _subproblem.reinitElemPhys(elem, _point_vec, 0); // Zero is for tid

        for (unsigned int j=0; j<_coupled_moose_vars.size(); j++)
          values[j] = _coupled_moose_vars[j]->sln()[0]; // The zero is for the "qp"

        _found_points[i] = true;
      }
    }
  }
}

void
FeatureFloodCount::FeatureData::updateBBoxExtremes(MeshTools::BoundingBox & bbox, const Point & node)
{
  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
  {
    bbox.min()(i) = std::min(bbox.min()(i), node(i));
    bbox.max()(i) = std::max(bbox.max()(i), node(i));
  }
}

MeshTools::BoundingBox
MultiApp::getBoundingBox(unsigned int app)
{
  if (!_has_an_app)
    mooseError("No app for " << name() << " on processor " << _orig_rank);

  FEProblem * problem = appProblem(app);

  MPI_Comm swapped = Moose::swapLibMeshComm(_my_comm);

  MooseMesh & mesh = problem->mesh();
  MeshTools::BoundingBox bbox = MeshTools::bounding_box(mesh);

  Moose::swapLibMeshComm(swapped);

  Point min = bbox.min();
  Point max = bbox.max();

  Point inflation_amount = (max-min)*_inflation;

  Point inflated_min = min - inflation_amount;
  Point inflated_max = max + inflation_amount;

  // This is where the app is located.  We need to shift by this amount.
  Point p = position(app);

  Point shifted_min = inflated_min;
  Point shifted_max = inflated_max;

  // If the problem is RZ then we're going to invent a box that would cover the whole "3D" app
  // FIXME: Assuming all subdomains are the same coordinate system type!
  if (problem->getCoordSystem(*(problem->mesh().meshSubdomains().begin())) == Moose::COORD_RZ)
  {
    shifted_min(0) = -inflated_max(0);
    shifted_min(1) = inflated_min(1);
    shifted_min(2) = -inflated_max(0);

    shifted_max(0) = inflated_max(0);
    shifted_max(1) = inflated_max(1);
    shifted_max(2) = inflated_max(0);
  }

  // Shift them to the position they're supposed to be
  shifted_min += p;
  shifted_max += p;

  return MeshTools::BoundingBox(shifted_min, shifted_max);
}

Teuchos::RCP<DataTransferKit::GeometryManager<DataTransferKit::Box,long unsigned int> >
MultiAppDTKUserObjectEvaluator::createSourceGeometry( const Teuchos::RCP<const Teuchos::Comm<int> >& comm )
{
  _boxes.resize(_multi_app.numLocalApps());
  _box_ids.resize(_multi_app.numLocalApps());

  // How much to inflate each bounding box by (helps with non-matching surfaces)
  unsigned int inflation = 0.01;

  comm->barrier();

  for(unsigned int app=0; app<_multi_app.numLocalApps(); app++)
  {
    unsigned int global_app = _multi_app.firstLocalApp() + app;

    MeshTools::BoundingBox bbox = _multi_app.getBoundingBox(global_app);

    _boxes[app] = DataTransferKit::Box(bbox.min()(0), bbox.min()(1), bbox.min()(2),
                                       bbox.max()(0), bbox.max()(1), bbox.max()(2));

    _box_ids[app] = global_app;
  }

  return Teuchos::rcp( new DataTransferKit::GeometryManager<DataTransferKit::Box,GlobalOrdinal>(_boxes, _box_ids, comm, 3));

}

void PostscriptIO::write (const std::string& fname)
{
  // We may need to gather a ParallelMesh to output it, making that
  // const qualifier in our constructor a dirty lie
  MeshSerializer serialize(const_cast<MeshBase&>(this->mesh()), !_is_parallel_format);

  if (this->mesh().processor_id() == 0)
    {
      // Get a constant reference to the mesh.
      const MeshBase& the_mesh = MeshOutput<MeshBase>::mesh();

      // Only works in 2D
      libmesh_assert_equal_to (the_mesh.mesh_dimension(), 2);

      // Create output file stream.
      // _out is now a private member of the class.
      _out.open(fname.c_str());

      // Make sure it opened correctly
      if (!_out.good())
        libmesh_file_error(fname.c_str());

      // The mesh bounding box gives us info about what the
      // Postscript bounding box should be.
      MeshTools::BoundingBox bbox = MeshTools::bounding_box(the_mesh);

      // Add a little extra padding to the "true" bounding box so
      // that we can still see the boundary
      const Real percent_padding = 0.01;
      const Real dx=bbox.second(0)-bbox.first(0); libmesh_assert_greater (dx, 0.0);
      const Real dy=bbox.second(1)-bbox.first(1); libmesh_assert_greater (dy, 0.0);

      const Real x_min = bbox.first(0)  - percent_padding*dx;
      const Real y_min = bbox.first(1)  - percent_padding*dy;
      const Real x_max = bbox.second(0) + percent_padding*dx;
      const Real y_max = bbox.second(1) + percent_padding*dy;

      // Width of the output as given in postscript units.
      // This usually is given by the strange unit 1/72 inch.
      // A width of 300 represents a size of roughly 10 cm.
      const Real width = 300;
      _scale = width / (x_max-x_min);
      _offset(0) = x_min;
      _offset(1) = y_min;

      // Header writing stuff stolen from Deal.II
      std::time_t  time1= std::time (0);
      std::tm     *time = std::localtime(&time1);
      _out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
	//<< "%!PS-Adobe-1.0" << '\n' // Lars' PS version
	  << "%%Filename: " << fname << '\n'
	  << "%%Title: LibMesh Output" << '\n'
	  << "%%Creator: LibMesh: A C++ finite element library" << '\n'
	  << "%%Creation Date: "
	  << time->tm_year+1900 << "/"
	  << time->tm_mon+1 << "/"
	  << time->tm_mday << " - "
	  << time->tm_hour << ":"
	  << std::setw(2) << time->tm_min << ":"
	  << std::setw(2) << time->tm_sec << '\n'
	  << "%%BoundingBox: "
	// lower left corner
	  << "0 0 "
	// upper right corner
	  << static_cast<unsigned int>( rint((x_max-x_min) * _scale ))
	  << ' '
	  << static_cast<unsigned int>( rint((y_max-y_min) * _scale ))
	  << '\n';

      // define some abbreviations to keep
      // the output small:
      // m=move turtle to
      // l=define a line
      // s=set rgb color
      // sg=set gray value
      // lx=close the line and plot the line
      // lf=close the line and fill the interior
      _out << "/m {moveto} bind def"      << '\n'
	  << "/l {lineto} bind def"      << '\n'
	  << "/s {setrgbcolor} bind def" << '\n'
	  << "/sg {setgray} bind def"    << '\n'
	  << "/cs {curveto stroke} bind def" << '\n'
	  << "/lx {lineto closepath stroke} bind def" << '\n'
	  << "/lf {lineto closepath fill} bind def"   << '\n';

      _out << "%%EndProlog" << '\n';
      //	  << '\n';

      // Set line width in the postscript file.
      _out << line_width << " setlinewidth" << '\n';

      // Set line cap and join options
      _out << "1 setlinecap" << '\n';
      _out << "1 setlinejoin" << '\n';

      // allow only five digits for output (instead of the default
      // six); this should suffice even for fine grids, but reduces
      // the file size significantly
      _out << std::setprecision (5);

//.........这里部分代码省略.........

void
ImageSampler::setupImageSampler(MooseMesh & mesh)
{
#ifdef LIBMESH_HAVE_VTK
  // Get access to the Mesh object
  MeshTools::BoundingBox bbox = MeshTools::bounding_box(mesh.getMesh());

  // Set the dimensions from the Mesh if not set by the User
  if (_is_pars.isParamValid("dimensions"))
    _physical_dims = _is_pars.get<Point>("dimensions");

  else
  {
    _physical_dims(0) = bbox.max()(0) - bbox.min()(0);
#if LIBMESH_DIM > 1
    _physical_dims(1) = bbox.max()(1) - bbox.min()(1);
#endif
#if LIBMESH_DIM > 2
    _physical_dims(2) = bbox.max()(2) - bbox.min()(2);
#endif
  }

  // Set the origin from the Mesh if not set in the input file
  if (_is_pars.isParamValid("origin"))
    _origin = _is_pars.get<Point>("origin");
  else
  {
    _origin(0) = bbox.min()(0);
#if LIBMESH_DIM > 1
    _origin(1) = bbox.min()(1);
#endif
#if LIBMESH_DIM > 2
    _origin(2) = bbox.min()(2);
#endif
  }

  // An array of filenames, to be filled in
  std::vector<std::string> filenames;

  // The file suffix, to be determined
  std::string file_suffix;

  // Try to parse our own file range parameters.  If that fails, then
  // see if the associated Mesh is an ImageMesh and use its.  If that
  // also fails, then we have to throw an error...
  //
  // The parseFileRange method sets parameters, thus a writable reference to the InputParameters
  // object must be obtained from the warehouse. Generally, this should be avoided, but
  // this is a special case.
  if (_status != 0)
  {
    // We don't have parameters, so see if we can get them from ImageMesh
    ImageMesh * image_mesh = dynamic_cast<ImageMesh*>(&mesh);
    if (!image_mesh)
      mooseError("No file range parameters were provided and the Mesh is not an ImageMesh.");

    // Get the ImageMesh's parameters.  This should work, otherwise
    // errors would already have been thrown...
    filenames = image_mesh->filenames();
    file_suffix = image_mesh->fileSuffix();
  }
  else
  {
    // Use our own parameters (using 'this' b/c of conflicts with filenames the local variable)
    filenames = this->filenames();
    file_suffix = fileSuffix();
  }

  // Storage for the file names
  _files = vtkSmartPointer<vtkStringArray>::New();

  for (unsigned i=0; i<filenames.size(); ++i)
    _files->InsertNextValue(filenames[i]);

  // Error if no files where located
  if (_files->GetNumberOfValues() == 0)
    mooseError("No image file(s) located");


  // Read the image stack.  Hurray for VTK not using polymorphism in a
  // smart way... we actually have to explicitly create the type of
  // reader based on the file extension, using an if-statement...
  if (file_suffix == "png")
    _image = vtkSmartPointer<vtkPNGReader>::New();
  else if (file_suffix == "tiff" || file_suffix == "tif")
    _image = vtkSmartPointer<vtkTIFFReader>::New();
  else
    mooseError("Un-supported file type '" << file_suffix << "'");

  // Now that _image is set up, actually read the images
  // Indicate that data read has started
  _is_console << "Reading image(s)..." << std::endl;

  // Extract the data
  _image->SetFileNames(_files);
  _image->Update();
  _data = _image->GetOutput();
  _algorithm = _image->GetOutputPort();

  // Set the image dimensions and voxel size member variable
//.........这里部分代码省略.........

void
ClosePackIC::computeCircleCenters()
{
  // Determine the extents of the mesh
  MeshTools::BoundingBox bbox = MeshTools::bounding_box(_fe_problem.mesh().getMesh());
  const Point & min = bbox.min();
  const Point & max = bbox.max();


  // Create the x,y,z limits for the while loops
  Real x_min = min(0);
  Real x_max = max(0) + 2*_radius;

  Real y_min = min(1) - 2*std::sqrt(3.0)*_radius + _radius;
  Real y_max = max(1) + 2*_radius;

  // Real z_min = min(2) - 2*std::sqrt(3.0)*_radius + _radius;
  Real z_max = 0;

  // Initialize the coordinates that will be used in looping
  Real x = x_min;
  Real y = y_min;
  Real z = 0;

  // Adjust the 3D z-dimension maximum
  if (_fe_problem.mesh().dimension() == 3)
    z_max = max(2) + 2*_radius;

  // Counters for offsetting every other row column in x,y dimensions
  unsigned int i = 0;
  unsigned int j = 0;

  while (z <= z_max)
  {

    // Offset the y-coordinate by sqrt(3)*r every other loop
    if (j % 2 != 0)
      y += std::sqrt(3)*_radius/2;

    while (y <= y_max)
    {

      // Offset the x-coordinate by r every other loop
      if (i % 2 == 0)
        x += _radius;

      while (x <= x_max)
      {
        Point p(x, y, z);
        _centers.push_back(p);
        _radii.push_back(_radius);
        x += 2*_radius;
      }

      // Reset x-coord and increment y-coord
      x = x_min;
      y += std::sqrt(3.)*_radius;
      i++;
    }

    // Reset y-coord and increment z-coord
    y = y_min;
    z += std::sqrt(3.) * _radius;
    j++;
  }
}

void PointLocatorTree::init (Trees::BuildType build_type)
{
  libmesh_assert (!this->_tree);

  if (this->_initialized)
    {
      // Warn that we are already initialized
      libMesh::err << "Warning: PointLocatorTree already initialized!  Will ignore this call..." << std::endl;

      // Further warn if we try to init() again with a different build_type
      if (_build_type != build_type)
        {
          libMesh::err << "Warning: PointLocatorTree is using build_type = " << _build_type << ".\n"
                       << "Your requested build_type, " << build_type << " will not be used!" << std::endl;
        }
    }

  else
    {
      // Let the requested build_type override the _build_type we were
      // constructed with.  This is no big deal since we have not been
      // initialized before.
      _build_type = build_type;

      if (this->_master == libmesh_nullptr)
        {
          LOG_SCOPE("init(no master)", "PointLocatorTree");

          if (this->_mesh.mesh_dimension() == 3)
            _tree = new Trees::OctTree (this->_mesh, get_target_bin_size(), _build_type);
          else
            {
              // A 1D/2D mesh in 3D space needs special consideration.
              // If the mesh is planar XY, we want to build a QuadTree
              // to search efficiently.  If the mesh is truly a manifold,
              // then we need an octree
#if LIBMESH_DIM > 2
              bool is_planar_xy = false;

              // Build the bounding box for the mesh.  If the delta-z bound is
              // negligibly small then we can use a quadtree.
              {
                MeshTools::BoundingBox bbox = MeshTools::bounding_box(this->_mesh);

                const Real
                  Dx = bbox.second(0) - bbox.first(0),
                  Dz = bbox.second(2) - bbox.first(2);

                if (std::abs(Dz/(Dx + 1.e-20)) < 1e-10)
                  is_planar_xy = true;
              }

              if (!is_planar_xy)
                _tree = new Trees::OctTree (this->_mesh, get_target_bin_size(), _build_type);
              else
#endif
#if LIBMESH_DIM > 1
                _tree = new Trees::QuadTree (this->_mesh, get_target_bin_size(), _build_type);
#else
              _tree = new Trees::BinaryTree (this->_mesh, get_target_bin_size(), _build_type);
#endif
            }
        }

      else
        {
          // We are _not_ the master.  Let our Tree point to
          // the master's tree.  But for this we first transform
          // the master in a state for which we are friends.
          // And make sure the master @e has a tree!
          const PointLocatorTree * my_master =
            cast_ptr<const PointLocatorTree *>(this->_master);

          if (my_master->initialized())
            this->_tree = my_master->_tree;
          else
            libmesh_error_msg("ERROR: Initialize master first, then servants!");
        }

      // Not all PointLocators may own a tree, but all of them
      // use their own element pointer.  Let the element pointer
      // be unique for every interpolator.
      // Suppose the interpolators are used concurrently
      // at different locations in the mesh, then it makes quite
      // sense to have unique start elements.
      this->_element = libmesh_nullptr;
    }

  // ready for take-off
  this->_initialized = true;
}

void
FeatureFloodCount::calculateBubbleVolumes()
{
  Moose::perf_log.push("calculateBubbleVolume()", "FeatureFloodCount");

  // Figure out which bubbles intersect the boundary if the user has enabled that capability.
  if (_compute_boundary_intersecting_volume)
  {
    // Create a std::set of node IDs which are on the boundary called all_boundary_node_ids.
    std::set<dof_id_type> all_boundary_node_ids;

    // Iterate over the boundary nodes, putting them into the std::set data structure
    MooseMesh::bnd_node_iterator
      boundary_nodes_it  = _mesh.bndNodesBegin(),
      boundary_nodes_end = _mesh.bndNodesEnd();
    for (; boundary_nodes_it != boundary_nodes_end; ++boundary_nodes_it)
    {
      BndNode * boundary_node = *boundary_nodes_it;
      all_boundary_node_ids.insert(boundary_node->_node->id());
    }

    // For each of the _maps_size BubbleData lists, determine if the set
    // of nodes includes any boundary nodes.
    for (unsigned int map_num = 0; map_num < _maps_size; ++map_num)
    {
      std::list<BubbleData>::iterator
        bubble_it = _bubble_sets[map_num].begin(),
        bubble_end = _bubble_sets[map_num].end();

      // Determine boundary intersection for each BubbleData object
      for (; bubble_it != bubble_end; ++bubble_it)
        bubble_it->_intersects_boundary = setsIntersect(all_boundary_node_ids.begin(), all_boundary_node_ids.end(),
                                                        bubble_it->_entity_ids.begin(), bubble_it->_entity_ids.end());
    }
  }

  // Size our temporary data structure
  std::vector<std::vector<Real> > bubble_volumes(_maps_size);
  for (unsigned int map_num = 0; map_num < _maps_size; ++map_num)
    bubble_volumes[map_num].resize(_bubble_sets[map_num].size());

  // Clear pre-existing values and allocate space to store the volume
  // of the boundary-intersecting grains for each variable.
  _total_volume_intersecting_boundary.clear();
  _total_volume_intersecting_boundary.resize(_maps_size);

  // Loop over the active local elements.  For each variable, and for
  // each BubbleData object, check whether a majority of the element's
  // nodes belong to that Bubble, and if so assign the element's full
  // volume to that bubble.
  const MeshBase::const_element_iterator el_end = _mesh.getMesh().active_local_elements_end();
  for (MeshBase::const_element_iterator el = _mesh.getMesh().active_local_elements_begin(); el != el_end; ++el)
  {
    Elem * elem = *el;
    unsigned int elem_n_nodes = elem->n_nodes();
    Real curr_volume = elem->volume();

    for (unsigned int map_num = 0; map_num < _maps_size; ++map_num)
    {
      std::list<BubbleData>::const_iterator
        bubble_it = _bubble_sets[map_num].begin(),
        bubble_end = _bubble_sets[map_num].end();

      for (unsigned int bubble_counter = 0; bubble_it != bubble_end; ++bubble_it, ++bubble_counter)
      {
        // Count the number of nodes on this element which are flooded.
        unsigned int flooded_nodes = 0;
        for (unsigned int node = 0; node < elem_n_nodes; ++node)
        {
          dof_id_type node_id = elem->node(node);
          if (bubble_it->_entity_ids.find(node_id) != bubble_it->_entity_ids.end())
            ++flooded_nodes;
        }

        // If a majority of the nodes for this element are flooded,
        // assign its volume to the current bubble_counter entry.
        if (flooded_nodes >= elem_n_nodes / 2)
        {
          bubble_volumes[map_num][bubble_counter] += curr_volume;

          // If the current bubble also intersects the boundary, also
          // accumlate the volume into the total volume of bubbles
          // which intersect the boundary.
          if (bubble_it->_intersects_boundary)
            _total_volume_intersecting_boundary[map_num] += curr_volume;
        }
      }
    }
  }

  // If we're calculating boundary-intersecting volumes, we have to normalize it by the
  // volume of the entire domain.
  if (_compute_boundary_intersecting_volume)
  {
    // Compute the total area using a bounding box.  FIXME: this
    // assumes the domain is rectangular and 2D, and is probably a
    // little expensive so we should only do it once if possible.
    MeshTools::BoundingBox bbox = MeshTools::bounding_box(_mesh);
    Real total_volume = (bbox.max()(0)-bbox.min()(0))*(bbox.max()(1)-bbox.min()(1));

//.........这里部分代码省略.........

void dataLoad(std::istream & stream, MeshTools::BoundingBox & bbox, void * context)
{
  loadHelper(stream, bbox.min(), context);
  loadHelper(stream, bbox.max(), context);
}

void dataStore(std::ostream & stream, MeshTools::BoundingBox & bbox, void * context)
{
  storeHelper(stream, bbox.min(), context);
  storeHelper(stream, bbox.max(), context);
}

void PointLocatorTree::init (const Trees::BuildType build_type)
{
  libmesh_assert (!this->_tree);

  if (this->_initialized)
    {
      libMesh::err << "ERROR: Already initialized!  Will ignore this call..."
                   << std::endl;
    }

  else

    {

      if (this->_master == NULL)
        {
          START_LOG("init(no master)", "PointLocatorTree");

          if (this->_mesh.mesh_dimension() == 3)
            _tree = new Trees::OctTree (this->_mesh, 200, build_type);
          else
            {
              // A 1D/2D mesh in 3D space needs special consideration.
              // If the mesh is planar XY, we want to build a QuadTree
              // to search efficiently.  If the mesh is truly a manifold,
              // then we need an octree
#if LIBMESH_DIM > 2
              bool is_planar_xy = false;

              // Build the bounding box for the mesh.  If the delta-z bound is
              // negligibly small then we can use a quadtree.
              {
                MeshTools::BoundingBox bbox = MeshTools::bounding_box(this->_mesh);

                const Real
                  Dx = bbox.second(0) - bbox.first(0),
                  Dz = bbox.second(2) - bbox.first(2);

                if (std::abs(Dz/(Dx + 1.e-20)) < 1e-10)
                  is_planar_xy = true;
              }

              if (!is_planar_xy)
                _tree = new Trees::OctTree (this->_mesh, 200, build_type);
              else
#endif
#if LIBMESH_DIM > 1
                _tree = new Trees::QuadTree (this->_mesh, 200, build_type);
#else
              _tree = new Trees::BinaryTree (this->_mesh, 200, build_type);
#endif
            }

          STOP_LOG("init(no master)", "PointLocatorTree");
        }

      else

        {
          // We are _not_ the master.  Let our Tree point to
          // the master's tree.  But for this we first transform
          // the master in a state for which we are friends.
          // And make sure the master @e has a tree!
          const PointLocatorTree* my_master =
            libmesh_cast_ptr<const PointLocatorTree*>(this->_master);

          if (my_master->initialized())
            this->_tree = my_master->_tree;
          else
            {
              libMesh::err << "ERROR: Initialize master first, then servants!"
                           << std::endl;
              libmesh_error();
            }
        }


      // Not all PointLocators may own a tree, but all of them
      // use their own element pointer.  Let the element pointer
      // be unique for every interpolator.
      // Suppose the interpolators are used concurrently
      // at different locations in the mesh, then it makes quite
      // sense to have unique start elements.
      this->_element = NULL;
    }


  // ready for take-off
  this->_initialized = true;
}

void
FeatureFloodCount::FeatureData::updateBBoxMax(MeshTools::BoundingBox & bbox, const Point & max)
{
  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
    bbox.max()(i) = std::max(bbox.max()(i), max(i));
}