pfs.c

Go to the documentation of this file.

/*
** ---------------------------------------------------------------
**
** @name Parametrized Faceted Surface module
**
** ---------------------------------------------------------------
**
** @date    20-Nov-2006
** @author  Luiz F. Martha
** @version 1.0
**
** ---------------------------------------------------------------
** pfs.c - Parametrized Faceted Surface module.
**
** Description:
** ---------------------------------------------------------------
**
** This module contains functions to represent a typical faceted 
** (triangulated) surface mesh data and functions to create a 
** two-dimensional parametric space associated to the faceted 
** surface.
** Basically, a pair of surface parametric values are associated 
** to each node in the mesh.  This is mainly based on the wrk
** of Levy, Petitjean, Ray, and Maillot, "Least squares conformal
** maps for automatic texture atlas generation," SIGRAPH 
** proceedings, 2002
** The local data structure consists of a node (vertex) vector
** with nodal coordinates and node-node adjacency relationships,
** and an element (face) vector with nodal connectivity.
**
** The API of this module is described in file "pfs.h".
**
** ---------------------------------------------------------------
**
** Created:    17-Jan-2005    Luiz F. Martha
**    Stolen from old pfs module based on the work of
**    Michael Floater: "Parametrization and smooth approximation 
**    of surface triangulation", Computer Aided Geometric Design, 
**    vol. 14, pp. 231-250, 1997.
**    Based on ArbSurfaceMapping.cpp by Wash.
**
** Modified:   17-Feb-2006    Alexandre A. Del Savio
**    Created the functions pfsProjectOn, pfsDistElem2Point,
**    pfsGetParPoint, pfsMapUV and pfsArea3d.
**    These functions were based on nSurfacePFS package, by 
**    William Lira and Antonio Miranda.
**
** Modified:   20-Feb-2006    Alexandre A. Del Savio
**    Created the functions pfsGetJacobian, pfsUnMapUV,
**    pfsGetJac, pfsArea2d, pfsIsParPtOnBound, pfsIsParPtInside, 
**    and pfsSnapParPtToBound.
**
** Modified:   22-Feb-2006    Alexandre A. Del Savio
**    Modified the function pfsGetJacobian. It passed to
**    return the components of the normal vector of  
**    element face where is the point.
**
** Modified:   22-Feb-2006    William W. M. Lira
**    Created the functions pfsEval and pfsGet3dPoint.
**    These functions were based on nSurfacePFS package, by 
**    William Lira and Antonio Miranda.
**
** Modified:   23-Feb-2006    Alexandre A. Del Savio
**    Modified the function pfsAddElem. It passed to have
**    a return. If there is no active surface it returns a
**    false (0) status, otherwise it returns a true (1) 
**    status.
**
** Modified:   27-Jul-2006    Luiz F. Martha
**    Created local function pfsNodeConstrNormal.
**    Modified function pfsAddAdjNode to order adjacent 
**    nodes around a reference node taking into account an
**    average normal vector that considers the the given 
**    element normal (if the element is not null) and
**    normal vectors of the adjacent elements in the list 
**    of adjacent nodes (it uses the created function 
**    pfsNodeConstrNormal).
**
** Modified:   30-Jul-2006    Luiz F. Martha
**    Created exported functions pfsItrNumLoops, pfsItrFirstLoopNode,
**    pfsItrNextLoopNode.
**    Created vector of first loop's node.
**    Modified functions pfsIsParPtOnBound, pfsIsParPtInside, and 
**    pfsSnapParPtToBound such that the use loop data structure 
**    instead of the dces package.
**    Created NODE_MARKED type.
**    Created local functions pfsMarkBdryNodes, pfsUnmarkLoopNodes,
**    and pfsBuildLoopList.
**
** Modified:   31-Aug-2006    Luiz F. Martha
**    Created exported function pfsCompleteParSpace to compute 
**    parametric bounding boxes of elements, to build the R-tree 
**    of elements indexed by parametric bounding boxes and to compute 
**    node Jacobians.
**    Modified function pfsSolveSurfPar to not perform these tasks.
**    Created local function pfsNormalizeParVals and got rid of local
**    functions pfsMapUV and pfsUnMapUV.
**    Replaced function pfsDistElem2Point by function pfsDistPntCurElem.
**    Replaced function pfsGetParPoint by function pfsParPntCurElem.
**    Modified function pfsProjectOn to use functions pfsDistPntCurElem
**    and pfsParPntCurElem.
**    Replaced function pfsGet3dPoint by function pfs3dPntCurElem.
**    Modified function pfsEval to use pfs3dPntCurElem.
**    Replaced function pfsGetJac by function pfsParDerivCurElem.
**    Modified function pfsGetJacobian to use pfsParDerivCurElem.
**    Created exported function pfsSetCurrElem.
**
** Modified:   20-Sep-2006    Gisele Cunha Holtz
**    Created exported function pfsShot and all the related 
**    local functions.
**
** Modified:   26-Sep-2006    Luiz F. Martha
**    Modified names of functions and did some clean up.
**
** Modified:   05-Oct-2006    Gisele Cunha Holtz
**    Criada a função pfsIsPntInElem. Modificada a função pfsArea3d,
**    que passou a receber a normal da area e a retornar a área com
**    sinal de acordo com a normal recebida.
**
** Modified:   20-Nov-2006    Luiz F. Martha
**    Created exported function pfsSnapToBound.
**    Modified functions pfsGetJacobian and pfsEval such that they
**    reject points outside current active surface.
**    Did some clean up.
**
** Modified:   09-July-2007    Gisele Cunha Holtz
**    Foram realizadas alterações para evitar problemas com a 
**    tolerância:
**    Modificada a função pfsParPntCurElem para retornar o status
**    de que um ponto está dentro de um elemento apenas quando
**    o ponto estiver dentro no espaço cartesiano e no paramétrico.
**    A função pfsUpdateShotLen passou a verificar se o 
**    comprimento restante do tiro é menor que a tolerância.
**    Na função pfsChkProxElemEdge foi feita uma pequena alteração
**    na verificação se o ponto está sobre um nó. Além disso,
**    adicionado um teste para certificar que a projeção do ponto
**    sobre a aresta está realmente sobre a aresta.
**    Criada a função pfsParPntCurElemAttract que é similiar a 
**    pfsParPntCurElem, exceto por corrigir a coordenada do ponto
**    qdo o ponto está dentro do triângulo. 
**    Removida a funçao pfsIsPntInElem e substituídas suas chamadas por
**    pfsParPntCurElem. Criadas as funções pfs3dGetAreaCoord e 
**    pfsParGetAreaCoord.
*/

#define _PFS_C

/*
** ---------------------------------------------------------------
** Global variables and symbols:
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "geo.h"
#include "pfs.h"
#include "amr3b.h"

/*
** ---------------------------------------------------------------
** Local definitions and variables:
*/

#define ITER_TOLER 0.0001

#define TOLER 1e-5

#define COL_TOLER 1e-5 

#ifndef ABS
#define ABS(x) (((x) < 0.0)? -(x): (x))
#endif

#ifndef MAX
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#endif

#ifndef Len2
#define Len2(x,y) sqrt((double)(((x)*(x))+((y)*(y))))
#endif

#ifndef Len3
#define Len3(x,y,z) sqrt((double)(((x)*(x))+((y)*(y))+((z)*(z))))
#endif

#define DET2D(a,b,c,d)   (((a)*(d))-((b)*(c)))

#define X_CROSS_LINE_LINE(x1,y1,x2,y2,x3,y3,x4,y4) \
        DET2D( DET2D(x1,y1,x2,y2), x1-x2, DET2D(x3,y3,x4,y4), x3-x4) / \
        DET2D(x1-x2,y1-y2,x3-x4,y3-y4)
#define Y_CROSS_LINE_LINE(x1,y1,x2,y2,x3,y3,x4,y4) \
        DET2D( DET2D(x1,y1,x2,y2), y1-y2, DET2D(x3,y3,x4,y4), y3-y4) / \
        DET2D(x1-x2,y1-y2,x3-x4,y3-y4)

#define X_CROSS_LINE_SHOT(x1,y1,x2,y2,x0,y0,xS,yS) \
        DET2D( DET2D(x1,y1,x2,y2), x1-x2, DET2D(xS,yS,x0,y0), xS) / \
        DET2D(x1-x2,y1-y2,xS,yS)
#define Y_CROSS_LINE_SHOT(x1,y1,x2,y2,x0,y0,xS,yS) \
        DET2D( DET2D(x1,y1,x2,y2), y1-y2, DET2D(xS,yS,x0,y0), yS) / \
        DET2D(x1-x2,y1-y2,xS,yS)


enum {
 NODE_INTERIOR,
 NODE_BOUNDARY,
 NODE_MARKED
 };

enum {
 NODE = 1,
 EDGE = 2,
 FACE = 4,
 LOOP = 8,
 FREE_FACE = 16
 };
typedef struct _pfselem {
  int           inc[3];
  PFSGeoPoint      normal;
  PFSGeoPoint      center;
  double          area;
  double       W[3][2];
  PFSGeoPoint    min3dbox;
  PFSGeoPoint    max3dbox;
  PFSGeoSurfPar minparbox;
  PFSGeoSurfPar maxparbox;
 } PfsElem;

typedef struct _pfsadjnode {
  struct _pfsadjnode  *next;
  int                 id;
  PfsElem             *elem;
  int                 refcnt;
  double           AtA[2][2];
 } PfsAdjNode;

typedef struct _pfsnode {
  PFSGeoPoint            coord;
  PFSGeoPoint            normal;
  PFSGeoPoint            ddu;
  PFSGeoPoint            ddv;
  PFSGeoSurfPar          par;
  int                 degree;
  double              AtA[2][2];
  struct _pfsadjnode  *adj_head;
  int                 type;
 } PfsNode;

typedef struct _pfssurf {
  int          numnodes;
  int          numelems;
  int          numloops;
  PfsNode      *nodelist;
  void         *elm3d_tree;
  void         *elm2d_tree;
  int          *looplist;
  int          *sol_order;
  int maxnodes;
  double box_xmin;
  double box_xmax;
  double box_ymin;
  double box_ymax;
  double box_zmin;
  double box_zmax;
  double box_umin;
  double box_umax;
  double box_vmin;
  double box_vmax;
  PfsElem *ref_elem;
  int     ref_node;
  double  iter_toler;

} PfsSurf;

/* Active surface data pointer.
 */
static PfsSurf *pfs;

/* Index of current element and current element node for traversal.
 */
static int cur_elemnode;

static PfsElem   *cur_elem = NULL;
static PfsElem *curr_ielem = NULL;
int _loop_ = 0;

/*
** ---------------------------------------------------------------
** Private functions: 
*/
static void pfs3dGetAreaCoord( double x, double y, double z,
                               double* xn, double* yn, double* zn,
                               double* b0, double* b1, double* b2 );
static void pfsParGetAreaCoord( double u, double v,
                                double* un, double* vn, 
                                double* b0, double* b1, double* b2 );
static void pfsNodeConstrNormal( int ref_id, PfsElem *elem, PFSGeoPoint *normal );
static int  pfsAddAdjNode( int ref_id, int adj_id, PfsElem* elem );
static int  pfsBuildNodeAdj( void );
static int  pfsChkBdryNode( int id );
static void pfsClassifyNodes( void );
static int  pfsOrderBdryNodeAdj( int id );
static int  pfsOrderAllBdryNodeAdj( void );
static void pfsMarkBdryNodes( void );
static void pfsUnmarkLoopNodes( int ref_id );
static void pfsBuildLoopList( void );
static void pfsDelNodeList( PfsSurf *surf );
static void pfsDelElemList( PfsSurf *surf );
static int  pfsGetEdgeAdjElems( int ni, int nj,
                                PfsElem **e_left, PfsElem **e_right );
static void pfsCompute3dBoundBox( void );
static void pfsElemNormal( PfsElem *elem );
static void pfsElemCenter( PfsElem *elem );
static void pfsElemLocalGradients( PfsElem *elem );
static void pfsElem3dBoundBox( PfsElem *elem );
static void pfsNodeNormal( int node );
static void pfsNodeNormals( void );
static void pfsFindRefElem( void );
static void pfsNodeConfMapCoeff( int ref_id );
static void pfsConfMapCoefficients( void );
static void pfsInitParams( void );
static int  pfsCreateNextSolOrderLevel( int lev_begin, int lev_end, int *pos,
                                        int *solorder, int *nodevisited );
static void pfsBuildSolOrderVec( void );
static void pfsDelSolOrderVec( PfsSurf *surf );
static void pfsSolve2x2( double a[2][2], double b[2], double *u, double *v );
static void pfsNormalizeParVals( void );
static void pfsComputeParBoundBox( void );
static void pfsElemParBoundBox ( PfsElem *elem );
static void pfsComputeElemParBoundBoxes( void );
static void pfsComputeElemJac  ( PfsElem *elem,
                                 double *xdu, double *ydu, double *zdu, 
                                 double *xdv, double *ydv, double *zdv );
static void   pfsNodeJacobian  ( int node );
static void   pfsNodeJacobians ( void );
static double pfsArea3d        ( double xa, double ya, double za, 
                                 double xb, double yb, double zb,
                                 double xc, double yc, double zc,
                                 double xn, double yn, double zn );
static double pfsArea2d        ( double u0, double v0,
                                 double u1, double v1,
                                 double u2, double v2 );
static double pfsDistPntCurElem( PFSGeoPoint *p, PFSGeoPoint *clst );
static int    pfsParPntCurElem ( double x, double y, double z,
                                 double *u, double *v );
static int    pfs3dPntCurElem  ( double u, double v,
                                 double *x, double *y, double *z );
static int    pfsParDerivCurElem( double u, double v,
                                 double *xdu, double *ydu, double *zdu,
                                  double *xdv, double *ydv, double *zdv,
                                  double *xdw, double *ydw, double *zdw );
static int    pfsClstPtParLine( double ua, double va, double ub, double vb,
                                double u, double v, double *up, double *vp,
                                double *dist );
static int    pfsIsParPtInside    ( PFSGeoSurfPar *p );
static int    pfsIsParPtOnBound   ( PFSGeoSurfPar *p );
static int    pfsSnapParPtToBound ( PFSGeoSurfPar *p );
static int pfsChkProxElemNode( PfsElem *elem, double *u, double *v, int *ni );
static int pfsChkProxElemEdge( PfsElem *elem, double *u, double *v, int *type, 
                       int *ni, int *nj );
static int    pfsShotXLine     ( PfsElem *elem, double u0, double v0, 
                                 double u_shot, double v_shot,
                         int na, int nb, double *u1, double *v1,
                         int *type, int *n_inter );
static int pfsNodeShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1, double *v1, 
                        int *type, int *ni, int *nj );
static int pfsEdgeShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1, double *v1, 
                        int *type, int *ni, int *nj );
static int pfsFaceShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1, double *v1, 
                        int *type, int *ni, int *nj );
static int pfsLoopShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1, double *v1, 
                        int *type, int *ni, int *nj );
static int pfsUpdateShotLen( double u0, double v0, PFSGeoPoint *p0, 
                             double *u1, double *v1, PFSGeoPoint *p1, 
                             double *len );
static void pfsFinishFreeFaceShotLen( double u0, double v0, PFSGeoPoint *p0, 
                                      double u_shot, double v_shot, double len,
                                      double *u1, double *v1, PFSGeoPoint *p1);
static int pfsShotSolve( PfsElem *init_elem, double u0, double v0, 
                         PFSGeoPoint p0, double u_shot, double v_shot, 
                         PfsElem **end_elem, double *u1, double *v1,
                         PFSGeoPoint *p1, int *type, int *ni, int *nj,
                         double *len );


/*=======================  pfs3dGetAreaCoord  ========================*/
static void pfs3dGetAreaCoord( double x, double y, double z,
                               double* xn, double* yn, double* zn,
                               double* b0, double* b1, double* b2 )
{
 double d;
 PFSGeoPoint *normal = &(cur_elem->normal);

 d = pfsArea3d( xn[0], yn[0], zn[0],
                xn[1], yn[1], zn[1],
                xn[2], yn[2], zn[2],
                normal->x, normal->y, normal->z );
 *b0 = pfsArea3d( x, y, z, xn[1], yn[1], zn[1], xn[2], yn[2], zn[2],
                  normal->x, normal->y, normal->z ) / d;
 *b1 = pfsArea3d( x, y, z, xn[2], yn[2], zn[2], xn[0], yn[0], zn[0],
                  normal->x, normal->y, normal->z ) / d;
 *b2 = pfsArea3d( x, y, z, xn[0], yn[0], zn[0], xn[1], yn[1], zn[1],
                  normal->x, normal->y, normal->z ) / d;
}

