This class provides some common primitives used by various propagation algorithms. More...

#include <Propagation_3.h>

Inheritance diagram for Propagation_3:

Collaboration diagram for Propagation_3:

Public Member Functions
	Propagation_3 ()

virtual	~Propagation_3 ()

	Propagation_3 (Manifold wm, COM::Window buf)
	Construct an object from a window manifold. More...

virtual double	time_stepping (const COM::Attribute spd, double dt, COM::Attribute disps, int *smoothed=NULL)=0
	Main entry for invoking a surface propagation algorithm. More...

virtual void	set_constraints (const COM::Attribute *cnstr_types)
	Set the types and directions of nodal constraints. More...

void	set_bounds (const COM::Attribute *bnd)
	Set the bounds. More...

void	bound_facial_speed (COM::Attribute *fa)

virtual void	enforce_nodal_constraints (COM::Attribute *du)
	Enforces the nodal constraints by projecting motion onto given direction. More...

virtual void	bound_nodal_motion (COM::Attribute *disps)

void	set_verbose (bool b)
	Set the verbose level. More...

Protected Member Functions
void	convert_constraints (const COM::Attribute ctypes_faces, COM::Attribute ctypes_nodes)
	Convert facial constraints to nodal constraints. More...

void	determine_constraint_boundary (const COM::Attribute ctypes_faces, COM::Attribute ctypes_bndry_edges, COM::Attribute *ctypes_bndry_nodes)
	Convert facial or panel constraints to nodal constraints. More...

Static Protected Member Functions
static void	get_constraint_directions (int type, int &ndirs, Vector_3 dirs[2])
	Get orthonormals of the constraints. More...

static void	enforce_nodal_constraint (int type, Vector_3 &du)
	Enforce constraint for a specific vector. More...

static void	bound_nodal_motion (const Point_3 &pnt, const double *bnd, Vector_3 &du, double eps=0)

static bool	check_spherical_bound (const Point_3 &pnt, const Point_3 &org, const double rad_min, const double rad_max, double eps=0)

static void	bound_spherical_disp (const Point_3 &pnt, const Point_3 &org, const double rad_min, const double rad_max, Vector_3 &disp, double eps=0)

static bool	check_radial_bound (const double x, const double y, const double bnd_min, const double bnd_max, double eps=0)

static void	bound_radial_disp (const double x, const double y, const double bnd_min, const double bnd_max, double &dx, double &dy, double eps=0)

static bool	check_axial_bound (const double x, const double bnd_min, const double bnd_max, double eps=0)

static void	bound_axial_disp (const double x, const double bnd_min, const double bnd_max, double &dx, double eps=0)

static bool	reached_nodal_bound (const Point_3 &pnt, const double *bnd, double eps=0)

static bool	in_bounding_box (const Point_3 &pnt, const Point_3 &lb, const Point_3 &ub)

static double	square (double x)

Protected Attributes
Manifold *	_surf

SURF::Access_Mode	_mode

COM::Window *	_buf

int	_rank

bool	_verb

bool	_cnstr_set

bool	_bnd_set

COM::Attribute *	_cnstr_nodes

COM::Attribute *	_cnstr_faces

COM::Attribute *	_cnstr_bndry_edges

COM::Attribute *	_cnstr_bndry_nodes

COM::Attribute *	_cnstr_bound

std::vector< COM::Pane * >	_panes

Detailed Description

This class provides some common primitives used by various propagation algorithms.

It is a pure virtual class.

Definition at line 40 of file Propagation_3.h.

Constructor & Destructor Documentation

Propagation_3 ( )

inline

Definition at line 42 of file Propagation_3.h.

     : _surf(NULL), _mode( SURF::ACROSS_PANE), _buf(NULL), _rank(0), 
       _verb(0), _cnstr_set(0), _cnstr_nodes(0), _cnstr_faces(0), 
       _cnstr_bndry_edges(0), _cnstr_bndry_nodes(), _cnstr_bound(0) {}

virtual ~Propagation_3 ( )

inlinevirtual

Definition at line 47 of file Propagation_3.h.

47 {}

PROP_BEGIN_NAMESPACE Propagation_3	(	Manifold *	wm,
		COM::Window *	buf
	)

Construct an object from a window manifold.

Definition at line 31 of file Propagation_3.C.

