CalculiX  2.8
A Free Software Three-Dimensional Structural Finite Element Program
 All Classes Files Functions Variables Macros
matrixstorage.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void matrixstorage (double *ad, double **aup, double *adb, double *aub, double *sigma, ITG *icol, ITG **irowp, ITG *neq, ITG *nzs, ITG *ntrans, ITG *inotr, double *trab, double *co, ITG *nk, ITG *nactdof, char *jobnamec, ITG *mi, ITG *ipkon, char *lakon, ITG *kon, ITG *ne, ITG *mei, ITG *nboun, ITG *nmpc, double *cs, ITG *mcs)
 

Function Documentation

void matrixstorage ( double *  ad,
double **  aup,
double *  adb,
double *  aub,
double *  sigma,
ITG icol,
ITG **  irowp,
ITG neq,
ITG nzs,
ITG ntrans,
ITG inotr,
double *  trab,
double *  co,
ITG nk,
ITG nactdof,
char *  jobnamec,
ITG mi,
ITG ipkon,
char *  lakon,
ITG kon,
ITG ne,
ITG mei,
ITG nboun,
ITG nmpc,
double *  cs,
ITG mcs 
)

Definition at line 25 of file matrixstorage.c.

31  {
32 
33  char fsti[132]="",fmas[132]="",fdof[132]="";
34  ITG i,j,l,*irow=NULL,*ai=NULL,*aj=NULL,kflag=2,ndim,jref,kstart,klen,
35  npoint_,npoint,neq3,index,i3l,i3c,i3lo,i3co,idof,n,il,
36  ic,id,itrans,ndim2,*ipoindex=NULL,mt=mi[1]+1,*nactdofinv=NULL,
37  *nodorig=NULL,inode,idir;
38  long long *ipoint=NULL,k;
39  double *au=NULL,*aa=NULL,*trans=NULL,*aa3=NULL,a[9];
40  FILE *f2,*f3,*f4;
41 
42  strcpy(fsti,jobnamec);
43  strcat(fsti,".sti");
44 
45  printf(" Storing the stiffness matrix in file %s \n\n",fsti);
46 
47  if((*mcs!=0)&&(cs[1]>=0)){
48  printf(" For cyclic symmetric calculations the complex\n");
49  printf(" Hermitian matrix is stored as a symmetric real\n");
50  printf(" matrix double its size; if R stands for the real\n");
51  printf(" part of the matrix and I for the imaginary part,\n");
52  printf(" the resulting matrix takes the form:\n");
53  printf(" _ _\n");
54  printf(" | |\n");
55  printf(" | R -I |\n");
56  printf(" | I R |\n");
57  printf(" |_ _|\n\n");
58  printf(" This applies to the stiffness and the mas matrix\n\n");
59  }
60 
61  if((f2=fopen(fsti,"wb"))==NULL){
62  printf("*ERROR in matrixstorage: cannot open %s for writing...\n",fsti);
63  FORTRAN(stop,());
64  }
65 
66  au=*aup;
67  irow=*irowp;
68 
69  ndim=*neq+*nzs;
70 
71  itrans=0;
72  if(*ntrans!=0){
73  for(i=0;i<*nk;i++){
74  if(inotr[2*i]!=0){
75  itrans=1;
76  break;
77  }
78  }
79  }
80 
81  /* stiffness matrix */
82 
83  if((itrans==0)||(mei[0]==1)){
84 
85  /* no transformation */
86 
87  NNEW(aa,double,ndim);
88  NNEW(ai,ITG,ndim);
89  NNEW(aj,ITG,ndim);
90 
91  k=0;
92  for(i=0;i<*neq;i++){
93  ai[k]=i+1;
94  aj[k]=i+1;
95  aa[k]=ad[i];
96  k++;
97  }
98  l=0;
99  for(i=0;i<*neq;i++){
100  for(j=0;j<icol[i];j++){
101  ai[k]=i+1;
102  aj[k]=irow[l];
103  aa[k]=au[l];
104  k++;l++;
105  }
106  }
107  }
108  else{
109 
110  /* transformation: storing the matrices in transformed coordinates
111  (needed by turboreduce (not part of CalculiX)) */
112 
113  if((*nboun!=0)||(*nmpc!=0)){
114  printf("*ERROR in matrixstorage: matrix storage in local\n");
115  printf(" coordinates is only possible in the absence\n");
116  printf(" of SPC's and MPC's\n\n");
117  FORTRAN(stop,());
118  }
119 
120  ndim2=*neq+2**nzs;
121 
122  /* aa contains the linear storage of the individual matrix elements */
123 
124  NNEW(aa,double,ndim2);
125 
126  /* aa3 contains the linear storage of the 3x3 matrices */
127 
128  NNEW(aa3,double,ndim2);
129  NNEW(ai,ITG,ndim2);
130  NNEW(aj,ITG,ndim2);
131 
132  k=0;
133  for(i=0;i<*neq;i++){
134  ai[k]=i+1;
135  aj[k]=i+1;
136  aa[k]=ad[i];
137  k++;
138  }
139  l=0;
140  for(i=0;i<*neq;i++){
141  for(j=0;j<icol[i];j++){
142  ai[k]=i+1;
143  aj[k]=irow[l];
144  aa[k]=au[l];
145  k++;
146  ai[k]=irow[l];
147  aj[k]=i+1;
148  aa[k]=au[l];
149  k++;
150  l++;
151  }
152  }
153 
154  FORTRAN(isortiid,(aj,ai,aa,&ndim2,&kflag));
155 
156  k=0;
157  for(i=0;i<*neq;i++){
158  jref=aj[k];
159  kstart=k;
160  do{
161  k++;
162  if(k==ndim2) break;
163  if(aj[k]!=jref) break;
164  }while(1);
165  klen=k-kstart;
166  FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
167  }
168 
169  /* npoint is the number of 3x3 matrices
170  ipoint contains the two-dimensional location of the matrices
171  (row and column stored as column x neq3 + row)
172  ipoindex contains the one-dimensional location of the matrices in
173  fields aa3 */
174 
175  npoint_=*neq;
176  npoint=0;
177  NNEW(ipoint,long long,npoint_);
178  NNEW(ipoindex,ITG,npoint_);
179 
180  neq3=*neq/3;
181  index=0;
182 
183  /* storing the matrix as a matrix of 3x3 submatrices */
184 
185  for(i=0;i<ndim2;i++){
186  il=ai[i];
187  ic=aj[i];
188  i3co=(ic-1)%3;
189  i3c=((ic-1)-i3co)/3;
190  i3lo=(il-1)%3;
191  i3l=((il-1)-i3lo)/3;
192  k=(long long)i3c*neq3+i3l;
193  FORTRAN(nidentll,(ipoint,&k,&npoint,&id));
194  if(npoint==0){
195  npoint++;
196  ipoint[npoint-1]=k;
197  ipoindex[npoint-1]=index;
198  }
199  else if(ipoint[id-1]!=k){
200  npoint++;
201  if(npoint>npoint_){
202  npoint_=(ITG)(1.1*npoint_);
203  RENEW(ipoint,long long,npoint_);
204  RENEW(ipoindex,ITG,npoint_);
205  }
206  index+=9;
207  ipoint[npoint-1]=k;
208  ipoindex[npoint-1]=index;
209  }
210  else{
211  index=ipoindex[id-1];
212  }
213  aa3[index+3*i3co+i3lo]=aa[i];
214  }
215 
216  /* defining the transformation matrix (diagonal matrix of
217  3x3 submatrices */
218 
219  NNEW(trans,double,9*neq3);
220  for (i=0;i<*nk;i++){
221  idof=nactdof[mt*i+1];
222  if(idof==0) continue;
223  itrans=inotr[2*i];
224  if(itrans==0){
225  for(j=0;j<9;j++){
226  trans[3*(idof-1)+j]=0.;
227  }
228  trans[3*(idof-1)]=1.;
229  trans[3*(idof-1)+4]=1.;
230  trans[3*(idof-1)+8]=1.;
231  }
232  else{
233  FORTRAN(transformatrix,(&trab[7*itrans-7],&co[3*i],a));
234  for(j=0;j<9;j++){
235  trans[3*(idof-1)+j]=a[j];
236  }
237  }
238  }
239 
240  /* postmultiplying the matrix with the transpose of the
241  transformation matrix */
242 
243  for(i=0;i<npoint;i++){
244  i3l=ipoint[i]%neq3;
245  i3c=(ipoint[i]-i3l)/neq3;
246  n=2;
247  FORTRAN(mult,(&aa3[9*i],&trans[9*i3c],&n));
248  }
249 
250  /* premultiplying the matrix with the transformation matrix */
251 
252  for(i=0;i<npoint;i++){
253  i3l=ipoint[i]%neq3;
254  n=1;
255  FORTRAN(mult,(&aa3[9*i],&trans[9*i3l],&n));
256  }
257 
258  /* storing the new matrix in a linear format */
259 
260  k=-1;
261  for(i=0;i<npoint;i++){
262  i3l=ipoint[i]%neq3;
263  i3c=(ipoint[i]-i3l)/neq3;
264  for(j=0;j<9;j++){
265  i3lo=j%3;
266  i3co=(j-i3lo)/3;
267  ic=3*i3c+i3co+1;
268  il=3*i3l+i3lo+1;
269  if(il>ic) continue;
270  k++;
271  ai[k]=il;
272  aj[k]=ic;
273  aa[k]=aa3[9*i+j];
274  }
275  }
276  SFREE(aa3);SFREE(ipoint);SFREE(ipoindex);SFREE(trans);
277  }
278 
279  FORTRAN(isortiid,(aj,ai,aa,&ndim,&kflag));
280 
281  k=0;
282  for(i=0;i<*neq;i++){
283  jref=aj[k];
284  kstart=k;
285  do{
286  k++;
287  if(aj[k]!=jref) break;
288  }while(1);
289  klen=k-kstart;
290  FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
291  }
292 
293  for(i=0;i<ndim;i++){
294  fprintf(f2,"%" ITGFORMAT " %" ITGFORMAT " %20.13e\n",ai[i],aj[i],aa[i]);
295  }
296 
297  fclose(f2);
298 
299  SFREE(ai);SFREE(aj);SFREE(aa);
300 
301  /* mass matrix */
302 
303  strcpy(fmas,jobnamec);
304  strcat(fmas,".mas");
305 
306  printf(" Storing the mass matrix in file %s \n\n",fmas);
307 
308  if((f3=fopen(fmas,"wb"))==NULL){
309  printf("*ERROR in matrixstorage: cannot open %s for writing...\n",fmas);
310  FORTRAN(stop,());
311  }
312 
313  if((itrans==0)||(mei[0]==1)){
314 
315  /* no transformation */
316 
317  NNEW(aa,double,ndim);
318  NNEW(ai,ITG,ndim);
319  NNEW(aj,ITG,ndim);
320 
321  k=0;
322  for(i=0;i<*neq;i++){
323  ai[k]=i+1;
324  aj[k]=i+1;
325  aa[k]=adb[i];
326  k++;
327  }
328  l=0;
329  for(i=0;i<*neq;i++){
330  for(j=0;j<icol[i];j++){
331  ai[k]=i+1;
332  aj[k]=irow[l];
333  aa[k]=aub[l];
334  k++;l++;
335  }
336  }
337  }
338  else{
339 
340 
341  /* transformation: storing the matrices in transformed coordinates
342  (needed by turboreduce (not part of CalculiX)) */
343 
344  ndim2=*neq+2**nzs;
345 
346  /* aa contains the linear storage of the individual matrix elements */
347 
348  NNEW(aa,double,ndim2);
349 
350  /* aa3 contains the linear storage of the 3x3 matrices */
351 
352  NNEW(aa3,double,ndim2);
353  NNEW(ai,ITG,ndim2);
354  NNEW(aj,ITG,ndim2);
355 
356  k=0;
357  for(i=0;i<*neq;i++){
358  ai[k]=i+1;
359  aj[k]=i+1;
360  aa[k]=adb[i];
361  k++;
362  }
363  l=0;
364  for(i=0;i<*neq;i++){
365  for(j=0;j<icol[i];j++){
366  ai[k]=i+1;
367  aj[k]=irow[l];
368  aa[k]=aub[l];
369  k++;
370  ai[k]=irow[l];
371  aj[k]=i+1;
372  aa[k]=aub[l];
373  k++;
374  l++;
375  }
376  }
377 
378  FORTRAN(isortiid,(aj,ai,aa,&ndim2,&kflag));
379 
380  k=0;
381  for(i=0;i<*neq;i++){
382  jref=aj[k];
383  kstart=k;
384  do{
385  k++;
386  if(k==ndim2) break;
387  if(aj[k]!=jref) break;
388  }while(1);
389  klen=k-kstart;
390  FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
391  }
392 
393  /* npoint is the number of 3x3 matrices
394  ipoint contains the two-dimensional location of the matrices
395  (row and column stored as column x neq3 + row)
396  ipoindex contains the one-dimensional location of the matrices in
397  fields aa3 */
398 
399  npoint_=*neq;
400  npoint=0;
401  NNEW(ipoint,long long,npoint_);
402  NNEW(ipoindex,ITG,npoint_);
403 
404  neq3=*neq/3;
405  index=0;
406 
407  /* storing the matrix as a matrix of 3x3 submatrices */
408 
409  for(i=0;i<ndim2;i++){
410  il=ai[i];
411  ic=aj[i];
412  i3co=(ic-1)%3;
413  i3c=((ic-1)-i3co)/3;
414  i3lo=(il-1)%3;
415  i3l=((il-1)-i3lo)/3;
416  k=(long long)i3c*neq3+i3l;
417  FORTRAN(nidentll,(ipoint,&k,&npoint,&id));
418  if(npoint==0){
419  npoint++;
420  ipoint[npoint-1]=k;
421  ipoindex[npoint-1]=index;
422  }
423  else if(ipoint[id-1]!=k){
424  npoint++;
425  if(npoint>npoint_){
426  npoint_=(ITG)(1.1*npoint_);
427  RENEW(ipoint,long long,npoint_);
428  RENEW(ipoindex,ITG,npoint_);
429  }
430  index+=9;
431  ipoint[npoint-1]=k;
432  ipoindex[npoint-1]=index;
433  }
434  else{
435  index=ipoindex[id-1];
436  }
437  aa3[index+3*i3co+i3lo]=aa[i];
438  }
439 
440  /* defining the transformation matrix (diagonal matrix of
441  3x3 submatrices */
442 
443  NNEW(trans,double,9*neq3);
444  for (i=0;i<*nk;i++){
445  idof=nactdof[mt*i+1];
446  if(idof==0) continue;
447  itrans=inotr[2*i];
448  if(itrans==0){
449  for(j=0;j<9;j++){
450  trans[3*(idof-1)+j]=0.;
451  }
452  trans[3*(idof-1)]=1.;
453  trans[3*(idof-1)+4]=1.;
454  trans[3*(idof-1)+8]=1.;
455  }
456  else{
457  FORTRAN(transformatrix,(&trab[7*itrans-7],&co[3*i],a));
458  for(j=0;j<9;j++){
459  trans[3*(idof-1)+j]=a[j];
460  }
461  }
462  }
463 
464  /* postmultiplying the matrix with the transpose of the
465  transformation matrix */
466 
467  for(i=0;i<npoint;i++){
468  i3l=ipoint[i]%neq3;
469  i3c=(ipoint[i]-i3l)/neq3;
470  n=2;
471  FORTRAN(mult,(&aa3[9*i],&trans[9*i3c],&n));
472  }
473 
474  /* premultiplying the matrix with the transformation matrix */
475 
476  for(i=0;i<npoint;i++){
477  i3l=ipoint[i]%neq3;
478  n=1;
479  FORTRAN(mult,(&aa3[9*i],&trans[9*i3l],&n));
480  }
481 
482  /* storing the new matrix in a linear format */
483 
484  k=-1;
485  for(i=0;i<npoint;i++){
486  i3l=ipoint[i]%neq3;
487  i3c=(ipoint[i]-i3l)/neq3;
488  for(j=0;j<9;j++){
489  i3lo=j%3;
490  i3co=(j-i3lo)/3;
491  ic=3*i3c+i3co+1;
492  il=3*i3l+i3lo+1;
493  if(il>ic) continue;
494  k++;
495  ai[k]=il;
496  aj[k]=ic;
497  aa[k]=aa3[9*i+j];
498  }
499  }
500  SFREE(aa3);SFREE(ipoint);SFREE(ipoindex);SFREE(trans);
501  }
502 
503  FORTRAN(isortiid,(aj,ai,aa,&ndim,&kflag));
504 
505  k=0;
506  for(i=0;i<*neq;i++){
507  jref=aj[k];
508  kstart=k;
509  do{
510  k++;
511  if(aj[k]!=jref) break;
512  }while(1);
513  klen=k-kstart;
514  FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
515  }
516 
517  for(i=0;i<ndim;i++){
518  fprintf(f3,"%" ITGFORMAT " %" ITGFORMAT " %20.13e\n",ai[i],aj[i],aa[i]);
519  }
520 
521  fclose(f3);
522 
523  SFREE(ai);SFREE(aj);SFREE(aa);
524 
525  *aup=au;
526  *irowp=irow;
527 
528  /* node and global direction for each degree of freedom in the
529  matrix (corresponding with a row or column number) */
530 
531  strcpy(fdof,jobnamec);
532  strcat(fdof,".dof");
533 
534  if((itrans==0)||(mei[0]==1)){
535  printf(" Storing the node and global direction per entry (row or column)\n in the stiffness and mass matrices in the form node.direction in file %s \n\n",fdof);
536  }else{
537  printf(" Storing the node and local direction per entry (row or column)\n in the stiffness and mass matrices in the form node.direction in file %s \n\n",fdof);
538  }
539 
540  if((f4=fopen(fdof,"wb"))==NULL){
541  printf("*ERROR in matrixstorage: cannot open %s for writing...\n",fdof);
542  FORTRAN(stop,());
543  }
544 
545  /* invert nactdof */
546 
547  NNEW(nactdofinv,ITG,mt**nk);
548  NNEW(nodorig,ITG,*nk);
549  FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
550  ipkon,lakon,kon,ne));
551  SFREE(nodorig);
552 
553  if((*mcs==0)||(cs[1]<0)){
554  for(i=0;i<*neq;i++){
555  inode=(ITG)((double)nactdofinv[(ITG)i]/mt)+1;
556  idir=nactdofinv[(ITG)i]-mt*(inode-1);
557  fprintf(f4,"%" ITGFORMAT ".%" ITGFORMAT "\n",inode,idir);
558  }
559  }else{
560  for(i=0;i<*neq/2;i++){
561  inode=(ITG)((double)nactdofinv[(ITG)i]/mt)+1;
562  idir=nactdofinv[(ITG)i]-mt*(inode-1);
563  fprintf(f4,"%" ITGFORMAT ".%" ITGFORMAT "\n",inode,idir);
564  }
565  }
566 
567  fclose(f4);
568 
569  FORTRAN(stop,());
570 
571  return;
572 }
subroutine nidentll(x, px, n, id)
Definition: nidentll.f:27
#define ITGFORMAT
Definition: CalculiX.h:52
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))
subroutine stop()
Definition: stop.f:19
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:19
#define RENEW(a, b, c)
Definition: CalculiX.h:40
#define SFREE(a)
Definition: CalculiX.h:41
subroutine gennactdofinv(nactdof, nactdofinv, nk, mi, nodorig, ipkon, lakon, kon, ne)
Definition: gennactdofinv.f:19
#define ITG
Definition: CalculiX.h:51
subroutine mult(matrix, trans, n)
Definition: mult.f:19
subroutine isortiid(ix, iy, dy, n, kflag)
Definition: isortiid.f:5
#define NNEW(a, b, c)
Definition: CalculiX.h:39
Hosted by OpenAircraft.com, (Michigan UAV, LLC)