_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"
#include "solversamg.h"


/*
** ------------------------------------------------------------------------
** 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 *);
static void ElmKp(double *, double *);
static void GCSolver(double *, double *);
static void SearchDirection(int i, int k, double *vector, double **w, double **v); 

static void ElmKpParalelo(int , int , int , int , double *, double *,
                          double *, double *, double *, double *);
static void ElmKprParalelo(int , int , int , int, double *,
                           double *, double *, double *, double *,
                           double *, double *);
static void ContourParalelo(int , int , double *, double *);
static void ContourKpParalelo(int , int , double *); 
static void ProdEscalarParalelo(int, int , double *, double *,
                                double *, double *, double *, double *);
static void UpDataParalelo(int , int , double *, double *, double ,
                           double *, double *, double *, double *,
                           double *, double *, double *);
static void ResParalelo(int , int , double ,
                        double *, double *, double *,
                        double *, double *, double *); 
static void PreParalelo(int , int , double *, double *,
                        double *, double *, double *,
                        double *, double *);                                                                                                            

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

/*
%F
*/
static void SearchDirection(int i, int k, double *vector, double **w, double **v)
{
/*
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);
  
  }
  else{
 
   /* chama o solver Samg para calcular o incremento dos deslocamentos
   */
   SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, b, vector );
   
  }
  
  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 );
}

/*INICIO TESTE*/

static void ProdutoEscalar(int n, int nn, double *w, double *v, double *a )
{
 int  i;
 
 *a = 0.0;
 
 for (i = n; i < nn; i++){
  *a += w[i]*v[i];
 }
}

static void SomaVector(int n, int nn, double *w, double *v, double *b, double a)
{
 int  i;
 
 for (i = n; i < nn; i++){
  b[i] = w[i] + v[i]*a;
  v[i] = b[i];
 }
}

static void SearchDirectionParalelo(int i, int k, double *vector, double **w, double **v,
                                    int ni[], int nf[], int gli[], int glf[])
{
/*
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;
  
  int    nd = 10;
  int    l;
  double ea;
  double *wi, *vi;
  
  n=0;
  a=0.0;
  m=k;

  b  = (double*)malloc(i*sizeof(double)); 
  wi = (double*)malloc(i*sizeof(double));
  vi = (double*)malloc(i*sizeof(double));
  
  for(j = 0; j < i; j++){
   wi[j] = w[j][m];
   vi[j] = v[j][m];
  } 
  
  ea = 0.0;
  
  for (l = 0; l < nd; l++){
   
   ProdutoEscalar(gli[l], glf[l], wi, vector, &a);
   
   ea += a;
  }
  
  for (l = 0; l < nd; l++){
   
   SomaVector(gli[l], glf[l], vector, vi, b, ea);
   
  }

  while(m>0){
   m -=1;
   ea = 0.0;
   
   for(j = 0; j < i; j++){
    wi[j] = w[j][m];
    vi[j] = v[j][m];
   }
   
   for (l = 0; l < nd; l++){
   
    ProdutoEscalar(gli[l], glf[l], wi, vector, &a);
   
    ea += a;
   }
   
   for (l = 0; l < nd; l++){
   
    SomaVector(gli[l], glf[l], vector, vi, b, ea);
   
   }
   
  }

  if(Config.solver == 0){

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

/*FIM TESTE*/

#if 0

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

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

/*
%F
*/
static void ElmKp(double *d, double *p)
{
 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    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);
 
 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);
  
  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  );
}

#else
/* TESTE PARA PARALELISMO*/                           
                                                      
/*                                                    
%F                                                    
*/                                                    
static void ElmKprParalelo(int n, int nn, int m, int mm, double *Kp,
                           double *d, double *p, double *z, double *rn,
                           double *u, double *R)                        
{                                                     
 int     i;                                           
                                                      
 for( i = n; i < nn; i++ ){           
  Kp[i] = avector[ia[i]];                             
 }                 
 
 ContourKpParalelo(m, mm, Kp); 
 
 PreParalelo(n, nn, Kp, d, p, z, rn, u, R); 
                                    
}                                                     
                                                      
