GRINS::RayfireMesh Class Reference

RayfireMesh. More...

#include <rayfire_mesh.h>

Public Member Functions
	RayfireMesh (libMesh::Point &origin, libMesh::Real theta)
	2D Constructor More...
	RayfireMesh (libMesh::Point &origin, libMesh::Real theta, libMesh::Real phi)
	3D Constructor More...
	RayfireMesh (const GetPot &input, const std::string &qoi_string)
	Input File Constructor. More...
	RayfireMesh (const RayfireMesh &original)
	Copy Constructor. More...
void	init (const libMesh::MeshBase &mesh_base)
	Initialization. More...
const libMesh::Elem *	map_to_rayfire_elem (const libMesh::dof_id_type elem_id)
	This function takes in an elem_id on the main mesh and returns an elem from the 1D rayfire mesh. More...
void	elem_ids_in_rayfire (std::vector< libMesh::dof_id_type > &id_vector) const
	Returns a vector of elem IDs that are currently along the rayfire. More...
void	reinit (const libMesh::MeshBase &mesh_base)
	Checks for refined and coarsened main mesh elements along the rayfire path. More...

Private Member Functions
	RayfireMesh ()
	User should never call the default constructor. More...
libMesh::Elem *	get_rayfire_elem (const libMesh::dof_id_type elem_id)
	Private function to get a rayfire elem from main_mesh elem ID. More...
const libMesh::Elem *	get_next_elem (const libMesh::Elem *cur_elem, libMesh::Point &start_point, libMesh::Point &end_point, bool same_parent=false)
	Calculate the intersection point. More...
bool	check_valid_point (libMesh::Point &intersection_point, libMesh::Point &start_point, const libMesh::Elem &edge_elem, libMesh::Point &next_point)
	Ensure the calculated intersection point is on the edge_elem and is not the start_point. More...
const libMesh::Elem *	get_correct_neighbor (libMesh::Point &start_point, libMesh::Point &end_point, const libMesh::Elem *cur_elem, unsigned int side, bool same_parent)
	Knowing the end_point, get the appropraite next elem along the path. More...
void	check_origin_on_boundary (const libMesh::Elem *start_elem)
	Ensure the supplied origin is on a boundary of the mesh. More...
const libMesh::Elem *	get_start_elem (const libMesh::MeshBase &mesh_base)
	Get the correct starting elem corresponding to _origin. More...
bool	validate_edge (const libMesh::Point &start_point, const libMesh::Point &end_point, const libMesh::Elem *side_elem, const libMesh::Elem *neighbor)
	Ensures the rayfire doesn't wander into the middle of an elem. More...
bool	rayfire_in_elem (const libMesh::Point &end_point, const libMesh::Elem *elem)
	Walks a short distance along the rayfire and checks if elem contains that point. More...
bool	newton_solve_intersection (libMesh::Point &initial_point, const libMesh::Elem *edge_elem, libMesh::Point &intersection_point)
	Iterative solver for calculating the intersection point of the rayfire on the side of an elem. More...
void	refine (const libMesh::Elem *main_elem, libMesh::Elem *rayfire_elem)
	Refinement of a rayfire element whose main mesh counterpart was refined. More...
void	coarsen (const libMesh::Elem *rayfire_elem)
	Coarsening of a rayfire element whose main mesh counterpart was coarsened. More...

Private Attributes
const unsigned int	_dim
	Dimension of the main mesh. More...
libMesh::Point	_origin
	Origin point. More...
libMesh::Real	_theta
	Rayfire Spherical polar angle (in radians) More...
libMesh::Real	_phi
	Rayfire Spherical azimuthal angle (in radians) More...
libMesh::UniquePtr< libMesh::Mesh >	_mesh
	Internal 1D mesh of EDGE2 elements. More...
std::map< libMesh::dof_id_type, libMesh::dof_id_type >	_elem_id_map
	Map of main mesh elem_id to rayfire mesh elem_id. More...

Detailed Description

RayfireMesh.

This class performs a rayfire across a given mesh. The user supplies an origin point on the mesh boundary, and spherical angle(s) theta (and phi for 3D) to describe the direction of the rayfire. As is standard in spherical coordinates, theta is the angle from the x-axis in the xy-plane, with counterclockwise being positive. -2π ≤ theta ≤ +2π. Phi is the angle from the xy-plane, with phi>0 being in the positive half of the z-plane. -π ≤ phi ≤ +π

Starting at the given origin point, this class will walk in a straight line across the mesh in the prescribed direction. On each element along the line, the class will calculate the point where it enters and leaves that element. Knowing the entering and exiting points, an EDGE2 element is created with those points as nodes, and is added to an internal 1D mesh. This repeats until a main mesh boundary is reached.

The intent of this class is to use the 1D mesh to integrate quantities along the rayfire path. Using the map_to_rayfire_elem() function, a separate integration class can access the 1D elements, prescribe a desired quadrature rule, and perform integration directly along the line, rather than using the entire 2D/3D elements of the main mesh.

