source: issm/trunk/src/c/objects/Inputs/TriaVertexForcing.cpp@ 8371

Last change on this file since 8371 was 8371, checked in by Mathieu Morlighem, 14 years ago

trunk: fixed marshall/demarshall of TriaVertecForcing
File size: 14.2 KB
Line 
1/*!\file TriaVertexForcing.c
2 * \brief: implementation of the TriaVertexForcing object
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include "config.h"
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11#include "stdio.h"
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/*TriaVertexForcing constructors and destructor*/
20/*FUNCTION TriaVertexForcing::TriaVertexForcing(){{{1*/
21TriaVertexForcing::TriaVertexForcing(){
22 enum_type=NoneEnum;
23 values=NULL;
24 numtimesteps=0;
25 parameters=NULL;
26 return;
27}
28/*}}}*/
29/*FUNCTION TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* time0values,double time0,int numtimesteps,Parameters* parameters){{{1*/
30TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* time0values,double time0,int in_numtimesteps,Parameters* parameters)
31 :TriaRef(1)
32{
33
34 /*Set TriaRef*/
35 this->SetElementType(P1Enum,0);
36 this->element_type=P1Enum;
37
38 /*Set Enum*/
39 enum_type=in_enum_type;
40
41 /*Allocate values and timesteps, and set first time step values: */
42 this->numtimesteps=in_numtimesteps;
43 this->values=(double*)xmalloc(this->numtimesteps*3*sizeof(double));
44 this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
45 this->values[0]=time0values[0];
46 this->values[1]=time0values[1];
47 this->values[2]=time0values[2];
48 this->timesteps[0]=time0;
49
50 this->parameters=parameters;
51
52}
53/*}}}*/
54/*FUNCTION TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* in_values,double* in_time,int numtimesteps,Parameters* parameters){{{1*/
55TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* in_values,double* in_time,int in_numtimesteps,Parameters* parameters)
56 :TriaRef(1)
57{
58
59 /*Set TriaRef*/
60 this->SetElementType(P1Enum,0);
61 this->element_type=P1Enum;
62
63 /*Set Enum*/
64 enum_type=in_enum_type;
65
66 /*Allocate values and timesteps, and copy: */
67 this->numtimesteps=in_numtimesteps;
68 this->values=(double*)xmalloc(this->numtimesteps*3*sizeof(double));
69 this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
70
71 memcpy(this->values,in_values,numtimesteps*3*sizeof(double));
72 memcpy(this->timesteps,in_time,numtimesteps*sizeof(double));
73
74 this->parameters=parameters;
75
76}
77/*}}}*/
78/*FUNCTION TriaVertexForcing::~TriaVertexForcing(){{{1*/
79TriaVertexForcing::~TriaVertexForcing(){
80 xfree((void**)&this->values);
81 xfree((void**)&this->timesteps);
82 parameters=NULL;
83 return;
84}
85/*}}}*/
86/*FUNCTION TriaVertexForcing::AddTimeValues(double* values,int step,double time);{{{1*/
87void TriaVertexForcing::AddTimeValues(double* values,int step,double time){
88
89 /*insert values at time step: */
90 if (step<0)_error_("timestep should start at least at time 0!");
91 if (step>numtimesteps-1)_error_("timestep cannot be more than number of time steps");
92
93 /*go ahead and plug: */
94 this->timesteps[step]=time;
95 this->values[3*step+0]=values[0];
96 this->values[3*step+1]=values[1];
97 this->values[3*step+2]=values[2];
98}
99/*}}}*/
100
101/*Object virtual functions definitions:*/
102/*FUNCTION TriaVertexForcing::Echo {{{1*/
103void TriaVertexForcing::Echo(void){
104 this->DeepEcho();
105}
106/*}}}*/
107/*FUNCTION TriaVertexForcing::DeepEcho{{{1*/
108void TriaVertexForcing::DeepEcho(void){
109
110 int i;
111
112 printf("TriaVertexForcing:\n");
113 printf(" enum: %i (%s)\n",this->enum_type,EnumToStringx(this->enum_type));
114 printf(" numtimesteps: %i\n",this->numtimesteps);
115 for(i=0;i<this->numtimesteps;i++){
116 printf(" time: %g values: [%g %g %g]\n",this->timesteps[i],this->values[3*i+0],this->values[3*i+1],this->values[3*i+2]);
117 }
118}
119/*}}}*/
120/*FUNCTION TriaVertexForcing::Id{{{1*/
121int TriaVertexForcing::Id(void){ return -1; }
122/*}}}*/
123/*FUNCTION TriaVertexForcing::MyRank{{{1*/
124int TriaVertexForcing::MyRank(void){
125 extern int my_rank;
126 return my_rank;
127}
128/*}}}*/
129/*FUNCTION TriaVertexForcing::Marshall{{{1*/
130void TriaVertexForcing::Marshall(char** pmarshalled_dataset){
131
132 char* marshalled_dataset=NULL;
133 int enum_value=0;
134
135 /*recover marshalled_dataset: */
136 marshalled_dataset=*pmarshalled_dataset;
137
138 /*get enum value of TriaVertexForcing: */
139 enum_value=TriaVertexForcingEnum;
140
141 /*marshall enum: */
142 memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
143
144 /*marshall TriaVertexForcing data: */
145 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
146 memcpy(marshalled_dataset,&numtimesteps,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps);
147 memcpy(marshalled_dataset,values,numtimesteps*3*sizeof(double));marshalled_dataset+=numtimesteps*3*sizeof(double);
148 memcpy(marshalled_dataset,timesteps,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double);
149
150 *pmarshalled_dataset=marshalled_dataset;
151}
152/*}}}*/
153/*FUNCTION TriaVertexForcing::MarshallSize{{{1*/
154int TriaVertexForcing::MarshallSize(){
155
156 return
157 +sizeof(enum_type)+
158 +sizeof(numtimesteps)+
159 3*numtimesteps*sizeof(double)+
160 numtimesteps*sizeof(double)+
161 +sizeof(int); //sizeof(int) for enum value
162}
163/*}}}*/
164/*FUNCTION TriaVertexForcing::Demarshall{{{1*/
165void TriaVertexForcing::Demarshall(char** pmarshalled_dataset){
166
167 char* marshalled_dataset=NULL;
168 int i;
169
170 /*recover marshalled_dataset: */
171 marshalled_dataset=*pmarshalled_dataset;
172
173 /*this time, no need to get enum type, the pointer directly points to the beginning of the
174 *object data (thanks to DataSet::Demarshall):*/
175 memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
176 memcpy(&numtimesteps,marshalled_dataset,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps);
177 /*allocate: */
178 timesteps=(double*)xmalloc(numtimesteps*sizeof(double));
179 values=(double*)xmalloc(numtimesteps*3*sizeof(double));
180
181 memcpy(values,marshalled_dataset,3*numtimesteps*sizeof(double));marshalled_dataset+=3*numtimesteps*sizeof(double);
182 memcpy(timesteps,marshalled_dataset,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double);
183
184 /*return: */
185 *pmarshalled_dataset=marshalled_dataset;
186 return;
187}
188/*}}}*/
189/*FUNCTION TriaVertexForcing::Enum{{{1*/
190int TriaVertexForcing::Enum(void){
191
192 return TriaVertexForcingEnum;
193
194}
195/*}}}*/
196/*FUNCTION TriaVertexForcing::copy{{{1*/
197Object* TriaVertexForcing::copy() {
198
199 return new TriaVertexForcing(this->enum_type,this->values,this->timesteps,this->numtimesteps,this->parameters);
200
201}
202/*}}}*/
203
204/*TriaVertexForcing management*/
205/*FUNCTION TriaVertexForcing::EnumType{{{1*/
206int TriaVertexForcing::EnumType(void){
207
208 return this->enum_type;
209
210}
211/*}}}*/
212/*FUNCTION TriaVertexForcing::SpawnTriaInput{{{1*/
213Input* TriaVertexForcing::SpawnTriaInput(int* indices){
214
215 /*output*/
216 TriaVertexForcing* outforcing=NULL;
217
218 /*Create new Tria forcing (copy of current forcing)*/
219 outforcing=new TriaVertexForcing(this->enum_type,this->values,this->timesteps,this->numtimesteps,this->parameters);
220
221 /*Assign output*/
222 return outforcing;
223
224}
225/*}}}*/
226/*FUNCTION TriaVertexForcing::SpawnResult{{{1*/
227ElementResult* TriaVertexForcing::SpawnResult(int step, double time){
228
229 double triavalues[3];
230
231 /*Ok, we want to spawn a TriaVertexElementResult. We have the time, just get
232 *the correct values at the three nodes: */
233 this->GetTimeValues(&triavalues[0],time);
234
235 return new TriaVertexElementResult(this->enum_type,&triavalues[0],step,time);
236
237}
238/*}}}*/
239
240
241/*Object functions*/
242/*FUNCTION TriaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
243void TriaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){
244
245 double time;
246 double triavalues[3];
247
248 /*First, recover current time from parameters: */
249 parameters->FindParam(&time,TimeEnum);
250
251 /*Retrieve interpolated values for this time step: */
252 this->GetTimeValues(&triavalues[0],time);
253
254 /*Call TriaRef function*/
255 TriaRef::GetParameterValue(pvalue,&triavalues[0],gauss);
256
257}
258/*}}}*/
259/*FUNCTION TriaVertexForcing::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/
260void TriaVertexForcing::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){
261
262 double time;
263 double triavalues[3];
264
265 /*First, recover current time from parameters: */
266 parameters->FindParam(&time,TimeEnum);
267
268 /*Retrieve interpolated values for this time step: */
269 this->GetTimeValues(&triavalues[0],time);
270
271 /*Call TriaRef function*/
272 TriaRef::GetParameterDerivativeValue(p,&triavalues[0],xyz_list,gauss);
273}
274/*}}}*/
275/*FUNCTION TriaVertexForcing::ChangeEnum{{{1*/
276void TriaVertexForcing::ChangeEnum(int newenumtype){
277 this->enum_type=newenumtype;
278}
279/*}}}*/
280/*FUNCTION TriaVertexForcing::GetParameterAverage{{{1*/
281void TriaVertexForcing::GetParameterAverage(double* pvalue){
282
283 double time;
284 double triavalues[3];
285
286 /*First, recover current time from parameters: */
287 parameters->FindParam(&time,TimeEnum);
288
289 /*Retrieve interpolated values for this time step: */
290 this->GetTimeValues(&triavalues[0],time);
291
292 *pvalue=1./3.*(triavalues[0]+triavalues[1]+triavalues[2]);
293}
294/*}}}*/
295
296/*Intermediary*/
297/*FUNCTION TriaVertexForcing::SquareMin{{{1*/
298void TriaVertexForcing::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
299
300 int i;
301 const int numnodes=3;
302 double valuescopy[numnodes];
303 double squaremin;
304 double time;
305
306 /*First, recover current time from parameters: */
307 parameters->FindParam(&time,TimeEnum);
308
309 /*Retrieve interpolated values for this time step: */
310 this->GetTimeValues(&valuescopy[0],time);
311
312 /*Process units if requested: */
313 if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type,parameters);
314
315 /*Now, figure out minimum of valuescopy: */
316 squaremin=pow(valuescopy[0],2);
317 for(i=1;i<numnodes;i++){
318 if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
319 }
320 /*Assign output pointers:*/
321 *psquaremin=squaremin;
322}
323/*}}}*/
324/*FUNCTION TriaVertexForcing::InfinityNorm{{{1*/
325double TriaVertexForcing::InfinityNorm(void){
326
327 /*Output*/
328 double norm=0;
329 const int numnodes=3;
330 double time;
331 double triavalues[3];
332
333 /*First, recover current time from parameters: */
334 parameters->FindParam(&time,TimeEnum);
335
336 /*Retrieve interpolated values for this time step: */
337 this->GetTimeValues(&triavalues[0],time);
338
339 for(int i=0;i<numnodes;i++) if(fabs(triavalues[i])>norm) norm=fabs(triavalues[i]);
340 return norm;
341}
342/*}}}*/
343/*FUNCTION TriaVertexForcing::Max{{{1*/
344double TriaVertexForcing::Max(void){
345
346 const int numnodes=3;
347 double max;
348 double time;
349 double triavalues[3];
350
351 /*First, recover current time from parameters: */
352 parameters->FindParam(&time,TimeEnum);
353
354 /*Retrieve interpolated values for this time step: */
355 this->GetTimeValues(&triavalues[0],time);
356
357 max=triavalues[0];
358
359 for(int i=1;i<numnodes;i++){
360 if(triavalues[i]>max) max=triavalues[i];
361 }
362 return max;
363}
364/*}}}*/
365/*FUNCTION TriaVertexForcing::MaxAbs{{{1*/
366double TriaVertexForcing::MaxAbs(void){
367
368 const int numnodes=3;
369 double max;
370 double time;
371 double triavalues[3];
372
373 /*First, recover current time from parameters: */
374 parameters->FindParam(&time,TimeEnum);
375
376 /*Retrieve interpolated values for this time step: */
377 this->GetTimeValues(&triavalues[0],time);
378
379
380 max=fabs(triavalues[0]);
381
382 for(int i=1;i<numnodes;i++){
383 if(fabs(triavalues[i])>max) max=fabs(triavalues[i]);
384 }
385 return max;
386}
387/*}}}*/
388/*FUNCTION TriaVertexForcing::Min{{{1*/
389double TriaVertexForcing::Min(void){
390
391 const int numnodes=3;
392 double min;
393 double time;
394 double triavalues[3];
395
396 /*First, recover current time from parameters: */
397 parameters->FindParam(&time,TimeEnum);
398
399 /*Retrieve interpolated values for this time step: */
400 this->GetTimeValues(&triavalues[0],time);
401
402 min=triavalues[0];
403
404 for(int i=1;i<numnodes;i++){
405 if(triavalues[i]<min) min=triavalues[i];
406 }
407 return min;
408}
409/*}}}*/
410/*FUNCTION TriaVertexForcing::MinAbs{{{1*/
411double TriaVertexForcing::MinAbs(void){
412
413 const int numnodes=3;
414 double min;
415 double time;
416 double triavalues[3];
417
418 /*First, recover current time from parameters: */
419 parameters->FindParam(&time,TimeEnum);
420
421 /*Retrieve interpolated values for this time step: */
422 this->GetTimeValues(&triavalues[0],time);
423
424 min=fabs(triavalues[0]);
425
426 for(int i=1;i<numnodes;i++){
427 if(fabs(triavalues[i])<min) min=fabs(triavalues[i]);
428 }
429 return min;
430}
431/*}}}*/
432/*FUNCTION TriaVertexForcing::GetVectorFromInputs{{{1*/
433void TriaVertexForcing::GetVectorFromInputs(Vec vector,int* doflist){
434
435 const int numvertices=3;
436 double time;
437 double triavalues[3];
438
439 /*First, recover current time from parameters: */
440 parameters->FindParam(&time,TimeEnum);
441
442 /*Retrieve interpolated values for this time step: */
443 this->GetTimeValues(&triavalues[0],time);
444
445 VecSetValues(vector,numvertices,doflist,(const double*)&triavalues[0],INSERT_VALUES);
446
447} /*}}}*/
448/*FUNCTION TriaVertexForcing::GetTimeValues{{{1*/
449void TriaVertexForcing::GetTimeValues(double* outvalues,double intime){
450
451 int i,j;
452 double deltat;
453 double alpha;
454 bool found=false;
455
456 /*Ok, we have the time, go through the timesteps, and figure out which interval we
457 *fall within. Then interpolate the values on this interval: */
458 if(intime<this->timesteps[0]){
459 /*get values for the first time: */
460 outvalues[0]=this->values[3*0+0];
461 outvalues[1]=this->values[3*0+1];
462 outvalues[2]=this->values[3*0+2];
463 found=true;
464 }
465 else if(intime>this->timesteps[this->numtimesteps-1]){
466 /*get values for the last time: */
467
468 outvalues[0]=this->values[3*this->numtimesteps-1+0];
469 outvalues[1]=this->values[3*this->numtimesteps-1+1];
470 outvalues[2]=this->values[3*this->numtimesteps-1+2];
471 found=true;
472 }
473 else{
474 /*Find which interval we fall within: */
475 for(i=0;i<this->numtimesteps;i++){
476 if(intime==this->timesteps[i]){
477 /*We are right on one step time: */
478 outvalues[0]=this->values[3*i+0];
479 outvalues[1]=this->values[3*i+1];
480 outvalues[2]=this->values[3*i+2];
481 found=true;
482 break; //we are done with the time interpolation.
483 }
484 else{
485 if(this->timesteps[i]<intime && intime<this->timesteps[i+1]){
486 /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
487 deltat=this->timesteps[i+1]-this->timesteps[i];
488 alpha=(intime-this->timesteps[i])/deltat;
489 for(j=0;j<3;j++){
490 outvalues[j]=this->values[3*i+j]+alpha*(this->values[3*(i+1)+j]-this->values[3*i+j]);
491 }
492 found=true;
493 break;
494 }
495 else continue; //keep looking on the next interval
496 }
497 }
498 }
499 if(!found)_error_("did not find time interval on which to interpolate forcing values!");
500
501}
502/*}}}*/
503/*FUNCTION TriaVertexForcing::Configure{{{1*/
504void TriaVertexForcing::Configure(Parameters* parameters){
505 this->parameters=parameters;
506}
507/*}}}*/
Note: See TracBrowser for help on using the repository browser.