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

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

Added InfinityNorm methods for inputs
File size: 14.9 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{{{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::GetParameterValues{{{1*/
187void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){
188 /*It is assumed that output values has been correctly allocated*/
189
190 int i,j;
191 double gauss[4];
192
193 for (i=0;i<numgauss;i++){
194
195 /*Get current Gauss point coordinates*/
196 for (j=0;j<4;j++) gauss[j]=gauss_pointers[i*4+j];
197
198 /*Assign parameter value*/
199 GetParameterValue(&values[i],&gauss[0]);
200 }
201}
202/*}}}*/
203/*FUNCTION PentaVertexInput::GetParameterDerivativeValue{{{1*/
204void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
205
206 /*Call PentaRef function*/
207 PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
208}
209/*}}}*/
210/*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
211void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
212 int i,j;
213
214 const int numgrids=6;
215 const int DOFVELOCITY=3;
216 double B[8][27];
217 double B_reduced[6][DOFVELOCITY*numgrids];
218 double velocity[numgrids][DOFVELOCITY];
219
220 /*Get B matrix: */
221 GetBStokes(&B[0][0], xyz_list, gauss);
222 /*Create a reduced matrix of B to get rid of pressure */
223 for (i=0;i<6;i++){
224 for (j=0;j<3;j++){
225 B_reduced[i][j]=B[i][j];
226 }
227 for (j=4;j<7;j++){
228 B_reduced[i][j-1]=B[i][j];
229 }
230 for (j=8;j<11;j++){
231 B_reduced[i][j-2]=B[i][j];
232 }
233 for (j=12;j<15;j++){
234 B_reduced[i][j-3]=B[i][j];
235 }
236 for (j=16;j<19;j++){
237 B_reduced[i][j-4]=B[i][j];
238 }
239 for (j=20;j<23;j++){
240 B_reduced[i][j-5]=B[i][j];
241 }
242 }
243
244 /*Here, we are computing the strain rate of (vx,0,0)*/
245 for(i=0;i<numgrids;i++){
246 velocity[i][0]=this->values[i];
247 velocity[i][1]=0.0;
248 velocity[i][2]=0.0;
249 }
250 /*Multiply B by velocity, to get strain rate: */
251 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
252
253}
254/*}}}*/
255/*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
256void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
257 int i,j;
258
259 const int numgrids=6;
260 const int DOFVELOCITY=3;
261 double B[8][27];
262 double B_reduced[6][DOFVELOCITY*numgrids];
263 double velocity[numgrids][DOFVELOCITY];
264
265 /*Get B matrix: */
266 GetBStokes(&B[0][0], xyz_list, gauss);
267 /*Create a reduced matrix of B to get rid of pressure */
268 for (i=0;i<6;i++){
269 for (j=0;j<3;j++){
270 B_reduced[i][j]=B[i][j];
271 }
272 for (j=4;j<7;j++){
273 B_reduced[i][j-1]=B[i][j];
274 }
275 for (j=8;j<11;j++){
276 B_reduced[i][j-2]=B[i][j];
277 }
278 for (j=12;j<15;j++){
279 B_reduced[i][j-3]=B[i][j];
280 }
281 for (j=16;j<19;j++){
282 B_reduced[i][j-4]=B[i][j];
283 }
284 for (j=20;j<23;j++){
285 B_reduced[i][j-5]=B[i][j];
286 }
287 }
288
289 /*Here, we are computing the strain rate of (0,vy,0)*/
290 for(i=0;i<numgrids;i++){
291 velocity[i][0]=0.0;
292 velocity[i][1]=this->values[i];
293 velocity[i][2]=0.0;
294 }
295 /*Multiply B by velocity, to get strain rate: */
296 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
297
298}
299/*}}}*/
300/*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
301void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
302 int i,j;
303
304 const int numgrids=6;
305 const int DOFVELOCITY=3;
306 double B[8][27];
307 double B_reduced[6][DOFVELOCITY*numgrids];
308 double velocity[numgrids][DOFVELOCITY];
309
310 /*Get B matrix: */
311 GetBStokes(&B[0][0], xyz_list, gauss);
312 /*Create a reduced matrix of B to get rid of pressure */
313 for (i=0;i<6;i++){
314 for (j=0;j<3;j++){
315 B_reduced[i][j]=B[i][j];
316 }
317 for (j=4;j<7;j++){
318 B_reduced[i][j-1]=B[i][j];
319 }
320 for (j=8;j<11;j++){
321 B_reduced[i][j-2]=B[i][j];
322 }
323 for (j=12;j<15;j++){
324 B_reduced[i][j-3]=B[i][j];
325 }
326 for (j=16;j<19;j++){
327 B_reduced[i][j-4]=B[i][j];
328 }
329 for (j=20;j<23;j++){
330 B_reduced[i][j-5]=B[i][j];
331 }
332 }
333
334 /*Here, we are computing the strain rate of (0,0,vz)*/
335 for(i=0;i<numgrids;i++){
336 velocity[i][0]=0.0;
337 velocity[i][1]=0.0;
338 velocity[i][2]=this->values[i];
339 }
340
341 /*Multiply B by velocity, to get strain rate: */
342 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
343
344}
345/*}}}*/
346/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
347void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){
348
349 int i;
350 const int numgrids=6;
351 const int NDOF2=2;
352 double B[5][NDOF2*numgrids];
353 double velocity[numgrids][NDOF2];
354
355 /*Get B matrix: */
356 GetBPattyn(&B[0][0], xyz_list, gauss);
357
358 /*Here, we are computing the strain rate of (vx,0)*/
359 for(i=0;i<numgrids;i++){
360 velocity[i][0]=this->values[i];
361 velocity[i][1]=0.0;
362 }
363
364 /*Multiply B by velocity, to get strain rate: */
365 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
366 &velocity[0][0],NDOF2*numgrids,1,0,
367 epsilonvx,0);
368
369}
370/*}}}*/
371/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
372void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){
373
374 int i;
375 const int numgrids=6;
376 const int NDOF2=2;
377 double B[5][NDOF2*numgrids];
378 double velocity[numgrids][NDOF2];
379
380 /*Get B matrix: */
381 GetBPattyn(&B[0][0], xyz_list, gauss);
382
383 /*Here, we are computing the strain rate of (0,vy)*/
384 for(i=0;i<numgrids;i++){
385 velocity[i][0]=0.0;
386 velocity[i][1]=this->values[i];
387 }
388
389 /*Multiply B by velocity, to get strain rate: */
390 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
391 &velocity[0][0],NDOF2*numgrids,1,0,
392 epsilonvy,0);
393
394}
395/*}}}*/
396/*FUNCTION PentaVertexInput::ChangeEnum{{{1*/
397void PentaVertexInput::ChangeEnum(int newenumtype){
398 this->enum_type=newenumtype;
399}
400/*}}}*/
401/*FUNCTION PentaVertexInput::GetParameterAverage{{{1*/
402void PentaVertexInput::GetParameterAverage(double* pvalue){
403 *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
404}
405/*}}}*/
406
407/*Intermediary*/
408/*FUNCTION PentaVertexInput::SquareMin{{{1*/
409void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
410
411 int i;
412 const int numnodes=6;
413 double valuescopy[numnodes];
414 double squaremin;
415
416 /*First, copy values, to process units if requested: */
417 for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];
418
419 /*Process units if requested: */
420 if(process_units)NodalValuesUnitConversion(&valuescopy[0],numnodes,enum_type,parameters);
421
422 /*Now, figure out minimum of valuescopy: */
423 squaremin=pow(valuescopy[0],2);
424 for(i=1;i<numnodes;i++){
425 if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
426 }
427 /*Assign output pointers:*/
428 *psquaremin=squaremin;
429}
430/*}}}*/
431/*FUNCTION PentaVertexInput::ConstrainMin{{{1*/
432void PentaVertexInput::ConstrainMin(double minimum){
433
434 int i;
435 const int numgrids=6;
436
437 for(i=0;i<numgrids;i++) if (values[i]<minimum) values[i]=minimum;
438}
439/*}}}*/
440/*FUNCTION PentaVertexInput::InfinityNorm{{{1*/
441double PentaVertexInput::InfinityNorm(void){
442
443 /*Output*/
444 const int numgrids=6;
445 double norm=0;
446
447 for(int i=0;i<numgrids;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
448 return norm;
449}
450/*}}}*/
451/*FUNCTION PentaVertexInput::Scale{{{1*/
452void PentaVertexInput::Scale(double scale_factor){
453
454 int i;
455 const int numgrids=6;
456
457 for(i=0;i<numgrids;i++)values[i]=values[i]*scale_factor;
458}
459/*}}}*/
460/*FUNCTION PentaVertexInput::AXPY{{{1*/
461void PentaVertexInput::AXPY(Input* xinput,double scalar){
462
463 int i;
464 const int numgrids=6;
465 PentaVertexInput* xpentavertexinput=NULL;
466
467 /*xinput is of the same type, so cast it: */
468 xpentavertexinput=(PentaVertexInput*)xinput;
469
470 /*Carry out the AXPY operation depending on type:*/
471 switch(xinput->Enum()){
472
473 case PentaVertexInputEnum:
474 for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
475 return;
476
477 default:
478 ISSMERROR("not implemented yet");
479 }
480
481}
482/*}}}*/
483/*FUNCTION PentaVertexInput::Constrain{{{1*/
484void PentaVertexInput::Constrain(double cm_min, double cm_max){
485
486 int i;
487 const int numgrids=6;
488
489 if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
490 if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
491
492}
493/*}}}*/
494/*FUNCTION PentaVertexInput::Extrude{{{1*/
495void PentaVertexInput::Extrude(void){
496
497 int i;
498
499 /*First 3 values copied on 3 last values*/
500 for(i=0;i<3;i++) this->values[3+i]=this->values[i];
501
502}
503/*}}}*/
504/*FUNCTION PentaVertexInput::VerticallyIntegrate{{{1*/
505void PentaVertexInput::VerticallyIntegrate(Input* thickness_input){
506
507 /*Intermediaries*/
508 int i;
509 const int numgrids = 6;
510 int num_thickness_values;
511 double *thickness_values = NULL;
512
513 /*Check that input provided is a thickness*/
514 if (thickness_input->EnumType()!=ThicknessEnum) ISSMERROR("Input provided is not a Thickness (enum_type is %s)",EnumToString(thickness_input->EnumType()));
515
516 /*Get Thickness value pointer*/
517 thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);
518
519 /*vertically integrate depending on type:*/
520 switch(thickness_input->Enum()){
521
522 case PentaVertexInputEnum:
523 for(i=0;i<3;i++){
524 this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
525 this->values[i+3]=this->values[i];
526 }
527 return;
528
529 default:
530 ISSMERROR("not implemented yet");
531 }
532}
533/*}}}*/
534/*FUNCTION PentaVertexInput::PointwiseDivide{{{1*/
535Input* PentaVertexInput::PointwiseDivide(Input* inputB){
536
537 /*Ouput*/
538 PentaVertexInput* outinput=NULL;
539
540 /*Intermediaries*/
541 int i;
542 PentaVertexInput *xinputB = NULL;
543 int B_numvalues;
544 double *B_values = NULL;
545 const int numgrids = 6;
546 double AdotBvalues[numgrids];
547
548 /*Check that inputB is of the same type*/
549 if (inputB->Enum()!=PentaVertexInputEnum) ISSMERROR("Operation not permitted because inputB is of type %s",EnumToString(inputB->Enum()));
550 xinputB=(PentaVertexInput*)inputB;
551
552 /*Create point wise sum*/
553 for(i=0;i<numgrids;i++){
554 ISSMASSERT(xinputB->values[i]!=0);
555 AdotBvalues[i]=this->values[i]/xinputB->values[i];
556 }
557
558 /*Create new Penta vertex input (copy of current input)*/
559 outinput=new PentaVertexInput(this->enum_type,&AdotBvalues[0]);
560
561 /*Return output pointer*/
562 return outinput;
563
564}
565/*}}}*/
566/*FUNCTION PentaVertexInput::GetVectorFromInputs{{{1*/
567void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
568
569 const int numvertices=6;
570 VecSetValues(vector,numvertices,doflist,(const double*)this->values,INSERT_VALUES);
571
572} /*}}}*/
573/*FUNCTION PentaVertexInput::GetValuesPtr{{{1*/
574void PentaVertexInput::GetValuesPtr(double** pvalues,int* pnum_values){
575
576 *pvalues=this->values;
577 *pnum_values=6;
578
579}
580/*}}}*/
Note: See TracBrowser for help on using the repository browser.