00001 #include "lfssolver.h" 00002 #include "global.h" 00003 #include "globmat.h" 00004 #include "loadcase.h" 00005 #include "gmatrix.h" 00006 #include "mechprint.h" 00007 #include "flsubdom.h" 00008 #include "elemswitch.h" 00009 00010 00011 /** 00012 function solves linear problem with floating subdomains 00013 it may be used for pull-out tests 00014 00015 JK, 11. 10. 2007, modified 3. 8. 2012 00016 */ 00017 void solve_linear_floating_subdomains () 00018 { 00019 long i; 00020 double *lhs,*rhs; 00021 00022 // stiffness matrix assembling 00023 stiffness_matrix (0); 00024 00025 // array is allocated in Lsrs 00026 lhs = Lsrs->give_lhs (0); 00027 rhs = Lsrs->give_rhs (0); 00028 00029 // loop over load cases 00030 for (i=0;i<Lsrs->nlc;i++){ 00031 // right hand side of system 00032 mefel_right_hand_side (i,rhs); 00033 } 00034 00035 // solution of equation system 00036 for (i=0;i<Lsrs->nlc;i++){ 00037 Mp->ssle->solve_system (Gtm,Smat,Lsrs->give_lhs(i),Lsrs->give_rhs(i),Out); 00038 } 00039 00040 print_init(-1, "wt"); 00041 for (i=0;i<Lsrs->nlc;i++){ 00042 // computes and prints required quantities 00043 //compute_ipstrains (i); 00044 //compute_ipstresses (i); 00045 compute_req_val (i); 00046 print_step(i, 0, 0.0, NULL); 00047 } 00048 print_close(); 00049 00050 } 00051 00052 00053 00054 00055 00056 /** 00057 function solves linear problem with floating subdomains 00058 it may be used for pull-out tests 00059 00060 JK, 11.10.2007 00061 */ 00062 void solve_linear_floating_subdomains2 () 00063 { 00064 00065 00066 //Out=fopen("dupnode","w"); 00067 00068 Gtm->edges (Out); 00069 Gtm->edgenode_sorting (); 00070 Gtm->prev_next_edges (); 00071 Gtm->edge_series (Out); 00072 Gtm->edge_dirvect (); 00073 Gtm->edge_normvect (); 00074 00075 Gtm->create_ltg (Out); 00076 00077 long i; 00078 00079 /* 00080 for (i=0;i<Gtm->nged;i++){ 00081 fprintf (Out,"\n %6ld %6ld %6ld %6ld",Gtm->gedges[i].nlist[0],Gtm->gedges[i].nlist[1],Gtm->gedges[i].nlist[2],Gtm->gedges[i].nlist[3]); 00082 } 00083 00084 fprintf (Out,"\n"); 00085 for (i=0;i<Gtm->nged;i++){ 00086 fprintf (Out,"\n %6ld %6ld",Gtm->gedges[i].nlist[0]+1,Gtm->gedges[i].nlist[1]+1); 00087 } 00088 i=Gtm->nged-1; 00089 fprintf (Out,"\n %6ld %6ld",Gtm->gedges[i].nlist[2]+1,Gtm->gedges[i].nlist[3]+1); 00090 00091 00092 00093 00094 00095 00096 00097 00098 00099 //fprintf (Out,"\n\n kontrola prev a next \n"); 00100 //for (i=0;i<Gtm->nged;i++){ 00101 //fprintf (Out,"\n hrana %6ld prev %6ld next %6ld",i,Gtm->gedges[i].prev,Gtm->gedges[i].next); 00102 //} 00103 //fprintf (Out,"\n\n"); 00104 00105 00106 fprintf (Out,"\n\n kontrola souslednosti\n"); 00107 long *av; 00108 av = new long [Gtm->nged]; 00109 00110 for (i=0;i<Gtm->nged;i++){ 00111 av[i]=-1; 00112 } 00113 for (i=0;i<Gtm->nged;i++){ 00114 if (Gtm->gedges[i].prev==-1){ 00115 a=i; 00116 n=Gtm->gedges[i].next; 00117 00118 fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00119 fprintf (Out,"\n %6ld %6ld",Gtm->gedges[a].nlist[0],Gtm->gedges[a].nlist[2]); 00120 //fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[0]); 00121 00122 av[i]=1; 00123 break; 00124 } 00125 } 00126 long ui=Gtm->nged; 00127 00128 for (i=1;i<ui;i++){ 00129 a=n; 00130 n=Gtm->gedges[n].next; 00131 av[a]=1; 00132 fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[0]); 00133 //fprintf (Out,"\n %6ld %6ld",Gtm->gedges[a].nlist[0],Gtm->gedges[a].nlist[2]); 00134 fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00135 if (n<0){ 00136 for (j=0;j<Gtm->nged;j++){ 00137 if (Gtm->gedges[j].prev==-1 && av[j]==-1){ 00138 a=j; 00139 n=Gtm->gedges[j].next; 00140 av[j]=1; 00141 ui--; 00142 fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[0]); 00143 /fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00144 break; 00145 } 00146 } 00147 } 00148 } 00149 fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[2]); 00150 00151 00152 fprintf (Out,"\n"); 00153 00154 00155 for (i=0;i<Gtm->nged;i++){ 00156 av[i]=-1; 00157 } 00158 for (i=0;i<Gtm->nged;i++){ 00159 if (Gtm->gedges[i].prev==-1){ 00160 a=i; 00161 n=Gtm->gedges[i].next; 00162 00163 //fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00164 //fprintf (Out,"\n %6ld %6ld ",Gtm->gedges[a].nlist[1]-205,Gtm->gedges[a].nlist[3]-205); 00165 fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[1]-205); 00166 00167 av[i]=1; 00168 break; 00169 } 00170 } 00171 ui=Gtm->nged; 00172 00173 for (i=1;i<ui;i++){ 00174 a=n; 00175 n=Gtm->gedges[n].next; 00176 av[a]=1; 00177 fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[1]-205); 00178 //fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00179 if (n<0){ 00180 for (j=0;j<Gtm->nged;j++){ 00181 if (Gtm->gedges[j].prev==-1 && av[j]==-1){ 00182 a=j; 00183 n=Gtm->gedges[j].next; 00184 av[j]=1; 00185 ui--; 00186 fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[1]-205); 00187 //fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00188 break; 00189 } 00190 } 00191 } 00192 } 00193 fprintf (Out,"\n %6ld",Gtm->gedges[a].nlist[3]-205); 00194 00195 00196 00197 00198 00199 fprintf (Out,"\n"); 00200 00201 00202 for (i=0;i<Gtm->nged;i++){ 00203 av[i]=-1; 00204 } 00205 for (i=0;i<Gtm->nged;i++){ 00206 if (Gtm->gedges[i].prev==-1){ 00207 a=i; 00208 n=Gtm->gedges[i].next; 00209 00210 fprintf (Out,"\n %lf %lf %lf %lf ",Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1]); 00211 //fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00212 //fprintf (Out,"\n %6ld %6ld ",Gtm->gedges[a].nlist[1]-205,Gtm->gedges[a].nlist[3]-205); 00213 00214 av[i]=1; 00215 break; 00216 } 00217 } 00218 ui=Gtm->nged; 00219 00220 for (i=1;i<ui;i++){ 00221 a=n; 00222 n=Gtm->gedges[n].next; 00223 av[a]=1; 00224 //fprintf (Out,"\n %6ld %6ld",Gtm->gedges[a].nlist[1]-205,Gtm->gedges[a].nlist[3]-205); 00225 //fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00226 fprintf (Out,"\n %lf %lf %lf %lf",Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1]); 00227 if (n<0){ 00228 for (j=0;j<Gtm->nged;j++){ 00229 if (Gtm->gedges[j].prev==-1 && av[j]==-1){ 00230 a=j; 00231 n=Gtm->gedges[j].next; 00232 av[j]=1; 00233 ui--; 00234 fprintf (Out,"\n %6ld %6ld",Gtm->gedges[a].nlist[1]-205,Gtm->gedges[a].nlist[3]-205); 00235 //fprintf (Out,"\n %6ld -> %6ld %lf %lf %lf %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00236 fprintf (Out,"\n %lf %lf %lf %lf",Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1]); 00237 break; 00238 } 00239 } 00240 } 00241 } 00242 00243 00244 00245 00246 00247 00248 00249 00250 00251 00252 00253 00254 00255 00256 00257 //a=0; n=Gtm->gedges[a].next; 00258 //fprintf (Out,"\n\n kontrola serie \n"); 00259 //for (i=0;i<Gtm->nged;i++){ 00260 //fprintf (Out,"\n%ld -> %ld",a,n); 00261 //a=n; 00262 //n=Gtm->gedges[a].next; 00263 //} 00264 00265 00266 00267 fprintf (Out,"\n\n kontrola souslednosti\n"); 00268 long *av; 00269 av = new long [Gtm->nged]; 00270 00271 for (i=0;i<Gtm->nged;i++){ 00272 av[i]=-1; 00273 } 00274 for (i=0;i<Gtm->nged;i++){ 00275 if (Gtm->gedges[i].prev==-1){ 00276 a=i; 00277 n=Gtm->gedges[i].next; 00278 00279 fprintf (Out,"\n %6ld %6ld %6ld %6ld",Gtm->gedges[a].nlist[0],Gtm->gedges[a].nlist[2],Gtm->gedges[a].nlist[1],Gtm->gedges[a].nlist[3]); 00280 fprintf (Out," %lf %lf %lf %lf %lf %lf",Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00281 00282 av[i]=1; 00283 break; 00284 } 00285 } 00286 long ui=Gtm->nged; 00287 00288 for (i=1;i<ui;i++){ 00289 a=n; 00290 n=Gtm->gedges[n].next; 00291 av[a]=1; 00292 fprintf (Out,"\n %6ld %6ld %6ld %6ld",Gtm->gedges[a].nlist[0],Gtm->gedges[a].nlist[2],Gtm->gedges[a].nlist[1],Gtm->gedges[a].nlist[3]); 00293 fprintf (Out," %lf %lf %lf %lf %lf %lf",Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00294 if (n<0){ 00295 for (j=0;j<Gtm->nged;j++){ 00296 if (Gtm->gedges[j].prev==-1 && av[j]==-1){ 00297 a=j; 00298 n=Gtm->gedges[j].next; 00299 av[j]=1; 00300 ui--; 00301 fprintf (Out,"\n %6ld %6ld %6ld %6ld",Gtm->gedges[a].nlist[0],Gtm->gedges[a].nlist[2],Gtm->gedges[a].nlist[1],Gtm->gedges[a].nlist[3]); 00302 fprintf (Out," %lf %lf %lf %lf %lf %lf",Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]); 00303 break; 00304 } 00305 } 00306 } 00307 } 00308 00309 00310 fprintf (Out,"\n\n\n"); 00311 00312 00313 00314 00315 00316 00317 for (i=0;i<Gtm->nged;i++){ 00318 fprintf (Out,"\n\n hrana cislo %ld",i); 00319 Gtm->gedges[i].print (Out); 00320 } 00321 00322 00323 fprintf (Out,"\n\n\n kontrola serii \n"); 00324 for (i=0;i<Gtm->nged;i++){ 00325 fprintf (Out,"\n edge %6ld je v serii %6ld",i,Gtm->edgeser[i]); 00326 } 00327 00328 //fclose (Out); 00329 00330 */ 00331 00332 00333 00334 double *lhs,*rhs; 00335 00336 // stiffness matrix assembling 00337 stiffness_matrix (0); 00338 00339 // array is allocated in Lsrs 00340 lhs = Lsrs->give_lhs (0); 00341 rhs = Lsrs->give_rhs (0); 00342 00343 // loop over load cases 00344 for (i=0;i<Lsrs->nlc;i++){ 00345 // right hand side of system 00346 mefel_right_hand_side (i,rhs); 00347 } 00348 00349 00350 00351 /* 00352 Smat->printmat(Out); 00353 FILE *p; 00354 p = fopen ("vektor.txt","w"); 00355 fprintf (p,"%ld\n",Ndofm); 00356 for (i=0;i<Ndofm;i++){ 00357 fprintf (p,"%20.16le\n",rhs[i]); 00358 } 00359 fclose (p); 00360 */ 00361 00362 00363 00364 00365 //Smat->printdiag(Out); 00366 00367 // solution of equation system 00368 for (i=0;i<Lsrs->nlc;i++){ 00369 Mp->ssle->solve_system (Gtm,Smat,Lsrs->give_lhs(i),Lsrs->give_rhs(i),Out); 00370 //Smat->solve_system (Lsrs->give_lhs(i),Lsrs->give_rhs(i)); 00371 } 00372 00373 print_init(-1, "wt"); 00374 for (i=0;i<Lsrs->nlc;i++){ 00375 // computes and prints required quantities 00376 //compute_ipstrains (i); 00377 //compute_ipstresses (i); 00378 compute_req_val (i); 00379 print_step(i, 0, 0.0, NULL); 00380 } 00381 print_close(); 00382 00383 00384 00385 00386 00387 00388 } 00389 00390 /** 00391 function solves linear problem with floating subdomains 00392 it may be used for pull-out tests 00393 00394 JK, 5.11.2005 00395 */ 00396 void solve_linear_floating_subdomains_old () 00397 { 00398 /* 00399 long i,lcid; 00400 double *lhs,*rhs; 00401 //flsubdom *fsd; 00402 00403 lcid=0; 00404 00405 Fsd = new flsubdom; 00406 00407 lhs=Lsrs->give_lhs (lcid); 00408 rhs=Lsrs->give_rhs (lcid); 00409 00410 mefel_right_hand_side (lcid,rhs); 00411 00412 stiffness_matrix (lcid); 00413 00414 Fsd->initialization (Ndofm,Mp->ense,rhs); 00415 00416 Fsd->solve_lin_alg_system (lhs,rhs); 00417 00418 print_init(-1, "wt"); 00419 for (i=0;i<Lsrs->nlc;i++){ 00420 compute_ipstrains (i); 00421 compute_ipstresses (i); 00422 compute_req_val (i); 00423 print_step(i, 0, 0.0, NULL); 00424 } 00425 print_close(); 00426 */ 00427 00428 00429 00430 00431 00432 long i; 00433 double *lhs,*rhs; 00434 00435 // stiffness matrix assembling 00436 stiffness_matrix (0); 00437 00438 // array is allocated in Lsrs 00439 lhs = Lsrs->give_lhs (0); 00440 rhs = Lsrs->give_rhs (0); 00441 00442 // loop over load cases 00443 for (i=0;i<Lsrs->nlc;i++){ 00444 // right hand side of system 00445 mefel_right_hand_side (i,rhs); 00446 } 00447 00448 00449 00450 00451 00452 /* 00453 Smat->printmat(Out); 00454 FILE *p; 00455 p = fopen ("vektor.txt","w"); 00456 fprintf (p,"%ld\n",Ndofm); 00457 for (i=0;i<Ndofm;i++){ 00458 fprintf (p,"%20.16le\n",rhs[i]); 00459 } 00460 fclose (p); 00461 */ 00462 00463 00464 00465 00466 //Smat->printdiag(Out); 00467 00468 double *condmat,*condvect; 00469 long nrdof=Ndofm-Gtm->nbdof; 00470 condmat = new double [nrdof*nrdof]; 00471 condvect = new double [nrdof]; 00472 00473 // solution of equation system 00474 for (i=0;i<Lsrs->nlc;i++){ 00475 //Mp->ssle.solve_system (Gtm,Smat,Lsrs->give_lhs(i),Lsrs->give_rhs(i),Out); 00476 //Smat->solve_system (Lsrs->give_lhs(i),Lsrs->give_rhs(i)); 00477 Smat->condense (Gtm,condmat,condvect,Lsrs->give_lhs(i),Lsrs->give_rhs(i),nrdof,1,Out); 00478 } 00479 00480 long j; 00481 FILE *out; 00482 out = fopen ("matice.txt","w"); 00483 for (i=0;i<nrdof;i++){ 00484 for (j=0;j<nrdof;j++){ 00485 fprintf (out,"%le ",condmat[i*nrdof+j]); 00486 } 00487 fprintf (out,"\n"); 00488 } 00489 fclose (out); 00490 out = fopen ("vektor.txt","w"); 00491 for (i=Gtm->nbdof;i<Ndofm;i++){ 00492 fprintf (out,"%le\n",rhs[i]); 00493 } 00494 fclose (out); 00495 00496 00497 out = fopen ("nezname.txt","w"); 00498 for (i=0;i<Mt->nn;i++){ 00499 fprintf (out,"%6ld %6ld\n",Gtm->gnodes[i].cn[0],Gtm->gnodes[i].cn[1]); 00500 } 00501 fclose (out); 00502 00503 00504 00505 print_init(-1, "wt"); 00506 for (i=0;i<Lsrs->nlc;i++){ 00507 // computes and prints required quantities 00508 //compute_ipstrains (i); 00509 //compute_ipstresses (i); 00510 compute_req_val (i); 00511 print_step(i, 0, 0.0, NULL); 00512 } 00513 print_close(); 00514 00515 00516 } 00517