fem.c

Go to the documentation of this file.

/*
%M This modules contains the Finite Element Method methods and definitions
%a Joao Luiz Elias Campos.
%d September 4th, 1998.
%r $Id: fem.c,v 1.3 2004/06/24 18:32:35 joaoluiz Exp $
%w (C) COPYRIGHT 1995-1996, Eduardo Nobre Lages.
   (C) COPYRIGHT 1997-1999, 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:    25-Aug-05    Alexandre A. Del Savio
 *     Modificada a função DisplacementVelocity onde implementou-se um 
 *     espaço paramétrico auxiliar para encontrar o novo ponto da 
 *     interação seguinte já sobre a restrição, sem introduzir
 *     desequilíbrio no sistema.
 *
 *   Modified:    30-Aug-05    Alexandre A. Del Savio
 *     Modificada a função DisplacementVelocity onde implementou-se o
 *     uso do nó piloto. Após a primeira interação os nós que possuem o 
 *     atributo DISPLACEMENT passam a ter o atributo OBJECT nas
 *     interações seguintes com exceção do nó piloto, para 
 *     "NodeVector[i].dof.{x,y,x}". 
 *
 *   Modified:    08-Sep-05    Alexandre A. Del Savio
 *     Modificada a função DisplacementVelocity onde corrigiu-se a 
 *     atualização do UVector e do VVector.
 *
 *   Modified:    23-Fev-06    Juan Pablo Ibañez
 *     Modificada a funcão DisplacementVelocity, que passa a retornar
 *     um valor int.
 *
 *   Modified:    23-Fev-06    Juan Pablo Ibañez
 *     Modificada a funcão DisplacementVelocity, onde a variável "id"
 *     é obtido segundo o valor de NumConstrains.     
 *
 *   Modified:    23-Fev-06    Juan Pablo Ibañez
 *     Modificada a funcão DisplacementVelocity, onde acessa-se o
 *     vetor ConstList em funcão do tipo de transformacão.
 *
 *   Modified:    24-Mai-06    Juan Pablo Ibañez
 *     Modificada a funcão DisplacementVelocity, onde se faz o 
 *     tratamento da "shear constrain" que é aplicada aos nós 
 *     no modo Pure Shear.
 *
 *   Modified:    23-Jun-06    Alexandre A. Del Savio
 *     Modified the DisplacementVelocity function when it was 
 *     corrected to treat curve and surface type constrains.
 *
 *   Modified:    5-Ago-06    Juan Pablo Ibanez
 *     Modificada a funcao DisplacementVelocity para permitir processar
 *     as restricoes independentemente do tipo de transformacao.
 *
 *   Modified     14-Set-06   Juan Pablo Ibanez
 *     Modificada a funcao DisplacementVelocity, onde adaptou-se a mesma
 *     as mudancas na struct sNode (NodeVector) no campo constrain.
 *
 *   Modified     18-Set-06   Márcio Santi & Gisele Cunha & Juan Pablo Ibañez
 *     When changed API for constrain methods that now gets the identifier
 *     for the current constrained node element. 
 *     Changed constrain methods name from GlobalParametric to Slide
 *
 *   Modified     27-Set-06   Márcio Santi e Juan Pablo Ibañez
 *     When reviewed all necessary routines for implementing the constrain
 *     modules. Here in fem was removed the COnstToGlobal call for any
 *     iteration excepts iteration 0 (zero).
 *     ConstSlide function now recieves vectors in local element constrain
 *     coordinates instead of global coordinates as it was used before.
 *     Reorganized code generaly.
 *
 *   Modified:    Dez-06 André Luis Muller
 *     Criadas as funções para analise com o método
 *     gradiente conjugado. (FirstSolution)- ElmKpr, ElmKp, GCSolver
 *
 *   Modified:    Fev-07 André Luis Muller
 *     Criadas as funções para analise com o método
 *     BFGS. (BFGS)- SearchDirection
 *
 *   Modified:    Jun-07 André Luis Muller
 *     Criadas as funções para analise com algoritmos implícitos
 *     LINEAR, BFGS e NRM.
 *
 *   Modified:    Jul-07 André Luis Muller
 *     Modificadas as funções ElmKpr e ElmKp para utilizar o
 *     formato CSR na armazernagem da matriz de rigidez.
 *     
 *     Possibilidade de utilizar o solver SAMG em todos os algoritmos
 *     implícitos.
 *     
 *
 */

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

