constcurve.c

Go to the documentation of this file.

/*
 * $Id: constcurve.c,v 1.4 2004/06/24 18:32:35 joaoluiz Exp $
 * (C) COPYRIGHT 2004, Joao Luiz Elias Campos.
 * All Rights Reserved
 * Duplication of this program or any part thereof without the express
 * written consent of the author is prohibited
 *
 *   Modified:    01-Aug-05    Alexandre A. Del Savio
 *     Modificada a funcao CheckPoint e criada a funcao LineBounds.
 *     Criada as definicoes: ZERO, MIN e MAX.
 *
 *   Modified:    03-Aug-05    Alexandre A. Del Savio
 *     Criada as funcoes VecLen, Dist, PolySurf e PolyInvSurf.
 *
 *   Modified:    18-Aug-05    Alexandre A. Del Savio
 *     Modificada as funcoes CurveToGlobal e CurveToLocal para usarem 
 *     diretamente o inverso do jacobiano e o jacobiano, respectivamente, 
 *     do elemento e nao mais a combinacao destes dos nos.
 *
 *   Modified:    25-Aug-05    Alexandre A. Del Savio
 *     Modificada a funcao ConstCurveInit e criada as funcoes 
 *     ConstCurveGlobalParametric e CurveToGlobalParametric. 
 *
 *   Modified:    01-Sep-05    Alexandre A. Del Savio
 *     Modificada as funcoes ConstCurveNew e ConstCurveRead de modo
 *     que as restricoes sejam estendidas. 
 *
 *   Modified:    25-Oct-05    Alexandre A. Del Savio
 *     Implementada uma tolerância adaptativa que se modifica em
 *     função das geometrias das curvas lidas por este pacote.
 *     Para tanto, criou-se uma função para calcular esta tolerância, 
 *     EvaluateTolerance, e adicionou-se uma chamada a esta função na
 *     função que efetua a leitura das curvas, ConstCurveRead.
 *
 *   Modified:    23-Fev-06    Juan Pablo Ibañez
 *     modificada a funcões CurveToLocal, ConstCurveLocal, CurveToGlobal, 
 *     CurveToGlobalParametric, ConstCurveGlobalParametric, ConstCurveGlobal,
 *     que passam a retornar 0 ou 1 em caso de falha ou sucesso respectivamente.
 *
 *  Modified:    29-Mar-06    Juan Pablo Ibañez
 *     Modificada a funcão ConstCurveNew que passa a receber mais um
 *     parâmetro, "i".
 *
 *   Modified:    11-Apr-06    Alexandre A. Del Savio
 *     Adicionada as funcoes: DiffVector, CrossProd e Invert3x3.
 *     Re-organizada as funcoes: EvaluateTolerance, CurveToLocal e 
 *     CurveToGlobal.
 *     Modificada as funcoes: LineBounds, VecLen, Dist, VecNormalize,
 *     DirNormalize, CheckPoint, PolySurf e ConstCurveRead.
 *     Estas modificacoes e novas funcoes tem como objetivo habilitar estas
 *     restricoes no espaco 3D. Ou seja, a partir de agora as restricoes de 
 *     linha podem ser dadas em qualquer posicao do espaco.
 *     Além disso, compatibilizou-se este arquivo com a versao do Juan Pablo.
 *     Funcoes compatibilizadas: CurveToLocal, ConstCurveLocal, CurveToGlobal, 
 *     CurveToGlobalParametric, ConstCurveGlobalParametric e ConstCurveGlobal.
 *
 *   Modified:    24-Mai-06    Juan Pablo Ibañez
 *     Modificada a funcão CurveToGlobalParametric, onde corrigui-se o 
 *     controle para saber se delta é negativo ou positivo.
 *
 *   Modified:    21-Jun-06    Alexandre A. Del Savio
 *     Modified the ConstCurveNew function when it was removed the first 
 *     parameter "int i" and it started to use the ContNumConstrainRead
 *     variable.
 *
 *    Modified:    12-Jul-06    Juan Pablo Iba
 *     Modificadas as funcoes PolySurf e PolyInvSurf onde foram definidas
 *     tolerancias para o cálculo.
 *
 *    Modified:    12-Ago-06    Juan Pablo Iba
 *     Criadas as funcoes ConstShotNew e ConstShotRead para o caso de constrain
 *     tipo SHOT (apoio inclinado).
 *
 *    Modified     14-Set-2006  Juan Pablo Ibanez
 *     Foram removidas as funcoes ConstShotNew e ConstShotRead que passaram
 *     ao modulo constshot.c
 *
 *    Modified     14-Dec-2006  Juan Pablo Ibanez
 *     Modificada a função _curveSlide que deixou de calcular o
 *     sinal do deslocamento porque o vetor de deslocamento 
 *     recebido já está no sistema local do elemento da curva.
 *
 *    Modified     02-Apr-2007  Gisele Holtz
 *     Modificada as funções _curvePolySlide, PolyInvSurf, _curveToLocal, 
 *     _curveToGlocal para prolongar o primeiro e o último
 *     segmento da curva, permitindo projetar sobre a curva um comprimento maior
 *     que o comprimento da curva.  
 *     Criada a função _checkPointInExtrem para verificar se um pto está sobre 
 *     os segmentos extremos da curva.
 *     
 */

