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

Last change on this file since 12530 was 12530, checked in by utke, 13 years ago

in the first round the assumption was made that the inputs remain passive but in the second round after discussion with M. Morlighem all inputs were consistently made active with the exception of the xyz_list arguments
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,IssmDouble* values){{{*/
26PentaP1Input::PentaP1Input(int in_enum_type,IssmDouble* 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 IssmDouble 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, IssmDouble time){
122
123 return new PentaP1ElementResult(this->enum_type,this->values,step,time);
124
125}
126/*}}}*/
127
128/*Object functions*/
129/*FUNCTION PentaP1Input::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/
130void PentaP1Input::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){
131
132 /*Call PentaRef function*/
133 PentaRef::GetInputValue(pvalue,&values[0],gauss);
134
135}
136/*}}}*/
137/*FUNCTION PentaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmPDouble* xyz_list, GaussPenta* gauss){{{*/
138void PentaP1Input::GetInputDerivativeValue(IssmDouble* 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(IssmDouble* epsilonvx,IssmPDouble* xyz_list, GaussPenta* gauss){
146 int i,j;
147
148 const int numnodes=6;
149 const int DOFVELOCITY=3;
150 IssmDouble B[8][27];
151 IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
152 IssmDouble 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(IssmDouble* epsilonvy,IssmPDouble* xyz_list, GaussPenta* gauss){
191 int i,j;
192
193 const int numnodes=6;
194 const int DOFVELOCITY=3;
195 IssmDouble B[8][27];
196 IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
197 IssmDouble 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(IssmDouble* epsilonvz,IssmPDouble* xyz_list, GaussPenta* gauss){
236 int i,j;
237
238 const int numnodes=6;
239 const int DOFVELOCITY=3;
240 IssmDouble B[8][27];
241 IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
242 IssmDouble 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(IssmDouble* epsilonvx,IssmPDouble* xyz_list, GaussPenta* gauss){
282
283 int i;
284 const int numnodes=6;
285 IssmDouble B[5][NDOF2*numnodes];
286 IssmDouble 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(IssmDouble* epsilonvy,IssmPDouble* xyz_list, GaussPenta* gauss){
306
307 int i;
308 const int numnodes=6;
309 IssmDouble B[5][NDOF2*numnodes];
310 IssmDouble 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(IssmDouble* 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(IssmDouble* psquaremin, bool process_units,Parameters* parameters){
342
343 int i;
344 const int numnodes=6;
345 IssmDouble valuescopy[numnodes];
346 IssmDouble 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(IssmDouble 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{{{*/
373IssmDouble PentaP1Input::InfinityNorm(void){
374
375 /*Output*/
376 const int numnodes=6;
377 IssmDouble 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{{{*/
384IssmDouble PentaP1Input::Max(void){
385
386 const int numnodes=6;
387 IssmDouble 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{{{*/
396IssmDouble PentaP1Input::MaxAbs(void){
397
398 const int numnodes=6;
399 IssmDouble 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{{{*/
408IssmDouble PentaP1Input::Min(void){
409
410 const int numnodes=6;
411 IssmDouble 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{{{*/
420IssmDouble PentaP1Input::MinAbs(void){
421
422 const int numnodes=6;
423 IssmDouble 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(IssmDouble 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,IssmDouble 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(IssmDouble cm_min, IssmDouble cm_max){
469
470 int i;
471 const int numnodes=6;
472
473 if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
474 if(!xIsNan<IssmDouble>(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 IssmDouble *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 IssmDouble 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 IssmDouble 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 IssmDouble 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(IssmDouble** 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.