#include <distance_function.h>

Inheritance diagram for GRINS::DistanceFunction:

Collaboration diagram for GRINS::DistanceFunction:

Public Member Functions
	DistanceFunction (libMesh::EquationSystems &es_in, const libMesh::UnstructuredMesh &bm_in)

	~DistanceFunction ()

virtual void	initialize ()

libMesh::Real	node_to_boundary (const libMesh::Node *node)

void	compute ()

libMesh::AutoPtr< libMesh::DenseVector< libMesh::Real > >	interpolate (const libMesh::Elem *elem, const std::vector< libMesh::Point > &qts) const

Private Attributes
libMesh::EquationSystems &	_equation_systems

const libMesh::UnstructuredMesh &	_boundary_mesh

libMesh::AutoPtr< libMesh::FEBase >	_dist_fe

Detailed Description

Definition at line 53 of file distance_function.h.

Constructor & Destructor Documentation

GRINS::DistanceFunction::DistanceFunction	(	libMesh::EquationSystems &	es_in,
		const libMesh::UnstructuredMesh &	bm_in
	)

Constructor

Definition at line 148 of file distance_function.C.

References _equation_systems.

                                                                                                         :
     _equation_systems (es_in),
     _boundary_mesh    (bm_in),
     _dist_fe          (libMesh::FEBase::build(_equation_systems.get_mesh().mesh_dimension(), libMesh::FEType(libMesh::FIRST, libMesh::LAGRANGE)))
   {
     // Ensure that libmesh is ready to roll
     libmesh_assert(libMesh::initialized());
 
     // Add distance function system
     _equation_systems.add_system<libMesh::System>("distance_function");
 
     // Get reference to distance function system we just added
     libMesh::System& sys = _equation_systems.get_system<libMesh::System>("distance_function");
 
     // Add distance function variable
     sys.add_variable("distance", libMesh::FIRST);
 
     // Attach initialization function
     sys.attach_init_object(*this);
 
   }

GRINS::DistanceFunction::~DistanceFunction ( )

inline

Destructor

Definition at line 65 of file distance_function.h.

65 {};

Member Function Documentation

void GRINS::DistanceFunction::compute ( )

Initialize "distance_function" equation system by computing distance from each mesh node to the nearest point in boundary_mesh

Definition at line 557 of file distance_function.C.

References _boundary_mesh, _equation_systems, and node_to_boundary().

Referenced by initialize().

   {
     // Get mesh
     const libMesh::MeshBase& mesh = _equation_systems.get_mesh();
 
     // Get reference to system and system number
     libMesh::System& system = _equation_systems.get_system<libMesh::System>("distance_function");
     const unsigned int sys_num = system.number();
 
     //Debugging code
     //system.print_info();
     //_equation_systems.print_info();
     //End debugging
 
     // The boundary mesh needs to all be on this processor for us to
     // calculate a correct distance function.  Since we don't need it to
     // be serial afterwards, we use a temporary serializer.
     {
       libMesh::MeshSerializer serialize(const_cast<libMesh::UnstructuredMesh&>(_boundary_mesh));
 
       // Loop over nodes in mesh
       libMesh::MeshBase::const_node_iterator node_it  = mesh.local_nodes_begin();
       const libMesh::MeshBase::const_node_iterator node_end = mesh.local_nodes_end();
 
       //Debugging code
       //const libMesh::Elem *elem = *mesh.elements_begin();
       //if (elem->default_order() == libMesh::FIRST)
       // std::cout<<"Mesh is using FIRST order elements."<<std::endl;
       //End debugging
 
 
       for ( ; node_it != node_end; ++node_it) {
 
         // Grab node
         const libMesh::Node* node = *node_it;
 
         // Compute distance to nearest point in boundary_mesh
         const libMesh::Real distance = DistanceFunction::node_to_boundary (node);
 
         // Stuff data into appropriate place in the system solution
         const unsigned int dof = node->dof_number(sys_num,0,0);
         system.solution->set (dof, distance);
 
       } // end loop over nodes
 
     } // end boundary mesh serialization
 
     system.solution->close();
     system.update();
 
     std::cout << "Distance Function computed." << std::endl;
   }

void GRINS::DistanceFunction::initialize ( )

virtual

Initializes the distance

Definition at line 173 of file distance_function.C.

References compute().

   {
     // Call the compute function
     this->compute();
     return;
   }