/*=====================  pfsParGetAreaCoord  ======================*/
static void pfsParGetAreaCoord( double u, double v,
                                double* un, double* vn, 
                                double* b0, double* b1, double* b2 )
{
 double d;

 d = pfsArea2d(un[0], vn[0], un[1], vn[1], un[2], vn[2]);
 *b0 = pfsArea2d( u, v, un[1], vn[1], un[2], vn[2]) / d;
 *b1 = pfsArea2d( u, v, un[2], vn[2], un[0], vn[0]) / d;
 *b2 = pfsArea2d( u, v, un[0], vn[0], un[1], vn[1]) / d;
}

/*=====================  pfsNodeConstrNormal  =====================*/
static void pfsNodeConstrNormal( int ref_id, PfsElem *elem, PFSGeoPoint *normal )
{
 PfsNode    *node;               /* pointer to given reference node */
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */
 double     size;                /* normal vector size */

/* Get the pointer to the reference node.
 */
 node = &pfs->nodelist[ref_id];

/* Initialize the normal vector as the normal of the given 
 * element or as null vector if the given element is null.
 */
 if( elem == NULL )
 {
  normal->x = 0.0;
  normal->y = 0.0;
  normal->z = 0.0;
 }
 else
 {
  *normal = elem->normal;
 }

/* Compute the average normal construction vector.
 */
 for( adjnode = node->adj_head;
      adjnode != NULL;
      adjnode = adjnode->next )
 {
  if( adjnode->elem != NULL )
  {
   normal->x += adjnode->elem->normal.x;
   normal->y += adjnode->elem->normal.y;
   normal->z += adjnode->elem->normal.z;
  }
 }

 size = PFSGeoVecLen( normal );
 if( size > 0.0 )
 {
  normal->x /= size;
  normal->y /= size;
  normal->z /= size;
 }
}

/*========================  pfsAddAdjNode  ========================*/
static int  pfsAddAdjNode( int ref_id, int adj_id, PfsElem* elem )
{
 PfsNode    *node;               /* pointer to given reference node */
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */
 PfsAdjNode *newadj;             /* new adjacent node created */
 PfsNode    *first_adjnode;      /* first adjacent node of ref. node */
 PfsAdjNode *prevadjnode;        /* pointer to previous adjacent node */
 double     cur_angle;           /* current adj.node angle around ref. node */
 double     new_angle;           /* new adj. node angle around ref. node */
 PFSGeoLine    nline;               /* normal line at reference point */

/* Get the pointer to the reference node.
 */
 node = &pfs->nodelist[ref_id];

/* Check to see whether the given node (given by its id) is already
 * in the list of adjacent nodes of the reference node.
 */
 for( adjnode = node->adj_head;
      adjnode != NULL;
      adjnode = adjnode->next )
 {
  if( adjnode->id == adj_id )
  {
/* If the given adjancent node was already in the list, 
 * check for more than two elements using a pair of nodes and
 * for reversed element incidence order.
 */
   if( adjnode->refcnt == 2 )
   {
    printf( "\n Found more than two elements framing into an edge!\n"
            " Ref. node: %d, Adj. node: %d, Adj. elem.:%p \n", 
            ref_id, adj_id, elem  );
    return( PFS_NON_MANIFOLD_EDGE );
   }
   if( adjnode->elem == NULL )
   {
    if( elem == NULL )
    {
     printf( "\n Found reversed element incidence order!\n"
             " Ref. node: %d, Adj. node: %d, Adj. elem.:%p \n", 
             ref_id, adj_id, elem  );
     return( PFS_INVERSE_INCIDENCE );
    }
    else
    {
     adjnode->elem = elem;
    }
   }
   else
   {
    if( elem != NULL )
    {
     printf( "\n Found reversed element incidence order!\n"
             " Ref. node: %d, Adj. node: %d, Adj. elem.:%p \n", 
             ref_id, adj_id, elem  );
     return( PFS_INVERSE_INCIDENCE );
    }
   }

/* If no problem is found, just increment its reference counter
 * and updates the index of the adjacent element.
 */
   adjnode->refcnt++;
   return( PFS_SUCCESS );
  }
 }

/* Update the number of adjacent nodes of the reference node.
 */
 node->degree += 1;

/* Get memory of new adjacent node and initialize it.
 */
 newadj = (PfsAdjNode *)malloc( sizeof( PfsAdjNode ) );
 newadj->id = adj_id;
 newadj->elem = elem;
 newadj->refcnt = 1;
 newadj->AtA[0][0] = 0.0;
 newadj->AtA[0][1] = 0.0;
 newadj->AtA[1][0] = 0.0;
 newadj->AtA[1][1] = 0.0;

/* For a new adjacent node, check to see whether the 
 * list is empty.  If that is the case, just 
 * initialize the list with the new adjacent node.
 * Otherwise, put it in the list ordering the adjacent
 * node in counter-clockwise order radially around the 
 * reference node, starting from the first adjacent node.
 * In this case, build a node construction normal based on
 * the given element normal (if the element is not null) and
 * on the normal vectors of the adjacent elements in the list
 * of adjacent nodes.
 */
 if( node->adj_head == NULL )
 {
  newadj->next = node->adj_head;
  node->adj_head = newadj;
 }
 else
 {
  first_adjnode = &pfs->nodelist[node->adj_head->id];
  nline.d = node->coord;
  pfsNodeConstrNormal( ref_id, elem, &nline.e );
  new_angle = PFSGeoAngPtAboutLine( &nline, 
                                 &pfs->nodelist[adj_id].coord, 
                                 &first_adjnode->coord );
  prevadjnode = node->adj_head;
  cur_angle = 0.0;
  adjnode = prevadjnode->next;
  while( adjnode != NULL )
  {
   cur_angle = PFSGeoAngPtAboutLine( &nline, 
                                  &pfs->nodelist[adjnode->id].coord, 
                                  &first_adjnode->coord );
   if( new_angle < cur_angle )
    break;
   prevadjnode = adjnode;
   adjnode = prevadjnode->next;
  }
  newadj->next = prevadjnode->next;
  prevadjnode->next = newadj;
 }
 return( PFS_SUCCESS );
}

/*=======================  pfsBuildNodeAdj  =======================*/
static int pfsBuildNodeAdj( void )
{
 int       node;                  /* current element node id     */
 int       prev;                  /* previous node on elem. bdry */
 int       next;                  /* next node on element boudry */
 int       j;                     /* counter                     */
 int       status;

 double xmin, xmax;
 double ymin, ymax;
 double zmin, zmax;

 void* e = NULL;

/* Traverse the list of elements in mesh.
 */
 PFSRtreeInitTraverse( pfs->elm3d_tree );

 while( (e = PFSRtreeTraverse( pfs->elm3d_tree,
                               &xmin, &xmax, 
                               &ymin, &ymax,
                               &zmin, &zmax )) != NULL )
 {
   PfsElem* elem = (PfsElem *)e;

/* Traverse the nodes of current element.
 */
  for( j = 0; j < 3; j++ )
  {
/* Get current, previous and next node on element boundary.
 */
   node = elem->inc[(j+1)%3];
   prev = elem->inc[j];
   next = elem->inc[(j+2)%3];

/* Add previous and next nodes to the linked list of adjacent
 * nodes of current node, or just increment their reference 
 * counters.
 * Set current element as adjacent element of current node 
 * only for next node.  By doing this, this element will
 * appear only once as an adjacent element of this node.
 */
   status = pfsAddAdjNode( node, prev, NULL );
   if( status != PFS_SUCCESS )
    return( status );
   status = pfsAddAdjNode( node, next, elem );
   if( status != PFS_SUCCESS )
    return( status );
  }
 }
 return( PFS_SUCCESS );
}

/*=======================  pfsChkBdryNode  ========================*/
static int pfsChkBdryNode( int id )
{
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */

 for( adjnode = pfs->nodelist[id].adj_head;
      adjnode != NULL;
      adjnode = adjnode->next )
 {
  if( adjnode->refcnt != 2 )
  {
   return( 1 );
  }
 }
 return( 0 );
}

/*======================  pfsClassifyNodes  =======================*/
static void pfsClassifyNodes( void )
{
 int        id;                    /* current node index           */

 for( id = 0; id < pfs->numnodes; id++ )
 {
  if( pfsChkBdryNode( id ) )
   pfs->nodelist[id].type = NODE_BOUNDARY;
  else
   pfs->nodelist[id].type = NODE_INTERIOR;
 }
}

/*=====================  pfsOrderBdryNodeAdj  =====================*/
static int pfsOrderBdryNodeAdj( int ref_id )
{
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */
 PfsAdjNode *curlast = NULL;     /* pointer to current last adj. node in list */
 PfsAdjNode *adjlast = NULL;     /* pointer to target last adj. node in list */

 for( adjnode = pfs->nodelist[ref_id].adj_head;
      adjnode != NULL;
      adjnode = adjnode->next )
 {
  if( adjnode->elem == NULL )
  {
   if( adjlast != NULL )
   {
    printf( "\n Found inconsistent boundary neighborhood!\n"
            " Ref. boundary node: %d\n", ref_id );
    return( PFS_NON_MANIFOLD_VERTEX );
   }
   else
   {
    adjlast = adjnode;
   }
  }
  if( adjnode->next == NULL )
   curlast = adjnode;
 }

/* If the given node is really a boundary node, there should be
 * an adjacent node with a field elem = NULL.  If no one is found
 * something wrong is going one.  In this case (which is really
 * a bug), send a message of bad adjacency neighborhood and return
 * an invalid status.
 */
 if( adjlast == NULL )
 {
  printf( "\n Found inconsistent boundary neighborhood!\n"
          " Ref. boundary node: %d\n", ref_id );
  return( PFS_NON_MANIFOLD_VERTEX );
 }

/* If target last adjacent node is already the last one,
 * there is nothing to do.  Just return a valid status.
 */
 if( adjlast == curlast )
  return( PFS_SUCCESS );

/* In the cannonical case, make the current last node precede 
 * the current head of list, and make next adjacent node be 
 * the new head of the list.  This will also consider the case
 * in which the target last adjacent node is the current 
 * head of the list.
 */
 curlast->next = pfs->nodelist[ref_id].adj_head;
 pfs->nodelist[ref_id].adj_head = adjlast->next;
 adjlast->next = NULL;
 
 return( PFS_SUCCESS );
}

/*===================  pfsOrderAllBdryNodeAdj  ====================*/
static int pfsOrderAllBdryNodeAdj( void )
{
 int       id;
 int       status;

 if( pfs->nodelist != NULL )
 {
  for( id = 0; id < pfs->numnodes; id++ )
  {
   if( pfs->nodelist[id].type != NODE_INTERIOR )
   {
    status = pfsOrderBdryNodeAdj( id );
    if( status != PFS_SUCCESS )
     return( status );
   }
  }
 }
 return( PFS_SUCCESS );
}

/*======================  pfsMarkBdryNodes  =======================*/
static void pfsMarkBdryNodes( void )
{
 int       id;

 if( pfs->nodelist != NULL )
 {
  for( id = 0; id < pfs->numnodes; id++ )
  {
   if( pfs->nodelist[id].type == NODE_BOUNDARY )
   {
    pfs->nodelist[id].type = NODE_MARKED;
   }
  }
 }
}

/*=====================  pfsUnmarkLoopNodes  ======================*/
static void pfsUnmarkLoopNodes( int ref_id )
{
 PfsNode    *ref_node;           /* pointer to given reference node */
 PfsNode    *node;               /* pointer to current node on loop */

 if( pfs->nodelist[ref_id].type != NODE_MARKED )
  return;

 ref_node = node = &pfs->nodelist[ref_id];

 do
 {
  node->type = NODE_BOUNDARY;
  node = &pfs->nodelist[node->adj_head->id];
 } while( node != ref_node );
}

/*======================  pfsBuildLoopList  =======================*/
static void pfsBuildLoopList( void )
{
 int       *looplist;   /* auxiliar loop list (stores first node) */
 int       found_loop;  /* flag for found next loop */
 int       id;          /* node id */
 int       j;           /* loop counter */

 pfs->numloops = 0;

 if( pfs->nodelist == NULL )
  return;

 if( pfs->numnodes == 0 )
  return;

 if( pfs->numelems == 0 )
  return;

 if( pfs->looplist != NULL )
  free( pfs->looplist );

 looplist = (int *)calloc( pfs->numnodes, sizeof( int ) );

 pfsMarkBdryNodes( );

 do
 {
  found_loop = 0;
  for( id = 0; id < pfs->numnodes; id++ )
  {
   if( pfs->nodelist[id].type == NODE_MARKED )
   {
    found_loop = 1;
    looplist[pfs->numloops++] = id;
    pfsUnmarkLoopNodes( id );
    break;
   }
  }
 } while( found_loop );

 pfs->looplist = (int *)calloc( pfs->numloops, sizeof( int ) );
 for( j = 0; j < pfs->numloops; j++ )
 {
  pfs->looplist[j] = looplist[j];
 }

 free( looplist );
}

/*=======================  pfsDelNodeList  ========================*/
static void pfsDelNodeList( PfsSurf *surf )
{
 PfsAdjNode *adjnode;    /* auxiliar adjacent node pointer   */
 int        i;           /* counter                          */

 if( surf->nodelist != NULL )
 {
  for( i = 0; i < surf->numnodes; i++ )
  {
   while( surf->nodelist[i].adj_head )
   {
    adjnode = surf->nodelist[i].adj_head;
    surf->nodelist[i].adj_head = adjnode->next;
    free( adjnode );
   }
  }
  free( surf->nodelist );
 }
 surf->nodelist = NULL;
}

/*=======================  pfsDelElemList  ========================*/
static void pfsDelElemList( PfsSurf *surf )
{
 double xmin, ymin, zmin;
 double xmax, ymax, zmax;
 void   *e = NULL;

 if( surf->elm2d_tree != NULL )
 {
  PFSRtreeDestroy( surf->elm2d_tree );
 }

 if( surf->elm3d_tree != NULL )
 {
  PFSRtreeInitTraverse( surf->elm3d_tree );

  while( (e = PFSRtreeTraverse( surf->elm3d_tree,
                                &xmin, &xmax, 
                                &ymin, &ymax,
                                &zmin, &zmax )) != NULL )
  {
   PfsElem* curr_elem = (PfsElem *)e;
   free( curr_elem );
  }

  PFSRtreeDestroy( surf->elm3d_tree );
 }

 surf->elm2d_tree = NULL;
 surf->elm3d_tree = NULL;
}

/*======================  pfsGetEdgeAdjElems  =====================*/
static int pfsGetEdgeAdjElems( int ni, int nj,
                               PfsElem **e_left, PfsElem **e_right )
{
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */
 int        found_node;

 found_node = 0;
 for( adjnode = pfs->nodelist[ni].adj_head;
      adjnode != NULL;
      adjnode = adjnode->next )
 {
  if( adjnode->id == nj )
  {
   *e_left = adjnode->elem;
   found_node = 1;
   break;
  }
 }

 if( !found_node )
  return( 0 );

 found_node = 0;
 for( adjnode = pfs->nodelist[nj].adj_head;
      adjnode != NULL;
      adjnode = adjnode->next )
 {
  if( adjnode->id == ni )
  {
   *e_right = adjnode->elem;
   found_node = 1;
   return( 1 );
  }
 }

 return( 0 );
}

/*=====================  pfsCompute3dBoundBox  ====================*/
static void pfsCompute3dBoundBox( void )
{
 int       i;                              /* counter   */

 if( pfs->numnodes < 1 )
  return;

 pfs->box_xmin = pfs->box_xmax = pfs->nodelist[0].coord.x;
 pfs->box_ymin = pfs->box_ymax = pfs->nodelist[0].coord.y;
 pfs->box_zmin = pfs->box_zmax = pfs->nodelist[0].coord.z;

 for( i = 1; i < pfs->numnodes; i++ )
 {
  if( pfs->box_xmin > pfs->nodelist[i].coord.x ) 
   pfs->box_xmin = pfs->nodelist[i].coord.x;
  if( pfs->box_xmax < pfs->nodelist[i].coord.x ) 
   pfs->box_xmax = pfs->nodelist[i].coord.x;
  if( pfs->box_ymin > pfs->nodelist[i].coord.y ) 
   pfs->box_ymin = pfs->nodelist[i].coord.y;
  if( pfs->box_ymax < pfs->nodelist[i].coord.y ) 
   pfs->box_ymax = pfs->nodelist[i].coord.y;
  if( pfs->box_zmin > pfs->nodelist[i].coord.z ) 
   pfs->box_zmin = pfs->nodelist[i].coord.z;
  if( pfs->box_zmax < pfs->nodelist[i].coord.z ) 
   pfs->box_zmax = pfs->nodelist[i].coord.z;
 }
}