/*
 * Global variables and symbols:
 */
#include <stdlib.h>
#include <math.h>
#include "constrain.h"
#include "rio.h"

typedef struct _curvedata 
{
  double jac[3][3];
  double ijac[3][3];
} sCCurveData;


/*
** ------------------------------------------------------------------------
** Local variables and symbols:
*/

/*
%K Tolerance definition
*/
static double Tolerance = 1.0e+10;

/*
%K Zero definition
*/
#ifndef ZERO
#define ZERO(v)     ((fabs(v)<Tolerance)?1:0) 
#endif

/*
%K Minimum definition
*/
#ifndef MIN
#define MIN(a,b)    ((a<b)?a:b)
#endif

/*
%K Maximum definition
*/
#ifndef MAX
#define MAX(a,b)    ((a>b)?a:b)
#endif


/*
** ------------------------------------------------------------------------
** Local functions:
*/

static void     DiffVec( sCoord *in0, sCoord *in1, sCoord *out )    
{
  out->x = in0->x - in1->x;
  out->y = in0->y - in1->y;
  out->z = in0->z - in1->z;
} 

static double ProdEscalar( sCoord *in0, sCoord* in1 )
{
  double prod;
  prod = (in0->x)*(in1->x) + (in0->y)*(in1->y) + (in0->z)*(in1->z);

  return (prod);
}

static void CrossProd( sCoord *in0, sCoord *in1, sCoord *out )
{
  out->x = (in0->y * in1->z) - (in0->z * in1->y); 
  out->y = (in0->z * in1->x) - (in0->x * in1->z); 
  out->z = (in0->x * in1->y) - (in0->y * in1->x); 
}  

static double VecLen( sCoord *vec )
{
  return(sqrt((vec->x * vec->x) + (vec->y * vec->y) + (vec->z * vec->z)));
}

static double Dist( sCoord *in0, sCoord *in1 )
{
  sCoord tmp;

  DiffVec(in1, in0, &tmp);
  
  return(VecLen(&tmp));
}

static void  VecNormalize( sCoord *vec, sCoord *n )
{
  double      len;

  len = VecLen(vec);

  n->x = vec->x / len;
  n->y = vec->y / len;
  n->z = vec->z / len;
}

static sCoord DirNormalize( sCoord* p1, sCoord* p2 )
{
  sCoord  dr;    /* Vetor direção da reta (p1,p2). */
  sCoord  nrm;

  dr.x = p2->x - p1->x;
  dr.y = p2->y - p1->y;
  dr.z = p2->z - p1->z;

  VecNormalize(&dr,&nrm);

  return nrm;
}

static int LineBounds( sCoord *c, sCoord *pi, sCoord *pf )
{
  double aux;

  aux = MIN(pi->x, pf->x);
  if(c->x < aux) if(!ZERO(c->x - aux)) return(0);
  aux = MAX(pi->x, pf->x);
  if(c->x > aux) if(!ZERO(c->x - aux)) return(0);

  aux = MIN(pi->y, pf->y);
  if(c->y < aux) if(!ZERO(c->y - aux)) return(0);
  aux = MAX(pi->y, pf->y);
  if(c->y > aux) if(!ZERO(c->y - aux)) return(0);

  aux = MIN(pi->z, pf->z);
  if(c->z < aux) if(!ZERO(c->z - aux)) return(0);
  aux = MAX(pi->z, pf->z);
  if(c->z > aux) if(!ZERO(c->z - aux)) return(0);

  return(1);
}