Refinement and coarsening of the rayfire mesh are supported through the reinit() function.

Definition at line 65 of file rayfire_mesh.h.

Constructor & Destructor Documentation

GRINS::RayfireMesh::RayfireMesh	(	libMesh::Point &	origin,
		libMesh::Real	theta
	)

2D Constructor

This is restricted to 2D meshes. The init() function must be called to perform the actual rayfire.

Parameters

origin	Origin point (x,y) of the rayfire on mesh boundary
theta	Spherical polar angle (in radians)

Definition at line 51 of file rayfire_mesh.C.

References _theta, and GRINS::Constants::pi.

                                                                  :
    _dim(2),
    _origin(origin),
    _theta(theta),
    _phi(-7.0) // bounds on angles are +/- 2pi
  {
    if (std::abs(_theta) > 2.0*Constants::pi)
      libmesh_error_msg("Please supply a theta value between -2*pi and 2*pi");
  }

GRINS::RayfireMesh::RayfireMesh	(	libMesh::Point &	origin,
		libMesh::Real	theta,
		libMesh::Real	phi
	)

3D Constructor

This is restricted to 3D meshes. The init() function must be called to perform the actual rayfire.

Parameters

origin	Origin point (x,y,z) of the rayfire on mesh boundary
theta	Spherical polar angle (in radians)
phi	Spherical azimuthal angle (in radians)

Definition at line 41 of file rayfire_mesh.C.

                                                                                   :
    _dim(3),
    _origin(origin),
    _theta(theta),
    _phi(phi)
  {
    libmesh_not_implemented();
  }

GRINS::RayfireMesh::RayfireMesh	(	const GetPot &	input,
		const std::string &	qoi_string
	)

Input File Constructor.

Creates rayfire from parameters specified in input file section: QoI/'qoi_string'/Rayfire/

Definition at line 61 of file rayfire_mesh.C.

References _origin, _theta, and GRINS::Constants::pi.

                                                                             :
    _dim(2),
    _phi(-7.0)
  {
    unsigned int rayfire_dim = input.vector_variable_size("QoI/"+qoi_string+"/Rayfire/origin");

    if (rayfire_dim != 2)
      libmesh_error_msg("ERROR: Only 2D Rayfires are currently supported");

    if (!input.have_variable("QoI/"+qoi_string+"/Rayfire/origin"))
      libmesh_error_msg("ERROR: No origin specified for Rayfire");

    if (input.vector_variable_size("QoI/"+qoi_string+"/Rayfire/origin") != 2)
      libmesh_error_msg("ERROR: Please specify a 2D point (x,y) for the rayfire origin");

    libMesh::Point origin;
    _origin(0) = input("QoI/"+qoi_string+"/Rayfire/origin", 0.0, 0);
    _origin(1) = input("QoI/"+qoi_string+"/Rayfire/origin", 0.0, 1);

    if (input.have_variable("QoI/"+qoi_string+"/Rayfire/theta"))
      _theta = input("QoI/"+qoi_string+"/Rayfire/theta", -7.0);
    else
      libmesh_error_msg("ERROR: Spherical polar angle theta must be given for Rayfire");

    if (std::abs(_theta) > 2.0*Constants::pi)
      libmesh_error_msg("Please supply a theta value between -2*pi and 2*pi");

    if (input.have_variable("QoI/"+qoi_string+"/Rayfire/phi"))
      libmesh_error_msg("ERROR: cannot specify spherical azimuthal angle phi for Rayfire, only 2D is currently supported");
  }

GRINS::RayfireMesh::RayfireMesh

(

const RayfireMesh &

original

)

Copy Constructor.

Required to deep-copy _mesh and _elem_id_map

Definition at line 92 of file rayfire_mesh.C.

References _elem_id_map, and _mesh.

                                                       :
    _dim(original._dim),
    _origin(original._origin),
    _theta(original._theta),
    _phi(original._phi)
  {
    if (original._mesh.get())
      this->_mesh.reset( new libMesh::Mesh( *((original._mesh).get()) ) );

    this->_elem_id_map = original._elem_id_map;
  }

GRINS::RayfireMesh::RayfireMesh ( )

private

User should never call the default constructor.

Member Function Documentation

void GRINS::RayfireMesh::check_origin_on_boundary ( const libMesh::Elem *  start_elem )

private

Ensure the supplied origin is on a boundary of the mesh.

Definition at line 247 of file rayfire_mesh.C.

References _origin.

Referenced by init().

  {
    libmesh_assert(start_elem);

    // first, make sure the elem is on a boundary
    if ( !(start_elem->on_boundary()) )
      libmesh_error_msg("The supplied origin point is not on a boundary element");

    // second, check all boundary sides of the elem
    // to see if one of them conatins the origin point
    bool valid = false;

    for (unsigned int s=0; s<start_elem->n_sides(); s++)
      {
        // neighbor_ptr() returns NULL on boundary elems
        if ( start_elem->neighbor_ptr(s) )
          continue;

        // we found a boundary elem, so make an edge and see if it contains the origin
        libMesh::UniquePtr<const libMesh::Elem> edge_elem = start_elem->build_edge_ptr(s);
        valid |= edge_elem->contains_point(_origin);
      }

    if (!valid)
      libmesh_error_msg("The supplied origin point is not on the boundary of the starting element");
  }