/*=========================  pfsElemNormal  =======================*/
static void pfsElemNormal( PfsElem *elem )
{
 int      i, first, curr, prev;
 double   size;

 elem->normal.x = 0.0;
 elem->normal.y = 0.0;
 elem->normal.z = 0.0;

 first = elem->inc[0];

 for( i = 2; i < 3; i++ )
 {
  curr = elem->inc[i];
  prev = elem->inc[i-1];

  elem->normal.x +=
        (pfs->nodelist[prev].coord.y - pfs->nodelist[first].coord.y) *
        (pfs->nodelist[curr].coord.z - pfs->nodelist[first].coord.z) -
        (pfs->nodelist[prev].coord.z - pfs->nodelist[first].coord.z) *
        (pfs->nodelist[curr].coord.y - pfs->nodelist[first].coord.y);
  elem->normal.y +=
       -(pfs->nodelist[prev].coord.x - pfs->nodelist[first].coord.x) *
        (pfs->nodelist[curr].coord.z - pfs->nodelist[first].coord.z) +
        (pfs->nodelist[prev].coord.z - pfs->nodelist[first].coord.z) *
        (pfs->nodelist[curr].coord.x - pfs->nodelist[first].coord.x);
  elem->normal.z +=
        (pfs->nodelist[prev].coord.x - pfs->nodelist[first].coord.x) *
        (pfs->nodelist[curr].coord.y - pfs->nodelist[first].coord.y) -
        (pfs->nodelist[prev].coord.y - pfs->nodelist[first].coord.y) *
        (pfs->nodelist[curr].coord.x - pfs->nodelist[first].coord.x);
 }

 size = PFSGeoVecLen( &elem->normal );

 if( size > 0.0 )
 {
  elem->normal.x /= size;
  elem->normal.y /= size;
  elem->normal.z /= size;
  elem->area = size * 0.5;
 }
}

/*=========================  pfsElemCenter  ========================*/
static void pfsElemCenter( PfsElem *elem )
{
 int node;
 int i;

 elem->center.x = 0.0;
 elem->center.y = 0.0;
 elem->center.z = 0.0;

 for( i = 0; i < 3; i++ )
 {
  node = elem->inc[i];
  elem->center.x += pfs->nodelist[node].coord.x;
  elem->center.y += pfs->nodelist[node].coord.y;
  elem->center.z += pfs->nodelist[node].coord.z;
 }

 elem->center.x /= 3.0;
 elem->center.y /= 3.0;
 elem->center.z /= 3.0;
}

/*====================  pfsElemLocalGradients  ====================*/
static void pfsElemLocalGradients( PfsElem *elem )
{
 int    v0;
 int    v1;
 int    v2;
 double r0;
 double r1;
 double r2;
 double r00;
 double r11;
 double r22;
 double x1;
 double x2;
 double y2;
 double sqrt_dt;

 v0 = elem->inc[0];
 v1 = elem->inc[1];
 v2 = elem->inc[2];

 r0 = Len3( pfs->nodelist[v0].coord.x - pfs->nodelist[v2].coord.x, 
            pfs->nodelist[v0].coord.y - pfs->nodelist[v2].coord.y, 
            pfs->nodelist[v0].coord.z - pfs->nodelist[v2].coord.z );
 r1 = Len3( pfs->nodelist[v2].coord.x - pfs->nodelist[v1].coord.x, 
            pfs->nodelist[v2].coord.y - pfs->nodelist[v1].coord.y, 
            pfs->nodelist[v2].coord.z - pfs->nodelist[v1].coord.z );
 r2 = Len3( pfs->nodelist[v1].coord.x - pfs->nodelist[v0].coord.x, 
            pfs->nodelist[v1].coord.y - pfs->nodelist[v0].coord.y, 
            pfs->nodelist[v1].coord.z - pfs->nodelist[v0].coord.z );
 r00 = r0*r0;
 r11 = r1*r1;
 r22 = r2*r2;

/* x0 = 0.0;
   y0 = 0.0; */
   x1 = r2;
/* y1 = 0.0; */
   x2 = 0.5 * (r22 + r00 - r11) / r2;
   y2 = 0.5 * sqrt(2.0*(r00*r22 + r22*r11 + r00*r11) -
                        r00*r00 - r11*r11 - r22*r22) / r2;

/* dt = (x0*y1 - y0*x1) + (x1*y2 - y1*x2) + (x2*y0 - y2*x0)
 */
 sqrt_dt = sqrt( x1*y2 );

/* [W]*sqrt_dt = |(x2-x1)   (y2-y1)|
 *               |(x0-x2)   (y0-y2)|
 *               |(x1-x0)   (y1-y0)|
 */

 elem->W[0][0] = (x2-x1) / sqrt_dt;
 elem->W[0][1] =  y2     / sqrt_dt;
 elem->W[1][0] = -x2     / sqrt_dt;
 elem->W[1][1] = -y2     / sqrt_dt;
 elem->W[2][0] =  x1     / sqrt_dt;
 elem->W[2][1] =  0.0;
}

/*=======================  pfsElem3dBoundBox  ========================*/
static void pfsElem3dBoundBox( PfsElem *elem )
{
 int    node;
 int    i;

 node = elem->inc[0];

 elem->min3dbox.x =
 elem->max3dbox.x = pfs->nodelist[node].coord.x;
 elem->min3dbox.y =
 elem->max3dbox.y = pfs->nodelist[node].coord.y;
 elem->min3dbox.z =
 elem->max3dbox.z = pfs->nodelist[node].coord.z;

 for( i = 1; i < 3; i++ )
 {
  node = elem->inc[i];

  if( pfs->nodelist[node].coord.x < elem->min3dbox.x )
      elem->min3dbox.x = pfs->nodelist[node].coord.x;

  if( pfs->nodelist[node].coord.y < elem->min3dbox.y )
      elem->min3dbox.y = pfs->nodelist[node].coord.y;

  if( pfs->nodelist[node].coord.z < elem->min3dbox.z )
      elem->min3dbox.z = pfs->nodelist[node].coord.z;

  if( pfs->nodelist[node].coord.x > elem->max3dbox.x )
      elem->max3dbox.x = pfs->nodelist[node].coord.x;

  if( pfs->nodelist[node].coord.y > elem->max3dbox.y )
      elem->max3dbox.y = pfs->nodelist[node].coord.y;

  if( pfs->nodelist[node].coord.z > elem->max3dbox.z )
      elem->max3dbox.z = pfs->nodelist[node].coord.z;
 }

 PFSRtreeInsert( pfs->elm3d_tree, (void *)elem, elem->min3dbox.x,
                                                elem->max3dbox.x,
                                                elem->min3dbox.y,
                                                elem->max3dbox.y,
                                                elem->min3dbox.z,
                                                elem->max3dbox.z );
}

/*=========================  pfsNodeNormal  =======================*/
static void pfsNodeNormal( int id )
{
 PfsNode    *node;               /* pointer to reference node */
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */
 PfsElem    *elem;               /* index of an adjacent element */
 double     size;                /* normal vector size */

 node = &pfs->nodelist[id];

 node->normal.x = 0.0;
 node->normal.y = 0.0;
 node->normal.z = 0.0;

 for( adjnode = node->adj_head; adjnode != NULL; adjnode = adjnode->next )
 {
  elem = adjnode->elem;
  if( elem != NULL )
  {
   node->normal.x += elem->normal.x;
   node->normal.y += elem->normal.y;
   node->normal.z += elem->normal.z;
  }
 }

 size = PFSGeoVecLen( &node->normal );
 if( size > 0.0 )
 {
  node->normal.x /= size;
  node->normal.y /= size;
  node->normal.z /= size;
 }
}

/*========================  pfsNodeNormals  =======================*/
static void pfsNodeNormals( void )
{
 int       i;

 if( pfs->nodelist != NULL )
 {
  for( i = 0; i < pfs->numnodes; i++ )
  {
   pfsNodeNormal( i );
  }
 }
}

/*========================  pfsFindRefElem  =======================*/
static void pfsFindRefElem( void )
{
 double   xmin, ymin, zmin;
 double   xmax, ymax, zmax;
 PfsElem  *elem;
 void     *e = NULL;
 double   area_avrg = 0.0;
 PFSGeoPoint normal_avrg;
 double   size;
 double   dotprod;
 double   max_dotprod = 0.0;

 pfs->ref_elem = NULL;
 pfs->ref_node = -1;

 if( pfs->elm3d_tree == NULL )
  return;

 normal_avrg.x = 0.0;
 normal_avrg.y = 0.0;
 normal_avrg.z = 0.0;

 PFSRtreeInitTraverse( pfs->elm3d_tree );
 while( (e = PFSRtreeTraverse( pfs->elm3d_tree,
                               &xmin, &xmax, 
                               &ymin, &ymax,
                               &zmin, &zmax )) != NULL )
 {
  elem = (PfsElem *)e;
  area_avrg += elem->area;
  normal_avrg.x += elem->area * elem->normal.x;
  normal_avrg.y += elem->area * elem->normal.y;
  normal_avrg.z += elem->area * elem->normal.z;
 }

 size = PFSGeoVecLen( &normal_avrg );

 if( size > 0.0 )
 {
  normal_avrg.x /= size;
  normal_avrg.y /= size;
  normal_avrg.z /= size;
 }
 else
 {
  return;
 }

 area_avrg /= (double)pfs->numelems;

 PFSRtreeInitTraverse( pfs->elm3d_tree );
 while( (e = PFSRtreeTraverse( pfs->elm3d_tree,
                               &xmin, &xmax, 
                               &ymin, &ymax,
                               &zmin, &zmax )) != NULL )
 {
  elem = (PfsElem *)e;
  if( elem->area > area_avrg*0.5 )
  {
   dotprod = PFSGeoDotProd( &elem->normal, &normal_avrg );
   if( dotprod > max_dotprod )
   {
    max_dotprod = dotprod;
    pfs->ref_elem = elem;
   }
  }
 }

 if( pfs->ref_elem != NULL )
 {
  pfs->ref_node = elem->inc[0];
  pfs->iter_toler = ITER_TOLER * Len3(
   pfs->nodelist[elem->inc[1]].coord.x - pfs->nodelist[pfs->ref_node].coord.x, 
   pfs->nodelist[elem->inc[1]].coord.y - pfs->nodelist[pfs->ref_node].coord.y, 
   pfs->nodelist[elem->inc[1]].coord.z - pfs->nodelist[pfs->ref_node].coord.z );
 }
}

/*=====================  pfsNodeConfMapCoeff  =====================*/
static void pfsNodeConfMapCoeff( int ref_id )
{
 PfsNode    *node;               /* pointer to reference node */
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */
 PfsElem    *e_prev;             /* CCW previous adjacent element */
 PfsElem    *e_next;             /* CCW next adjacent element */
 int        adj_id;              /* id of adjacent node */
 int        ref_previd = 0;      /* local id of ref. node in previous elem. */
 int        adj_previd = 0;      /* local id of adj. node in previous elem. */
 int        ref_nextid = 0;      /* local id of ref. node in next elem. */
 int        adj_nextid = 0;      /* local id of adj. node in next elem. */
 int        j;

 node = &pfs->nodelist[ref_id];

/* Reset the coefficient values.
 */
 node->AtA[0][0] = 0.0;
 node->AtA[0][1] = 0.0;
 node->AtA[1][0] = 0.0;
 node->AtA[1][1] = 0.0;

/* Get the last adjacent node in counter-clockwise order.
 * The adjacent element of the last adjacent node (if it exists) is 
 * the previous element of the first adjacent node.
 * It might be null if the reference node is a boundary node.
 */
 for( adjnode = node->adj_head; adjnode->next != NULL; adjnode = adjnode->next )
 ;
 e_next = adjnode->elem;

/* Traverse the adjacent nodes and elements in counter-clockwise order
 * around reference node and add the contributions of the previous
 * and next adjacent elements to the coefficients of the reference node
 * and of each adjacent node.
 */
 for( adjnode = node->adj_head; adjnode != NULL; adjnode = adjnode->next )
 {
  adjnode->AtA[0][0] = 0.0;
  adjnode->AtA[0][1] = 0.0;
  adjnode->AtA[1][0] = 0.0;
  adjnode->AtA[1][1] = 0.0;
  e_prev = e_next;
  e_next = adjnode->elem;
  adj_id = adjnode->id;

/* Find the position of reference and adjacent nodes in both
 * previous and next adjacent elements.
 */
  for( j = 0; j < 3; j++ )
  {
   if( e_prev != NULL )
   {
    if( e_prev->inc[j] == ref_id )
     ref_previd = j;
    if( e_prev->inc[j] == adj_id )
     adj_previd = j;
   }
   if( e_next != NULL )
   {
    if( e_next->inc[j] == ref_id )
     ref_nextid = j;
    if( e_next->inc[j] == adj_id )
     adj_nextid = j;
   }
  }

/* Add the contribution of the next element (if it exists) to the 
 * coefficients of the reference node.
 */
  if( e_next != NULL )
  {
   node->AtA[0][0] += e_next->W[ref_nextid][0] * e_next->W[ref_nextid][0] +
                      e_next->W[ref_nextid][1] * e_next->W[ref_nextid][1];
   node->AtA[0][1] += e_next->W[ref_nextid][0] * e_next->W[ref_nextid][1] -
                      e_next->W[ref_nextid][1] * e_next->W[ref_nextid][0];
   node->AtA[1][0] += e_next->W[ref_nextid][1] * e_next->W[ref_nextid][0] -
                      e_next->W[ref_nextid][0] * e_next->W[ref_nextid][1];
   node->AtA[1][1] += e_next->W[ref_nextid][1] * e_next->W[ref_nextid][1] +
                      e_next->W[ref_nextid][0] * e_next->W[ref_nextid][0];
  }

/* Add the contribution of the previous and next elements (if 
 * they exist) to the coefficients of the adjacent node.
 */
  if( e_prev != NULL )
  {
   adjnode->AtA[0][0] += e_prev->W[adj_previd][0] * e_prev->W[ref_previd][0] +
                         e_prev->W[adj_previd][1] * e_prev->W[ref_previd][1];
   adjnode->AtA[0][1] += e_prev->W[adj_previd][0] * e_prev->W[ref_previd][1] -
                         e_prev->W[adj_previd][1] * e_prev->W[ref_previd][0];
   adjnode->AtA[1][0] += e_prev->W[adj_previd][1] * e_prev->W[ref_previd][0] -
                         e_prev->W[adj_previd][0] * e_prev->W[ref_previd][1];
   adjnode->AtA[1][1] += e_prev->W[adj_previd][1] * e_prev->W[ref_previd][1] +
                         e_prev->W[adj_previd][0] * e_prev->W[ref_previd][0];
  }

  if( e_next != NULL )
  {
   adjnode->AtA[0][0] += e_next->W[adj_nextid][0] * e_next->W[ref_nextid][0] +
                         e_next->W[adj_nextid][1] * e_next->W[ref_nextid][1];
   adjnode->AtA[0][1] += e_next->W[adj_nextid][0] * e_next->W[ref_nextid][1] -
                         e_next->W[adj_nextid][1] * e_next->W[ref_nextid][0];
   adjnode->AtA[1][0] += e_next->W[adj_nextid][1] * e_next->W[ref_nextid][0] -
                         e_next->W[adj_nextid][0] * e_next->W[ref_nextid][1];
   adjnode->AtA[1][1] += e_next->W[adj_nextid][1] * e_next->W[ref_nextid][1] +
                         e_next->W[adj_nextid][0] * e_next->W[ref_nextid][0];
  }
 }
}

/*===================  pfsConfMapCoefficients  ====================*/
static void pfsConfMapCoefficients( void )
{
 int       i;

 if( pfs->ref_elem == NULL )
  return;

 if( pfs->nodelist != NULL )
 {
  for( i = 0; i < pfs->numnodes; i++ )
  {
   pfsNodeConfMapCoeff( i );
  }
 }
}

