00001 #include "pmetrinit.h"
00002 #include "metrinit.h"
00003 #include "mefelinit.h"
00004 #include "trfelinit.h"
00005 #include "pglobalc.h"
00006 #include "pglobalt.h"
00007 #include "pglobal.h"
00008 #include "xfile.h"
00009 #include "globalc.h"
00010 #include "globalt.h"
00011 #include "globalc.h"
00012 #include "mpi.h"
00013 #include <stdio.h>
00014 #include <string.h>
00015 #include <time.h>
00016 #include "stacktrace.h"
00017 #include "globmatc.h"
00018
00019 void pmetr_init (int argc,const char *argv[], stochdriver *stochdm, stochdrivert *stochdt)
00020 {
00021 long i, merr, tmp,n;
00022 char name[1100];
00023 XFILE *inm,*intr,*inc;
00024
00025 set_prgname(argv[0]);
00026
00027 Pcp = new pprobdescc;
00028
00029
00030 Psolm = new psolver (Nproc,Myrank,proc_namec,nameLengthc);
00031
00032
00033 Psolt = new psolver (Nproc,Myrank,proc_namec,nameLengthc);
00034
00035
00036
00037 if (Bar2d==NULL) Bar2d = new barel2d;
00038 if (Beam2d==NULL) Bar2d = new barel2d;
00039 if (Pelq==NULL) Pelq = new planeelemlq;
00040 if (Peqq==NULL) Peqq = new planeelemqq;
00041 if (Asymlq==NULL) Asymlq = new axisymlq;
00042 if (Lhex==NULL) Lhex = new linhex;
00043
00044 if (Ltet==NULL) Ltet = new lintet;
00045 if (Bar3d==NULL) Bar3d = new barel3d;
00046
00047
00048
00049 Cp = new probdescc;
00050 if (process_argsc(argc, argv)){
00051 delete Cp;
00052 abort();
00053 }
00054
00055 Mp = new probdesc;
00056
00057 Tp = new probdesct;
00058
00059 Gtm = new gtopology;
00060
00061 Gtt = new gtopology;
00062
00063 Gtu = new gtopology;
00064
00065 Mt = new mechtop;
00066
00067 Tt = new transtop;
00068
00069 Ct = new couptop;
00070
00071 Mm = new mechmat;
00072
00073 Tm = new transmat;
00074
00075 Cmu = new coupmatu;
00076
00077 Cml = new coupmatl;
00078
00079 Mc = new mechcrsec;
00080
00081 Tc = new transcrsec;
00082
00083 Cb = new coupbclc;
00084
00085 Mb = new mechbclc;
00086
00087 Tb = new transbclc;
00088
00089 Outdm = new outdriverm;
00090
00091 Outdt = new outdrivert;
00092
00093 Outdc = new outdriverc;
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 strcpy (name,argv[1]);
00105 sprintf (&name[strlen (argv[1])],"%d.in",Myrank+1);
00106 inc = xfopen (name,"r");
00107 if (inc==NULL){
00108 par_print_err(Myrank+1, proc_namec, "input file cannot be opened.",__FILE__, __LINE__, __func__);
00109 metr_dealloc();
00110 abort();
00111 }
00112
00113 inc->warning = 1;
00114 inc->ignorecase = 1;
00115 if ((Cp->kwdsw == nokwd) || (Cp->kwdsw == kwd_outd))
00116 inc->kwdmode = ignore;
00117 else
00118 inc->kwdmode = sequent_mode;
00119
00120
00121 xfscanf (inc,"%ld",&Ndom);
00122 Ndom--;
00123
00124
00125 Cp->read (inc);
00126
00127
00128 strcpy (name,Cp->minfile);
00129 sprintf (&name[strlen (Cp->minfile)],"%d.in",Myrank+1);
00130 inm = xfopen(name,"r");
00131 if (inm==NULL){
00132 par_print_err(Myrank+1, proc_namec, "input file for mechanical part cannot be opened.",__FILE__, __LINE__, __func__);
00133 abort ();
00134 }
00135
00136 strcpy (name,Cp->tinfile);
00137 sprintf (&name[strlen (Cp->tinfile)],"%d.in",Myrank+1);
00138 intr = xfopen(name,"r");
00139 if (intr==NULL){
00140 par_print_err(Myrank+1, proc_namec, "input file for thermo-hydro part cannot be opened.",__FILE__, __LINE__, __func__);
00141 abort ();
00142 }
00143 strcpy (name,argv[1]);
00144 sprintf (&name[strlen (argv[1])],"%d.log",Myrank+1);
00145
00146 Outc = fopen (name,"w");
00147 if (Outc==NULL){
00148 par_print_err(Myrank+1, proc_namec, "log file for coupled part cannot be opened.",__FILE__, __LINE__, __func__);
00149 abort ();
00150 }
00151 strcpy (name,Cp->minfile);
00152 sprintf (&name[strlen (Cp->minfile)],"%d.log",Myrank+1);
00153 Out = fopen (name,"w");
00154
00155 if (Out==NULL){
00156 par_print_err(Myrank+1, proc_namec, "log file for mechanical part cannot be opened.",__FILE__, __LINE__, __func__);
00157 abort ();
00158 }
00159 strcpy (name,Cp->tinfile);
00160 sprintf (&name[strlen (Cp->tinfile)],"%d.log",Myrank+1);
00161 Outt = fopen (name,"w");
00162
00163 if (Outt==NULL){
00164 par_print_err(Myrank+1, proc_namec, "log file for thermo-hydro part cannot be opened.",__FILE__, __LINE__, __func__);
00165 abort ();
00166 }
00167
00168 inm->warning = 1;
00169 inm->ignorecase = 1;
00170 if ((Mp->kwdsw == nokwd) || (Mp->kwdsw == kwd_outd))
00171 inm->kwdmode = ignore;
00172 else
00173 inm->kwdmode = sequent_mode;
00174
00175 intr->warning = 1;
00176 intr->ignorecase = 1;
00177 if ((Tp->kwdsw == nokwd) || (Tp->kwdsw == kwd_outd))
00178 intr->kwdmode = ignore;
00179 else
00180 intr->kwdmode = sequent_mode;
00181
00182
00183 xfscanf (inm,"%ld", &tmp);
00184
00185 Mp->read (inm);
00186
00187
00188 xfscanf (intr,"%ld", &tmp);
00189
00190 Tp->read (intr);
00191
00192
00193 Psolm->read (inm,Gtm,Ndom,Mp->tstorsm,Mespr);
00194
00195
00196 Psolt->read (intr,Gtt,Ndom,Tp->tstorkm,Mesprt);
00197
00198 if ((Cp->tprob==growing_par_coupl_mech_trans) || (Cp->tprob==par_coupl_mech_trans)){
00199
00200 Cp->time=Cp->timecon.starttime ();
00201 Tp->time=Tp->timecont.starttime ();
00202 Mp->time=Mp->timecon.starttime ();
00203 }
00204
00205 if (Mp->kwdsw == kwd_full)
00206 inm->kwdmode = sequent_mode;
00207 else
00208 inm->kwdmode = ignore;
00209
00210 if (Tp->kwdsw == kwd_full)
00211 intr->kwdmode = sequent_mode;
00212 else
00213 intr->kwdmode = ignore;
00214
00215 if (Cp->kwdsw == kwd_full)
00216 inc->kwdmode = sequent_mode;
00217 else
00218 inc->kwdmode = ignore;
00219
00220
00221 Psolm->procdomcorr ();
00222 Psolt->procdomcorr ();
00223
00224
00225
00226 Mt->read (inm);
00227 Tt->read (intr);
00228 Ct->read (inc);
00229
00230
00231 if (Mp->tprob == growing_mech_structure){
00232 if (Mespr==1){
00233 merr=Mt->mesh_check();
00234 if (merr > 0){
00235 fprintf (stdout,"MEFEL Mesh error = %ld\n\n\n",merr);
00236 exit(0);
00237 }
00238 }
00239 }
00240
00241
00242
00243 if (Tp->tprob == growing_np_problem || Tp->tprob == growing_np_problem_nonlin){
00244 if (Mesprt==1){
00245 merr=Tt->mesh_check();
00246 if (merr > 0){
00247 fprintf (stdout,"TRFEL Mesh error = %ld\n\n\n",merr);
00248 exit(0);
00249 }
00250 }
00251 }
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266 if (Cp->tprob==fully_coupled_mech_trans){
00267
00268 Gtu->comptop (Gtm,Gtt);
00269 }
00270
00271
00272 Psolm->initiate (Gtm,argc,argv);
00273
00274 Psolt->initiate (Gtt,argc,argv);
00275
00276
00277
00278 Psolm->read_ltg (inm,Gtm,Out);
00279 Psolt->read_ltg (intr,Gtt,Outt);
00280
00281
00282
00283 Ndofm=Psolm->ordering (Gtm,Out,proc_namec);
00284 if (Mesprc==1)
00285 fprintf (stdout,"\n number of degrees of freedom in mechanics %ld",Ndofm);
00286
00287 Ndoft=Psolt->ordering (Gtt,Outt,proc_namec);
00288 if (Mesprc==1)
00289 fprintf (stdout,"\n number of degrees of freedom in transports %ld",Ndoft);
00290 if (Cp->tprob==fully_coupled_mech_trans){
00291
00292 Ndofc=Gtu->gencodnum ();
00293
00294 Gtu->cndistr (Gtm,Gtt);
00295 Ndofm=Ndofc;
00296 Ndoft=Ndofc;
00297 if (Mesprc==1)
00298 fprintf (stdout,"\n total number of degrees of freedom %ld",Ndofc);
00299 }
00300
00301
00302 if (Tp->tprob == growing_np_problem || Tp->tprob == growing_np_problem_nonlin){
00303 Gtt->read_gf (intr);
00304 }
00305
00306
00307 Mm->read (inm);
00308
00309
00310 Mm->init_ip_1 ();
00311
00312
00313 Tm->read (intr);
00314
00315
00316 Tm->intpointalloc ();
00317 Tm->intpointinit ();
00318
00319
00320 if (Cp->tprob==fully_coupled_mech_trans){
00321 Cmu->read (inc);
00322 Cml->read (inc);
00323 }
00324
00325
00326 Mc->read (inm);
00327 Tc->read (intr);
00328
00329
00330
00331 Mb->read (inm);
00332 Tb->read (intr);
00333
00334
00335 if (Mp->tprob == growing_mech_structure){
00336 Gtm->read_gf (inm);
00337 Gtm->alloc_growstr();
00338 Mt->alloc_growstr();
00339 }
00340
00341
00342
00343
00344 Tb->elemsource ();
00345
00346
00347
00348 n = Mm->search_reqnmq(Mm->nmqo);
00349 if (n){
00350 for (i=0; i<n; i++)
00351 Mm->nmqid[Mm->nmqo[i]-1] = i;
00352
00353 Mm->alloc_nonmechq(n);
00354 }
00355
00356 n = Tm->search_reqntq(Tm->ntqo);
00357 if (n){
00358 for (i=0; i<n; i++)
00359 Tm->ntqid[Tm->ntqo[i]-1] = i;
00360
00361 Tm->alloc_nontransq(n);
00362 }
00363
00364
00365 TM_ip = new long[Tm->tnip];
00366 MT_ip = new long[Mm->tnip];
00367 mefel_trfel_ip_mapping(TM_ip, MT_ip);
00368
00369 if (Cp->tprob==fully_coupled_mech_trans){
00370 if (Mp->pore_press)
00371 Mp->pore_press = 3;
00372
00373
00374 Lsrsc = new lhsrhsc;
00375 Lsrsc->alloc ();
00376
00377
00378 Lsrs = new lhsrhs;
00379 Lsrs->lhsi=Lsrsc->lhsi;
00380 Lsrs->initcond (inm);
00381
00382 Lsrst = new lhsrhst;
00383 Lsrst->lhsi=Lsrsc->lhsi;
00384 Lsrst->initcond (intr);
00385 }
00386
00387 if ((Cp->tprob==par_coupl_mech_trans) || (Cp->tprob==growing_par_coupl_mech_trans)){
00388 if (Mp->pore_press)
00389 Mp->pore_press = 2;
00390
00391
00392 Lsrs = new lhsrhs;
00393 Lsrs->alloc ();
00394 Lsrs->initcond (inc);
00395
00396
00397 Lsrst = new lhsrhst;
00398 Lsrst->alloc ();
00399 Lsrst->initcond (intr);
00400
00401 }
00402
00403 Mb->alloc_sumcomp ();
00404 if (Mm->csol)
00405 Gtm->comp_domain_sizes ();
00406
00407
00408
00409 if ((Cp->kwdsw == nokwd) || (Cp->kwdsw == kwd_pd))
00410 inc->kwdmode = ignore;
00411 else
00412 inc->kwdmode = sequent_mode;
00413
00414 Outdc->read(inc);
00415
00416
00417 if ((Mp->kwdsw == nokwd) || (Mp->kwdsw == kwd_pd))
00418 inm->kwdmode = ignore;
00419 else
00420 inm->kwdmode = sequent_mode;
00421
00422 Outdm->read(inm);
00423
00424
00425 if ((Tp->kwdsw == nokwd) || (Tp->kwdsw == kwd_pd))
00426 intr->kwdmode = ignore;
00427 else
00428 intr->kwdmode = sequent_mode;
00429
00430 Outdt->read(intr);
00431
00432
00433
00434 if (Mp->hdbcont.restore_stat())
00435 {
00436 strcpy (name,Mp->hdbcont.hdbnamer);
00437 sprintf (&name[strlen (Mp->hdbcont.hdbnamer)],"%d.out",Myrank+1);
00438 strcpy (Mp->hdbcont.hdbnamer, name);
00439 }
00440 if (Mp->hdbcont.save_stat())
00441 {
00442 strcpy (name,Mp->hdbcont.hdbnames);
00443 sprintf (&name[strlen (Mp->hdbcont.hdbnames)],"%d.out",Myrank+1);
00444 strcpy (Mp->hdbcont.hdbnames, name);
00445 }
00446
00447
00448 strcpy (name,Outdm->outfn);
00449 sprintf (&name[strlen (Outdm->outfn)],"%d.out",Myrank+1);
00450 strcpy (Outdm->outfn, name);
00451 if (Outdm->gf != grfmt_no)
00452 {
00453 strcpy (name,Outdm->outgrfn);
00454 sprintf (&name[strlen (Outdm->outgrfn)],"%d",Myrank+1);
00455 strcpy (Outdm->outgrfn, name);
00456 }
00457 if (Outdm->ndiag > 0)
00458 {
00459 strcpy (name,Outdm->outdiagfn);
00460 sprintf (&name[strlen (Outdm->outdiagfn)],"%d.dat",Myrank+1);
00461 strcpy (Outdm->outdiagfn, name);
00462 }
00463
00464
00465 if (Tp->hdbcont.restore_stat())
00466 {
00467 strcpy (name,Tp->hdbcont.hdbnamer);
00468 sprintf (&name[strlen (Tp->hdbcont.hdbnamer)],"%d.out",Myrank+1);
00469 strcpy (Tp->hdbcont.hdbnamer, name);
00470 }
00471 if (Tp->hdbcont.save_stat())
00472 {
00473 strcpy (name,Tp->hdbcont.hdbnames);
00474 sprintf (&name[strlen (Tp->hdbcont.hdbnames)],"%d.out",Myrank+1);
00475 strcpy (Tp->hdbcont.hdbnames, name);
00476 }
00477
00478
00479 strcpy (name,Outdt->outfn);
00480 sprintf (&name[strlen (Outdt->outfn)],"%d.out",Myrank+1);
00481 strcpy (Outdt->outfn, name);
00482 if (Outdt->gf != grftt_no)
00483 {
00484 strcpy (name,Outdt->outgrfn);
00485 sprintf (&name[strlen (Outdt->outgrfn)],"%d",Myrank+1);
00486 strcpy (Outdt->outgrfn, name);
00487 }
00488 if (Outdt->ndiag > 0)
00489 {
00490 strcpy (name,Outdt->outdiagfn);
00491 sprintf (&name[strlen (Outdt->outdiagfn)],"%d.dat",Myrank+1);
00492 strcpy (Outdt->outdiagfn, name);
00493 }
00494
00495
00496
00497
00498 Mm->init_ip_2 ();
00499
00500
00501
00502 Tt->alloc_nodes ();
00503
00504 Mt->alloc_nodes_arrays ();
00505
00506
00507 Mt->give_maxncompstr(Mm->max_ncompstrn, Mm->max_ncompstre);
00508
00509 Mt->give_maxncompo(Mm->max_nncompo, Mm->max_encompo);
00510
00511
00512
00513
00514 Mt->elemprescdisp ();
00515
00516 if (Mp->reactcomp==1){
00517 Mt->comreac ();
00518 }
00519
00520
00521 if (Mb->ico > 0)
00522 Mb->inicipval();
00523
00524 if (Mp->matmodel == nonlocal){
00525
00526 Gtm->adjacelem (Out);
00527
00528
00529 if (restore_adjacip() == 0)
00530 adjacip ();
00531 save_adjacip();
00532 ipvolume ();
00533 }
00534
00535
00536
00537 Mm->alloceigstrains ();
00538 Mm->alloctempstrains ();
00539 Mm->alloceigstresses ();
00540
00541
00542
00543 xfclose (inc);
00544 xfclose (inm);
00545 xfclose (intr);
00546
00547 Pcp->shift_indices ();
00548 }
00549
00550 void metr_dealloc()
00551 {
00552 delete Pcp;
00553 delete Psolm;
00554 delete Psolt;
00555
00556 delete Pelq;
00557 delete Asymlq;
00558 delete Bar2d;
00559
00560
00561 delete Cp;
00562 delete Mp;
00563 delete Tp;
00564 delete Gtm;
00565 delete Gtt;
00566 delete Gtu;
00567 delete Mt;
00568 delete Tt;
00569 delete Ct;
00570 delete Mm;
00571 delete Tm;
00572 delete Cmu;
00573 delete Cml;
00574 delete Mc;
00575 delete Tc;
00576 delete Cb;
00577 delete Mb;
00578 delete Tb;
00579 delete Outdm;
00580 delete Outdt;
00581 delete Outdc;
00582 }