References _buf, _cnstr_bndry_edges, _cnstr_bndry_nodes, _cnstr_bound, _cnstr_faces, _cnstr_nodes, _mode, _panes, _rank, ACROSS_PANE, BOUND_LEN, COM_CHAR, COM_DOUBLE, COM_INT, COMMPI_Comm_rank(), COMMPI_Initialized(), and WHOLE_PANE.

   : _surf(wm), _buf(buf), _cnstr_set(false), _bnd_set(0), 
     _cnstr_nodes(NULL), _cnstr_bound(NULL) {
 
   // Depending on whether we have pconn for ghost nodes and elements,
   // run in either WHOLE_PANE or ACROSS_PANE mode.
   if ( wm->pconn_nblocks() >=3)
     _mode = SURF::WHOLE_PANE;
   else
     _mode = SURF::ACROSS_PANE;
 
   if ( _buf) {
     _cnstr_nodes = _buf->new_attribute( "PROP_cnstr_nodes", 'n', COM_INT, 1, "");
     _cnstr_faces = _buf->new_attribute( "PROP_cnstr_faces", 'e', COM_INT, 1, "");
     _cnstr_bndry_edges = _buf->new_attribute( "PROP_cnstr_bndry_edges", 'e', COM_CHAR, 1, "");
     _cnstr_bndry_nodes = _buf->new_attribute( "PROP_cnstr_bndry_nodes", 'n', COM_CHAR, 1, "");
     _cnstr_bound = _buf->new_attribute( "PROP_cnstr_bound", 'p', COM_DOUBLE, BOUND_LEN, "");
 
     _buf->panes( _panes); 
   }
 
   MPI_Comm comm = _buf->get_communicator();
   if ( COMMPI_Initialized())  _rank = COMMPI_Comm_rank(comm);
   else _rank = 0;
 }

Here is the call graph for this function:

Member Function Documentation

void bound_axial_disp	(	const double	x,
		const double	bnd_min,
		const double	bnd_max,
		double &	dx,
		double	eps = `0`
	)

staticprotected

Definition at line 377 of file Propagation_3.C.

References dx, max(), min(), and x.

Referenced by bound_nodal_motion().

                                                   {
 
   double xnew = x+dx;
   if ( std::max( x, xnew) >= bnd_max-eps && std::min(x, xnew) <= bnd_max+eps) {
     dx = bnd_max - x;
   }
   else if ( std::min( x, xnew) <= bnd_min+eps && std::max(x, xnew) >= bnd_min-eps) {
     dx = bnd_min - x;
   }
   else if ( x >= bnd_max+eps && xnew > x || x<=bnd_min-eps && xnew < x)
     dx = 0; // Do not alter dx if moves inwards
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void bound_facial_speed ( COM::Attribute * fa )

Definition at line 550 of file Propagation_3.C.

References _bnd_set, _buf, _cnstr_bound, _panes, BOUND_LEN, COM_assertion, COM_assertion_msg, COM_compatible_types(), COM_DOUBLE, COM_NC, i, iend, j, k, min(), Element_node_enumerator::next(), nj, reached_nodal_bound(), and Element_node_enumerator::size_of_edges().

                                       {
   COM_assertion( fa->window() == _buf); 
 
   COM_assertion_msg( _bnd_set, "Bound must be set first before calling bound_facial_speed");
 
   COM_assertion_msg( fa->is_elemental() && fa->size_of_components()==1 &&
                      COM_compatible_types( fa->data_type(), COM_DOUBLE),
                      "Propagation_3::bound_facial_speed expects a scalar double-precision elemental attribute.");
 
   std::vector< COM::Pane*>::iterator it = _panes.begin(), iend=_panes.end();
   // Loop through the panes and its real nodes
   for ( it=_panes.begin(); it!=iend; ++it) {
     COM::Pane *pane = *it;
     
     double *val_dbl = (double*)pane->attribute(fa->id())->pointer();
     // Skip the panes without the attribute.
     if ( val_dbl==NULL) continue;
 
     COM_assertion( pane->attribute(COM_NC)->stride()==3);
     const Point_3 *pnts = reinterpret_cast<const Point_3*>(pane->coordinates()); 
     int nbnd = _bnd_set ? pane->attribute( _cnstr_bound->id())->size_of_items() : 0;
     if ( nbnd==0) continue;
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1);
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       // Skip the faces that are already initialized to 0.
       if ( val_dbl[j]==0) continue;
       int ne = ene.size_of_edges();
 
       const double *bnd = reinterpret_cast<const double*>
         ( pane->attribute( _cnstr_bound->id())->pointer());
       bool out[] = {false, false, false, false};
       // Bound all nodes
       for ( int i=0; i<nbnd; ++i, bnd+=BOUND_LEN) {
         double eps = 1.e-5 * std::min( std::min( bnd[9]-bnd[6], bnd[10]-bnd[7]),
                                        bnd[11]-bnd[8]);
         if ( eps>1 || eps <=0 ) eps = 1.e-10;
 
         for ( int k=0; k<ne; ++k) {
           out[k] |= reached_nodal_bound( pnts[ene[k]-1], bnd, eps);
         }
       }
       
       if ( out[0] && out[1] && out[2] && (ne==3 || out[3])) val_dbl[j] = 0;
     }
   }
 }