/*========================  pfsInitParams  ========================*/
static void pfsInitParams( void )
{
 PFSGeoNormPlane refplane;         /* plane that contains reference element */
 PfsNode      *refnode;         /* pointer to reference node */
 PfsNode      *node;            /* pointer to current node */
 int          nodeid;           /* index of a node */
 PFSGeoPoint     proj_pt;          /* proj. of node on plane */
 PFSGeoPoint     x_axis = { 1.0, 0.0, 0.0 };
 PFSGeoPoint     y_axis = { 0.0, 1.0, 0.0 };
 PFSGeoPoint     z_axis = { 0.0, 0.0, 1.0 };

 if( !pfs )
  return;

 refnode = &pfs->nodelist[pfs->ref_node];

/* Make the normal to the reference plane be the normal of the ref. elem.
 */
 refplane.n.x = pfs->ref_elem->normal.x;
 refplane.n.y = pfs->ref_elem->normal.y;
 refplane.n.z = pfs->ref_elem->normal.z;

/* Make the origin of parametric plane space be the first ref. point.
 */
 refplane.a.x = refnode->coord.x;
 refplane.a.y = refnode->coord.y;
 refplane.a.z = refnode->coord.z;

/* Try to lay the first plane parametric vector be at the projection
 * of the major x-axis with the plane.  If the x-axis is almost (less
 * then 30 degrees) parallel to the plane normal vector, try to do 
 * the same with the y-axis and z-axis (in this order).
 */
 PFSGeoCrossProd( &refplane.n, &x_axis, &refplane.v );
 if( ABS( PFSGeoVecLen( &refplane.v ) ) < 0.5 )
 {
  PFSGeoCrossProd( &refplane.n, &y_axis, &refplane.v );
  if( ABS( PFSGeoVecLen( &refplane.v ) ) < 0.5 )
  {
   PFSGeoCrossProd( &refplane.n, &z_axis, &refplane.v );
  }
 }
 PFSGeoVecNormalize( &refplane.v, &refplane.v );
 PFSGeoCrossProd( &refplane.v, &refplane.n, &refplane.u );

/* Compute the initial parametric coordinates of all nodes 
 * from the projection of their points on the reference plane.
 */
 for( nodeid = 0; nodeid < pfs->numnodes; nodeid++ )
 {
  node = &pfs->nodelist[nodeid];
  PFSGeoClstPtNormPlane( &refplane, &node->coord, &proj_pt, &node->par );
 }

/* To avoid round off errors, set the fixed null parametric values of
 * the reference point.
 */
 refnode->par.u = 0.0;
 refnode->par.v = 0.0;
}

/*==================  pfsCreateSolOrderLevel  ==================*/
static int pfsCreateNextSolOrderLevel( int lev_begin, int lev_end, int *pos,
                                       int *solorder, int *nodevisited )
{
 int        node_id;             /* id of reference node */
 PfsNode    *node;               /* pointer to reference node */
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */
 int        found_node = 0;
 int        i;

 for( i = lev_begin; i <= lev_end; i++ )
 {
  node_id = solorder[i];
  node = &pfs->nodelist[node_id];

  for( adjnode = node->adj_head; adjnode != NULL; adjnode = adjnode->next )
  {
   if( ! nodevisited[adjnode->id] )
   {
    found_node = 1;
    nodevisited[adjnode->id] = 1;
    solorder[*pos] = adjnode->id;
    (*pos)++;
   }
  }
 }

 return( found_node );
}

/*====================  pfsBuildSolOrderVec  ======================*/
static void pfsBuildSolOrderVec( void )
{
 int    *nodevisited;  /* array with flags for already visited nodes */
 int    lev_begin;     /* id of first node of current level in sol. order vec.*/
 int    lev_end;       /* id of last node of current level in sol. order vec.*/
 int    pos;           /* position of next node in solution order vector */
 int    i;

 if( pfs->sol_order != NULL )
  pfsDelSolOrderVec( pfs );

 pfs->sol_order = (int *)calloc( pfs->numnodes, sizeof( int ) );
 nodevisited = (int *)calloc( pfs->numnodes, sizeof( int ) );

 pfs->sol_order[0] = pfs->ref_node;
 for( i = 0; i < pfs->numnodes; i++ )
  nodevisited[i] = 0;
 nodevisited[pfs->ref_node] = 1;

 lev_begin = 0;
 lev_end = 0;
 pos = 1;
 while( pfsCreateNextSolOrderLevel( lev_begin, lev_end, &pos,
        pfs->sol_order, nodevisited ) )
 {
  lev_begin = lev_end + 1;
  lev_end = pos - 1;
 }

 free( nodevisited );
}

/*======================  pfsDelSolOrderVec  ======================*/
static void pfsDelSolOrderVec( PfsSurf *surf )
{
 if( surf->sol_order != NULL )
 {
  free( surf->sol_order );
  surf->sol_order = NULL;
 }
}

/*========================  pfsSolve2x2  ==========================*/
static void pfsSolve2x2( double a[2][2], double b[2], double *u, double *v )
{
 double det = a[0][0]*a[1][1] - a[0][1]*a[1][0];

 if( det > 0.0 )
 {
  *u = (a[1][1]*b[0] - a[0][1]*b[1]) / det;
  *v = (a[0][0]*b[1] - a[1][0]*b[0]) / det;
 }
 else
 {
  *u = 0.0;
  *v = 0.0;
 }
}

/*=====================  pfsSolveSurfPar1Pass  ====================*/
static int pfsSolveSurfPar1Pass( void )
{
 int        updated = 0;
 PfsNode    *node;
 PfsAdjNode *cur_adj;
 double     b[2];
 PFSGeoSurfPar newpar;
 int        i;

 for( i = 1; i < pfs->numnodes; i++ )
 {
  node = &pfs->nodelist[pfs->sol_order[i]];
  b[0] = b[1] = 0.0;
  newpar.u = newpar.v = 0.0;
  for( cur_adj = node->adj_head; cur_adj != NULL; cur_adj = cur_adj->next )
  {
   b[0] -= cur_adj->AtA[0][0] * pfs->nodelist[cur_adj->id].par.u + 
           cur_adj->AtA[0][1] * pfs->nodelist[cur_adj->id].par.v;
   b[1] -= cur_adj->AtA[1][0] * pfs->nodelist[cur_adj->id].par.u + 
           cur_adj->AtA[1][1] * pfs->nodelist[cur_adj->id].par.v;
   pfsSolve2x2( node->AtA, b, &newpar.u, &newpar.v );
  }
  if( ABS( newpar.u - node->par.u ) > pfs->iter_toler )
  {
   node->par.u = newpar.u;
   updated = 1;
  }
  if( ABS( newpar.v - node->par.v ) > pfs->iter_toler )
  {
   node->par.v = newpar.v;
   updated = 1;
  }
 }

 return( updated );
}

/*====================  pfsNormalizeParVals  ====================*/
static void pfsNormalizeParVals( void )
{
 double    u_size;                   /* current u size of parametric box */
 double    v_size;                   /* current u size of parametric box */
 double    scale;                    /* scaling factor */
 int       i;                        /* counter   */
 
 if( pfs->numnodes < 1 )
  return;

/* Compute scaling factor for normalization.
 */
 u_size = pfs->box_umax - pfs->box_umin;
 v_size = pfs->box_vmax - pfs->box_vmin;
 scale = MAX( u_size, v_size );
 scale = 1.0 / scale;

/* Now normalize parametric coordinates using same scale factor.
 */
 for( i = 0; i < pfs->numnodes; i++ )
 {
  pfs->nodelist[i].par.u = (pfs->nodelist[i].par.u - pfs->box_umin) * scale;
  pfs->nodelist[i].par.v = (pfs->nodelist[i].par.v - pfs->box_vmin) * scale;
 }

/* Finally update parametric bounding box.
 */
 pfs->box_umin = 0.0;
 pfs->box_umax = u_size * scale;
 pfs->box_vmin = 0.0;
 pfs->box_vmax = v_size * scale;
}

/*====================  pfsComputeParBoundBox  ====================*/
static void pfsComputeParBoundBox( void )
{
 int       i;                              /* counter   */

 if( pfs->numnodes < 1 )
  return;

 pfs->box_umin = pfs->box_umax = pfs->nodelist[0].par.u;
 pfs->box_vmin = pfs->box_vmax = pfs->nodelist[0].par.v;

 for( i = 1; i < pfs->numnodes; i++ )
 {
  if( pfs->box_umin > pfs->nodelist[i].par.u ) 
   pfs->box_umin = pfs->nodelist[i].par.u;
  if( pfs->box_umax < pfs->nodelist[i].par.u ) 
   pfs->box_umax = pfs->nodelist[i].par.u;
  if( pfs->box_vmin > pfs->nodelist[i].par.v ) 
   pfs->box_vmin = pfs->nodelist[i].par.v;
  if( pfs->box_vmax < pfs->nodelist[i].par.v ) 
   pfs->box_vmax = pfs->nodelist[i].par.v;
 }
}

/*=======================  pfsElemParBoundBox  ========================*/
static void pfsElemParBoundBox( PfsElem *elem )
{
 int    node;
 int    i;

 node = elem->inc[0];

 elem->minparbox.u =
 elem->maxparbox.u = pfs->nodelist[node].par.u;
 elem->minparbox.v =
 elem->maxparbox.v = pfs->nodelist[node].par.v;

 for( i = 1; i < 3; i++ )
 {
  node = elem->inc[i];

  if( pfs->nodelist[node].par.u < elem->minparbox.u )
      elem->minparbox.u = pfs->nodelist[node].par.u;

  if( pfs->nodelist[node].par.v < elem->minparbox.v )
      elem->minparbox.v = pfs->nodelist[node].par.v;

  if( pfs->nodelist[node].par.u > elem->maxparbox.u )
      elem->maxparbox.u = pfs->nodelist[node].par.u;

  if( pfs->nodelist[node].par.v > elem->maxparbox.v )
      elem->maxparbox.v = pfs->nodelist[node].par.v;
 }

 PFSRtreeInsert( pfs->elm2d_tree, (void *)elem, elem->minparbox.u,
                                                elem->maxparbox.u,
                                                elem->minparbox.v,
                                                elem->maxparbox.v, 0.0, 0.0 );
}

/*=================  pfsComputeElemParBoundBoxes  =====================*/
static void pfsComputeElemParBoundBoxes( void )
{
 double  xmin, ymin, zmin;
 double  xmax, ymax, zmax;
 PfsElem *elem;
 void    *e = NULL;

 if( pfs->elm3d_tree != NULL )
 {
  PFSRtreeInitTraverse( pfs->elm3d_tree );

  while( (e = PFSRtreeTraverse( pfs->elm3d_tree,
                                &xmin, &xmax, 
                                &ymin, &ymax,
                                &zmin, &zmax )) != NULL )
  {
   elem = (PfsElem *)e;
   pfsElemParBoundBox( elem );
  }
 }
}

/*=====================  pfsComputeElemJac  ======================*/
static void pfsComputeElemJac( PfsElem *elem, 
                               double *xdu, double *ydu, double *zdu, 
                               double *xdv, double *ydv, double *zdv )
{
 PFSGeoPoint   p1 = pfs->nodelist[elem->inc[0]].coord;
 PFSGeoPoint   p2 = pfs->nodelist[elem->inc[1]].coord;
 PFSGeoPoint   p3 = pfs->nodelist[elem->inc[2]].coord;
 PFSGeoSurfPar l1 = pfs->nodelist[elem->inc[0]].par;
 PFSGeoSurfPar l2 = pfs->nodelist[elem->inc[1]].par;
 PFSGeoSurfPar l3 = pfs->nodelist[elem->inc[2]].par;

 double area = ((l2.u - l1.u) * (l3.v - l1.v) -
                (l2.v - l1.v) * (l3.u - l1.u)) * 0.5;

 double fact = 0.5 / area;

 *xdu = fact * (p1.x * (l2.v - l3.v) + 
                p2.x * (l3.v - l1.v) +
                p3.x * (l1.v - l2.v));
 *ydu = fact * (p1.y * (l2.v - l3.v) + 
                p2.y * (l3.v - l1.v) +
                p3.y * (l1.v - l2.v));
 *zdu = fact * (p1.z * (l2.v - l3.v) + 
                p2.z * (l3.v - l1.v) +
                p3.z * (l1.v - l2.v));

 *xdv = fact * (p1.x * (l3.u - l2.u) + 
                p2.x * (l1.u - l3.u) +
                p3.x * (l2.u - l1.u));
 *ydv = fact * (p1.y * (l3.u - l2.u) + 
                p2.y * (l1.u - l3.u) +
                p3.y * (l2.u - l1.u));
 *zdv = fact * (p1.z * (l3.u - l2.u) + 
                p2.z * (l1.u - l3.u) +
                p3.z * (l2.u - l1.u));
}

/*=========================  pfsNodeJacobian  =======================*/
static void pfsNodeJacobian( int id )
{
 PfsElem    *elem;               /* index of an adjacent element */
 PfsNode    *node;               /* pointer to reference node */
 PfsAdjNode *adjnode;            /* pointer to current adjacent node */

 double xdu, ydu, zdu;
 double xdv, ydv, zdv;

 int nelem = 0;

 node = &pfs->nodelist[id];

 node->ddu.x = 0.0;
 node->ddu.y = 0.0;
 node->ddu.z = 0.0;

 node->ddv.x = 0.0;
 node->ddv.y = 0.0;
 node->ddv.z = 0.0;

 for( adjnode = node->adj_head; adjnode != NULL; adjnode = adjnode->next )
 {
  elem = adjnode->elem;

  if( elem != NULL )
  {
   pfsComputeElemJac( elem, &xdu, &ydu, &zdu, &xdv, &ydv, &zdv );

   node->ddu.x += xdu;
   node->ddu.y += ydu;
   node->ddu.z += zdu;

   node->ddv.x += xdv;
   node->ddv.y += ydv;
   node->ddv.z += zdv;

   nelem++;
  }
 }

 node->ddu.x /= (double)nelem;
 node->ddu.y /= (double)nelem;
 node->ddu.z /= (double)nelem;

 node->ddv.x /= (double)nelem;
 node->ddv.y /= (double)nelem;
 node->ddv.z /= (double)nelem;
}

/*========================  pfsNodeJacobians  =======================*/
static void pfsNodeJacobians( void )
{
 int i;

 if(pfs->nodelist != NULL)
 {
  for(i = 0; i < pfs->numnodes; i++)
  {
   pfsNodeJacobian(i);
  }
 }
}

/*============================  pfsArea3d  ==========================*/
static double pfsArea3d( double xa, double ya, double za, 
                         double xb, double yb, double zb,
                         double xc, double yc, double zc,
                         double xn, double yn, double zn )
{
 PFSGeoPoint vba, vca, vd, normal;
 double area;

 vba.x = xb - xa;
 vba.y = yb - ya;
 vba.z = zb - za; 
  
 vca.x = xc - xa; 
 vca.y = yc - ya;
 vca.z = zc - za; 
 
 normal.x = xn; 
 normal.y = yn;
 normal.z = zn;
 
 PFSGeoCrossProd(&vba, &vca, &vd);
 area = PFSGeoVecLen( &vd )/2;
 
 if( PFSGeoDotProd(&vd, &normal) < 0 )
  area *= -1;         
                                                 
 return( area );
}

/*============================  pfsArea2d  ==========================*/
static double pfsArea2d( double u0, double v0,
                         double u1, double v1,
                         double u2, double v2 )
{
 double u01, v01, u02, v02;

 u01 = u1 - u0;
 v01 = v1 - v0;
 u02 = u2 - u0;
 v02 = v2 - v0;

 return((u01*v02 - v01*u02) * 0.5);
}

/*========================  pfsDistPntCurElem  ======================*/
static double pfsDistPntCurElem( PFSGeoPoint *p, PFSGeoPoint *clst )
{
 PFSGeoPoint     p0;           /* coordinates of first element node */
 PFSGeoPoint     p1;           /* coordinates of second element node */
 PFSGeoNormPlane plane;        /* normalized plane that contains element */
 PFSGeoSurfPar   spar;         /* parametric coordinates of closest pt on plane */
 double       dist;         /* distance between given point and closest pt */
 
 pfsRItrNodeCoords( cur_elem->inc[0], &p0.x, &p0.y, &p0.z );
 pfsRItrNodeCoords( cur_elem->inc[1], &p1.x, &p1.y, &p1.z );

/* Make the normal to the plane be the normal of the current element.
 */
 plane.n = cur_elem->normal;

/* Make the origin of parametric plane space be the first node of element.
 */
 plane.a = p0;

/* Make the u axis of parametric plane space be from first to second element
 * node.  The v axis is computed by the cross-product between the normal
 * vector and the u axis.
 */
 PFSGeoDiffVec( &p1, &p0, &plane.u );
 PFSGeoVecNormalize( &plane.u, &plane.u );
 PFSGeoCrossProd( &plane.n, &plane.u, &plane.v );

/* Project given point on element plane to get the closest point.
 */
 PFSGeoClstPtNormPlane( &plane, p, clst, &spar );

/* Get distance between given point and closest point.
 */
 dist = PFSGeoDist( p, clst );

 return( dist );
}