bool GRINS::RayfireMesh::check_valid_point ( libMesh::Point &  intersection_point, libMesh::Point &  start_point, const libMesh::Elem &  edge_elem, libMesh::Point &  next_point  )

private

Ensure the calculated intersection point is on the edge_elem and is not the start_point.

Definition at line 355 of file rayfire_mesh.C.

Referenced by get_next_elem().

  {
    bool is_not_start = !(intersection_point.absolute_fuzzy_equals(start_point));
    bool is_on_edge = edge_elem.contains_point(intersection_point);

    bool is_valid = is_not_start && is_on_edge;

    if ( is_valid )
      {
        next_point(0) = intersection_point(0);
        next_point(1) = intersection_point(1);
      }

    return is_valid;
  }

void GRINS::RayfireMesh::coarsen ( const libMesh::Elem *  rayfire_elem )

private

Coarsening of a rayfire element whose main mesh counterpart was coarsened.

Definition at line 731 of file rayfire_mesh.C.

References _elem_id_map, _mesh, _origin, and get_rayfire_elem().

Referenced by reinit().

  {
    libmesh_assert(child_elem);

    if (this->get_rayfire_elem(child_elem->id()))
      {
        const libMesh::Elem* parent_elem = child_elem->parent();
        libmesh_assert(parent_elem);

        const libMesh::Node* start_node = NULL;
        const libMesh::Node* end_node = NULL;

        for (unsigned int c=0; c<parent_elem->n_children(); c++)
          {
            libMesh::Elem* rayfire_child = this->get_rayfire_elem(parent_elem->child_ptr(c)->id());

            if (rayfire_child)
              {
                if (!start_node)
                  {
                    start_node = rayfire_child->node_ptr(0);
                    end_node = rayfire_child->node_ptr(1);
                  }
                else
                  {
                    if ( (this->_origin - *(rayfire_child->node_ptr(0))).norm() < (this->_origin - *start_node).norm() )
                      start_node = rayfire_child->node_ptr(0);

                    if ( (this->_origin - *(rayfire_child->node_ptr(1))).norm() > (this->_origin - *end_node).norm() )
                      end_node = rayfire_child->node_ptr(1);
                  }

                rayfire_child->set_refinement_flag(libMesh::Elem::RefinementState::INACTIVE);
              }
          } // for c

        // make sure we found nodes
        libmesh_assert(start_node);
        libmesh_assert(end_node);

        // add a new rayfire elem
        libMesh::Elem* elem = _mesh->add_elem(new libMesh::Edge2);
        elem->set_node(0) = _mesh->node_ptr(start_node->id());
        elem->set_node(1) = _mesh->node_ptr(end_node->id());

        libmesh_assert_less( (*(elem->node_ptr(0))-_origin).norm(),  (*(elem->node_ptr(1))-_origin).norm());

        // add new rayfire elem to the map
        _elem_id_map[parent_elem->id()] = elem->id();
      }

  }

void GRINS::RayfireMesh::elem_ids_in_rayfire

(

std::vector< libMesh::dof_id_type > &

id_vector

)

const

Returns a vector of elem IDs that are currently along the rayfire.

Definition at line 187 of file rayfire_mesh.C.

References _elem_id_map, and _mesh.

  {
    std::map<libMesh::dof_id_type,libMesh::dof_id_type>::const_iterator it = _elem_id_map.begin();
    for(; it != _elem_id_map.end(); it++)
      {
        if (_mesh->elem_ptr(it->second)->active())
          id_vector.push_back(it->first);
      }
  }

const libMesh::Elem * GRINS::RayfireMesh::get_correct_neighbor ( libMesh::Point &  start_point, libMesh::Point &  end_point, const libMesh::Elem *  cur_elem, unsigned int  side, bool  same_parent  )

private

Knowing the end_point, get the appropraite next elem along the path.

Definition at line 432 of file rayfire_mesh.C.

References rayfire_in_elem(), and validate_edge().

