GRINS-0.8.0
List of all members | Public Member Functions | Private Attributes
GRINS::DistanceFunction Class Reference

#include <distance_function.h>

Inheritance diagram for GRINS::DistanceFunction:
Inheritance graph
[legend]
Collaboration diagram for GRINS::DistanceFunction:
Collaboration graph
[legend]

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::UniquePtr< 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::UniquePtr< 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.

148  :
149  _equation_systems (es_in),
150  _boundary_mesh (bm_in),
151  _dist_fe (libMesh::FEBase::build(_equation_systems.get_mesh().mesh_dimension(), libMesh::FEType(libMesh::FIRST, libMesh::LAGRANGE)))
152  {
153  // Ensure that libmesh is ready to roll
154  libmesh_assert(libMesh::initialized());
155 
156  // Add distance function system
157  _equation_systems.add_system<libMesh::System>("distance_function");
158 
159  // Get reference to distance function system we just added
160  libMesh::System& sys = _equation_systems.get_system<libMesh::System>("distance_function");
161 
162  // Add distance function variable
163  sys.add_variable("distance", libMesh::FIRST);
164 
165  // Attach initialization function
166  sys.attach_init_object(*this);
167 
168  }
libMesh::UniquePtr< libMesh::FEBase > _dist_fe
libMesh::EquationSystems & _equation_systems
const libMesh::UnstructuredMesh & _boundary_mesh
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().

558  {
559  // Get mesh
560  const libMesh::MeshBase& mesh = _equation_systems.get_mesh();
561 
562  // Get reference to system and system number
563  libMesh::System& system = _equation_systems.get_system<libMesh::System>("distance_function");
564  const unsigned int sys_num = system.number();
565 
566  //Debugging code
567  //system.print_info();
568  //_equation_systems.print_info();
569  //End debugging
570 
571  // The boundary mesh needs to all be on this processor for us to
572  // calculate a correct distance function. Since we don't need it to
573  // be serial afterwards, we use a temporary serializer.
574  {
575  libMesh::MeshSerializer serialize(const_cast<libMesh::UnstructuredMesh&>(_boundary_mesh));
576 
577  // Loop over nodes in mesh
578  libMesh::MeshBase::const_node_iterator node_it = mesh.local_nodes_begin();
579  const libMesh::MeshBase::const_node_iterator node_end = mesh.local_nodes_end();
580 
581  //Debugging code
582  //const libMesh::Elem *elem = *mesh.elements_begin();
583  //if (elem->default_order() == libMesh::FIRST)
584  // std::cout<<"Mesh is using FIRST order elements."<<std::endl;
585  //End debugging
586 
587 
588  for ( ; node_it != node_end; ++node_it) {
589 
590  // Grab node
591  const libMesh::Node* node = *node_it;
592 
593  // Compute distance to nearest point in boundary_mesh
594  const libMesh::Real distance = DistanceFunction::node_to_boundary (node);
595 
596  // Stuff data into appropriate place in the system solution
597  const unsigned int dof = node->dof_number(sys_num,0,0);
598  system.solution->set (dof, distance);
599 
600  } // end loop over nodes
601 
602  } // end boundary mesh serialization
603 
604  system.solution->close();
605  system.update();
606 
607  std::cout << "Distance Function computed." << std::endl;
608  }
libMesh::EquationSystems & _equation_systems
libMesh::Real node_to_boundary(const libMesh::Node *node)
const libMesh::UnstructuredMesh & _boundary_mesh
void GRINS::DistanceFunction::initialize ( )
virtual

Initializes the distance

Definition at line 173 of file distance_function.C.

References compute().

174  {
175  // Call the compute function
176  this->compute();
177  return;
178  }
libMesh::UniquePtr< 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.