libMesh::AutoPtr< libMesh::DenseVector< libMesh::Real > > GRINS::DistanceFunction::interpolate	(	const libMesh::Elem *	elem,
		const std::vector< libMesh::Point > &	qts
	)		const

Interpolate distance function to points qpts (in reference space) for element *elem

Definition at line 615 of file distance_function.C.

References _dist_fe, and _equation_systems.

   {
     libmesh_assert( elem != NULL );    // can't interpolate in NULL elem
     libmesh_assert( qpts.size() > 0 ); // can't interpolate if no points requested
 
     // grab basis functions (evaluated at qpts)
     const std::vector<std::vector<libMesh::Real> > &phi = _dist_fe->get_phi();
 
     // reinitialize finite element data at qpts
     _dist_fe->reinit(elem, &qpts);
 
     // number of basis functions
     const unsigned int n_dofs = phi.size();
 
     // number of points
     const unsigned int n_pts = qpts.size();
 
     // instantiate auto_ptr to dense vector to hold results
     libMesh::AutoPtr< DenseVector<libMesh::Real> > ap( new libMesh::DenseVector<libMesh::Real>(qpts.size()) );
     (*ap).zero();
 
     // pull off distance function at nodes on this element
     libMesh::System& sys = _equation_systems.get_system<libMesh::System>("distance_function");
     const libMesh::DofMap& dof_map = sys.get_dof_map();
 
     std::vector<unsigned int> dof_ind;
     dof_map.dof_indices(elem, dof_ind);
 
     libMesh::DenseVector<libMesh::Real> nodal_dist;
     nodal_dist.resize(n_dofs);
     //dof_map.extract_local_vector( *(sys.solution), dof_ind, nodal_dist);
     dof_map.extract_local_vector( *(sys.current_local_solution), dof_ind, nodal_dist);
 
     for ( unsigned int idof=0; idof<n_dofs; idof++ ) {
       for ( unsigned int iqpt=0; iqpt<n_pts; iqpt++ ) {
         (*ap)(iqpt) += nodal_dist(idof) * phi[idof][iqpt];
       }
     }
 
     return ap;
   }

libMesh::Real GRINS::DistanceFunction::node_to_boundary ( const libMesh::Node * node )

Compute distance from input node to boundary_mesh

Definition at line 184 of file distance_function.C.

References _boundary_mesh, and _equation_systems.