static int Invert3x3( double a[3][3], double b[3][3] )
{
  double det;
  int    i, j;
        
  /* Limpa a matriz inversa. */
  for(i = 0; i < 3; i++)
    for(j = 0; j < 3; j++) 
      b[i][j] = 0.0 ;

  /* Calcula o determinante da matriz. */
  det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
        a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
        a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]) ;

  /* Se não foi possível inverter a matriz retorna 0. */
  if(ZERO(det)) 
    return 0;

  /* Calcula os coeficientes da matriz inversa. */
  b[0][0] =  (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det ;
  b[0][1] = -(a[0][1] * a[2][2] - a[0][2] * a[2][1]) / det ;
  b[0][2] =  (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det ;
  b[1][0] = -(a[1][0] * a[2][2] - a[1][2] * a[2][0]) / det ;
  b[1][1] =  (a[0][0] * a[2][2] - a[0][2] * a[2][0]) / det ;
  b[1][2] = -(a[0][0] * a[1][2] - a[0][2] * a[1][0]) / det ;
  b[2][0] =  (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det ;
  b[2][1] = -(a[0][0] * a[2][1] - a[0][1] * a[2][0]) / det ;
  b[2][2] =  (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det ;

  return 1;
}

static int CheckPoint( sCoord *c, sCoord *pi, sCoord *pf )
{
  sCoord du, dv, dw;
  double coef[3][3];
  double cinv[3][3];

  /* Direções. */
  DiffVec(pf, pi, &du);
  DiffVec(c,  pi, &dv);

  /* Direção da linha perpendicular. */
  CrossProd(&du, &dv, &dw);

  /* Constroi os coeficientes da matriz. */
  coef[0][0] = du.x;  coef[0][1] = dv.x;  coef[0][2] = dw.x;
  coef[1][0] = du.y;  coef[1][1] = dv.y;  coef[1][2] = dw.y;
  coef[2][0] = du.z;  coef[2][1] = dv.z;  coef[2][2] = dw.z;
 
  /* Inverte a matriz dos coeficientes. */
  if(!Invert3x3(coef, cinv)) 
  {
    /* Como nao foi possível inverter a matriz significa que o determinante é
     * nulo. Ou seja, os vetores sao paralelos ou colineares. Desta forma,
     * basta testar se o ponto está sobre a reta testando a sua bounding box.
     */
    /*printf("lineBounds...\n");*/
    if(!LineBounds(c, pi, pf))
      return(0);
    else
      return(1);
  }

  return(0);
}


static int _checkPointInExtrem(sConstrain *c, sCoord *p, int *id )
{
  sCoord d, dp, n;
  double area;
  sCoord *pi;
  sCoord *pf;
  int    inc[2];
  double TOL = 0.001;
  
  /* Verifica primeiro segmento  da curva
   *
   *          d                           
   *      i  <--   f                     
   *      .----------.       .........    . ----------. 
   */
  ElmConnect(c->elems[0], inc);
  pi = &c->nodes[inc[0]-1].coord;
  pf = &c->nodes[inc[1]-1].coord;
  d = DirNormalize(pi, pf);
  dp = DirNormalize( p, pf );
  
  /* Verifica se o ponto está na direção do segmento */
  CrossProd(&d, &dp, &n);
  area = VecLen( &n );
  
  /* Verifica se o pto está no mesmo sentido do segmento */
  if( area < TOL )
  {
        if( ProdEscalar( &d, &dp ) > 0 )
        {
         /* ponto está antes do primeiro ponto da curva */
    *id = 0;
    return( 1 );
   }
  }
  
   /* Verifica último segmento  da curva   
   *                                             d
   *                                        i    -->   f
   *      .----------.       .........      . ----------.
   */
  ElmConnect(c->elems[c->numElems-1], inc);
  pi = &c->nodes[inc[0]-1].coord;
  pf = &c->nodes[inc[1]-1].coord;
  d = DirNormalize(pf, pi);
  dp = DirNormalize( p, pi );

  /* Direção da linha perpendicular. */
  CrossProd(&d, &dp, &n);
  area = VecLen( &n );

  if( area < TOL )
  {
        /* Ponto está na mesma direção do último segmento */
        if( ProdEscalar( &d, &dp ) > 0 )
        {
         /* ponto está depois do último ponto da curva */
    *id = c->numElems-1;
    return( 1 );
   }
  }
  return( 0 );
}

static int PolyInvSurf( sConstrain *c, sCoord *p, double *dist )
{
  int  flag = 0;
  int  e;
  int  inc[2];
  int  id;
  double tol = 0.000001;
  
  /* Inicializa a distancia.
   */
  *dist = 0.0;

  /* Testa se o ponto p é o primeiro ponto inicial da curva (poligonal).
   */
  ElmConnect(c->elems[0], inc);
  
  if((Dist(&c->nodes[inc[0]-1].coord, p)) < tol)
  {
    id = 0;
    *dist = 0.0;
    return(id);
  }

  /* Encontra o segmento onde esta o ponto p, acumulando a distancia
   * ate chegar a ele.
   */
  for(e = 0; e < c->numElems; e++)
  {
    ElmConnect(c->elems[e], inc);
    if(CheckPoint(p, &c->nodes[inc[0]-1].coord, &c->nodes[inc[1]-1].coord))
    {
      id = e;
      flag = 1;
      *dist += Dist(&c->nodes[inc[0]-1].coord, p);
      break;
    }
    *dist += Dist(&c->nodes[inc[0]-1].coord, &c->nodes[inc[1]-1].coord);
  }

  /* O ponto p não está sobre a poligonal.
  */
  if(!flag)
  {
    flag = _checkPointInExtrem(c, p, &id );
    
    if( !flag )
    {
     id = -1;
     *dist = 0.0;
     return(id);
    }
    ElmConnect(c->elems[id], inc);
    
    if( id == 0 )
        *dist = 0.0;
    else 
      *dist += Dist(&c->nodes[inc[1]-1].coord, p);
  }
   
  return(id);
}

static int _curvePolySlide( sConstrain *c, double dist, void **elm, sCoord *p )
{
  double  dist_pol = 0.0;
  int     flag     = 0;
  double  D, Dx, Dy, Dz, d;
  int     e;
  int     inc[2];
  double tol = 0.000001;
  int     first = 0; /* flag q indica q o pto está sobre a prolongação do primeiro
                      * segmento */

  /* A distância passada como parâmetro é negativa.
  */
  if(dist < 0.0)
  {
    /* printf( "_curvePolySlide: distância negativa\n" ); */
    /* mantains current element node constrain */
    ElmConnect(c->elems[0], inc);
    *elm = c->elems[0];
    first = 1;
    dist_pol = Dist(&c->nodes[inc[1]-1].coord, &c->nodes[inc[0]-1].coord);
  }

  /* Distância passada como parâmetro igual a zero.
  ** Retorna o primeiro ponto da curva (poligonal).
  */
  if(fabs( dist ) < tol)
  {
   /* printf("_curvePolySlide: distancia nula\n");*/
    ElmConnect(c->elems[0], inc);
    p->x = c->nodes[inc[0]-1].coord.x;
    p->y = c->nodes[inc[0]-1].coord.y;
    p->z = c->nodes[inc[0]-1].coord.z;
    *elm = c->elems[0];

    return(0);
  }
  if( !first )
  {
    /* Encontra o segmento onde esta o ponto p.
     */
    for(e = 0; e < c->numElems; e++)
    {
      ElmConnect(c->elems[e], inc);
      dist_pol += Dist(&c->nodes[inc[0]-1].coord, &c->nodes[inc[1]-1].coord);
      if((dist_pol > dist) || (fabs(dist_pol - dist)<0.00001))
      {
        /*printf("distancia achada\n");*/
        flag = 1;
        *elm = c->elems[e];   
        break;
      }
    }
  }

  /* A distância passada como parâmetro é maior que o comprimento
  ** de toda a curva (polyline).
  */
  if( !first && !flag )
  {
    /*printf("_curvePolySlide: distancia maior do que a curva\n");*/
    /* mantains current element node constrain */
    *elm = c->elems[c->numElems-1];
    ElmConnect(c->elems[c->numElems-1], inc);
  }
  
  D  = Dist(&c->nodes[inc[1]-1].coord, &c->nodes[inc[0]-1].coord);
  Dx = c->nodes[inc[1]-1].coord.x - c->nodes[inc[0]-1].coord.x;
  Dy = c->nodes[inc[1]-1].coord.y - c->nodes[inc[0]-1].coord.y;
  Dz = c->nodes[inc[1]-1].coord.z - c->nodes[inc[0]-1].coord.z;
  d  = dist - dist_pol;
  
  p->x = c->nodes[inc[1]-1].coord.x + (d * Dx) / D;
  p->y = c->nodes[inc[1]-1].coord.y + (d * Dy) / D;
  p->z = c->nodes[inc[1]-1].coord.z + (d * Dz) / D;
  
  return(1);
}

static void EvaluateTolerance( sConstrain *c )
{
  int inc[3];
  double dist = 0.0;
  double auxTol = 0.0;
  int i;

  for(i = 0; i < c->numElems; i++)
  {
    ElmConnect(c->elems[i],inc);
    dist = Dist(&c->nodes[inc[1]-1].coord, &c->nodes[inc[0]-1].coord);
    auxTol = dist * 0.1;
    if(auxTol < Tolerance)
      Tolerance = auxTol;
  }

} /* End of EvaluateTolerance */

static int _curveToLocal(sConstrain *c, sCoord *p, void **elm, double jac[3][3])
{
  int    e, i, j;
  int    inc[2];
  double jace[3][3], ijace[3][3];
  int    found = 0;

  if( !(*elm) )
  {
    for(e = 0; e < c->numElems; e++) 
    {
      ElmConnect(c->elems[e],inc);
      if(CheckPoint(p, &c->nodes[inc[0]-1].coord, &c->nodes[inc[1]-1].coord)) 
      {
        ElmJacobian(c->elems[e],jace,ijace);
        *elm = c->elems[e];
        found = 1;
        break;
      }
    }

    if(!found)
    {
      found = _checkPointInExtrem(c, p, &e );

      if( !found )
      {
        fprintf(stderr, "Could not find element(1).\n\n");
        *elm = NULL; /* not necessary */
        return 0;
      }
      ElmJacobian(c->elems[e],jace,ijace);
      *elm = c->elems[e];
    }
  }
  else
  {
    sElement *_elm = (sElement*)(*elm);

    ElmJacobian( _elm, jace, ijace );
  }

  for( i = 0; i < 3; i++ ) 
    for( j = 0; j < 3; j++ ) 
      jac[i][j] = jace[i][j];

  return 1;

} /* End of _curveToLocal */

static int _curveToGlobal( sConstrain *c, sCoord *p, 
                           void **elm, double ijac[3][3] )
{ 
  int    i, j;
  int    inc[2];
  double jace[3][3], ijace[3][3]; 
  int    found = 0;

  if( !(*elm) )
  {
    for(i = 0; i < c->numElems; i++) 
    {
      ElmConnect(c->elems[i],inc);
      if(CheckPoint(p, &c->nodes[inc[0]-1].coord, &c->nodes[inc[1]-1].coord))
      {
        ElmJacobian(c->elems[i],jace,ijace);
        *elm = c->elems[i];
        found = 1;
        break;
      } 
    }

    if(!found)
    {
      found = _checkPointInExtrem(c, p, &i );

      if( !found )
      {
        fprintf(stderr, "Could not find element(2).\n\n");
        *elm = NULL; 
        return 0;
      }
      ElmJacobian(c->elems[i],jace,ijace);
      *elm = c->elems[i];
    }
  }
  else
  {
    sElement *_elm = ((sElement*)(*elm));

    ElmJacobian( _elm, jace, ijace );
  }

  for( i = 0; i < 3; i++ )
    for( j = 0; j < 3; j++ ) 
      ijac[i][j] = ijace[i][j];

  return 1;

} /* End of _curveToGlobal */

static int _curveSlide( sConstrain *c, sCoord *p, void **elm, sCoord *d,
                        double ijac[3][3], sCoord *pl )
{
  int    i, j;
  int    inc[2];
  double jace[3][3], ijace[3][3];
  double dist_p;  /* Distancia do inicio da curva ate o ponto p.            */
  double dist_pl; /* Distancia do inicio da curva ate o ponto pl.           */
  int    id_p;    /* Identificacao do segmento onde encontra-se o ponto p.  */
  int    status;  /* buleano: encontrou o  ponto pl em algum elemento?     */

  /* Calcula as distancias. */
  id_p = PolyInvSurf(c, p, &dist_p);
  if(id_p == -1)
  {
    fprintf(stderr, "Could not find element(3).\n\n");
    return 0;
  }

  /* calcula distância no espaço da restrição */
  dist_pl = dist_p + d->x;

  /* Encontra o segmento em funcao da distancia ate o ponto pl e
   * encontra o ponto correspondente percorrendo o espaco da restricao.
   */
  status = _curvePolySlide(c, dist_pl, elm, pl);
  if(!status)
  { 
    pl->x = p->x;
    pl->y = p->y;
    pl->z = p->z;
  }

  /* Calcula os parametros em funcao do segmento encontrado. */
  {
    sElement *_elm = (sElement*)(*elm);

    ElmConnect( _elm, inc );
    ElmJacobian( _elm, jace, ijace );
  }
 
  for(i = 0; i < 3; i++)
    for(j = 0; j < 3; j++) 
      ijac[i][j] = ijace[i][j];

  return 1;

} /* End of _curveSlide */

/*
**--------------------------------------------------------------------------
*/

/*
 * Subclass methods:
 */
static void ConstCurveNew(int, int, int, sConstrain **, sConstrain **);
static void ConstCurveFree(sConstrain *);
static void ConstCurveRead(sConstrain *);
static void ConstCurveBuild(sConstrain *);

static void ConstCurveNew( int label, int num_nodes, int num_elems,
                           sConstrain **c, sConstrain **clst )
{
 sCCurveData *data = NULL;

 /* A partir de agora sera adicionado a curva de restricao mais dois
  * segmentos que correspondem a extensao do primeiro e do ultimo
  * segmento da restricao, na direcao dada pelos proprios segmentos que
  * estao sendo estendidos.
  */
 num_nodes += 2;
 num_elems += 2;

 data = (sCCurveData *)calloc(num_nodes, sizeof(sCCurveData));

 (*c) = (sConstrain *)malloc(sizeof(sConstrain));
 (*c)->type = CURVE;
 (*c)->label = label;
 (*c)->numNodes = num_nodes;
 (*c)->numElems = num_elems;
 (*c)->nodes = (sNode *)calloc(num_nodes, sizeof(sNode));
 (*c)->elems = (sElement **)calloc(num_elems, sizeof(sElement *));
 (*c)->data = data;
 
 clst[ContNumConstrainRead] = (*c);
 ContNumConstrainRead++;

} /* End of ConstCurveNew */

static void ConstCurveFree(sConstrain *c)
{
 int i;

 if( c->data ) {
  free(c->data);
 }

 if( c->numNodes != 0 ) {
  free(c->nodes);
 }

 if( c->numElems != 0 ) {
  for( i = 0; i < c->numElems; i++ ) {
   ElmFree(c->elems[i]);
  }
  free(c->elems);
 }
 c->numNodes = 0;
 c->numElems = 0;
 c->nodes = NULL;
 c->elems = NULL;
 c->data = NULL;

} /* End of ConstCurveFree */

static void ConstCurveRead(sConstrain *c)
{
 int i, id;
 double x, y, z;
 int matid, tckid, ordid;
 sElement *elm = NULL;

#if 0
 static int first = 0;
 if(!first)
 {
  for( i = 0; i < c->numNodes; i++ ) 
  {
    fscanf(nf, "%d %lf %lf %lf", &id, &x, &y, &z);
    c->nodes[i].id = id;
    c->nodes[i].coord.x = x;
    c->nodes[i].coord.y = y;
    c->nodes[i].coord.z = z;
    c->nodes[i].dof.x = FORCE;
    c->nodes[i].dof.y = FORCE;
    c->nodes[i].dof.z = FORCE;
    c->nodes[i].dof.psi = 0.0;
    c->nodes[i].rezone = NONE;
    c->nodes[i].curve = 0;
    c->nodes[i].dspIteration = DSP_NO;
    c->nodes[i].dspTime = DSP_NO;
  }

  for( i = 0; i < c->numElems; i++ ) 
  {
    fscanf(nf, "%d%d%d%d", &id, &matid, &tckid, &ordid);
    ElmNew(LINE2,id,matid,ordid,tckid,&elm,c->elems,c->nodes);
    ElmRead(elm);
  }
  first = 1;
 }
 else
 {
#else
  /* A partir de agora sera adicionado a curva de restricao mais dois
   * segmentos que correspondem a extensao do primeiro e do ultimo
   * segmento da restricao, na direcao dada pelos proprios segmentos que
   * estao sendo estendidos.
   */

  /* Le as coordenadas dos nos "internos".
   */ 
  for(i = 1; i < c->numNodes - 1; i++) 
  {
    fscanf(nf, "%d %lf %lf %lf", &id, &x, &y, &z);
    c->nodes[i].id = id+1;
    c->nodes[i].coord.x = x;
    c->nodes[i].coord.y = y;
    c->nodes[i].coord.z = z;
    c->nodes[i].dof.x = FORCE;
    c->nodes[i].dof.y = FORCE;
    c->nodes[i].dof.z = FORCE;
    c->nodes[i].dof.psi = 0.0;
    c->nodes[i].rezone = NONE;
    c->nodes[i].curve = 0;
    c->nodes[i].dspIteration = DSP_NO;
    c->nodes[i].dspTime = DSP_NO;
  }  
 
  /* Le os segmentos "internos".
   */ 
  for(i = 1; i < c->numElems - 1; i++) 
  {
    fscanf(nf, "%d%d%d%d", &id, &matid, &tckid, &ordid);
    ElmNew(LINE2, id+1, matid, ordid, tckid, &elm, c->elems, c->nodes);
    ElmRead(elm);
  }

  /* Corrige as incidencias dos elementos para permitir a adicao de um
   * novo segmento no inicio e no fim.
   */
  for(i = 1; i < c->numElems - 1; i++) 
  {
    ElmUpdateConnect(c->elems[i]);
  }

  /* Calcula a extensao do primeiro e do ultimo segmento,
   * consequentemente ocorre a criacao de mais dois nos e dois
   * segmentos.
   */
  {
    sCoord dirn;   /* Direcao normalizada do vetor dado pelo segmento. */
    int    inc[2]; /* Vetor auxiliar para capturar a incidencia do segmento. */
    sCoord pe;     /* Ponto estendido. */
    double extension = 10.0; /* Fator de extensao. */
    int    lastElem;         /* Ultimo segmento.   */
    int    lastNode;         /* Ultimo no.         */
   
    /* Primeiro segmento. 
     */ 
    ElmConnect(c->elems[1], inc);
    dirn = DirNormalize(&c->nodes[inc[1]-1].coord, &c->nodes[inc[0]-1].coord);
    pe.x = c->nodes[inc[0]-1].coord.x + extension * dirn.x;
    pe.y = c->nodes[inc[0]-1].coord.y + extension * dirn.y;
    pe.z = c->nodes[inc[0]-1].coord.z + extension * dirn.z;
    c->nodes[0].id = 1;
    c->nodes[0].coord.x = pe.x;
    c->nodes[0].coord.y = pe.y;
    c->nodes[0].coord.z = pe.z;
    c->nodes[0].dof.x = FORCE;
    c->nodes[0].dof.y = FORCE;
    c->nodes[0].dof.z = FORCE;
    c->nodes[0].dof.psi = 0.0;
    c->nodes[0].rezone = NONE;
    c->nodes[0].curve = 0;
    c->nodes[0].dspIteration = DSP_NO;
    c->nodes[0].dspTime = DSP_NO;
    ElmNew(LINE2, 1, 1, 1, 1, &elm, c->elems, c->nodes);
    inc[0] = 1;
    inc[1] = 2;
    ElmSetConnect(elm, inc);

    /* Ultimo segmento. 
     */
    lastElem = c->numElems - 1;
    lastNode = c->numNodes - 1;
    ElmConnect(c->elems[lastElem-1], inc);
    dirn = DirNormalize(&c->nodes[inc[0]-1].coord, &c->nodes[inc[1]-1].coord);
    pe.x = c->nodes[inc[1]-1].coord.x + extension * dirn.x;
    pe.y = c->nodes[inc[1]-1].coord.y + extension * dirn.y;
    pe.z = c->nodes[inc[1]-1].coord.z + extension * dirn.z;
    c->nodes[lastNode].id = c->numNodes;
    c->nodes[lastNode].coord.x = pe.x;
    c->nodes[lastNode].coord.y = pe.y;
    c->nodes[lastNode].coord.z = pe.z;
    c->nodes[lastNode].dof.x = FORCE;
    c->nodes[lastNode].dof.y = FORCE;
    c->nodes[lastNode].dof.z = FORCE;
    c->nodes[lastNode].dof.psi = 0.0;
    c->nodes[lastNode].rezone = NONE;
    c->nodes[lastNode].curve = 0;
    c->nodes[lastNode].dspIteration = DSP_NO;
    c->nodes[lastNode].dspTime = DSP_NO;
    ElmNew(LINE2, c->numElems, 1, 1, 1, &elm, c->elems, c->nodes);
    inc[0] = inc[0] + 1;
    inc[1] = inc[1] + 1;
    ElmSetConnect(elm, inc); 
 }
#endif

 EvaluateTolerance(c);

} /* End of ConstCurveRead */

static void ConstCurveBuild(sConstrain *c)
{
 int i, j, n, k, l;
 int inc[3];
 double jac[3][3], ijac[3][3];
 int *count = NULL;
 sCCurveData *data = NULL;

 count = (int *)calloc(c->numNodes, sizeof(int));
 data = (sCCurveData *)c->data;

 for( i = 0; i < c->numElems; i++ ) {
  ElmNumNodes(c->elems[i],&n);
  ElmConnect(c->elems[i],inc);
  ElmJacobian(c->elems[i],jac,ijac);
  for( j = 0; j < n; j++ ) {
   for( k = 0; k < 3; k++ ) {
    for( l = 0; l < 3; l++ ) {
     data[inc[j]-1].jac[k][l] += jac[k][l];
     data[inc[j]-1].ijac[k][l] += ijac[k][l];
    }
   }
   count[inc[j]-1]++;
  }
 }

 for( i = 0; i < c->numNodes; i++ ) {
  for( k = 0; k < 3; k++ ) {
   for( l = 0; l < 3; l++ ) {
    data[i].jac[k][l] /= count[i];
    data[i].ijac[k][l] /= count[i];
   }
  }
 }
 
 free(count);

} /* End of ConstCurveBuild */

static int ConstCurveLocal(sConstrain *c, sCoord *p, void **elm,
                           double *g, double *l)
{
  int i, j;
  double jac[3][3];

  if(!_curveToLocal(c, p, elm, jac)) 
    return 0;

  for(i = 0; i < NDof; i++) 
  {
    l[i] = 0.0;
    for(j = 0; j < NDof; j++) 
      l[i] += jac[i][j] * g[j];
  }

  return 1;

} /* End of ConstCurveLocal */

static int ConstCurveGlobal(sConstrain *c, sCoord *p, void **elm, 
                            double *l, double *g)
{
  int i, j;
  double ijac[3][3];

  if(!_curveToGlobal(c, p, elm, ijac)) 
    return 0;

  for(i = 0; i < NDof; i++) 
  {
    g[i] = 0.0;
    for( j = 0; j < NDof; j++ )
      g[i] += ijac[i][j] * l[j];
  }

  return 1;

} /* End of ConstCurveGlobal */

static int ConstCurveSlide( sConstrain *c, sCoord *p, void **elm, sCoord *d,
                            double *l, double *g, sCoord *pl )
{
  int i, j;
  double ijac[3][3];

  if(!_curveSlide(c, p, elm, d, ijac, pl)) 
    return 0;

  for(i = 0; i < NDof; i++)
  {
    g[i] = 0.0;
    for(j = 0; j < NDof; j++)
      g[i] += ijac[i][j] * l[j];
  }

  return 1;

} /* End of ConstCurveSlide */

/*
 * Public functions:
 */

void ConstCurveInit(void);
void ConstCurveInit(void)
{
  ConstClass[CURVE].new    = ConstCurveNew;
  ConstClass[CURVE].free   = ConstCurveFree;
  ConstClass[CURVE].read   = ConstCurveRead;
  ConstClass[CURVE].build  = ConstCurveBuild;
  ConstClass[CURVE].local  = ConstCurveLocal;
  ConstClass[CURVE].global = ConstCurveGlobal;
  ConstClass[CURVE].slide  = ConstCurveSlide;

  Tolerance = 1.0e+10;

} /* End of ConstCurveInit */

/*
 * End of File.
 */

---