#include "constrain.h"
#include "load.h"
#include "elm.h"
#include "node.h"
#include "alg.h"
#include "xgplib.h"
#include "fem.h"
#include "csr.h"
#include "pfs.h"

#ifdef USE_SAMG
  #include "solversamg.h"
#endif



/*
** ------------------------------------------------------------------------
** Local variables and symbols:
*/
#ifndef maxcor
#define maxcor 30
#endif

#ifndef maxiter
#define maxiter 5000
#endif

#ifndef cgtol
#define cgtol 1.0e-10
#endif

#ifndef ngr
#define ngr 10e20
#endif


/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void ElmKpr(double *, int *, double *);
static void ElmKp(double *, double *, int *, int *, double *);
static void GCSolver(double *, double *, int *, int *, double *);
static void SearchDirection(int i, int k, double *vector, double **w, double **v,
                            int nna, int nnu, int *ia, int *ja, double *avector); 

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

/*
%F
*/
static void SearchDirection(int i, int k, double *vector, double **w, double **v,
                            int nna, int nnu, int *ia, int *ja, double *avector)
{
/*
Essa função calcula a direção de busca para
o método BFGS
*/      
  int n;
  int j;
  int m;
  double a;
  double *b;
  int    matrix = 12;
  int    iversion = 4;
  int    nsys = 1;
  int    npnt = 0;
  
  n=0;
  a=0.0;
  m=k;

  b =  (double*)malloc(i*sizeof(double)); 

  for (j = 0; j < i; j++){
   a += w[j][m]*vector[j];
  }
  for (j = 0; j < i; j++){
   b[j] = vector[j]+v[j][m]*a;
  }
  for (j = 0; j < i; j++){
   vector[j] = b[j];
  }
  while(m>0){
   m -=1;
   a=0.0;
   for (j = 0; j < i; j++){
    a += w[j][m]*vector[j];
   }
   for (j = 0; j < i; j++){
    b[j] = vector[j]+v[j][m]*a;
   }
   for (j = 0; j < i; j++){
    vector[j] = b[j];
   }
  }

  if(Config.solver == 0){

   /* usa o método do gradiente conjugado
   */
   GCSolver(b, vector, ia, ja, avector);
  
  }
  else{
 
   /* chama o solver Samg para calcular o incremento dos deslocamentos*/
   #ifdef USE_SAMG
    SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, b, vector );
   #endif
  }
  
  a=0.0;
  for (j = 0; j < i; j++){
   a += v[j][n]*vector[j];
  }
  for (j = 0; j < i; j++){
   b[j] = vector[j]+w[j][n]*a;
  }
  for (j = 0; j < i; j++){
   vector[j] = b[j];
  }
  while(n<k){
   n +=1;
   a=0.0;
   for (j = 0; j < i; j++){
    a += v[j][n]*vector[j];
   }
   for (j = 0; j < i; j++){
    b[j] = vector[j]+w[j][n]*a;
   }
   for (j = 0; j < i; j++){
    vector[j] = b[j];
   }
  }
  free ( b );
}


/*
%F
*/
static void ElmKpr(double *Kp, int *ia, double *avector)
{
 int     i;           

 for( i = 0; i < NumNodes * NDof; i++ ){         
  Kp[i] = avector[ia[i]];
 }
}

/*
%F
*/
static void ElmKp(double *d, double *p, int *ia, int *ja, double *avector)
{
 int     i,j,k,l,m;           

 for( i = 0; i < NumNodes * NDof; i++ ){  
  j = ia[i];
  l = ia[i+1];  
  for (k = j; k < l; k++){
   m=ja[k]; 
   d[i] += avector[k] * p[m];
  }
 }
}

