alg.c

Go to the documentation of this file.

/*
%M This module contains the funtions to manipulate the solution algorithm
   for the Relax program.
%a Joao Luiz Elias Campos.
%d September 3rd, 1998.
%r $Id: alg.c,v 1.1 2004/06/22 05:29:58 joaoluiz Exp $
%w (C) COPYRIGHT 1995-1996, Eduardo Nobre Lajes.
   (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:    08-Aug-05    Alexandre A. Del Savio
 *     Modificada a funcao DRSolver, onde adicionou-se "_count_it".
 *
 *   Modified:    13-Sep-05    Alexandre A. Del Savio
 *     Criada a funcao AlgReset.
 *
 *   Modified:    08-Dec-05    Juan Pablo Ibañez
 *     Criadas as variaveis locais vetError e delta para controle de 
 *     convergencia.
 *     Criada a funcão ConvergeError, para controle de convergencia.
 *     Modificada a funcão NeedIteration, que carrega o vetor de 
 *     forcas desequilibradas e chama a funcão ConvergeError.
 *     Modificada a funcão DRSolver, onde incluiu-se a alocacão
 *     dinamica de vetError, assim como sua liberacão.
 *     Acrescentou-se o parâmetro delta na funcão _count_it.
 *
 *   Modified:    29-Mar-06    Juan Pablo Ibañez
 *     Modificada a funcão DRSolver, onde a variável delta é zerada no 
 *     final do processo de calculo. Também é chamada a funcão 
 *     ConstrainFree para liberar memória da lista ConstList.
 *     
 *   Modified:    Dez-06       André Luis Muller
 *     Modificação na função DRSolver, possibilitando uma primeira
 *     solução, linear elástica.
 *
 *   Modified:    07       André Luis Muller
 *     Modificação na função DRSolver, possibilitando uma 
 *     solução híbrida.
 *
 *     Novo critério de convergência.
 *
 *     Passos de carga.
 *
 *     Adição do um novo item para verificar parada da análise
 *     função (stop).
 */

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

#include "load.h"
#include "elm.h"
#include "node.h"
#include "fem.h"
#include "alg.h"
#include "nfi.h"
#include "rio.h"
#include "drv.h"
#include "damp.h"
#include "constrain.h"
#include "csr.h"

/*
** ------------------------------------------------------------------------
**definitions of minimum and maximus
*/
#ifndef MAX
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#endif

#ifndef MIN
#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#endif

/*
%V Configuration parameters
*/
sConfig Config = { 0, 0, 0, 1, 0, 0, 0.0, 0, 0.0, 0.0, 0.0, 0 }; 


/*
%V Flag to see if is a coupled problem
*/
eCoupState State = UNCOUPLED;


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

#ifndef PI
#define PI 3.141592654
#endif

/* local variables for convergence control */
double *vetError;   
double delta = 0.0;  

/*
** ------------------------------------------------------------------------
** Local functions:
*/
static double _TimeIncrement ( void );
static void   _UpdateGeometry( double, double * );
static int    _ConvergeError ( );
static double InitialError( double, double *, double *, double *, double *);
/*
%F
*/
static double _TimeIncrement( void )
{
 double dt  = 1.0e6;
 double edt = 1.0e6;
 int    i;

/* For each element computes the time step
 */
 for( i = 0; i < NumElements; i++ )
 {
 /* Check to see if the element was rezoned
  */
  if( ElmList[i]->rezone != NONE )
   continue;

 /* Compute the time step for the element
  */
  ElmTimeStep( ElmList[i], &edt );

 /* Get the minimum time step
  */
  if( edt < dt )
   dt = edt;
 }

 return( dt );

} /* End of _TimeIncrement */


static int _ConvergeError ( )
{
 int i;
 int flag = 1;
 int dimVetError = Config.numsteptol;   
 
 for(i = 0; i < dimVetError; i++)
 {
   if(vetError[i] >= Config.tolerance)
    flag = 0;   
 }

 return( flag );
        
} /* End of _ConvergeError */


/*
%F This function updates the model geometry for the larger displaments case.
*/
static void _UpdateGeometry( double dtime, double *VVector )
{
 int i;

/* Update geometry
 */
 for( i = 0; i < NumNodes; i++ )
 {
  NodeVector[i].coord.x += dtime * VVector[NDof*i];
  NodeVector[i].coord.y += dtime * VVector[NDof*i+1];
  if( NDof == 3 )
   NodeVector[i].coord.z += dtime * VVector[NDof*i+2];
 }

} /* End of _UpdateGeometry */


static double InitialError( double dtime, double *F, double *U, double *V, double *U0)
{
 int       i, j;
 double    error;
 double vp[3];
 eRestType dof[3]; 
 
 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;
   }
  }
 }
 
 /* Compute the internal forces
 */
 InternalForces( dtime, F, U, V );
 
/* Initialize error
 */
 error = 0.0;

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

 /* Loop over the degree of freedown
  */
  for( j = 0; j < NDof; j++ )
  {
   /* Compute the error for this dof
    */
   error += fabs(F[NDof*i+j]);
  }
 }
 
 for( i = 0; i < NumNodes; i++ ) {
  for( j = 0; j < NDof; j++ ) {
   U[NDof*i+j] = U0[NDof*i+j];
  }
 }
 
 return (error);
}

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


/*
%F
*/
double ComputeError( double *FVector, double *UVector, double error0 )
{
 int       i, j;
 double    err;
 
/* Initialize error
 */
 err = 0.0;

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

 /* Loop over the degree of freedown
  */
  for( j = 0; j < NDof; j++ )
  {
   /* Compute the error for this dof
    */
   err += fabs(FVector[NDof*i+j] * UVector[NDof*i+j]);
  }
 }

 return( err/error0 );

} /* End of ComputeError */

/*
%F
*/
void _dUVector( double *UVector, double *UVector0, double *dUVector )
{
 int       i, j;

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

 /* Loop over the degree of freedown
  */
  for( j = 0; j < NDof; j++ )
  {
   /* Compute the dUVector and UVector0
    */
   dUVector[NDof*i+j] = UVector[NDof*i+j] - UVector0[NDof*i+j]; 
   UVector0[NDof*i+j] = UVector[NDof*i+j];
  }
 }

} /* End of _dUVector */

/*
%F
*/
int NeedIteration( double *error, int *iteration, int stop )
{
 int i;
 int flag;
 int dimVetError = Config.numsteptol;
    
 /* load error vector (iteration < dimVetError) */
 if ( (*iteration < dimVetError)  && (stop == 0) )
 {
   vetError[*iteration] = *error;
   flag = 1;
   (*iteration)++; 
   return( flag );
 }  
    
 /* load error vector (iteration > dimVetError) */
 for(i = 0; i < dimVetError-1; i++)
 {
   vetError[i] = vetError[i+1];
 } 
 vetError[dimVetError-1] = *error;

 /* convergence control */
 if( (_ConvergeError( ) != 0) || ((*iteration + 1) > Config.max_iter) ||
     (stop == 1)) 
 {
  flag = 0;
 }
 else
 {
  flag = 1;
  (*iteration)++;
 }

 return( flag );

} /* End of NeedIteration */

/*
%F
*/
int DRSolver(UI_State *R, double *FVector, double *UVector, double *VVector,
                          double *MVector, int loadstep)
{
 int    iteration, needite;
 double dtime, alpha, CinEng, error;
 double *dUVector, *UVector0;
 double error0;
 int    stop = 0;
 unsigned long int i_time = time( NULL );

/* Initialize variables
 */
 iteration = 0;
 needite   = 0;
 dtime     = 0.0;
 alpha     = 0.0;
 error     = 0.0;
 CinEng    = 0.0;
 error0    = 1.0;
 
 dUVector  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 UVector0  = (double *)calloc( NDof*NumNodes, sizeof(double) );

 if(R)
  UISetDisplacement(R, NDof, UVector);

/* Compute the time increment according Cundall (1989)
 */
 dtime  = Config.dtfrac * _TimeIncrement();

/* Mount the mass vector
 */
 MassVector( MVector );
 
/* allocate error vector */ 
 vetError = (double*)malloc(Config.numsteptol*sizeof(double));

 /* Faz analise para cada passo de carga
  */

 iteration = 0;
 
 Config.loadfactor = (loadstep + 1.0) / Config.num_load_step;
 
 /* Compute the initial error
 */
 error0 = InitialError(dtime, FVector, UVector, VVector, UVector0);
 
 /* Process iteration until convergence or maximum iteration number
 */
 do {
  
  /* Compute the rigidity coefficient
  */
  DampCalcCoeff(DampObj,iteration,dtime,MVector,VVector,&alpha,&CinEng);
  
  /* Compute displacements and veclocities
  */
  if( !DisplacementVelocity( iteration, alpha, dtime, 
                            FVector, VVector, UVector, MVector ) ) 
   return 0;
  /*
  */
  if( R && (iteration%Config.draw_step) )
   UIDraw(R);
  
  /* Update geometry for the larger displacements case
  */
  if( IsGeomNonLinear )
   _UpdateGeometry( dtime, VVector );
  
  /* Compute the percolation forces
  */
  if( IsCoupled )
   PercolationForces( FVector );
  
  /* Compute the internal forces
  */
  InternalForces( dtime, FVector, UVector, VVector );
  
  /* Compute the dUVector
  */
  _dUVector(UVector, UVector0, dUVector);
   
  /* Compute the error
  */
  error = ComputeError( FVector , dUVector, error0);
  
  /* For larger displaments correct the time step and compute a new
  * mass vector for the model
  */
  if( IsGeomNonLinear ) {
   if( error > 1.0 ) {
    dtime = Config.dtfrac * _TimeIncrement();
    MassVector( MVector );
   }
  }
  
  if( _count_it != NULL )
    (*_count_it) ( iteration, error, loadstep);
    
  if( _fstop != NULL )
    (*_fstop) ( &stop );  
    
  /* Check for another iteration
  */
  needite = NeedIteration( &error, &iteration, stop );    
  
  fprintf( stdout, "\tStep: %d\tError: %3.3e\r", iteration+1, error);

  fflush( stdout );
  
  if(Config.num_load_step == 1){
   /* Save results on a temporary file
   */
   DrvPrintResult( DrvObj, iteration, UVector, VVector );
  }
 } while( needite );

 if(Config.num_load_step > 1){
  /* Save results on a temporary file
  */
  DrvPrintResult( DrvObj, loadstep, UVector, VVector );
 }
 

 /* free error vector */
 free ( vetError );

 fprintf( stdout, "\n\tTotal solver time....................: %lu segs\n\n", (time(NULL)-i_time) );

 free( dUVector );
 free( UVector0 );
 
 return( 1 );

} /* End of DRSolver */

/*
%F
*/
int IMPLINEARSolver(UI_State *R, double *FVector, double *UVector, double *VVector,
                    double *MVector, int loadstep)
{
 int    iteration, needite;
 double dtime, alpha, CinEng, error;
 double *dUVector, *UVector0;
 double error0;
 int    stop = 0;
 unsigned long int i_time = time( NULL );
 int i;
 int    *ia;
 int    *ja = NULL;
 double *avector;
 int    nna = 0;
 int    nnu = 0;
 CsrNode   *node_list = 0L;      /* lista de nós Csr*/

/* Initialize variables
 */
 iteration = 0;
 needite   = 0;
 dtime     = 0.0;
 alpha     = 0.0;
 error     = 0.0;
 CinEng    = 0.0;
 error0    = 1.0;
 Config.numsteptol    = 1;
 
 dUVector  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 UVector0  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 
 nnu = NDof*NumNodes;

 /* aloca memória para lista de nós
 */
 node_list = (CsrNode*) calloc(NumNodes + 1, sizeof(CsrNode));
 
 /* constroi a lista de nós
 */
 BuildNodeList(node_list);

 /* constroi a lista de nós
 */
 BuildNodeList(node_list);
 
 /* aloca memória para ia*/
 ia  = (int*)calloc( nnu + 1, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildIaVec( node_list, ia, nnu, &nna );

 /* aloca memória para ia*/
 avector  = (double*)calloc( nna, sizeof(double) );
 
 /* aloca memória para ja*/
 ja  = (int*)calloc( nna, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildJaVec( node_list, ja );
 
 /* libera memória da lista de nós
 */
 FreeNodeList(node_list);
 
 /* cria o vetor avector*/
 BuildAVector( ia, ja, avector );
 
 if(Config.solver == 1){
 
 /* Mudança de índices necessárias
    para utilização do solver SAMG*/
  for( i = 0; i < nnu+1; i++ ) {
   ia[i] += 1;
  }
  for( i = 0; i < nna; i++ ) {
   ja[i] += 1;
  }
  
 }
 
 if(R)
  UISetDisplacement(R, NDof, UVector);

/* Compute the time increment according Cundall (1989)
 */
 dtime  = Config.dtfrac * _TimeIncrement();

/* Mount the mass vector
 */
 MassVector( MVector );
 
/* allocate error vector */ 
 vetError = (double*)malloc(Config.numsteptol*sizeof(double));

 /* Faz analise para cada passo de carga
  */

 iteration = 0;
 
 Config.loadfactor = (loadstep + 1.0) / Config.num_load_step;
 
 /* Compute the initial error
 */
 error0 = InitialError(dtime, FVector, UVector, VVector, UVector0);
 
 do {
        
  LINEAR( dtime, FVector, UVector, VVector, &iteration, &error, error0, loadstep, nna, nnu, ia, ja, avector );
   
  if( _count_it != NULL )
    (*_count_it) ( iteration, error, loadstep);
    
  if( _fstop != NULL )
    (*_fstop) ( &stop );  
    
  needite = 0;
  
  fprintf( stdout, "\tStep: %d             Error: %3.3e\r", iteration+1, error );
  fflush( stdout );
  
  if(Config.num_load_step == 1){
   /* Save results on a temporary file
   */
   DrvPrintResult( DrvObj, iteration, UVector, VVector );
  }
 } while( needite );

 if(Config.num_load_step > 1){
  /* Save results on a temporary file
  */
  DrvPrintResult( DrvObj, loadstep, UVector, VVector );
 }
 

 /* free error vector */
 free ( vetError );

 fprintf( stdout, "\n\tTotal solver time....................: %lu segs\n\n", (time(NULL)-i_time) );

 free( dUVector );
 free( UVector0 );
 free( ia );
 free( ja );
 free( avector );
 
 return( 1 );

} /* End of IMPLINEARSolver */

/*
%F
*/
int IMPBFGSSolver(UI_State *R, double *FVector, double *UVector, double *VVector,
               double *MVector, int loadstep)
{
 int    iteration, needite;
 double dtime, alpha, CinEng, error;
 double *dUVector, *UVector0;
 double error0;
 int    stop = 0;
 unsigned long int i_time = time( NULL );
 int i;
 int    *ia;
 int    *ja;
 double *avector;
 int    nna = 0;
 int    nnu = 0;
 CsrNode   *node_list = 0L;      /* lista de nós Csr*/
 
/* Initialize variables
 */
 iteration = 0;
 needite   = 0;
 dtime     = 0.0;
 alpha     = 0.0;
 error     = 0.0;
 CinEng    = 0.0;
 error0    = 1.0;
 Config.numsteptol    = 1;
 
 dUVector  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 UVector0  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 
 nnu = NDof*NumNodes;
 
 /* aloca memória para lista de nós
 */
 node_list = (CsrNode*) calloc(NumNodes + 1, sizeof(CsrNode));
 
 /* constroi a lista de nós
 */
 BuildNodeList(node_list);

 /* constroi a lista de nós
 */
 BuildNodeList(node_list);
 
 /* aloca memória para ia*/
 ia  = (int*)calloc( nnu + 1, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildIaVec( node_list, ia, nnu, &nna );
 
 /* aloca memória para ja*/
 ja  = (int*)calloc( nna, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildJaVec( node_list, ja );
 
 /* libera memória da lista de nós
 */
 FreeNodeList(node_list);
 
 /* aloca memória para ia*/
 avector  = (double*)calloc( nna, sizeof(double) );
 
 /* cria o vetor avector*/
 BuildAVector( ia, ja, avector );
 
 if(Config.solver == 1){
 
 /* Mudança de índices necessárias
    para utilização do solver SAMG*/
  for( i = 0; i < nnu+1; i++ ) {
   ia[i] += 1;
  }
  for( i = 0; i < nna; i++ ) {
   ja[i] += 1;
  }
  
 }
 
 if(R)
  UISetDisplacement(R, NDof, UVector);

/* Compute the time increment according Cundall (1989)
 */
 dtime  = Config.dtfrac * _TimeIncrement();

/* Mount the mass vector
 */
 MassVector( MVector );
 
/* allocate error vector */ 
 vetError = (double*)malloc(Config.numsteptol*sizeof(double));

 /* Faz analise para cada passo de carga
  */

 iteration = 0;
 
 Config.loadfactor = (loadstep + 1.0) / Config.num_load_step;
 
 /* Compute the initial error
 */
 error0 = InitialError(dtime, FVector, UVector, VVector, UVector0);
 
 do {
  
  BFGS( dtime, FVector, UVector, VVector, &iteration, &error, error0, loadstep, nna, nnu, ia, ja, avector );
  
  if( _count_it != NULL )
    (*_count_it) ( iteration, error, loadstep);
    
  if( _fstop != NULL )
    (*_fstop) ( &stop );  
    
  /* Check for another iteration
  */
  needite = NeedIteration( &error, &iteration, stop );
  
  if(Config.num_load_step == 1){
   /* Save results on a temporary file
   */
   DrvPrintResult( DrvObj, iteration, UVector, VVector );
  }
 } while( needite );

 if(Config.num_load_step > 1){
  /* Save results on a temporary file
  */
  DrvPrintResult( DrvObj, loadstep, UVector, VVector );
 }

 /* free error vector */
 free ( vetError );

 fprintf( stdout, "\n\tTotal solver time....................: %lu segs\n\n", (time(NULL)-i_time) );

 printf( "\n" );

 free( dUVector );
 free( UVector0 );
 free( ia );
 free( ja );
 free( avector );
 
 return( 1 );

} /* End of BFGSSolver */

/*
%F
*/
int IMPNRMSolver(UI_State *R, double *FVector, double *UVector, double *VVector,
              double *MVector, int loadstep)
{
 int    iteration, needite;
 double dtime, alpha, CinEng, error;
 double *dUVector, *UVector0;
 double error0;
 int    stop = 0;
 unsigned long int i_time = time( NULL );
 int i;
 int    *ia;
 int    *ja;
 double *avector;
 int    nna = 0;
 int    nnu = 0;
 CsrNode   *node_list = 0L;      /* lista de nós Csr*/
 
/* Initialize variables
 */
 iteration = 0;
 needite   = 0;
 dtime     = 0.0;
 alpha     = 0.0;
 error     = 0.0;
 CinEng    = 0.0;
 error0    = 1.0;
 Config.numsteptol    = 1;
 
 dUVector  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 UVector0  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 
 nnu = NDof*NumNodes;
 
 /* aloca memória para lista de nós
 */
 node_list = (CsrNode*) calloc(NumNodes + 1, sizeof(CsrNode));
 
 /* constroi a lista de nós
 */
 BuildNodeList(node_list);

 /* constroi a lista de nós
 */
 BuildNodeList(node_list);
 
 /* aloca memória para ia*/
 ia  = (int*)calloc( nnu + 1, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildIaVec( node_list, ia, nnu, &nna );
 
 /* aloca memória para ja*/
 ja  = (int*)calloc( nna, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildJaVec( node_list, ja );
 
 /* libera memória da lista de nós
 */
 FreeNodeList(node_list);
 
 /* aloca memória para ia*/
 avector  = (double*)calloc( nna, sizeof(double) );
 
 /* cria o vetor avector*/
 BuildAVector( ia, ja, avector );
 
 if(Config.solver == 1){
 
 /* Mudança de índices necessárias
    para utilização do solver SAMG*/
  for( i = 0; i < nnu+1; i++ ) {
   ia[i] += 1;
  }
  for( i = 0; i < nna; i++ ) {
   ja[i] += 1;
  }
  
 }

 if(R)
  UISetDisplacement(R, NDof, UVector);

/* Compute the time increment according Cundall (1989)
 */
 dtime  = Config.dtfrac * _TimeIncrement();

/* Mount the mass vector
 */
 MassVector( MVector );
 
/* allocate error vector */ 
 vetError = (double*)malloc(Config.numsteptol*sizeof(double));

 /* Faz analise para cada passo de carga
  */

 iteration = 0;
 
 Config.loadfactor = (loadstep + 1.0) / Config.num_load_step;
 
 /* Compute the initial error
 */
 error0 = InitialError(dtime, FVector, UVector, VVector, UVector0);
 
 do {
  
  NRM( dtime, FVector, UVector, VVector, &iteration, &error, error0, loadstep, nna, nnu, ia, ja, avector );
  
  if( _count_it != NULL )
    (*_count_it) ( iteration, error, loadstep);
    
  if( _fstop != NULL )
    (*_fstop) ( &stop );  
    
  /* Check for another iteration
  */
  needite = NeedIteration( &error, &iteration, stop );
  
  if(Config.num_load_step == 1){
   /* Save results on a temporary file
   */
   DrvPrintResult( DrvObj, iteration, UVector, VVector );
  }
 } while( needite );

 if(Config.num_load_step > 1){
  /* Save results on a temporary file
  */
  DrvPrintResult( DrvObj, loadstep, UVector, VVector );
 }

 /* free error vector */
 free ( vetError );

 fprintf( stdout, "\n\tTotal solver time....................: %lu segs\n\n", (time(NULL)-i_time) );

 printf( "\n" );

 free( dUVector );
 free( UVector0 );
 free( ia );
 free( ja );
 free( avector );
 
 return( 1 );

} /* End of NRMSolver */

/*
%F
*/
int HYBRIDSolver(UI_State *R, double *FVector, double *UVector, double *VVector,
               double *MVector, int loadstep)
{
 int    iteration, needite;
 double dtime, alpha, CinEng, error;
 double *dUVector, *UVector0;
 double error0;
 int    stop = 0;
 int    hs = 1;
 unsigned long int i_time = time( NULL );
 int i;
 int    *ia;
 int    *ja;
 double *avector;
 int    nna = 0;
 int    nnu = 0;
 CsrNode   *node_list = 0L;      /* lista de nós Csr*/
 
/* Initialize variables
 */
 iteration = 0;
 needite   = 0;
 dtime     = 0.0;
 alpha     = 0.0;
 error     = 0.0;
 CinEng    = 0.0;
 error0    = 1.0;
 Config.max_iter      = 20;
 Config.numsteptol    = 1;
 
 dUVector  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 UVector0  = (double *)calloc( NDof*NumNodes, sizeof(double) );
 
 nnu = NDof*NumNodes;
 
 /* aloca memória para lista de nós
 */
 node_list = (CsrNode*) calloc(NumNodes + 1, sizeof(CsrNode));
 
 /* constroi a lista de nós
 */
 BuildNodeList(node_list);

 /* constroi a lista de nós
 */
 BuildNodeList(node_list);
 
 /* aloca memória para ia*/
 ia  = (int*)calloc( nnu + 1, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildIaVec( node_list, ia, nnu, &nna );
 
 /* aloca memória para ja*/
 ja  = (int*)calloc( nna, sizeof(int) );
 
 /* cria o vetor ia*/
 BuildJaVec( node_list, ja );
 
 /* libera memória da lista de nós
 */
 FreeNodeList(node_list);
 
 /* aloca memória para ia*/
 avector  = (double*)calloc( nna, sizeof(double) );
 
 /* cria o vetor avector*/
 BuildAVector( ia, ja, avector );
 
 if(Config.solver == 1){
 
 /* Mudança de índices necessárias
    para utilização do solver SAMG*/
  for( i = 0; i < nnu+1; i++ ) {
   ia[i] += 1;
  }
  for( i = 0; i < nna; i++ ) {
   ja[i] += 1;
  }
  
 }


 if(R)
  UISetDisplacement(R, NDof, UVector);

/* Compute the time increment according Cundall (1989)
 */
 dtime  = Config.dtfrac * _TimeIncrement();

/* Mount the mass vector
 */
 MassVector( MVector );
 
/* allocate error vector */ 
 vetError = (double*)malloc(Config.numsteptol*sizeof(double));

 /* Faz analise para cada passo de carga
  */

 iteration = 0;
 
 Config.loadfactor = (loadstep + 1.0) / Config.num_load_step;
 
 /* Compute the initial error
 */
 error0 = InitialError(dtime, FVector, UVector, VVector, UVector0);
 
 do {
  
  if(hs){
   BFGS( dtime, FVector, UVector, VVector, &iteration, &error, error0, loadstep, nna, nnu, ia, ja, avector );
   if(error >= Config.tolerance ){
    /* nao convergiu com BFGS, passa para RD */
    hs = 0;      
    Config.max_iter   = 10000;
    Config.numsteptol = 50;
    iteration = 0;
    free ( vetError );
    /* allocate error vector */                                   
    vetError = (double*)malloc(Config.numsteptol*sizeof(double));
   }
  }
  if(!hs){   
  
   /* Compute the rigidity coefficient
   */
   DampCalcCoeff(DampObj,iteration,dtime,MVector,VVector,&alpha,&CinEng);
   
   /* Compute displacements and veclocities
   */
   if( !DisplacementVelocity( iteration, alpha, dtime, 
                             FVector, VVector, UVector, MVector ) ) 
    return 0;
   /*
   */
   if( R && (iteration%Config.draw_step) )
    UIDraw(R);
   
   /* Update geometry for the larger displacements case
   */
   if( IsGeomNonLinear )
    _UpdateGeometry( dtime, VVector );
   
   /* Compute the percolation forces
   */
   if( IsCoupled )
    PercolationForces( FVector );
   
   /* Compute the internal forces
   */
   InternalForces( dtime, FVector, UVector, VVector );
   
   /* Compute the dUVector
   */
   _dUVector(UVector, UVector0, dUVector);
    
   /* Compute the error
   */
   error = ComputeError( FVector , dUVector, error0);
   
   /* For larger displaments correct the time step and compute a new
   * mass vector for the model
   */
   if( IsGeomNonLinear ) {
    if( error > 1.0 ) {
     dtime = Config.dtfrac * _TimeIncrement();
     MassVector( MVector );
    }
   }
   
   if( _count_it != NULL )
     (*_count_it) ( iteration, error, loadstep);
     
   if( _fstop != NULL )
     (*_fstop) ( &stop );  
     
   /* Check for another iteration
   */
   needite = NeedIteration( &error, &iteration, stop );
   
   fprintf( stdout, "\tStep: %d             Error: %3.3e\r", iteration+1, error );
   fflush( stdout );
   
   if(Config.num_load_step == 1){
    /* Save results on a temporary file
    */
    DrvPrintResult( DrvObj, iteration, UVector, VVector );
   }
  }
 } while( needite );

 if(Config.num_load_step > 1){
  /* Save results on a temporary file
  */
  DrvPrintResult( DrvObj, loadstep, UVector, VVector );
 }

 /* free error vector */
 free ( vetError );

 fprintf( stdout, "\n\tTotal solver time....................: %lu segs\n\n", (time(NULL)-i_time) );

 printf( "\n" );

 free( dUVector );
 free( UVector0 );
 free( ia );
 free( ja );
 free( avector );
 
 return( 1 );

} /* End of HSSolver */


/*
%F This function reset this packtage.
*/
void AlgReset( void )
{
  Config.max_iter        = 0;
  Config.nonlinear       = 0;
  Config.print_step      = 0;
  Config.draw_step       = 1;
  Config.num_step        = 0;
  Config.num_time_step   = 0;
  Config.tolerance       = 0.0;
  Config.numsteptol      = 50;
  Config.dtfrac          = 0.0;
  Config.totalTime       = 0.0;
  Config.timeStep        = 0.0;  
  Config.num_load_step   = 1;
  Config.loadfactor      = 1;
  Config.algtype         = 0;
  Config.numgaussresults =-1;
  Config.gaussresults    = NULL;
  Config.solver          = 0;

} /* End of AlgReset */


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

---