1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 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
prediction_em(double * uam,ITG * nmethod,double * bet,double * gam,double * dtime,ITG * ithermal,ITG * nk,double * veold,double * v,ITG * iinc,ITG * idiscon,double * vold,ITG * nactdof,ITG * mi)33 void prediction_em(double *uam, ITG *nmethod, double *bet, double *gam,
34 double *dtime,
35 ITG *ithermal, ITG *nk, double *veold, double *v,
36 ITG *iinc, ITG *idiscon, double *vold, ITG *nactdof, ITG *mi){
37
38 ITG j,k,mt=mi[1]+1,jstart;
39 double dextrapol;
40
41 uam[0]=0.;
42 uam[1]=0.;
43
44 if(*ithermal<2){
45 jstart=1;
46 }else{
47 jstart=0;
48 }
49
50 if(*nmethod==4){
51 for(k=0;k<*nk;++k){
52 for(j=jstart;j<mt;j++){
53 dextrapol=*dtime*veold[mt*k+j];
54 if(fabs(dextrapol)>100.) dextrapol=100.*dextrapol/fabs(dextrapol);
55 if(j==0){
56 if((fabs(dextrapol)>uam[1])&&(nactdof[mt*k]>0)) {uam[1]=fabs(dextrapol);}
57 }
58 v[mt*k+j]=vold[mt*k+j]+dextrapol;
59 }
60 }
61 }
62
63 /* for the static case: extrapolation of the previous increment
64 (if any within the same step) */
65
66 else{
67 if(*iinc>1){
68 for(k=0;k<*nk;++k){
69 for(j=jstart;j<mt;j++){
70 if(*idiscon==0){
71 dextrapol=*dtime*veold[mt*k+j];
72 if(fabs(dextrapol)>100.) dextrapol=100.*dextrapol/fabs(dextrapol);
73 if(j==0){
74 if((fabs(dextrapol)>uam[1])&&(nactdof[mt*k]>0)) {uam[1]=fabs(dextrapol);}
75 }
76 v[mt*k+j]=vold[mt*k+j]+dextrapol;
77 }else{
78 v[mt*k+j]=vold[mt*k+j];
79 }
80 }
81 }
82 }
83 else{
84 for(k=0;k<*nk;++k){
85 for(j=jstart;j<mt;j++){
86 v[mt*k+j]=vold[mt*k+j];
87 }
88 }
89 }
90 }
91 *idiscon=0;
92
93 return;
94 }
95