/*
%F
*/
static void GCSolver(double *R, double *u, int *ia, int *ja, double *avector)
{
 int    i,j,m;                      
 double alpha,beta,a,b;                
 double *rn,*p,*d,*Kp,*z;
 double error;
 int    iter;
 
 m = NumNodes * NDof;
 
 rn =  (double*)malloc(m*sizeof(double));
 p  =  (double*)malloc(m*sizeof(double));
 d  =  (double*)malloc(m*sizeof(double));
 Kp =  (double*)malloc(m*sizeof(double));
 z  =  (double*)malloc(m*sizeof(double));
 
 for( i = 0; i < m; i++ ){
  Kp[i] = 0.0;
 }
 
 /* monta matriz de precondicionamento diagonal
 */
 ElmKpr(Kp,ia,avector);
 
 for( i = 0; i < NumNodes; i++ ){
  if( NodeVector[i].dof.x == DISPLACEMENT ){
   Kp[NDof*i]+= ngr;
  }
  if( NodeVector[i].dof.y == DISPLACEMENT ){
   Kp[NDof*i+1]+= ngr;
  }
  if( NDof == 3 ) {
   if( NodeVector[i].dof.z == DISPLACEMENT ){
    Kp[NDof*i+2]+= ngr;
   }
  }
 }      
 
 for( j = 0; j < m; j++ ){
  if(Kp[j] != 0.0){
   Kp[j] = 1.0/Kp[j];
  }
  else{
   Kp[j] = 1.0;
  }
  rn[j] = R[j];
  z[j] = p[j] = Kp[j]*rn[j];
  d[j] = u[j] = 0.0;
 }
 
 error = 1.0; iter = 0;
 
 while((error>=cgtol)&&(iter<maxiter)){
  iter ++;
  
  error = 0.0;
  
  ElmKp(d,p,ia,ja,avector);
  
  for( i = 0; i < NumNodes; i++ ){
   if( NodeVector[i].dof.x == DISPLACEMENT ){
    d[NDof*i] = ngr*p[NDof*i];
   }
   if( NodeVector[i].dof.y == DISPLACEMENT ){
    d[NDof*i+1] = ngr*p[NDof*i+1];
   }
   if( NDof == 3 ) {
    if( NodeVector[i].dof.z == DISPLACEMENT ){
     d[NDof*i+2] = ngr*p[NDof*i+2];
    }
   }
  }
  
  a = b = 0.0;
  
  for (j = 0; j < m; j++){
   a+= z[j]*rn[j]; 
   b+= p[j]*d[j];
  }
  
  alpha = a/b;
  a = b = 0.0;
  
  for( j = 0; j < m; j++ ){
   u[j]+= alpha*p[j];
   rn[j]= R[j]-alpha*d[j];
   b+= z[j]*R[j];
   z[j]= Kp[j]*rn[j];
   a+= z[j]*rn[j];
  }
  
  beta=a/b;
  
  for( j = 0; j < m; j++ ){
   p[j]= z[j]+beta*p[j];
   R[j]= rn[j]; 
   d[j]= 0.0;
   error+= rn[j]*rn[j];
  }
  
  error=sqrt(error);
  
 }
 
 free ( rn );
 free ( p  );
 free ( d  );
 free ( Kp );
 free ( z  );
}

/*
** ------------------------------------------------------------------------
** Public functions:
*/

/*
%F
*/
void BuildAVector( int *ia, int *ja, double *avector )
{
 int    i,j,k,l,elm;                      
 double ke[24][24];
 int    u[24];
 int    index;
 
 /* monta matriz de rigidez dos elementos
    no formato csr
 */
 for( elm = 0; elm < NumElements; elm++ ){
  /* pega os graus de liberdade do elemento
  */
  ElmGetDof( ElmList[elm], u, &index );
  /* monta matriz de rigidez do elemento
  */
  ElmKMatrix( ElmList[elm], ke );
  /* coloca os termos da matriz do elemento
     no vetor de coeficientes avector
  */
  for(i = 0; i < index; i++){
   for(j = 0; j < index; j++){
    k = u[i];  
    l = u[j];
    BuildaVector( k, l, ke[i][j], ia, ja, avector );
   }
  }
 }
 /* aplica condições de contorno em avector
 */
 for( i = 0; i < NumNodes; i++ ){
  if( NodeVector[i].dof.x == DISPLACEMENT ){
   avector[ia[NDof*i]]+= ngr;
  }
  if( NodeVector[i].dof.y == DISPLACEMENT ){
   avector[ia[NDof*i+1]]+= ngr;
  }
  if( NDof == 3 ) {
   if( NodeVector[i].dof.z == DISPLACEMENT ){
    avector[ia[NDof*i+2]]+= ngr;
   }
  }
 }
}       