Referenced by get_next_elem().

  {
    libmesh_assert(cur_elem);
    libmesh_assert(cur_elem->active());

    const libMesh::Elem * neighbor = cur_elem->neighbor_ptr(side);

    // a valid neighbor is an ACTIVE elem or null
    bool is_valid = true;
    if (neighbor)
      is_valid = neighbor->active();

    if (is_valid)
      {
        // check if the intersection point is a vertex
        bool is_vertex = false;
        const libMesh::Node * vertex = NULL;
        for(unsigned int n=0; n<cur_elem->n_nodes(); n++)
          {
            if ((cur_elem->node_ptr(n))->absolute_fuzzy_equals(end_point))
              {
                is_vertex = true;
                vertex = cur_elem->node_ptr(n);
                break;
              }
          }

        if (is_vertex)
          {
            // rayfire goes through vertex

            // check elem neighbors first
            for (unsigned int s=0; s<cur_elem->n_sides(); s++)
              {
                libMesh::UniquePtr<const libMesh::Elem> side_elem = cur_elem->build_side_ptr(s);
                if (side_elem->contains_point(end_point))
                  {
                    const libMesh::Elem * neighbor = cur_elem->neighbor_ptr(s);
                    if (!neighbor)
                      continue;

                    if (same_parent)
                      {
                        if (neighbor->parent())
                          {
                            if (neighbor->parent()->id() != cur_elem->parent()->id())
                              continue;
                            else
                              if (this->rayfire_in_elem(end_point,neighbor))
                                if (this->validate_edge(start_point,end_point,side_elem.get(),neighbor))
                                  return neighbor;
                          }
                        else
                          continue;
                      }
                    else
                      if (this->rayfire_in_elem(end_point,neighbor))
                        if (this->validate_edge(start_point,end_point,side_elem.get(),neighbor))
                          return neighbor;
                  }
              }

            // next elem is not a neighbor,
            // so get all elems that share this vertex
            std::set<const libMesh::Elem*> elem_set;
            cur_elem->find_point_neighbors(*vertex,elem_set);
            std::set<const libMesh::Elem *>::const_iterator       it  = elem_set.begin();
            const std::set<const libMesh::Elem *>::const_iterator end = elem_set.end();

            // iterate over each elem
            for (; it != end; ++it)
              {
                const libMesh::Elem* elem = *it;

                if (elem == cur_elem) // skip the current elem
                  continue;

                if (this->rayfire_in_elem(end_point,elem))
                  return elem;
              }

          }
        else
          {
            // check if side is a boundary
            if( !(cur_elem->neighbor_ptr(side)) )
              return NULL;

            // not a vertex or boundary, so just get the elem on that side
            return cur_elem->neighbor_ptr(side);
          }
      } // if (is_valid)
    else
      {
        for (unsigned int c=0; c<neighbor->n_children(); c++)
          {
            const libMesh::Elem * child = neighbor->child_ptr(c);

            if (child->contains_point(end_point))
              if (this->rayfire_in_elem(end_point,child))
                return child;
          }
        // could not find a valid child elem, so we return NULL below
      }

    return NULL;
  }

const libMesh::Elem * GRINS::RayfireMesh::get_next_elem ( const libMesh::Elem *  cur_elem, libMesh::Point &  start_point, libMesh::Point &  end_point, bool  same_parent = false  )

private

Calculate the intersection point.

Parameters

[out]

end_point

The intersection point in (x,y) coordinates

Returns

Elem* if intersection is found

NULL if no intersection point found (i.e. start_point is on a boundary)

Definition at line 323 of file rayfire_mesh.C.

References check_valid_point(), get_correct_neighbor(), and newton_solve_intersection().

Referenced by init(), and refine().

  {
    libmesh_assert(cur_elem);

    libMesh::Point intersection_point;

    // loop over all sides of the elem and check each one for intersection
    for (unsigned int s=0; s<cur_elem->n_sides(); s++)
      {
        libMesh::UniquePtr<const libMesh::Elem> edge_elem = cur_elem->build_edge_ptr(s);

        // Using the default tol can cause a false positive when start_point is near a node,
        // causing this loop to skip over an otherwise valid edge to check
        if (edge_elem->contains_point(start_point,libMesh::TOLERANCE*0.1))
          continue;

        bool converged = this->newton_solve_intersection(start_point,edge_elem.get(),intersection_point);

        if (converged)
          {
            if ( this->check_valid_point(intersection_point,start_point,*edge_elem,next_point) )
              return this->get_correct_neighbor(start_point,intersection_point,cur_elem,s,same_parent);
          }
        else
          continue;

      } // for s

    return NULL; // no intersection
  }

libMesh::Elem * GRINS::RayfireMesh::get_rayfire_elem ( const libMesh::dof_id_type  elem_id )

private

Private function to get a rayfire elem from main_mesh elem ID.

Does not return a const pointer, and is used within this class to simplify rayfire elem access.

Also used by map_to_rayfire_elem(), which adds a const to prevent modification outside this class.

Definition at line 308 of file rayfire_mesh.C.

References _elem_id_map, and _mesh.

Referenced by coarsen(), and map_to_rayfire_elem().

  {
    // return value; set if valid rayfire elem is found
    libMesh::Elem* retval = NULL;

    std::map<libMesh::dof_id_type,libMesh::dof_id_type>::iterator it;
    it = _elem_id_map.find(elem_id);
    if (it != _elem_id_map.end())
      if (_mesh->elem_ptr(it->second)->refinement_flag() != libMesh::Elem::RefinementState::INACTIVE)
        retval = _mesh->elem_ptr(it->second);

    return retval;
  }

