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

Go to the source code of this file.

Functions

void expand (double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc, ITG *nmpc, ITG *nodeforc, ITG *ndirforc, double *xforc, ITG *nforc, ITG *nelemload, char *sideload, double *xload, ITG *nload, ITG *nactdof, ITG *neq, ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun, double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, double *alcon, ITG *nalcon, double *alzero, ITG *ielmat, ITG *ielorien, ITG *norien, double *orab, ITG *ntmat_, double *t0, ITG *ithermal, double *prestr, ITG *iprestr, double *vold, ITG *iperturb, double *sti, ITG *nzs, double *adb, double *aub, char *filab, double *eme, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double *xstate, ITG *npmat_, char *matname, ITG *mi, ITG *ics, double *cs, ITG *mpcend, ITG *ncmat_, ITG *nstate_, ITG *mcs, ITG *nkon, double *ener, char *jobnamec, char *output, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, ITG *nprint, char *prlab, char *prset, ITG *nener, double *trab, ITG *inotr, ITG *ntrans, double *ttime, double *fmpc, ITG *nev, double **zp, ITG *iamboun, double *xbounold, ITG *nsectors, ITG *nm, ITG *icol, ITG *irow, ITG *nzl, ITG *nam, ITG *ipompcold, ITG *nodempcold, double *coefmpcold, char *labmpcold, ITG *nmpcold, double *xloadold, ITG *iamload, double *t1old, double *t1, ITG *iamt1, double *xstiff, ITG **icolep, ITG **jqep, ITG **irowep, ITG *isolver, ITG *nzse, double **adbep, double **aubep, ITG *iexpl, ITG *ibody, double *xbody, ITG *nbody, double *cocon, ITG *ncocon, char *tieset, ITG *ntie, ITG *imddof, ITG *nmddof, ITG *imdnode, ITG *nmdnode, ITG *imdboun, ITG *nmdboun, ITG *imdmpc, ITG *nmdmpc, ITG **izdofp, ITG *nzdof, ITG *nherm, double *xmr, double *xmi, char *typeboun, ITG *ielprop, double *prop)
 

Function Documentation

void expand ( double *  co,
ITG nk,
ITG kon,
ITG ipkon,
char *  lakon,
ITG ne,
ITG nodeboun,
ITG ndirboun,
double *  xboun,
ITG nboun,
ITG ipompc,
ITG nodempc,
double *  coefmpc,
char *  labmpc,
ITG nmpc,
ITG nodeforc,
ITG ndirforc,
double *  xforc,
ITG nforc,
ITG nelemload,
char *  sideload,
double *  xload,
ITG nload,
ITG nactdof,
ITG neq,
ITG nmethod,
ITG ikmpc,
ITG ilmpc,
ITG ikboun,
ITG ilboun,
double *  elcon,
ITG nelcon,
double *  rhcon,
ITG nrhcon,
double *  alcon,
ITG nalcon,
double *  alzero,
ITG ielmat,
ITG ielorien,
ITG norien,
double *  orab,
ITG ntmat_,
double *  t0,
ITG ithermal,
double *  prestr,
ITG iprestr,
double *  vold,
ITG iperturb,
double *  sti,
ITG nzs,
double *  adb,
double *  aub,
char *  filab,
double *  eme,
double *  plicon,
ITG nplicon,
double *  plkcon,
ITG nplkcon,
double *  xstate,
ITG npmat_,
char *  matname,
ITG mi,
ITG ics,
double *  cs,
ITG mpcend,
ITG ncmat_,
ITG nstate_,
ITG mcs,
ITG nkon,
double *  ener,
char *  jobnamec,
char *  output,
char *  set,
ITG nset,
ITG istartset,
ITG iendset,
ITG ialset,
ITG nprint,
char *  prlab,
char *  prset,
ITG nener,
double *  trab,
ITG inotr,
ITG ntrans,
double *  ttime,
double *  fmpc,
ITG nev,
double **  zp,
ITG iamboun,
double *  xbounold,
ITG nsectors,
ITG nm,
ITG icol,
ITG irow,
ITG nzl,
ITG nam,
ITG ipompcold,
ITG nodempcold,
double *  coefmpcold,
char *  labmpcold,
ITG nmpcold,
double *  xloadold,
ITG iamload,
double *  t1old,
double *  t1,
ITG iamt1,
double *  xstiff,
ITG **  icolep,
ITG **  jqep,
ITG **  irowep,
ITG isolver,
ITG nzse,
double **  adbep,
double **  aubep,
ITG iexpl,
ITG ibody,
double *  xbody,
ITG nbody,
double *  cocon,
ITG ncocon,
char *  tieset,
ITG ntie,
ITG imddof,
ITG nmddof,
ITG imdnode,
ITG nmdnode,
ITG imdboun,
ITG nmdboun,
ITG imdmpc,
ITG nmdmpc,
ITG **  izdofp,
ITG nzdof,
ITG nherm,
double *  xmr,
double *  xmi,
char *  typeboun,
ITG ielprop,
double *  prop 
)

Definition at line 33 of file expand.c.

