source: issm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.cpp@ 15100

Last change on this file since 15100 was 15100, checked in by Eric.Larour, 12 years ago

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