/*========================  pfsParPntCurElemAttract  =======================*/
static int pfsParPntCurElemAttract( double *x, double *y, double *z,
                           double *u, double *v )
{
 double tol = COL_TOLER;
 double xn[3], yn[3], zn[3], un[3], vn[3];
 double b0, b1, b2, A;
 double b0c, b1c, b2c, Ac;
 int status_cartesiano = 1;
 int status_param = 1;

 pfsRItrNodeCoords( cur_elem->inc[0], &xn[0], &yn[0], &zn[0] );
 pfsRItrNodeCoords( cur_elem->inc[1], &xn[1], &yn[1], &zn[1] );
 pfsRItrNodeCoords( cur_elem->inc[2], &xn[2], &yn[2], &zn[2] );
 pfsRItrNodeParVals( cur_elem->inc[0], &un[0], &vn[0] );
 pfsRItrNodeParVals( cur_elem->inc[1], &un[1], &vn[1] );
 pfsRItrNodeParVals( cur_elem->inc[2], &un[2], &vn[2] );

 pfs3dGetAreaCoord( *x, *y, *z, xn, yn, zn, &b0c, &b1c, &b2c );

 Ac = b0c + b1c + b2c;

/* Check for point outside given element
 */
 Ac = b0c + b1c + b2c;
 if( ((Ac + tol) < 1.0) || ((Ac - tol) > 1.0) )
  status_cartesiano = 0;
 if( ((b0c + tol) < 0.0) || ((b1c + tol) < 0.0) || ((b2c + tol) < 0.0) )
  status_cartesiano = 0;

 *u = b0c*un[0] + b1c*un[1] + b2c*un[2];
 *v = b0c*vn[0] + b1c*vn[1] + b2c*vn[2];

 if( status_cartesiano )
 { 
  pfsParGetAreaCoord( *u, *v, un, vn, &b0, &b1, &b2 );
 
 /* Check for point outside given element
  */
  A = b0 + b1 + b2;
  if( ((A + tol) < 1.0) || ((A - tol) > 1.0) )
   status_param = 0;
  if( ((b0 + tol) < 0.0) || ((b1 + tol) < 0.0) || ((b2 + tol) < 0.0) )
   status_param = 0;
 
  if( status_param )
  {
   /* Se o ponto estiver dentro do triângulo no espaço cartesiano
    * e no espaço paramétrico, recalcula suas coordenadas para
    * evitar erros devido a tolerância
    */
   if( b0c < 0.0 ) b0c = 0.0;
   if( b1c < 0.0 ) b1c = 0.0;
   if( b2c < 0.0 ) b2c = 0.0;
 
   Ac = b0c + b1c + b2c;
   b0c /= Ac; 
   b1c /= Ac; 
   b2c /= Ac; 
   *u = b0c*un[0] + b1c*un[1] + b2c*un[2];
   *v = b0c*vn[0] + b1c*vn[1] + b2c*vn[2];
 
   *x = b0c*xn[0] + b1c*xn[1] + b2c*xn[2];
   *y = b0c*yn[0] + b1c*yn[1] + b2c*yn[2];
   *z = b0c*zn[0] + b1c*zn[1] + b2c*zn[2];
   return 1;
  }
 }
 return 0;
}   

/*========================  pfsParPntCurElem  =======================*/
static int pfsParPntCurElem ( double x, double y, double z,
                           double *u, double *v )
{
 double tol = COL_TOLER;
 double xn[3], yn[3], zn[3], un[3], vn[3];
 double b0, b1, b2, A;
 double b0c, b1c, b2c, Ac;
 int status_cartesiano = 1;
 int status_param = 1;

 pfsRItrNodeCoords( cur_elem->inc[0], &xn[0], &yn[0], &zn[0] );
 pfsRItrNodeCoords( cur_elem->inc[1], &xn[1], &yn[1], &zn[1] );
 pfsRItrNodeCoords( cur_elem->inc[2], &xn[2], &yn[2], &zn[2] );
 pfsRItrNodeParVals( cur_elem->inc[0], &un[0], &vn[0] );
 pfsRItrNodeParVals( cur_elem->inc[1], &un[1], &vn[1] );
 pfsRItrNodeParVals( cur_elem->inc[2], &un[2], &vn[2] );

 pfs3dGetAreaCoord( x, y, z, xn, yn, zn, &b0c, &b1c, &b2c );

 Ac = b0c + b1c + b2c;

/* Check for point outside given element
 */
 Ac = b0c + b1c + b2c;
 if( ((Ac + tol) < 1.0) || ((Ac - tol) > 1.0) )
  status_cartesiano = 0;
 if( ((b0c + tol) < 0.0) || ((b1c + tol) < 0.0) || ((b2c + tol) < 0.0) )
  status_cartesiano = 0;

 *u = b0c*un[0] + b1c*un[1] + b2c*un[2];
 *v = b0c*vn[0] + b1c*vn[1] + b2c*vn[2];
 
 if( status_cartesiano )
 {
  pfsParGetAreaCoord( *u, *v, un, vn, &b0, &b1, &b2 );
 
 /* Check for point outside given element
  */
  A = b0 + b1 + b2;
  if( ((A + tol) < 1.0) || ((A - tol) > 1.0) )
   status_param = 0;
  if( ((b0 + tol) < 0.0) || ((b1 + tol) < 0.0) || ((b2 + tol) < 0.0) )
   status_param = 0;

  if( status_param )
  {
   /* Point is inside triangle
  */
   return 1;
  }
 }
 return 0;
}   

/*=========================  pfs3dPntCurElem  =========================*/
static int pfs3dPntCurElem( double u, double v,
                           double *x, double *y, double *z )
{
 double un[3], vn[3];
 double xn[3], yn[3], zn[3];
 double b0, b1, b2, A;

 pfsRItrNodeCoords( cur_elem->inc[0], &xn[0], &yn[0], &zn[0] );
 pfsRItrNodeCoords( cur_elem->inc[1], &xn[1], &yn[1], &zn[1] );
 pfsRItrNodeCoords( cur_elem->inc[2], &xn[2], &yn[2], &zn[2] );
 pfsRItrNodeParVals( cur_elem->inc[0], &un[0], &vn[0] );
 pfsRItrNodeParVals( cur_elem->inc[1], &un[1], &vn[1] );
 pfsRItrNodeParVals( cur_elem->inc[2], &un[2], &vn[2] );

 pfsParGetAreaCoord( u, v, un, vn, &b0, &b1, &b2 );

 *x = b0*xn[0] + b1*xn[1] + b2*xn[2];
 *y = b0*yn[0] + b1*yn[1] + b2*yn[2];
 *z = b0*zn[0] + b1*zn[1] + b2*zn[2];

/* Check for point outside given element
 */
 A = b0 + b1 + b2;
 if((A + TOLER) < 1 || (A - TOLER) > 1)
  return( 0 );
 if((b0 + TOLER) < 0 || (b1 + TOLER) < 0 || (b2 + TOLER) < 0)
  return( 0 );

/* Point is inside triangle
 */
 return( 1 );
}

/*=======================  pfsParDerivCurElem  ========================*/
static int pfsParDerivCurElem( double u, double v,
                       double *xdu, double *ydu, double *zdu,
                               double *xdv, double *ydv, double *zdv,
                               double *xdw, double *ydw, double *zdw )
{
 double un[3], vn[3];
 double xdun[3], ydun[3], zdun[3];
 double xdvn[3], ydvn[3], zdvn[3];
 double xdwn[3], ydwn[3], zdwn[3];
 double b0, b1, b2, A;

 pfsRItrNodeParVals( cur_elem->inc[0], &un[0], &vn[0] );
 pfsRItrNodeParVals( cur_elem->inc[1], &un[1], &vn[1] );
 pfsRItrNodeParVals( cur_elem->inc[2], &un[2], &vn[2] );

 pfsRItrNodeJacobian( cur_elem->inc[0], &xdun[0], &ydun[0], &zdun[0],
                                       &xdvn[0], &ydvn[0], &zdvn[0]);
 pfsRItrNodeJacobian( cur_elem->inc[1], &xdun[1], &ydun[1], &zdun[1],
                                       &xdvn[1], &ydvn[1], &zdvn[1]);
 pfsRItrNodeJacobian( cur_elem->inc[2], &xdun[2], &ydun[2], &zdun[2],
                                       &xdvn[2], &ydvn[2], &zdvn[2]);

 pfsRItrNodeNormal( cur_elem->inc[0], &xdwn[0], &ydwn[0], &zdwn[0]);
 pfsRItrNodeNormal( cur_elem->inc[1], &xdwn[1], &ydwn[1], &zdwn[1]);
 pfsRItrNodeNormal( cur_elem->inc[2], &xdwn[2], &ydwn[2], &zdwn[2]);

 pfsParGetAreaCoord( u, v, un, vn, &b0, &b1, &b2 );

 *xdu = b0*xdun[0] + b1*xdun[1] + b2*xdun[2];
 *ydu = b0*ydun[0] + b1*ydun[1] + b2*ydun[2];
 *zdu = b0*zdun[0] + b1*zdun[1] + b2*zdun[2];

 *xdv = b0*xdvn[0] + b1*xdvn[1] + b2*xdvn[2];
 *ydv = b0*ydvn[0] + b1*ydvn[1] + b2*ydvn[2];
 *zdv = b0*zdvn[0] + b1*zdvn[1] + b2*zdvn[2];
  
 *xdw = b0*xdwn[0] + b1*xdwn[1] + b2*xdwn[2];
 *ydw = b0*ydwn[0] + b1*ydwn[1] + b2*ydwn[2];
 *zdw = b0*zdwn[0] + b1*zdwn[1] + b2*zdwn[2];

/* Check for point outside given element
 */
 A = b0 + b1 + b2;
 if((A + TOLER) < 1 || (A - TOLER) > 1)
  return( 0 );
 if((b0 + TOLER) < 0 || (b1 + TOLER) < 0 || (b2 + TOLER) < 0)
  return( 0 );

/* Point is inside triangle
 */
 return( 1 );
}

/*======================  pfsClstPtParLine  =======================*/
static int pfsClstPtParLine( double ua, double va, double ub, double vb,
                             double u, double v, double *up, double *vp,
                             double *dist )
{
 double u_perp, v_perp;   /* vector perpendicular to edge */
 double size;
 double perp_pa, perp_pb;

/* Get normalized vector perpendicular to edge.
 */
 u_perp = vb - va;
 v_perp = ua - ub;
 size = sqrt( u_perp*u_perp + v_perp*v_perp );
 if( size == 0.0 )
  return( 0 );
 u_perp /= size;
 v_perp /= size;

/* The perpendicular line passing at the given point WILL intercept the
 * line (ab) if the cross product (perp x pa) and (perp x pb) have 
 * different signs (or if the product of the cross products is negative).
 * If (perp x pa) or (perp x pb) is equal zero, the perpendicular line 
 * intersepts the node "a" or "b", respectively. 
 */
 perp_pa = ( u_perp * (va - v) ) - ( v_perp * (ua - u) );
 perp_pb = ( u_perp * (vb - v) ) - ( v_perp * (ub - u) );
 if( perp_pa == 0.0 )
 {
  *up = ua;
  *vp = va;
 }
 else if( perp_pb == 0.0 )
 { 
  *up = ub;
  *vp = vb;
 }
 else if( perp_pa * perp_pb < 0.0 )
 {
/* Compute the intersection point of the perpendicular line passing at  
 * (u,v) with line (ab).
 */
  *up = X_CROSS_LINE_SHOT( ua, va, ub, vb, u, v, u_perp, v_perp );
  *vp = Y_CROSS_LINE_SHOT( ua, va, ub, vb, u, v, u_perp, v_perp );
 }
 else /* if( perp_pa * perp_pb > 0.0 ) */
 {
  return( 0 );
 }

/* Compute the distance between the given point and the closest point
 * on edge.
 */
 *dist = Len2( u - (*up), v - (*vp) );

 return( 1 );
}

/*============================  pfsIsParPtInside  ==========================*/
static int pfsIsParPtInside( PFSGeoSurfPar *p )
{
 int ni = 0;
 int numloops;
 int loop;
 int node;

 numloops = pfsRItrNumLoops( );

 for( loop = 0; loop < numloops; loop++ )
 {
  int stop = 0;

  if( pfsRItrFirstLoopNode( loop, &node ) == PFS_SUCCESS )
  {
   int frst = node;
   int curr;
   int next;

   do
   {
    double p1u, p1v, p2u, p2v;
    curr = node;

    if( pfsRItrNextLoopNode( loop, curr, &node ) != PFS_SUCCESS )
     stop = 1;

    if(!stop)
     next = node;
    else
     next = frst;

    pfsRItrNodeParVals(curr, &p1u, &p1v);
    pfsRItrNodeParVals(next, &p2u, &p2v);

       /* discards horizontal edges */ 
    if(!(p1v == p2v) &&
       /* discards edges entirely above */
       !((p1v > p->v) && (p2v > p->v)) &&
       /* discards edges entirely below */
       !((p1v < p->v) && (p2v < p->v)) &&
       /* discards edges entirely to the left */
       !((p1u < p->u) && (p2u < p->u)))
    {
     /* if the initial point is on the same cote */
     if(p1v == p->v)
     {
      /* if the point is on the left and the edge is above it */
      if((p1u > p->u) && (p2v > p->v))
       ni++;
     }
     else
     {
      /* if the last point is on the same cote */
      if(p2v == p->v)
      {
       /* if the point is on the left and the edge is above it */
       if((p2u > p->u) && (p1v > p->v))
        ni++;
      }
      else
      {
       /* if the edge is entirely on the right */
       if((p1u > p->u) && (p2u > p->u))
        ni++;
       else
       {
        /* find the intersection point */
        double dx = p1u-p2u;
        double xc = p1u;

        /*if ( dx != 0 )*/
        if(fabs(dx) > TOLER)
         xc += (p->v-p1v)*dx / (p1v-p2v);

        /* if it is on the right */
        if(xc > p->u)
         ni++;
       }
      }
     }
    }
   } while(!stop);
  }
 }

 /* TRUE if the number of intersections is odd */
 return(ni%2 != 0);
}

/*=========================  pfsIsParPtOnBound  ======================*/
static int pfsIsParPtOnBound( PFSGeoSurfPar *p )
{
 int numloops;
 int loop;
 int node;
 double uc, vc;
 double un, vn;
 double up, vp;
 double dist;

 numloops = pfsRItrNumLoops( );

 for( loop = 0; loop < numloops; loop++ )
 {
  int stop = 0;

  if( pfsRItrFirstLoopNode( loop, &node ) == PFS_SUCCESS )
  {
   int frst = node;
   int curr;
   int next;

   do
   {
    curr = node;

    if( pfsRItrNextLoopNode( loop, curr, &node ) != PFS_SUCCESS )
     stop = 1;

    if(!stop)
     next = node;
    else
     next = frst;

    pfsRItrNodeParVals(curr, &uc, &vc);
    pfsRItrNodeParVals(next, &un, &vn);

    if( pfsClstPtParLine( uc, vc, un, vn, p->u, p->v, &up, &vp, &dist ) )
    {
     if( dist < TOLER )
    {
      p->u = up;
      p->v = vp;
      return( 1 );
     }
    }
   } while(!stop);
  }
 }
 return( 0 );
}

/*=========================  pfsSnapParPtToBound  ======================*/
static int pfsSnapParPtToBound( PFSGeoSurfPar *p )
{
 PFSGeoSurfPar clst;
 double dmin = 1e20;
 double dist;
 double ua, va, ub, vb;
 double up, vp;
 int found = 0;
 int numloops;
 int loop;
 int node;

 clst = *p;

 numloops = pfsRItrNumLoops( );

 for( loop = 0; loop < numloops; loop++ )
 {
  int stop = 0;

  if( pfsRItrFirstLoopNode( loop, &node ) == PFS_SUCCESS )
  {
   int frst = node;
   int curr;
   int next;

   do
   {
    curr = node;

    if( pfsRItrNextLoopNode( loop, curr, &node ) != PFS_SUCCESS )
     stop = 1;

    if(!stop)
     next = node;
    else
     next = frst;

    pfsRItrNodeParVals( curr, &ua, &va );
    pfsRItrNodeParVals( next, &ub, &vb );

    if( pfsClstPtParLine( ua, va, ub, vb, p->u, p->v, &up, &vp, &dist ) )
    {
     if( dist < dmin )
     {
      found = 1;
      clst.u = up;
      clst.v = vp;
      dmin = dist;
     }
    }
   } while(!stop);
  }
 }

 if( !found )
  return( 0 );

 *p = clst;

 return( 1 );
}

