CalculiX  2.8
A Free Software Three-Dimensional Structural Finite Element Program
 All Classes Files Functions Variables Macros
dynboun.c File Reference
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "CalculiX.h"
Include dependency graph for dynboun.c:

Go to the source code of this file.

Functions

void dynboun (double *amta, ITG *namta, ITG *nam, double *ampli, double *time, double *ttime, double *dtime, double *xbounold, double *xboun, double *xbounact, ITG *iamboun, ITG *nboun, ITG *nodeboun, ITG *ndirboun, double *ad, double *au, double *adb, double *aub, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, double *sigma, double *b, ITG *isolver, double *alpham, double *betam, ITG *nzl, ITG *init, double *bact, double *bmin, ITG *jq, char *amname, double *bv, double *bprev, double *bdiff, ITG *nactmech, ITG *icorrect, ITG *iprev)
 

Function Documentation

void dynboun ( double *  amta,
ITG namta,
ITG nam,
double *  ampli,
double *  time,
double *  ttime,
double *  dtime,
double *  xbounold,
double *  xboun,
double *  xbounact,
ITG iamboun,
ITG nboun,
ITG nodeboun,
ITG ndirboun,
double *  ad,
double *  au,
double *  adb,
double *  aub,
ITG icol,
ITG irow,
ITG neq,
ITG nzs,
double *  sigma,
double *  b,
ITG isolver,
double *  alpham,
double *  betam,
ITG nzl,
ITG init,
double *  bact,
double *  bmin,
ITG jq,
char *  amname,
double *  bv,
double *  bprev,
double *  bdiff,
ITG nactmech,
ITG icorrect,
ITG iprev 
)

Definition at line 37 of file dynboun.c.

46  {
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 }
void spooles_solve(double *b, ITG *neq)
subroutine init(nktet, inodfa, ipofa, netet_)
Definition: init.f:19
void FORTRAN(addimdnodecload,(ITG *nodeforc, ITG *i, ITG *imdnode, ITG *nmdnode, double *xforc, ITG *ikmpc, ITG *ilmpc, ITG *ipompc, ITG *nodempc, ITG *nmpc, ITG *imddof, ITG *nmddof, ITG *nactdof, ITG *mi, ITG *imdmpc, ITG *nmdmpc, ITG *imdboun, ITG *nmdboun, ITG *ikboun, ITG *nboun, ITG *ilboun, ITG *ithermal))
void sgi_solve(double *b, ITG token)
subroutine op(n, x, y, ad, au, jq, irow)
Definition: op.f:25
#define DMEMSET(a, b, c, d)
Definition: CalculiX.h:45
subroutine temploadmodal(amta, namta, nam, ampli, time, ttime, dtime, xbounold, xboun, xbounact, iamboun, nboun, nodeboun, ndirboun, amname)
Definition: temploadmodal.f:19
#define SFREE(a)
Definition: CalculiX.h:41
void tau_solve(double *b, ITG *neq)
#define ITG
Definition: CalculiX.h:51
void pardiso_solve(double *b, ITG *neq, ITG *symmetryflag)
#define NNEW(a, b, c)
Definition: CalculiX.h:39
Hosted by OpenAircraft.com, (Michigan UAV, LLC)