const libMesh::Elem * GRINS::RayfireMesh::get_start_elem ( const libMesh::MeshBase &  mesh_base )

private

Get the correct starting elem corresponding to _origin.

Returns

NULL Origin is not on the mesh

Elem* First elem on the rayfire

Definition at line 275 of file rayfire_mesh.C.

References _origin, and rayfire_in_elem().

Referenced by init().

  {
    const libMesh::Elem* start_elem = NULL;

    libMesh::UniquePtr<libMesh::PointLocatorBase> locator = mesh_base.sub_point_locator();
    const libMesh::Elem* elem = (*locator)(_origin);

    // elem would be NULL if origin is not on mesh
    if (elem)
      {
        if (this->rayfire_in_elem(_origin,elem))
          start_elem = elem;
        else
          {
            for (unsigned int i=0; i<elem->n_neighbors(); i++)
              {
                const libMesh::Elem* neighbor_elem = elem->neighbor_ptr(i);
                if (!neighbor_elem)
                  continue;

                if (this->rayfire_in_elem(_origin,neighbor_elem))
                  {
                    start_elem = neighbor_elem;
                    break;
                  }
              }
          }
      }

    return start_elem; // might be NULL if _origin is not on mesh
  }

void GRINS::RayfireMesh::init

(

const libMesh::MeshBase &

mesh_base

)

Initialization.

This function performs the rayfire and assembles the 1D mesh. Call this immediately after the constructor. Must be called before any refinements of the main mesh.

Parameters

mesh_base

Reference to the main mesh

Definition at line 104 of file rayfire_mesh.C.

References _dim, _elem_id_map, _mesh, _origin, check_origin_on_boundary(), get_next_elem(), and get_start_elem().

  {
    // consistency check
    if(mesh_base.mesh_dimension() != _dim)
      {
        std::stringstream ss;
        ss <<"The supplied mesh object is " <<mesh_base.mesh_dimension()
           <<"D, but the RayfireMesh object was created with the "
           <<_dim <<"D constructor";

        libmesh_error_msg(ss.str());
      }

    _mesh.reset( new libMesh::Mesh(mesh_base.comm(),(unsigned char)1) );

    libMesh::Point start_point(_origin);

    // get first element
    const libMesh::Elem* start_elem = this->get_start_elem(mesh_base);

    if (!start_elem)
      libmesh_error_msg("Origin is not on mesh");

    // ensure the origin is on a boundary element
    // AND on the boundary of said element
    this->check_origin_on_boundary(start_elem);

    // add the origin point to the point list
    libMesh::Node * start_node = _mesh->add_point(start_point);
    libMesh::Node * end_node = NULL;

    libMesh::Point end_point;

    const libMesh::Elem* next_elem;
    const libMesh::Elem* prev_elem = start_elem;

    do
      {
        // calculate the end point and
        // get the next elem in the rayfire
        next_elem = this->get_next_elem(prev_elem,start_point,end_point);

#ifndef NDEBUG
        // make sure we are only picking up active elements
        if (next_elem)
          libmesh_assert( next_elem->active() );
#endif

        // add end point as node on the rayfire mesh
        end_node = _mesh->add_point(end_point);
        libMesh::Elem* elem = _mesh->add_elem(new libMesh::Edge2);
        elem->set_node(0) = start_node;
        elem->set_node(1) = end_node;

        // warn if rayfire elem is shorter than TOLERANCE
        if ( (start_point-end_point).norm() < libMesh::TOLERANCE)
          {
            std::stringstream ss;
            ss  <<"********\n"
                <<"WARNING\n"
                <<"Detected rayfire element shorter than TOLERANCE\n"
                <<"Element ID: " <<prev_elem->id() <<", rayfire element length: " <<(start_point-end_point).norm();

            libmesh_warning(ss.str());
          }

        // add new rayfire elem to the map
        _elem_id_map[prev_elem->id()] = elem->id();

        start_point = end_point;
        start_node = end_node;
        prev_elem = next_elem;
      } while(next_elem);

  }

const libMesh::Elem * GRINS::RayfireMesh::map_to_rayfire_elem

(

const libMesh::dof_id_type

elem_id

)

This function takes in an elem_id on the main mesh and returns an elem from the 1D rayfire mesh.

Parameters

elem_id

The ID of the elem on the main mesh

Returns

Elem* to 1D elem from the rayfire mesh if the elem_id falls along the rayfire path

NULL if the given elem_id does not correspond to an elem through which the rayfire passes

Definition at line 181 of file rayfire_mesh.C.

References get_rayfire_elem().

  {
    return this->get_rayfire_elem(elem_id);
  }

bool GRINS::RayfireMesh::newton_solve_intersection ( libMesh::Point &  initial_point, const libMesh::Elem *  edge_elem, libMesh::Point &  intersection_point  )

private

Iterative solver for calculating the intersection point of the rayfire on the side of an elem.

The rayfire line can be represented in point-slope form:

y-y0 = m(x-x0)

