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

Go to the source code of this file.

Macros

#define min(a, b)   ((a) <= (b) ? (a) : (b))
 
#define max(a, b)   ((a) >= (b) ? (a) : (b))
 

Functions

void cascade (ITG *ipompc, double **coefmpcp, ITG **nodempcp, ITG *nmpc, ITG *mpcfree, ITG *nodeboun, ITG *ndirboun, ITG *nboun, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun, ITG *mpcend, ITG *mpcmult, char *labmpc, ITG *nk, ITG *memmpc_, ITG *icascade, ITG *maxlenmpc, ITG *callfrommain, ITG *iperturb, ITG *ithermal)
 

Macro Definition Documentation

#define max (   a,
 
)    ((a) >= (b) ? (a) : (b))

Definition at line 32 of file cascade.c.

#define min (   a,
 
)    ((a) <= (b) ? (a) : (b))

Definition at line 31 of file cascade.c.

Function Documentation

void cascade ( ITG ipompc,
double **  coefmpcp,
ITG **  nodempcp,
ITG nmpc,
ITG mpcfree,
ITG nodeboun,
ITG ndirboun,
ITG nboun,
ITG ikmpc,
ITG ilmpc,
ITG ikboun,
ITG ilboun,
ITG mpcend,
ITG mpcmult,
char *  labmpc,
ITG nk,
ITG memmpc_,
ITG icascade,
ITG maxlenmpc,
ITG callfrommain,
ITG iperturb,
ITG ithermal 
)

Definition at line 34 of file cascade.c.