/*
%F
*/
void BFGS(double dtime, double *F, double *U, double *V,
          int *iteration, double *error, double error0, int loadstep,
          int nna, int nnu, int *ia, int *ja, double *avector )
{
/*
Essa função busca minimizar o vetor resíduo
com o método BFGS
*/
  int    i,j,k,n;
  double vp[3];
  eRestType dof[3];                     
  double b, c;                                                 
  double *q, *F0;
  double **w, **v;
  int    needite;
  int    stop = 0;
  int    matrix = 12;
  int    iversion = 4;
  int    nsys = 1;
  int    npnt = 0;
  
  n = NumNodes*NDof;
  
  q  = (double*)malloc(n*sizeof(double));
  F0 = (double*)malloc(n*sizeof(double));
  w  = (double**)malloc(n*sizeof(double*));
  v  = (double**)malloc(n*sizeof(double*));
  for(i=0; i < n; i++){
   w[i] = (double*)malloc(maxcor*sizeof(double));
   v[i] = (double*)malloc(maxcor*sizeof(double));
  }
 
  for( i = 0; i < NumNodes; i++ ) {
   dof[0] = NodeVector[i].dof.x;
   dof[1] = NodeVector[i].dof.y;
   dof[2] = NodeVector[i].dof.z;
   vp[0]  = NodeVector[i].dof.vpx;
   vp[1]  = NodeVector[i].dof.vpy;
   vp[2]  = NodeVector[i].dof.vpz;
   for( j = 0; j < NDof; j++ ) {
   /* Verifica a condição de contorno de deslocamento
    */
    if(dof[j] == DISPLACEMENT){
     U[NDof*i+j] = vp[j] * Config.loadfactor;
    }
   }
  }
 
  /* Calcula o vetor de forças desequilibradas
  */
  InternalForces( dtime, F, U, V );

  for( j = 0; j < n; j++ ){
   q[j] = F0[j] = V[j] = F[j];
  }
  k=0; 
 
  /* Calcula o vetor de incrementos de deslocamentos
  */
  if(Config.solver == 0){

   /* usa o método do gradiente conjugado
   */
   GCSolver(F, V, ia, ja, avector);
   
  }
  else{
 
   /* chama o solver Samg para calcular o incremento dos deslocamentos*/
   #ifdef USE_SAMG
    SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, F, V );
   #endif
  }
  
  *error = Config.tolerance;  needite   = 0;
  
  do {
   /* atualiza o vetor de deslocamentos
   */
   for( i = 0; i < NumNodes; i++ ) {
    U[NDof*i]   +=  V[NDof*i];    
    U[NDof*i+1] +=  V[NDof*i+1];  
    if( NDof == 3 ) {
     U[NDof*i+2] += V[NDof*i+2];  
    }
   }
   
   /* Calcula o vetor de forças desequilibradas
   */
   InternalForces( dtime, F, U, V); 
   
   *error = ComputeError( F, V, error0 );
  
   if( _count_it != NULL )
    (*_count_it) ( *iteration, *error, loadstep); 
    
   if( _fstop != NULL )
     (*_fstop) ( &stop ); 
     
   /* Check for another iteration
   */
   needite = NeedIteration( error, iteration, stop );
   
   printf( "\tStep: %d             Error: %5.3E\r", *iteration, *error );
   
   if(!needite) break;
   if(*error < Config.tolerance) break;

   b=c=0.0;
   for (j = 0; j < n; j++){
    q[j] -= F[j];
    b += V[j]*q[j];
    c += V[j]*F0[j];
   }
   for (j = 0; j < n; j++){
    v[j][k]= -(pow(fabs(b/c),0.5))*F0[j]-q[j];   
    w[j][k]= V[j]/b;                            
    V[j] = q[j] = F0[j] = F[j];
   }
   /* calcula a direção de busca*/
   SearchDirection(n, k, V, w, v, nna, nnu, ia, ja, avector);    
   k ++;
   if( k == maxcor ) k = 0;
  }while( needite); 
 
  free ( q );
  free ( F0 );
  for(i=0; i < n; i++){
   free( w[i] );
   free( v[i] );
  }
  free( w ); 
  free( v );    
 
}