65  {
66 
67  /* calls the Arnoldi Package (ARPACK) for cyclic symmetry calculations */
68 
69  char *filabt,*tchar1=NULL,*tchar2=NULL,*tchar3=NULL,lakonl[2]=" \0";
70 
71  ITG *inum=NULL,k,idir,lfin,j,iout=0,index,inode,id,i,idof,im,
72  ielas,icmd,kk,l,nkt,icntrl,imag=1,icomplex,kkv,kk6,iterm,
73  lprev,ilength,ij,i1,i2,iel,ielset,node,indexe,nope,ml1,nelem,
74  *inocs=NULL,*ielcs=NULL,jj,l1,l2,is,nlabel,*nshcon=NULL,
75  nodeleft,*noderight=NULL,numnodes,ileft,kflag=2,itr,locdir,
76  neqh,j1,nodenew,mt=mi[1]+1,istep=1,iinc=1,
77  tint=-1,tnstart=-1,tnend=-1,tint2=-1,
78  noderight_,*izdof=*izdofp,iload,iforc,*iznode=NULL,nznode,ll,ne0,
79  icfd=0,*inomat=NULL,mortar=0,*islavact=NULL,
80  *islavnode=NULL,*nslavnode=NULL,*islavsurf=NULL;
81 
82  long long lint;
83 
84  double *stn=NULL,*v=NULL,*temp_array=NULL,*vini=NULL,*csmass=NULL,
85  *een=NULL,cam[5],*f=NULL,*fn=NULL,qa[3],*epn=NULL,summass,
86  *stiini=NULL,*emn=NULL,*emeini=NULL,*clearini=NULL,
87  *xstateini=NULL,theta,pi,*coefmpcnew=NULL,t[3],ctl,stl,
88  *stx=NULL,*enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,
89  *qfx=NULL,*qfn=NULL,xreal,ximag,*vt=NULL,sum,
90  *coefright=NULL,coef,a[9],ratio,reltime,
91  *shcon=NULL,*springarea=NULL,*z=*zp, *zdof=NULL, *thicke=NULL,
92  atrab[9],acs[9],diff,fin[3],fout[3],*sumi=NULL,
93  *vti=NULL,*pslavsurf=NULL,*pmastsurf=NULL,*cdn=NULL;
94 
95  /* dummy arguments for the results call */
96 
97  double *veold=NULL,*accold=NULL,bet,gam,dtime,time;
98 
99  pi=4.*atan(1.);
100  neqh=neq[1]/2;
101 
102  noderight_=10;
103  NNEW(noderight,ITG,noderight_);
104  NNEW(coefright,double,noderight_);
105 
106  NNEW(v,double,2*mt**nk);
107  NNEW(vt,double,mt**nk**nsectors);
108 
109  NNEW(fn,double,2*mt**nk);
110  NNEW(stn,double,12**nk);
111  NNEW(inum,ITG,*nk);
112  NNEW(stx,double,6*mi[0]**ne);
113 
114  nlabel=46;
115  NNEW(filabt,char,87*nlabel);
116  for(i=1;i<87*nlabel;i++) filabt[i]=' ';
117  filabt[0]='U';
118 
119  NNEW(temp_array,double,neq[1]);
120  NNEW(coefmpcnew,double,*mpcend);
121 
122  nkt=*nsectors**nk;
123 
124  /* assigning nodes and elements to sectors */
125 
126  NNEW(inocs,ITG,*nk);
127  NNEW(ielcs,ITG,*ne);
128  ielset=cs[12];
129  if((*mcs!=1)||(ielset!=0)){
130  for(i=0;i<*nk;i++) inocs[i]=-1;
131  for(i=0;i<*ne;i++) ielcs[i]=-1;
132  }
133  NNEW(csmass,double,*mcs);
134  if(*mcs==1) csmass[0]=1.;
135 
136  for(i=0;i<*mcs;i++){
137  is=cs[17*i];
138  // if(is==1) continue;
139  ielset=cs[17*i+12];
140  if(ielset==0) continue;
141  for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
142  if(ialset[i1]>0){
143  iel=ialset[i1]-1;
144  if(ipkon[iel]<0) continue;
145  ielcs[iel]=i;
146  indexe=ipkon[iel];
147  if(*mcs==1){
148  if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
149  else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
150  else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
151  else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
152  else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
153  else if (strcmp1(&lakon[8*iel+3],"6")==0)nope=6;
154  else if (strcmp1(&lakon[8*iel],"ES")==0){
155  lakonl[0]=lakon[8*iel+7];
156  nope=atoi(lakonl)+1;}
157  else continue;
158  }else{
159  nelem=iel+1;
160  FORTRAN(calcmass,(ipkon,lakon,kon,co,mi,&nelem,ne,thicke,
161  ielmat,&nope,t0,t1,rhcon,nrhcon,ntmat_,
162  ithermal,&csmass[i]));
163  }
164  for(i2=0;i2<nope;++i2){
165  node=kon[indexe+i2]-1;
166  inocs[node]=i;
167  }
168  }
169  else{
170  iel=ialset[i1-2]-1;
171  do{
172  iel=iel-ialset[i1];
173  if(iel>=ialset[i1-1]-1) break;
174  if(ipkon[iel]<0) continue;
175  ielcs[iel]=i;
176  indexe=ipkon[iel];
177  if(*mcs==1){
178  if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
179  else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
180  else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
181  else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
182  else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
183  else {nope=6;}
184  }else{
185  nelem=iel+1;
186  FORTRAN(calcmass,(ipkon,lakon,kon,co,mi,&nelem,ne,thicke,
187  ielmat,&nope,t0,t1,rhcon,nrhcon,ntmat_,
188  ithermal,&csmass[i]));
189  }
190  for(i2=0;i2<nope;++i2){
191  node=kon[indexe+i2]-1;
192  inocs[node]=i;
193  }
194  }while(1);
195  }
196  }
197 // printf("expand.c mass = %" ITGFORMAT ",%e\n",i,csmass[i]);
198  }
199 
200  /* copying imdnode into iznode
201  iznode contains the nodes in which output is requested and
202  the nodes in which loading is applied */
203 
204  NNEW(iznode,ITG,*nk);
205  for(j=0;j<*nmdnode;j++){iznode[j]=imdnode[j];}
206  nznode=*nmdnode;
207 
208 /* expanding imddof, imdnode, imdboun and imdmpc */
209 
210  for(i=1;i<*nsectors;i++){
211  for(j=0;j<*nmddof;j++){
212  imddof[i**nmddof+j]=imddof[j]+i*neqh;
213  }
214  for(j=0;j<*nmdnode;j++){
215  imdnode[i**nmdnode+j]=imdnode[j]+i**nk;
216  }
217  for(j=0;j<*nmdboun;j++){
218  imdboun[i**nmdboun+j]=imdboun[j]+i**nboun;
219  }
220  for(j=0;j<*nmdmpc;j++){
221  imdmpc[i**nmdmpc+j]=imdmpc[j]+i**nmpc;
222  }
223  }
224  (*nmddof)*=(*nsectors);
225  (*nmdnode)*=(*nsectors);
226  (*nmdboun)*=(*nsectors);
227  (*nmdmpc)*=(*nsectors);
228 
229 /* creating a field with the degrees of freedom in which the eigenmodes
230  are needed:
231  1. all dofs in which the solution is needed (=imddof)
232  2. all dofs in which loading was applied
233  */
234 
235  NNEW(izdof,ITG,neqh**nsectors);
236  for(j=0;j<*nmddof;j++){izdof[j]=imddof[j];}
237  *nzdof=*nmddof;
238 
239  /* generating the coordinates for the other sectors */
240 
241  icntrl=1;
242 
243  FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filabt,&imag,mi,emn));
244 
245  for(jj=0;jj<*mcs;jj++){
246  is=(ITG)(cs[17*jj]+0.5);
247  for(i=1;i<is;i++){
248 
249  theta=i*2.*pi/cs[17*jj];
250 
251  for(l=0;l<*nk;l++){
252  if(inocs[l]==jj){
253  co[3*l+i*3**nk]=co[3*l];
254  co[1+3*l+i*3**nk]=co[1+3*l]+theta;
255  co[2+3*l+i*3**nk]=co[2+3*l];
256  if(*ntrans>0) inotr[2*l+i*2**nk]=inotr[2*l];
257  }
258  }
259  for(l=0;l<*nkon;l++){kon[l+i**nkon]=kon[l]+i**nk;}
260  for(l=0;l<*ne;l++){
261  if(ielcs[l]==jj){
262  if(ipkon[l]>=0){
263  ipkon[l+i**ne]=ipkon[l]+i**nkon;
264  ielmat[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
265  if(*norien>0) ielorien[l+i**ne]=ielorien[l];
266  for(l1=0;l1<8;l1++){
267  l2=8*l+l1;
268  lakon[l2+i*8**ne]=lakon[l2];
269  }
270  }else{
271  ipkon[l+i**ne]=ipkon[l];
272  }
273  }
274  }
275  }
276  }
277 
278  icntrl=-1;
279 
280  FORTRAN(rectcyl,(co,vt,fn,stn,qfn,een,cs,&nkt,&icntrl,t,filabt,&imag,mi,emn));
281 
282 /* expand nactdof */
283 
284  for(i=1;i<*nsectors;i++){
285  lint=i*mt**nk;
286  for(j=0;j<mt**nk;j++){
287  if(nactdof[j]!=0){
288  nactdof[lint+j]=nactdof[j]+i*neqh;
289  }else{
290  nactdof[lint+j]=0;
291  }
292  }
293  }
294 
295 /* copying the boundary conditions
296  (SPC's must be defined in cylindrical coordinates) */
297 
298  for(i=1;i<*nsectors;i++){
299  for(j=0;j<*nboun;j++){
300  nodeboun[i**nboun+j]=nodeboun[j]+i**nk;
301  ndirboun[i**nboun+j]=ndirboun[j];
302  xboun[i**nboun+j]=xboun[j];
303  xbounold[i**nboun+j]=xbounold[j];
304  if(*nam>0) iamboun[i**nboun+j]=iamboun[j];
305  ikboun[i**nboun+j]=ikboun[j]+8*i**nk;
306  ilboun[i**nboun+j]=ilboun[j]+i**nboun;
307  }
308  }
309 
310  /* distributed loads */
311 
312  for(i=0;i<*nload;i++){
313  if(nelemload[2*i+1]<*nsectors){
314  nelemload[2*i]+=*ne*nelemload[2*i+1];
315  }else{
316  nelemload[2*i]+=*ne*(nelemload[2*i+1]-(*nsectors));
317  }
318  iload=i+1;
319  FORTRAN(addizdofdload,(nelemload,sideload,ipkon,kon,lakon,
320  nactdof,izdof,nzdof,mi,&iload,iznode,&nznode,nk,
321  imdnode,nmdnode));
322  }
323 
324  /* body loads */
325 
326  if(*nbody>0){
327  printf("*ERROR in expand: body loads are not allowed for modal dynamics\n and steady state dynamics calculations in cyclic symmetric structures\n\n");
328  FORTRAN(stop,());
329  }
330 
331  /* sorting the elements with distributed loads */
332 
333  if(*nload>0){
334  if(*nam>0){
335  FORTRAN(isortiiddc,(nelemload,iamload,xload,xloadold,sideload,nload,&kflag));
336  }else{
337  FORTRAN(isortiddc,(nelemload,xload,xloadold,sideload,nload,&kflag));
338  }
339  }
340 
341  /* point loads */
342 
343 /* for(i=0;i<*nforc;i++){
344  if(nodeforc[2*i+1]<*nsectors){
345  nodeforc[2*i]+=*nk*nodeforc[2*i+1];
346  }else{
347  nodeforc[2*i]+=*nk*(nodeforc[2*i+1]-(*nsectors));
348  }
349  iforc=i+1;
350  FORTRAN(addizdofcload,(nodeforc,ndirforc,nactdof,mi,izdof,
351  nzdof,&iforc,iznode,&nznode,nk,imdnode,nmdnode,xforc));
352  }*/
353 
354  i=0;
355  while(i<*nforc){
356  node=nodeforc[2*i];
357 
358  /* checking for a cylindrical transformation;
359  comparison with the cyclic symmetry system */
360 
361  itr=inotr[2*node-2];
362 
363  if(itr==0){
364 
365  /* carthesian coordinate system */
366 
367  if(nodeforc[2*i+1]<*nsectors){
368  nodeforc[2*i]+=*nk*nodeforc[2*i+1];
369  }else{
370  nodeforc[2*i]+=*nk*(nodeforc[2*i+1]-(*nsectors));
371  }
372  i++;iforc=i;
373  FORTRAN(addizdofcload,(nodeforc,ndirforc,nactdof,mi,izdof,
374  nzdof,&iforc,iznode,&nznode,nk,imdnode,nmdnode,xforc));
375  }else{
376 
377  /* cylindrical coordinate system */
378 
379  FORTRAN(transformatrix,(&trab[7*(itr-1)],&co[3*(node-1)],atrab));
380  FORTRAN(transformatrix,(&cs[5],&co[3*(node-1)],acs));
381  diff=0.; for(j=0;j<9;j++) diff+=(atrab[j]-acs[j])*(atrab[j]-acs[j]);
382 
383  if((ndirforc[i]!=1)||
384  (nodeforc[2*i+2]!=node)||(ndirforc[i+1]!=2)||
385  (nodeforc[2*i+4]!=node)||(ndirforc[i+2]!=3)||
386  ((diff>1.e-10)&&(fabs(diff-8.)>1.e-10))){
387  printf("*ERROR: forces in a modal dynamic or steady state dynamics\n");
388  printf(" calculation with cyclic symmetry must be defined in\n");
389  printf(" the cyclic symmetric cylindrical coordinate system\n");
390  printf(" force at fault is applied in node %" ITGFORMAT "\n",node);
391  FORTRAN(stop,());
392  }
393 
394  /* changing the complete force in the node in the basis sector from
395  the global rectangular system into the cylindrical system */
396 
397  fin[0]=xforc[i];
398  fin[1]=xforc[i+1];
399  fin[2]=xforc[i+2];
400  icntrl=2;
401  FORTRAN(rectcyltrfm,(&node,co,cs,&icntrl,fin,fout));
402 
403  /* new node number (= node number in the target sector) */
404 
405  if(nodeforc[2*i+1]<*nsectors){
406  nodeforc[2*i]+=*nk*nodeforc[2*i+1];
407  }else{
408  nodeforc[2*i]+=*nk*(nodeforc[2*i+1]-(*nsectors));
409  }
410  nodeforc[2*i+2]=nodeforc[2*i];
411  nodeforc[2*i+4]=nodeforc[2*i];
412 
413  /* changing the complete force in the node in the target sector from
414  the cylindrical system into the global rectangular system */
415 
416  node=nodeforc[2*i];
417  fin[0]=fout[0];
418  fin[1]=fout[1];
419  fin[2]=fout[2];
420  icntrl=-2;
421  FORTRAN(rectcyltrfm,(&node,co,cs,&icntrl,fin,fout));
422  xforc[i]=fout[0];
423  xforc[i+1]=fout[1];
424  xforc[i+2]=fout[2];
425 
426  /* storing the node and the dof into iznode and izdof */
427 
428  for(j=0;j<3;j++){
429  i++;iforc=i;
430  FORTRAN(addizdofcload,(nodeforc,ndirforc,nactdof,mi,izdof,
431  nzdof,&iforc,iznode,&nznode,nk,imdnode,nmdnode,xforc));
432  }
433  }
434  }
435 
436  /* loop over all eigenvalues; the loop starts from the highest eigenvalue
437  so that the reuse of z is not a problem
438  z before: real and imaginary part for a segment for all eigenvalues
439  z after: real part for all segments for all eigenvalues */
440 
441  if(*nherm==1){
442  NNEW(zdof,double,(long long)*nev**nzdof);
443  }else{
444  NNEW(zdof,double,(long long)2**nev**nzdof);
445  NNEW(sumi,double,*nev);
446  }
447 
448  lfin=0;
449  for(j=*nev-1;j>-1;--j){
450  lint=2*j*neqh;
451 
452  /* calculating the cosine and sine of the phase angle */
453 
454  for(jj=0;jj<*mcs;jj++){
455  theta=nm[j]*2.*pi/cs[17*jj];
456  cs[17*jj+14]=cos(theta);
457  cs[17*jj+15]=sin(theta);
458  }
459 
460  /* generating the cyclic MPC's (needed for nodal diameters
461  different from 0 */
462 
463  NNEW(eei,double,6*mi[0]**ne);
464 
465  DMEMSET(v,0,2*mt**nk,0.);
466 
467  for(k=0;k<2*neqh;k+=neqh){
468 
469  for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
470 
471  if(k==0) {kk=0;kkv=0;kk6=0;}
472  else {kk=*nk;kkv=mt**nk;kk6=6**nk;}
473  for(i=0;i<*nmpc;i++){
474  index=ipompc[i]-1;
475  /* check whether thermal mpc */
476  if(nodempc[3*index+1]==0) continue;
477  coefmpcnew[index]=coefmpc[index];
478  while(1){
479  index=nodempc[3*index+2];
480  if(index==0) break;
481  index--;
482 
483  icomplex=0;
484  inode=nodempc[3*index];
485  if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
486  icomplex=atoi(&labmpc[20*i+6]);}
487  else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
488  for(ij=0;ij<*mcs;ij++){
489  lprev=cs[ij*17+13];
490  ilength=cs[ij*17+3];
491  FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
492  if(id!=0){
493  if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
494  }
495  }
496  }
497 
498  if(icomplex!=0){
499  idir=nodempc[3*index+1];
500  idof=nactdof[mt*(inode-1)+idir]-1;
501  if(idof==-1){xreal=1.;ximag=1.;}
502  else{xreal=z[lint+idof];ximag=z[lint+idof+neqh];}
503  if(k==0) {
504  if(fabs(xreal)<1.e-30)xreal=1.e-30;
505  coefmpcnew[index]=coefmpc[index]*
506  (cs[17*(icomplex-1)+14]+
507  ximag/xreal*cs[17*(icomplex-1)+15]);}
508  else {
509  if(fabs(ximag)<1.e-30)ximag=1.e-30;
510  coefmpcnew[index]=coefmpc[index]*
511  (cs[17*(icomplex-1)+14]-
512  xreal/ximag*cs[17*(icomplex-1)+15]);}
513  }
514  else{coefmpcnew[index]=coefmpc[index];}
515  }
516  }
517 
518  results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
519  stx,elcon,
520  nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
521  norien,orab,ntmat_,t0,t0,ithermal,
522  prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
523  f,&fn[kkv],nactdof,&iout,qa,vold,&z[lint+k],
524  nodeboun,ndirboun,xboun,nboun,ipompc,
525  nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neqh,veold,accold,
526  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
527  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
528  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
529  xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
530  ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
531  nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,&reltime,
532  &ne0,xforc,nforc,thicke,shcon,nshcon,
533  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
534  &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
535  islavsurf,ielprop,prop);
536 
537  }
538  SFREE(eei);
539 
540  /* mapping the results to the other sectors */
541 
542  if(*nherm!=1)NNEW(vti,double,mt**nk**nsectors);
543 
544  icntrl=2;imag=1;
545 
546  FORTRAN(rectcylexp,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filabt,&imag,mi,
547  iznode,&nznode,nsectors,nk,emn));
548 
549  /* basis sector */
550 
551  for(ll=0;ll<nznode;ll++){
552  l1=iznode[ll]-1;
553  for(l2=0;l2<mt;l2++){
554  l=mt*l1+l2;
555  vt[l]=v[l];
556  if(*nherm!=1)vti[l]=v[l+mt**nk];
557  }
558  }
559 
560  /* other sectors */
561 
562  for(jj=0;jj<*mcs;jj++){
563  ilength=cs[17*jj+3];
564  lprev=cs[17*jj+13];
565  for(i=1;i<*nsectors;i++){
566 
567  theta=i*nm[j]*2.*pi/cs[17*jj];
568  ctl=cos(theta);
569  stl=sin(theta);
570 
571  for(ll=0;ll<nznode;ll++){
572  l1=iznode[ll]-1;
573  if(inocs[l1]==jj){
574 
575  /* check whether node lies on axis */
576 
577  ml1=-l1-1;
578  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
579  if(id!=0){
580  if(ics[lprev+id-1]==ml1){
581  for(l2=0;l2<mt;l2++){
582  l=mt*l1+l2;
583  vt[l+mt**nk*i]=v[l];
584  if(*nherm!=1)vti[l+mt**nk*i]=v[l+mt**nk];
585  }
586  continue;
587  }
588  }
589  for(l2=0;l2<mt;l2++){
590  l=mt*l1+l2;
591  vt[l+mt**nk*i]=ctl*v[l]-stl*v[l+mt**nk];
592  if(*nherm!=1)vti[l+mt**nk*i]=stl*v[l]+ctl*v[l+mt**nk];
593  }
594  }
595  }
596  }
597  }
598 
599  icntrl=-2;imag=0;
600 
601  FORTRAN(rectcylexp,(co,vt,fn,stn,qfn,een,cs,&nkt,&icntrl,t,filabt,
602  &imag,mi,iznode,&nznode,nsectors,nk,emn));
603 
604 /* storing the displacements into the expanded eigenvectors */
605 
606  for(ll=0;ll<nznode;ll++){
607  i=iznode[ll]-1;
608 // for(i=0;i<*nk;i++){
609  for(j1=0;j1<mt;j1++){
610 
611  for(k=0;k<*nsectors;k++){
612  /* C-convention */
613  idof=nactdof[mt*(k**nk+i)+j1]-1;
614  if(idof!=-1){
615  FORTRAN(nident,(izdof,&idof,nzdof,&id));
616  if(id!=0){
617  if(izdof[id-1]==idof){
618  if(*nherm==1){
619  zdof[(long long)j**nzdof+id-1]=vt[k*mt**nk+mt*i+j1];
620  }else{
621  zdof[(long long)2*j**nzdof+id-1]=vt[k*mt**nk+mt*i+j1];
622  zdof[(long long)(2*j+1)**nzdof+id-1]=vti[k*mt**nk+mt*i+j1];
623  }
624  }
625  }
626  }
627  }
628  }
629  }
630 
631  if(*nherm!=1) SFREE(vti);
632 
633 /* normalizing the eigenvectors with the mass */
634 
635 /* if (nm[j]==0||(nm[j]==(ITG)((cs[0]/2))&&(fmod(cs[0],2.)==0.)))
636  {sum=sqrt(cs[0]);}
637  else{sum=sqrt(cs[0]/2);}*/
638 
639  sum=0.;
640  summass=0.;
641  for(i=0;i<*mcs;i++){
642  if (nm[j]==0||(nm[j]==(ITG)((cs[17*i]/2))&&(fmod(cs[17*i],2.)==0.))){
643  sum+=cs[17*i]*csmass[i];
644  }else{
645  sum+=cs[17*i]*csmass[i]/2.;
646  }
647  summass+=csmass[i];
648  }
649  if(fabs(summass)>1.e-20){
650  sum=sqrt(sum/summass);
651  }else{
652  printf("*ERROR in expand.c: total mass of structure is zero\n");
653  printf(" maybe no element sets were specified for the\n");
654  printf(" cyclic symmetry ties\n");
655  FORTRAN(stop,());
656  }
657 
658  if(*nherm==1){
659  for(i=0;i<*nzdof;i++){zdof[(long long)j**nzdof+i]/=sum;}
660  }else{
661  for(i=0;i<*nzdof;i++){zdof[(long long)(2*j)**nzdof+i]/=sum;}
662  for(i=0;i<*nzdof;i++){zdof[(long long)(2*j+1)**nzdof+i]/=sum;}
663  sumi[j]=sqrt(sum);
664  }
665  }
666 
667 /* copying zdof into z */
668 
669  if(*nherm==1){
670  RENEW(z,double,(long long)*nev**nzdof);
671  memcpy(&z[0],&zdof[0],(long long)sizeof(double)**nev**nzdof);
672  }else{
673  RENEW(z,double,(long long)2**nev**nzdof);
674  memcpy(&z[0],&zdof[0],(long long)sizeof(double)*2**nev**nzdof);
675  for(i=0;i<*nev;i++){
676  for(j=0;j<*nev;j++){
677  xmr[i**nev+j]/=(sumi[i]*sumi[j]);
678  xmi[i**nev+j]/=(sumi[i]*sumi[j]);
679  }
680  }
681  SFREE(sumi);
682  }
683  SFREE(zdof);
684 
685 /* copying the multiple point constraints */
686 
687  *nmpc=0;
688  *mpcend=0;
689  for(i=0;i<*nsectors;i++){
690  if(i==0){
691  ileft=*nsectors-1;
692  }else{
693  ileft=i-1;
694  }
695  for(j=0;j<*nmpcold;j++){
696  if(noderight_>10){
697  noderight_=10;
698  RENEW(noderight,ITG,noderight_);
699  RENEW(coefright,double,noderight_);
700  }
701  ipompc[*nmpc]=*mpcend+1;
702  ikmpc[*nmpc]=ikmpc[j]+8*i**nk;
703  ilmpc[*nmpc]=ilmpc[j]+i**nmpcold;
704  strcpy1(&labmpc[20**nmpc],&labmpcold[20*j],20);
705  if(strcmp1(&labmpcold[20*j],"CYCLIC")==0){
706  index=ipompcold[j]-1;
707  nodeleft=nodempcold[3*index];
708  idir=nodempcold[3*index+1];
709  index=nodempcold[3*index+2]-1;
710  numnodes=0;
711  do{
712  node=nodempcold[3*index];
713  if(nodempcold[3*index+1]==idir){
714  noderight[numnodes]=node;
715  coefright[numnodes]=coefmpcold[index];
716  numnodes++;
717  if(numnodes>=noderight_){
718  noderight_=(ITG)(1.5*noderight_);
719  RENEW(noderight,ITG,noderight_);
720  RENEW(coefright,double,noderight_);
721  }
722  }
723  index=nodempcold[3*index+2]-1;
724  if(index==-1) break;
725  }while(1);
726  if(numnodes>0){
727  sum=0.;
728  for(k=0;k<numnodes;k++){
729  sum+=coefright[k];
730  }
731  for(k=0;k<numnodes;k++){
732  coefright[k]/=sum;
733  }
734  }else{coefright[0]=1.;}
735  nodempc[3**mpcend]=nodeleft+i**nk;
736  nodempc[3**mpcend+1]=idir;
737  nodempc[3**mpcend+2]=*mpcend+2;
738  coefmpc[*mpcend]=1.;
739  for(k=0;k<numnodes;k++){
740  (*mpcend)++;
741  nodempc[3**mpcend]=noderight[k]+ileft**nk;
742  nodempc[3**mpcend+1]=idir;
743  nodempc[3**mpcend+2]=*mpcend+2;
744  coefmpc[*mpcend]=-coefright[k];
745  }
746  nodempc[3**mpcend+2]=0;
747  (*mpcend)++;
748  }else{
749  index=ipompcold[j]-1;
750  iterm=0;
751  do{
752  iterm++;
753  node=nodempcold[3*index];
754  idir=nodempcold[3*index+1];
755  coef=coefmpcold[index];
756 
757  /* check whether SUBCYCLIC MPC: if the current node
758  is an independent node of a CYCLIC MPC, the
759  node in the new MPC should be the cylic previous
760  one */
761 
762  nodenew=node+i**nk;
763  if(strcmp1(&labmpcold[20*j],"SUBCYCLIC")==0){
764  for(ij=0;ij<*mcs;ij++){
765  lprev=cs[ij*17+13];
766  ilength=cs[ij*17+3];
767  FORTRAN(nident,(&ics[lprev],&node,&ilength,&id));
768  if(id!=0){
769  if(ics[lprev+id-1]==node){
770  nodenew=node+ileft**nk;
771  break;
772  }
773  }
774  }
775  }
776 
777  /* modification of the MPC coefficients if
778  cylindrical coordinate system is active
779  it is assumed that all terms in the MPC are
780  either in the radial, the circumferential
781  or axial direction */
782 
783  if(*ntrans<=0){itr=0;}
784  else if(inotr[2*node-2]==0){itr=0;}
785  else{itr=inotr[2*node-2];}
786 
787  if(iterm==1) locdir=-1;
788 
789  if((itr!=0)&&(idir!=0)){
790  if(trab[7*itr-1]<0){
791  FORTRAN(transformatrix,(&trab[7*itr-7],
792  &co[3*node-3],a));
793  if(iterm==1){
794  for(k=0;k<3;k++){
795  if(fabs(a[3*k+idir-1]-coef)<1.e-10){
796  FORTRAN(transformatrix,(&trab[7*itr-7],
797  &co[3*nodenew-3],a));
798  coef=a[3*k+idir-1];
799  locdir=k;
800  break;
801  }
802  if(fabs(a[3*k+idir-1]+coef)<1.e-10){
803  FORTRAN(transformatrix,(&trab[7*itr-7],
804  &co[3*nodenew-3],a));
805  coef=-a[3*k+idir-1];
806  locdir=k;
807  break;
808  }
809  }
810  }else{
811  if(locdir!=-1){
812  if(fabs(a[3*locdir+idir-1])>1.e-10){
813  ratio=coef/a[3*locdir+idir-1];
814  }else{ratio=0.;}
815  FORTRAN(transformatrix,(&trab[7*itr-7],
816  &co[3*nodenew-3],a));
817  coef=ratio*a[3*locdir+idir-1];
818  }
819  }
820  }
821  }
822 
823  nodempc[3**mpcend]=nodenew;
824  nodempc[3**mpcend+1]=idir;
825  coefmpc[*mpcend]=coef;
826  index=nodempcold[3*index+2]-1;
827  if(index==-1) break;
828  nodempc[3**mpcend+2]=*mpcend+2;
829  (*mpcend)++;
830  }while(1);
831  nodempc[3**mpcend+2]=0;
832  (*mpcend)++;
833  }
834  (*nmpc)++;
835  }
836  }
837 
838  /* copying the temperatures */
839 
840  if(*ithermal!=0){
841  for(i=1;i<*nsectors;i++){
842  lint=i**nk;
843  for(j=0;j<*nk;j++){
844  t0[lint+j]=t0[j];
845  t1old[lint+j]=t1old[j];
846  t1[lint+j]=t1[j];
847  }
848  }
849  if(*nam>0){
850  for(i=1;i<*nsectors;i++){
851  lint=i**nk;
852  for(j=0;j<*nk;j++){
853  iamt1[lint+j]=iamt1[j];
854  }
855  }
856  }
857  }
858 
859  /* copying the contact definition */
860 
861  if(*nmethod==4){
862 
863  /* first find the startposition to append the expanded contact fields*/
864 
865  for(j=0; j<*nset; j++){
866  if(iendset[j]>tint){
867  tint=iendset[j];
868  }
869  }
870  tint++;
871  /* now append and expand the contact definitons*/
872  NNEW(tchar1,char,81);
873  NNEW(tchar2,char,81);
874  NNEW(tchar3,char,81);
875  for(i=0; i<*ntie; i++){
876  if(tieset[i*(81*3)+80]=='C'){
877  memcpy(tchar2,&tieset[i*(81*3)+81],81);
878  tchar2[80]='\0';
879  memcpy(tchar3,&tieset[i*(81*3)+81+81],81);
880  tchar3[80]='\0';
881  //a contact constraint was found, so append and expand the information
882  for(j=0; j<*nset; j++){
883  memcpy(tchar1,&set[j*81],81);
884  tchar1[80]='\0';
885  if(strcmp(tchar1,tchar2)==0){
886  /* dependent nodal surface was found,copy the original information first */
887  tnstart=tint;
888  for(k=0; k<iendset[j]-istartset[j]+1; k++){
889  ialset[tint-1]=ialset[istartset[j]-1+k];
890  tint++;
891  }
892  /* now append the expanded information */
893  for(l=1; l<*nsectors; l++){
894  for(k=0; k<iendset[j]-istartset[j]+1; k++){
895  ialset[tint-1]=(ialset[istartset[j]-1+k]!=-1)?ialset[istartset[j]-1+k]+*nk*l:-1;
896  tint++;
897  }
898  }
899  tnend=tint-1;
900  /* now replace the information in istartset and iendset*/
901  istartset[j]=tnstart;
902  iendset[j]=tnend;
903  }
904  else if(strcmp(tchar1,tchar3)==0){
905  /* independent element face surface was found */
906  tnstart=tint;
907  for(k=0; k<iendset[j]-istartset[j]+1; k++){
908  ialset[tint-1]=ialset[istartset[j]-1+k];
909  tint++;
910  }
911  /* now append the expanded information*/
912  for(l=1; l<*nsectors; l++){
913  for(k=0; k<iendset[j]-istartset[j]+1; k++){
914  tint2=((ITG)(ialset[istartset[j]-1+k]))/10;
915  ialset[tint-1]=(ialset[istartset[j]-1+k]!=-1)?(tint2+*ne*l)*10+(ialset[istartset[j]-1+k]-(tint2*10)):-1;
916  tint++;
917  }
918  }
919  tnend=tint-1;
920  /* now replace the information in istartset and iendset*/
921  istartset[j]=tnstart;
922  iendset[j]=tnend;
923  }
924  }
925  }
926  }
927  SFREE(tchar1);
928  SFREE(tchar2);
929  SFREE(tchar3);
930  }
931 
932  *nk=nkt;
933  (*ne)*=(*nsectors);
934  (*nkon)*=(*nsectors);
935  (*nboun)*=(*nsectors);
936  neq[1]=neqh**nsectors;
937 
938  *zp=z;*izdofp=izdof;
939 
940  SFREE(temp_array);SFREE(coefmpcnew);SFREE(noderight);SFREE(coefright);
941  SFREE(v);SFREE(vt);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
942  SFREE(inocs);SFREE(ielcs);SFREE(filabt);SFREE(iznode);SFREE(csmass);
943 
944  return;
945 }
#define ITGFORMAT
Definition: CalculiX.h:52
subroutine calcmass(ipkon, lakon, kon, co, mi, nelem, ne, thicke, ielmat, nope, t0, t1, rhcon, nrhcon, ntmat_, ithermal, csmasstot)
Definition: calcmass.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))
ITG strcpy1(char *s1, const char *s2, ITG length)
Definition: strcpy1.c:24
ITG strcmp1(const char *s1, const char *s2)
Definition: strcmp1.c:24
#define DMEMSET(a, b, c, d)
Definition: CalculiX.h:45
subroutine stop()
Definition: stop.f:19
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:19
subroutine rectcylexp(co, v, fn, stn, qfn, een, cs, nkt, icntrl, t, filab, imag, mi, iznode, nznode, nsectors, nk, emn)
Definition: rectcylexp.f:19
#define RENEW(a, b, c)
Definition: CalculiX.h:40
#define SFREE(a)
Definition: CalculiX.h:41
void results(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, double *v, double *stn, ITG *inum, double *stx, double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, double *alcon, ITG *nalcon, double *alzero, ITG *ielmat, ITG *ielorien, ITG *norien, double *orab, ITG *ntmat_, double *t0, double *t1, ITG *ithermal, double *prestr, ITG *iprestr, char *filab, double *eme, double *emn, double *een, ITG *iperturb, double *f, double *fn, ITG *nactdof, ITG *iout, double *qa, double *vold, double *b, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc, ITG *nmpc, ITG *nmethod, double *vmax, ITG *neq, double *veold, double *accold, double *beta, double *gamma, double *dtime, double *time, double *ttime, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double *xstateini, double *xstiff, double *xstate, ITG *npmat_, double *epl, char *matname, ITG *mi, ITG *ielas, ITG *icmd, ITG *ncmat_, ITG *nstate_, double *stiini, double *vini, ITG *ikboun, ITG *ilboun, double *ener, double *enern, double *emeini, double *xstaten, double *eei, double *enerini, double *cocon, ITG *ncocon, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, ITG *nprint, char *prlab, char *prset, double *qfx, double *qfn, double *trab, ITG *inotr, ITG *ntrans, double *fmpc, ITG *nelemload, ITG *nload, ITG *ikmpc, ITG *ilmpc, ITG *istep, ITG *iinc, double *springarea, double *reltime, ITG *ne0, double *xforc, ITG *nforc, double *thicke, double *shcon, ITG *nshcon, char *sideload, double *xload, double *xloadold, ITG *icfd, ITG *inomat, double *pslavsurf, double *pmastsurf, ITG *mortar, ITG *islavact, double *cdn, ITG *islavnode, ITG *nslavnode, ITG *ntie, double *clearini, ITG *islavsurf, ITG *ielprop, double *prop)
Definition: results.c:42
subroutine addizdofcload(nodeforc, ndirforc, nactdof, mi, izdof, nzdof, iforc, iznode, nznode, nk, imdnode, nmdnode, xforc)
Definition: addizdofcload.f:19
subroutine nident(x, px, n, id)
Definition: nident.f:25
subroutine isortiiddc(ix1, ix2, dy1, dy2, cy, n, kflag)
Definition: isortiiddc.f:5
subroutine isortiddc(ix, dy1, dy2, cy, n, kflag)
Definition: isortiddc.f:5
subroutine rectcyltrfm(node, co, cs, icntrl, fin, fout)
Definition: rectcyltrfm.f:19
#define ITG
Definition: CalculiX.h:51
subroutine addizdofdload(nelemload, sideload, ipkon, kon, lakon, nactdof, izdof, nzdof, mi, iload, iznode, nznode, nk, imdnode, nmdnode)
Definition: addizdofdload.f:19
#define NNEW(a, b, c)
Definition: CalculiX.h:39
subroutine rectcyl(co, v, fn, stn, qfn, een, cs, n, icntrl, t, filab, imag, mi, emn)
Definition: rectcyl.f:19
Hosted by OpenAircraft.com, (Michigan UAV, LLC)