/*                                                    
%F                                                    
*/                                                    
static void ElmKpParalelo(int n, int nn, int ni, int nf, double *d, double *p,
                          double *z, double *rn, double *a, double *b)                          
{                                                     
 int     i,j,k,l,m;                                   
                                                      
 for( i = n; i < nn; i++ ){                            
  j = ia[i];                                      
  l = ia[i+1];                                        
  for (k = j; k < l; k++){                            
   m=ja[k];                                           
   d[i] += avector[k] * p[m];                         
  }                                                   
 }         
     
 ContourParalelo(ni, nf, d, p);
   
 ProdEscalarParalelo(n, nn, a, b, d, p, z, rn);  
                                      
}          

static void ContourParalelo(int n, int nn, double *d, double *p) 
{
 int i;
 
 for( i = n; i < nn; 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];                  
   }                                                
  } 
 }  
   
}  

static void ContourKpParalelo(int n, int nn, double *Kp) 
{
 int i;
 
 for( i = n; i < nn; 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;                               
   }                                                  
  }                                                   
 }
   
}    

static void ProdEscalarParalelo(int n, int nn, double *a, double *b,
                                double *d, double *p, double *z, double *rn) 
{
 int i;
 
 *a = *b = 0.0;  

 for (i = n; i < nn; i++){                           
  *a += z[i]*rn[i];                                   
  *b += p[i]*d[i];                                    
 }
 
}   

static void UpDataParalelo(int n, int nn, double *a, double *b, double alpha,
                           double *u, double *d, double *p, double *z,
                           double *rn, double *R, double *Kp) 
                             
{
 int i;
 
 *a = *b = 0.0;  
 
 for (i = n; i < nn; i++){                           
  u[i]+= alpha*p[i];                               
  rn[i]= R[i]-alpha*d[i];                          
  *b += z[i]*R[i];                                   
  z[i]= Kp[i]*rn[i];                               
  *a += z[i]*rn[i];                                   
 }
 
}

static void ResParalelo(int n, int nn, double beta,
                        double *d, double *p, double *z,
                        double *rn, double *R, double *error) 
{
 int i;
 
 for( i = n; i < nn; i++ ){                         
  p[i]= z[i]+beta*p[i];                            
  R[i]= rn[i];                                     
  d[i]= 0.0;                                       
  *error+= rn[i]*rn[i];                             
 } 
 
}

static void PreParalelo(int n, int nn, double *Kp, double *d,
                        double *p, double *z, double *rn,
                        double *u, double *R) 
{
 int i;
 
 for( i = n; i < nn; i++ ){                                                    
  if(Kp[i] != 0.0){                                   
   Kp[i] = 1.0/Kp[i];                                 
  }                                                   
  else{                                               
   Kp[i] = 1.0;                                       
  }                                                   
  rn[i] = R[i];                                       
  z[i] = p[i] = Kp[i]*rn[i];                          
  d[i] = u[i] = 0.0;                                  
 }
 
}    

                       
                                                      
