source: issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp@ 5660

Last change on this file since 5660 was 5660, checked in by Mathieu Morlighem, 15 years ago

Deleted ComputePressure
Deleted GetParameterValues, now use GetParameterOnVertices more robust
File size: 20.5 KB
Line 
1/*!\file PentaVertexInput.c
2 * \brief: implementation of the PentaVertexInput 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/*PentaVertexInput constructors and destructor*/
20/*FUNCTION PentaVertexInput::PentaVertexInput(){{{1*/
21PentaVertexInput::PentaVertexInput(){
22 return;
23}
24/*}}}*/
25/*FUNCTION PentaVertexInput::PentaVertexInput(int in_enum_type,double* values){{{1*/
26PentaVertexInput::PentaVertexInput(int in_enum_type,double* 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 PentaVertexInput::~PentaVertexInput(){{{1*/
44PentaVertexInput::~PentaVertexInput(){
45 return;
46}
47/*}}}*/
48
49/*Object virtual functions definitions:*/
50/*FUNCTION PentaVertexInput::Echo {{{1*/
51void PentaVertexInput::Echo(void){
52 this->DeepEcho();
53}
54/*}}}*/
55/*FUNCTION PentaVertexInput::DeepEcho{{{1*/
56void PentaVertexInput::DeepEcho(void){
57
58 printf("PentaVertexInput:\n");
59 printf(" enum: %i (%s)\n",this->enum_type,EnumToString(this->enum_type));
60 printf(" values: [%g %g %g %g %g %g]\n",this->values[0],this->values[1],this->values[2],this->values[3],this->values[4],this->values[5]);
61}
62/*}}}*/
63/*FUNCTION PentaVertexInput::Id{{{1*/
64int PentaVertexInput::Id(void){ return -1; }
65/*}}}*/
66/*FUNCTION PentaVertexInput::MyRank{{{1*/
67int PentaVertexInput::MyRank(void){
68 extern int my_rank;
69 return my_rank;
70}
71/*}}}*/
72/*FUNCTION PentaVertexInput::Marshall{{{1*/
73void PentaVertexInput::Marshall(char** pmarshalled_dataset){
74
75 char* marshalled_dataset=NULL;
76 int enum_value=0;
77
78 /*recover marshalled_dataset: */
79 marshalled_dataset=*pmarshalled_dataset;
80
81 /*get enum value of PentaVertexInput: */
82 enum_value=PentaVertexInputEnum;
83
84 /*marshall enum: */
85 memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
86
87 /*marshall PentaVertexInput data: */
88 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
89 memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values);
90
91 *pmarshalled_dataset=marshalled_dataset;
92}
93/*}}}*/
94/*FUNCTION PentaVertexInput::MarshallSize{{{1*/
95int PentaVertexInput::MarshallSize(){
96
97 return sizeof(values)+
98 +sizeof(enum_type)+
99 +sizeof(int); //sizeof(int) for enum value
100}
101/*}}}*/
102/*FUNCTION PentaVertexInput::Demarshall{{{1*/
103void PentaVertexInput::Demarshall(char** pmarshalled_dataset){
104
105 char* marshalled_dataset=NULL;
106 int i;
107
108 /*recover marshalled_dataset: */
109 marshalled_dataset=*pmarshalled_dataset;
110
111 /*this time, no need to get enum type, the pointer directly points to the beginning of the
112 *object data (thanks to DataSet::Demarshall):*/
113 memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
114 memcpy(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values);
115
116 /*return: */
117 *pmarshalled_dataset=marshalled_dataset;
118 return;
119}
120/*}}}*/
121/*FUNCTION PentaVertexInput::Enum{{{1*/
122int PentaVertexInput::Enum(void){
123
124 return PentaVertexInputEnum;
125
126}
127/*}}}*/
128
129/*PentaVertexInput management*/
130/*FUNCTION PentaVertexInput::copy{{{1*/
131Object* PentaVertexInput::copy() {
132
133 return new PentaVertexInput(this->enum_type,this->values);
134
135}
136/*}}}*/
137/*FUNCTION PentaVertexInput::EnumType{{{1*/
138int PentaVertexInput::EnumType(void){
139
140 return this->enum_type;
141
142}
143/*}}}*/
144/*FUNCTION PentaVertexInput::SpawnTriaInput{{{1*/
145Input* PentaVertexInput::SpawnTriaInput(int* indices){
146
147 /*output*/
148 TriaVertexInput* outinput=NULL;
149 double newvalues[3];
150
151 /*Loop over the new indices*/
152 for(int i=0;i<3;i++){
153
154 /*Check index value*/
155 ISSMASSERT(indices[i]>=0 && indices[i]<6);
156
157 /*Assign value to new input*/
158 newvalues[i]=this->values[indices[i]];
159 }
160
161 /*Create new Tria input*/
162 outinput=new TriaVertexInput(this->enum_type,&newvalues[0]);
163
164 /*Assign output*/
165 return outinput;
166
167}
168/*}}}*/
169/*FUNCTION PentaVertexInput::SpawnResult{{{1*/
170ElementResult* PentaVertexInput::SpawnResult(int step, double time){
171
172 return new PentaVertexElementResult(this->enum_type,this->values,step,time);
173
174}
175/*}}}*/
176
177/*Object functions*/
178/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
179void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
180
181 /*Call PentaRef function*/
182 PentaRef::GetParameterValue(pvalue,&values[0],gauss);
183
184}
185/*}}}*/
186/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
187void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){
188
189 /*Call PentaRef function*/
190 PentaRef::GetParameterValue(pvalue,&values[0],gauss);
191
192}
193/*}}}*/
194/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/
195void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
196
197 /*Call PentaRef function*/
198 PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
199}
200/*}}}*/
201/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){{{1*/
202void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){
203
204 /*Call PentaRef function*/
205 PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
206}
207/*}}}*/
208/*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
209void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
210 int i,j;
211
212 const int numgrids=6;
213 const int DOFVELOCITY=3;
214 double B[8][27];
215 double B_reduced[6][DOFVELOCITY*numgrids];
216 double velocity[numgrids][DOFVELOCITY];
217
218 /*Get B matrix: */
219 GetBStokes(&B[0][0], xyz_list, gauss);
220 /*Create a reduced matrix of B to get rid of pressure */
221 for (i=0;i<6;i++){
222 for (j=0;j<3;j++){
223 B_reduced[i][j]=B[i][j];
224 }
225 for (j=4;j<7;j++){
226 B_reduced[i][j-1]=B[i][j];
227 }
228 for (j=8;j<11;j++){
229 B_reduced[i][j-2]=B[i][j];
230 }
231 for (j=12;j<15;j++){
232 B_reduced[i][j-3]=B[i][j];
233 }
234 for (j=16;j<19;j++){
235 B_reduced[i][j-4]=B[i][j];
236 }
237 for (j=20;j<23;j++){
238 B_reduced[i][j-5]=B[i][j];
239 }
240 }
241
242 /*Here, we are computing the strain rate of (vx,0,0)*/
243 for(i=0;i<numgrids;i++){
244 velocity[i][0]=this->values[i];
245 velocity[i][1]=0.0;
246 velocity[i][2]=0.0;
247 }
248 /*Multiply B by velocity, to get strain rate: */
249 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
250
251}
252/*}}}*/
253/*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
254void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
255 int i,j;
256
257 const int numgrids=6;
258 const int DOFVELOCITY=3;
259 double B[8][27];
260 double B_reduced[6][DOFVELOCITY*numgrids];
261 double velocity[numgrids][DOFVELOCITY];
262
263 /*Get B matrix: */
264 GetBStokes(&B[0][0], xyz_list, gauss);
265 /*Create a reduced matrix of B to get rid of pressure */
266 for (i=0;i<6;i++){
267 for (j=0;j<3;j++){
268 B_reduced[i][j]=B[i][j];
269 }
270 for (j=4;j<7;j++){
271 B_reduced[i][j-1]=B[i][j];
272 }
273 for (j=8;j<11;j++){
274 B_reduced[i][j-2]=B[i][j];
275 }
276 for (j=12;j<15;j++){
277 B_reduced[i][j-3]=B[i][j];
278 }
279 for (j=16;j<19;j++){
280 B_reduced[i][j-4]=B[i][j];
281 }
282 for (j=20;j<23;j++){
283 B_reduced[i][j-5]=B[i][j];
284 }
285 }
286
287 /*Here, we are computing the strain rate of (0,vy,0)*/
288 for(i=0;i<numgrids;i++){
289 velocity[i][0]=0.0;
290 velocity[i][1]=this->values[i];
291 velocity[i][2]=0.0;
292 }
293 /*Multiply B by velocity, to get strain rate: */
294 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
295
296}
297/*}}}*/
298/*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
299void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
300 int i,j;
301
302 const int numgrids=6;
303 const int DOFVELOCITY=3;
304 double B[8][27];
305 double B_reduced[6][DOFVELOCITY*numgrids];
306 double velocity[numgrids][DOFVELOCITY];
307
308 /*Get B matrix: */
309 GetBStokes(&B[0][0], xyz_list, gauss);
310 /*Create a reduced matrix of B to get rid of pressure */
311 for (i=0;i<6;i++){
312 for (j=0;j<3;j++){
313 B_reduced[i][j]=B[i][j];
314 }
315 for (j=4;j<7;j++){
316 B_reduced[i][j-1]=B[i][j];
317 }
318 for (j=8;j<11;j++){
319 B_reduced[i][j-2]=B[i][j];
320 }
321 for (j=12;j<15;j++){
322 B_reduced[i][j-3]=B[i][j];
323 }
324 for (j=16;j<19;j++){
325 B_reduced[i][j-4]=B[i][j];
326 }
327 for (j=20;j<23;j++){
328 B_reduced[i][j-5]=B[i][j];
329 }
330 }
331
332 /*Here, we are computing the strain rate of (0,0,vz)*/
333 for(i=0;i<numgrids;i++){
334 velocity[i][0]=0.0;
335 velocity[i][1]=0.0;
336 velocity[i][2]=this->values[i];
337 }
338
339 /*Multiply B by velocity, to get strain rate: */
340 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
341
342}
343/*}}}*/
344/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
345void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){
346
347 int i;
348 const int numgrids=6;
349 const int NDOF2=2;
350 double B[5][NDOF2*numgrids];
351 double velocity[numgrids][NDOF2];
352
353 /*Get B matrix: */
354 GetBPattyn(&B[0][0], xyz_list, gauss);
355
356 /*Here, we are computing the strain rate of (vx,0)*/
357 for(i=0;i<numgrids;i++){
358 velocity[i][0]=this->values[i];
359 velocity[i][1]=0.0;
360 }
361
362 /*Multiply B by velocity, to get strain rate: */
363 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
364 &velocity[0][0],NDOF2*numgrids,1,0,
365 epsilonvx,0);
366
367}
368/*}}}*/
369/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
370void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){
371
372 int i;
373 const int numgrids=6;
374 const int NDOF2=2;
375 double B[5][NDOF2*numgrids];
376 double velocity[numgrids][NDOF2];
377
378 /*Get B matrix: */
379 GetBPattyn(&B[0][0], xyz_list, gauss);
380
381 /*Here, we are computing the strain rate of (0,vy)*/
382 for(i=0;i<numgrids;i++){
383 velocity[i][0]=0.0;
384 velocity[i][1]=this->values[i];
385 }
386
387 /*Multiply B by velocity, to get strain rate: */
388 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
389 &velocity[0][0],NDOF2*numgrids,1,0,
390 epsilonvy,0);
391
392}
393/*}}}*/
394/*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
395void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){
396 int i,j;
397
398 const int numgrids=6;
399 const int DOFVELOCITY=3;
400 double B[8][27];
401 double B_reduced[6][DOFVELOCITY*numgrids];
402 double velocity[numgrids][DOFVELOCITY];
403
404 /*Get B matrix: */
405 GetBStokes(&B[0][0], xyz_list, gauss);
406 /*Create a reduced matrix of B to get rid of pressure */
407 for (i=0;i<6;i++){
408 for (j=0;j<3;j++){
409 B_reduced[i][j]=B[i][j];
410 }
411 for (j=4;j<7;j++){
412 B_reduced[i][j-1]=B[i][j];
413 }
414 for (j=8;j<11;j++){
415 B_reduced[i][j-2]=B[i][j];
416 }
417 for (j=12;j<15;j++){
418 B_reduced[i][j-3]=B[i][j];
419 }
420 for (j=16;j<19;j++){
421 B_reduced[i][j-4]=B[i][j];
422 }
423 for (j=20;j<23;j++){
424 B_reduced[i][j-5]=B[i][j];
425 }
426 }
427
428 /*Here, we are computing the strain rate of (vx,0,0)*/
429 for(i=0;i<numgrids;i++){
430 velocity[i][0]=this->values[i];
431 velocity[i][1]=0.0;
432 velocity[i][2]=0.0;
433 }
434 /*Multiply B by velocity, to get strain rate: */
435 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
436
437}
438/*}}}*/
439/*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
440void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){
441 int i,j;
442
443 const int numgrids=6;
444 const int DOFVELOCITY=3;
445 double B[8][27];
446 double B_reduced[6][DOFVELOCITY*numgrids];
447 double velocity[numgrids][DOFVELOCITY];
448
449 /*Get B matrix: */
450 GetBStokes(&B[0][0], xyz_list, gauss);
451 /*Create a reduced matrix of B to get rid of pressure */
452 for (i=0;i<6;i++){
453 for (j=0;j<3;j++){
454 B_reduced[i][j]=B[i][j];
455 }
456 for (j=4;j<7;j++){
457 B_reduced[i][j-1]=B[i][j];
458 }
459 for (j=8;j<11;j++){
460 B_reduced[i][j-2]=B[i][j];
461 }
462 for (j=12;j<15;j++){
463 B_reduced[i][j-3]=B[i][j];
464 }
465 for (j=16;j<19;j++){
466 B_reduced[i][j-4]=B[i][j];
467 }
468 for (j=20;j<23;j++){
469 B_reduced[i][j-5]=B[i][j];
470 }
471 }
472
473 /*Here, we are computing the strain rate of (0,vy,0)*/
474 for(i=0;i<numgrids;i++){
475 velocity[i][0]=0.0;
476 velocity[i][1]=this->values[i];
477 velocity[i][2]=0.0;
478 }
479 /*Multiply B by velocity, to get strain rate: */
480 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
481
482}
483/*}}}*/
484/*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
485void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){
486 int i,j;
487
488 const int numgrids=6;
489 const int DOFVELOCITY=3;
490 double B[8][27];
491 double B_reduced[6][DOFVELOCITY*numgrids];
492 double velocity[numgrids][DOFVELOCITY];
493
494 /*Get B matrix: */
495 GetBStokes(&B[0][0], xyz_list, gauss);
496 /*Create a reduced matrix of B to get rid of pressure */
497 for (i=0;i<6;i++){
498 for (j=0;j<3;j++){
499 B_reduced[i][j]=B[i][j];
500 }
501 for (j=4;j<7;j++){
502 B_reduced[i][j-1]=B[i][j];
503 }
504 for (j=8;j<11;j++){
505 B_reduced[i][j-2]=B[i][j];
506 }
507 for (j=12;j<15;j++){
508 B_reduced[i][j-3]=B[i][j];
509 }
510 for (j=16;j<19;j++){
511 B_reduced[i][j-4]=B[i][j];
512 }
513 for (j=20;j<23;j++){
514 B_reduced[i][j-5]=B[i][j];
515 }
516 }
517
518 /*Here, we are computing the strain rate of (0,0,vz)*/
519 for(i=0;i<numgrids;i++){
520 velocity[i][0]=0.0;
521 velocity[i][1]=0.0;
522 velocity[i][2]=this->values[i];
523 }
524
525 /*Multiply B by velocity, to get strain rate: */
526 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
527
528}
529/*}}}*/
530/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
531void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){
532
533 int i;
534 const int numgrids=6;
535 const int NDOF2=2;
536 double B[5][NDOF2*numgrids];
537 double velocity[numgrids][NDOF2];
538
539 /*Get B matrix: */
540 GetBPattyn(&B[0][0], xyz_list, gauss);
541
542 /*Here, we are computing the strain rate of (vx,0)*/
543 for(i=0;i<numgrids;i++){
544 velocity[i][0]=this->values[i];
545 velocity[i][1]=0.0;
546 }
547
548 /*Multiply B by velocity, to get strain rate: */
549 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
550 &velocity[0][0],NDOF2*numgrids,1,0,
551 epsilonvx,0);
552
553}
554/*}}}*/
555/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
556void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){
557
558 int i;
559 const int numgrids=6;
560 const int NDOF2=2;
561 double B[5][NDOF2*numgrids];
562 double velocity[numgrids][NDOF2];
563
564 /*Get B matrix: */
565 GetBPattyn(&B[0][0], xyz_list, gauss);
566
567 /*Here, we are computing the strain rate of (0,vy)*/
568 for(i=0;i<numgrids;i++){
569 velocity[i][0]=0.0;
570 velocity[i][1]=this->values[i];
571 }
572
573 /*Multiply B by velocity, to get strain rate: */
574 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
575 &velocity[0][0],NDOF2*numgrids,1,0,
576 epsilonvy,0);
577
578}
579/*}}}*/
580/*FUNCTION PentaVertexInput::ChangeEnum{{{1*/
581void PentaVertexInput::ChangeEnum(int newenumtype){
582 this->enum_type=newenumtype;
583}
584/*}}}*/
585/*FUNCTION PentaVertexInput::GetParameterAverage{{{1*/
586void PentaVertexInput::GetParameterAverage(double* pvalue){
587 *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
588}
589/*}}}*/
590
591/*Intermediary*/
592/*FUNCTION PentaVertexInput::SquareMin{{{1*/
593void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
594
595 int i;
596 const int numnodes=6;
597 double valuescopy[numnodes];
598 double squaremin;
599
600 /*First, copy values, to process units if requested: */
601 for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];
602
603 /*Process units if requested: */
604 if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type,parameters);
605
606 /*Now, figure out minimum of valuescopy: */
607 squaremin=pow(valuescopy[0],2);
608 for(i=1;i<numnodes;i++){
609 if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
610 }
611 /*Assign output pointers:*/
612 *psquaremin=squaremin;
613}
614/*}}}*/
615/*FUNCTION PentaVertexInput::ConstrainMin{{{1*/
616void PentaVertexInput::ConstrainMin(double minimum){
617
618 int i;
619 const int numgrids=6;
620
621 for(i=0;i<numgrids;i++) if (values[i]<minimum) values[i]=minimum;
622}
623/*}}}*/
624/*FUNCTION PentaVertexInput::InfinityNorm{{{1*/
625double PentaVertexInput::InfinityNorm(void){
626
627 /*Output*/
628 const int numgrids=6;
629 double norm=0;
630
631 for(int i=0;i<numgrids;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
632 return norm;
633}
634/*}}}*/
635/*FUNCTION PentaVertexInput::Max{{{1*/
636double PentaVertexInput::Max(void){
637
638 const int numgrids=6;
639 double max=values[0];
640
641 for(int i=1;i<numgrids;i++){
642 if(values[i]>max) max=values[i];
643 }
644 return max;
645}
646/*}}}*/
647/*FUNCTION PentaVertexInput::MaxAbs{{{1*/
648double PentaVertexInput::MaxAbs(void){
649
650 const int numgrids=6;
651 double max=fabs(values[0]);
652
653 for(int i=1;i<numgrids;i++){
654 if(fabs(values[i])>max) max=fabs(values[i]);
655 }
656 return max;
657}
658/*}}}*/
659/*FUNCTION PentaVertexInput::Min{{{1*/
660double PentaVertexInput::Min(void){
661
662 const int numgrids=6;
663 double min=values[0];
664
665 for(int i=1;i<numgrids;i++){
666 if(values[i]<min) min=values[i];
667 }
668 return min;
669}
670/*}}}*/
671/*FUNCTION PentaVertexInput::MinAbs{{{1*/
672double PentaVertexInput::MinAbs(void){
673
674 const int numgrids=6;
675 double min=fabs(values[0]);
676
677 for(int i=1;i<numgrids;i++){
678 if(fabs(values[i])<min) min=fabs(values[i]);
679 }
680 return min;
681}
682/*}}}*/
683/*FUNCTION PentaVertexInput::Scale{{{1*/
684void PentaVertexInput::Scale(double scale_factor){
685
686 int i;
687 const int numgrids=6;
688
689 for(i=0;i<numgrids;i++)values[i]=values[i]*scale_factor;
690}
691/*}}}*/
692/*FUNCTION PentaVertexInput::AXPY{{{1*/
693void PentaVertexInput::AXPY(Input* xinput,double scalar){
694
695 int i;
696 const int numgrids=6;
697 PentaVertexInput* xpentavertexinput=NULL;
698
699 /*xinput is of the same type, so cast it: */
700 xpentavertexinput=(PentaVertexInput*)xinput;
701
702 /*Carry out the AXPY operation depending on type:*/
703 switch(xinput->Enum()){
704
705 case PentaVertexInputEnum:
706 for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
707 return;
708
709 default:
710 ISSMERROR("not implemented yet");
711 }
712
713}
714/*}}}*/
715/*FUNCTION PentaVertexInput::Constrain{{{1*/
716void PentaVertexInput::Constrain(double cm_min, double cm_max){
717
718 int i;
719 const int numgrids=6;
720
721 if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
722 if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
723
724}
725/*}}}*/
726/*FUNCTION PentaVertexInput::Extrude{{{1*/
727void PentaVertexInput::Extrude(void){
728
729 int i;
730
731 /*First 3 values copied on 3 last values*/
732 for(i=0;i<3;i++) this->values[3+i]=this->values[i];
733
734}
735/*}}}*/
736/*FUNCTION PentaVertexInput::VerticallyIntegrate{{{1*/
737void PentaVertexInput::VerticallyIntegrate(Input* thickness_input){
738
739 /*Intermediaries*/
740 int i;
741 const int numgrids = 6;
742 int num_thickness_values;
743 double *thickness_values = NULL;
744
745 /*Check that input provided is a thickness*/
746 if (thickness_input->EnumType()!=ThicknessEnum) ISSMERROR("Input provided is not a Thickness (enum_type is %s)",EnumToString(thickness_input->EnumType()));
747
748 /*Get Thickness value pointer*/
749 thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);
750
751 /*vertically integrate depending on type:*/
752 switch(thickness_input->Enum()){
753
754 case PentaVertexInputEnum:
755 for(i=0;i<3;i++){
756 this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
757 this->values[i+3]=this->values[i];
758 }
759 return;
760
761 default:
762 ISSMERROR("not implemented yet");
763 }
764}
765/*}}}*/
766/*FUNCTION PentaVertexInput::PointwiseDivide{{{1*/
767Input* PentaVertexInput::PointwiseDivide(Input* inputB){
768
769 /*Ouput*/
770 PentaVertexInput* outinput=NULL;
771
772 /*Intermediaries*/
773 int i;
774 PentaVertexInput *xinputB = NULL;
775 int B_numvalues;
776 double *B_values = NULL;
777 const int numgrids = 6;
778 double AdotBvalues[numgrids];
779
780 /*Check that inputB is of the same type*/
781 if (inputB->Enum()!=PentaVertexInputEnum) ISSMERROR("Operation not permitted because inputB is of type %s",EnumToString(inputB->Enum()));
782 xinputB=(PentaVertexInput*)inputB;
783
784 /*Create point wise sum*/
785 for(i=0;i<numgrids;i++){
786 ISSMASSERT(xinputB->values[i]!=0);
787 AdotBvalues[i]=this->values[i]/xinputB->values[i];
788 }
789
790 /*Create new Penta vertex input (copy of current input)*/
791 outinput=new PentaVertexInput(this->enum_type,&AdotBvalues[0]);
792
793 /*Return output pointer*/
794 return outinput;
795
796}
797/*}}}*/
798/*FUNCTION PentaVertexInput::GetVectorFromInputs{{{1*/
799void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
800
801 const int numvertices=6;
802 VecSetValues(vector,numvertices,doflist,(const double*)this->values,INSERT_VALUES);
803
804} /*}}}*/
805/*FUNCTION PentaVertexInput::GetValuesPtr{{{1*/
806void PentaVertexInput::GetValuesPtr(double** pvalues,int* pnum_values){
807
808 *pvalues=this->values;
809 *pnum_values=6;
810
811}
812/*}}}*/
Note: See TracBrowser for help on using the repository browser.