where m is the slope and (x0,y0) is the known initial_point. To find the intersection point, we will use a parametric representation of the edge utilizing the shape functions as such:

X(ξ) = sum( x_j φ_j(ξ) )

Y(ξ) = sum( y_j φ_j(ξ) )

where x_j and y_j are the node x- and y-coordinates, respectively, and φ_j is the value of shape function j at reference coordinate ξ

The intersection point is then where the x and y coordinates from the point-slope equation equal the X and Y coordinates from the parametric equations for the given edge, respectively. We can then form a residual by substituting the parametric coordinates into the point slope form and bringing all terms to one side of the equation as:

f = 0 = m(X(ξ)-x0) - (Y(ξ)-y0)

The residual is then a function of a single parameter, namely the reference coordinate ξ. The derivative is rather trivial and can be expressed as

df = m*dX - dY

where

dX = sum( x_j dφ_j/dξ )

dY = sum( y_j dφ_j/dξ )

Parameters

[out]

intersection_point

The calculated intersection point (only set if convergence is achieved)

Returns

Whether or not the solver converged before hitting the iteration limit

Definition at line 541 of file rayfire_mesh.C.

References _theta.

Referenced by get_next_elem().

  {
    libmesh_assert(edge_elem);

    unsigned int iter_max = 20; // max iterations

    // the number of shape functions needed for the edge_elem
    unsigned int n_sf = libMesh::FE<1,libMesh::LAGRANGE>::n_shape_functions(edge_elem->type(),edge_elem->default_order());

    // starting point on the elem
    libMesh::Real x0 = initial_point(0);
    libMesh::Real y0 = initial_point(1);

    // shape functions and derivatives w.r.t reference coordinate
    std::vector<libMesh::Real> phi(n_sf);
    std::vector<libMesh::Real> dphi(n_sf);

    // Newton iteration step
    libMesh::Real d_xi;

    // tan(theta) is the slope, so precompute since it is used repeatedly
    libMesh::Real tan_theta = std::tan(_theta);

    // Initial guess is center of the edge_elem
    libMesh::Real xi = 0.0;

    // Newton iteration
    for(unsigned int it=0; it<iter_max; it++)
      {
        // Get the shape function and derivative values at the reference coordinate
        // phi.size() == dphi.size()
        for(unsigned int i=0; i<phi.size(); i++)
          {
            phi[i] = libMesh::FE<1,libMesh::LAGRANGE>::shape(edge_elem->type(),
                                                             edge_elem->default_order(),
                                                             i,
                                                             xi);

            dphi[i] = libMesh::FE<1,libMesh::LAGRANGE>::shape_deriv(edge_elem->type(),
                                                                    edge_elem->default_order(),
                                                                    i,
                                                                    0, // const unsigned int libmesh_dbg_varj
                                                                    xi);
          } // for i

        libMesh::Real X=0.0, Y=0.0, dX=0.0, dY=0.0;

        for(unsigned int i=0; i<phi.size(); i++)
          {
            X  += (*(edge_elem->node_ptr(i)))(0) * phi[i];
            dX += (*(edge_elem->node_ptr(i)))(0) * dphi[i];
            Y  += (*(edge_elem->node_ptr(i)))(1) * phi[i];
            dY += (*(edge_elem->node_ptr(i)))(1) * dphi[i];
          }

        libMesh::Real  f = tan_theta*(X-x0) - (Y-y0);
        libMesh::Real df = tan_theta*(dX) - dY;

        d_xi = f/df;

        if(std::abs(d_xi) < libMesh::TOLERANCE)
          {
            // convergence
            intersection_point(0) = X;
            intersection_point(1) = Y;
            return true;
          }
        else
          xi -= d_xi;

      } // for it

    // no convergence
    return false;
  }

bool GRINS::RayfireMesh::rayfire_in_elem ( const libMesh::Point &  end_point, const libMesh::Elem *  elem  )

private

Walks a short distance along the rayfire and checks if elem contains that point.

Definition at line 372 of file rayfire_mesh.C.

References _theta.

Referenced by get_correct_neighbor(), get_start_elem(), and refine().

  {
    libmesh_assert(elem);

    // move a little bit along the rayfire and see if we are still in the elem
    // need to move more than TOLERANCE to avoid false positive on contains_point()
    libMesh::Real L = 2*libMesh::TOLERANCE;

    // If the elem is too small, need to shorten L so we stay within the elem
    if ( elem->hmin() < libMesh::TOLERANCE )
      L = elem->hmin() * 0.1;

    // parametric representation of rayfire line
    libMesh::Real x = end_point(0) + L*std::cos(_theta);
    libMesh::Real y = end_point(1) + L*std::sin(_theta);

    return elem->contains_point(libMesh::Point(x,y));
  }

void GRINS::RayfireMesh::refine ( const libMesh::Elem *  main_elem, libMesh::Elem *  rayfire_elem  )

private

Refinement of a rayfire element whose main mesh counterpart was refined.

Definition at line 618 of file rayfire_mesh.C.