/*=======================  pfsChkProxElemNode  ========================*/
static int pfsChkProxElemNode( PfsElem *elem, double *u, double *v, int *ni )
{
 double un[3], vn[3];
 double b0, b1, b2;
 int    na;
 int    found = 1;

 pfsRItrNodeParVals( elem->inc[0], &un[0], &vn[0] );
 pfsRItrNodeParVals( elem->inc[1], &un[1], &vn[1] );
 pfsRItrNodeParVals( elem->inc[2], &un[2], &vn[2] );

 pfsParGetAreaCoord( *u, *v, un, vn, &b0, &b1, &b2 );

 /* Caso uma das áreas seja igual a um, considerando uma certa
  * tolerância, o pto está sobre um nó ou muito próximo a ele 
  */
 if( fabs( b0-1.0 ) < COL_TOLER )      na = 0;
 else if( fabs( b1-1.0 ) < COL_TOLER ) na = 1;
 else if( fabs( b2-1.0 ) < COL_TOLER ) na = 2;
 else found = 0;
   
 /* o pto é atraído para a posição do nó 
  */
 if( found )
 {
  *u = un[na];
  *v = vn[na];
  *ni = elem->inc[na];
  return( 1 );
 }
 return( 0 );
}

/*=========================  pfsChkProxElemEdge  ==========================*/
static int pfsChkProxElemEdge( PfsElem *elem, double *u, double *v, 
                               int *type, int *ni, int *nj )
{
 double un[3], vn[3];
 double b0, b1, b2;
 int    na, nb;
 double du, dv;
 double tp;
 int    found = 1;

 pfsRItrNodeParVals( elem->inc[0], &un[0], &vn[0] );
 pfsRItrNodeParVals( elem->inc[1], &un[1], &vn[1] );
 pfsRItrNodeParVals( elem->inc[2], &un[2], &vn[2] );

 pfsParGetAreaCoord( *u, *v, un, vn, &b0, &b1, &b2 );

 /* Caso uma das áreas seja igual a um, considerando uma certa
  * tolerância, o pto está sobre um nó ou muito próximo a ele 
  */
 if( fabs( b1 ) < COL_TOLER && fabs( b2 ) < COL_TOLER )      na = 0;
 else if( fabs( b0 ) < COL_TOLER && fabs( b2 ) < COL_TOLER ) na = 1;
 else if( fabs( b0 ) < COL_TOLER && fabs( b1 ) < COL_TOLER ) na = 2;
 else found = 0;
   
 /* o pto é atraído para a posição do nó 
  */
 if( found )
 {
  *u = un[na];
  *v = vn[na];
  *ni = elem->inc[na];
  *type = NODE;
  return( 1 );
 }

 /* Caso uma das áreas seja praticamente nula, o pto está sobre uma
  * aresta ou muito próximo a ela 
  */
 if( fabs( b0 ) < COL_TOLER && 
     fabs( b1 ) > COL_TOLER && fabs( b2) > COL_TOLER ) 
 {
  na = 1;
  nb = 2;
 }
 else if( fabs( b1 ) < COL_TOLER && 
          fabs( b0 ) > COL_TOLER && fabs( b2 ) > COL_TOLER ) 
 { 
  na = 0;
  nb = 2;
 }
 else if( fabs( b2 ) < COL_TOLER &&
          fabs( b0 ) > COL_TOLER && fabs( b1 ) > COL_TOLER )
 { 
  na = 0;
  nb = 1;
 }
 else return( 0 ); 

 *ni = elem->inc[na];
 *nj = elem->inc[nb];
 *type = EDGE;

 /* Atrai o pto para a posição da projeção do pto na aresta. 
  */
 du = un[nb] - un[na];
 dv = vn[nb] - vn[na];
 if( ( fabs(du) < TOLER ) && ( fabs(dv) < TOLER ) ) 
  return( 0 );

 /* Faz uma verificação para ver se a  projeção do ponto está
  * realmente sobre a aresta*/
 tp = (du*(*u - un[na]) + dv*(*v - vn[na])) / (du*du + dv*dv);
 if( tp <= 0.0 )
 {
  tp = 0.0;
  *ni = elem->inc[na];
  *type = NODE;
 } 
 else if( tp >= 1.0 )
 {
  tp = 1.0;
  *ni = elem->inc[nb];
  *type = NODE;
 }
 *u = un[na] + tp*du;
 *v = vn[na] + tp*dv;
 
 return( 1 );
}

/*========================  pfsShotXLine  =========================*/
static int pfsShotXLine( PfsElem *elem, double u0, double v0,
                         double u_shot, double v_shot,
                        int na, int nb, double *u1, double *v1,
                        int *type, int *n_inter )
{
 double size;
 double shot_0a, shot_0b;
 double ua, va, ub, vb;
 int    aux_type;

 pfsRItrNodeParVals( na, &ua, &va );
 pfsRItrNodeParVals( nb, &ub, &vb );

/* Normalize shot vector.
 */
 size = sqrt( u_shot*u_shot + v_shot*v_shot );
 if( size == 0.0 )
  return( 0 );
 u_shot /= size;
 v_shot /= size;

/* The shot vector WILL intercept line (ab) if the cross product 
 * (shot x 0a) and (shot x 0b) have diferent signs (or if the product
 * of the cross products is negative). If (shot x 0a) or (shot x 0b) 
 * is equal zero, the shot intersepts the node "a" or "b", 
 * respectively. 
 */
 shot_0a = ( u_shot * (va - v0) ) - ( v_shot * (ua - u0) );
 shot_0b = ( u_shot * (vb - v0) ) - ( v_shot * (ub - u0) );
 if( shot_0a == 0.0 )
 {
  *u1 = ua;
  *v1 = va;
  aux_type = NODE;
  *n_inter = na;
 }
 else if( shot_0b == 0.0 )
 { 
  *u1 = ub;
  *v1 = vb;
  aux_type = NODE;
  *n_inter = nb;
 }
 else if( shot_0a * shot_0b < 0.0 )
 {
/* Compute the intersection point of the shot line passing at  
 * (u0,v0) with line (ab).
 */
 *u1 = X_CROSS_LINE_SHOT( ua, va, ub, vb, u0, v0, u_shot, v_shot );
 *v1 = Y_CROSS_LINE_SHOT( ua, va, ub, vb, u0, v0, u_shot, v_shot );
 aux_type = EDGE;
 }
 else if( shot_0a * shot_0b > 0.0 )
 {
  return( 0 );
 }

/* Check to see if the intersection point, as seem from the starting 
 * point, is in the same direction of the shot vector.  This is known
 * by checking the sign of the dot product between a vector from the
 * starting point to the inserction point and the shot vector.
 */
 if( (((*u1 - u0) * u_shot) + ((*v1 - v0) * v_shot) ) <= 0.0 )
 {
  return( 0 );
 }

 if( ( aux_type == EDGE ) && ( pfsChkProxElemNode( elem, u1, v1, n_inter ) ) )
 {
  *type = NODE;
 }
 else
 {
  *type = aux_type;
 }
 return( 1 );
}

/*=========================  pfsNodeShot  =========================*/
static int pfsNodeShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1, double *v1, 
                        int *type, int *ni, int *nj )
{
 PfsNode    *node = &pfs->nodelist[*ni];
 PfsAdjNode *adjnode; /* ponteiro para os nós adjacentes */
 PfsElem    *elem;       /* ponteiro para os triângulos adjacentes */
 int        na, nb;
 int        i;
 int        n_inter;

 /* Verifica se o tiro intercepta as arestas opostas ao nó ou se
  * intercepta os nós das arestas.
  */
 for( adjnode = node->adj_head;
      ( adjnode!= NULL ) && ( adjnode->elem != NULL );
      adjnode = adjnode->next )
 {
  elem = adjnode->elem;

  if( ( elem != NULL ) && ( elem == init_elem ) )
   continue;

  /* Pego os nós das arestas opostas */
  na = adjnode->id;
 
  for( i=0; 
       i<3 && ( ( elem->inc[i] == *ni ) || ( elem->inc[i] == na ) );
       i++ ) ;
  nb = elem->inc[i];

  if( pfsShotXLine( elem, u0, v0, u_shot, v_shot, na, nb,
                    u1, v1, type, &n_inter ) )
  {
   *end_elem = elem;
   if(*type == NODE ) 
   {
    *ni = n_inter;
   }
   else if( *type == EDGE ) 
   {
    *ni = na;
    *nj = nb;
   }
    return( 1 ); 
   }
 }

 /* caso o ponto esteja inicialmente sobre um nó da fronteira
  * da malha, não se conheça o triangulo de origem (init_elem) 
  * e não ocorra interseção com os triângulos adjacentes ao nó é 
  * porque o tiro está na direção da face externa. Sendo assim, 
  * toma-se como triangulo final o primeiro triangulo adjacente ao nó.  
  *             ______ _______
  *             \    / \    / \
  *              \  /   \  /  <--- end_elem
  *               \/_____\/_____\  <--- fronteira
  *                      |
  *                      |
  *                     \|/ tiro
  */
  

 if( ( adjnode != NULL ) && ( adjnode->elem == NULL ) )
 {
  if(init_elem == NULL ) 
   *end_elem =  node->adj_head->elem;
  else
   *end_elem = init_elem;
  *type |= LOOP;
  return( 1 ); 
 }
 printf("\n%d\tPfsNodeShot - Nenhuma aresta ou no foram interceptados\n", pfs_node );
 return( 0 );
}

/*=========================  pfsEdgeShot  =========================*/
static int pfsEdgeShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1, double *v1, 
                        int *type, int *ni, int *nj )
{
 PfsElem *elem;
 PfsElem *e_left, *e_right;
 int     na, nb;  
 int     boundary = 0;
 int     first = 1;
 int     i;
 int     n_inter;  

 /* pega os triangulos adjacentes a aresta */
 pfsGetEdgeAdjElems( *ni, *nj, &e_left, &e_right );

 /* verifica se a aresta está na fronteira da malha */
 if( ( e_left == NULL ) || ( e_right == NULL ) )
  boundary = 1;

 /* A busca é feita no triangulo a esquerda e em seguida no 
  * triangulo a direita. Caso a aresta esteja na fronteira da malha
  * ou seja conhecido o triangulo de onde o tiro se originou 
  * (init_elem não nulo) a busca é feita em apenas um dos 
  * triângulos. 
  */
 if( ( e_left != NULL ) && ( e_left != init_elem ) )
 {
  elem = e_left;
 }
 else
 {
  elem = e_right;
 } 

 if( elem )
 {
  /* verifica se ocorre interseção do tiro com as arestas dos 
   * triângulos adjacentes com exceção da aresta sobre a qual o pto 
   * se encontra.
   */
  do{
   na = elem->inc[0]; 
   for( i=0; i<3; i++ )
   {
    nb = elem->inc[(i+1)%3];
  
    /* se for a aresta sobre a qual o pto está, não faz o teste */
    if( ( ( na == *ni ) && (nb == *nj ) ) ||
        ( ( na == *nj ) && (nb == *ni ) ) )
    {
     na = nb;
     continue;
    }
     
    if( pfsShotXLine( elem, u0, v0, u_shot, v_shot, na, nb,
                      u1, v1, type, &n_inter ) )
    {
     *end_elem = elem;
     if(*type == NODE ) 
     {
      *ni = n_inter;
     }
     else if( *type == EDGE ) 
     {
      *ni = na;
      *nj = nb;
     }
     return( 1 ); 
    }
    na = nb;
   }
   elem = ( elem == e_left ) ? e_right : NULL;
  }while( !boundary && ( elem == e_right )  && ( elem != init_elem ) );
 }

 /* caso o ponto esteja inicialmente sobre uma aresta da fronteira
  * da malha, não se conheça o triangulo de origem (init_elem) 
  * e não ocorra interseção com o triângulo adjacente à aresta é 
  * porque o tiro está na direção da face externa e o triangulo 
  * origem é o único triângulo adjacente. Sendo assim, o triangulo
  * final será igual ao único adjacente a aresta.  
  *             ______ _______
  *             \    / \    / \
  *              \  /   \  /   \
  *               \/_____\/_____\  <--- fronteira
  *                   |
  *                   |
  *                  \|/ tiro
  */
 if( boundary )
 {
  if(init_elem == NULL ) 
   *end_elem = ( e_left != NULL ) ? e_left : e_right;
  else
   *end_elem = init_elem;

  *u1 = u0;
  *v1 = v0;
  *type |= LOOP;
  return( 1 ); 
 }
 printf("\n%d PfsEdgeShot - Nenhuma aresta ou no foram interceptados\n", pfs_node );
 return( 0 ); 
}

/*=========================  pfsFaceShot  =========================*/
static int pfsFaceShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1, double *v1, 
                        int *type, int *ni, int *nj )
{
 int     na, nb;
 int     i;
 int     n_inter;

 /* Verifica se o tiro intercepta as arestas 
  */
 for( i=0; i<3; i++ )
 {
  /* Pego os nós das arestas opostas */
  na = init_elem->inc[i];
  nb = init_elem->inc[(i+1)%3];
 
  if( pfsShotXLine( init_elem, u0, v0, u_shot, v_shot, na, nb,
                    u1, v1, type, &n_inter ) )
  {
   *end_elem = init_elem;
   if(*type == NODE ) 
   {
    *ni = n_inter;
   }
   else if( *type == EDGE ) 
   {
    *ni = na;
    *nj = nb;
   }
   return( 1 ); 
  }
 }
 printf("\n%d\tPfsFaceShot - Nenhuma aresta ou nó foram interceptados\n", pfs_node );
 return( 0 );
}

/*===========================  pfsLoopShot ========================*/
static int pfsLoopShot( PfsElem *init_elem, double u0, double v0, 
                        double u_shot, double v_shot,
                        PfsElem **end_elem, double *u1,
                        double *v1, int *type, int *ni, int *nj )
{
 int           prev_id, cur_id;
 double        cur_len, len; 
 int           found = 0;
 int           aux_type, cur_type;
 PFSGeoPoint   p0, p_cur;
 double        u_inter, v_inter;
 int           ni_ref, nj_ref;
 int           n_inter;
 int           first;
 PfsNode       *prev_node;    
 PfsElem       *elem;

 ni_ref = *ni;
 nj_ref = *nj;

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1)
            || (pfs->numloops < 1) )
  return( 0 );

 if( ( pfs->nodelist[*ni].type != NODE_BOUNDARY ) || 
     ( ( *type & EDGE ) && 
       ( pfs->nodelist[*nj].type != NODE_BOUNDARY ) ) )
  return( 0 );

/* Calcula coordenada cartesiana do pto de partida do tiro */
 if( *type & NODE ) 
  pfsRItrNodeCoords( ni_ref, &p0.x, &p0.y, &p0.z );
 else
  pfs3dPntCurElem( u0, v0, &p0.x, &p0.y, &p0.z );

/* Percorre o loop por onde o tiro saiu e verifico se o tiro está 
 * interceptando outra aresta do loop.
 * Quando o tipo do ponto de interseção no loop for NODE,  
 * não verifica-se a interseção do tiro com as arestas do loop
 * adjacentes ao nó ni.
 * Qdo o tipo é EDGE, não testa a interseção do tiro com a aresta  
 * ninj. Qdo o tipo é FREE_FACE, verifica-se a interseção com todas
 * as arestas do loop. 
 */
 if( *type & NODE ) 
  prev_id = pfs->nodelist[ni_ref].adj_head->id; 
 else if( *type & EDGE ) 
  prev_id = nj_ref;
 else if( *type & FREE_FACE )
 {
  prev_id = ni_ref;
  first = 1; /*ao invés do for, o loop seria "do ... while" */
 }

 for( prev_node = &pfs->nodelist[prev_id], 
      cur_id = prev_node->adj_head->id; 
      ( ( *type & NODE ) && ( cur_id != ni_ref ) ) ||
      ( ( *type & EDGE ) && ( prev_id != ni_ref ) ) ||
      ( ( *type & FREE_FACE ) && 
        ( first || ( prev_id != ni_ref ) ) ); 
      prev_node = &pfs->nodelist[prev_id], 
      cur_id = prev_node->adj_head->id )
 {
  first = 0;
  elem = prev_node->adj_head->elem;
  if( pfsShotXLine( elem, u0, v0, u_shot, v_shot, 
                    prev_id, cur_id,  
                    &u_inter, &v_inter, &cur_type, &n_inter ) )
  {
   pfsRSetCurrElem( elem );
/* Coordenada cartesiana do pto do loop interceptado pelo tiro 
 */
   pfs3dPntCurElem( u_inter, v_inter, &p_cur.x, &p_cur.y, &p_cur.z );

/* Distância entre o nó interceptado e o pto de partida do tiro 
 */
   cur_len = PFSGeoDist( &p_cur, &p0 );
   
/* Mantém o pto de interseção mais próximo ao pto de partida 
 */ 
   if( !found || ( cur_len < len ) )
   {
    if( cur_type == NODE ) 
    {
     *ni = n_inter;
    }
    else if( cur_type == EDGE ) 
    {
     *ni = prev_id;
     *nj = cur_id;
    }
    len = cur_len;
    *u1 = u_inter;
    *v1 = v_inter;
    aux_type = cur_type;
    found = 1;
   }
  }
  prev_id = cur_id;
 }
 *end_elem = init_elem;
 
 if( found )
  *type = aux_type;
 else
  *type = FREE_FACE;
 
 return( 1 ); 
}

