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

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

New UnitConversion routine, to be used everywhere for transparent conversion of units.
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)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,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.