root/src/prediction.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. prediction

   1 /*     CalculiX - A 3-dimensional finite element program                 */
   2 /*              Copyright (C) 1998-2015 Guido Dhondt                          */
   3 
   4 /*     This program is free software; you can redistribute it and/or     */
   5 /*     modify it under the terms of the GNU General Public License as    */
   6 /*     published by the Free Software Foundation(version 2);    */
   7 /*                    */
   8 
   9 /*     This program is distributed in the hope that it will be useful,   */
  10 /*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
  11 /*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
  12 /*     GNU General Public License for more details.                      */
  13 
  14 /*     You should have received a copy of the GNU General Public License */
  15 /*     along with this program; if not, write to the Free Software       */
  16 /*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
  17 
  18 #include <stdio.h>
  19 #include <math.h>
  20 #include <stdlib.h>
  21 #include "CalculiX.h"
  22 #ifdef SPOOLES
  23    #include "spooles.h"
  24 #endif
  25 #ifdef SGI
  26    #include "sgi.h"
  27 #endif
  28 #ifdef TAUCS
  29    #include "tau.h"
  30 #endif
  31 
  32 
  33 void prediction(double *uam, ITG *nmethod, double *bet, double *gam, 
  34                double *dtime,
  35                ITG *ithermal, ITG *nk, double *veold, double *accold, double *v,
  36                ITG *iinc, ITG *idiscon, double *vold, ITG *nactdof, ITG *mi){
  37 
  38     ITG j,k,mt=mi[1]+1;
  39     double dextrapol,scal1,scal2;
  40 
  41     uam[0]=0.;
  42     uam[1]=0.;
  43     if(*nmethod==4){
  44       
  45         scal1=0.5*(1.-2.**bet)**dtime**dtime;
  46         scal2=(1.-*gam)**dtime;
  47         
  48         if(*ithermal<2){
  49             for(k=0;k<*nk;++k){
  50                 for(j=0;j<mt;j++){
  51                     dextrapol=*dtime*veold[mt*k+j]+scal1*accold[mt*k+j];
  52                     if((fabs(dextrapol)>uam[0])&&(nactdof[mt*k+j]>0)) {uam[0]=fabs(dextrapol);}
  53                     v[mt*k+j]=vold[mt*k+j]+dextrapol;
  54                     veold[mt*k+j]=veold[mt*k+j]+scal2*accold[mt*k+j];
  55                     accold[mt*k+j]=0.;
  56                 }
  57             }
  58         }else if(*ithermal==2){
  59             for(k=0;k<*nk;++k){
  60                 for(j=0;j<mt;j++){
  61                     v[mt*k+j]=vold[mt*k+j];
  62                 }
  63             }
  64             for(k=0;k<*nk;++k){
  65                 dextrapol=*dtime*veold[mt*k];
  66                 if(fabs(dextrapol)>100.) dextrapol=100.*dextrapol/fabs(dextrapol);
  67                 if((fabs(dextrapol)>uam[1])&&(nactdof[mt*k]>0)) {uam[1]=fabs(dextrapol);}
  68                 v[mt*k]+=dextrapol;
  69             }
  70         }else{
  71             for(k=0;k<*nk;++k){
  72                 for(j=0;j<mt;++j){
  73                     dextrapol=*dtime*veold[mt*k+j]+scal1*accold[mt*k+j];
  74                     if((j==0)&&fabs(dextrapol)>100.) dextrapol=100.*dextrapol/fabs(dextrapol);
  75                     if(j==0){
  76                         if((fabs(dextrapol)>uam[1])&&(nactdof[mt*k]>0)) {uam[1]=fabs(dextrapol);}
  77                     }else{
  78                         if((fabs(dextrapol)>uam[0])&&(nactdof[mt*k+j]>0)) {uam[0]=fabs(dextrapol);}
  79                     }
  80                     v[mt*k+j]=vold[mt*k+j]+dextrapol;
  81                     veold[mt*k+j]=veold[mt*k+j]+scal2*accold[mt*k+j];
  82                     accold[mt*k+j]=0.;
  83                 }
  84             }
  85         }
  86     }
  87     
  88     /* for the static case: extrapolation of the previous increment
  89        (if any within the same step) */
  90     
  91     else{
  92         if(*iinc>1){
  93             if(*ithermal<2){
  94                 for(k=0;k<*nk;++k){
  95                     for(j=0;j<mt;++j){
  96                         if(*idiscon==0){
  97                             dextrapol=*dtime*veold[mt*k+j];
  98                             if((fabs(dextrapol)>uam[0])&&(nactdof[mt*k+j]>0)) {uam[0]=fabs(dextrapol);} 
  99                             v[mt*k+j]=vold[mt*k+j]+dextrapol;
 100                         }else{
 101                             v[mt*k+j]=vold[mt*k+j];
 102                         }
 103                     }
 104                 }
 105             }else if(*ithermal==2){
 106                 for(k=0;k<*nk;++k){
 107                     for(j=0;j<mt;++j){
 108                         v[mt*k+j]=vold[mt*k+j];
 109                     }
 110                 }
 111                 for(k=0;k<*nk;++k){
 112                     if(*idiscon==0){
 113                         dextrapol=*dtime*veold[mt*k];
 114                         if(fabs(dextrapol)>100.) dextrapol=100.*dextrapol/fabs(dextrapol);
 115                         if((fabs(dextrapol)>uam[1])&&(nactdof[mt*k]>0)) {uam[1]=fabs(dextrapol);}       
 116                         v[mt*k]+=dextrapol;
 117                     }
 118                 }
 119             }else{
 120                 for(k=0;k<*nk;++k){
 121                     for(j=0;j<mt;++j){
 122                         if(*idiscon==0){
 123                             dextrapol=*dtime*veold[mt*k+j];
 124                             if((j==0)&&fabs(dextrapol)>100.) dextrapol=100.*dextrapol/fabs(dextrapol);
 125                             if(j==0){
 126                                 if((fabs(dextrapol)>uam[1])&&(nactdof[mt*k+j]>0)) {uam[1]=fabs(dextrapol);}
 127                             }else{
 128                                 if((fabs(dextrapol)>uam[0])&&(nactdof[mt*k+j]>0)) {uam[0]=fabs(dextrapol);}
 129                             }   
 130                             v[mt*k+j]=vold[mt*k+j]+dextrapol;
 131                         }else{
 132                             v[mt*k+j]=vold[mt*k+j];
 133                         }
 134                     }
 135                 }
 136             }
 137         }
 138         else{
 139             for(k=0;k<*nk;++k){
 140                 for(j=0;j<mt;++j){
 141                     v[mt*k+j]=vold[mt*k+j];
 142                 }
 143             }
 144         }
 145     }
 146     *idiscon=0;
 147 
 148   return;
 149 }

/* [<][>][^][v][top][bottom][index][help] */