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(){{{*/
|
---|
21 | PentaP1Input::PentaP1Input(){
|
---|
22 | return;
|
---|
23 | }
|
---|
24 | /*}}}*/
|
---|
25 | /*FUNCTION PentaP1Input::PentaP1Input(int in_enum_type,IssmPDouble* values){{{*/
|
---|
26 | PentaP1Input::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(){{{*/
|
---|
44 | PentaP1Input::~PentaP1Input(){
|
---|
45 | return;
|
---|
46 | }
|
---|
47 | /*}}}*/
|
---|
48 |
|
---|
49 | /*Object virtual functions definitions:*/
|
---|
50 | /*FUNCTION PentaP1Input::Echo {{{*/
|
---|
51 | void PentaP1Input::Echo(void){
|
---|
52 | this->DeepEcho();
|
---|
53 | }
|
---|
54 | /*}}}*/
|
---|
55 | /*FUNCTION PentaP1Input::DeepEcho{{{*/
|
---|
56 | void 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{{{*/
|
---|
64 | int PentaP1Input::Id(void){ return -1; }
|
---|
65 | /*}}}*/
|
---|
66 | /*FUNCTION PentaP1Input::MyRank{{{*/
|
---|
67 | int PentaP1Input::MyRank(void){
|
---|
68 | extern int my_rank;
|
---|
69 | return my_rank;
|
---|
70 | }
|
---|
71 | /*}}}*/
|
---|
72 | /*FUNCTION PentaP1Input::ObjectEnum{{{*/
|
---|
73 | int PentaP1Input::ObjectEnum(void){
|
---|
74 |
|
---|
75 | return PentaP1InputEnum;
|
---|
76 |
|
---|
77 | }
|
---|
78 | /*}}}*/
|
---|
79 |
|
---|
80 | /*PentaP1Input management*/
|
---|
81 | /*FUNCTION PentaP1Input::copy{{{*/
|
---|
82 | Object* PentaP1Input::copy() {
|
---|
83 |
|
---|
84 | return new PentaP1Input(this->enum_type,this->values);
|
---|
85 |
|
---|
86 | }
|
---|
87 | /*}}}*/
|
---|
88 | /*FUNCTION PentaP1Input::InstanceEnum{{{*/
|
---|
89 | int PentaP1Input::InstanceEnum(void){
|
---|
90 |
|
---|
91 | return this->enum_type;
|
---|
92 |
|
---|
93 | }
|
---|
94 | /*}}}*/
|
---|
95 | /*FUNCTION PentaP1Input::SpawnTriaInput{{{*/
|
---|
96 | Input* 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{{{*/
|
---|
121 | ElementResult* 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){{{*/
|
---|
130 | void 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){{{*/
|
---|
138 | void 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{{{*/
|
---|
145 | void 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{{{*/
|
---|
190 | void 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{{{*/
|
---|
235 | void 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{{{*/
|
---|
281 | void 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{{{*/
|
---|
305 | void 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{{{*/
|
---|
329 | void PentaP1Input::ChangeEnum(int newenumtype){
|
---|
330 | this->enum_type=newenumtype;
|
---|
331 | }
|
---|
332 | /*}}}*/
|
---|
333 | /*FUNCTION PentaP1Input::GetInputAverage{{{*/
|
---|
334 | void 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{{{*/
|
---|
341 | void 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{{{*/
|
---|
364 | void 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{{{*/
|
---|
373 | IssmPDouble 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{{{*/
|
---|
384 | IssmPDouble 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{{{*/
|
---|
396 | IssmPDouble 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{{{*/
|
---|
408 | IssmPDouble 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{{{*/
|
---|
420 | IssmPDouble 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{{{*/
|
---|
432 | void 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{{{*/
|
---|
441 | void 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{{{*/
|
---|
468 | void 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{{{*/
|
---|
479 | void 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{{{*/
|
---|
488 | void 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{{{*/
|
---|
518 | Input* 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{{{*/
|
---|
549 | Input* 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{{{*/
|
---|
580 | Input* 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{{{*/
|
---|
611 | void 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{{{*/
|
---|
618 | void PentaP1Input::GetValuesPtr(IssmPDouble** pvalues,int* pnum_values){
|
---|
619 |
|
---|
620 | *pvalues=this->values;
|
---|
621 | *pnum_values=6;
|
---|
622 |
|
---|
623 | }
|
---|
624 | /*}}}*/
|
---|
625 | /*FUNCTION PentaP1Input::Configure{{{*/
|
---|
626 | void PentaP1Input::Configure(Parameters* parameters){
|
---|
627 | /*do nothing: */
|
---|
628 | }
|
---|
629 | /*}}}*/
|
---|