00001 #include "phsolvert.h"
00002 #include "pglobalt.h"
00003 #include "seqfilest.h"
00004 #include <string.h>
00005 #include <math.h>
00006 #include "mpi.h"
00007
00008 void par_homogenization ()
00009 {
00010 long i,j,n;
00011 double s,zero,alpha,*d,*p,*lhs,*lhsi,*tdlhs,*rhs;
00012
00013 double dt,end_time;
00014
00015
00016 n=Ndoft;
00017
00018 lhs = Lsrst->give_lhs (0);
00019 tdlhs = Lsrst->give_tdlhs (0);
00020 lhsi = Lsrst->lhsi;
00021 rhs = Lsrst->give_rhs (0);
00022
00023 d = new double [n];
00024 p = new double [n];
00025
00026
00027
00028 nullv (lhs,n);
00029 nullv (tdlhs,n);
00030 nullv (d,n);
00031 nullv (p,n);
00032
00033
00034
00035
00036 alpha=Tp->alpha;
00037
00038 zero=Tp->zero;
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 i=0;
00049
00050
00051 Tp->time = Tp->timecont.starttime ();
00052
00053 dt = Tp->timecont.initialtimeincr ();
00054
00055 end_time = Tp->timecont.endtime ();
00056
00057
00058
00059 approximation ();
00060 actual_previous_change ();
00061
00062 Tm->initmaterialmodels();
00063
00064
00065 print_initt(-1, "wt",Ptp->fni,Ptp->fei);
00066 print_stept(0,i,Tp->time,NULL);
00067 print_flusht();
00068
00069
00070 do{
00071
00072
00073 if (Myrank==0) fprintf (stdout,"\n iteration number %ld",i);
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089 if (Myrank!=0){
00090
00091
00092 }
00093
00094
00095 conductivity_matrix (0);
00096
00097 capacity_matrix (0);
00098
00099 if (i==0){
00100 Psolt->computation_stat (Gtt,Kmat);
00101 }
00102
00103
00104
00105
00106
00107
00108
00109 for (j=0;j<n;j++){
00110 p[j]=lhs[j]+(1.0-alpha)*dt*tdlhs[j];
00111 }
00112
00113
00114 Kmat->gmxv (p,d);
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 Kmat->scalgm (alpha*dt);
00126 Kmat->addgm (1.0,*Cmat);
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145 trfel_right_hand_side (0,rhs,n);
00146
00147
00148
00149
00150
00151
00152
00153
00154 for (j=0;j<n;j++){
00155 rhs[j] = rhs[j] - d[j];
00156 d[j]=tdlhs[j];
00157 }
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 Psolt->par_linear_solver (Gtt,Kmat,tdlhs,rhs,Outt,Mesprt);
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 for (j=0;j<n;j++){
00190 s=(1.0-alpha)*d[j]+alpha*tdlhs[j];
00191 lhs[j]+=dt*s;
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
00201 approximation ();
00202
00203
00204
00205
00206 Tp->time = Tp->timecont.newtime ();
00207 dt = Tp->timecont.actualbacktimeincr ();
00208 i++;
00209
00210 print_stept(0,i,Tp->time,NULL);
00211 print_flusht();
00212
00213 printf ("\n Tp->time %e",Tp->time);
00214 printf ("\n i %ld",i);
00215
00216 }while(Tp->time<end_time);
00217
00218
00219
00220 delete [] p; delete [] d;
00221
00222
00223 print_closet();
00224
00225 Psolt->computation_stat_print (Outt);
00226
00227
00228 }
00229