Here is the call graph for this function:

void bound_nodal_motion ( COM::Attribute * disps )

virtual

Definition at line 518 of file Propagation_3.C.

References _bnd_set, _cnstr_bound, _panes, BOUND_LEN, COM_NC, i, iend, min(), and n.

Referenced by FaceOffset_3::time_stepping().

                                       {
   if ( !_bnd_set) return;
 
   std::vector< COM::Pane*>::iterator it = _panes.begin(), iend=_panes.end();
   // Loop through the panes and its real nodes
   for ( it=_panes.begin(); it!=iend; ++it) {
     COM::Pane *pane = *it;
     
     Vector_3 *ds = reinterpret_cast<Vector_3*>
       (pane->attribute( du->id())->pointer());
 
     const Point_3 *pnts = reinterpret_cast<const Point_3*>
       (pane->attribute(COM_NC)->pointer());
     const double *bnd = _bnd_set ? reinterpret_cast<const double*>
       ( pane->attribute( _cnstr_bound->id())->pointer()) : NULL;
     int nbnd = _bnd_set ? pane->attribute( _cnstr_bound->id())->size_of_items() : 0;
 
     // Bound all nodes
     for ( int i=0; i<nbnd; ++i, bnd+=BOUND_LEN) {
       double eps = 1.e-5 * std::min( std::min( bnd[9]-bnd[6], bnd[10]-bnd[7]),
                                      bnd[11]-bnd[8]);
       if ( eps>1 || eps <= 0) eps = 1.e-10;
 
       for (int i=0, n=pane->size_of_real_nodes(); i<n; ++i)
         bound_nodal_motion( pnts[i], bnd, ds[i]);
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void bound_nodal_motion	(	const Point_3 &	pnt,
		const double *	bnd,
		Vector_3 &	du,
		double	eps = `0`
	)

staticprotected

Definition at line 393 of file Propagation_3.C.

References bound_axial_disp(), bound_radial_disp(), bound_spherical_disp(), COM_assertion_msg, and in_bounding_box().

                                               {
   if ( !in_bounding_box( pnt, (const Point_3&)bnd[6], (const Point_3&)bnd[9]))
     return;
 
   switch ( (int)bnd[0]) {
   case 's':  bound_spherical_disp( pnt, (const Point_3&)bnd[3],
                                    bnd[1], bnd[2], du, eps); break;
 
   case 'x':  bound_axial_disp( pnt[0]-bnd[3], bnd[1], bnd[2], du[0], eps); break;
   case 'y':  bound_axial_disp( pnt[1]-bnd[4], bnd[1], bnd[2], du[1], eps); break;
   case 'z':  bound_axial_disp( pnt[2]-bnd[5], bnd[1], bnd[2], du[2], eps); break;
 
   case -'x':  bound_radial_disp( pnt[1]-bnd[4], pnt[2]-bnd[5], bnd[1], bnd[2], 
                                  du[1], du[2], eps); break;
   case -'y':  bound_radial_disp( pnt[0]-bnd[3], pnt[2]-bnd[5], bnd[1], bnd[2], 
                                  du[0], du[2], eps); break;
   case -'z':  bound_radial_disp( pnt[0]-bnd[3], pnt[1]-bnd[4], bnd[1], bnd[2], 
                                  du[0], du[1], eps); break;
     
   default: COM_assertion_msg( false, "Unknown type of bnd");
   }
 }

Here is the call graph for this function:

void bound_radial_disp	(	const double	x,
		const double	y,
		const double	bnd_min,
		const double	bnd_max,
		double &	dx,
		double &	dy,
		double	eps = `0`
	)

staticprotected

Definition at line 343 of file Propagation_3.C.

References max(), min(), sqrt(), square(), x, and y.

Referenced by bound_nodal_motion().

                                                                            {
 
   double sq_xy = square( x)+square( y);
   double sq_rad_actual=square( x+dx)+square( y+dy);
   double t;
 
   if ( std::max( sq_xy, sq_rad_actual) >= square((t=bnd_max)-eps) && 
        std::min( sq_xy, sq_rad_actual) <= square(bnd_max+eps) ||
        std::min( sq_xy, sq_rad_actual) <= square((t=bnd_min)+eps) && 
        std::max( sq_xy, sq_rad_actual) >= square(bnd_min-eps)) {
 
     if ( sq_xy == 0)
       dx = dy = 0;
     else {
       // Try to propagate toward the cylinder
       double ratio = t / std::sqrt( sq_xy);
       dx = ratio*x - x;
       dy = ratio*y - y;
     }
   }
   else if ( sq_xy >= square(bnd_max-eps) && sq_rad_actual >= sq_xy || 
             bnd_min>=0 && sq_xy <= square(bnd_min+eps) && sq_rad_actual<=sq_xy) {
     dx = dy = 0; // Do not alter dx if moves inwards
   }
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void bound_spherical_disp	(	const Point_3 &	pnt,
		const Point_3 &	org,
		const double	rad_min,
		const double	rad_max,
		Vector_3 &	disp,
		double	eps = `0`
	)

staticprotected

Definition at line 305 of file Propagation_3.C.

References Vector_3< Type >::is_null(), max(), min(), square(), and Vector_3< Type >::squared_norm().

Referenced by bound_nodal_motion().

                                                   {
   Vector_3 dir = pnt - org;
   double sqnrm = dir.squared_norm(), sqnrm1 = (dir+disp).squared_norm(), t;
 
   if ( std::max( sqnrm, sqnrm1) >= square((t=rad_max)-eps) && 
        std::min( sqnrm, sqnrm1) <= square(rad_max+eps) ||
        (t=rad_min) >= 0 && std::min( sqnrm, sqnrm1) <= square(rad_min+eps) && 
        std::max( sqnrm, sqnrm1) >= square(rad_min-eps)) {
 
     if ( dir.is_null()) 
       disp[0] = disp[1] = disp[2] = 0;
     else {
       Vector_3 dir_unit; (dir_unit = dir).normalize();
       disp += ( t - disp*dir_unit) * dir_unit - dir;
     }
   }
   else if ( sqnrm >= square(rad_max) && sqnrm1 > sqnrm || 
             rad_min >= 0 && sqnrm <= square(rad_min) && sqnrm1 < sqnrm) {
     if ( dir.is_null()) 
       disp[0] = disp[1] = disp[2] = 0;
     else
       // Only move tangentially
       disp -= disp * dir / dir.squared_norm();
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

bool check_axial_bound	(	const double	x,
		const double	bnd_min,
		const double	bnd_max,
		double	eps = `0`
	)

staticprotected

Definition at line 372 of file Propagation_3.C.

Referenced by reached_nodal_bound().

                                                                          {
   return x >= bnd_max-eps || x <= bnd_min+eps;
 }

Here is the caller graph for this function:

bool check_radial_bound	(	const double	x,
		const double	y,
		const double	bnd_min,
		const double	bnd_max,
		double	eps = `0`
	)

staticprotected

Definition at line 333 of file Propagation_3.C.

References square().

Referenced by reached_nodal_bound().

                                                     {
 
   double sq_xy = square( x)+square( y);
   return  sq_xy >= square(bnd_max-eps) || 
     bnd_min >=0 && sq_xy <= square(bnd_min+eps);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

bool check_spherical_bound	(	const Point_3 &	pnt,
		const Point_3 &	org,
		const double	rad_min,
		const double	rad_max,
		double	eps = `0`
	)

staticprotected

Definition at line 295 of file Propagation_3.C.

References square().

Referenced by reached_nodal_bound().

                                    {
   double sqnrm = (pnt - org).squared_norm();
 
   return sqnrm >= square( rad_max-eps) || 
     rad_min >=0 && sqnrm <= square( rad_min+eps);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void convert_constraints	(	const COM::Attribute *	ctypes_faces,
		COM::Attribute *	ctypes_nodes
	)

protected

Convert facial constraints to nodal constraints.

Definition at line 118 of file Propagation_3.C.

References _panes, _surf, NTS::abs(), COM_assertion, COM_assertion_msg, COM_NC, Rocblas::copy_scalar(), Vector_3< Type >::cross_product(), iend, j, k, matchall(), matchany(), Element_node_enumerator::next(), nj, Vector_3< Type >::normalize(), Element_node_enumerator::size_of_edges(), and Element_node_enumerator::size_of_nodes().

Referenced by set_constraints().

                                                  {
   int zero = 0;
   Rocblas::copy_scalar( &zero, ctypes_nodes);
 
   // Loop through all the nodes to compute their displacements
   enum { FIXED=1, XDIR = 2, YDIR = 4, ZDIR = 8, 
          NOX = 16, NOY = 32, NOZ = 64};
 
   std::vector< COM::Pane*>::iterator it = _panes.begin(), iend=_panes.end();
   // Loop through the panes and its real nodes
   for ( it=_panes.begin(); it!=iend; ++it) {
     COM::Pane *pane = *it;
 
     COM_assertion( pane->attribute(COM_NC)->stride()==3);
     const Point_3 *pnts = 
       reinterpret_cast<const Point_3*>(pane->coordinates()); 
     
     const int *val_faces=(const int*)pane->attribute(ctypes_faces->id())->pointer();
     int *val_nodes=(int*)pane->attribute(ctypes_nodes->id())->pointer();
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1); 
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int type = val_faces[j];
 
       int ne = ene.size_of_edges();
       // Loop through all vertices of current face.
       for ( int k=0; k<ne; ++k) {
         int vindex = ene[k]-1;
 
         switch ( type) {
         case 0: break;
         case 2:    val_nodes[vindex] |= FIXED; break;
         case 'x':  val_nodes[vindex] |= XDIR; break;
         case -'x': val_nodes[vindex] |= NOX;   break;
         case 'y':  val_nodes[vindex] |= YDIR; break;
         case -'y': val_nodes[vindex] |= NOY;   break;
         case 'z':  val_nodes[vindex] |= ZDIR; break;
         case -'z': val_nodes[vindex] |= NOZ;   break;
         case 't':
         case -'t': {
           SURF::Generic_element_2 e(ene.size_of_edges(), ene.size_of_nodes());
           COM::Element_node_vectors_k_const<Point_3 > ps;
           ps.set( pnts, ene, 1);
 
           Vector_2 nc(0.5,0.5);
           Vector_3 J[2], nrm;
           e.Jacobian( ps, nc, J);
           nrm = Vector_3::cross_product( J[0], J[1]).normalize();
           
           double eps = 1.e-1;
           if ( std::abs(std::abs(nrm[0])-1) < eps) {
             if ( type > 0) val_nodes[vindex] |= NOX;
             else val_nodes[vindex] |= XDIR;
           }
           else if ( std::abs(std::abs(nrm[1])-1) < eps) {
             if ( type > 0) val_nodes[vindex] |= NOY;
             else val_nodes[vindex] |= YDIR;
           }
           else if ( std::abs(std::abs(nrm[2])-1) < eps) {
             if ( type > 0) val_nodes[vindex] |= NOZ;
             else val_nodes[vindex] |= ZDIR;
           }
           else
             COM_assertion_msg(false, "Plane must be along axial direction");
           break;
         }
         default: 
           COM_assertion_msg( false, "Unknown/unsupported type of constraint");
         } // switch
       } // for k
     } // for j
   } // for it
 
   _surf->reduce_on_shared_nodes( ctypes_nodes, Manifold::OP_BOR);
 
   // Loop through the panes and its real nodes
   for ( it=_panes.begin(); it!=iend; ++it) {
     COM::Pane *pane = *it;
     int *val_nodes = (int*)pane->attribute( ctypes_nodes->id())->pointer();
 
     for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
       int &type = val_nodes[j];
 
       if ( type == 0) continue; // no constrained
 
       if ( (type & FIXED) || matchall(type, NOX|NOY|NOZ) ||
            matchall(type, XDIR|YDIR) || matchall(type, XDIR|ZDIR) ||
            matchall(type, YDIR|ZDIR) )
         type = 2;
       else if ( matchall(type, NOX|NOY)) {
         if ( matchany(type, XDIR|YDIR)) { type = 2; }
         else { type = 'z'; }
       }
       else if ( matchall(type, NOX|NOZ)) {
         if ( matchany(type, XDIR|ZDIR)) { type = 2; }
         else { type = 'y'; }
       }
       else if ( matchall(type, NOY|NOZ)) {
         if ( matchany(type, YDIR|ZDIR)) { type = 2; }
         else { type = 'x'; }
       }
       else if ( type & NOX) {
         if ( type & XDIR) { type = 2; }
         else if ( type & YDIR) { type = 'y'; }
         else if ( type & ZDIR) { type = 'z'; }
         else type = -'x';
       }
       else if ( type & NOY) {
         if ( type & YDIR) { type = 2; }
         else if ( type & XDIR) { type = 'x'; }
         else if ( type & ZDIR) { type = 'z'; }
         else type = -'y';
       }
       else if ( type & NOZ) {
         if ( type & ZDIR) { type = 2; }
         else if ( type & XDIR) { type = 'x'; }
         else if ( type & YDIR) { type = 'y'; }
         else type = -'z';
       }
       else if ( type & XDIR)
         type = 'x';
       else if ( type & YDIR)
         type = 'y';
       else {
         COM_assertion( type & ZDIR);
         type = 'z'; 
       }
     } // for j
   } // for it
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void determine_constraint_boundary	(	const COM::Attribute *	ctypes_faces,
		COM::Attribute *	ctypes_bndry_edges,
		COM::Attribute *	ctypes_bndry_nodes
	)

protected

Convert facial or panel constraints to nodal constraints.

Definition at line 253 of file Propagation_3.C.

References _panes, _surf, Rocblas::copy_scalar(), iend, j, k, Element_node_enumerator::next(), nj, and Element_node_enumerator::size_of_edges().

Referenced by set_constraints().

                                                                  {
   _surf->update_bd_flags( ctypes_faces); // Update facal flags along boundaries
 
   char zero = 0;
   Rocblas::copy_scalar( &zero, ctypes_bndry_edges);
   Rocblas::copy_scalar( &zero, ctypes_bndry_nodes);
 
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   std::vector< COM::Pane*>::iterator it = _panes.begin(), iend=_panes.end();
 
   // Loop through the panes and its real elements
   for ( it=_panes.begin(); it!=iend; ++it, ++pm_it) {
     COM::Pane *pane = *it;
 
     const int *val_faces=(const int*)pane->attribute(ctypes_faces->id())->pointer();
     char *val_bndry_edges=(char*)pane->attribute(ctypes_bndry_edges->id())->pointer();
     char *val_bndry_nodes=(char*)pane->attribute(ctypes_bndry_nodes->id())->pointer();
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1);
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       for ( int k=0, ne=ene.size_of_edges(); k<ne; ++k) {
         Edge_ID eid(j+1,k), eid_opp = pm_it->get_opposite_real_edge( eid);
         
         if (!pm_it->is_physical_border_edge( eid_opp)) {
           int flag_opp = (eid_opp.is_border() ? 
                           pm_it->get_bd_flag( eid_opp) : val_faces[eid_opp.eid()-1]);
           if ( val_faces[j] != flag_opp) {
             val_bndry_edges[j] |= (1<<k);
             val_bndry_nodes[ene[k]-1] = 1;
           }
         }
       } // for k
     } // for j
   } // for it
 
   _surf->reduce_on_shared_nodes( ctypes_bndry_nodes, Manifold::OP_MAX);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void enforce_nodal_constraint	(	int	type,
		Vector_3 &	du
	)

staticprotected

Enforce constraint for a specific vector.

Definition at line 441 of file Propagation_3.C.

References COM_assertion_msg.

Referenced by enforce_nodal_constraints().

                                                   {
   switch ( type) {
   case 0: break;
   case 2:    du = Vector_3(0.); break;
   case 'x':  du[1] = du[2] = 0.; break;
   case -'x': du[0] = 0.; break;
   case 'y':  du[0] = du[2] = 0.; break;
   case -'y': du[1] = 0.; break;
   case 'z':  du[0] = du[1] = 0.; break;
   case -'z': du[2] = 0.; break;
   default: COM_assertion_msg( false, "Unknown type of constraint");
   }
 }

Here is the caller graph for this function:

void enforce_nodal_constraints ( COM::Attribute * du )

virtual

Enforces the nodal constraints by projecting motion onto given direction.

Definition at line 497 of file Propagation_3.C.

References _cnstr_nodes, _cnstr_set, _panes, enforce_nodal_constraint(), i, iend, and n.

                                              {
   if ( !_cnstr_set) return;
 
   std::vector< COM::Pane*>::iterator it = _panes.begin(), iend=_panes.end();
   // Loop through the panes and its real nodes
   for ( it=_panes.begin(); it!=iend; ++it) {
     COM::Pane *pane = *it;
     
     const int *types = reinterpret_cast<const int*>
       (pane->attribute(_cnstr_nodes->id())->pointer());
     Vector_3 *ds = reinterpret_cast<Vector_3*>
       (pane->attribute( du->id())->pointer());
 
     for (int i=0, n=pane->size_of_real_nodes(); i<n; ++i) {
       enforce_nodal_constraint( types[i], ds[i]);
     }
   }
 }

Here is the call graph for this function:

void get_constraint_directions	(	int	type,
		int &	ndirs,
		Vector_3	dirs[2]
	)

staticprotected

Get orthonormals of the constraints.

If nnrms is 0, then not constrained. If nnrms is 3, then the point is fixed. Otherwise, the point is constrained either in a plane or a line.

Definition at line 456 of file Propagation_3.C.

References COM_assertion_msg.

Referenced by FaceOffset_3::obtain_constrained_directions().

                                                                    {
   switch ( type) {
   case 0: 
     ndirs = 0;
     break;
   case 2: 
     ndirs = 3;
     break;
   case 'x':  
     ndirs = 1; 
     dirs[0] = Vector_3(1,0,0); 
     break;
   case -'x': 
     ndirs = 2; 
     dirs[0] = Vector_3(0,1,0); 
     dirs[1] = Vector_3(0,0,1); 
     break;
   case 'y': 
     ndirs = 1; 
     dirs[0] = Vector_3(0,1,0); 
     break;
   case -'y':
     ndirs = 2;
     dirs[0] = Vector_3(1,0,0); 
     dirs[1] = Vector_3(0,0,1); 
     break;
   case 'z': 
     ndirs = 1; 
     dirs[0] = Vector_3(0,0,1); 
     break;
   case -'z':
     ndirs = 2;
     dirs[0] = Vector_3(1,0,0); 
     dirs[1] = Vector_3(0,1,0); 
     break;
   default: COM_assertion_msg( false, "Unknown type of constraint");
   }
 }

Here is the caller graph for this function:

static bool in_bounding_box	(	const Point_3 &	pnt,
		const Point_3 &	lb,
		const Point_3 &	ub
	)

inlinestaticprotected

Definition at line 147 of file Propagation_3.h.

Referenced by bound_nodal_motion(), and reached_nodal_bound().

                                                                      {
     return pnt[0] >= lb[0] && pnt[0] <= ub[0] &&
       pnt[1] >= lb[1] && pnt[1] <= ub[1] &&
       pnt[2] >= lb[2] && pnt[2] <= ub[2];
   }

Here is the caller graph for this function:

bool reached_nodal_bound	(	const Point_3 &	pnt,
		const double *	bnd,
		double	eps = `0`
	)

staticprotected

Definition at line 418 of file Propagation_3.C.

References check_axial_bound(), check_radial_bound(), check_spherical_bound(), COM_assertion_msg, and in_bounding_box().

Referenced by bound_facial_speed().

                                                                         {
   if ( !in_bounding_box( pnt, (const Point_3&)bnd[6], (const Point_3&)bnd[9]))
     return false;
 
   switch ( (int)bnd[0]) {
   case 's':  return check_spherical_bound( pnt, (const Point_3&)bnd[3],
                                            bnd[1], bnd[2], eps);
 
   case 'x':  return check_axial_bound( pnt[0]-bnd[3], bnd[1], bnd[2], eps); 
   case 'y':  return check_axial_bound( pnt[1]-bnd[4], bnd[1], bnd[2], eps); 
   case 'z':  return check_axial_bound( pnt[2]-bnd[5], bnd[1], bnd[2], eps); 
 
   case -'x': return check_radial_bound( pnt[1]-bnd[4], pnt[2]-bnd[5], bnd[1], bnd[2], eps); 
   case -'y': return check_radial_bound( pnt[0]-bnd[3], pnt[2]-bnd[5], bnd[1], bnd[2], eps); 
   case -'z': return check_radial_bound( pnt[0]-bnd[3], pnt[1]-bnd[4], bnd[1], bnd[2], eps); 
     
   default: COM_assertion_msg( false, "Unknown type of bnd");
   }
 
   return false;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void set_bounds ( const COM::Attribute * bnd )

Set the bounds.

Definition at line 58 of file Propagation_3.C.

References _bnd_set, _buf, BOUND_LEN, COM_assertion_msg, COM_compatible_types(), and COM_DOUBLE.

                                      {
   if ( _buf == NULL) return;
 
   if ( bnd == NULL) {
     // Deallocate cnstr_bound if constraints are unset
     _buf->dealloc_array( "PROP_cnstr_bound", 0);
     _bnd_set = false;
   }
   else {
     COM_assertion_msg( COM_compatible_types( bnd->data_type(), COM_DOUBLE) &&
                        bnd->size_of_components()==BOUND_LEN && 
                        bnd->is_panel() || bnd->is_windowed(),
                        "Propagation_3::set_bounds expects double-precision 3-vector");
 
     _buf->inherit( const_cast<COM::Attribute*>(bnd), "PROP_cnstr_bound",
                    COM::Pane::INHERIT_CLONE, true, NULL, 0);
     _bnd_set = true;
   }
   
   _buf->init_done( false);
 }

Here is the call graph for this function:

void set_constraints ( const COM::Attribute * cnstr_types )

virtual

Set the types and directions of nodal constraints.

Parameters

cnstr_types is an integer-type attribute. If it is a nodal attribute, then the value for each node is one of the following: 0: No constraint. 1: Move along the direction to be given in the displacement vector when calling

See Also: time_stepping(). -1: Move orthogonal to the direction to be given in the displacement vector when calling; time_stepping(). 2: Fixed point. 'x': Move along x direction. 'y': Move along y direction. 'z': Move along z direction. -'x': Move orthogonal to x direction (in the yz plane). -'y': Move orthogonal to y direction (in the xz plane). -'z': Move orthogonal to z direction (in the xy plane). 't': Move orthogonal to normal direction.

If the attribute is facial, panel, or windowed, then its value cannot be 1 or -1 but can be any other value as above. These constraints will be converted into nodal constraints by

See Also: convert_constraints().

Definition at line 81 of file Propagation_3.C.

References _buf, _cnstr_bndry_edges, _cnstr_bndry_nodes, _cnstr_faces, _cnstr_nodes, _cnstr_set, COM_assertion_msg, convert_constraints(), Rocblas::copy(), and determine_constraint_boundary().

                                             {
   if ( _buf == NULL) return;
 
   if ( ct_in == NULL) {
     // Deallocate cnstr_types if constraints are unset
     _buf->dealloc_array( "PROP_cnstr_nodes", 0);
     _buf->dealloc_array( "PROP_cnstr_faces", 0);
     _buf->dealloc_array( "PROP_cnstr_bndry_edges", 0);
     _buf->dealloc_array( "PROP_cnstr_bndry_nodes", 0);
     _cnstr_set = false;
   }
   else {
     COM_assertion_msg( ct_in->is_panel() || ct_in->is_elemental(), 
                        "Constraints must be associated with panes or faces.");
 
     // Copy the constraints onto the window
     _buf->resize_array( _cnstr_faces, 0);
     Rocblas::copy( ct_in, _cnstr_faces);   // Set facal flags
 
     // Determine the boundary of patches with constraints
     _buf->resize_array( _cnstr_bndry_edges, 0);
     _buf->resize_array( _cnstr_bndry_nodes, 0);
     determine_constraint_boundary( _cnstr_faces, _cnstr_bndry_edges, _cnstr_bndry_nodes);
 
     // Determine nodal constraints.
     _buf->resize_array( _cnstr_nodes, 0);
 
     // Convert constraints from facial/panel to nodal
     convert_constraints( _cnstr_faces, _cnstr_nodes);
     _cnstr_set = true;
   }
 }

Here is the call graph for this function:

void set_verbose ( bool b )

inline

Set the verbose level.

Definition at line 102 of file Propagation_3.h.

References _verb.

102 { _verb = b; }

Propagation_3::_verb

bool _verb

Definition: Propagation_3.h:162

static double square ( double x )

inlinestaticprotected

Definition at line 155 of file Propagation_3.h.

References x.

Referenced by bound_radial_disp(), bound_spherical_disp(), check_radial_bound(), check_spherical_bound(), FaceOffset_3::reset_default_parameters(), and FaceOffset_3::set_fangle_weak().

155 { return x*x; }

x

void int int REAL * x

Definition: read.cpp:74

Here is the caller graph for this function:

virtual double time_stepping	(	const COM::Attribute *	spd,
		double	dt,
		COM::Attribute *	disps,
		int *	smoothed = `NULL`
	)

pure virtual

Main entry for invoking a surface propagation algorithm.

Note that disps serves as both input and output. At input, it contains the (normalized) constrained nodal direction for the vertices with a customized directional or planar constraint (