616  {
617  libmesh_assert( elem != NULL ); // can't interpolate in NULL elem
618  libmesh_assert( qpts.size() > 0 ); // can't interpolate if no points requested
619 
620  // grab basis functions (evaluated at qpts)
621  const std::vector<std::vector<libMesh::Real> > &phi = _dist_fe->get_phi();
622 
623  // reinitialize finite element data at qpts
624  _dist_fe->reinit(elem, &qpts);
625 
626  // number of basis functions
627  const unsigned int n_dofs = phi.size();
628 
629  // number of points
630  const unsigned int n_pts = qpts.size();
631 
632  // instantiate auto_ptr to dense vector to hold results
633  libMesh::UniquePtr< DenseVector<libMesh::Real> > ap( new libMesh::DenseVector<libMesh::Real>(qpts.size()) );
634  (*ap).zero();
635 
636  // pull off distance function at nodes on this element
637  libMesh::System& sys = _equation_systems.get_system<libMesh::System>("distance_function");
638  const libMesh::DofMap& dof_map = sys.get_dof_map();
639 
640  std::vector<unsigned int> dof_ind;
641  dof_map.dof_indices(elem, dof_ind);
642 
643  libMesh::DenseVector<libMesh::Real> nodal_dist;
644  nodal_dist.resize(n_dofs);
645  //dof_map.extract_local_vector( *(sys.solution), dof_ind, nodal_dist);
646  dof_map.extract_local_vector( *(sys.current_local_solution), dof_ind, nodal_dist);
647 
648  for ( unsigned int idof=0; idof<n_dofs; idof++ ) {
649  for ( unsigned int iqpt=0; iqpt<n_pts; iqpt++ ) {
650  (*ap)(iqpt) += nodal_dist(idof) * phi[idof][iqpt];
651  }
652  }
653 
654  return ap;
655  }
libMesh::UniquePtr< libMesh::FEBase > _dist_fe
libMesh::EquationSystems & _equation_systems
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().

