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

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

Use GaussTria

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