00001 #include "homogtrans.h"
00002 #include "spsolvert.h"
00003 #include "nspsolvert.h"
00004 #include "globalt.h"
00005 #include "globmatt.h"
00006 #include "transprint.h"
00007 #include <string.h>
00008 #include "npsolvert.h"
00009
00010
00011
00012
00013
00014
00015 void transport_homogenization ()
00016 {
00017 long i,j;
00018 long ntm,ncomp;
00019 double puc_area;
00020 double *rhs,*lhs;
00021 matrix ltm,lm,kpuc,km,kM,kinv,khelp,cM;
00022
00023
00024 Tm->initmaterialmodels();
00025
00026 approximation_puc ();
00027
00028
00029
00030
00031 switch (Tp->tprob){
00032 case stationary_problem:{
00033
00034 solve_stationary_problem ();
00035 break;
00036 }
00037 case nonlinear_stationary_problem:{
00038 solve_nonlinear_stationary_problem_pokus ();
00039
00040 break;
00041 }
00042 default:{
00043 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00044 }
00045 }
00046
00047
00048
00049
00050 approximation_puc ();
00051
00052
00053
00054
00055
00056 ntm = Tp->ntm;
00057 ncomp = Tt->nodes[0].ncompgrad;
00058 allocm (Ndoft,ntm*ncomp,ltm);
00059 allocm (ntm*ncomp,Ndoft,lm);
00060 allocm (Ndoft,Ndoft,kpuc);
00061 allocm (Ndoft,Ndoft,kinv);
00062 allocm (Ndoft,ntm*ncomp,khelp);
00063
00064 assemble_l_matrix (lm,ltm);
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087 conductivity_matrix (0);
00088
00089
00090 allocm (ntm*ncomp,ntm*ncomp,km);
00091
00092 assemble_average_d_matrix(km,puc_area);
00093
00094
00095 allocm (ntm*ncomp,ntm*ncomp,kM);
00096
00097
00098
00099
00100 lhs = new double [kinv.m];
00101 nullv (lhs,kinv.m);
00102
00103 rhs = new double [kinv.m];
00104 nullv (rhs,kinv.m);
00105
00106 for(i=0;i<kinv.m;i++){
00107 nullv (rhs,kinv.m);
00108 nullv (lhs,kinv.m);
00109 rhs[i] = 1.0;
00110 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
00111 for(j=0;j<kinv.n;j++){
00112 kinv[j][i]=lhs[j];
00113 }
00114 }
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 mxm(kinv,ltm,khelp);
00126 mxm(lm,khelp,kM);
00127 cmulm(-1.0,kM);
00128 puc_area = 1.0/puc_area;
00129 cmulm(puc_area,kM);
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 addm(km,kM,kM);
00151
00152
00153 allocm (ntm,ntm,cM);
00154
00155
00156 assemble_average_c_matrix(cM);
00157
00158
00159 if (Mesprt != 0){
00160 fprintf(Outt,"\n\n vysledna matice D^M\n");
00161 for (i=0;i<ntm*ncomp;i++){
00162 for (j=0;j<ntm*ncomp;j++)
00163 fprintf(Outt,"%e \n",kM[i][j]);
00164 fprintf(Outt,"\n");
00165 }
00166
00167 fprintf(Outt,"\n\n vysledna matice C^M\n");
00168 for (i=0;i<ntm;i++){
00169 for (j=0;j<ntm;j++)
00170 fprintf(Outt,"%e \n",cM[i][j]);
00171 fprintf(Outt,"\n");
00172 }
00173 }
00174
00175
00176 destrm (lm);
00177 destrm (ltm);
00178 destrm (kpuc);
00179 destrm (kinv);
00180 destrm (km);
00181 destrm (kM);
00182 destrm (cM);
00183 destrm (khelp);
00184 delete [] rhs;
00185 delete [] lhs;
00186 }
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197 void paral_transport_homogenization (double *unkn_r,double *matrix_k,double *matrix_c)
00198 {
00199 long i,j,k;
00200 long ntm,ncomp;
00201 double puc_area;
00202 double *rhs,*lhs;
00203 matrix ltm,lm,kpuc,km,kM,kinv,khelp,cM;
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 Tm->initmaterialmodels();
00214
00215
00216 ntm = Tp->ntm;
00217 ncomp = Tt->nodes[0].ncompgrad;
00218 k = -1;
00219 for (i=0;i<ntm;i++){
00220 k++;
00221 Tb->lc[i].masterval = unkn_r[k];
00222 for (j=0;j<ncomp;j++){
00223 k++;
00224 Tb->lc[i].mastergrad[j] = unkn_r[k];
00225 }
00226 }
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 approximation_puc ();
00239
00240
00241
00242 switch (Tp->tprob){
00243 case stationary_problem:{
00244
00245 solve_stationary_problem ();
00246 break;
00247 }
00248 case nonlinear_stationary_problem:{
00249 solve_nonlinear_stationary_problem_pokus ();
00250
00251 break;
00252 }
00253 default:{
00254 print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
00255 }
00256 }
00257
00258
00259
00260
00261 approximation_puc ();
00262
00263
00264
00265
00266
00267 allocm (Ndoft,ntm*ncomp,ltm);
00268 allocm (ntm*ncomp,Ndoft,lm);
00269 allocm (Ndoft,Ndoft,kpuc);
00270 allocm (Ndoft,Ndoft,kinv);
00271 allocm (Ndoft,ntm*ncomp,khelp);
00272
00273 assemble_l_matrix (lm,ltm);
00274
00275
00276
00277 conductivity_matrix (0);
00278
00279
00280 allocm (ntm*ncomp,ntm*ncomp,km);
00281
00282 assemble_average_d_matrix(km,puc_area);
00283
00284
00285 allocm (ntm*ncomp,ntm*ncomp,kM);
00286
00287
00288
00289
00290 lhs = new double [kinv.m];
00291 nullv (lhs,kinv.m);
00292
00293 rhs = new double [kinv.m];
00294 nullv (rhs,kinv.m);
00295
00296 for(i=0;i<kinv.m;i++){
00297 nullv (rhs,kinv.m);
00298 nullv (lhs,kinv.m);
00299 rhs[i] = 1.0;
00300 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
00301 for(j=0;j<kinv.n;j++){
00302 kinv[j][i]=lhs[j];
00303 }
00304 }
00305
00306 mxm(kinv,ltm,khelp);
00307 mxm(lm,khelp,kM);
00308 cmulm(-1.0,kM);
00309 puc_area = 1.0/puc_area;
00310 cmulm(puc_area,kM);
00311
00312
00313 addm(km,kM,kM);
00314
00315
00316 allocm (ntm,ntm,cM);
00317
00318
00319 assemble_average_c_matrix(cM);
00320
00321
00322
00323 k = 0;
00324 for (i=0;i<ntm*ncomp;i++){
00325 for (j=0;j<ntm*ncomp;j++){
00326 matrix_k[k] = kM[i][j];
00327 k++;
00328 }
00329 }
00330
00331
00332 k = 0;
00333 for (i=0;i<ntm;i++){
00334 for (j=0;j<ntm;j++){
00335 matrix_c[k] = cM[i][j];
00336 k++;
00337 }
00338 }
00339
00340 destrm (lm);
00341 destrm (ltm);
00342 destrm (kpuc);
00343 destrm (kinv);
00344 destrm (km);
00345 destrm (kM);
00346 destrm (cM);
00347 destrm (khelp);
00348 delete [] rhs;
00349 delete [] lhs;
00350
00351 }