/*
%F
*/
void NRM(double dtime, double *F, double *U, double *V,
         int *iteration, double *error, double error0, int loadstep,
         int nna, int nnu, int *ia, int *ja, double *avector )
{
/*
Essa função busca minimizar o vetor resíduo
com o método NRM
*/
  int    i,j,n;
  double vp[3];
  eRestType dof[3];                       
  int    needite = 0;
  int    stop = 0;
  int    matrix = 12;
  int    iversion = 4;
  int    nsys = 1;
  int    npnt = 0;
  
  n = NumNodes*NDof;
  
  do {
        
   for( i = 0; i < NumNodes; i++ ) {
    dof[0] = NodeVector[i].dof.x;
    dof[1] = NodeVector[i].dof.y;
    dof[2] = NodeVector[i].dof.z;
    vp[0]  = NodeVector[i].dof.vpx;
    vp[1]  = NodeVector[i].dof.vpy;
    vp[2]  = NodeVector[i].dof.vpz;
    for( j = 0; j < NDof; j++ ) {
    /* Verifica a condição de contorno de deslocamento
     */
     if(dof[j] == DISPLACEMENT){
      U[NDof*i+j] = vp[j] * Config.loadfactor;
     }
    }
   }
 
   /* Calcula o vetor de forças desequilibradas
   */
   InternalForces( dtime, F, U, V );

   for( j = 0; j < n; j++ ){
    V[j] = F[j];
   }
 
   /* Calcula o vetor de incrementos de deslocamentos
   */
  if(Config.solver == 0){

   /* usa o método do gradiente conjugado
   */
   GCSolver(F, V, ia, ja, avector);
  
  }
  else{
 
   /* chama o solver Samg para calcular o incremento dos deslocamentos*/
   #ifdef USE_SAMG
    SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, F, V );
   #endif
  }
 
   /* atualiza o vetor de deslocamentos
   */
   for( i = 0; i < NumNodes; i++ ) {
    U[NDof*i]   +=  V[NDof*i];    
    U[NDof*i+1] +=  V[NDof*i+1];  
    if( NDof == 3 ) {
     U[NDof*i+2] += V[NDof*i+2];  
    }
   }
   
   /* Calcula o vetor de forças desequilibradas
   */
   InternalForces( dtime, F, U, V); 
   
   *error = ComputeError( F, V, error0 );
  
   if( _count_it != NULL )
    (*_count_it) ( *iteration, *error, loadstep); 
    
   if( _fstop != NULL )
     (*_fstop) ( &stop ); 
     
   /* Check for another iteration
   */
   needite = NeedIteration( error, iteration, stop );
   
   printf( "\tStep: %d             Error: %5.3E\r", *iteration, *error );            
                 
  }while(needite); 
 
}

/*
%F
*/
void LINEAR(double dtime, double *F, double *U, double *V,
            int *iteration, double *error, double error0, int loadstep,
            int nna, int nnu, int *ia, int *ja, double *avector )
{
  int    i,j;
  double vp[3];
  eRestType dof[3];                     
  int    stop = 0;
  int    matrix = 12;
  int    iversion = 4;
  int    nsys = 1;
  int    npnt = 0;
 
  for( i = 0; i < NumNodes; i++ ) {
   dof[0] = NodeVector[i].dof.x;
   dof[1] = NodeVector[i].dof.y;
   dof[2] = NodeVector[i].dof.z;
   vp[0]  = NodeVector[i].dof.vpx;
   vp[1]  = NodeVector[i].dof.vpy;
   vp[2]  = NodeVector[i].dof.vpz;
   for( j = 0; j < NDof; j++ ) {
   /* Verifica a condição de contorno de deslocamento
    */
    if(dof[j] == DISPLACEMENT){
     U[NDof*i+j] = vp[j] * Config.loadfactor;
    }
   }
  }
 
  /* Calcula o vetor de forças desequilibradas
  */
  InternalForces( dtime, F, U, V );
  
  for( j = 0; j < NumNodes*NDof; j++ ){
   V[j] = F[j];
  }
 
  /* Calcula o vetor de incrementos de deslocamentos
  */
  if(Config.solver == 0){

   /* usa o método do gradiente conjugado
   */
   GCSolver(F, V, ia, ja, avector);
   
  }
  else{
  
   /* chama o solver Samg para calcular o incremento dos deslocamentos*/
   #ifdef USE_SAMG
    SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, F, V );
   #endif
  }
 
  /* atualiza o vetor de deslocamentos
  */
  for( i = 0; i < NumNodes; i++ ) {
   U[NDof*i]   +=  V[NDof*i];    
   U[NDof*i+1] +=  V[NDof*i+1];  
   if( NDof == 3 ) {
    U[NDof*i+2] += V[NDof*i+2];  
   }
  }
  
  /* Calcula o vetor de forças desequilibradas
   */
  InternalForces( dtime, F, U, V); 
   
  *error = ComputeError( F, V, error0 );
  
  if( _count_it != NULL )
    (*_count_it) ( *iteration, *error, loadstep); 
    
  if( _fstop != NULL )
     (*_fstop) ( &stop );

}