Referenced by compute().

   {
     // Ensure that node is not NULL
     libmesh_assert( node != NULL );
 
     // Initialize distance to infinity
     libMesh::Real distance = std::numeric_limits<libMesh::Real>::infinity();
 
     // Get dimension
     const unsigned int dim = _equation_systems.get_mesh().mesh_dimension();
     libmesh_assert( (dim==2) || (dim==3) );
 
     // Get coordinates of node
     std::vector<libMesh::Real> xnode(dim);
     for( unsigned int ii=0; ii<dim; ii++ ) xnode[ii] = (*node)(ii);
 
     // This function will work on a distributed interior mesh, but won't
     // give correct results on a distributed boundary mesh.
     libmesh_assert(_boundary_mesh.is_serial());
 
     // First loop over the boundary mesh and compute the distance to the
     // nodes.  This will give us an idea of what elements to check more
     // closely.
 
     libMesh::MeshBase::const_element_iterator       el     = _boundary_mesh.active_elements_begin();
     const libMesh::MeshBase::const_element_iterator end_el = _boundary_mesh.active_elements_end();
 
     if (el==end_el)
       {
         std::cout << "There are no elements to iterate through!!!  Returning..." << std::endl;
         return distance;
       }
 
     const libMesh::Elem* elem_containing_min = NULL;
     libMesh::Real dmin = std::numeric_limits<Real>::infinity();
 
     for ( ; el != end_el; ++el) {
 
       const libMesh::Elem* belem = *el;
 
       for (unsigned int bnode=0; bnode<belem->n_nodes(); ++bnode)
         {
           const libMesh::Point& pbndry = belem->point(bnode);
 
           libMesh::Real dnode = 0.0;
           for (unsigned int idim=0; idim<dim; ++idim)
             dnode += (xnode[idim] - pbndry(idim))*(xnode[idim] - pbndry(idim));
 
           dnode = std::sqrt(dnode);
 
           if (dnode<dmin)
             {
               dmin = dnode;
               elem_containing_min = belem;
             }
         }
 
     } // end element loop
 
     // grab the elements around the minimum
     std::set<const libMesh::Elem*> near_min_elems;
     elem_containing_min->find_point_neighbors (near_min_elems);
 
     // insert the element itself into the set... we must check it also
     near_min_elems.insert(elem_containing_min);
 
     std::set<const libMesh::Elem*>::iterator nm_el;
 
     // loop over the near_min_elems
     for (nm_el=near_min_elems.begin(); nm_el!=near_min_elems.end(); ++nm_el)
       {
 
         const libMesh::Elem* belem = *nm_el; //*el;
 
         // Ensure that elem defined by edge/face is linear
         libmesh_assert( belem->default_order() == FIRST );
         // For now we assume that the boundary_elems are linear, but give a warning
         //libmesh_warning( "Moving ahead, but need to make sure that the boundary elems are linear");
 
         // Initialize distance to this edge/face to infinity
         libMesh::Real dedge = std::numeric_limits<libMesh::Real>::infinity();
 
         if ( dim==2 )
           { // 2d
 
             // For 2d mesh, boundary elems must be 1d.  Due to the way
             // find_point_neighbors works, there may be some 2d elements
             // in the near_min_elems set.  These elements are clearly
             // not part of the boundary mesh, and thus we skip them.
             //
             if (belem->dim()!=1) continue;
 
             //std::cout<<"Number of nodes: "<<belem->n_nodes()<<std::endl;
 
             libmesh_assert( belem->n_nodes() == 2 );
             libmesh_assert( belem->type() == EDGE2 );
 
             // Points defining the edge
             const libMesh::Point& p0 = belem->point(0);
             const libMesh::Point& p1 = belem->point(1);
 
             Line line;
             line.p0 = p0;
             line.p1 = p1;
 
             DistanceToSegment<2>(*node, line, dedge);
 
             if ( dedge < distance ) distance = dedge;
 
           }
         else
           { // 3d
 
             // For 3d mesh, boundary elems must be 2d.  Due to the way
             // find_point_neighbors works, there may be some 3d elements
             // in the near_min_elems set.  These elements are clearly
             // not part of the boundary mesh, and thus we skip them.
             //
             if (belem->dim()!=2) continue;
 
             libmesh_assert( (belem->type() == TRI3) || (belem->type() == QUAD4) );
 
             if ( belem->type() == TRI3 )
               { // triangular boundary faces
 
                 // Find the point in the plane defined by the boundary nodes that
                 // minimizes the distance between the plane and the input node.
                 //
                 // Done in terms of the reference coordinates of the boundary element.
                 //
                 const libMesh::Point& p0 = belem->point(0);
                 const libMesh::Point& p1 = belem->point(1);
                 const libMesh::Point& p2 = belem->point(2);
 
                 const libMesh::Real x_xi = (p1(0) - p0(0)), x_et = (p2(0) - p0(0));
                 const libMesh::Real y_xi = (p1(1) - p0(1)), y_et = (p2(1) - p0(1));
                 const libMesh::Real z_xi = (p1(2) - p0(2)), z_et = (p2(2) - p0(2));
 
                 libMesh::Real A[2][2], b[2];
 
                 A[0][0] = x_xi*x_xi + y_xi*y_xi + z_xi*z_xi;
                 A[0][1] = x_xi*x_et + y_xi*y_et + z_xi*z_et;
                 A[1][0] = x_xi*x_et + y_xi*y_et + z_xi*z_et;
                 A[1][1] = x_et*x_et + y_et*y_et + z_et*z_et;
 
                 b[0] = (xnode[0] - p0(0))*x_xi + (xnode[1] - p0(1))*y_xi + (xnode[2] - p0(2))*z_xi;
                 b[1] = (xnode[0] - p0(0))*x_et + (xnode[1] - p0(1))*y_et + (xnode[2] - p0(2))*z_et;
 
                 const libMesh::Real detA = A[0][0]*A[1][1] - A[1][0]*A[0][1];
                 libmesh_assert( fabs(detA) > 0.0 ); // assert that A is not singular
 
                 const libMesh::Real xi = ( A[1][1]*b[0] - A[0][1]*b[1])/detA;
                 const libMesh::Real et = (-A[1][0]*b[0] + A[0][0]*b[1])/detA;
 
 
                 // If projection of node onto plane defined by boundary face
                 // (i.e., what we just computed) is not inside boundary face
                 // then we need to figure out closest point that is inside the
                 // boundary face
                 //
                 if ( (xi<0.0) || (et<0.0) || (et>1.0-xi) ) {
 
                   libMesh::Real dtmp = std::numeric_limits<libMesh::Real>::infinity();
 
                   // for each edge of boundary face, find distance from
                   // input node to the segment defined by that edge
                   //
                   // the minimum of these distances minimizes the distance
                   // to the node over the points in the boundary face
                   //
                   for ( unsigned int iedge=0; iedge<3; iedge++ ){
 
                     Line line;
 
                     // Get the endpoints of this edge
                     if      (iedge==0) { line.p0 = p1; line.p1 = p2; }
                     else if (iedge==1) { line.p0 = p0; line.p1 = p2; }
                     else   /*iedge==2*/{ line.p0 = p0; line.p1 = p1; }
 
                     DistanceToSegment<3>(*node, line, dtmp);
 
                     if( dtmp < dedge ) dedge = dtmp;
 
                   }
 
                 } else { // projection is inside face, so we're good to go
 
                   // Map from reference to physical space
                   libMesh::Real xint[3];
                   for ( unsigned int ii=0; ii<3; ii++ ) {
                     xint[ii] = p0(ii)*(1.0 - xi - et) + p1(ii)*xi + p2(ii)*et;
                   }
 
                   // compute distance
                   dedge = 0.0;
                   for ( unsigned int ii=0; ii<3; ii++ ) {
                     dedge += (xnode[ii] - xint[ii])*(xnode[ii] - xint[ii]);
                   }
                   dedge = sqrt(dedge);
 
                 }
 
               }
             else if ( belem->type() == QUAD4 )
               {
                 //std::cout << "WARNING: this functionality is not well-tested.  Sorry." << std::endl;
 
                 libMesh::Real RTOL = 1e-10;
                 libMesh::Real ATOL = 1e-20;
                 const unsigned int ITER_MAX=1000;
                 unsigned int iter=0;
 
                 ComputeDistanceResidual res(belem, node);
 
                 ComputeDistanceJacobian jac;
 
                 libMesh::DenseVector<libMesh::Real> X(2); X(0) = X(1) = 0.0;
                 libMesh::DenseVector<libMesh::Real> dX(2);
                 libMesh::Real det;
 
                 libMesh::DenseVector<libMesh::Real> R(2);
                 libMesh::DenseMatrix<libMesh::Real> dRdX(2,2);
                 libMesh::DenseMatrix<libMesh::Real> dRdXinv(2,2);
 
                 // evaluate residual... maybe we don't have to iterate
                 res(X,R);
 
                 libMesh::Real R0 = R.l2_norm();
 
                 // iterate
                 while ( (R.l2_norm() > ATOL) && (R.l2_norm()/R0 > RTOL) && (iter<ITER_MAX) )
                   {
                     // compute jacobian
                     jac(X, R, res, dRdX);
 
                     // invert jacobian
                     det = dRdX(0,0)*dRdX(1,1) - dRdX(0,1)*dRdX(1,0);
 
                     // protect against divide by zero
                     if(std::abs(det)<=1e-20)
                       {
                         std::cout << "WARNING: about to divide by zero." << std::endl;
                       }
 
                     dRdXinv(0,0) =  dRdX(1,1)/det;
                     dRdXinv(0,1) = -dRdX(0,1)/det;
                     dRdXinv(1,0) = -dRdX(1,0)/det;
                     dRdXinv(1,1) =  dRdX(0,0)/det;
 
                     // compute delta
                     dX(0) = -( dRdXinv(0,0)*R(0) + dRdXinv(0,1)*R(1) );
                     dX(1) = -( dRdXinv(1,0)*R(0) + dRdXinv(1,1)*R(1) );
 
                     // update
                     X += dX;
 
                     // if ( std::abs(X(0)) > 1.0 || std::abs(X(1)) > 1.0 )
                     //   {
                     //     std::cout << "WARNING: on iter = " << iter << ", xi is leaving the element!" << std::endl;
                     //   }
 
                     // recompute Residual
                     res(X,R);
 
                     // increment counter
                     iter++;
                   }
 
                 // check that we converged
                 if (iter==ITER_MAX)
                   {
                     // std::cerr << "Failed to converge distance function!!!" << std::endl;
                     // std::cerr << "Started with ||R|| = " << R0 << std::endl;
                     // std::cerr << "Finished " << iter << " iterations with ||R|| = " << R.l2_norm() << std::endl;
                     // std::cout << "Final location was xi = (" << X(0) << ", " << X(1) << ")." << std::endl;
                     // libmesh_error();
 
                     std::cout << "WARNING: Failed to converge distance function iteration.  Assuming that min is outside this face.  Use closest distance to edges to proceed." << std::endl;
                   }
                 else
                   {
                     //std::cout << "Converged!" << std::endl;
                   }
 
 
                 if ( std::abs(X(0)) > 1.0 || std::abs(X(1)) > 1.0 )
                   {
                     // converged to point outside the face... so min must
                     // be on edge, and luckily the edges are linear
 
                     libMesh::Real dtmp = std::numeric_limits<libMesh::Real>::infinity();
 
                     // for each edge of boundary face, find distance from
                     // input node to the segment defined by that edge
                     //
                     // the minimum of these distances minimizes the distance
                     // to the node over the points in the boundary face
                     //
                     for ( unsigned int iedge=0; iedge<4; iedge++ ){
 
                       Line line;
 
                       // Get the endpoints of this edge
                       if      (iedge==0) { line.p0 = belem->point(0); line.p1 = belem->point(1); }
                       else if (iedge==1) { line.p0 = belem->point(1); line.p1 = belem->point(2); }
                       else if (iedge==2) { line.p0 = belem->point(2); line.p1 = belem->point(3); }
                       else   /*iedge==3*/{ line.p0 = belem->point(3); line.p1 = belem->point(0); }
 
                       DistanceToSegment<3>(*node, line, dtmp);
 
                       if( dtmp < dedge ) dedge = dtmp;
                     }
                   }
                 else
                   {
                     // converged to point inside, so compute point in
                     // physical space and use it to evaluate the distance
 
                     // assuming first order lagrange basis here
                     libMesh::AutoPtr<libMesh::FEBase> fe( libMesh::FEBase::build(2, libMesh::FEType(libMesh::FIRST, libMesh::LAGRANGE)) );
 
                     std::vector<libMesh::Point> xi(1);
                     xi[0](0) = X(0);
                     xi[0](1) = X(1);
 
                     // grab basis functions (evaluated at qpts)
                     const std::vector<std::vector<libMesh::Real> > &basis = fe->get_phi();
 
                     // reinitialize finite element data at xi
                     fe->reinit(belem, &xi);
 
                     // interpolate location
                     libMesh::DenseVector<libMesh::Real> xx(3);
                     xx.zero();
 
                     for (unsigned int inode=0; inode<belem->n_nodes(); ++inode)
                       for (unsigned int idim=0; idim<3; ++idim)
                         xx(idim) += belem->point(inode)(idim) * basis[inode][0];
 
                     // compute distance
                     dedge = 0.0;
                     for ( unsigned int ii=0; ii<3; ii++ ) {
                       dedge += (xnode[ii] - xx(ii))*(xnode[ii] - xx(ii));
                     }
                     dedge = sqrt(dedge);
 
                   }
 
               }
             else
               { // higher-order faces not supported yet
                 std::cout << "My type is " << belem->type() << std::endl;
                 libmesh_not_implemented();
               }
 
             if ( dedge < distance ) distance = dedge;
 
           } // end if(2d)
 
       } // end loop over elements
 
     // make sure we are returning valid distance---i.e., 0 <= distance < infinity
     // But this may fail if we try to get the distance on a mesh with no
     // boundary
     //  libmesh_assert( (std::isfinite(distance)) && (distance>=0.0) );
 
     return distance;
   }

Member Data Documentation

const libMesh::UnstructuredMesh& GRINS::DistanceFunction::_boundary_mesh

private

Pointer to boundary mesh object

Definition at line 100 of file distance_function.h.

Referenced by compute(), and node_to_boundary().

libMesh::AutoPtr<libMesh::FEBase> GRINS::DistanceFunction::_dist_fe

private

Finite element to use for interpolation of distance. For internal use only, so don't provide any access to the outside world.

Currently type is hardcoded to first order, Lagrange (see constructor), but this could be easily changed.

Definition at line 110 of file distance_function.h.

Referenced by interpolate().

libMesh::EquationSystems& GRINS::DistanceFunction::_equation_systems

private

Pointer to EquationSystems object

Definition at line 95 of file distance_function.h.

Referenced by compute(), DistanceFunction(), interpolate(), and node_to_boundary().

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

src/utilities/include/grins/distance_function.h
src/utilities/src/distance_function.C

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation