xalm.cpp

Zobrazit dokumentaci tohoto souboru.

#include "xalm.h"
#include "gelib.h"
#include "arrays.h"

#include <stdio.h>
#include <stdlib.h>


namespace xalm {
  
// *****************************************************************************************
// *************************             VEŘEJNÉ METODY             ************************
// *****************************************************************************************
XALM :: XALM (void)
{
  step = 0;
  loadLevel = 0.0;
  status = 0;
  
  xalm_NR_Mode = xalm_NR_OldMode = xalm_modifiedNRM;
  xalm_NR_ModeTick = -1;  // do not swith to xalm_NR_OldMode
  xalm_MANRMSteps = 0;
  
  nsteps = 0;
  minStepLength = maxStepLength = initialStepLength = 0.0;
  nsMin =  0;
  nsReq =  5;
  nsMax = 60;
  
  Psi    = 0.0;       // displacement control on
  //Psi    = 1.0;     // load control on
  
  verbose = 1;
  rtolf = rtold = 1.e-3;
  
  neq = -1;
  
  // Variables for Hyper Plane Control
  xalm_Control = xalm_hpc_off; // HPControl is not default
}
  
XALM :: ~XALM()
{
  
}


// *****************************************************************************************
// *************************       CHRÁNĚNÉ METODY A ATRIBUTY       ************************
// *****************************************************************************************
void XALM :: initialize_and_check_consistency (void)
{
  if (neq == -1)  _errorr("Number of equations was not set.");
  
  internalForces.resize_ignore_vals(neq);
  internalForces.zero();
  
  totalDisplacement.resize_ignore_vals(neq);
  totalDisplacement.zero();
  
  incrementOfDisplacement.resize_ignore_vals(neq);
  incrementOfDisplacement.zero();
  
  //internalForcesEBENorm.resize_ignore_vals(neq);
  //internalForcesEBENorm.zero();
  
  xalm_NR_OldMode = xalm_NR_Mode;
}



int XALM :: solve (void)
{
  // zde se inicializuji pocatecni nastaveni tridy XALM, alokuji a resizuji pole
  // tato fce musi byt v kazde XALM_interface tride
  this->initialize();
  
  // Nastavení soukromých atributů a kontrola konzistence zadání.
  this->initialize_and_check_consistency();
  
  status = 1;
  
  // Iterační smyčka.
  for (step = 0; step <nsteps; step++) {
    this->solve_step (&incrementalLoadVector, &initialLoadVector,
                   &totalDisplacement, &incrementOfDisplacement,
                   &internalForces, step);
    
    if (status == 0) break;
    
    this->update_step();
  }
  
  return status;
}

void XALM :: solve_step (Dvctr *R, Dvctr *R0, Dvctr *X, Dvctr *dX, Dvctr *F,  long step)
{
  Dvctr rhs, deltaXt, deltaX_, dXm1,  ddX;
  Dvctr XInitial;

  double XX, RR, RR0, XR, p;
  double deltaLambda, Lambda, eta, DeltaLambdam1, DeltaLambda;
  double deltaL;
  double drProduct;
  int irest = 0;
  int HPsize, i, ind;
  double _RR, _XX;
  bool converged, errorOutOfRangeFlag;
  
  //
  deltaL = this->initialStepLength;
  
  if (verbose)  fprintf(stdout, "XALM:       Initial step length: %-15e\n", deltaL);
  if (verbose)  fprintf(stdout, "XALM:       Iteration       LoadLevel       ForceError      DisplError    \n");
  if (verbose)  fprintf(stdout, "----------------------------------------------------------------------------\n");
  
  
  //
  if ( xalm_NR_ModeTick == 0 )
    xalm_NR_Mode = xalm_NR_OldMode;
  else
    if ( xalm_NR_ModeTick > 0 )
      xalm_NR_ModeTick--;
  
  // zaloha pro pripad restartu
  XInitial.beCopyOf(X);

  ddX    .resize_ignore_vals(neq);  ddX    .zero();
  deltaXt.resize_ignore_vals(neq);  deltaXt.zero();
  

  
  if ( R0 )  RR0 = R0->compute_squared_norm();
  else       RR0 = 0.0;
  
  //
  // A  initial step (predictor)
  //
  //
  
  //
  // A.2.   We assume positive-definitive (0)Kt (tangent stiffness mtrx).
  //
  RR = R->compute_squared_norm();
  
 restart:
  //
  // A.1. calculation of (0)VarRt
  //
  dX->zero();
  this->update_stiffness_matrix(&totalDisplacement);
  this->lineq_solve(&deltaXt, R);
  
  p = 0.0;
  if ( xalm_Control == xalm_hpc_off ) {
    XX = deltaXt.compute_squared_norm();
    p = sqrt(XX + Psi * Psi * RR);
  }
  else if ( xalm_Control == xalm_hpc_on ) {
    HPsize = xalm_HPCIndirectDofMask.give_size();
    _XX = 0;
    _RR = 0;
    for (i=0; i < HPsize; i++) {
      if ( ( ind = xalm_HPCIndirectDofMask[i] - 1 ) != 0 ) {
        _XX += deltaXt[ind] * deltaXt[ind];
        _RR += (*R)[ind] * (*R)[ind];
      }
    }
    
    p = sqrt(_XX + Psi * Psi * _RR);
  }
  else if ( xalm_Control == xalml_hpc ) {
    HPsize = xalm_HPCIndirectDofMask.give_size();
    p = 0.;
    for (i=0; i < HPsize; i++) {
      if ( ( ind = xalm_HPCIndirectDofMask[i] -1 ) != 0 ) {
        p += deltaXt[ind] * xalm_HPCWeights[i];
      }
    }
    
  }
  
  XR = deltaXt.give_dotproduct(*R);
  
  /* XR is unscaled Bergan's param of current stiffness XR = deltaXt^T k deltaXt
   * this is used to test whether k has negative or positive slope */
  
  Lambda = loadLevel;
  DeltaLambda = deltaLambda = sgn(XR) * deltaL / p;
  Lambda += DeltaLambda;
  //
  // A.3.
  //
  
  rhs.beCopyOf(R);
  rhs.tmsby(DeltaLambda);
  if ( R0 ) {
    rhs.add(*R0);
  }
  this->lineq_solve(dX, &rhs);
  
  X->add(*dX);
  
  this->update_internal_forces(&internalForces, &totalDisplacement);
  
  int nite = 0;
  
  do {
    nite++;
    dXm1.beCopyOf(dX);
    DeltaLambdam1 = DeltaLambda;
    
    // Sestavení nové tečné matice tuhosti a deltaXt.
    if ( ( xalm_NR_Mode == xalm_fullNRM ) || ( xalm_NR_Mode == xalm_accelNRM  &&  nite % xalm_MANRMSteps == 0 ) ) {
      this->update_stiffness_matrix(&totalDisplacement);
      this->lineq_solve(&deltaXt, R);
    }
    
    rhs.beCopyOf(R);
    rhs.tmsby(Lambda);
    if (R0)  rhs.add(*R0);
    
    rhs.sbt(*F);
    deltaX_.resize_ignore_vals(neq);
    
    this->lineq_solve(&deltaX_, &rhs);
    eta = 1.0;
    //
    if ( this->computeDeltaLambda(deltaLambda, * dX, deltaXt, deltaX_, * R, RR, eta, deltaL, DeltaLambda, neq) ) {
      irest++;
      
      // Problem konvergence, snížení délky kroku a restart.
      if ( irest <= XALM_MAX_RESTARTS ) {
    status = 2;
    deltaL =  deltaL * XALM_RESET_STEP_REDUCE;
    if ( deltaL < minStepLength )
      deltaL = minStepLength;
    
        // Znovunastavení počátečních hodnot.
    X->beCopyOf(XInitial);
    dX->zero();
        this->update_stiffness_matrix(&totalDisplacement);
    
    if (verbose)  fprintf(stdout, "XALM:       Iteration Reset ...\n");
    
    xalm_NR_OldMode  = xalm_NR_Mode;
    xalm_NR_Mode     = xalm_fullNRM;
    xalm_NR_ModeTick = XALM_DEFAULT_NRM_TICKS;
    goto restart;
      }
      else {
    if (verbose) _warningg("XALM :: solve - can't continue further");
    status = 0;
        return;
      }
    }
    
    //
    ddX.resize_ignore_vals(neq);
    ddX.addtms(deltaXt, eta*deltaLambda);
    ddX.addtms(deltaX_, eta);
    dX->beCopyOf(dXm1);
    dX->add(ddX);
    X->beCopyOf(XInitial);
    X->add(*dX);
    
    drProduct = ddX.compute_squared_norm();
    
    //
    DeltaLambda = DeltaLambdam1 + deltaLambda;
    Lambda = loadLevel + DeltaLambda;
    
    this->update_internal_forces(&internalForces, &totalDisplacement);
    
    // Kontrola konvergence.
    converged = this->checkConvergence(* R, R0, * F, * X, ddX, Lambda, RR0, RR, drProduct,
                            nite, errorOutOfRangeFlag);
    
    if ( ( nite >= nsMax ) || errorOutOfRangeFlag ) {
      irest++;
      
      // Problem konvergence, snížení délky kroku a restart.
      if ( irest <= XALM_MAX_RESTARTS ) {
    status = 2;
    deltaL =  deltaL * XALM_RESET_STEP_REDUCE;
    if ( deltaL < minStepLength )
      deltaL = minStepLength;
    
        // Znovunastavení počátečních hodnot.
    X->beCopyOf(XInitial);
    dX->zero();
        this->update_stiffness_matrix(&totalDisplacement);
    
    if (verbose)  fprintf(stdout, "XALM:       Iteration Reset ...\n");
    
    xalm_NR_OldMode  = xalm_NR_Mode;
    xalm_NR_Mode     = xalm_fullNRM;
    xalm_NR_ModeTick = XALM_DEFAULT_NRM_TICKS;
    goto restart;
      }
      else {
    if (verbose) _warningg("XALM :: solve - can't continue further");
    status = 0;
        return;
      }
    }
  } while ( !converged || ( nite < nsMin ) );
  
  // Upravení délky kroku, aby bylo dosaženo požadovaného počtu iterací.
  if ( nite > nsReq ) {
    deltaL =  deltaL * nsReq / nite;
  }
  else {
    deltaL =  deltaL * sqrt( sqrt( ( double ) nsReq / ( double ) nite ) );
  }
  
  if (deltaL > maxStepLength)  deltaL = maxStepLength;
  if (deltaL < minStepLength)  deltaL = minStepLength;
  
  if (verbose)  fprintf(stdout, "XALM:       Adjusted step length: %-15e\n", deltaL);
  
  loadLevel = Lambda;
} // end of XALM :: solve

bool XALM :: checkConvergence (const Dvctr &R, const Dvctr *R0, const Dvctr &F,
                   const Dvctr &X, const Dvctr &ddX,
                   double Lambda, double RR0, double RR, double drProduct,
                    int nite, bool &errorOutOfRange)
{
  double forceErr, dispErr;
  Dvctr rhs; // residual of momentum balance eq (unbalanced nodal forces)
  bool answer;
  
  answer = true;
  errorOutOfRange = false;
  
  // compute residual vector
  rhs.beCopyOf(R);
  
  rhs.tmsby(Lambda);
  if ( R0 ) {
    rhs.add(*R0);
  }
  
  rhs.sbt(F);
  
  //
  // compute force error(s)
  //
  double dXX;
  
  // err is relative error of unbalanced forces
  forceErr = rhs.compute_squared_norm();
  // err is relative displacement change
  dXX = X.compute_squared_norm();
  
  // double eNorm = internalForcesEBENorm.give_sum();
  // we compute a relative error norm
  if ( ( RR0 + RR * Lambda * Lambda ) > XALM_SMALL_ERROR_NUM ) {
    forceErr = sqrt( forceErr / ( RR0 + RR * Lambda * Lambda ) );
  }
  // else if ( eNorm > XALM_SMALL_ERROR_NUM ) {
  //   forceErr = sqrt(forceErr / eNorm );
  // }
  else {
    forceErr = sqrt(forceErr);
  }
  
  //
  // compute displacement error
  //
  if ( dXX < XALM_SMALL_ERROR_NUM ) {
    dispErr = drProduct;
  } else {
    dispErr = drProduct / dXX;
    dispErr = sqrt(dispErr);
  }
  
  if ( ( fabs(forceErr) > rtolf * XALM_MAX_REL_ERROR_BOUND ) ||
       ( fabs(dispErr)  > rtold * XALM_MAX_REL_ERROR_BOUND ) ) {
    errorOutOfRange = true;
  }
    
  if ( ( fabs(forceErr) > rtolf ) || ( fabs(dispErr) > rtold ) ) {
    answer = false;
  }
  
  if (verbose)  fprintf (stdout, "XALM:       %-15d %-15e %-15e %-15e\n", nite, Lambda, forceErr, dispErr);
  
  
  return answer;
}



int XALM :: computeDeltaLambda (double &deltaLambda, const Dvctr &dX, const Dvctr &deltaXt,
                const Dvctr &deltaX_, const Dvctr &R, double RR, double eta,
                double deltaL, double DeltaLambda0, int neq)
{
  double a1 = 0.0, a2 = 0.0, a3 = 0.0, a4, a5;
  double _RR, _rr, _a2, _a3, _pr;
  double lam1, lam2, cos1, cos2;
  int i, ind, HPsize = 0;
  Dvctr help;
  double XX;
  
  //
  // B.3.
  //
  if ( ( xalm_Control == xalm_hpc_off ) || ( xalm_Control == xalm_hpc_on ) ) {
    if ( xalm_Control == xalm_hpc_off ) {
      // this two lines are necessary if NRM is used
      // (for MNRM they can be computed at startup A1).
      
      XX = deltaXt.compute_squared_norm();
      a1 = eta * eta * XX + Psi * Psi * RR;
      a2 = RR * Psi * Psi * DeltaLambda0 * 2.0;
      a2 += 2.0 * eta * dX.give_dotproduct(deltaXt);
      a2 += 2.0 * eta * eta * deltaX_.give_dotproduct(deltaXt);
      a3 = dX.compute_squared_norm();
      a3 += 2.0 *eta * dX.give_dotproduct(deltaX_);
      a3 += eta * eta * deltaX_.compute_squared_norm();
      a3 += DeltaLambda0 * DeltaLambda0 * RR * Psi * Psi - deltaL * deltaL;
      
    }
    else if ( xalm_Control == xalm_hpc_on ) {
      HPsize = xalm_HPCIndirectDofMask.give_size();
      _rr = 0.;
      _RR = 0.;
      _a2 = _a3 = 0.;
      for ( i = 0; i < HPsize; i++ ) {
    if ( ( ind = xalm_HPCIndirectDofMask[i] - 1 ) != 0 ) {
      _rr += deltaXt[ind] * deltaXt[ind];
      _RR += R[ind] * R[ind];
      _pr   = ( dX[ind] + eta * deltaX_[ind] );
      _a2 += eta * deltaXt[ind]  * _pr;
      _a3 += _pr * _pr;
    }
      }
      
      a1 = eta * eta * _rr + Psi * Psi * _RR;
      a2 = _RR * Psi * Psi * DeltaLambda0 * 2.0;
      a2 += 2.0 * _a2;
      a3 = _a3 - deltaL * deltaL + DeltaLambda0 * DeltaLambda0 * _RR * Psi * Psi;
      
      //printf ("\na1=%e, a2=%e, a3=%e", a1,a2,a3);
    }
    
    // solution of quadratic eqn.
    double discr = a2 * a2 - 4.0 * a1 * a3;
    
    if ( discr < 0.0 ) {
      if (verbose) _warningg("XALM :: computeDeltaLambda: discriminant is negative, solution failed");
      status = 0;
    }
    
    discr = sqrt(discr);
    lam1 = ( -a2 + discr ) / 2. / a1;
    lam2 = ( -a2 - discr ) / 2. / a1;
    
    // select better lam (according to angle between deltar0 and deltar1(2).
    //
    // up to there rewritten
    //
    if ( xalm_Control == xalm_hpc_off ) {
      a4 = eta * dX.give_dotproduct(deltaX_) + dX.compute_squared_norm();
      a5 = eta * dX.give_dotproduct(deltaXt);
    }
    else {
      a4 = 0.;
      a5 = 0.;
      for ( i = 0; i < HPsize; i++ ) {
    if ( ( ind = xalm_HPCIndirectDofMask[i] - 1 ) != 0 ) {
      a4 += eta * dX[ind] * deltaX_[ind] + dX[ind] * dX[ind];
      a5 += eta * dX[ind] * deltaXt[ind];
    }
      }
    }
    
    cos1 = ( a4 + a5 * lam1 ) / deltaL / deltaL;
    cos2 = ( a4 + a5 * lam2 ) / deltaL / deltaL;
    
    if ( cos1 > cos2 )  deltaLambda = lam1;
    else                deltaLambda = lam2;
    
    //printf ("eta=%e, lam1=%e, lam2=%e", eta, lam1, lam2);
  }
  else if ( xalm_Control == xalml_hpc ) {
    // linearized control
    double nom = 0., denom = 0.;
    
    HPsize = xalm_HPCIndirectDofMask.give_size();
    for ( i = 0; i < HPsize; i++ ) {
      if ( ( ind = xalm_HPCIndirectDofMask[i] - 1 ) != 0 ) {
    nom += ( dX[ind] + eta * deltaX_[ind] ) * xalm_HPCWeights[i];
    denom += eta * deltaXt[ind] * xalm_HPCWeights[i];
      }
    }
    
    if ( fabs(denom) < XALM_SMALL_NUM ) {
      if (verbose) _warningg("XALM :: \nxalm: zero denominator in linearized control");
      status = 0;
    }
    
    deltaLambda = ( deltaL - nom ) / denom;
  }
  
  return 0;
}



} // namespace xalm