/*
%F
*/
void InternalForces( double dtime, double *F, double *U, double *V )
{
 int     i;
 sTensor stress[8], strain[8];
 double  iforce[24];

/* Initialize the internal forces
 */
 for( i = 0; i < NumNodes; i++ ) {
  if( NodeVector[i].rezone == NONE ) {
   if( NodeVector[i].dof.x == FORCE ) F[NDof*i] = NodeVector[i].dof.vpx * Config.loadfactor;
   else                               F[NDof*i] = 0.0;

   if( NodeVector[i].dof.y == FORCE ) F[NDof*i+1] = NodeVector[i].dof.vpy * Config.loadfactor;
   else                               F[NDof*i+1] = 0.0;

   if( NDof == 3 ) {
    if( NodeVector[i].dof.z == FORCE ) F[NDof*i+2] = NodeVector[i].dof.vpz * Config.loadfactor;
    else                               F[NDof*i+2] = 0.0;
   }
  }
 }

/* Loop over the elements
 */
 for( i = 0; i < NumElements; i++ ) {
 /* Check to see if the element was not rezoned
  */
  if( ElmList[i]->rezone != NONE )
   continue;

 /* Compute element stress strain
  */
  ElmStressStrain( ElmList[i], dtime, U, 0L, stress, strain );

 /* Correct stresses for the larger displacements case
  */
  if( IsGeomNonLinear && (ElmList[i]->rezone != DONE) )
   ElmUpdateStress( ElmList[i], dtime, V, stress );

 /* Compute element internal forces
  */
  ElmInterForce( ElmList[i], stress, iforce );

 /* Assemble in the global vector
  */
  ElmAssVector( ElmList[i], F, iforce );
 }

} /* End of InternalForces */

