root/src/dynboun.c

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

DEFINITIONS

This source file includes following definitions.
  1. dynboun

   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 <string.h>
  22 #include "CalculiX.h"
  23 
  24 #ifdef SPOOLES
  25    #include "spooles.h"
  26 #endif
  27 #ifdef SGI
  28    #include "sgi.h"
  29 #endif
  30 #ifdef TAUCS
  31    #include "tau.h"
  32 #endif
  33 #ifdef PARDISO
  34    #include "pardiso.h"
  35 #endif
  36 
  37 void dynboun(double *amta,ITG *namta,ITG *nam,double *ampli, double *time,
  38              double *ttime,double *dtime,double *xbounold,double *xboun,
  39              double *xbounact,ITG *iamboun,ITG *nboun,ITG *nodeboun,
  40              ITG *ndirboun, double *ad, double *au, double *adb,
  41              double *aub, ITG *icol, ITG *irow, ITG *neq, ITG *nzs,
  42              double *sigma, double *b, ITG *isolver,
  43              double *alpham, double *betam, ITG *nzl,
  44              ITG *init,double *bact, double *bmin, ITG *jq, 
  45              char *amname,double *bv, double *bprev, double *bdiff,
  46              ITG *nactmech, ITG *icorrect, ITG *iprev){
  47 
  48     ITG idiff[3],i,j,ic,ir,im,symmetryflag=0;
  49 
  50     double *xbounmin=NULL,*xbounplus=NULL,*bplus=NULL,
  51         *ba=NULL,deltatime,deltatime2,deltatimesq,timemin,ttimemin,
  52         timeplus,ttimeplus,*aux=NULL,*b1=NULL,*b2=NULL,*bnew=NULL;
  53 
  54 #ifdef SGI
  55   ITG token=1;
  56 #endif
  57     
  58   NNEW(xbounmin,double,*nboun);
  59   NNEW(xbounplus,double,*nboun);
  60 
  61       /* time increment for the calculation of the change of the
  62          particular solution (needed to account for nonzero
  63          SPC's) */
  64 
  65   deltatime=*dtime;
  66   deltatime2=2.*deltatime;
  67   deltatimesq=deltatime*deltatime;
  68 
  69       /* the SPC value at timemin is stored in xbounmin */
  70 
  71   if(*init==1){
  72 
  73       /* at the start of a new step it is assumed that the previous step
  74          has reached steady state (at least for the SPC conditions) */
  75 
  76       for(i=0;i<*nboun;i++){
  77           xbounmin[i]=xbounold[i];
  78           xbounact[i]=xbounold[i];
  79       }
  80   }
  81   else{
  82       timemin=*time-deltatime;
  83       ttimemin=*ttime-deltatime;
  84       FORTRAN(temploadmodal,(amta,namta,nam,ampli,&timemin,&ttimemin,dtime,
  85            xbounold,xboun,xbounmin,iamboun,nboun,nodeboun,ndirboun,
  86            amname));
  87   }
  88 
  89       /* the SPC value at timeplus is stored in xbounplus */
  90 
  91   timeplus=*time+deltatime;
  92   ttimeplus=*ttime+deltatime;
  93   FORTRAN(temploadmodal,(amta,namta,nam,ampli,&timeplus,&ttimeplus,dtime,
  94           xbounold,xboun,xbounplus,iamboun,nboun,nodeboun,ndirboun,
  95           amname));
  96 
  97   NNEW(bplus,double,neq[1]);
  98   NNEW(ba,double,neq[1]);
  99   NNEW(b1,double,neq[1]);
 100   NNEW(b2,double,neq[1]);
 101 
 102       /* check whether boundary conditions changed 
 103          comparision of min with prev */
 104 
 105   if(*init==1){
 106       for(i=0;i<*nboun;i++){
 107           ic=neq[1]+i;
 108           for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
 109               ir=irow[j]-1;
 110               bmin[ir]=bmin[ir]-au[j]*xbounmin[i];
 111           }
 112       }
 113       if(*isolver==0){
 114 #ifdef SPOOLES
 115           spooles_solve(bmin,&neq[1]);
 116 #endif
 117       }
 118       else if(*isolver==4){
 119 #ifdef SGI
 120           sgi_solve(bmin,token);
 121 #endif
 122       }
 123       else if(*isolver==5){
 124 #ifdef TAUCS
 125           tau_solve(bmin,&neq[1]);
 126 #endif
 127       }
 128       if(*isolver==7){
 129 #ifdef PARDISO
 130           pardiso_solve(bmin,&neq[1],&symmetryflag);
 131 #endif
 132       }
 133   }
 134 
 135   /* check whether boundary conditions changed 
 136      comparision of act with min */
 137   
 138   idiff[1]=0;
 139   for(i=0;i<*nboun;i++){
 140       if(fabs(xbounact[i]-xbounmin[i])>1.e-10){
 141           idiff[1]=1;
 142           break;
 143       }
 144   }
 145   if(*init==1){
 146       for(i=0;i<*nboun;i++){
 147           ic=neq[1]+i;
 148           for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
 149               ir=irow[j]-1;
 150               bact[ir]=bact[ir]-au[j]*xbounact[i];
 151           }
 152       }
 153       if(*isolver==0){
 154 #ifdef SPOOLES
 155           spooles_solve(bact,&neq[1]);
 156 #endif
 157       }
 158       else if(*isolver==4){
 159 #ifdef SGI
 160           sgi_solve(bact,token);
 161 #endif
 162       }
 163       else if(*isolver==5){
 164 #ifdef TAUCS
 165           tau_solve(bact,&neq[1]);
 166 #endif
 167       }
 168       if(*isolver==7){
 169 #ifdef PARDISO
 170           pardiso_solve(bact,&neq[1],&symmetryflag);
 171 #endif
 172       }
 173   }
 174 
 175       /* check whether boundary conditions changed 
 176          comparision of plus with act */
 177   
 178   idiff[2]=0;
 179   for(i=0;i<*nboun;i++){
 180       if(fabs(xbounplus[i]-xbounact[i])>1.e-10){
 181           idiff[2]=1;
 182           break;
 183       }
 184   }
 185   if(idiff[2]==1){
 186       for(i=0;i<*nboun;i++){
 187           ic=neq[1]+i;
 188           for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
 189               ir=irow[j]-1;
 190               bplus[ir]=bplus[ir]-au[j]*xbounplus[i];
 191           }
 192       }
 193       if(*isolver==0){
 194 #ifdef SPOOLES
 195           spooles_solve(bplus,&neq[1]);
 196 #endif
 197       }
 198       else if(*isolver==4){
 199 #ifdef SGI
 200           sgi_solve(bplus,token);
 201 #endif
 202       }
 203       else if(*isolver==5){
 204 #ifdef TAUCS
 205           tau_solve(bplus,&neq[1]);
 206 #endif
 207       }
 208       if(*isolver==7){
 209 #ifdef PARDISO
 210           pardiso_solve(bplus,&neq[1],&symmetryflag);
 211 #endif
 212       }
 213   }
 214   
 215   if((idiff[1]!=0)||(idiff[2]!=0)){
 216 
 217       /* present value is not zero */
 218 
 219       if(idiff[2]==0){
 220           for(i=0;i<neq[1];i++){bplus[i]=bact[i];}
 221       }
 222       for(i=0;i<neq[1];i++){
 223           
 224           /* bv is the velocity */
 225           
 226           bv[i]=(bplus[i]-bmin[i])/deltatime2;
 227           
 228           /* ba is the acceleration */
 229           
 230           ba[i]=(bmin[i]-2.*bact[i]+bplus[i])/deltatimesq;
 231           
 232           b1[i]=ba[i]+*alpham*bv[i];
 233           b2[i]=*betam*bv[i];
 234 
 235           bmin[i]=bact[i];
 236           bact[i]=bplus[i];
 237       }
 238       NNEW(bnew,double,neq[1]);
 239       FORTRAN(op,(&neq[1],b1,bplus,adb,aub,jq,irow));
 240       for(i=0;i<neq[1];i++){bnew[i]=-bplus[i];}
 241       FORTRAN(op,(&neq[1],b2,bplus,ad,au,jq,irow));
 242       if(*icorrect==2){
 243           for(i=0;i<neq[1];i++){
 244               bnew[i]-=bplus[i];
 245               b[i]+=bnew[i];
 246           }
 247       }else if(*icorrect==0){
 248           for(i=0;i<neq[1];i++){
 249               bnew[i]-=bplus[i];
 250               bdiff[i]=bnew[i]-bprev[i];
 251               b[i]+=bdiff[i];
 252 //            printf("dynboun %e,%e,%e,%e\n",bprev[i],bnew[i],bdiff[i],b[i]);
 253           }
 254           memcpy(&bprev[0],&bnew[0],sizeof(double)*neq[1]);
 255       }else{
 256           for(i=0;i<neq[1];i++){
 257               bnew[i]-=bplus[i];
 258               bdiff[i]+=bnew[i]-bprev[i];
 259               b[i]+=bdiff[i];
 260           }
 261           memcpy(&bprev[0],&bnew[0],sizeof(double)*neq[1]);
 262       }
 263       SFREE(bnew);
 264       *nactmech=neq[1];
 265       *iprev=1;
 266   }else if((*iprev!=0)&&(*icorrect!=2)){
 267 
 268       /* present value of b is zero, previous value was not zero */
 269 
 270       if(*icorrect==0){
 271           for(i=0;i<neq[1];i++){
 272               bdiff[i]=-bprev[i];
 273               b[i]+=bdiff[i];
 274 //            printf("dynboun %e,%e,%e,%e\n",bprev[i],bdiff[i],b[i]);
 275           }
 276 //        memset(&bprev[0],0.,sizeof(double)*neq[1]);
 277           DMEMSET(bprev,0,neq[1],0.);
 278       }else{
 279           for(i=0;i<neq[1];i++){
 280               bdiff[i]+=-bprev[i];
 281               b[i]+=bdiff[i];
 282           }
 283 //        memset(&bprev[0],0.,sizeof(double)*neq[1]);
 284           DMEMSET(bprev,0,neq[1],0.);
 285       }
 286       *nactmech=neq[1];
 287       *iprev=0;
 288   }
 289   
 290   SFREE(xbounmin);SFREE(xbounplus);
 291   SFREE(bplus);SFREE(ba);SFREE(b1);SFREE(b2);
 292   
 293   return;
 294 }

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