source: issm/trunk/src/c/objects/Inputs/PentaVertexForcing.cpp@ 8376

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

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