References _elem_id_map, _mesh, _origin, get_next_elem(), and rayfire_in_elem().

Referenced by reinit().

  {
    libmesh_assert(main_elem);
    libmesh_assert(rayfire_elem);

    libmesh_assert_equal_to(main_elem->refinement_flag(),libMesh::Elem::RefinementState::INACTIVE);

    // these nodes cannot change
    libMesh::Node* start_node = rayfire_elem->node_ptr(0);
    libMesh::Node* end_node   = rayfire_elem->node_ptr(1);

    // set the rayfire_elem as INACTIVE
    rayfire_elem->set_refinement_flag(libMesh::Elem::RefinementState::INACTIVE);

    // find which child elem we start with
    libMesh::dof_id_type start_child = libMesh::invalid_uint;
    for(unsigned int i=0; i<main_elem->n_children(); i++)
      {
        if ( (main_elem->child_ptr(i))->contains_point(*start_node) )
          {
            // check if the start point is a vertex
            bool is_vertex = false;
            for(unsigned int n=0; n<main_elem->child_ptr(i)->n_nodes(); n++)
              {
                if ((main_elem->child_ptr(i)->node_ptr(n))->absolute_fuzzy_equals(*start_node))
                  {
                    is_vertex = true;
                    break;
                  }
              }

            if (is_vertex)
              {
                if ( this->rayfire_in_elem(*start_node, main_elem->child_ptr(i)) )
                  {
                    start_child = i;
                    break;
                  }
              }
            else
              {
                start_child = i;
                break;
              }
          }
      }

    // make sure we found a starting child elem
    libmesh_assert_not_equal_to(start_child,libMesh::invalid_uint);

    // we found the starting element, so perform the rayfire
    // until we reach the stored end_node
    libMesh::Point start_point = *start_node;
    libMesh::Point end_point;

    const libMesh::Elem* next_elem;
    const libMesh::Elem* prev_elem = main_elem->child_ptr(start_child);

    // if prev_elem is INACTIVE, then more than one refinement
    // has taken place between reinit() calls and will
    // break this
    libmesh_assert_equal_to( prev_elem->refinement_flag(), libMesh::Elem::RefinementState::JUST_REFINED );

    libMesh::Node * prev_node = start_node;
    libMesh::Node * new_node = NULL;

    // iterate until we reach the stored end_node
    do
      {
        next_elem = this->get_next_elem(prev_elem,start_point,end_point,true);

#ifndef NDEBUG
        // make sure we are only picking up active elements
        if (next_elem)
          libmesh_assert( next_elem->active() );
#endif

        // add end point as node on the rayfire mesh
        new_node = _mesh->add_point(end_point);
        libMesh::Elem* elem = _mesh->add_elem(new libMesh::Edge2);

        elem->set_node(0) = prev_node;
        elem->set_node(1) = new_node;

        libmesh_assert_less( (*(elem->node_ptr(0))-_origin).norm(),  (*(elem->node_ptr(1))-_origin).norm());

        // warn if rayfire elem is shorter than TOLERANCE
        if ( (start_point-end_point).norm() < libMesh::TOLERANCE)
          {
            std::stringstream ss;
            ss  <<"********\n"
                <<"WARNING\n"
                <<"Detected rayfire element shorter than TOLERANCE on refinement\n"
                <<"Element ID: " <<prev_elem->id() <<", rayfire element length: " <<(start_point-end_point).norm();

            libmesh_warning(ss.str());
          }

        // set rayfire_elem as the parent of this new elem
        // in case it gets coarsened
        elem->set_parent(rayfire_elem);

        // add new rayfire elem to the map
        _elem_id_map[prev_elem->id()] = elem->id();
        start_point = end_point;
        prev_elem = next_elem;
        prev_node = new_node;

      } while(!(new_node->absolute_fuzzy_equals(*end_node)));

  }

void GRINS::RayfireMesh::reinit

(

const libMesh::MeshBase &

mesh_base

)

Checks for refined and coarsened main mesh elements along the rayfire path.

They are then passed to refine() and coarsen(), respectively, to update the rayfire mesh.

Only 1 level of refinement and/or coarsening can be done between reinit() calls. Note that this is not limited to doing either refinement or coarsening between reinits. Both can be done on different elements, as long as they are only 1 level in either direction.

Parameters

mesh	reference to main mesh, needed to get Elem* from the stored elem_id's

Definition at line 198 of file rayfire_mesh.C.

