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

Last change on this file since 15128 was 15128, checked in by Mathieu Morlighem, 12 years ago

DEL: removed ProcessUnits from core, moved everything to m/py
File size: 15.1 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 /*Now, figure out minimum of valuescopy: */
341 squaremin=pow(valuescopy[0],2);
342 for(i=1;i<numnodes;i++){
343 if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
344 }
345 /*Assign output pointers:*/
346 *psquaremin=squaremin;
347}
348/*}}}*/
349/*FUNCTION PentaP1Input::ConstrainMin{{{*/
350void PentaP1Input::ConstrainMin(IssmDouble minimum){
351
352 int i;
353 const int numnodes=6;
354
355 for(i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
356}
357/*}}}*/
358/*FUNCTION PentaP1Input::InfinityNorm{{{*/
359IssmDouble PentaP1Input::InfinityNorm(void){
360
361 /*Output*/
362 const int numnodes=6;
363 IssmDouble norm=0;
364
365 for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
366 return norm;
367}
368/*}}}*/
369/*FUNCTION PentaP1Input::Max{{{*/
370IssmDouble PentaP1Input::Max(void){
371
372 const int numnodes=6;
373 IssmDouble max=values[0];
374
375 for(int i=1;i<numnodes;i++){
376 if(values[i]>max) max=values[i];
377 }
378 return max;
379}
380/*}}}*/
381/*FUNCTION PentaP1Input::MaxAbs{{{*/
382IssmDouble PentaP1Input::MaxAbs(void){
383
384 const int numnodes=6;
385 IssmDouble max=fabs(values[0]);
386
387 for(int i=1;i<numnodes;i++){
388 if(fabs(values[i])>max) max=fabs(values[i]);
389 }
390 return max;
391}
392/*}}}*/
393/*FUNCTION PentaP1Input::Min{{{*/
394IssmDouble PentaP1Input::Min(void){
395
396 const int numnodes=6;
397 IssmDouble min=values[0];
398
399 for(int i=1;i<numnodes;i++){
400 if(values[i]<min) min=values[i];
401 }
402 return min;
403}
404/*}}}*/
405/*FUNCTION PentaP1Input::MinAbs{{{*/
406IssmDouble PentaP1Input::MinAbs(void){
407
408 const int numnodes=6;
409 IssmDouble min=fabs(values[0]);
410
411 for(int i=1;i<numnodes;i++){
412 if(fabs(values[i])<min) min=fabs(values[i]);
413 }
414 return min;
415}
416/*}}}*/
417/*FUNCTION PentaP1Input::Scale{{{*/
418void PentaP1Input::Scale(IssmDouble scale_factor){
419
420 int i;
421 const int numnodes=6;
422
423 for(i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
424}
425/*}}}*/
426/*FUNCTION PentaP1Input::AXPY{{{*/
427void PentaP1Input::AXPY(Input* xinput,IssmDouble scalar){
428
429 int i;
430 const int numnodes=6;
431
432 /*xinput is of the same type, so cast it: */
433
434 /*Carry out the AXPY operation depending on type:*/
435 switch(xinput->ObjectEnum()){
436
437 case PentaP1InputEnum:{
438 PentaP1Input* cast_input=(PentaP1Input*)xinput;
439 for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
440 return;
441 case ControlInputEnum:{
442 ControlInput* cont_input=(ControlInput*)xinput;
443 if(cont_input->values->ObjectEnum()!=PentaP1InputEnum) _error_("not supported yet");
444 PentaP1Input* cast_input=(PentaP1Input*)cont_input->values;
445 for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
446 return;
447 default:
448 _error_("not implemented yet");
449 }
450
451}
452/*}}}*/
453/*FUNCTION PentaP1Input::Constrain{{{*/
454void PentaP1Input::Constrain(IssmDouble cm_min, IssmDouble cm_max){
455
456 int i;
457 const int numnodes=6;
458
459 if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
460 if(!xIsNan<IssmDouble>(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
461
462}
463/*}}}*/
464/*FUNCTION PentaP1Input::Extrude{{{*/
465void PentaP1Input::Extrude(void){
466
467 int i;
468
469 /*First 3 values copied on 3 last values*/
470 for(i=0;i<3;i++) this->values[3+i]=this->values[i];
471}
472/*}}}*/
473/*FUNCTION PentaP1Input::VerticallyIntegrate{{{*/
474void PentaP1Input::VerticallyIntegrate(Input* thickness_input){
475
476 IssmDouble thickness;
477
478 /*Check that input provided is a thickness*/
479 if (thickness_input->InstanceEnum()!=ThicknessEnum) _error_("Input provided is not a Thickness (enum_type is " << EnumToStringx(thickness_input->InstanceEnum()) << ")");
480
481 /*vertically integrate depending on type:*/
482 switch(thickness_input->ObjectEnum()){
483
484 case PentaP1InputEnum:{
485 GaussPenta *gauss=new GaussPenta();
486 for(int iv=0;iv<3;iv++){
487 gauss->GaussVertex(iv);
488 thickness_input->GetInputValue(&thickness,gauss);
489 this->values[iv]=0.5*(this->values[iv]+this->values[iv+3]) * thickness;
490 this->values[iv+3]=this->values[iv];
491 }
492 delete gauss;
493 return; }
494
495 default:
496 _error_("not implemented yet");
497 }
498}
499/*}}}*/
500/*FUNCTION PentaP1Input::PointwiseDivide{{{*/
501Input* PentaP1Input::PointwiseDivide(Input* inputB){
502
503 /*Ouput*/
504 PentaP1Input* outinput=NULL;
505
506 /*Intermediaries*/
507 PentaP1Input *xinputB = NULL;
508 const int numnodes = 6;
509 IssmDouble AdotBvalues[numnodes];
510
511 /*Check that inputB is of the same type*/
512 if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
513 xinputB=(PentaP1Input*)inputB;
514
515 /*Create point wise sum*/
516 for(int i=0;i<numnodes;i++){
517 _assert_(xinputB->values[i]!=0);
518 AdotBvalues[i]=this->values[i]/xinputB->values[i];
519 }
520
521 /*Create new Penta vertex input (copy of current input)*/
522 outinput=new PentaP1Input(this->enum_type,&AdotBvalues[0]);
523
524 /*Return output pointer*/
525 return outinput;
526
527}
528/*}}}*/
529/*FUNCTION PentaP1Input::PointwiseMin{{{*/
530Input* PentaP1Input::PointwiseMin(Input* inputB){
531
532 /*Ouput*/
533 PentaP1Input* outinput=NULL;
534
535 /*Intermediaries*/
536 int i;
537 PentaP1Input *xinputB = NULL;
538 const int numnodes = 6;
539 IssmDouble minvalues[numnodes];
540
541 /*Check that inputB is of the same type*/
542 if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
543 xinputB=(PentaP1Input*)inputB;
544
545 /*Create point wise min*/
546 for(i=0;i<numnodes;i++){
547 if(this->values[i] > xinputB->values[i]) minvalues[i]=xinputB->values[i];
548 else minvalues[i]=this->values[i];
549 }
550
551 /*Create new Penta vertex input (copy of current input)*/
552 outinput=new PentaP1Input(this->enum_type,&minvalues[0]);
553
554 /*Return output pointer*/
555 return outinput;
556
557}
558/*}}}*/
559/*FUNCTION PentaP1Input::PointwiseMax{{{*/
560Input* PentaP1Input::PointwiseMax(Input* inputB){
561
562 /*Ouput*/
563 PentaP1Input* outinput=NULL;
564
565 /*Intermediaries*/
566 int i;
567 PentaP1Input *xinputB = NULL;
568 const int numnodes = 6;
569 IssmDouble maxvalues[numnodes];
570
571 /*Check that inputB is of the same type*/
572 if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
573 xinputB=(PentaP1Input*)inputB;
574
575 /*Create point wise max*/
576 for(i=0;i<numnodes;i++){
577 if(this->values[i] < xinputB->values[i]) maxvalues[i]=xinputB->values[i];
578 else maxvalues[i]=this->values[i];
579 }
580
581 /*Create new Penta vertex input (copy of current input)*/
582 outinput=new PentaP1Input(this->enum_type,&maxvalues[0]);
583
584 /*Return output pointer*/
585 return outinput;
586
587}
588/*}}}*/
589/*FUNCTION PentaP1Input::GetVectorFromInputs{{{*/
590void PentaP1Input::GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist){
591
592 const int numvertices=6;
593 vector->SetValues(numvertices,doflist,this->values,INS_VAL);
594
595} /*}}}*/
596/*FUNCTION PentaP1Input::Configure{{{*/
597void PentaP1Input::Configure(Parameters* parameters){
598 /*do nothing: */
599}
600/*}}}*/
Note: See TracBrowser for help on using the repository browser.