CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions
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, char *orname)
 

Function Documentation

◆ expand()

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