/*
%F
*/
int DisplacementVelocity( int iteration, double alpha, double dtime,
                                         double *FVector, double *VVector,
                                         double *UVector, double *MVector )
{
 int         i, j;
 double      a, b, df, disp, veloc;
 double      vp[3];
 eRestType   dof[3];
 
 /* Cumpute the auxiliares constantes
  */
 a = 2.0 - (alpha * dtime);
 b = 2.0 + (alpha * dtime);
 /*teste adicionado por Gisele*/
 pfs_iter = iteration;

 /* Loop over the nodes
  */
 for( i = 0; i < NumNodes; i++ ) 
 {
  pfs_node = i;
  /* Check to see if the node was not rezoned
   */
  if( NodeVector[i].rezone != NONE )
   continue;

  /* Fill the auxiliare vector for the prescribed value
   */
  dof[0] = NodeVector[i].dof.x;
  dof[1] = NodeVector[i].dof.y;
  dof[2] = NodeVector[i].dof.z;
  vp[0]  = NodeVector[i].dof.vpx;
  vp[1]  = NodeVector[i].dof.vpy;
  vp[2]  = NodeVector[i].dof.vpz;

  /* Inicialize total displacement and total velocity
   */ 
  disp  = 0.0;
  veloc = 0.0;

  /* case: node with any constrain */  
  if((dof[0] == OBJECT) || (dof[1] == OBJECT) || (dof[2] == OBJECT))
  {
    /* id node constrain */
    int idNodeConst = NodeVector[i].constrain.id - 1;
    /* pointer to id element constrain */
    void* ConstrElm = &(NodeVector[i].constrain.curr_elm);

    double fl[3];
    double vl[3];
    double ml = MVector[NDof*i];

    sCoord p;

    if( iteration == 0 ) 
    {
      if( !ConstToLocal(ConstList[idNodeConst], &NodeVector[i].coord,
                        ConstrElm, vp, fl) ) 
      {
        return 0;
      } 

      for(j = 0; j < NDof; j++) 
        vl[j] = 0.01 * dtime * (fl[j] / (2.0 * ml));

      if( !ConstToGlobal(ConstList[idNodeConst], &NodeVector[i].coord,
                         ConstrElm, vl, &VVector[NDof*i]) )
        return 0;

      /* Usado somente para a implementação ESPAÇO PARAMÉTRICO */
      for(j = 0; j < NDof; j++) 
        UVector[NDof*i+j] += dtime * VVector[NDof*i+j];
    }
    else
    {
      /* O ponto p calculado abaixo tem que estar sobre a curva -> ConstList.
       * Se nao estiver nao sera possivel: ConstToLocal e ConstToGlobal.
       */
      p.x = NodeVector[i].coord.x + UVector[NDof*i];
      p.y = NodeVector[i].coord.y + UVector[NDof*i+1];
      if(NDof == 3) 
        p.z = NodeVector[i].coord.z + UVector[NDof*i+2];
      else
        p.z = 0.0;
      if( !ConstToLocal(ConstList[idNodeConst], &p, ConstrElm,
                        &FVector[NDof*i], fl) ) return 0;
      if( !ConstToLocal(ConstList[idNodeConst], &p, ConstrElm,
                        &VVector[NDof*i], vl) ) return 0;
      /* Integração.
       */
      for(j = 0; j < NDof; j++) 
        vl[j] = ((a / b) * vl[j]) + (2.0 * dtime * (fl[j] / (b * ml)));

      /* Projetando as grandezas no elemento corrente da restricao
       */
      fl[2] = vl[2] = 0.0;
      if(NDof <= 2)
        fl[1] = vl[1] = 0.0;

      {
        double UVecPrev[3];   /* UVector no passo anterior.             */
        sCoord pl;            /* Ponto sobre a restrição.               */
        sCoord delta_u;       /* deslocamento no elemento da restrição. */

        /* Calcula o deslocamento ("delta") no espaço paramétrico.
         */
        delta_u.x = dtime * vl[0];

        if(NDof == 3)
          delta_u.y = dtime * vl[1];
        else
          delta_u.y = 0.0;       /* dtime * vl[1] */

        delta_u.z = 0.0;         /* dtime * vl[2] */
           
        if( !ConstSlide(ConstList[idNodeConst], &p, ConstrElm, 
                        &delta_u, vl, &VVector[NDof*i], &pl) ) 
        {
         /* ERROR */
          printf( "DisplacementVelocity: Node didn't slide\n" );
         
        }

        /* Guarda o UVector antes de atualizar ele para atualizar o 
         * VVector.
         */
        for(j = 0; j < NDof; j++) 
          UVecPrev[j] = UVector[NDof*i+j];

        /* Calcula o UVector. 
         */
        UVector[NDof*i]   = pl.x - NodeVector[i].coord.x;
        UVector[NDof*i+1] = pl.y - NodeVector[i].coord.y; 
        if(NDof == 3) 
          UVector[NDof*i+2] = pl.z - NodeVector[i].coord.z; 

        /* Calcula o VVector.
         */    
        for(j = 0; j < NDof; j++) 
          VVector[NDof*i+j] = (UVector[NDof*i+j] - UVecPrev[j]) / dtime;
      }
    }
  }

 /* Loop over the D.O.F. for each node
  */
  for( j = 0; j < NDof; j++ ) {
  /* Check the D.O.F. condition
   */
   if(MVector[NDof*i+j] != 0.0){
           
    switch( dof[j] ) {
     case FORCE:
      if( iteration == 0 )
       df = vp[j];
      else
       df = FVector[NDof*i+j];
    
      if( iteration == 0 )
        VVector[NDof*i+j] = 0.01 * dtime * (df / (2.0 * MVector[NDof*i+j]));
      else
        VVector[NDof*i+j] = ((a / b) * VVector[NDof*i+j]) +
                           (2.0 * dtime * (df / (b * MVector[NDof*i+j])));
    
      UVector[NDof*i+j] += dtime * VVector[NDof*i+j];
      break;
    
     case DISPLACEMENT:
      VVector[NDof*i+j] = 0.0;
      UVector[NDof*i+j] = vp[j] * Config.loadfactor;
      /* Com o nó piloto. */
      if(iteration == 0)
      {
        /* O nó tem que estar sobre uma das restrições.
         */
        if(NodeVector[i].constrain.id) 
        {
          /* O nó não pode ser o nó piloto.
           */
          if(!NodeVector[i].pilot) 
          {
            NodeVector[i].dof.x = OBJECT;
            NodeVector[i].dof.y = OBJECT;
            NodeVector[i].dof.z = OBJECT;
          }
        }
      }
      break;
    
     case VELOCITY:
      VVector[NDof*i+j]  = vp[j];
      UVector[NDof*i+j] += dtime * VVector[NDof*i+j];
      break;
    
     default:
      break;
    }
   }

   if( (NodeVector[i].dspIteration == DSP_YES) && (GpType == 0) )
    disp += (UVector[NDof*i+j] * UVector[NDof*i+j]);
   else if( (NodeVector[i].dspIteration == DSP_YES) && (GpType == 1) )
    veloc += (VVector[NDof*i+j] * VVector[NDof*i+j]);
  }

  if( (NodeVector[i].dspIteration == DSP_YES) && (GpType == 0) )
   XGPSendPoint( NodeVector[i].curve, (double)iteration, sqrt( disp ) );
  else if( (NodeVector[i].dspIteration == DSP_YES) && (GpType == 1) )
   XGPSendPoint( NodeVector[i].curve, (double)iteration, sqrt( veloc ) );
 }
   return 1;
} /* End of DisplacementVelocity */