/*                                                    
%F                                                    
*/                                                    
static void GCSolver(double *R, double *u)            
{                                                     
 int    i,m;                                        
 double alpha,beta;                               
 double *rn,*p,*d,*Kp,*z;                             
 double error;                                        
 int    iter;                                         
 int    gli[10], glf[10];
 int    ni[10], nf[10];
 double ea, eb;
 double a, b;
 int    nno;
 
 /*define número de divisões*/
 int    nd = 10;
 
 /*define nó inicial e nó final de cada divisão*/ 
 nno = NumNodes/nd;
                                                      
 ni[0] = 0;
 nf[0] = nno;
 for( i = 1; i < nd-1; i++ ){  
  ni[i] = nf[i-1] + 1;                          
  nf[i] = nf[i-1] + nno;                                        
 }
 ni[nd-1] = nf[nd-2] + 1;
 nf[nd-1] = NumNodes - 1;
 
 /*define dof inicial e dof final de cada divisão*/
 gli[0] = 0;
 glf[0] = nf[0] * NDof;
 for( i = 1; i < nd-1; i++ ){  
  gli[i] = glf[i-1];                          
  glf[i] = nf[i] * NDof;                                        
 }
 gli[nd-1] = glf[nd-2];
 glf[nd-1] = nf[nd-1] * NDof + NDof;
 
 /* dof total */                                                
 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));             
                                                                                                         
 /* monta matriz de precondicionamento diagonal       
 */      
 for(i = 0; i < nd; i++){
                
  ElmKprParalelo(gli[i], glf[i], ni[i], nf[i], Kp, d, p, z, rn, u, R);
  
 }                                             
                                                                                                                                                                                                                                                      
                                                      
 error = 1.0; iter = 0;                               
                                                      
 while((error>=cgtol)&&(iter<maxiter)){
                       
  iter ++;                                            
                                                      
  error = ea = eb = 0.0;                                   
  
  /*divide o processo em nd*/
  for(i = 0; i < nd; i++){
                
   ElmKpParalelo(gli[i], glf[i], ni[i], nf[i], d, p, z, rn, &a, &b);
   
   ea += a;
   eb += b; 
   
  }  
  
  alpha = ea/eb;
  ea = eb = 0.0; 
  
  /*divide o processo em nd*/
  for(i = 0; i < nd; i++){ 
        
   UpDataParalelo(gli[i], glf[i], &a, &b, alpha, u, d, p, z, rn, R, Kp);                                            
   
   ea += a;
   eb += b;
   
  }
  
  beta = ea/eb;                                         
  
  /*divide o processo em nd*/
  for(i = 0; i < nd; i++){ 
        
   ResParalelo(gli[i], glf[i], beta, d, p, z, rn, R, &error);
                                                                                                
  }       
  error = sqrt(error);                         
 }                                                    
                                                      
 free ( rn );                                         
 free ( p  );                                         
 free ( d  );                                         
 free ( Kp );                                         
 free ( z  );                                         
}                                                     
                                                      
/* FIM DO TESTE*/
#endif

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

/*
%F
*/
void BuildAVector( void )
{
 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]);
   }
  }
 }
 /* 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 )
{
/*
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);
   
  }
  else{
 
   /* chama o solver Samg para calcular o incremento dos deslocamentos*/
   SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, F, V );
  }
  
  *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);    
   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 )
{
/*
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);
  
  }
  else{
 
   /* chama o solver Samg para calcular o incremento dos deslocamentos*/
   SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, F, V );
   
  }
 
   /* 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    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);
   
  }
  else{
  
   /* chama o solver Samg para calcular o incremento dos deslocamentos*/
   SamgSolver( iversion, nna, nnu, matrix, nsys, npnt, ia, ja, avector, F, V );
  }
 
  /* 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 );

}


/*INICIO TESTE*/
#if 0
/*
%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 */

#else

/*
%F
*/
void InternalForces( double dtime, double *F, double *U, double *V )
{
 int     i,j;
 sTensor stress[8], strain[8];
 double  iforce[24];
 /*número de divisões*/
 int nd = 10;
 int ei[10], ef[10];
 int nelm;
 double *Faux;
 
 nelm = NumElements/nd;
 
 /*elemento inicial e elemento final de cada divisão*/
 ei[0] = 0;
 ef[0] = nelm;
 
 for(i = 1; i < (nd-1); i++){
  ei[i] = ef[i-1];
  ef[i] = ei[i] + nelm;
 }
 ei[nd-1] = ef[nd-2];
 ef[nd-1] = NumElements;
 
 /*vetor auxiliar. pega o resíduo para cada divisão*/
 Faux  = (double *)calloc( NDof*NumNodes, sizeof(double) );

/* 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;
   }
  }
 }
 
 for( j = 0; j < nd; j++ ) { 
 
  for( i = 0; i < NDof*NumNodes; i++ ) {
   Faux[i] = 0.0;
  }
 
  for( i = ei[j]; i < ef[j]; i++ ) {
 /* Loop over the elements
  */
  /* 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], Faux, iforce );
  }
  
  for( i = 0; i < NDof*NumNodes; i++ ) {
   F[i] += Faux[i];
  }
 
 }
 
 free(Faux);
 

} /* End of InternalForces */

#endif

/*
%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  */

---