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

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

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