00001 #define EXTERN
00002 #include "globalt.h"
00003 #include "globmatt.h"
00004
00005 #include "stochdrivert.h"
00006 #include "trfelinit.h"
00007 #include "ctlinterface.h"
00008 #include "npglvec.h"
00009 #include "npsolvert.h"
00010
00011
00012 #define DEBUG
00013 #define CTL_Connect
00014 #include "../../ctl/include/ctl/ci/simuTRFEL.ci"
00015
00016 Ctlinterface::Ctlinterface()
00017 {
00018 }
00019
00020
00021
00022
00023
00024 Ctlinterface::Ctlinterface(const std::string& input)
00025 {
00026 stochdrivert *stochd;
00027 stochd = new stochdrivert;
00028
00029 Stt = NULL;
00030 Kmat = NULL;
00031 Cmat = NULL;
00032 Jmat = NULL;
00033 Bmat = NULL;
00034 Lbt = NULL;
00035 Lbat = NULL;
00036 Qbt = NULL;
00037 Qbat = NULL;
00038 Ltt = NULL;
00039 Ltat = NULL;
00040 Lqt = NULL;
00041 Qqt = NULL;
00042 Qqat = NULL;
00043 Lqat = NULL;
00044 Ltett = NULL;
00045 Lht = NULL;
00046 Qht = NULL;
00047 G2d = NULL;
00048
00049 static const char *argv[3];
00050
00051 argv[0] = "trfel";
00052 argv[1] = input.c_str();
00053 argv[2] = NULL;
00054
00055
00056 trfel_init (2,argv,stochd);
00057 }
00058
00059
00060
00061
00062
00063
00064
00065
00066 Ctlinterface::Ctlinterface( const std::string& input, const std::vector<double>& params )
00067 {
00068 }
00069
00070
00071
00072
00073
00074
00075
00076
00077 Ctlinterface::~Ctlinterface()
00078 {
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 fclose(Outt);
00108 }
00109
00110
00111
00112
00113
00114
00115 void Ctlinterface::set_params( const std::vector<double>& params )
00116 {
00117
00118 }
00119
00120
00121
00122
00123
00124
00125
00126 int Ctlinterface::solve()
00127 {
00128 long i,li,ret, j;
00129 double newtime,dt,end_time;
00130 np_glob_vec np_gv;
00131 long lcid = 0;
00132
00133
00134 approximation();
00135 nonstat_solver_init (lcid, np_gv);
00136
00137 dt = Tp->timecont.initialtimeincr ();
00138
00139 end_time = Tp->timecont.endtime ();
00140
00141 li = i = np_gv.istep;
00142 j = 0;
00143
00144 bool breakloop = false;
00145
00146
00147
00148 do{
00149 i++;
00150
00151 newtime = Tp->time = Tp->timecont.newtime();
00152
00153 dt = Tp->timecont.actualbacktimeincr ();
00154
00155 if(Tp->tprob == nonstationary_problem)
00156
00157
00158 ret = one_step_linear(lcid, newtime, dt, i, li, np_gv);
00159
00160 if(Tp->tprob == nonlinear_nonstationary_problem)
00161
00162
00163 ret = one_step_nonlinear(lcid, newtime, dt, i, li, np_gv);
00164
00165
00166 if (Tp->adaptivityflag && Tp->time < end_time)
00167 if (i%2 == 0)
00168 breakloop = Adat->run (2, true);
00169 j++;
00170
00171 }while(Tp->time < end_time && breakloop == false);
00172
00173 return (Tp->time-end_time < 0.0);
00174 }
00175
00176
00177
00178
00179
00180
00181 void Ctlinterface::get_state( std::vector<double>& state ) const
00182 {
00183 long nt;
00184 nt = Tp->timecont.starttime() - Tp->timecont.endtime();
00185 nt /=Tp->timecont.initialtimeincr ();
00186
00187
00188 state.resize(Ndoft*(nt+1));
00189 double *statep= new double[Ndoft*(nt+1)];
00190
00191 Stt->exportvalues((double*) statep);
00192 for(int i=0; i<Ndoft*(nt+1); i++)
00193 state[i] = statep[i];
00194 }
00195
00196 struct connectDetail
00197 {
00198 CTL_Constructor(1, (const std::string&), 1);
00199 CTL_Constructor(2, (const std::string&, const std::vector<double>&), 2);
00200 };
00201
00202 void CTL_connect()
00203 {
00204 ctl::connect<simuCI, Ctlinterface, connectDetail>();
00205 }
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 void Ctlinterface::get_residual(const std::vector<double>& params,const std::vector<double>& state,
00217 const unsigned int niter, std::vector<double>& residual ) const
00218 {
00219 long i, j, n=Ndoft, nt;
00220 double *aux, *lhs, *du, dt;
00221 aux = new double [n];
00222 du = new double [n];
00223
00224
00225 lhs = Lsrst->give_lhs(0);
00226
00227 Stt->importvalues((double*) ¶ms[0]);
00228
00229 dt = Tp->timecont.initialtimeincr ();
00230
00231 nt = Tp->timecont.starttime() - Tp->timecont.endtime();
00232 nt /= dt;
00233
00234 for (i=0; i<nt+1; i++)
00235 {
00236
00237
00238 for (j=0;j<n;j++){
00239 lhs[j]=state[n*i+j];
00240 }
00241
00242 subv(lhs, (double*) &state[n*(i+1)], du, n);
00243 residuum((double*) &residual[n*i], du, aux, dt, n, 0);
00244 }
00245
00246
00247
00248
00249
00250
00251 delete [] du;
00252 delete [] aux;
00253 }
00254