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

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

New results API
File size: 29.0 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 "./InputLocal.h"
14#include "../objects.h"
15#include "../../EnumDefinitions/EnumDefinitions.h"
16#include "../../shared/shared.h"
17#include "../../DataSet/DataSet.h"
18#include "../../include/include.h"
19
20/*Object constructors and destructor*/
21/*FUNCTION PentaVertexInput::PentaVertexInput(){{{1*/
22PentaVertexInput::PentaVertexInput(){
23 return;
24}
25/*}}}*/
26/*FUNCTION PentaVertexInput::PentaVertexInput(int in_enum_type,double* values){{{1*/
27PentaVertexInput::PentaVertexInput(int in_enum_type,double* in_values){
28
29 enum_type=in_enum_type;
30 values[0]=in_values[0];
31 values[1]=in_values[1];
32 values[2]=in_values[2];
33 values[3]=in_values[3];
34 values[4]=in_values[4];
35 values[5]=in_values[5];
36}
37/*}}}*/
38/*FUNCTION PentaVertexInput::~PentaVertexInput(){{{1*/
39PentaVertexInput::~PentaVertexInput(){
40 return;
41}
42/*}}}*/
43
44/*Object management*/
45/*FUNCTION PentaVertexInput::copy{{{1*/
46Object* PentaVertexInput::copy() {
47
48 return new PentaVertexInput(this->enum_type,this->values);
49
50}
51/*}}}*/
52/*FUNCTION PentaVertexInput::DeepEcho{{{1*/
53void PentaVertexInput::DeepEcho(void){
54
55 printf("PentaVertexInput:\n");
56 printf(" enum: %i (%s)\n",this->enum_type,EnumAsString(this->enum_type));
57 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]);
58}
59/*}}}*/
60/*FUNCTION PentaVertexInput::Demarshall{{{1*/
61void PentaVertexInput::Demarshall(char** pmarshalled_dataset){
62
63 char* marshalled_dataset=NULL;
64 int i;
65
66 /*recover marshalled_dataset: */
67 marshalled_dataset=*pmarshalled_dataset;
68
69 /*this time, no need to get enum type, the pointer directly points to the beginning of the
70 *object data (thanks to DataSet::Demarshall):*/
71 memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
72 memcpy(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values);
73
74 /*return: */
75 *pmarshalled_dataset=marshalled_dataset;
76 return;
77}
78/*}}}*/
79/*FUNCTION PentaVertexInput::Echo {{{1*/
80void PentaVertexInput::Echo(void){
81 this->DeepEcho();
82}
83/*}}}*/
84/*FUNCTION PentaVertexInput::Enum{{{1*/
85int PentaVertexInput::Enum(void){
86
87 return PentaVertexInputEnum;
88
89}
90/*}}}*/
91/*FUNCTION PentaVertexInput::EnumType{{{1*/
92int PentaVertexInput::EnumType(void){
93
94 return this->enum_type;
95
96}
97/*}}}*/
98/*FUNCTION PentaVertexInput::Id{{{1*/
99int PentaVertexInput::Id(void){ return -1; }
100/*}}}*/
101/*FUNCTION PentaVertexInput::Marshall{{{1*/
102void PentaVertexInput::Marshall(char** pmarshalled_dataset){
103
104 char* marshalled_dataset=NULL;
105 int enum_value=0;
106
107 /*recover marshalled_dataset: */
108 marshalled_dataset=*pmarshalled_dataset;
109
110 /*get enum value of PentaVertexInput: */
111 enum_value=PentaVertexInputEnum;
112
113 /*marshall enum: */
114 memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
115
116 /*marshall PentaVertexInput data: */
117 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
118 memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values);
119
120 *pmarshalled_dataset=marshalled_dataset;
121}
122/*}}}*/
123/*FUNCTION PentaVertexInput::MarshallSize{{{1*/
124int PentaVertexInput::MarshallSize(){
125
126 return sizeof(values)+
127 +sizeof(enum_type)+
128 +sizeof(int); //sizeof(int) for enum value
129}
130/*}}}*/
131/*FUNCTION PentaVertexInput::MyRank{{{1*/
132int PentaVertexInput::MyRank(void){
133 extern int my_rank;
134 return my_rank;
135}
136/*}}}*/
137/*FUNCTION PentaVertexInput::SpawnSingInput{{{1*/
138Input* PentaVertexInput::SpawnSingInput(int index){
139
140 /*output*/
141 SingVertexInput* outinput=NULL;
142
143 /*Create new Sing input (copy of current input)*/
144 ISSMASSERT(index<6 && index>=0);
145 outinput=new SingVertexInput(this->enum_type,this->values[index]);
146
147 /*Assign output*/
148 return outinput;
149
150}
151/*}}}*/
152/*FUNCTION PentaVertexInput::SpawnBeamInput{{{1*/
153Input* PentaVertexInput::SpawnBeamInput(int* indices){
154
155 /*output*/
156 BeamVertexInput* outinput=NULL;
157 double newvalues[2];
158
159 /*Loop over the new indices*/
160 for(int i=0;i<2;i++){
161
162 /*Check index value*/
163 ISSMASSERT(indices[i]>=0 && indices[i]<6);
164
165 /*Assign value to new input*/
166 newvalues[i]=this->values[indices[i]];
167 }
168
169 /*Create new Beam input*/
170 outinput=new BeamVertexInput(this->enum_type,&newvalues[0]);
171
172 /*Assign output*/
173 return outinput;
174
175}
176/*}}}*/
177/*FUNCTION PentaVertexInput::SpawnTriaInput{{{1*/
178Input* PentaVertexInput::SpawnTriaInput(int* indices){
179
180 /*output*/
181 TriaVertexInput* outinput=NULL;
182 double newvalues[3];
183
184 /*Loop over the new indices*/
185 for(int i=0;i<3;i++){
186
187 /*Check index value*/
188 ISSMASSERT(indices[i]>=0 && indices[i]<6);
189
190 /*Assign value to new input*/
191 newvalues[i]=this->values[indices[i]];
192 }
193
194 /*Create new Tria input*/
195 outinput=new TriaVertexInput(this->enum_type,&newvalues[0]);
196
197 /*Assign output*/
198 return outinput;
199
200}
201/*}}}*/
202
203/*Object functions*/
204/*FUNCTION PentaVertexInput::GetParameterValue(bool* pvalue) {{{1*/
205void PentaVertexInput::GetParameterValue(bool* pvalue){ISSMERROR(" not supported yet!");}
206/*}}}*/
207/*FUNCTION PentaVertexInput::GetParameterValue(int* pvalue){{{1*/
208void PentaVertexInput::GetParameterValue(int* pvalue){ISSMERROR(" not supported yet!");}
209/*}}}*/
210/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue){{{1*/
211void PentaVertexInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");}
212/*}}}*/
213/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,Node* node){{{1*/
214void PentaVertexInput::GetParameterValue(double* pvalue,Node* node){ISSMERROR(" not supported yet!");}
215/*}}}*/
216/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord){{{1*/
217void PentaVertexInput::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord){ISSMERROR(" not supported yet!");}
218/*}}}*/
219/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
220void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
221 /*P1 interpolation on Gauss point*/
222
223 /*intermediary*/
224 double l1l6[6];
225
226 /*nodal functions: */
227 GetNodalFunctionsP1(&l1l6[0],gauss);
228
229 /*Assign output pointers:*/
230 *pvalue=l1l6[0]*values[0]+l1l6[1]*values[1]+l1l6[2]*values[2]+l1l6[3]*values[3]+l1l6[4]*values[4]+l1l6[5]*values[5];
231
232}
233/*}}}*/
234/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/
235void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){ISSMERROR(" not supported yet!");}
236/*}}}*/
237/*FUNCTION PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){{{1*/
238void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){
239 /*It is assumed that output values has been correctly allocated*/
240
241 int i,j;
242 double gauss[4];
243
244 for (i=0;i<numgauss;i++){
245
246 /*Get current Gauss point coordinates*/
247 for (j=0;j<4;j++) gauss[j]=gauss_pointers[i*4+j];
248
249 /*Assign parameter value*/
250 GetParameterValue(&values[i],&gauss[0]);
251 }
252}
253/*}}}*/
254/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
255void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
256 /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
257 * dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
258 * dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
259 * dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
260 *
261 * p is a vector of size 3x1 already allocated.
262 */
263
264 const int NDOF3=3;
265 const int numgrids=6;
266 double dh1dh6[NDOF3][numgrids];
267
268 /*Get nodal funnctions derivatives in actual coordinate system: */
269 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
270
271 p[0]=this->values[0]*dh1dh6[0][0]+this->values[1]*dh1dh6[0][1]+this->values[2]*dh1dh6[0][2]+this->values[3]*dh1dh6[0][3]+this->values[4]*dh1dh6[0][4]+this->values[5]*dh1dh6[0][5];
272 p[1]=this->values[0]*dh1dh6[1][0]+this->values[1]*dh1dh6[1][1]+this->values[2]*dh1dh6[1][2]+this->values[3]*dh1dh6[1][3]+this->values[4]*dh1dh6[1][4]+this->values[5]*dh1dh6[1][5];
273 p[2]=this->values[0]*dh1dh6[2][0]+this->values[1]*dh1dh6[2][1]+this->values[2]*dh1dh6[2][2]+this->values[3]*dh1dh6[2][3]+this->values[4]*dh1dh6[2][4]+this->values[5]*dh1dh6[2][5];
274
275}
276/*}}}*/
277/*FUNCTION PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/
278void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
279 int i,j;
280
281 const int numgrids=6;
282 const int DOFVELOCITY=3;
283 double B[8][27];
284 double B_reduced[6][DOFVELOCITY*numgrids];
285 double velocity[numgrids][DOFVELOCITY];
286
287 /*Get B matrix: */
288 GetBStokes(&B[0][0], xyz_list, gauss);
289 /*Create a reduced matrix of B to get rid of pressure */
290 for (i=0;i<6;i++){
291 for (j=0;j<3;j++){
292 B_reduced[i][j]=B[i][j];
293 }
294 for (j=4;j<7;j++){
295 B_reduced[i][j-1]=B[i][j];
296 }
297 for (j=8;j<11;j++){
298 B_reduced[i][j-2]=B[i][j];
299 }
300 for (j=12;j<15;j++){
301 B_reduced[i][j-3]=B[i][j];
302 }
303 for (j=16;j<19;j++){
304 B_reduced[i][j-4]=B[i][j];
305 }
306 for (j=20;j<23;j++){
307 B_reduced[i][j-5]=B[i][j];
308 }
309 }
310
311 /*Here, we are computing the strain rate of (vx,0,0)*/
312 for(i=0;i<numgrids;i++){
313 velocity[i][0]=this->values[i];
314 velocity[i][1]=0.0;
315 velocity[i][2]=0.0;
316 }
317 /*Multiply B by velocity, to get strain rate: */
318 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
319
320}
321/*}}}*/
322/*FUNCTION PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss) {{{1*/
323void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
324 int i,j;
325
326 const int numgrids=6;
327 const int DOFVELOCITY=3;
328 double B[8][27];
329 double B_reduced[6][DOFVELOCITY*numgrids];
330 double velocity[numgrids][DOFVELOCITY];
331
332 /*Get B matrix: */
333 GetBStokes(&B[0][0], xyz_list, gauss);
334 /*Create a reduced matrix of B to get rid of pressure */
335 for (i=0;i<6;i++){
336 for (j=0;j<3;j++){
337 B_reduced[i][j]=B[i][j];
338 }
339 for (j=4;j<7;j++){
340 B_reduced[i][j-1]=B[i][j];
341 }
342 for (j=8;j<11;j++){
343 B_reduced[i][j-2]=B[i][j];
344 }
345 for (j=12;j<15;j++){
346 B_reduced[i][j-3]=B[i][j];
347 }
348 for (j=16;j<19;j++){
349 B_reduced[i][j-4]=B[i][j];
350 }
351 for (j=20;j<23;j++){
352 B_reduced[i][j-5]=B[i][j];
353 }
354 }
355
356 /*Here, we are computing the strain rate of (0,vy,0)*/
357 for(i=0;i<numgrids;i++){
358 velocity[i][0]=0.0;
359 velocity[i][1]=this->values[i];
360 velocity[i][2]=0.0;
361 }
362 /*Multiply B by velocity, to get strain rate: */
363 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
364
365}
366/*}}}*/
367/*FUNCTION PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss) {{{1*/
368void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
369 int i,j;
370
371 const int numgrids=6;
372 const int DOFVELOCITY=3;
373 double B[8][27];
374 double B_reduced[6][DOFVELOCITY*numgrids];
375 double velocity[numgrids][DOFVELOCITY];
376
377 /*Get B matrix: */
378 GetBStokes(&B[0][0], xyz_list, gauss);
379 /*Create a reduced matrix of B to get rid of pressure */
380 for (i=0;i<6;i++){
381 for (j=0;j<3;j++){
382 B_reduced[i][j]=B[i][j];
383 }
384 for (j=4;j<7;j++){
385 B_reduced[i][j-1]=B[i][j];
386 }
387 for (j=8;j<11;j++){
388 B_reduced[i][j-2]=B[i][j];
389 }
390 for (j=12;j<15;j++){
391 B_reduced[i][j-3]=B[i][j];
392 }
393 for (j=16;j<19;j++){
394 B_reduced[i][j-4]=B[i][j];
395 }
396 for (j=20;j<23;j++){
397 B_reduced[i][j-5]=B[i][j];
398 }
399 }
400
401 /*Here, we are computing the strain rate of (0,0,vz)*/
402 for(i=0;i<numgrids;i++){
403 velocity[i][0]=0.0;
404 velocity[i][1]=0.0;
405 velocity[i][2]=this->values[i];
406 }
407
408 /*Multiply B by velocity, to get strain rate: */
409 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
410
411}
412/*}}}*/
413/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/
414void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){
415
416 int i;
417 const int numgrids=6;
418 const int NDOF2=2;
419 double B[5][NDOF2*numgrids];
420 double velocity[numgrids][NDOF2];
421
422 /*Get B matrix: */
423 GetBPattyn(&B[0][0], xyz_list, gauss);
424
425 /*Here, we are computing the strain rate of (vx,0)*/
426 for(i=0;i<numgrids;i++){
427 velocity[i][0]=this->values[i];
428 velocity[i][1]=0.0;
429 }
430
431 /*Multiply B by velocity, to get strain rate: */
432 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
433 &velocity[0][0],NDOF2*numgrids,1,0,
434 epsilonvx,0);
435
436}
437/*}}}*/
438/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss) {{{1*/
439void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){
440
441 int i;
442 const int numgrids=6;
443 const int NDOF2=2;
444 double B[5][NDOF2*numgrids];
445 double velocity[numgrids][NDOF2];
446
447 /*Get B matrix: */
448 GetBPattyn(&B[0][0], xyz_list, gauss);
449
450 /*Here, we are computing the strain rate of (0,vy)*/
451 for(i=0;i<numgrids;i++){
452 velocity[i][0]=0.0;
453 velocity[i][1]=this->values[i];
454 }
455
456 /*Multiply B by velocity, to get strain rate: */
457 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
458 &velocity[0][0],NDOF2*numgrids,1,0,
459 epsilonvy,0);
460
461}
462/*}}}*/
463/*FUNCTION PentaVertexInput::ChangeEnum(int newenumtype){{{1*/
464void PentaVertexInput::ChangeEnum(int newenumtype){
465 this->enum_type=newenumtype;
466}
467/*}}}*/
468/*FUNCTION PentaVertexInput::GetParameterAverage(double* pvalue){{{1*/
469void PentaVertexInput::GetParameterAverage(double* pvalue){
470 *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
471}
472/*}}}*/
473
474/*Intermediary*/
475/*FUNCTION PentaVertexInput::GetNodalFunctionsP1 {{{1*/
476void PentaVertexInput::GetNodalFunctionsP1(double* l1l6, double* gauss_coord){
477
478 /*This routine returns the values of the nodal functions at the gaussian point.*/
479
480 l1l6[0]=gauss_coord[0]*(1-gauss_coord[3])/2.0;
481
482 l1l6[1]=gauss_coord[1]*(1-gauss_coord[3])/2.0;
483
484 l1l6[2]=gauss_coord[2]*(1-gauss_coord[3])/2.0;
485
486 l1l6[3]=gauss_coord[0]*(1+gauss_coord[3])/2.0;
487
488 l1l6[4]=gauss_coord[1]*(1+gauss_coord[3])/2.0;
489
490 l1l6[5]=gauss_coord[2]*(1+gauss_coord[3])/2.0;
491
492}
493/*}}}*/
494/*FUNCTION PentaVertexInput::GetNodalFunctionsMINI{{{1*/
495void PentaVertexInput::GetNodalFunctionsMINI(double* l1l7, double* gauss_coord){
496
497 /*This routine returns the values of the nodal functions at the gaussian point.*/
498
499 /*First nodal function: */
500 l1l7[0]=gauss_coord[0]*(1.0-gauss_coord[3])/2.0;
501
502 /*Second nodal function: */
503 l1l7[1]=gauss_coord[1]*(1.0-gauss_coord[3])/2.0;
504
505 /*Third nodal function: */
506 l1l7[2]=gauss_coord[2]*(1.0-gauss_coord[3])/2.0;
507
508 /*Fourth nodal function: */
509 l1l7[3]=gauss_coord[0]*(1.0+gauss_coord[3])/2.0;
510
511 /*Fifth nodal function: */
512 l1l7[4]=gauss_coord[1]*(1.0+gauss_coord[3])/2.0;
513
514 /*Sixth nodal function: */
515 l1l7[5]=gauss_coord[2]*(1.0+gauss_coord[3])/2.0;
516
517 /*Seventh nodal function: */
518 l1l7[6]=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(1.0+gauss_coord[3])*(1.0-gauss_coord[3]);
519
520}
521/*}}}*/
522/*FUNCTION PentaVertexInput::GetNodalFunctionsP1Derivatives {{{1*/
523void PentaVertexInput::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss_coord){
524
525 /*This routine returns the values of the nodal functions derivatives (with respect to the actual coordinate system: */
526 int i;
527 const int NDOF3=3;
528 const int numgrids=6;
529
530 double dh1dh6_ref[NDOF3][numgrids];
531 double Jinv[NDOF3][NDOF3];
532
533 /*Get derivative values with respect to parametric coordinate system: */
534 GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss_coord);
535
536 /*Get Jacobian invert: */
537 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
538
539 /*Build dh1dh3:
540 *
541 * [dhi/dx]= Jinv*[dhi/dr]
542 * [dhi/dy] [dhi/ds]
543 * [dhi/dz] [dhi/dn]
544 */
545
546 for (i=0;i<numgrids;i++){
547 *(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
548 *(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
549 *(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
550 }
551
552}
553/*}}}*/
554/*FUNCTION PentaVertexInput::GetNodalFunctionsMINIDerivatives{{{1*/
555void PentaVertexInput::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss_coord){
556
557 /*This routine returns the values of the nodal functions derivatives (with respect to the
558 * actual coordinate system: */
559
560 int i;
561
562 const int numgrids=7;
563 double dh1dh7_ref[3][numgrids];
564 double Jinv[3][3];
565
566
567 /*Get derivative values with respect to parametric coordinate system: */
568 GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss_coord);
569
570 /*Get Jacobian invert: */
571 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
572
573 /*Build dh1dh6:
574 *
575 * [dhi/dx]= Jinv'*[dhi/dr]
576 * [dhi/dy] [dhi/ds]
577 * [dhi/dz] [dhi/dzeta]
578 */
579
580 for (i=0;i<numgrids;i++){
581 *(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
582 *(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
583 *(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
584 }
585
586}
587/*}}}*/
588/*FUNCTION PentaVertexInput::GetNodalFunctionsP1DerivativesReference {{{1*/
589void PentaVertexInput::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss_coord){
590
591 /*This routine returns the values of the nodal functions derivatives (with respect to the
592 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
593
594 const int numgrids=6;
595 double A1,A2,A3,z;
596
597 A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
598 A2=gauss_coord[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
599 A3=gauss_coord[2]; //third area coordinate value In term of xi and eta: A3=y/SQRT3;
600 z=gauss_coord[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
601
602
603 /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
604 *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
605 *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
606 *(dl1dl6+numgrids*2+0)=-0.5*A1;
607
608 /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
609 *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
610 *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
611 *(dl1dl6+numgrids*2+1)=-0.5*A2;
612
613 /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
614 *(dl1dl6+numgrids*0+2)=0.0;
615 *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
616 *(dl1dl6+numgrids*2+2)=-0.5*A3;
617
618 /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
619 *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
620 *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
621 *(dl1dl6+numgrids*2+3)=0.5*A1;
622
623 /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
624 *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
625 *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
626 *(dl1dl6+numgrids*2+4)=0.5*A2;
627
628 /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
629 *(dl1dl6+numgrids*0+5)=0.0;
630 *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
631 *(dl1dl6+numgrids*2+5)=0.5*A3;
632}
633/*}}}*/
634/*FUNCTION PentaVertexInput::GetNodalFunctionsMINIDerivativesReference{{{1*/
635void PentaVertexInput::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss_coord){
636
637 /*This routine returns the values of the nodal functions derivatives (with respect to the
638 * natural coordinate system) at the gaussian point. */
639
640 int numgrids=7; //six plus bubble grids
641
642 double r=gauss_coord[1]-gauss_coord[0];
643 double s=-3.0/SQRT3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
644 double zeta=gauss_coord[3];
645
646 /*First nodal function: */
647 *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
648 *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
649 *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
650
651 /*Second nodal function: */
652 *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
653 *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
654 *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
655
656 /*Third nodal function: */
657 *(dl1dl7+numgrids*0+2)=0;
658 *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
659 *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
660
661 /*Fourth nodal function: */
662 *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
663 *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
664 *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
665
666 /*Fith nodal function: */
667 *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
668 *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
669 *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
670
671 /*Sixth nodal function: */
672 *(dl1dl7+numgrids*0+5)=0;
673 *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
674 *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
675
676 /*Seventh nodal function: */
677 *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
678 *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
679 *(dl1dl7+numgrids*2+6)=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(-2.0*zeta);
680
681}
682/*}}}*/
683/*FUNCTION PentaVertexInput::GetJacobian {{{1*/
684void PentaVertexInput::GetJacobian(double* J, double* xyz_list,double* gauss_coord){
685
686 const int NDOF3=3;
687 int i,j;
688
689 /*The Jacobian is constant over the element, discard the gaussian points.
690 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
691
692 double A1,A2,A3; //area coordinates
693 double xi,eta,zi; //parametric coordinates
694
695 double x1,x2,x3,x4,x5,x6;
696 double y1,y2,y3,y4,y5,y6;
697 double z1,z2,z3,z4,z5,z6;
698
699 /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
700 A1=gauss_coord[0];
701 A2=gauss_coord[1];
702 A3=gauss_coord[2];
703
704 xi=A2-A1;
705 eta=SQRT3*A3;
706 zi=gauss_coord[3];
707
708 x1=*(xyz_list+3*0+0);
709 x2=*(xyz_list+3*1+0);
710 x3=*(xyz_list+3*2+0);
711 x4=*(xyz_list+3*3+0);
712 x5=*(xyz_list+3*4+0);
713 x6=*(xyz_list+3*5+0);
714
715 y1=*(xyz_list+3*0+1);
716 y2=*(xyz_list+3*1+1);
717 y3=*(xyz_list+3*2+1);
718 y4=*(xyz_list+3*3+1);
719 y5=*(xyz_list+3*4+1);
720 y6=*(xyz_list+3*5+1);
721
722 z1=*(xyz_list+3*0+2);
723 z2=*(xyz_list+3*1+2);
724 z3=*(xyz_list+3*2+2);
725 z4=*(xyz_list+3*3+2);
726 z5=*(xyz_list+3*4+2);
727 z6=*(xyz_list+3*5+2);
728
729 *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
730 *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
731 *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
732
733 *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
734 *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
735 *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
736
737 *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
738 *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
739 *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
740
741}
742/*}}}*/
743/*FUNCTION PentaVertexInput::GetJacobianInvert {{{1*/
744void PentaVertexInput::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_coord){
745
746 double Jdet;
747 const int NDOF3=3;
748
749 /*Call Jacobian routine to get the jacobian:*/
750 GetJacobian(Jinv, xyz_list, gauss_coord);
751
752 /*Invert Jacobian matrix: */
753 MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet);
754}
755/*}}}*/
756/*FUNCTION PentaVertexInput::GetBPattyn {{{1*/
757void PentaVertexInput::GetBPattyn(double* B, double* xyz_list, double* gauss_coord){
758 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
759 * For grid i, Bi can be expressed in the actual coordinate system
760 * by:
761 * Bi=[ dh/dx 0 ]
762 * [ 0 dh/dy ]
763 * [ 1/2*dh/dy 1/2*dh/dx ]
764 * [ 1/2*dh/dz 0 ]
765 * [ 0 1/2*dh/dz ]
766 * where h is the interpolation function for grid i.
767 *
768 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
769 */
770
771 int i;
772 const int numgrids=6;
773 const int NDOF3=3;
774 const int NDOF2=2;
775
776 double dh1dh6[NDOF3][numgrids];
777
778 /*Get dh1dh6 in actual coordinate system: */
779 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
780
781 /*Build B: */
782 for (i=0;i<numgrids;i++){
783 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
784 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
785
786 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
787 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
788
789 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
790 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
791
792 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];
793 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
794
795 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
796 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
797 }
798
799}
800/*}}}*/
801/*FUNCTION PentaVertexInput::GetBStokes {{{1*/
802void PentaVertexInput::GetBStokes(double* B, double* xyz_list, double* gauss_coord){
803
804 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
805 * For grid i, Bi can be expressed in the actual coordinate system
806 * by: Bi=[ dh/dx 0 0 0 ]
807 * [ 0 dh/dy 0 0 ]
808 * [ 0 0 dh/dy 0 ]
809 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ]
810 * [ 1/2*dh/dz 0 1/2*dh/dx 0 ]
811 * [ 0 1/2*dh/dz 1/2*dh/dy 0 ]
812 * [ 0 0 0 h ]
813 * [ dh/dx dh/dy dh/dz 0 ]
814 * where h is the interpolation function for grid i.
815 * Same thing for Bb except the last column that does not exist.
816 */
817
818 int i;
819 const int calculationdof=3;
820 const int numgrids=6;
821 int DOFPERGRID=4;
822
823 double dh1dh7[calculationdof][numgrids+1];
824 double l1l6[numgrids];
825
826
827 /*Get dh1dh7 in actual coordinate system: */
828 GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss_coord);
829 GetNodalFunctionsP1(l1l6, gauss_coord);
830
831 /*Build B: */
832 for (i=0;i<numgrids+1;i++){
833 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
834 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
835 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
836 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
837 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
838 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
839 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
840 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
841 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
842 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];
843 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];
844 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
845 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
846 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
847 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
848 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
849 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
850 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
851 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
852 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
853 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
854 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
855 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
856 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
857 }
858
859 for (i=0;i<numgrids;i++){ //last column not for the bubble function
860 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
861 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
862 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
863 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
864 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
865 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
866 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
867 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
868 }
869
870}
871/*}}}*/
872/*FUNCTION PentaVertexInput::PatchSize(void);{{{1*/
873int PentaVertexInput::PatchSize(void){
874
875 /*Return the number of nodal values this input holds, so that
876 * results can be correctl dimensionned. See InputToResultsx
877 * module for more explanations: */
878 return 6;
879}
880/*}}}*/
881/*FUNCTION PentaVertexInput::PatchFill(double* patches, int max_vertices,Parameters* parameters);{{{1*/
882void PentaVertexInput::PatchFill(double* patches, int max_vertices,Parameters* parameters){
883
884 /*A patch is made of the following information:
885 * element_id interpolation_type vertex_ids values.
886 * For example:
887
888 1 P0 1 2 4 11 12 14 4.5 NaN NaN NaN NaN NaN
889 2 P1 2 4 5 12 14 15 4.5 23.3 23.3 4.2 4.2 3.2
890 3 P0 5 2 1 15 12 11 5.5 NaN NaN NaN NaN NaN
891 4 P1 2 3 5 12 13 15 4.5 30.2 322.2 4.2 3.2 8.3
892 ...
893
894 Here, we fill the info relevant to the input, ie interpolation_type and nodal values: */
895
896 int i;
897
898
899 patches[1]=P1Enum;
900 for(i=0;i<6;i++)patches[2+max_vertices+i]=values[i]; //start of nodal values is at position 2+max_vertices (2 for id and interpolation_type) and max_vertices for vertices ids.
901
902 /*Now, post-processing (essentially, unit conversion): */
903 ProcessResults(patches+2+max_vertices,6,this->enum_type,parameters);
904
905}
906/*}}}*/
907
908
Note: See TracBrowser for help on using the repository browser.