root/src/calcresidual.c

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

DEFINITIONS

This source file includes following definitions.
  1. calcresidual

   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 calcresidual(ITG *nmethod, ITG *neq, double *b, double *fext, double *f,
  34         ITG *iexpl, ITG *nactdof, double *aux2, double *vold,
  35         double *vini, double *dtime, double *accold, ITG *nk, double *adb,
  36         double *aub, ITG *jq, ITG *irow, ITG *nzl, double *alpha,
  37         double *fextini, double *fini, ITG *islavnode, ITG *nslavnode,
  38         ITG *mortar, ITG *ntie,double *f_cm,
  39         double* f_cs, ITG *mi,ITG *nzs,ITG *nasym,ITG *idamping,
  40         double *veold,double *adc,double *auc,double *cvini,double *cv){
  41 
  42     ITG j,k,mt=mi[1]+1;
  43     double scal1;
  44       
  45     /* residual for a static analysis */
  46       
  47     if(*nmethod!=4){
  48         for(k=0;k<neq[1];++k){
  49             b[k]=fext[k]-f[k];
  50         }
  51     }
  52       
  53     /* residual for implicit dynamics */
  54       
  55     else if(*iexpl<=1){
  56         for(k=0;k<*nk;++k){
  57             if(nactdof[mt*k]!=0){
  58                 aux2[nactdof[mt*k]-1]=(vold[mt*k]-vini[mt*k])/(*dtime);}
  59             for(j=1;j<mt;++j){
  60                 if(nactdof[mt*k+j]!=0){aux2[nactdof[mt*k+j]-1]=accold[mt*k+j];}
  61             }
  62         }
  63         if(*nasym==0){
  64             FORTRAN(op,(&neq[1],aux2,b,adb,aub,jq,irow)); 
  65         }else{
  66             FORTRAN(opas,(&neq[1],aux2,b,adb,aub,jq,irow,nzs)); 
  67         }
  68         scal1=1.+*alpha;
  69         for(k=0;k<neq[0];++k){
  70             b[k]=scal1*(fext[k]-f[k])-*alpha*(fextini[k]-fini[k])-b[k];
  71         } 
  72         for(k=neq[0];k<neq[1];++k){
  73             b[k]=fext[k]-f[k]-b[k];
  74         } 
  75 
  76         /* correction for damping */
  77 
  78         if(*idamping==1){
  79             for(k=0;k<*nk;++k){
  80                 if(nactdof[mt*k]!=0){aux2[nactdof[mt*k]-1]=0.;}
  81                 for(j=1;j<mt;++j){
  82                     if(nactdof[mt*k+j]!=0){
  83                         aux2[nactdof[mt*k+j]-1]=veold[mt*k+j];}
  84                 }
  85             }
  86             if(*nasym==0){
  87                 FORTRAN(op,(&neq[1],aux2,cv,adc,auc,jq,irow));
  88             }else{
  89                 FORTRAN(opas,(&neq[1],aux2,cv,adc,auc,jq,irow,nzs)); 
  90             }
  91             for(k=0;k<neq[0];++k){
  92                 b[k]-=scal1*cv[k]-*alpha*cvini[k];
  93             }
  94         }
  95     }
  96 
  97     /* residual for explicit dynamics */
  98     
  99     else{
 100         for(k=0;k<*nk;++k){
 101             if(nactdof[mt*k]!=0){
 102                 aux2[nactdof[mt*k]-1]=(vold[mt*k]-vini[mt*k])/(*dtime);}
 103             for(j=1;j<mt;++j){
 104                 if(nactdof[mt*k+j]!=0){aux2[nactdof[mt*k+j]-1]=accold[mt*k+j];}
 105             }
 106         }
 107         scal1=1.+*alpha;
 108         for(k=0;k<neq[0];++k){
 109             b[k]=scal1*(fext[k]-f[k])-*alpha*(fextini[k]-fini[k])
 110                 -adb[k]*aux2[k];
 111         } 
 112         for(k=neq[0];k<neq[1];++k){
 113             b[k]=fext[k]-f[k]-adb[k]*aux2[k];
 114         } 
 115     }
 116 
 117     return;
 118 }

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