185  {
186  // Ensure that node is not NULL
187  libmesh_assert( node != NULL );
188 
189  // Initialize distance to infinity
190  libMesh::Real distance = std::numeric_limits<libMesh::Real>::infinity();
191 
192  // Get dimension
193  const unsigned int dim = _equation_systems.get_mesh().mesh_dimension();
194  libmesh_assert( (dim==2) || (dim==3) );
195 
196  // Get coordinates of node
197  std::vector<libMesh::Real> xnode(dim);
198  for( unsigned int ii=0; ii<dim; ii++ ) xnode[ii] = (*node)(ii);
199 
200  // This function will work on a distributed interior mesh, but won't
201  // give correct results on a distributed boundary mesh.
202  libmesh_assert(_boundary_mesh.is_serial());
203 
204  // First loop over the boundary mesh and compute the distance to the
205  // nodes. This will give us an idea of what elements to check more
206  // closely.
207 
208  libMesh::MeshBase::const_element_iterator el = _boundary_mesh.active_elements_begin();
209  const libMesh::MeshBase::const_element_iterator end_el = _boundary_mesh.active_elements_end();
210 
211  if (el==end_el)
212  {
213  std::cout << "There are no elements to iterate through!!! Returning..." << std::endl;
214  return distance;
215  }
216 
217  const libMesh::Elem* elem_containing_min = NULL;
218  libMesh::Real dmin = std::numeric_limits<Real>::infinity();
219 
220  for ( ; el != end_el; ++el) {
221 
222  const libMesh::Elem* belem = *el;
223 
224  for (unsigned int bnode=0; bnode<belem->n_nodes(); ++bnode)
225  {
226  const libMesh::Point& pbndry = belem->point(bnode);
227 
228  libMesh::Real dnode = 0.0;
229  for (unsigned int idim=0; idim<dim; ++idim)
230  dnode += (xnode[idim] - pbndry(idim))*(xnode[idim] - pbndry(idim));
231 
232  dnode = std::sqrt(dnode);
233 
234  if (dnode<dmin)
235  {
236  dmin = dnode;
237  elem_containing_min = belem;
238  }
239  }
240 
241  } // end element loop
242 
243  // grab the elements around the minimum
244  std::set<const libMesh::Elem*> near_min_elems;
245  elem_containing_min->find_point_neighbors (near_min_elems);
246 
247  // insert the element itself into the set... we must check it also
248  near_min_elems.insert(elem_containing_min);
249 
250  std::set<const libMesh::Elem*>::iterator nm_el;
251 
252  // loop over the near_min_elems
253  for (nm_el=near_min_elems.begin(); nm_el!=near_min_elems.end(); ++nm_el)
254  {
255 
256  const libMesh::Elem* belem = *nm_el; //*el;
257 
258  // Ensure that elem defined by edge/face is linear
259  libmesh_assert( belem->default_order() == FIRST );
260  // For now we assume that the boundary_elems are linear, but give a warning
261  //libmesh_warning( "Moving ahead, but need to make sure that the boundary elems are linear");
262 
263  // Initialize distance to this edge/face to infinity
264  libMesh::Real dedge = std::numeric_limits<libMesh::Real>::infinity();
265 
266  if ( dim==2 )
267  { // 2d
268 
269  // For 2d mesh, boundary elems must be 1d. Due to the way
270  // find_point_neighbors works, there may be some 2d elements
271  // in the near_min_elems set. These elements are clearly
272  // not part of the boundary mesh, and thus we skip them.
273  //
274  if (belem->dim()!=1) continue;
275 
276  //std::cout<<"Number of nodes: "<<belem->n_nodes()<<std::endl;
277 
278  libmesh_assert( belem->n_nodes() == 2 );
279  libmesh_assert( belem->type() == EDGE2 );
280 
281  // Points defining the edge
282  const libMesh::Point& p0 = belem->point(0);
283  const libMesh::Point& p1 = belem->point(1);
284 
285  Line line;
286  line.p0 = p0;
287  line.p1 = p1;
288 
289  DistanceToSegment<2>(*node, line, dedge);
290 
291  if ( dedge < distance ) distance = dedge;
292 
293  }
294  else
295  { // 3d
296 
297  // For 3d mesh, boundary elems must be 2d. Due to the way
298  // find_point_neighbors works, there may be some 3d elements
299  // in the near_min_elems set. These elements are clearly
300  // not part of the boundary mesh, and thus we skip them.
301  //
302  if (belem->dim()!=2) continue;
303 
304  libmesh_assert( (belem->type() == TRI3) || (belem->type() == QUAD4) );
305 
306  if ( belem->type() == TRI3 )
307  { // triangular boundary faces
308 
309  // Find the point in the plane defined by the boundary nodes that
310  // minimizes the distance between the plane and the input node.
311  //
312  // Done in terms of the reference coordinates of the boundary element.
313  //
314  const libMesh::Point& p0 = belem->point(0);
315  const libMesh::Point& p1 = belem->point(1);
316  const libMesh::Point& p2 = belem->point(2);
317 
318  const libMesh::Real x_xi = (p1(0) - p0(0)), x_et = (p2(0) - p0(0));
319  const libMesh::Real y_xi = (p1(1) - p0(1)), y_et = (p2(1) - p0(1));
320  const libMesh::Real z_xi = (p1(2) - p0(2)), z_et = (p2(2) - p0(2));
321 
322  libMesh::Real A[2][2], b[2];
323 
324  A[0][0] = x_xi*x_xi + y_xi*y_xi + z_xi*z_xi;
325  A[0][1] = x_xi*x_et + y_xi*y_et + z_xi*z_et;
326  A[1][0] = x_xi*x_et + y_xi*y_et + z_xi*z_et;
327  A[1][1] = x_et*x_et + y_et*y_et + z_et*z_et;
328 
329  b[0] = (xnode[0] - p0(0))*x_xi + (xnode[1] - p0(1))*y_xi + (xnode[2] - p0(2))*z_xi;
330  b[1] = (xnode[0] - p0(0))*x_et + (xnode[1] - p0(1))*y_et + (xnode[2] - p0(2))*z_et;
331 
332  const libMesh::Real detA = A[0][0]*A[1][1] - A[1][0]*A[0][1];
333  libmesh_assert( fabs(detA) > 0.0 ); // assert that A is not singular
334 
335  const libMesh::Real xi = ( A[1][1]*b[0] - A[0][1]*b[1])/detA;
336  const libMesh::Real et = (-A[1][0]*b[0] + A[0][0]*b[1])/detA;
337 
338 
339  // If projection of node onto plane defined by boundary face
340  // (i.e., what we just computed) is not inside boundary face
341  // then we need to figure out closest point that is inside the
342  // boundary face
343  //
344  if ( (xi<0.0) || (et<0.0) || (et>1.0-xi) ) {
345 
346  libMesh::Real dtmp = std::numeric_limits<libMesh::Real>::infinity();
347 
348  // for each edge of boundary face, find distance from
349  // input node to the segment defined by that edge
350  //
351  // the minimum of these distances minimizes the distance
352  // to the node over the points in the boundary face
353  //
354  for ( unsigned int iedge=0; iedge<3; iedge++ ){
355 
356  Line line;
357 
358  // Get the endpoints of this edge
359  if (iedge==0) { line.p0 = p1; line.p1 = p2; }
360  else if (iedge==1) { line.p0 = p0; line.p1 = p2; }
361  else /*iedge==2*/{ line.p0 = p0; line.p1 = p1; }
362 
363  DistanceToSegment<3>(*node, line, dtmp);
364 
365  if( dtmp < dedge ) dedge = dtmp;
366 
367  }
368 
369  } else { // projection is inside face, so we're good to go
370 
371  // Map from reference to physical space
372  libMesh::Real xint[3];
373  for ( unsigned int ii=0; ii<3; ii++ ) {
374  xint[ii] = p0(ii)*(1.0 - xi - et) + p1(ii)*xi + p2(ii)*et;
375  }
376 
377  // compute distance
378  dedge = 0.0;
379  for ( unsigned int ii=0; ii<3; ii++ ) {
380  dedge += (xnode[ii] - xint[ii])*(xnode[ii] - xint[ii]);
381  }
382  dedge = sqrt(dedge);
383 
384  }
385 
386  }
387  else if ( belem->type() == QUAD4 )
388  {
389  //std::cout << "WARNING: this functionality is not well-tested. Sorry." << std::endl;
390 
391  libMesh::Real RTOL = 1e-10;
392  libMesh::Real ATOL = 1e-20;
393  const unsigned int ITER_MAX=1000;
394  unsigned int iter=0;
395 
396  ComputeDistanceResidual res(belem, node);
397 
398  ComputeDistanceJacobian jac;
399 
400  libMesh::DenseVector<libMesh::Real> X(2); X(0) = X(1) = 0.0;
401  libMesh::DenseVector<libMesh::Real> dX(2);
402  libMesh::Real det;
403 
404  libMesh::DenseVector<libMesh::Real> R(2);
405  libMesh::DenseMatrix<libMesh::Real> dRdX(2,2);
406  libMesh::DenseMatrix<libMesh::Real> dRdXinv(2,2);
407 
408  // evaluate residual... maybe we don't have to iterate
409  res(X,R);
410 
411  libMesh::Real R0 = R.l2_norm();
412 
413  // iterate
414  while ( (R.l2_norm() > ATOL) && (R.l2_norm()/R0 > RTOL) && (iter<ITER_MAX) )
415  {
416  // compute jacobian
417  jac(X, R, res, dRdX);
418 
419  // invert jacobian
420  det = dRdX(0,0)*dRdX(1,1) - dRdX(0,1)*dRdX(1,0);
421 
422  // protect against divide by zero
423  if(std::abs(det)<=1e-20)
424  {
425  std::cout << "WARNING: about to divide by zero." << std::endl;
426  }
427 
428  dRdXinv(0,0) = dRdX(1,1)/det;
429  dRdXinv(0,1) = -dRdX(0,1)/det;
430  dRdXinv(1,0) = -dRdX(1,0)/det;
431  dRdXinv(1,1) = dRdX(0,0)/det;
432 
433  // compute delta
434  dX(0) = -( dRdXinv(0,0)*R(0) + dRdXinv(0,1)*R(1) );
435  dX(1) = -( dRdXinv(1,0)*R(0) + dRdXinv(1,1)*R(1) );
436 
437  // update
438  X += dX;
439 
440  // if ( std::abs(X(0)) > 1.0 || std::abs(X(1)) > 1.0 )
441  // {
442  // std::cout << "WARNING: on iter = " << iter << ", xi is leaving the element!" << std::endl;
443  // }
444 
445  // recompute Residual
446  res(X,R);
447 
448  // increment counter
449  iter++;
450  }
451 
452  // check that we converged
453  if (iter==ITER_MAX)
454  {
455  // std::cerr << "Failed to converge distance function!!!" << std::endl;
456  // std::cerr << "Started with ||R|| = " << R0 << std::endl;
457  // std::cerr << "Finished " << iter << " iterations with ||R|| = " << R.l2_norm() << std::endl;
458  // std::cout << "Final location was xi = (" << X(0) << ", " << X(1) << ")." << std::endl;
459  // libmesh_error();
460 
461  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;
462  }
463  else
464  {
465  //std::cout << "Converged!" << std::endl;
466  }
467 
468 
469  if ( std::abs(X(0)) > 1.0 || std::abs(X(1)) > 1.0 )
470  {
471  // converged to point outside the face... so min must
472  // be on edge, and luckily the edges are linear
473 
474  libMesh::Real dtmp = std::numeric_limits<libMesh::Real>::infinity();
475 
476  // for each edge of boundary face, find distance from
477  // input node to the segment defined by that edge
478  //
479  // the minimum of these distances minimizes the distance
480  // to the node over the points in the boundary face
481  //
482  for ( unsigned int iedge=0; iedge<4; iedge++ ){
483 
484  Line line;
485 
486  // Get the endpoints of this edge
487  if (iedge==0) { line.p0 = belem->point(0); line.p1 = belem->point(1); }
488  else if (iedge==1) { line.p0 = belem->point(1); line.p1 = belem->point(2); }
489  else if (iedge==2) { line.p0 = belem->point(2); line.p1 = belem->point(3); }
490  else /*iedge==3*/{ line.p0 = belem->point(3); line.p1 = belem->point(0); }
491 
492  DistanceToSegment<3>(*node, line, dtmp);
493 
494  if( dtmp < dedge ) dedge = dtmp;
495  }
496  }
497  else
498  {
499  // converged to point inside, so compute point in
500  // physical space and use it to evaluate the distance
501 
502  // assuming first order lagrange basis here
503  libMesh::UniquePtr<libMesh::FEBase> fe( libMesh::FEBase::build(2, libMesh::FEType(libMesh::FIRST, libMesh::LAGRANGE)) );
504 
505  std::vector<libMesh::Point> xi(1);
506  xi[0](0) = X(0);
507  xi[0](1) = X(1);
508 
509  // grab basis functions (evaluated at qpts)
510  const std::vector<std::vector<libMesh::Real> > &basis = fe->get_phi();
511 
512  // reinitialize finite element data at xi
513  fe->reinit(belem, &xi);
514 
515  // interpolate location
516  libMesh::DenseVector<libMesh::Real> xx(3);
517  xx.zero();
518 
519  for (unsigned int inode=0; inode<belem->n_nodes(); ++inode)
520  for (unsigned int idim=0; idim<3; ++idim)
521  xx(idim) += belem->point(inode)(idim) * basis[inode][0];
522 
523  // compute distance
524  dedge = 0.0;
525  for ( unsigned int ii=0; ii<3; ii++ ) {
526  dedge += (xnode[ii] - xx(ii))*(xnode[ii] - xx(ii));
527  }
528  dedge = sqrt(dedge);
529 
530  }
531 
532  }
533  else
534  { // higher-order faces not supported yet
535  std::cout << "My type is " << belem->type() << std::endl;
536  libmesh_not_implemented();
537  }
538 
539  if ( dedge < distance ) distance = dedge;
540 
541  } // end if(2d)
542 
543  } // end loop over elements
544 
545  // make sure we are returning valid distance---i.e., 0 <= distance < infinity
546  // But this may fail if we try to get the distance on a mesh with no
547  // boundary
548  // libmesh_assert( (std::isfinite(distance)) && (distance>=0.0) );
549 
550  return distance;
551  }
libMesh::EquationSystems & _equation_systems
const libMesh::UnstructuredMesh & _boundary_mesh

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::UniquePtr<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:

Generated on Tue Dec 19 2017 12:47:30 for GRINS-0.8.0 by  doxygen 1.8.9.1