38  {
39 
40  /* detects cascaded mpc's and decascades them; checks multiple
41  occurrence of the same dependent DOF's in different mpc/spc's
42 
43  data structure of ipompc,coefmpc,nodempc:
44  for each mpc, e.g. i,
45  -the nodes are stored in nodempc(1,ipompc(i)),
46  nodempc(1,nodempc(3,ipompc(i))),
47  nodempc(1,nodempc(3,nodempc(3,ipompc(i))))... till
48  nodempc(3,nodempc(3,nodempc(3,.......))))))=0;
49  -the corresponding directions in nodempc(2,ipompc(i)),
50  nodempc(2,nodempc(3,ipompc(i))),.....
51  -the corresponding coefficient in coefmpc(ipompc(i)),
52  coefmpc(nodempc(3,ipompc(i))),.....
53  the mpc is written as a(1)u(i1,j1)+a(2)u(i2,j2)+...
54  +....a(k)u(ik,jk)=0, the first term is the dependent term,
55  the others are independent, at least after execution of the
56  present routine. The mpc's must be homogeneous, otherwise a
57  error message is generated and the program stops. */
58 
59  ITG i,j,index,id,idof,nterm,idepend,*nodempc=NULL,
60  ispooles,iexpand,ichange,indexold,ifluidmpc,
61  mpc,indexnew,index1,index2,index1old,index2old,*jmpc=NULL,nl;
62 
63  double coef,*coefmpc=NULL;
64 
65 #ifdef SPOOLES
66 
67  ITG irow,icolumn,node,idir,irownl,icolnl,*ipointer=NULL,*icoef=NULL,
68  ifree,*indepdof=NULL,nindep;
69 
70  double *xcoef=NULL,b;
71 
72  DenseMtx *mtxB, *mtxX ;
73  Chv *rootchv ;
74  ChvManager *chvmanager ;
75  SubMtxManager *mtxmanager ;
76  FrontMtx *frontmtx ;
77  InpMtx *mtxA ;
78  double tau = 100.;
79  double cpus[10] ;
80  ETree *frontETree ;
81  FILE *msgFile ;
82  Graph *graph ;
83  ITG jrhs, msglvl=0, nedges,error,
84  nent, neqns, nrhs, pivotingflag=1, seed=389,
85  symmetryflag=2, type=1,maxdomainsize,maxzeros,maxsize;
86  ITG *oldToNew ;
87  ITG stats[20] ;
88  IV *newToOldIV, *oldToNewIV ;
89  IVL *adjIVL, *symbfacIVL ;
90 #endif
91 
92  nodempc=*nodempcp;
93  coefmpc=*coefmpcp;
94 
95  /* for(i=0;i<*nmpc;i++){
96  j=i+1;
97  FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
98  }*/
99 
100  NNEW(jmpc,ITG,*nmpc);
101  idepend=0;
102 
103 /* check whether a node is used as a dependent node in a MPC
104  and in a SPC */
105 
106  for(i=0;i<*nmpc;i++){
107  if(*nboun>0){
108  FORTRAN(nident,(ikboun,&ikmpc[i],nboun,&id));}
109  else{id=0;}
110  if(id>0){
111  if(ikboun[id-1]==ikmpc[i]){
112  if(strcmp1(&labmpc[20*i],"FLUID")!=0){
113  printf("*ERROR in cascade: the DOF corresponding to \n node %" ITGFORMAT " in direction %" ITGFORMAT " is detected on the \n dependent side of a MPC and a SPC\n",
114  (ikmpc[i])/8+1,ikmpc[i]-8*((ikmpc[i])/8));
115  }else{
116  printf("*ERROR in cascade: the DOF corresponding to \n face %" ITGFORMAT " in direction %" ITGFORMAT " is detected on the \n dependent side of a MPC and a SPC\n",
117  (-ikmpc[i])/8+1,-ikmpc[i]-8*((-ikmpc[i])/8));
118  }
119  FORTRAN(stop,());
120  }
121  }
122  }
123 
124 /* check whether there are user mpc's: in user MPC's the
125  dependent DOF can change, however, the number of terms
126  cannot change */
127 
128  for(i=0;i<*nmpc;i++){
129 
130  /* linear mpc */
131 
132  /* because of the next line the size of field labmpc
133  has to be defined as 20*nmpc+1: without "+1" an
134  undefined field is accessed */
135 
136  if((strcmp1(&labmpc[20*i]," ")==0) ||
137  (strcmp1(&labmpc[20*i],"CYCLIC")==0) ||
138  (strcmp1(&labmpc[20*i],"SUBCYCLIC")==0)||
139  (strcmp1(&labmpc[20*i],"CONTACT")==0)||
140  (strcmp1(&labmpc[20*i],"FLUID")==0)||
141  (*iperturb<2)) jmpc[i]=0;
142 
143  /* nonlinear mpc */
144 
145  else if((strcmp1(&labmpc[20*i],"RIGID")==0) ||
146  (strcmp1(&labmpc[20*i],"KNOT")==0) ||
147  (strcmp1(&labmpc[20*i],"PLANE")==0) ||
148  (strcmp1(&labmpc[20*i],"STRAIGHT")==0)||
149  (strcmp1(&labmpc[20*i],"ISOCHORIC")==0)) jmpc[i]=1;
150 
151  /* user mpc */
152 
153  else{
154  jmpc[i]=1;
155  if(*icascade==0) *icascade=1;
156  }
157  }
158 
159 /* decascading */
160 
161  ispooles=0;
162 
163  /* decascading using simple substitution */
164 
165  do{
166  ichange=0;
167  for(i=0;i<*nmpc;i++){
168 
169  if(strcmp1(&labmpc[20*i],"FLUID")!=0){
170  ifluidmpc=0;
171  }else{
172  ifluidmpc=1;
173  }
174 
175  if(jmpc[i]==1) nl=1;
176  else nl=0;
177  iexpand=0;
178  index=nodempc[3*ipompc[i]-1];
179  if(index==0) continue;
180  do{
181  if(ifluidmpc==0){
182  /* MPC on node */
183  idof=(nodempc[3*index-3]-1)*8+nodempc[3*index-2];
184  }else{
185  if(nodempc[3*index-3]>0){
186  /* MPC on face */
187  idof=-((nodempc[3*index-3]-1)*8+nodempc[3*index-2]);
188  }else{
189  /* MPC on node
190  SPC number: -nodempc[3*index-3]
191  node: nodeboun[-nodempc[3*index-3]-1] */
192 
193  idof=(nodeboun[-nodempc[3*index-3]-1]-1)*8+nodempc[3*index-2];
194  }
195  }
196 
197  FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
198  if((id>0)&&(ikmpc[id-1]==idof)){
199 
200  /* a term on the independent side of the MPC is
201  detected as dependent node in another MPC */
202 
203  indexold=nodempc[3*index-1];
204  coef=coefmpc[index-1];
205  mpc=ilmpc[id-1];
206 
207  /* no expansion if there is a dependence of a
208  nonlinear MPC on another linear or nonlinear MPC
209  and the call is from main */
210 
211  if((jmpc[mpc-1]==1)||(nl==1)){
212  *icascade=2;
213  if(idepend==0){
214  printf("*INFO in cascade: linear MPCs and\n");
215  printf(" nonlinear MPCs depend on each other\n");
216  printf(" common node: %" ITGFORMAT " in direction %" ITGFORMAT "\n\n",nodempc[3*index-3],nodempc[3*index-2]);
217  idepend=1;}
218  if(*callfrommain==1){
219  index=nodempc[3*index-1];
220  if(index!=0) continue;
221  else break;}
222  }
223 
224 /* printf("*INFO in cascade: DOF %" ITGFORMAT " of node %" ITGFORMAT " is expanded\n",
225  nodempc[3*index-2],nodempc[3*index-3]);*/
226 
227  /* collecting terms corresponding to the same DOF */
228 
229  index1=ipompc[i];
230  do{
231  index2old=index1;
232  index2=nodempc[3*index1-1];
233  if(index2==0) break;
234  do{
235  if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
236  (nodempc[3*index1-2]==nodempc[3*index2-2])){
237  coefmpc[index1-1]+=coefmpc[index2-1];
238  nodempc[3*index2old-1]=nodempc[3*index2-1];
239  nodempc[3*index2-1]=*mpcfree;
240  *mpcfree=index2;
241  index2=nodempc[3*index2old-1];
242  if(index2==0) break;
243  }
244  else{
245  index2old=index2;
246  index2=nodempc[3*index2-1];
247  if(index2==0) break;
248  }
249  }while(1);
250  index1=nodempc[3*index1-1];
251  if(index1==0) break;
252  }while(1);
253 
254  /* check for zero coefficients on the dependent side */
255 
256  index1=ipompc[i];
257  if(fabs(coefmpc[index1-1])<1.e-10){
258  printf("*ERROR in cascade: zero coefficient on the\n");
259  printf(" dependent side of an equation\n");
260  printf(" dependent node: %" ITGFORMAT "",nodempc[3*index1-3]);
261  printf(" direction: %" ITGFORMAT "",nodempc[3*index1-2]);
262  FORTRAN(stop,());
263  }
264 
265  ichange=1;iexpand=1;
266  if((strcmp1(&labmpc[20*i]," ")==0)&&
267  (strcmp1(&labmpc[20*(mpc-1)],"CYCLIC")==0))
268  strcpy1(&labmpc[20*i],"SUBCYCLIC",9);
269  indexnew=ipompc[mpc-1];
270  coef=-coef/coefmpc[indexnew-1];
271  indexnew=nodempc[3*indexnew-1];
272  do{
273  coefmpc[index-1]=coef*coefmpc[indexnew-1];
274  nodempc[3*index-3]=nodempc[3*indexnew-3];
275  nodempc[3*index-2]=nodempc[3*indexnew-2];
276  indexnew=nodempc[3*indexnew-1];
277  if(indexnew!=0){
278  nodempc[3*index-1]=*mpcfree;
279  index=*mpcfree;
280  *mpcfree=nodempc[3**mpcfree-1];
281  if(*mpcfree==0){
282  *mpcfree=*memmpc_+1;
283  nodempc[3*index-1]=*mpcfree;
284  *memmpc_=(ITG)(1.1**memmpc_);
285  printf("*INFO in cascade: reallocating nodempc; new size = %" ITGFORMAT "\n\n",*memmpc_);
286  RENEW(nodempc,ITG,3**memmpc_);
287  RENEW(coefmpc,double,*memmpc_);
288  for(j=*mpcfree;j<*memmpc_;j++){
289  nodempc[3*j-1]=j+1;
290  }
291  nodempc[3**memmpc_-1]=0;
292  }
293  continue;
294  }
295  else{
296  nodempc[3*index-1]=indexold;
297  break;
298  }
299  }while(1);
300  break;
301  }
302  else{
303  index=nodempc[3*index-1];
304  if(index!=0) continue;
305  else break;
306  }
307  }while(1);
308  if(iexpand==0) continue;
309 
310  /* one term of the mpc was expanded
311  collecting terms corresponding to the same DOF */
312 
313  index1=ipompc[i];
314  do{
315  index2old=index1;
316  index2=nodempc[3*index1-1];
317  if(index2==0) break;
318  do{
319  if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
320  (nodempc[3*index1-2]==nodempc[3*index2-2])){
321  coefmpc[index1-1]+=coefmpc[index2-1];
322  nodempc[3*index2old-1]=nodempc[3*index2-1];
323  nodempc[3*index2-1]=*mpcfree;
324  *mpcfree=index2;
325  index2=nodempc[3*index2old-1];
326  if(index2==0) break;
327  }
328  else{
329  index2old=index2;
330  index2=nodempc[3*index2-1];
331  if(index2==0) break;
332  }
333  }while(1);
334  index1=nodempc[3*index1-1];
335  if(index1==0) break;
336  }while(1);
337 
338  /* check for zero coefficients on the dependent and
339  independent side */
340 
341  index1=ipompc[i];
342  index1old=0;
343  do {
344  if(fabs(coefmpc[index1-1])<1.e-10){
345  if(index1old==0){
346  printf("*ERROR in cascade: zero coefficient on the\n");
347  printf(" dependent side of an equation\n");
348  printf(" dependent node: %" ITGFORMAT "",nodempc[3*index1-3]);
349  printf(" direction: %" ITGFORMAT "",nodempc[3*index1-2]);
350  FORTRAN(stop,());
351  }
352  else{
353  nodempc[3*index1old-1]=nodempc[3*index1-1];
354  nodempc[3*index1-1]=*mpcfree;
355  *mpcfree=index1;
356  index1=nodempc[3*index1old-1];
357  }
358  }
359  else{
360  index1old=index1;
361  index1=nodempc[3*index1-1];
362  }
363  if(index1==0) break;
364  }while(1);
365  }
366  if(ichange==0) break;
367  }while(1);
368 
369  /* decascading using spooles */
370 
371 #ifdef SPOOLES
372  if((*icascade==1)&&(ispooles==1)){
373  if ( (msgFile = fopen("spooles.out", "a")) == NULL ) {
374  fprintf(stderr, "\n fatal error in spooles.c"
375  "\n unable to open file spooles.out\n") ;
376  }
377  NNEW(ipointer,ITG,7**nk);
378  NNEW(indepdof,ITG,7**nk);
379  NNEW(icoef,ITG,2**memmpc_);
380  NNEW(xcoef,double,*memmpc_);
381  ifree=0;
382  nindep=0;
383 
384  for(i=*nmpc-1;i>-1;i--){
385  index=ipompc[i];
386  while(1){
387  idof=8*(nodempc[3*index-3]-1)+nodempc[3*index-2]-1;
388 
389 /* check whether idof is a independent dof which has not yet been
390  stored in indepdof */
391 
392  FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
393  if((id==0)||(ikmpc[id-1]!=idof)){
394  FORTRAN(nident,(indepdof,&idof,&nindep,&id));
395  if((id==0)||(indepdof[id-1]!=idof)){
396  for(j=nindep;j>id;j--){
397  indepdof[j]=indepdof[j-1];
398  }
399  indepdof[id]=idof;
400  nindep++;
401  }
402  }
403 
404  icoef[2*ifree]=i+1;
405  icoef[2*ifree+1]=ipointer[idof];
406  xcoef[ifree]=coefmpc[index-1];
407  ipointer[idof]=++ifree;
408  index=nodempc[3*index-1];
409  if(index==0) break;
410  }
411  }
412 
413 /* filling the left hand side */
414 
415  nent=*memmpc_;
416  neqns=*nmpc;
417  mtxA = InpMtx_new() ;
418  InpMtx_init(mtxA, INPMTX_BY_ROWS, type, nent, neqns) ;
419 
420  for(i=0;i<*nmpc;i++){
421  idof=ikmpc[i];
422  icolumn=ilmpc[i]-1;
423  if(strcmp1(&labmpc[20*icolumn],"RIGID")==0) icolnl=1;
424  else icolnl=0;
425  index=ipointer[idof-1];
426  while(1){
427  irow=icoef[2*index-2]-1;
428  if(irow!=icolumn){
429  if(strcmp1(&labmpc[20*irow],"RIGID")==0)irownl=1;
430  else irownl=0;
431  if((irownl==1)||(icolnl==1)){
432  *icascade=2;
433  InpMtx_free(mtxA);
434  printf("*ERROR in cascade: linear and nonlinear MPCs depend on each other");
435  FORTRAN(stop,());
436  }
437  }
438  if((strcmp1(&labmpc[20*irow]," ")==0)&&
439  (strcmp1(&labmpc[20*icolumn],"CYCLIC")==0)){
440  strcpy1(&labmpc[20*irow],"SUBCYCLIC",9);}
441  coef=xcoef[index-1];
442  InpMtx_inputRealEntry(mtxA,irow,icolumn,coef);
443  index=icoef[2*index-1];
444  if(index==0) break;
445  }
446  ipointer[idof-1]=0;
447  }
448 
449  InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
450  if ( msglvl > 1 ) {
451  fprintf(msgFile, "\n\n input matrix") ;
452  InpMtx_writeForHumanEye(mtxA, msgFile) ;
453  fflush(msgFile) ;
454  }
455 /*--------------------------------------------------------------------*/
456 /*
457  -------------------------------------------------
458  STEP 2 : find a low-fill ordering
459  (1) create the Graph object
460  (2) order the graph using multiple minimum degree
461  -------------------------------------------------
462 */
463  graph = Graph_new() ;
464  adjIVL = InpMtx_fullAdjacency(mtxA) ;
465  nedges = IVL_tsize(adjIVL) ;
466  Graph_init2(graph, 0, neqns, 0, nedges, neqns, nedges, adjIVL,
467  NULL, NULL) ;
468  if ( msglvl > 1 ) {
469  fprintf(msgFile, "\n\n graph of the input matrix") ;
470  Graph_writeForHumanEye(graph, msgFile) ;
471  fflush(msgFile) ;
472  }
473  maxdomainsize=800;maxzeros=1000;maxsize=64;
474  /*maxdomainsize=neqns/100;*/
475  /*frontETree = orderViaMMD(graph, seed, msglvl, msgFile) ;*/
476  /*frontETree = orderViaND(graph,maxdomainsize,seed,msglvl,msgFile); */
477  /*frontETree = orderViaMS(graph,maxdomainsize,seed,msglvl,msgFile);*/
478  frontETree=orderViaBestOfNDandMS(graph,maxdomainsize,maxzeros,
479  maxsize,seed,msglvl,msgFile);
480  if ( msglvl > 1 ) {
481  fprintf(msgFile, "\n\n front tree from ordering") ;
482  ETree_writeForHumanEye(frontETree, msgFile) ;
483  fflush(msgFile) ;
484  }
485 /*--------------------------------------------------------------------*/
486 /*
487  -----------------------------------------------------
488  STEP 3: get the permutation, permute the matrix and
489  front tree and get the symbolic factorization
490  -----------------------------------------------------
491 */
492  oldToNewIV = ETree_oldToNewVtxPerm(frontETree) ;
493  oldToNew = IV_entries(oldToNewIV) ;
494  newToOldIV = ETree_newToOldVtxPerm(frontETree) ;
495  ETree_permuteVertices(frontETree, oldToNewIV) ;
496  InpMtx_permute(mtxA, oldToNew, oldToNew) ;
497 /* InpMtx_mapToUpperTriangle(mtxA) ;*/
498  InpMtx_changeCoordType(mtxA,INPMTX_BY_CHEVRONS);
499  InpMtx_changeStorageMode(mtxA,INPMTX_BY_VECTORS);
500  symbfacIVL = SymbFac_initFromInpMtx(frontETree, mtxA) ;
501  if ( msglvl > 1 ) {
502  fprintf(msgFile, "\n\n old-to-new permutation vector") ;
503  IV_writeForHumanEye(oldToNewIV, msgFile) ;
504  fprintf(msgFile, "\n\n new-to-old permutation vector") ;
505  IV_writeForHumanEye(newToOldIV, msgFile) ;
506  fprintf(msgFile, "\n\n front tree after permutation") ;
507  ETree_writeForHumanEye(frontETree, msgFile) ;
508  fprintf(msgFile, "\n\n input matrix after permutation") ;
509  InpMtx_writeForHumanEye(mtxA, msgFile) ;
510  fprintf(msgFile, "\n\n symbolic factorization") ;
511  IVL_writeForHumanEye(symbfacIVL, msgFile) ;
512  fflush(msgFile) ;
513  }
514 /*--------------------------------------------------------------------*/
515 /*
516  ------------------------------------------
517  STEP 4: initialize the front matrix object
518  ------------------------------------------
519 */
520  frontmtx = FrontMtx_new() ;
521  mtxmanager = SubMtxManager_new() ;
522  SubMtxManager_init(mtxmanager, NO_LOCK, 0) ;
523  FrontMtx_init(frontmtx, frontETree, symbfacIVL, type, symmetryflag,
524  FRONTMTX_DENSE_FRONTS, pivotingflag, NO_LOCK, 0, NULL,
525  mtxmanager, msglvl, msgFile) ;
526 /*--------------------------------------------------------------------*/
527 /*
528  -----------------------------------------
529  STEP 5: compute the numeric factorization
530  -----------------------------------------
531 */
532  chvmanager = ChvManager_new() ;
533  ChvManager_init(chvmanager, NO_LOCK, 1) ;
534  DVfill(10, cpus, 0.0) ;
535  IVfill(20, stats, 0) ;
536  rootchv = FrontMtx_factorInpMtx(frontmtx, mtxA, tau, 0.0, chvmanager,
537  &error,cpus, stats, msglvl, msgFile) ;
538  ChvManager_free(chvmanager) ;
539  if ( msglvl > 1 ) {
540  fprintf(msgFile, "\n\n factor matrix") ;
541  FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
542  fflush(msgFile) ;
543  }
544  if ( rootchv != NULL ) {
545  fprintf(msgFile, "\n\n matrix found to be singular\n") ;
546  exit(-1) ;
547  }
548  if(error>=0){
549  fprintf(msgFile,"\n\nerror encountered at front %" ITGFORMAT "",error);
550  exit(-1);
551  }
552 /*--------------------------------------------------------------------*/
553 /*
554  --------------------------------------
555  STEP 6: post-process the factorization
556  --------------------------------------
557 */
558  FrontMtx_postProcess(frontmtx, msglvl, msgFile) ;
559  if ( msglvl > 1 ) {
560  fprintf(msgFile, "\n\n factor matrix after post-processing") ;
561  FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
562  fflush(msgFile) ;
563  }
564 
565 /* reinitialize nodempc */
566 
567  *mpcfree=1;
568  for(j=0;j<*nmpc;j++){
569  ipompc[j]=0;}
570 
571 /* filling the RHS */
572 
573  jrhs=0;
574  nrhs=1;
575  mtxB=DenseMtx_new();
576  mtxX=DenseMtx_new();
577 
578  for(i=nindep;i>0;i--){
579  idof=indepdof[i-1];
580  if(ipointer[idof]>0){
581 
582 /* new RHS column */
583 
584  DenseMtx_init(mtxB, type, 0, 0, neqns, nrhs, 1, neqns) ;
585  DenseMtx_zero(mtxB) ;
586 
587  index=ipointer[idof];
588  while(1){
589  irow=icoef[2*index-2]-1;
590  coef=xcoef[index-1];
591  DenseMtx_setRealEntry(mtxB,irow,jrhs,coef);
592  index=icoef[2*index-1];
593  if(index==0) break;
594  }
595 
596  if ( msglvl > 1 ) {
597  fprintf(msgFile, "\n\n rhs matrix in original ordering") ;
598  DenseMtx_writeForHumanEye(mtxB, msgFile) ;
599  fflush(msgFile) ;
600  }
601 
602 /*--------------------------------------------------------------------*/
603 /*
604  ---------------------------------------------------------
605  STEP 8: permute the right hand side into the new ordering
606  ---------------------------------------------------------
607 */
608  DenseMtx_permuteRows(mtxB, oldToNewIV) ;
609  if ( msglvl > 1 ) {
610  fprintf(msgFile, "\n\n right hand side matrix in new ordering") ;
611  DenseMtx_writeForHumanEye(mtxB, msgFile) ;
612  fflush(msgFile) ;
613  }
614 /*--------------------------------------------------------------------*/
615 /*
616  -------------------------------
617  STEP 9: solve the linear system
618  -------------------------------
619 */
620  DenseMtx_init(mtxX, type, 0, 0, neqns, nrhs, 1, neqns) ;
621  DenseMtx_zero(mtxX) ;
622  FrontMtx_solve(frontmtx, mtxX, mtxB, mtxmanager,cpus, msglvl, msgFile) ;
623  if ( msglvl > 1 ) {
624  fprintf(msgFile, "\n\n solution matrix in new ordering") ;
625  DenseMtx_writeForHumanEye(mtxX, msgFile) ;
626  fflush(msgFile) ;
627  }
628 /*--------------------------------------------------------------------*/
629 /*
630  --------------------------------------------------------
631  STEP 10: permute the solution into the original ordering
632  --------------------------------------------------------
633 */
634  DenseMtx_permuteRows(mtxX, newToOldIV) ;
635  if ( msglvl > 1 ) {
636  fprintf(msgFile, "\n\n solution matrix in original ordering") ;
637  DenseMtx_writeForHumanEye(mtxX, msgFile) ;
638  fflush(msgFile) ;
639  }
640 
641 
642  for(j=0;j<*nmpc;j++){
643  b=DenseMtx_entries(mtxX)[j];
644  if(fabs(b)>1.e-10){
645  nodempc[3**mpcfree-1]=ipompc[j];
646  node=(ITG)((idof+8)/8);
647  idir=idof+1-8*(node-1);
648  nodempc[3**mpcfree-3]=node;
649  nodempc[3**mpcfree-2]=idir;
650  coefmpc[*mpcfree-1]=b;
651  ipompc[j]=(*mpcfree)++;
652  if(*mpcfree>*memmpc_){
653  *memmpc_=(ITG)(1.1**memmpc_);
654  RENEW(nodempc,ITG,3**memmpc_);
655  RENEW(coefmpc,double,*memmpc_);
656  }
657  }
658  }
659  }
660  }
661 /*--------------------------------------------------------------------*/
662 /*
663  -----------
664  free memory
665  -----------
666 */
667  FrontMtx_free(frontmtx) ;
668  DenseMtx_free(mtxB) ;
669  DenseMtx_free(mtxX) ;
670  IV_free(newToOldIV) ;
671  IV_free(oldToNewIV) ;
672  InpMtx_free(mtxA) ;
673  ETree_free(frontETree) ;
674  IVL_free(symbfacIVL) ;
675  SubMtxManager_free(mtxmanager) ;
676  Graph_free(graph) ;
677 
678 /* diagonal terms */
679 
680  for(i=0;i<*nmpc;i++){
681  j=ilmpc[i]-1;
682  idof=ikmpc[i];
683  node=(ITG)((idof+7)/8);
684  idir=idof-8*(node-1);
685  nodempc[3**mpcfree-1]=ipompc[j];
686  nodempc[3**mpcfree-3]=node;
687  nodempc[3**mpcfree-2]=idir;
688  coefmpc[*mpcfree-1]=1.;
689  ipompc[j]=(*mpcfree)++;
690  if(*mpcfree>*memmpc_){
691  *memmpc_=(ITG)(1.1**memmpc_);
692  RENEW(nodempc,ITG,3**memmpc_);
693  RENEW(coefmpc,double,*memmpc_);
694  }
695  }
696 
697  SFREE(ipointer);SFREE(indepdof);SFREE(icoef);SFREE(xcoef);
698 
699  fclose(msgFile);
700 
701  }
702 #endif
703 
704 /* determining the effective size of nodempc and coefmpc for
705  the reallocation*/
706 
707  *mpcend=0;
708  *mpcmult=0;
709  *maxlenmpc=0;
710  for(i=0;i<*nmpc;i++){
711  index=ipompc[i];
712  *mpcend=max(*mpcend,index);
713  nterm=1;
714  while(1){
715  index=nodempc[3*index-1];
716  if(index==0){
717  *mpcmult+=nterm*(nterm-1);
718  *maxlenmpc=max(*maxlenmpc,nterm);
719  break;
720  }
721  *mpcend=max(*mpcend,index);
722  nterm++;
723  }
724  }
725 
726  SFREE(jmpc);
727 
728  *nodempcp=nodempc;
729  *coefmpcp=coefmpc;
730 
731 /* for(i=0;i<*nmpc;i++){
732  j=i+1;
733  FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
734  }*/
735 
736  return;
737 }
#define ITGFORMAT
Definition: CalculiX.h:52
#define max(a, b)
Definition: cascade.c:32
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
subroutine stop()
Definition: stop.f:19
void tau(double *ad, double **aup, double *adb, double *aubp, double *sigma, double *b, ITG *icol, ITG **irowp, ITG *neq, ITG *nzs)
#define RENEW(a, b, c)
Definition: CalculiX.h:40
#define SFREE(a)
Definition: CalculiX.h:41
subroutine nident(x, px, n, id)
Definition: nident.f:25
#define ITG
Definition: CalculiX.h:51
#define NNEW(a, b, c)
Definition: CalculiX.h:39
Hosted by OpenAircraft.com, (Michigan UAV, LLC)