/*
%F This function computes the mass vector for the model.
*/
void MassVector( double *MVector )
{
 int    i;
 double mass[24];

/* Initialize the mass vector
 */
 for( i = 0; i < NDof*NumNodes; i++ )
  MVector[i] = 0.0;

/* Loop over the elements
 */
 for( i = 0; i < NumElements; i++ ) {
 /* Check to see if the element was rezoned
  */
  if( ElmList[i]->rezone == DONE )
   continue;

 /* Compute the element mass matrix
  */
  ElmMass( ElmList[i], mass );

 /* Assemble in to the global vector
  */
  ElmAssVector( ElmList[i], MVector, mass );
 }

} /* End of MassVector */


/*
%F
*/
int DoRezone( int step, double *FVector )
{
 int i, j, id, numnodes;
 int inc[8];

/* Check to see if the rezone descriptor exists
 */
 if( Rezone == 0L )
  return( 0 );

/* Update element rezone status for the current step
 */
 for( i = 0; i < Rezone[step-1].numelm; i++ ) {
  id = Rezone[step-1].elm[i] - 1;
  if( ElmList[id]->rezone == NONE ) ElmList[id]->rezone = TODO;
  else                              ElmList[id]->rezone = DONE;
 }

/* Update node flag
 */
 for( i = 0; i < NumNodes; i++ )
  NodeVector[i].rezone = DONE;

 for( i = 0; i < NumElements; i++ ) {
  if( ElmList[i]->rezone == NONE ) {
   ElmConnect( ElmList[i], inc );
   ElmNumNodes( ElmList[i], &numnodes );
   for( j = 0; j < numnodes; j++ )
    NodeVector[inc[j]-1].rezone = NONE;
  }
 }

/* Update prescribed values
 */
 UpdatePrescribedValues( FVector );

/* Set element rezone state to DONE
 */
 for( i = 0; i < NumElements; i++ )
  if( ElmList[i]->rezone == TODO )
   ElmList[i]->rezone = DONE;

 return( 1 );

} /* End of DoRezone */


/*
%F
*/
void PercolationForces( double *F )
{
 int    i;
 double pforce[24];

/* Loop over the elements
 */
 for( i = 0; i < NumElements; i++ ) {
 /* Check to see if the element was not rezoned
  */
  if( ElmList[i]->rezone != NONE )
   continue;

 /* Compute element percolation forces
  */
  ElmPercForce( ElmList[i], pforce );

 /* Assemble in the global vector
  */
  ElmAssVector( ElmList[i], F, pforce );
 }

} /* End of PercolationForces */


/* =======================================================  End of File  */

---