00001 #include "solvert.h"
00002 #include "spsolvert.h"
00003 #include "nspsolvert.h"
00004 #include "npsolvert.h"
00005 #include "nnpsolvert.h"
00006 #include "cpnpsolvert.h"
00007 #include "cpnnpsolvert.h"
00008 #include "dnpsolvert.h"
00009 #include "dnnpsolvert.h"
00010 #include "homogtrans.h"
00011 #include "globalt.h"
00012 #include "globmatt.h"
00013 #include "libtrace.h"
00014 #include "trfelinit.h"
00015
00016
00017 void solve_trfel_problem (stochdrivert *stochd)
00018 {
00019
00020 if (Tp->stochasticcalc==0 && Tp->adaptivityflag==0)
00021 solve_trfel_deterministic_problem ();
00022 if (Tp->stochasticcalc && Tp->adaptivityflag==0)
00023 solve_trfel_stochastic_problem (stochd);
00024 if (Tp->stochasticcalc==0 && Tp->adaptivityflag)
00025 solve_trfel_adaptivity_problem ();
00026 if (Tp->stochasticcalc && Tp->adaptivityflag)
00027 print_err("stochastic and adaptivity mix not supported",__FILE__,__LINE__,__func__);
00028 }
00029
00030 void solve_trfel_deterministic_problem ()
00031 {
00032 switch (Tp->tprob){
00033 case stationary_problem:{
00034 if(Tp->homogt == 1)
00035 transport_homogenization();
00036 else{
00037 solve_stationary_problem ();
00038
00039 }
00040 break;
00041 }
00042 case nonlinear_stationary_problem:{
00043 if(Tp->homogt == 1)
00044 transport_homogenization();
00045 else{
00046
00047
00048 solve_nonlinear_stationary_problem_old ();
00049 }
00050 break;
00051 }
00052 case nonstationary_problem:{
00053 solve_nonstationary_problem ();
00054 break;
00055 }
00056 case nonlinear_nonstationary_problem:{
00057
00058 solve_nonstationary_problem ();
00059 break;
00060 }
00061 case growing_np_problem:{
00062
00063 solve_nonstationary_growing_vform ();
00064 break;
00065 }
00066 case growing_np_problem_nonlin:{
00067 solve_nonstationary_growing_problem_nonlin ();
00068 break;
00069 }
00070 case discont_nonstat_problem:{
00071 solve_discont_nonstationary_problem ();
00072 break;
00073 }
00074 case discont_nonlin_nonstat_problem:{
00075 solve_discont_nonlin_nonstationary_problem ();
00076 break;
00077 }
00078 default:{
00079 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00080 }
00081 }
00082 }
00083
00084 void solve_trfel_stochastic_problem (stochdrivert *stochd)
00085 {
00086 print_err("stochastic transport problems are not implemented at this time",__FILE__,__LINE__,__func__);
00087 }
00088
00089
00090
00091
00092
00093
00094
00095
00096 void print_VTK_nodedata (char *filename, gtopology *gt, const double *data)
00097 {
00098 long i;
00099
00100 FILE *out = fopen (filename,"w"); if (out==NULL) { print_err("test file has not been opened", __FILE__, __LINE__, __func__); exit (1); }
00101
00102
00103 fprintf (out, "# vtk DataFile Version 3.0\n");
00104 fprintf (out, "testing VTK file written by Termitovo fci in SIFEL\n");
00105 fprintf (out, "ASCII\n");
00106 fprintf (out, "DATASET POLYDATA\n");
00107 fprintf (out, "POINTS %ld float\n", gt->nn);
00108
00109 for (i=0; i<gt->nn; i++)
00110 fprintf (out, "%13.6f %13.6f %13.6f\n", gt->gnodes[i].x, gt->gnodes[i].y, gt->gnodes[i].z);
00111
00112 fprintf (out, "POLYGONS %ld %ld\n", gt->ne, 4*gt->ne);
00113
00114 for (i=0; i<gt->ne; i++)
00115 fprintf (out, "3 %ld %ld %ld\n", gt->gelements[i].nodes[0], gt->gelements[i].nodes[1], gt->gelements[i].nodes[2]);
00116
00117 fprintf (out, "POINT_DATA %ld\n", gt->nn);
00118 fprintf (out, "SCALARS displacement float 2\n");
00119 fprintf (out, "LOOKUP_TABLE default\n");
00120
00121 for (i=0; i<gt->nn; i++)
00122 fprintf (out, "%13.6f %13.6f\n", data[i], data[i+gt->nn]);
00123
00124 fclose (out);
00125 }
00126
00127 void delete_glob (void)
00128 {
00129 delete Kmat; Kmat = NULL;
00130 delete Cmat; Cmat = NULL;
00131 delete Jmat; Jmat = NULL;
00132 delete Bmat; Bmat = NULL;
00133 delete Lbt; Lbt = NULL;
00134 delete Lbat; Lbat = NULL;
00135 delete Qbt; Qbt = NULL;
00136 delete Qbat; Qbat = NULL;
00137 delete Ltt; Ltt = NULL;
00138 delete Ltat; Ltat = NULL;
00139 delete Lqt; Lqt = NULL;
00140 delete Qqt; Qqt = NULL;
00141 delete Qqat; Qqat = NULL;
00142 delete Lqat; Lqat = NULL;
00143 delete Ltett; Ltett = NULL;
00144 delete Lht; Lht = NULL;
00145 delete Qht; Qht = NULL;
00146 delete G2d; G2d = NULL;
00147 delete Lsrst; Lsrst = NULL;
00148 delete Adat; Adat = NULL;
00149 delete Outdt; Outdt = NULL;
00150 delete Tb; Tb = NULL;
00151 delete Tc; Tc = NULL;
00152 delete Tm; Tm = NULL;
00153 delete Tt; Tt = NULL;
00154 delete Gtt; Gtt = NULL;
00155 delete Tp; Tp = NULL;
00156
00157 Ndoft = Mesprt = 0;
00158 fclose (Outt);
00159 Outt = Outt1 = Outt2 = NULL;
00160 }
00161
00162 void newmeshgen (long i)
00163 {
00164 int ret;
00165 char procname[1023];
00166 char genbin[255]; sprintf (genbin, "/home/dr/Bin/T3d");
00167 char prepbin[255]; sprintf (prepbin, "./transprep");
00168 double d = 0.6;
00169
00170 const char *filename = Adat->give_filename();
00171 const char *ni = Adat->give_ni();
00172 int w = Adat->give_niwidth();
00173
00174
00175 if (Mesprt) fprintf (stdout," ADAPTIVITY: T3D generation\n");
00176 sprintf (procname, "%s -i %s.T3d.in -o %s.T3d.out -m %s.T3d.bgm%s -d %g -p 264 -r 1 -k %d", genbin, filename, filename, filename, ni, d, 1);
00177 fprintf (stdout,"\n %s\n",procname);
00178 if (Mesprt<2) strcat(procname, " >> adaptivity.log 2>&1");
00179 if (Mesprt>1) fprintf (stdout, "%s\n", procname);
00180 ret = system (procname); if (ret != 0) { print_err("Mesh generation failed", __FILE__, __LINE__, __func__); exit (1); }
00181
00182
00183 if (Mesprt) fprintf (stdout," ADAPTIVITY: sifel.in generation\n");
00184 sprintf (procname, "%s %s.sifel.pr %s.sifel.in.%0*ld", prepbin, filename, filename, w, i);
00185 if (Mesprt<2) strcat(procname, " >> adaptivity.log 2>&1");
00186 if (Mesprt>1) fprintf (stdout, "%s\n", procname);
00187 ret = system (procname); if (ret != 0) { print_err("TRFEL PREP failed", __FILE__, __LINE__, __func__); exit (1); }
00188
00189
00190
00191 sprintf (procname,"mv %s.T3d.out %s.T3d.out.%0*ld", filename, filename, w, i); system (procname);
00192
00193 }
00194
00195 void solve_trfel_adaptivity_problem ()
00196 {
00197
00198 fprintf (stdout,"\n\n ADAPTIVITY: start with mesh 0\n");
00199 solve_trfel_deterministic_problem ();
00200
00201 switch (Tp->tprob) {
00202
00203 case stationary_problem:
00204 if (Adat->answer) fprintf (stdout,"\n ADAPTIVITY: new mesh required\n");
00205 else fprintf (stdout,"\n ADAPTIVITY: new mesh NOT required\n");
00206 break;
00207
00208 case nonstationary_problem: {
00209 gtopology *Gtt_old;
00210 adaptivityt *Adat_old;
00211
00212 int argc = 2;
00213 const char *argv[2];
00214 argv[0] = "trfel";
00215
00216 char infile[255];
00217 argv[1] = infile;
00218
00219 int i = 1;
00220 while (Adat->answer) {
00221 fprintf (stdout,"\n ADAPTIVITY: switch form mesh %d to mesh %d\n", i-1, i);
00222
00223
00224 newmeshgen (i);
00225 sprintf (infile, "%s.sifel.in.%0*d", Adat->give_filename(), Adat->give_niwidth(), i);
00226 if (Gtt->nadjelnod == NULL)
00227 Gtt->adjacelem (Outt);
00228
00229
00230
00231 Adat->statedata_backup();
00232
00233 Gtt_old = Gtt; Gtt = NULL;
00234 Adat_old = Adat; Adat = NULL;
00235
00236
00237
00238 delete_glob ();
00239
00240
00241
00242 trfel_init (argc,argv,NULL);
00243 if (Gtt->nadjelnod == NULL)
00244 Gtt->adjacelem (Outt);
00245
00246
00247 Tp->timecont.be_copy_of(*Adat_old->tctrl);
00248 Tp->time = Tp->timecont.time;
00249 Adat->initialize(i);
00250
00251
00252
00253 Adat->statedata_transfer(Adat_old, Gtt_old);
00254
00255
00256 char tmpfile[255];
00257 sprintf (tmpfile, "%s.r_old.vtk", infile); print_VTK_nodedata (tmpfile, Gtt_old, Adat_old->give_r());
00258 sprintf (tmpfile, "%s.r_new.vtk", infile); print_VTK_nodedata (tmpfile, Gtt, Adat ->give_r());
00259 sprintf (tmpfile, "%s.rdr_old.vtk", infile); print_VTK_nodedata (tmpfile, Gtt_old, Adat_old->give_rdr());
00260 sprintf (tmpfile, "%s.rdr_new.vtk", infile); print_VTK_nodedata (tmpfile, Gtt, Adat ->give_rdr());
00261
00262
00263 delete Gtt_old; Gtt_old = NULL;
00264 delete Adat_old; Adat_old = NULL;
00265
00266
00267
00268
00269
00270
00271
00272
00273 solve_trfel_deterministic_problem ();
00274
00275
00276 i++;
00277
00278
00279
00280 }
00281
00282 break;
00283 }
00284 default:
00285 print_err("unknown problem type is required",__FILE__,__LINE__,__func__); exit (1);
00286 }
00287 }