source: issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.cpp@ 12511

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

replaced all printf to _printString and _printLine
File size: 15.9 KB
Line 
1/*!\file PentaP1Input.c
2 * \brief: implementation of the PentaP1Input 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/*PentaP1Input constructors and destructor*/
20/*FUNCTION PentaP1Input::PentaP1Input(){{{*/
21PentaP1Input::PentaP1Input(){
22 return;
23}
24/*}}}*/
25/*FUNCTION PentaP1Input::PentaP1Input(int in_enum_type,IssmPDouble* values){{{*/
26PentaP1Input::PentaP1Input(int in_enum_type,IssmPDouble* in_values)
27 :PentaRef(1)
28{
29
30 /*Set PentaRef*/
31 this->SetElementType(P1Enum,0);
32 this->element_type=P1Enum;
33
34 enum_type=in_enum_type;
35 values[0]=in_values[0];
36 values[1]=in_values[1];
37 values[2]=in_values[2];
38 values[3]=in_values[3];
39 values[4]=in_values[4];
40 values[5]=in_values[5];
41}
42/*}}}*/
43/*FUNCTION PentaP1Input::~PentaP1Input(){{{*/
44PentaP1Input::~PentaP1Input(){
45 return;
46}
47/*}}}*/
48
49/*Object virtual functions definitions:*/
50/*FUNCTION PentaP1Input::Echo {{{*/
51void PentaP1Input::Echo(void){
52 this->DeepEcho();
53}
54/*}}}*/
55/*FUNCTION PentaP1Input::DeepEcho{{{*/
56void PentaP1Input::DeepEcho(void){
57
58 _printLine_("PentaP1Input:");
59 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")");
60 _printLine_(" values: [" << this->values[0] << " " << this->values[1] << " " << this->values[2] << " " << this->values[3] << " " << this->values[4] << " " << this->values[5] << "]");
61}
62/*}}}*/
63/*FUNCTION PentaP1Input::Id{{{*/
64int PentaP1Input::Id(void){ return -1; }
65/*}}}*/
66/*FUNCTION PentaP1Input::MyRank{{{*/
67int PentaP1Input::MyRank(void){
68 extern int my_rank;
69 return my_rank;
70}
71/*}}}*/
72/*FUNCTION PentaP1Input::ObjectEnum{{{*/
73int PentaP1Input::ObjectEnum(void){
74
75 return PentaP1InputEnum;
76
77}
78/*}}}*/
79
80/*PentaP1Input management*/
81/*FUNCTION PentaP1Input::copy{{{*/
82Object* PentaP1Input::copy() {
83
84 return new PentaP1Input(this->enum_type,this->values);
85
86}
87/*}}}*/
88/*FUNCTION PentaP1Input::InstanceEnum{{{*/
89int PentaP1Input::InstanceEnum(void){
90
91 return this->enum_type;
92
93}
94/*}}}*/
95/*FUNCTION PentaP1Input::SpawnTriaInput{{{*/
96Input* PentaP1Input::SpawnTriaInput(int* indices){
97
98 /*output*/
99 TriaP1Input* outinput=NULL;
100 IssmPDouble newvalues[3];
101
102 /*Loop over the new indices*/
103 for(int i=0;i<3;i++){
104
105 /*Check index value*/
106 _assert_(indices[i]>=0 && indices[i]<6);
107
108 /*Assign value to new input*/
109 newvalues[i]=this->values[indices[i]];
110 }
111
112 /*Create new Tria input*/
113 outinput=new TriaP1Input(this->enum_type,&newvalues[0]);
114
115 /*Assign output*/
116 return outinput;
117
118}
119/*}}}*/
120/*FUNCTION PentaP1Input::SpawnResult{{{*/
121ElementResult* PentaP1Input::SpawnResult(int step, IssmPDouble time){
122
123 return new PentaP1ElementResult(this->enum_type,this->values,step,time);
124
125}
126/*}}}*/
127
128/*Object functions*/
129/*FUNCTION PentaP1Input::GetInputValue(IssmPDouble* pvalue,GaussPenta* gauss){{{*/
130void PentaP1Input::GetInputValue(IssmPDouble* pvalue,GaussPenta* gauss){
131
132 /*Call PentaRef function*/
133 PentaRef::GetInputValue(pvalue,&values[0],gauss);
134
135}
136/*}}}*/
137/*FUNCTION PentaP1Input::GetInputDerivativeValue(IssmPDouble* p, IssmPDouble* xyz_list, GaussPenta* gauss){{{*/
138void PentaP1Input::GetInputDerivativeValue(IssmPDouble* p, IssmPDouble* xyz_list, GaussPenta* gauss){
139
140 /*Call PentaRef function*/
141 PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
142}
143/*}}}*/
144/*FUNCTION PentaP1Input::GetVxStrainRate3d{{{*/
145void PentaP1Input::GetVxStrainRate3d(IssmPDouble* epsilonvx,IssmPDouble* xyz_list, GaussPenta* gauss){
146 int i,j;
147
148 const int numnodes=6;
149 const int DOFVELOCITY=3;
150 IssmPDouble B[8][27];
151 IssmPDouble B_reduced[6][DOFVELOCITY*numnodes];
152 IssmPDouble velocity[numnodes][DOFVELOCITY];
153
154 /*Get B matrix: */
155 GetBStokes(&B[0][0], xyz_list, gauss);
156 /*Create a reduced matrix of B to get rid of pressure */
157 for (i=0;i<6;i++){
158 for (j=0;j<3;j++){
159 B_reduced[i][j]=B[i][j];
160 }
161 for (j=4;j<7;j++){
162 B_reduced[i][j-1]=B[i][j];
163 }
164 for (j=8;j<11;j++){
165 B_reduced[i][j-2]=B[i][j];
166 }
167 for (j=12;j<15;j++){
168 B_reduced[i][j-3]=B[i][j];
169 }
170 for (j=16;j<19;j++){
171 B_reduced[i][j-4]=B[i][j];
172 }
173 for (j=20;j<23;j++){
174 B_reduced[i][j-5]=B[i][j];
175 }
176 }
177
178 /*Here, we are computing the strain rate of (vx,0,0)*/
179 for(i=0;i<numnodes;i++){
180 velocity[i][0]=this->values[i];
181 velocity[i][1]=0.0;
182 velocity[i][2]=0.0;
183 }
184 /*Multiply B by velocity, to get strain rate: */
185 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvx,0);
186
187}
188/*}}}*/
189/*FUNCTION PentaP1Input::GetVyStrainRate3d{{{*/
190void PentaP1Input::GetVyStrainRate3d(IssmPDouble* epsilonvy,IssmPDouble* xyz_list, GaussPenta* gauss){
191 int i,j;
192
193 const int numnodes=6;
194 const int DOFVELOCITY=3;
195 IssmPDouble B[8][27];
196 IssmPDouble B_reduced[6][DOFVELOCITY*numnodes];
197 IssmPDouble velocity[numnodes][DOFVELOCITY];
198
199 /*Get B matrix: */
200 GetBStokes(&B[0][0], xyz_list, gauss);
201 /*Create a reduced matrix of B to get rid of pressure */
202 for (i=0;i<6;i++){
203 for (j=0;j<3;j++){
204 B_reduced[i][j]=B[i][j];
205 }
206 for (j=4;j<7;j++){
207 B_reduced[i][j-1]=B[i][j];
208 }
209 for (j=8;j<11;j++){
210 B_reduced[i][j-2]=B[i][j];
211 }
212 for (j=12;j<15;j++){
213 B_reduced[i][j-3]=B[i][j];
214 }
215 for (j=16;j<19;j++){
216 B_reduced[i][j-4]=B[i][j];
217 }
218 for (j=20;j<23;j++){
219 B_reduced[i][j-5]=B[i][j];
220 }
221 }
222
223 /*Here, we are computing the strain rate of (0,vy,0)*/
224 for(i=0;i<numnodes;i++){
225 velocity[i][0]=0.0;
226 velocity[i][1]=this->values[i];
227 velocity[i][2]=0.0;
228 }
229 /*Multiply B by velocity, to get strain rate: */
230 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvy,0);
231
232}
233/*}}}*/
234/*FUNCTION PentaP1Input::GetVzStrainRate3d{{{*/
235void PentaP1Input::GetVzStrainRate3d(IssmPDouble* epsilonvz,IssmPDouble* xyz_list, GaussPenta* gauss){
236 int i,j;
237
238 const int numnodes=6;
239 const int DOFVELOCITY=3;
240 IssmPDouble B[8][27];
241 IssmPDouble B_reduced[6][DOFVELOCITY*numnodes];
242 IssmPDouble velocity[numnodes][DOFVELOCITY];
243
244 /*Get B matrix: */
245 GetBStokes(&B[0][0], xyz_list, gauss);
246 /*Create a reduced matrix of B to get rid of pressure */
247 for (i=0;i<6;i++){
248 for (j=0;j<3;j++){
249 B_reduced[i][j]=B[i][j];
250 }
251 for (j=4;j<7;j++){
252 B_reduced[i][j-1]=B[i][j];
253 }
254 for (j=8;j<11;j++){
255 B_reduced[i][j-2]=B[i][j];
256 }
257 for (j=12;j<15;j++){
258 B_reduced[i][j-3]=B[i][j];
259 }
260 for (j=16;j<19;j++){
261 B_reduced[i][j-4]=B[i][j];
262 }
263 for (j=20;j<23;j++){
264 B_reduced[i][j-5]=B[i][j];
265 }
266 }
267
268 /*Here, we are computing the strain rate of (0,0,vz)*/
269 for(i=0;i<numnodes;i++){
270 velocity[i][0]=0.0;
271 velocity[i][1]=0.0;
272 velocity[i][2]=this->values[i];
273 }
274
275 /*Multiply B by velocity, to get strain rate: */
276 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvz,0);
277
278}
279/*}}}*/
280/*FUNCTION PentaP1Input::GetVxStrainRate3dPattyn{{{*/
281void PentaP1Input::GetVxStrainRate3dPattyn(IssmPDouble* epsilonvx,IssmPDouble* xyz_list, GaussPenta* gauss){
282
283 int i;
284 const int numnodes=6;
285 IssmPDouble B[5][NDOF2*numnodes];
286 IssmPDouble velocity[numnodes][NDOF2];
287
288 /*Get B matrix: */
289 GetBPattyn(&B[0][0], xyz_list, gauss);
290
291 /*Here, we are computing the strain rate of (vx,0)*/
292 for(i=0;i<numnodes;i++){
293 velocity[i][0]=this->values[i];
294 velocity[i][1]=0.0;
295 }
296
297 /*Multiply B by velocity, to get strain rate: */
298 MatrixMultiply( &B[0][0],5,NDOF2*numnodes,0,
299 &velocity[0][0],NDOF2*numnodes,1,0,
300 epsilonvx,0);
301
302}
303/*}}}*/
304/*FUNCTION PentaP1Input::GetVyStrainRate3dPattyn{{{*/
305void PentaP1Input::GetVyStrainRate3dPattyn(IssmPDouble* epsilonvy,IssmPDouble* xyz_list, GaussPenta* gauss){
306
307 int i;
308 const int numnodes=6;
309 IssmPDouble B[5][NDOF2*numnodes];
310 IssmPDouble velocity[numnodes][NDOF2];
311
312 /*Get B matrix: */
313 GetBPattyn(&B[0][0], xyz_list, gauss);
314
315 /*Here, we are computing the strain rate of (0,vy)*/
316 for(i=0;i<numnodes;i++){
317 velocity[i][0]=0.0;
318 velocity[i][1]=this->values[i];
319 }
320
321 /*Multiply B by velocity, to get strain rate: */
322 MatrixMultiply( &B[0][0],5,NDOF2*numnodes,0,
323 &velocity[0][0],NDOF2*numnodes,1,0,
324 epsilonvy,0);
325
326}
327/*}}}*/
328/*FUNCTION PentaP1Input::ChangeEnum{{{*/
329void PentaP1Input::ChangeEnum(int newenumtype){
330 this->enum_type=newenumtype;
331}
332/*}}}*/
333/*FUNCTION PentaP1Input::GetInputAverage{{{*/
334void PentaP1Input::GetInputAverage(IssmPDouble* pvalue){
335 *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
336}
337/*}}}*/
338
339/*Intermediary*/
340/*FUNCTION PentaP1Input::SquareMin{{{*/
341void PentaP1Input::SquareMin(IssmPDouble* psquaremin, bool process_units,Parameters* parameters){
342
343 int i;
344 const int numnodes=6;
345 IssmPDouble valuescopy[numnodes];
346 IssmPDouble squaremin;
347
348 /*First, copy values, to process units if requested: */
349 for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];
350
351 /*Process units if requested: */
352 if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type);
353
354 /*Now, figure out minimum of valuescopy: */
355 squaremin=pow(valuescopy[0],2);
356 for(i=1;i<numnodes;i++){
357 if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
358 }
359 /*Assign output pointers:*/
360 *psquaremin=squaremin;
361}
362/*}}}*/
363/*FUNCTION PentaP1Input::ConstrainMin{{{*/
364void PentaP1Input::ConstrainMin(IssmPDouble minimum){
365
366 int i;
367 const int numnodes=6;
368
369 for(i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
370}
371/*}}}*/
372/*FUNCTION PentaP1Input::InfinityNorm{{{*/
373IssmPDouble PentaP1Input::InfinityNorm(void){
374
375 /*Output*/
376 const int numnodes=6;
377 IssmPDouble norm=0;
378
379 for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
380 return norm;
381}
382/*}}}*/
383/*FUNCTION PentaP1Input::Max{{{*/
384IssmPDouble PentaP1Input::Max(void){
385
386 const int numnodes=6;
387 IssmPDouble max=values[0];
388
389 for(int i=1;i<numnodes;i++){
390 if(values[i]>max) max=values[i];
391 }
392 return max;
393}
394/*}}}*/
395/*FUNCTION PentaP1Input::MaxAbs{{{*/
396IssmPDouble PentaP1Input::MaxAbs(void){
397
398 const int numnodes=6;
399 IssmPDouble max=fabs(values[0]);
400
401 for(int i=1;i<numnodes;i++){
402 if(fabs(values[i])>max) max=fabs(values[i]);
403 }
404 return max;
405}
406/*}}}*/
407/*FUNCTION PentaP1Input::Min{{{*/
408IssmPDouble PentaP1Input::Min(void){
409
410 const int numnodes=6;
411 IssmPDouble min=values[0];
412
413 for(int i=1;i<numnodes;i++){
414 if(values[i]<min) min=values[i];
415 }
416 return min;
417}
418/*}}}*/
419/*FUNCTION PentaP1Input::MinAbs{{{*/
420IssmPDouble PentaP1Input::MinAbs(void){
421
422 const int numnodes=6;
423 IssmPDouble min=fabs(values[0]);
424
425 for(int i=1;i<numnodes;i++){
426 if(fabs(values[i])<min) min=fabs(values[i]);
427 }
428 return min;
429}
430/*}}}*/
431/*FUNCTION PentaP1Input::Scale{{{*/
432void PentaP1Input::Scale(IssmPDouble scale_factor){
433
434 int i;
435 const int numnodes=6;
436
437 for(i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
438}
439/*}}}*/
440/*FUNCTION PentaP1Input::AXPY{{{*/
441void PentaP1Input::AXPY(Input* xinput,IssmPDouble scalar){
442
443 int i;
444 const int numnodes=6;
445
446 /*xinput is of the same type, so cast it: */
447
448 /*Carry out the AXPY operation depending on type:*/
449 switch(xinput->ObjectEnum()){
450
451 case PentaP1InputEnum:{
452 PentaP1Input* cast_input=(PentaP1Input*)xinput;
453 for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
454 return;
455 case ControlInputEnum:{
456 ControlInput* cont_input=(ControlInput*)xinput;
457 if(cont_input->values->ObjectEnum()!=PentaP1InputEnum) _error2_("not supported yet");
458 PentaP1Input* cast_input=(PentaP1Input*)cont_input->values;
459 for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
460 return;
461 default:
462 _error2_("not implemented yet");
463 }
464
465}
466/*}}}*/
467/*FUNCTION PentaP1Input::Constrain{{{*/
468void PentaP1Input::Constrain(IssmPDouble cm_min, IssmPDouble cm_max){
469
470 int i;
471 const int numnodes=6;
472
473 if(!isnan(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
474 if(!isnan(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
475
476}
477/*}}}*/
478/*FUNCTION PentaP1Input::Extrude{{{*/
479void PentaP1Input::Extrude(void){
480
481 int i;
482
483 /*First 3 values copied on 3 last values*/
484 for(i=0;i<3;i++) this->values[3+i]=this->values[i];
485}
486/*}}}*/
487/*FUNCTION PentaP1Input::VerticallyIntegrate{{{*/
488void PentaP1Input::VerticallyIntegrate(Input* thickness_input){
489
490 /*Intermediaries*/
491 int i;
492 const int numnodes = 6;
493 int num_thickness_values;
494 IssmPDouble *thickness_values = NULL;
495
496 /*Check that input provided is a thickness*/
497 if (thickness_input->InstanceEnum()!=ThicknessEnum) _error2_("Input provided is not a Thickness (enum_type is " << EnumToStringx(thickness_input->InstanceEnum()) << ")");
498
499 /*Get Thickness value pointer*/
500 thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);
501
502 /*vertically integrate depending on type:*/
503 switch(thickness_input->ObjectEnum()){
504
505 case PentaP1InputEnum:
506 for(i=0;i<3;i++){
507 this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
508 this->values[i+3]=this->values[i];
509 }
510 return;
511
512 default:
513 _error2_("not implemented yet");
514 }
515}
516/*}}}*/
517/*FUNCTION PentaP1Input::PointwiseDivide{{{*/
518Input* PentaP1Input::PointwiseDivide(Input* inputB){
519
520 /*Ouput*/
521 PentaP1Input* outinput=NULL;
522
523 /*Intermediaries*/
524 int i;
525 PentaP1Input *xinputB = NULL;
526 int B_numvalues;
527 const int numnodes = 6;
528 IssmPDouble AdotBvalues[numnodes];
529
530 /*Check that inputB is of the same type*/
531 if (inputB->ObjectEnum()!=PentaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
532 xinputB=(PentaP1Input*)inputB;
533
534 /*Create point wise sum*/
535 for(i=0;i<numnodes;i++){
536 _assert_(xinputB->values[i]!=0);
537 AdotBvalues[i]=this->values[i]/xinputB->values[i];
538 }
539
540 /*Create new Penta vertex input (copy of current input)*/
541 outinput=new PentaP1Input(this->enum_type,&AdotBvalues[0]);
542
543 /*Return output pointer*/
544 return outinput;
545
546}
547/*}}}*/
548/*FUNCTION PentaP1Input::PointwiseMin{{{*/
549Input* PentaP1Input::PointwiseMin(Input* inputB){
550
551 /*Ouput*/
552 PentaP1Input* outinput=NULL;
553
554 /*Intermediaries*/
555 int i;
556 PentaP1Input *xinputB = NULL;
557 int B_numvalues;
558 const int numnodes = 6;
559 IssmPDouble minvalues[numnodes];
560
561 /*Check that inputB is of the same type*/
562 if (inputB->ObjectEnum()!=PentaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
563 xinputB=(PentaP1Input*)inputB;
564
565 /*Create point wise min*/
566 for(i=0;i<numnodes;i++){
567 if(this->values[i] > xinputB->values[i]) minvalues[i]=xinputB->values[i];
568 else minvalues[i]=this->values[i];
569 }
570
571 /*Create new Penta vertex input (copy of current input)*/
572 outinput=new PentaP1Input(this->enum_type,&minvalues[0]);
573
574 /*Return output pointer*/
575 return outinput;
576
577}
578/*}}}*/
579/*FUNCTION PentaP1Input::PointwiseMax{{{*/
580Input* PentaP1Input::PointwiseMax(Input* inputB){
581
582 /*Ouput*/
583 PentaP1Input* outinput=NULL;
584
585 /*Intermediaries*/
586 int i;
587 PentaP1Input *xinputB = NULL;
588 int B_numvalues;
589 const int numnodes = 6;
590 IssmPDouble maxvalues[numnodes];
591
592 /*Check that inputB is of the same type*/
593 if (inputB->ObjectEnum()!=PentaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
594 xinputB=(PentaP1Input*)inputB;
595
596 /*Create point wise max*/
597 for(i=0;i<numnodes;i++){
598 if(this->values[i] < xinputB->values[i]) maxvalues[i]=xinputB->values[i];
599 else maxvalues[i]=this->values[i];
600 }
601
602 /*Create new Penta vertex input (copy of current input)*/
603 outinput=new PentaP1Input(this->enum_type,&maxvalues[0]);
604
605 /*Return output pointer*/
606 return outinput;
607
608}
609/*}}}*/
610/*FUNCTION PentaP1Input::GetVectorFromInputs{{{*/
611void PentaP1Input::GetVectorFromInputs(Vector* vector,int* doflist){
612
613 const int numvertices=6;
614 vector->SetValues(numvertices,doflist,this->values,INS_VAL);
615
616} /*}}}*/
617/*FUNCTION PentaP1Input::GetValuesPtr{{{*/
618void PentaP1Input::GetValuesPtr(IssmPDouble** pvalues,int* pnum_values){
619
620 *pvalues=this->values;
621 *pnum_values=6;
622
623}
624/*}}}*/
625/*FUNCTION PentaP1Input::Configure{{{*/
626void PentaP1Input::Configure(Parameters* parameters){
627 /*do nothing: */
628}
629/*}}}*/
Note: See TracBrowser for help on using the repository browser.