/*======================  pfsUpdateShotLen  =======================*/
static int pfsUpdateShotLen( double u0, double v0, PFSGeoPoint *p0, 
                             double *u1, double *v1, PFSGeoPoint *p1, 
                             double *len )
{
 double cur_len;
 PFSGeoPoint shot;

/* consulta a coordenada cartesiana da posição final do tiro 
 */
 pfs3dPntCurElem( *u1, *v1, &p1->x, &p1->y, &p1->z );

/* calcula a distância entre a posição final e inicial do tiro 
 */
 PFSGeoDiffVec( p1, p0, &shot );
 cur_len = PFSGeoVecLen( &shot );
  
/* caso o tiro tenha ultrapassado o comprimento necessário, 
 * recalcula a nova posiçao
 */
 if( cur_len > *len )
 {
#if 1
  double ratio = *len/cur_len;

  *u1 = u0 + ( *u1 - u0 ) * ratio;
  *v1 = v0 + ( *v1 - v0 ) * ratio;

 /* consulta a coordenada cartesiana da posição final do tiro 
  */
  if( !pfs3dPntCurElem( *u1, *v1, &p1->x, &p1->y, &p1->z ) )
  {
   if( !_loop_ ) 
    printf( "\n%d pfsUpdateShotLen\n", pfs_node );
  }
  return( 1 );
#else
 /* normaliza o vetor que representa a direção do tiro no espaço
  * cartesiano e escala para ele ter apenas o tamanho que falta
  * percorrer.
  */
  PFSGeoVecNormalize( &shot , &shot );
  PFSGeoProdScalVec( &shot, len, &shot );

 /* calcula a posição final do tiro no espaço cartesiano 
  */
  PFSGeoSumVec( p0, &shot, p1 );

 /* calcula a posição final do tiro no espaço paramétrico 
  */
  pfsParPntCurElemAttract( &p1->x, &p1->y, &p1->z, u1, v1 );
  *len = 0.0;
  return( 1 );
#endif
 }
 else
 {
 /* calcula o comprimento q ainda falta percorrer 
  */
  *len -= cur_len; 
  if( *len < TOLER )
   return( 1 );
  else
   return( 0 );
 }
}

/*==================  pfsFinishFreeFaceShotLen  ===================*/
static void pfsFinishFreeFaceShotLen( double u0, double v0, PFSGeoPoint *p0, 
                                      double u_shot, double v_shot, double len,
                                      double *u1, double *v1, PFSGeoPoint *p1 )  
{
 PFSGeoPoint shot;

/* calcula a coordenada paramétrica de um pto qualquer na direção do
 * tiro. 
 */
 *u1 = u0+u_shot; 
 *v1 = v0+v_shot; 

/* calculo a coord cartesiana desse ponto 
 */
 pfs3dPntCurElem( *u1, *v1, &p1->x, &p1->y, &p1->z ); 
 
/* calculo o vetor do tiro normalizado no espaço cartesiano e 
 * escala pro comprimento q ainda falta 
 */
 PFSGeoDiffVec( p1, p0, &shot );
 PFSGeoVecNormalize( &shot , &shot );
 PFSGeoProdScalVec( &shot, &len, &shot );
 
/* calcula a posição final do tiro no espaço cartesiano 
*/
 PFSGeoSumVec( p0, &shot, p1 ); 

/* calcula a posição final do tiro no espaço paramétrico 
 */
 pfsParPntCurElemAttract( &p1->x, &p1->y, &p1->z, u1, v1 );
}

 /*========================  pfsShotSolve  =========================*/
 static int pfsShotSolve( PfsElem *init_elem, double u0, double v0, 
                          PFSGeoPoint p0, double u_shot, double v_shot,
                          PfsElem **end_elem, double *u1, double *v1,
                          PFSGeoPoint *p1, int *type, int *ni, int *nj,
                          double *len )
 {
  int status = 0;
  int loop_in = 0; /* indica se o tiro está entrando na malha */
  _loop_ = 0;

  switch( *type )
  {
   case NODE:
    status = pfsNodeShot( init_elem, u0, v0, u_shot, v_shot, 
                          end_elem, u1, v1, type, ni, nj );
    break;
   case EDGE:
    status = pfsEdgeShot( init_elem, u0, v0, u_shot, v_shot, 
                          end_elem, u1, v1, type, ni, nj );
    break;
   case FACE:
    status = pfsFaceShot( init_elem, u0, v0, u_shot, v_shot, 
                          end_elem, u1, v1, type, ni, nj );
    break;
  default:
   if( *type & LOOP )
   {
    _loop_ = 1;
    status = pfsLoopShot( init_elem, u0, v0, u_shot, v_shot,
                          end_elem, u1, v1, type, ni, nj );
    loop_in = 1;
   }
   break;
 }
 
 /* Verifica se ocorreu algum erro e, caso não tenha ocorrido,
  * atualiza o elemento corrente. 
  */
 if( !status )
  return( 0 );
 else if( *end_elem != NULL )
  pfsRSetCurrElem( *end_elem );
  
  if( init_elem == NULL )
   pfs3dPntCurElem( u0, v0, &p0.x, &p0.y, &p0.z );
  
 /* Se não encontrou interseção é pque o tiro sai da malha 
  */
 if( *type == FREE_FACE )
 {
  pfsFinishFreeFaceShotLen( u0, v0, &p0, u_shot, v_shot, *len, u1, v1, p1 );  
 }
 else if(  ( *type & LOOP ) || 
           !pfsUpdateShotLen( u0, v0, &p0, u1, v1, p1, len ) )
 {
  if( *type & LOOP ) 
  {
   p1->x = p0.x;
   p1->y = p0.y;
   p1->z = p0.z;
  } 
  else if( loop_in ) 
  {
/* Qdo o tiro atravessa um loop, ao entrar novamente na malha
 * o tipo será sempre nó ou aresta e a busca deve ser feita em
 * todos os elementos adjacentes. Para isso, o elemento inicial
 * do próximimo passo tem q ser nulo
 */
   *end_elem = NULL;
  }

  pfsShotSolve( *end_elem, *u1, *v1, *p1,
                u_shot, v_shot, end_elem, u1, v1, p1, type, 
                ni, nj, len );  
 }
 return( 1 );
}

/*
** ---------------------------------------------------------------
** Entry points begin here:
*/

/*=========================  pfsRInitSurf  =========================*/

void *pfsRInitSurf( int mxnds, int mxels )
{
 PfsSurf *newsurf;

/* Allocate memory for a new surface.
 */
 newsurf = (PfsSurf *)malloc( sizeof( PfsSurf ) );
 newsurf->box_xmin = newsurf->box_xmin = 0.0;
 newsurf->box_ymin = newsurf->box_ymin = 0.0;
 newsurf->box_zmin = newsurf->box_zmin = 0.0;
 newsurf->box_umin = newsurf->box_umin = 0.0;
 newsurf->box_vmin = newsurf->box_vmin = 0.0;

/* Allocate memory for node list vector.
 */
 newsurf->maxnodes = mxnds;
 newsurf->nodelist = (PfsNode *)calloc( newsurf->maxnodes, sizeof( PfsNode ) );

/* Create two R-trees: one indexed by 3D element bounding box and
 * the other indexed by parametric element bounding box.
 */
 newsurf->elm2d_tree = PFSRtreeCreate( );
 newsurf->elm3d_tree = PFSRtreeCreate( );

/* Reset node solution order vector.
 */
 newsurf->sol_order = NULL;

/* Initialize total number of nodes and elements.
 */
 newsurf->numnodes = 0;
 newsurf->numelems = 0;

/* Initialize loops data.
 */
 newsurf->numloops = 0;
 newsurf->looplist = NULL;

/* Initialize reference element and reference node.
 */
 newsurf->ref_elem = NULL;
 newsurf->ref_node = -1;

 return( newsurf );
}

/*=======================  pfsRActivateSurf  =======================*/

void pfsRActivateSurf( void *surf )
{
 if( !surf )
  return;

 pfs = (PfsSurf *)surf;
 cur_elem = NULL;
 curr_ielem = NULL;
}

/*=========================  pfsRFreeSurf  =========================*/

void pfsRFreeSurf( void *surf )
{
 PfsSurf *delsurf;

 if( !surf )
  return;

 delsurf = (PfsSurf *)surf;

 pfsDelNodeList( delsurf );
 pfsDelElemList( delsurf );
 pfsDelSolOrderVec( delsurf );

 cur_elem = NULL;
 curr_ielem = NULL;

 free( delsurf );
}

/*==========================  pfsRAddNode  =========================*/

int pfsRAddNode( double x, double y, double z )
{
 int    node_id;

 if( !pfs )
  return( PFS_NO_SURFACE_DEFINED );

 if( pfs->numnodes == pfs->maxnodes )
 {
  pfs->maxnodes *= 2;
  pfs->nodelist = (PfsNode *)realloc( pfs->nodelist, 
                                      pfs->maxnodes*sizeof( PfsNode ) );
 }

 node_id = pfs->numnodes;
 pfs->nodelist[node_id].coord.x = x;
 pfs->nodelist[node_id].coord.y = y;
 pfs->nodelist[node_id].coord.z = z;
 pfs->nodelist[node_id].par.u = 0.0;
 pfs->nodelist[node_id].par.v = 0.0;
 pfs->nodelist[node_id].AtA[0][0] = 0.0;
 pfs->nodelist[node_id].AtA[0][1] = 0.0;
 pfs->nodelist[node_id].AtA[1][0] = 0.0;
 pfs->nodelist[node_id].AtA[1][1] = 0.0;
 pfs->nodelist[node_id].degree = 0;
 pfs->nodelist[node_id].adj_head = NULL;
 pfs->nodelist[node_id].type = NODE_INTERIOR;

 pfs->numnodes++;

 return( node_id );
}

/*==========================  pfsRAddElem  =========================*/

int pfsRAddElem( int v0, int v1, int v2 )
{
 if( !pfs )
  return 0;

 curr_ielem = (PfsElem *)calloc( 1, sizeof( PfsElem ) );

 curr_ielem->inc[0] = v0;
 curr_ielem->inc[1] = v1;
 curr_ielem->inc[2] = v2;

 pfsElemNormal( curr_ielem );
 pfsElemCenter( curr_ielem );
 pfsElemLocalGradients( curr_ielem );
 pfsElem3dBoundBox( curr_ielem );

 curr_ielem->minparbox.u = 0.0;
 curr_ielem->minparbox.v = 0.0;
 curr_ielem->maxparbox.u = 0.0;
 curr_ielem->maxparbox.v = 0.0;

 pfs->numelems++;

 return 1;
}

/*=======================  pfsRCompleteSurf  =======================*/

int pfsRCompleteSurf( void )
{
 int  status;

 if( !pfs )
  return( PFS_NO_SURFACE_DEFINED );

/* Build node adjacency data structure.
 */
 status = pfsBuildNodeAdj( );
 if( status != PFS_SUCCESS )
 {
  pfsRFreeSurf( pfs );
  pfs = NULL;
  return( status );
 }
 pfsClassifyNodes( );

/* Order the adjacent node list of all nodes such that the first 
 * adjacent node in list is the next boundary node when 
 * traversing the surface boundary.
 */
 status = pfsOrderAllBdryNodeAdj( );
 if( status != PFS_SUCCESS )
 {
  pfsRFreeSurf( pfs );
  pfs = NULL;
  return( status );
 }

/* Create loop list vector that contains the id of the first
 * node of each loop.
 */
 pfsBuildLoopList( );

/* Compute mesh 3D bounding box.
 */
 pfsCompute3dBoundBox( );

/* Compute node normal vectors.
 */
 pfsNodeNormals( );

/* Find reference element (which contains the origin of parametric space).
 */
 pfsFindRefElem( );

/* Computed conformal map coefficients for all nodes in surface.
 */
 pfsConfMapCoefficients( );

/* Initialize nodal parametric values.
 */
 pfsInitParams( );

/* Compute mesh parametric bounding box.
 */
 pfsComputeParBoundBox( );

/* Allocated and build node solution order vector.
 */
 pfsBuildSolOrderVec( );

 return( PFS_SUCCESS );
}

/*========================  pfsRSolveSurfPar  ======================*/

void pfsRSolveSurfPar( int maxnumiter )
{
 int        numiter = 0;
 int        not_done = 1;

 if( !pfs )
  return;

 if( maxnumiter == 0 )
 {
  while( not_done )
  {
   if( ! pfsSolveSurfPar1Pass( ) )
    not_done = 0;
  }
 }
 else
 {
  while( not_done && (numiter < maxnumiter) )
  {
   numiter++;
   if( ! pfsSolveSurfPar1Pass( ) )
    not_done = 0;
  }
 }

/* Compute mesh parametric bounding box.
 */
 pfsComputeParBoundBox( );
}

/*======================  pfsCompleteParSpace  ====================*/

void pfsCompleteParSpace( void )
{
/* Normalized nodal parametric coordinates.
 */
 pfsNormalizeParVals( );

/* Compute element parametric bounding boxes and store the 
 * elements in the R-tree indexed by parametric bounding boxes.
 */
 pfsComputeElemParBoundBoxes( );

/* Compute avarage Jacobian (mapping derivatves) for all nodes.
 */
 pfsNodeJacobians( );
}

/*=======================  pfsRItr3dBoundBox  ======================*/

int pfsRItr3dBoundBox( double *xmin, double *xmax,
                      double *ymin, double *ymax,
                      double *zmin, double *zmax )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return( 0 );

 if( xmin ) *xmin = pfs->box_xmin;
 if( xmax ) *xmax = pfs->box_xmax;
 if( ymin ) *ymin = pfs->box_ymin;
 if( ymax ) *ymax = pfs->box_ymax;
 if( zmin ) *zmin = pfs->box_zmin;
 if( zmax ) *zmax = pfs->box_zmax;

 return( 1 );
}

/*======================  pfsRItrParBoundBox  ======================*/

int pfsRItrParBoundBox( double *umin, double *umax,
                       double *vmin, double *vmax )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return( 0 );

 if( umin ) *umin = pfs->box_umin;
 if( umax ) *umax = pfs->box_umax;
 if( vmin ) *vmin = pfs->box_vmin;
 if( vmax ) *vmax = pfs->box_vmax;

 return( 1 );
}

/*=======================  pfsRItrFirstElem  =======================*/

int pfsRItrFirstElem( void )
{
 double xmin, ymin, zmin;
 double xmax, ymax, zmax;

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1) )
  return( PFS_ITERATION_END );

 PFSRtreeInitTraverse( pfs->elm3d_tree );

 cur_elem = PFSRtreeTraverse( pfs->elm3d_tree,
                              &xmin, &xmax, 
                              &ymin, &ymax,
                              &zmin, &zmax );

 cur_elemnode = 0;

 if( cur_elem == NULL )
  return( PFS_ITERATION_END );

 return( PFS_SUCCESS );
}

/*=======================  pfsRItrNextElem  ========================*/

int pfsRItrNextElem( void )
{
 double xmin, ymin, zmin;
 double xmax, ymax, zmax;

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1) )
  return( PFS_ITERATION_END );

 cur_elemnode = 0;

 cur_elem = PFSRtreeTraverse( pfs->elm3d_tree,
                              &xmin, &xmax, 
                              &ymin, &ymax,
                              &zmin, &zmax );

 if( cur_elem == NULL )
  return( PFS_ITERATION_END );

 return( PFS_SUCCESS );
}

/*=======================  pfsRItrFirst3dElem  =======================*/

int pfsRItrFirst3dElem( double xmin, double xmax,
                       double ymin, double ymax,
                       double zmin, double zmax )
{
 double xmn, ymn, zmn;
 double xmx, ymx, zmx;

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1) )
  return( PFS_ITERATION_END );

 PFSRtreeInitSearchBox( pfs->elm3d_tree, xmin, xmax, ymin, ymax, zmin, zmax );

 cur_elem = PFSRtreeSearchBox( pfs->elm3d_tree,
                               &xmn, &xmx, 
                               &ymn, &ymx,
                               &zmn, &zmx );

 cur_elemnode = 0;

 if( cur_elem == NULL )
  return( PFS_ITERATION_END );

 return( PFS_SUCCESS );
}