References _elem_id_map, _mesh, coarsen(), and refine().

  {
    // store the elems to be refined until later
    // so we don't mess with the _elem_id_map while we
    // iterate over it
    std::vector<std::pair<const libMesh::Elem*,libMesh::Elem*> > elems_to_refine;

    // same with elems to coarsen
    // store the main_mesh elem
    std::vector<const libMesh::Elem*> elems_to_coarsen;

    // iterate over all main elems along the rayfire and look for
    // refinement: INACTIVE parent with JUST_REFINED children
    // coarsening: JUST_COARSENED parent
    std::map<libMesh::dof_id_type,libMesh::dof_id_type>::iterator it = _elem_id_map.begin();
    for(; it != _elem_id_map.end(); it++)
      {
        const libMesh::Elem* main_elem = mesh_base.elem_ptr(it->first);
        libmesh_assert(main_elem);

        if (main_elem->parent())
          {
            if (main_elem->parent()->refinement_flag() == libMesh::Elem::RefinementState::JUST_COARSENED)
              elems_to_coarsen.push_back(main_elem);
          }

        libMesh::Elem::RefinementState state = main_elem->refinement_flag();

        if (state == libMesh::Elem::RefinementState::INACTIVE)
          {
            if (main_elem->has_children())
              if (main_elem->child_ptr(0)->refinement_flag() == libMesh::Elem::RefinementState::JUST_REFINED)
                elems_to_refine.push_back(std::pair<const libMesh::Elem*, libMesh::Elem*>(main_elem,_mesh->elem_ptr(it->second)));
          }
      }

    // refine the elements that need it
    for (unsigned int i=0; i<elems_to_refine.size(); i++)
      this->refine(elems_to_refine[i].first, elems_to_refine[i].second);

    // coarsen the elements that need it
    for (unsigned int i=0; i<elems_to_coarsen.size(); i++)
      this->coarsen(elems_to_coarsen[i]);

  }

bool GRINS::RayfireMesh::validate_edge ( const libMesh::Point &  start_point, const libMesh::Point &  end_point, const libMesh::Elem *  side_elem, const libMesh::Elem *  neighbor  )

private

Ensures the rayfire doesn't wander into the middle of an elem.

Definition at line 392 of file rayfire_mesh.C.

Referenced by get_correct_neighbor().

  {
    bool is_valid = true;

    const libMesh::Node * node0 = side_elem->node_ptr(0);
    const libMesh::Node * node1 = side_elem->node_ptr(1);

    unsigned int side = libMesh::invalid_uint;

    for (unsigned int s=0; s<neighbor->n_sides(); s++)
      {
        libMesh::UniquePtr<const libMesh::Elem> side_elem = neighbor->build_side_ptr(s);

        if ( (side_elem->contains_point(*node0)) && (side_elem->contains_point(*node1)) )
          {
            side = s;
            break;
          }
      }

    if ( side != libMesh::invalid_uint )
      {
        libMesh::UniquePtr<const libMesh::Elem> edge_elem = neighbor->build_side_ptr(side);

        bool start_point_on_edge = edge_elem->contains_point(start_point);
        bool end_point_on_edge = edge_elem->contains_point(end_point);

        if ( end_point_on_edge && start_point_on_edge )
          {
            bool l_to_r = ( start_point.absolute_fuzzy_equals(*(edge_elem->node_ptr(0))) ) && ( end_point.absolute_fuzzy_equals(*(edge_elem->node_ptr(1))) );
            bool r_to_l = ( end_point.absolute_fuzzy_equals(*(edge_elem->node_ptr(0))) )   && ( start_point.absolute_fuzzy_equals(*(edge_elem->node_ptr(1))) );

            is_valid &= ( l_to_r || r_to_l );
          }
      }

    return is_valid;
  }

Member Data Documentation

const unsigned int GRINS::RayfireMesh::_dim

private

Dimension of the main mesh.

Definition at line 136 of file rayfire_mesh.h.

Referenced by init().

std::map<libMesh::dof_id_type,libMesh::dof_id_type> GRINS::RayfireMesh::_elem_id_map

private

Map of main mesh elem_id to rayfire mesh elem_id.

Definition at line 151 of file rayfire_mesh.h.

Referenced by coarsen(), elem_ids_in_rayfire(), get_rayfire_elem(), init(), RayfireMesh(), refine(), and reinit().

libMesh::UniquePtr<libMesh::Mesh> GRINS::RayfireMesh::_mesh

private

Internal 1D mesh of EDGE2 elements.

Definition at line 148 of file rayfire_mesh.h.

Referenced by coarsen(), elem_ids_in_rayfire(), get_rayfire_elem(), init(), RayfireMesh(), refine(), and reinit().

libMesh::Point GRINS::RayfireMesh::_origin

private

Origin point.

Definition at line 139 of file rayfire_mesh.h.

Referenced by check_origin_on_boundary(), coarsen(), get_start_elem(), init(), RayfireMesh(), and refine().

libMesh::Real GRINS::RayfireMesh::_phi

private

Rayfire Spherical azimuthal angle (in radians)

Definition at line 145 of file rayfire_mesh.h.

libMesh::Real GRINS::RayfireMesh::_theta

private

Rayfire Spherical polar angle (in radians)

Definition at line 142 of file rayfire_mesh.h.

Referenced by newton_solve_intersection(), rayfire_in_elem(), and RayfireMesh().

---

The documentation for this class was generated from the following files:

• src/qoi/include/grins/rayfire_mesh.h
• src/qoi/src/rayfire_mesh.C