source: issm/trunk-jpl/src/c/objects/Inputs/TransientInput.cpp@ 12326

Last change on this file since 12326 was 12326, checked in by Mathieu Morlighem, 13 years ago

merged trunk and trunk-jpl

File size: 11.2 KB
RevLine 
[8757]1/*!\file TransientInput.c
2 * \brief: implementation of the TransientInput object
3 */
[12326]4/*Headers{{{*/
[8757]5#ifdef HAVE_CONFIG_H
[9320]6 #include <config.h>
[8757]7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
[9320]11#include <stdio.h>
[8757]12#include <string.h>
13#include "../objects.h"
14#include "../../EnumDefinitions/EnumDefinitions.h"
15#include "../../shared/shared.h"
16#include "../../Container/Container.h"
17#include "../../include/include.h"
18/*}}}*/
19
20/*TransientInput constructors and destructor*/
[12326]21/*FUNCTION TransientInput::TransientInput(){{{*/
[8757]22TransientInput::TransientInput(){
[8791]23
[8967]24 enum_type=UNDEF;
[8757]25 inputs=NULL;
[8791]26 this->numtimesteps=0;
27 this->parameters=NULL;
28 this->timesteps=NULL;
29
[8757]30}
31/*}}}*/
[12326]32/*FUNCTION TransientInput::TransientInput(int in_enum_type){{{*/
[8757]33TransientInput::TransientInput(int in_enum_type)
34{
35 /*Set Enum*/
36 enum_type=in_enum_type;
37
38 /*Allocate values and timesteps, and copy: */
39 this->numtimesteps=0;
40 this->timesteps=NULL;
41 inputs = new Inputs();
42 this->parameters=NULL;
43
44}
45/*}}}*/
[12326]46/*FUNCTION TransientInput::~TransientInput{{{*/
[8757]47TransientInput::~TransientInput(){
48 xfree((void**)&this->timesteps);
[8791]49 this->timesteps=NULL;
[8757]50 this->numtimesteps=0;
51 parameters=NULL;
52 delete this->inputs;
53 return;
54}
55/*}}}*/
56
57/*Object virtual functions definitions:*/
[12326]58/*FUNCTION TransientInput::Echo {{{*/
[8757]59void TransientInput::Echo(void){
60 this->DeepEcho();
61}
62/*}}}*/
[12326]63/*FUNCTION TransientInput::DeepEcho{{{*/
[8757]64void TransientInput::DeepEcho(void){
65
66 int i;
67
68 printf("TransientInput:\n");
69 printf(" enum: %i (%s)\n",this->enum_type,EnumToStringx(this->enum_type));
70 printf(" numtimesteps: %i\n",this->numtimesteps);
71 printf("---inputs: \n");
72 for(i=0;i<this->numtimesteps;i++){
73 printf(" time: %g \n",this->timesteps[i]);
74 ((Input*)this->inputs->GetObjectByOffset(i))->Echo();
75 }
76}
77/*}}}*/
[12326]78/*FUNCTION TransientInput::Id{{{*/
[8757]79int TransientInput::Id(void){ return -1; }
80/*}}}*/
[12326]81/*FUNCTION TransientInput::MyRank{{{*/
[8757]82int TransientInput::MyRank(void){
83 extern int my_rank;
84 return my_rank;
85}
86/*}}}*/
[12326]87/*FUNCTION TransientInput::ObjectEnum{{{*/
[9883]88int TransientInput::ObjectEnum(void){
[8757]89
90 return TransientInputEnum;
91
92}
93/*}}}*/
[12326]94/*FUNCTION TransientInput::copy{{{*/
[8757]95Object* TransientInput::copy() {
96
97 TransientInput* output=NULL;
98
99 output = new TransientInput();
100 output->enum_type=this->enum_type;
101 output->numtimesteps=this->numtimesteps;
[8791]102 output->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
[8757]103 memcpy(output->timesteps,this->timesteps,this->numtimesteps*sizeof(double));
104 output->inputs=(Inputs*)this->inputs->Copy();
105 output->parameters=this->parameters;
106
107 return output;
108
109}
110/*}}}*/
111
112/*TransientInput management*/
[12326]113/*FUNCTION TransientInput::InstanceEnum{{{*/
[9883]114int TransientInput::InstanceEnum(void){
[8757]115
116 return this->enum_type;
117
118}
119/*}}}*/
[12326]120/*FUNCTION TransientInput::SpawnTriaInput{{{*/
[8757]121Input* TransientInput::SpawnTriaInput(int* indices){
122
123 /*output*/
124 TransientInput* outinput=NULL;
125
126 /*Create new Transientinput (copy of current input)*/
127 outinput=new TransientInput();
128 outinput->enum_type=this->enum_type;
[8791]129 outinput->numtimesteps=this->numtimesteps;
130 outinput->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
131 memcpy(outinput->timesteps,this->timesteps,this->numtimesteps*sizeof(double));
[8757]132 outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices);
[8791]133 outinput->parameters=this->parameters;
[8757]134
135 /*Assign output*/
136 return outinput;
137
138}
139/*}}}*/
[12326]140/*FUNCTION TransientInput::SpawnResult{{{*/
[8757]141ElementResult* TransientInput::SpawnResult(int step, double time){
142
143 ElementResult* elementresult=NULL;
144
145 /*Ok, we want to spawn an ElementResult. We have the time, just get
146 *the correct values: */
147 Input* input=GetTimeInput(time);
148
149 elementresult=input->SpawnResult(step,time);
150
[8791]151 delete input;
[8757]152
153 return elementresult;
154}
155/*}}}*/
156
157/*Object functions*/
[12326]158/*FUNCTION TransientInput::GetInputValue(double* pvalue,GaussTria* gauss){{{*/
[10135]159void TransientInput::GetInputValue(double* pvalue,GaussTria* gauss){
[8757]160 double time;
161
162 /*First, recover current time from parameters: */
[8791]163 this->parameters->FindParam(&time,TimeEnum);
[8757]164
165 /*Retrieve interpolated values for this time step: */
166 Input* input=GetTimeInput(time);
167
168 /*Call input function*/
[10135]169 input->GetInputValue(pvalue,gauss);
[8757]170
[8791]171 delete input;
[12326]172}
173/*}}}*/
174/*FUNCTION TransientInput::GetInputValue(double* pvalue,GaussTria* gauss,double time){{{*/
175void TransientInput::GetInputValue(double* pvalue,GaussTria* gauss,double time){
[8757]176
[12326]177 /*Retrieve interpolated values for this time step: */
178 Input* input=GetTimeInput(time);
179
180 /*Call input function*/
181 input->GetInputValue(pvalue,gauss);
182
183 delete input;
[8757]184}
185/*}}}*/
[12326]186/*FUNCTION TransientInput::GetInputDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{*/
[10135]187void TransientInput::GetInputDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){
[8757]188
189 double time;
190
191 /*First, recover current time from parameters: */
192 parameters->FindParam(&time,TimeEnum);
193
194 /*Retrieve interpolated values for this time step: */
195 Input* input=GetTimeInput(time);
196
197 /*Call input function*/
[10135]198 input->GetInputDerivativeValue(p,xyz_list,gauss);
[8757]199
[8791]200 delete input;
201
[8757]202}
203/*}}}*/
[12326]204/*FUNCTION TransientInput::ChangeEnum{{{*/
[8757]205void TransientInput::ChangeEnum(int newenumtype){
206 this->enum_type=newenumtype;
207}
208/*}}}*/
[12326]209/*FUNCTION TransientInput::GetInputAverage{{{*/
[10135]210void TransientInput::GetInputAverage(double* pvalue){
[8757]211
212 double time;
213
214 /*First, recover current time from parameters: */
215 parameters->FindParam(&time,TimeEnum);
216
217 /*Retrieve interpolated values for this time step: */
218 Input* input=GetTimeInput(time);
219
220 /*Call input function*/
[10135]221 input->GetInputAverage(pvalue);
[8757]222
[8791]223 delete input;
[8757]224
225}
226/*}}}*/
227
228/*Intermediary*/
[12326]229/*FUNCTION TransientInput::AddTimeInput{{{*/
230void TransientInput::AddTimeInput(Input* input,double time){
231
232 /*insert values at time step: */
233 if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _assert_("timestep values must increase sequentially");
234
235 //copy timesteps, add the new time, delete previous timesteps, and add the new input: inputs->AddObject(input);
236 double* old_timesteps=NULL;
237
238 if (this->numtimesteps > 0){
239 old_timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
240 memcpy(old_timesteps,this->timesteps,this->numtimesteps*sizeof(double));
241 xfree((void**)&this->timesteps);
242 }
243
244 this->numtimesteps=this->numtimesteps+1;
245 this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
246
247 if (this->numtimesteps > 1){
248 memcpy(this->timesteps,old_timesteps,(this->numtimesteps-1)*sizeof(double));
249 xfree((void**)&old_timesteps);
250 }
251
252 /*go ahead and plug: */
253 this->timesteps[this->numtimesteps-1]=time;
254 inputs->AddObject(input);
255
256}
257/*}}}*/
258/*FUNCTION TransientInput::SquareMin{{{*/
[8757]259void TransientInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
260
261 double time;
262
263 /*First, recover current time from parameters: */
264 parameters->FindParam(&time,TimeEnum);
265
266 /*Retrieve interpolated values for this time step: */
267 Input* input=GetTimeInput(time);
268
269 /*Call input function*/
270 input->SquareMin(psquaremin,process_units,parameters);
271
[8791]272 delete input;
[8757]273
274}
275/*}}}*/
[12326]276/*FUNCTION TransientInput::InfinityNorm{{{*/
[8757]277double TransientInput::InfinityNorm(void){
278
279 double time;
280 double infnorm;
281
282 /*First, recover current time from parameters: */
283 parameters->FindParam(&time,TimeEnum);
284
285 /*Retrieve interpolated values for this time step: */
286 Input* input=GetTimeInput(time);
287
288 /*Call input function*/
289 infnorm=input->InfinityNorm();
290
[11197]291 /*Clean-up and return*/
[8791]292 delete input;
[11197]293 return infnorm;
[8757]294}
295/*}}}*/
[12326]296/*FUNCTION TransientInput::Max{{{*/
[8757]297double TransientInput::Max(void){
298
299 double time;
300 double max;
301
302 /*First, recover current time from parameters: */
303 parameters->FindParam(&time,TimeEnum);
304
305 /*Retrieve interpolated values for this time step: */
306 Input* input=GetTimeInput(time);
307
308 /*Call input function*/
309 max=input->Max();
310
[8791]311 delete input;
[8757]312
313 return max;
314}
315/*}}}*/
[12326]316/*FUNCTION TransientInput::MaxAbs{{{*/
[8757]317double TransientInput::MaxAbs(void){
318
319 double time;
320 double maxabs;
321
322 /*First, recover current time from parameters: */
323 parameters->FindParam(&time,TimeEnum);
324
325 /*Retrieve interpolated values for this time step: */
326 Input* input=GetTimeInput(time);
327
328 /*Call input function*/
329 maxabs=input->MaxAbs();
330
[11197]331 /*Clean-up and return*/
[8791]332 delete input;
[8757]333 return maxabs;
[8791]334
[8757]335}
336/*}}}*/
[12326]337/*FUNCTION TransientInput::Min{{{*/
[8757]338double TransientInput::Min(void){
339
340 double time;
341 double min;
342
343 /*First, recover current time from parameters: */
344 parameters->FindParam(&time,TimeEnum);
345
346 /*Retrieve interpolated values for this time step: */
347 Input* input=GetTimeInput(time);
348
349 /*Call input function*/
350 min=input->Min();
351
[11197]352 /*Clean-up and return*/
[8791]353 delete input;
[8757]354 return min;
[8791]355
[8757]356}
357/*}}}*/
[12326]358/*FUNCTION TransientInput::MinAbs{{{*/
[8757]359double TransientInput::MinAbs(void){
360
361 double time;
362 double minabs;
363
364 /*First, recover current time from parameters: */
365 parameters->FindParam(&time,TimeEnum);
366
367 /*Retrieve interpolated values for this time step: */
368 Input* input=GetTimeInput(time);
369
370 /*Call input function*/
371 minabs=input->MinAbs();
372
[11197]373 /*Clean-up and return*/
[8791]374 delete input;
[8757]375 return minabs;
376}
377/*}}}*/
[12326]378/*FUNCTION TransientInput::GetVectorFromInputs{{{*/
[11695]379void TransientInput::GetVectorFromInputs(Vector* vector,int* doflist){
[8757]380
381 double time;
382
383 /*First, recover current time from parameters: */
384 parameters->FindParam(&time,TimeEnum);
385
386 /*Retrieve interpolated values for this time step: */
387 Input* input=GetTimeInput(time);
388
389 /*Call input function*/
390 input->GetVectorFromInputs(vector,doflist);
391
[8791]392 delete input;
[8757]393
394} /*}}}*/
[12326]395/*FUNCTION TransientInput::GetTimeInput{{{*/
[8757]396Input* TransientInput::GetTimeInput(double intime){
397
398 int i,j;
399 double deltat;
400 double alpha1,alpha2;
401 bool found=false;
402 Input* input=NULL;
403 Input* input1=NULL;
404 Input* input2=NULL;
405
406 /*Ok, we have the time, go through the timesteps, and figure out which interval we
407 *fall within. Then interpolate the values on this interval: */
408 if(intime<this->timesteps[0]){
409 /*get values for the first time: */
410 input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
411 found=true;
412 }
413 else if(intime>this->timesteps[this->numtimesteps-1]){
414 /*get values for the last time: */
415 input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
416 found=true;
417 }
418 else{
419 /*Find which interval we fall within: */
420 for(i=0;i<this->numtimesteps;i++){
421 if(intime==this->timesteps[i]){
422 /*We are right on one step time: */
423 input=(Input*)((Input*)this->inputs->GetObjectByOffset(i))->copy();
424 found=true;
425 break; //we are done with the time interpolation.
426 }
427 else{
428 if(this->timesteps[i]<intime && intime<this->timesteps[i+1]){
429 /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
430 deltat=this->timesteps[i+1]-this->timesteps[i];
431 alpha2=(intime-this->timesteps[i])/deltat;
432 alpha1=(1-alpha2);
433
434 input1=(Input*)this->inputs->GetObjectByOffset(i);
435 input2=(Input*)this->inputs->GetObjectByOffset(i+1);
436
437 input=(Input*)input1->copy();
438 input->Scale(alpha1);
439 input->AXPY(input2,alpha2);
440
441 found=true;
442 break;
443 }
444 else continue; //keep looking on the next interval
445 }
446 }
447 }
448 if(!found)_error_("did not find time interval on which to interpolate forcing values!");
449
450 /*Assign output pointer*/
451 return input;
452}
453/*}}}*/
[12326]454/*FUNCTION TransientInput::Configure{{{*/
[8757]455void TransientInput::Configure(Parameters* parameters){
456 this->parameters=parameters;
457}
458/*}}}*/
Note: See TracBrowser for help on using the repository browser.