/*=======================  pfsRItrNext3dElem  ========================*/

int pfsRItrNext3dElem( double xmin, double xmax,
                      double ymin, double ymax,
                      double zmin, double zmax )
{
 double xmn, ymn, zmn;
 double xmx, ymx, zmx;

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1) )
  return( PFS_ITERATION_END );

 cur_elemnode = 0;

 cur_elem = PFSRtreeSearchBox( pfs->elm3d_tree,
                               &xmn, &xmx, 
                               &ymn, &ymx,
                               &zmn, &zmx );

 if( cur_elem == NULL )
  return( PFS_ITERATION_END );

 return( PFS_SUCCESS );
}

/*=======================  pfsRItrFirstParElem  =======================*/

int pfsRItrFirstParElem( double umin, double umax,
                        double vmin, double vmax )
{
 double xmn, ymn, zmn;
 double xmx, ymx, zmx;

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1) )
  return( PFS_ITERATION_END );

 PFSRtreeInitSearchBox( pfs->elm2d_tree, umin, umax, vmin, vmax, 0, 0 );

 cur_elem = PFSRtreeSearchBox( pfs->elm2d_tree,
                               &xmn, &xmx, 
                               &ymn, &ymx,
                               &zmn, &zmx );

 cur_elemnode = 0;

 if( cur_elem == NULL )
  return( PFS_ITERATION_END );

 return( PFS_SUCCESS );
}

/*=======================  pfsRItrNextParElem  ========================*/

int pfsRItrNextParElem( double umin, double umax,
                       double vmin, double vmax )
{
 double xmn, ymn, zmn;
 double xmx, ymx, zmx;

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1) )
  return( PFS_ITERATION_END );

 cur_elemnode = 0;

 cur_elem = PFSRtreeSearchBox( pfs->elm2d_tree,
                               &xmn, &xmx, 
                               &ymn, &ymx,
                               &zmn, &zmx );

 if( cur_elem == NULL )
  return( PFS_ITERATION_END );

 return( PFS_SUCCESS );
}

/*=======================  pfsRSetCurrElem  ========================*/

void pfsRSetCurrElem( void *elem )
{
 if( elem == NULL )
  return;

 cur_elem = (PfsElem *)elem;
 cur_elemnode = 0;
}

/*=========================  pfsRGetElem  ==========================*/

int pfsRGetElem( double x, double y, double z, void **elem )
{
 int itr;
 double u, v;

 for(itr  = pfsRItrFirstElem( );
     itr != PFS_ITERATION_END;
     itr  = pfsRItrNextElem( ) )
 {
   if( pfsParPntCurElem( x, y, z, &u, &v ))
   {
     *elem = cur_elem;
     return(1);
   }
 }
 *elem = NULL;

 return(0);
}

/*=========================  pfsRGetLoopElem  ==========================*/
/* Procura o elemento de um loop externo mais próximo ao ponto dado pelas
 * coordenadas xyz
 */
int pfsRGetLoopElem( double x, double y, double z, void **elem )
{
 int itr;
 int nloop = pfsRItrNumLoops( );
 int ni;   
 int i;
 double xc, yc, zc;
 int first = 1;
 PfsNode *node;
 double dist, dist_cur;
 
 for( i=0; i<nloop; i++ )
 {
  for(itr  = pfsRItrFirstLoopNode( i, &ni );
      itr != PFS_ITERATION_END;
      itr  = pfsRItrNextLoopNode(  i, ni, &ni ) )
  {
        node = &pfs->nodelist[ni];
        cur_elem = node->adj_head->elem;
   pfsRItrElemCenter( &xc, &yc, &zc );
   dist_cur = sqrt( ( x - xc ) * ( x - xc ) +
                    ( y - yc ) * ( y - yc ) +
                    ( z - zc ) * ( z - zc ) );
        /* calcula a distância do pto ao centro do elemento */
   if( first || ( dist_cur < dist ) )
   {
    dist = dist_cur;
    *elem = cur_elem;
    first = 0;
   }
  }
 }
 
 if( first )
 {
  *elem = NULL;
  return(0);
 }
 return(1);
}

/*=======================  pfsRItrElemNode  ========================*/

int pfsRItrElemNode( int *id )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return( PFS_ITERATION_END );

 if( cur_elem == NULL )
  return( PFS_ITERATION_END );

 if( cur_elemnode == 3 )
  return( PFS_ITERATION_END );

 *id = cur_elem->inc[cur_elemnode];

 cur_elemnode += 1;

 return( PFS_SUCCESS );
}

/*======================  pfsRItrElemCenter  =======================*/

void pfsRItrElemCenter( double *x, double *y, double *z )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 if( cur_elem == NULL )
  return;

 *x = cur_elem->center.x;
 *y = cur_elem->center.y;
 *z = cur_elem->center.z;
}

/*======================  pfsRItrElemNormal  =======================*/

void pfsRItrElemNormal( double *nx, double *ny, double *nz )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 if( cur_elem == NULL )
  return;

 *nx = cur_elem->normal.x;
 *ny = cur_elem->normal.y;
 *nz = cur_elem->normal.z;
}

/*======================  pfsRItrElem3dBox  ======================*/

void pfsRItrElem3dBox( double *xmin, double *ymin, double *zmin,
                      double *xmax, double *ymax, double *zmax )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 if( cur_elem == NULL )
  return;

 *xmin = cur_elem->min3dbox.x;
 *ymin = cur_elem->min3dbox.y;
 *zmin = cur_elem->min3dbox.z;
 *xmax = cur_elem->max3dbox.x;
 *ymax = cur_elem->max3dbox.y;
 *zmax = cur_elem->max3dbox.z;
}

/*=====================  pfsRItrElemParBox  ======================*/

void pfsRItrElemParBox( double *umin, double *vmin,
                       double *umax, double *vmax )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 if( cur_elem == NULL )
  return;

 *umin = cur_elem->minparbox.u;
 *vmin = cur_elem->minparbox.v;
 *umax = cur_elem->maxparbox.u;
 *vmax = cur_elem->maxparbox.v;
}

/*=======================  pfsRItrNumNodes  ========================*/

int pfsRItrNumNodes( void )
{
 if( !pfs )
  return( 0 );

 return( pfs->numnodes );
}

/*======================  pfsRItrNodeCoords  =======================*/

void pfsRItrNodeCoords( int id, double *x, double *y, double *z )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 *x = pfs->nodelist[id].coord.x;
 *y = pfs->nodelist[id].coord.y;
 *z = pfs->nodelist[id].coord.z;
}

/*======================  pfsRItrNodeNormal  =======================*/

void pfsRItrNodeNormal( int id, double *nx, double *ny, double *nz )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 *nx = pfs->nodelist[id].normal.x;
 *ny = pfs->nodelist[id].normal.y;
 *nz = pfs->nodelist[id].normal.z;
}

/*======================  pfsRItrNodeJacobian  =======================*/

void pfsRItrNodeJacobian( int id, double *xdu, double *ydu, double *zdu,
                                 double *xdv, double *ydv, double *zdv )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 *xdu = pfs->nodelist[id].ddu.x;
 *ydu = pfs->nodelist[id].ddu.y;
 *zdu = pfs->nodelist[id].ddu.z;
 *xdv = pfs->nodelist[id].ddv.x;
 *ydv = pfs->nodelist[id].ddv.y;
 *zdv = pfs->nodelist[id].ddv.z;
}

/*=====================  pfsRItrNodeParVals  =======================*/

void pfsRItrNodeParVals( int id, double *u, double *v )
{
 if( (!pfs) || (pfs->numnodes < 1) )
  return;

 *u = pfs->nodelist[id].par.u;
 *v = pfs->nodelist[id].par.v;
}

/*=======================  pfsRItrNumLoops  ========================*/

int pfsRItrNumLoops( void )
{
 if( !pfs )
  return( 0 );

 return( pfs->numloops );
}

/*=====================  pfsRItrFirstLoopNode  =====================*/

int pfsRItrFirstLoopNode( int loop_id, int *id )
{
 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1)
            || (pfs->numloops < 1) || (loop_id >= pfs->numloops) )
  return( PFS_ITERATION_END );

 *id = pfs->looplist[loop_id];

 return( PFS_SUCCESS );
}

/*=====================  pfsRItrNextLoopNode  ======================*/

int pfsRItrNextLoopNode( int loop_id, int ref_id, int *id )
{
 PfsNode    *ref_node;           /* pointer to given reference node */

 if( (!pfs) || (pfs->numnodes < 1) || (pfs->numelems < 1)
            || (pfs->numloops < 1) )
  return( PFS_ITERATION_END );

 if( pfs->nodelist[ref_id].type != NODE_BOUNDARY )
  return( PFS_ITERATION_END );

 ref_node = &pfs->nodelist[ref_id];
 *id = ref_node->adj_head->id;

 if( *id == pfs->looplist[loop_id] )
  return( PFS_ITERATION_END );

 return( PFS_SUCCESS );
}

/*========================= pfsRProjectOn ==========================*/

int pfsRProjectOn( double x, double y, double z,
                  double *u, double *v )
{
 double dmin, xmin, xmax, ymin, ymax, zmin, zmax, dist;
 int kmin = 0;
 int itr;
 PFSGeoPoint p;
 PFSGeoPoint clst;
 double tol = TOLER;

 if(pfs == NULL)
  return( 0 );

 dmin = 1e20;

 p.x = x;
 p.y = y;
 p.z = z;

 xmin = x - tol;
 xmax = x + tol;
 ymin = y - tol;
 ymax = y + tol;
 zmin = z - tol;
 zmax = z + tol;

 for(itr  = pfsRItrFirst3dElem(xmin, xmax, ymin, ymax, zmin, zmax);
     itr != PFS_ITERATION_END;
     itr  = pfsRItrNext3dElem (xmin, xmax, ymin, ymax, zmin, zmax))
 {
  dist = pfsDistPntCurElem( &p, &clst );

  if(dist < dmin)
  {
   if( !pfsParPntCurElem( x, y, z, u, v ) )
    continue;

   dmin = dist;
   kmin++;
  }
 }

 if(kmin != 0)
  return( 1 );

 return 0;
}

/*======================== pfsSnapToBound =========================*/

int pfsSnapToBound( double *u, double *v, double toler )
{
 PFSGeoSurfPar p;

 if(pfs == NULL)
  return( 0 );

 p.u = *u;
 p.v = *v;

 if( !pfsSnapParPtToBound( &p ) )
  return( 0 );

 if( Len2( *u - p.u, *v - p.v ) <= toler )
 {
  *u = p.u;
  *v = p.v;
  return( 1 );
 }

 return( 0 );
}

/*======================== pfsRGetJacobian =========================*/

int pfsRGetJacobian( double u, double v,
                    double *xdu, double *ydu, double *zdu,
                    double *xdv, double *ydv, double *zdv,
                    double *xdw, double *ydw, double *zdw )
{
 PFSGeoSurfPar p;
 double umin, vmin, umax, vmax;
 double tol = 1e-3;
 int itr;

 if(pfs == NULL)
  return( 0 );

/* In case given point is outside of surface parametric borders,
 * reject point.  If it is close to the boundary, update the
 * parametric coordinates to the coordinates of the closest point.
  */
 p.u = u;
 p.v = v;
 if( (! pfsIsParPtInside( &p )) && (! pfsIsParPtOnBound( &p )) )
 {
  return( 0 );
 }
 u = p.u;
 v = p.v;
 
 umin = u;
 vmin = v;
 umax = u;
 vmax = v;

 umin -= tol;
 vmin -= tol;
 umax += tol;
 vmax += tol;

 for(itr  = pfsRItrFirstParElem(umin, umax, vmin, vmax);
     itr != PFS_ITERATION_END;
     itr  = pfsRItrNextParElem(umin, umax, vmin, vmax))
 {
  if( pfsParDerivCurElem( u, v, xdu, ydu, zdu, xdv, ydv, zdv, xdw, ydw, zdw ) )
  {
   return( 1 );
  }
 }

 return( 0 );
}

/*============================ pfsREval ============================*/

int pfsREval( double u, double v, double *x, double *y, double *z )
{
 PFSGeoSurfPar p;
 double umin, vmin, umax, vmax;
 double tol = 1e-3;
 int itr;

 if(pfs == NULL)
  return( 0 );

/* In case given point is outside of surface parametric borders,
 * reject point.  If it is close to the boundary, update the
 * parametric coordinates to the coordinates of the closest point.
  */
 p.u = u;
 p.v = v;

 if( (! pfsIsParPtInside( &p )) && (! pfsIsParPtOnBound( &p )) )
 {
  return( 0 );
 }
 u = p.u;
 v = p.v;

 umin = u;
 vmin = v;
 umax = u;
 vmax = v;

 umin -= tol;
 vmin -= tol;
 umax += tol;
 vmax += tol;

 for(itr  = pfsRItrFirstParElem(umin, umax, vmin, vmax);
     itr != PFS_ITERATION_END;
     itr  = pfsRItrNextParElem(umin, umax, vmin, vmax))
 {
  if( pfs3dPntCurElem( u, v, x, y, z ) )
  {
   return( 1 );
  }
 }

 return( 0 );
}

/*===========================  pfsShot  ===========================*/
int pfsRShot( void *init_elem, double x0, double y0, double z0,
             double u_shot, double v_shot, double len,
             void **end_elem, double *x1, double *y1, double *z1 )
{
 int ni, nj;     
 int type;
 double u0, v0;  /* coord paramétrica da posição inicial do tiro */
 double u1, v1;  /* coord paramétrica da posição final do tiro   */
 double u_aux, v_aux;
 PFSGeoPoint p0, p1; /* posição inicial e final, respect.     */
 void *init_elem_bck = init_elem;
 int  status;
 int  in;
 
 if( init_elem == NULL )
 {
  printf( "\nERRO: Shot nao pode ter elemento inicial nulo\n" ); 
  return( 0 );
 }
 if( len < TOLER )
 {
  *x1 = x0;
  *y1 = y0;
  *z1 = z0;
  *end_elem = init_elem;
  return( 1 );
 }

/* atualiza o elemento corrente para q ele seja o elemento inicial
 */
  pfsRSetCurrElem( (PfsElem *)init_elem );

/* consulta a coordenada paramétrica da posição inicial do tiro 
 */
 in = pfsParPntCurElemAttract( &x0, &y0, &z0, &u0, &v0 );

/* Verifica se o ponto está muito próximo a um nó ou aresta e, caso
 * positivo, atrai o pto para a posição do nó ou da sua projeção 
 * na aresta. Neste caso, considera-se o triangulo inicial, que é 
 * dado por init_elem, desconhecido para que a busca seja feita em
 * todos os triângulos adjacentes ao nó ou aresta para onde o pto 
 * foi atraido.
 */
 u_aux = u0;
 v_aux = v0;
 if( pfsChkProxElemEdge( (PfsElem *)init_elem, &u_aux, &v_aux, &type, &ni, &nj ) )
 {
  /* atualiza a direção do tiro, considerando a nova posição do pt atraído */
  u_shot = (u0 + u_shot) - u_aux;
  v_shot = (v0 + v_shot) - v_aux;

  u0 = u_aux;
  v0 = v_aux;

 /* atualiza a coord cartesiana do pto, considerando sua nova posição */
  pfs3dPntCurElem( u0, v0, &x0, &y0, &z0 ); 
  init_elem = NULL;
 }
 else
 {
  if( in ) 
   type = FACE;
  else
  {
   /* pega um nó do loop por onde o tiro saiu */
   PfsElem *elem = (PfsElem *)init_elem;
   int i;

   status = 0;
   for( i=0; i<3; i++ )
   {
    if( pfs->nodelist[elem->inc[i]].type == NODE_BOUNDARY )
    { 
     ni = elem->inc[i];
     type = LOOP;
     type |= FREE_FACE;
     status = 1;
     break;
    }
   }
   if( !status ) 
   {
    printf("\nERRO: %d Não foi encontrado o nó  do loop por onde o tiro saiu.\n", pfs_node);
    return( 0 );
   }
  }     
 }

 p0.x = x0;
 p0.y = y0;
 p0.z = z0;
 status = pfsShotSolve( (PfsElem *)init_elem, u0, v0, p0, u_shot, v_shot, 
                        (PfsElem **)end_elem, &u1, &v1, &p1, &type, &ni, &nj, 
                        &len );
 
 if( status )
 {
  *x1 = p1.x;
  *y1 = p1.y;
  *z1 = p1.z;
  
  if( *end_elem == NULL )
   *end_elem = init_elem_bck;

  return( 1 );
 }
 return( 0 );
}

---