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

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

Split results between element results and external results.
element results are held in the Results* dataset of each element.
external results are held in the femmodel->Results* dataset, and hold anything that cannot fit in the elements.
File size: 29.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 "../../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(int in_enum_type,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 (%s)\n",this->enum_type,EnumAsString(this->enum_type));
56 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]);
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/*FUNCTION PentaVertexInput::SpawnSingInput{{{1*/
137Input* PentaVertexInput::SpawnSingInput(int index){
138
139 /*output*/
140 SingVertexInput* outinput=NULL;
141
142 /*Create new Sing input (copy of current input)*/
143 ISSMASSERT(index<6 && index>=0);
144 outinput=new SingVertexInput(this->enum_type,this->values[index]);
145
146 /*Assign output*/
147 return outinput;
148
149}
150/*}}}*/
151/*FUNCTION PentaVertexInput::SpawnBeamInput{{{1*/
152Input* PentaVertexInput::SpawnBeamInput(int* indices){
153
154 /*output*/
155 BeamVertexInput* outinput=NULL;
156 double newvalues[2];
157
158 /*Loop over the new indices*/
159 for(int i=0;i<2;i++){
160
161 /*Check index value*/
162 ISSMASSERT(indices[i]>=0 && indices[i]<6);
163
164 /*Assign value to new input*/
165 newvalues[i]=this->values[indices[i]];
166 }
167
168 /*Create new Beam input*/
169 outinput=new BeamVertexInput(this->enum_type,&newvalues[0]);
170
171 /*Assign output*/
172 return outinput;
173
174}
175/*}}}*/
176/*FUNCTION PentaVertexInput::SpawnTriaInput{{{1*/
177Input* PentaVertexInput::SpawnTriaInput(int* indices){
178
179 /*output*/
180 TriaVertexInput* outinput=NULL;
181 double newvalues[3];
182
183 /*Loop over the new indices*/
184 for(int i=0;i<3;i++){
185
186 /*Check index value*/
187 ISSMASSERT(indices[i]>=0 && indices[i]<6);
188
189 /*Assign value to new input*/
190 newvalues[i]=this->values[indices[i]];
191 }
192
193 /*Create new Tria input*/
194 outinput=new TriaVertexInput(this->enum_type,&newvalues[0]);
195
196 /*Assign output*/
197 return outinput;
198
199}
200/*}}}*/
201/*FUNCTION PentaVertexInput::SpawnResult{{{1*/
202ElementResult* PentaVertexInput::SpawnResult(int step, double time){
203
204 return new PentaVertexElementResult(this->enum_type,this->values,step,time);
205
206}
207/*}}}*/
208
209/*Object functions*/
210/*FUNCTION PentaVertexInput::GetParameterValue(bool* pvalue) {{{1*/
211void PentaVertexInput::GetParameterValue(bool* pvalue){ISSMERROR(" not supported yet!");}
212/*}}}*/
213/*FUNCTION PentaVertexInput::GetParameterValue(int* pvalue){{{1*/
214void PentaVertexInput::GetParameterValue(int* pvalue){ISSMERROR(" not supported yet!");}
215/*}}}*/
216/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue){{{1*/
217void PentaVertexInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");}
218/*}}}*/
219/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,Node* node){{{1*/
220void PentaVertexInput::GetParameterValue(double* pvalue,Node* node){ISSMERROR(" not supported yet!");}
221/*}}}*/
222/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord){{{1*/
223void PentaVertexInput::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord){ISSMERROR(" not supported yet!");}
224/*}}}*/
225/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
226void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
227 /*P1 interpolation on Gauss point*/
228
229 /*intermediary*/
230 double l1l6[6];
231
232 /*nodal functions: */
233 GetNodalFunctionsP1(&l1l6[0],gauss);
234
235 /*Assign output pointers:*/
236 *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];
237
238}
239/*}}}*/
240/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/
241void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){ISSMERROR(" not supported yet!");}
242/*}}}*/
243/*FUNCTION PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){{{1*/
244void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){
245 /*It is assumed that output values has been correctly allocated*/
246
247 int i,j;
248 double gauss[4];
249
250 for (i=0;i<numgauss;i++){
251
252 /*Get current Gauss point coordinates*/
253 for (j=0;j<4;j++) gauss[j]=gauss_pointers[i*4+j];
254
255 /*Assign parameter value*/
256 GetParameterValue(&values[i],&gauss[0]);
257 }
258}
259/*}}}*/
260/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
261void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
262 /*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:
263 * 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;
264 * 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;
265 * 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;
266 *
267 * p is a vector of size 3x1 already allocated.
268 */
269
270 const int NDOF3=3;
271 const int numgrids=6;
272 double dh1dh6[NDOF3][numgrids];
273
274 /*Get nodal funnctions derivatives in actual coordinate system: */
275 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
276
277 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];
278 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];
279 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];
280
281}
282/*}}}*/
283/*FUNCTION PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/
284void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
285 int i,j;
286
287 const int numgrids=6;
288 const int DOFVELOCITY=3;
289 double B[8][27];
290 double B_reduced[6][DOFVELOCITY*numgrids];
291 double velocity[numgrids][DOFVELOCITY];
292
293 /*Get B matrix: */
294 GetBStokes(&B[0][0], xyz_list, gauss);
295 /*Create a reduced matrix of B to get rid of pressure */
296 for (i=0;i<6;i++){
297 for (j=0;j<3;j++){
298 B_reduced[i][j]=B[i][j];
299 }
300 for (j=4;j<7;j++){
301 B_reduced[i][j-1]=B[i][j];
302 }
303 for (j=8;j<11;j++){
304 B_reduced[i][j-2]=B[i][j];
305 }
306 for (j=12;j<15;j++){
307 B_reduced[i][j-3]=B[i][j];
308 }
309 for (j=16;j<19;j++){
310 B_reduced[i][j-4]=B[i][j];
311 }
312 for (j=20;j<23;j++){
313 B_reduced[i][j-5]=B[i][j];
314 }
315 }
316
317 /*Here, we are computing the strain rate of (vx,0,0)*/
318 for(i=0;i<numgrids;i++){
319 velocity[i][0]=this->values[i];
320 velocity[i][1]=0.0;
321 velocity[i][2]=0.0;
322 }
323 /*Multiply B by velocity, to get strain rate: */
324 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
325
326}
327/*}}}*/
328/*FUNCTION PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss) {{{1*/
329void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
330 int i,j;
331
332 const int numgrids=6;
333 const int DOFVELOCITY=3;
334 double B[8][27];
335 double B_reduced[6][DOFVELOCITY*numgrids];
336 double velocity[numgrids][DOFVELOCITY];
337
338 /*Get B matrix: */
339 GetBStokes(&B[0][0], xyz_list, gauss);
340 /*Create a reduced matrix of B to get rid of pressure */
341 for (i=0;i<6;i++){
342 for (j=0;j<3;j++){
343 B_reduced[i][j]=B[i][j];
344 }
345 for (j=4;j<7;j++){
346 B_reduced[i][j-1]=B[i][j];
347 }
348 for (j=8;j<11;j++){
349 B_reduced[i][j-2]=B[i][j];
350 }
351 for (j=12;j<15;j++){
352 B_reduced[i][j-3]=B[i][j];
353 }
354 for (j=16;j<19;j++){
355 B_reduced[i][j-4]=B[i][j];
356 }
357 for (j=20;j<23;j++){
358 B_reduced[i][j-5]=B[i][j];
359 }
360 }
361
362 /*Here, we are computing the strain rate of (0,vy,0)*/
363 for(i=0;i<numgrids;i++){
364 velocity[i][0]=0.0;
365 velocity[i][1]=this->values[i];
366 velocity[i][2]=0.0;
367 }
368 /*Multiply B by velocity, to get strain rate: */
369 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
370
371}
372/*}}}*/
373/*FUNCTION PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss) {{{1*/
374void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
375 int i,j;
376
377 const int numgrids=6;
378 const int DOFVELOCITY=3;
379 double B[8][27];
380 double B_reduced[6][DOFVELOCITY*numgrids];
381 double velocity[numgrids][DOFVELOCITY];
382
383 /*Get B matrix: */
384 GetBStokes(&B[0][0], xyz_list, gauss);
385 /*Create a reduced matrix of B to get rid of pressure */
386 for (i=0;i<6;i++){
387 for (j=0;j<3;j++){
388 B_reduced[i][j]=B[i][j];
389 }
390 for (j=4;j<7;j++){
391 B_reduced[i][j-1]=B[i][j];
392 }
393 for (j=8;j<11;j++){
394 B_reduced[i][j-2]=B[i][j];
395 }
396 for (j=12;j<15;j++){
397 B_reduced[i][j-3]=B[i][j];
398 }
399 for (j=16;j<19;j++){
400 B_reduced[i][j-4]=B[i][j];
401 }
402 for (j=20;j<23;j++){
403 B_reduced[i][j-5]=B[i][j];
404 }
405 }
406
407 /*Here, we are computing the strain rate of (0,0,vz)*/
408 for(i=0;i<numgrids;i++){
409 velocity[i][0]=0.0;
410 velocity[i][1]=0.0;
411 velocity[i][2]=this->values[i];
412 }
413
414 /*Multiply B by velocity, to get strain rate: */
415 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
416
417}
418/*}}}*/
419/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/
420void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){
421
422 int i;
423 const int numgrids=6;
424 const int NDOF2=2;
425 double B[5][NDOF2*numgrids];
426 double velocity[numgrids][NDOF2];
427
428 /*Get B matrix: */
429 GetBPattyn(&B[0][0], xyz_list, gauss);
430
431 /*Here, we are computing the strain rate of (vx,0)*/
432 for(i=0;i<numgrids;i++){
433 velocity[i][0]=this->values[i];
434 velocity[i][1]=0.0;
435 }
436
437 /*Multiply B by velocity, to get strain rate: */
438 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
439 &velocity[0][0],NDOF2*numgrids,1,0,
440 epsilonvx,0);
441
442}
443/*}}}*/
444/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss) {{{1*/
445void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){
446
447 int i;
448 const int numgrids=6;
449 const int NDOF2=2;
450 double B[5][NDOF2*numgrids];
451 double velocity[numgrids][NDOF2];
452
453 /*Get B matrix: */
454 GetBPattyn(&B[0][0], xyz_list, gauss);
455
456 /*Here, we are computing the strain rate of (0,vy)*/
457 for(i=0;i<numgrids;i++){
458 velocity[i][0]=0.0;
459 velocity[i][1]=this->values[i];
460 }
461
462 /*Multiply B by velocity, to get strain rate: */
463 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
464 &velocity[0][0],NDOF2*numgrids,1,0,
465 epsilonvy,0);
466
467}
468/*}}}*/
469/*FUNCTION PentaVertexInput::ChangeEnum(int newenumtype){{{1*/
470void PentaVertexInput::ChangeEnum(int newenumtype){
471 this->enum_type=newenumtype;
472}
473/*}}}*/
474/*FUNCTION PentaVertexInput::GetParameterAverage(double* pvalue){{{1*/
475void PentaVertexInput::GetParameterAverage(double* pvalue){
476 *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
477}
478/*}}}*/
479
480/*Intermediary*/
481/*FUNCTION PentaVertexInput::GetNodalFunctionsP1 {{{1*/
482void PentaVertexInput::GetNodalFunctionsP1(double* l1l6, double* gauss_coord){
483
484 /*This routine returns the values of the nodal functions at the gaussian point.*/
485
486 l1l6[0]=gauss_coord[0]*(1-gauss_coord[3])/2.0;
487
488 l1l6[1]=gauss_coord[1]*(1-gauss_coord[3])/2.0;
489
490 l1l6[2]=gauss_coord[2]*(1-gauss_coord[3])/2.0;
491
492 l1l6[3]=gauss_coord[0]*(1+gauss_coord[3])/2.0;
493
494 l1l6[4]=gauss_coord[1]*(1+gauss_coord[3])/2.0;
495
496 l1l6[5]=gauss_coord[2]*(1+gauss_coord[3])/2.0;
497
498}
499/*}}}*/
500/*FUNCTION PentaVertexInput::GetNodalFunctionsMINI{{{1*/
501void PentaVertexInput::GetNodalFunctionsMINI(double* l1l7, double* gauss_coord){
502
503 /*This routine returns the values of the nodal functions at the gaussian point.*/
504
505 /*First nodal function: */
506 l1l7[0]=gauss_coord[0]*(1.0-gauss_coord[3])/2.0;
507
508 /*Second nodal function: */
509 l1l7[1]=gauss_coord[1]*(1.0-gauss_coord[3])/2.0;
510
511 /*Third nodal function: */
512 l1l7[2]=gauss_coord[2]*(1.0-gauss_coord[3])/2.0;
513
514 /*Fourth nodal function: */
515 l1l7[3]=gauss_coord[0]*(1.0+gauss_coord[3])/2.0;
516
517 /*Fifth nodal function: */
518 l1l7[4]=gauss_coord[1]*(1.0+gauss_coord[3])/2.0;
519
520 /*Sixth nodal function: */
521 l1l7[5]=gauss_coord[2]*(1.0+gauss_coord[3])/2.0;
522
523 /*Seventh nodal function: */
524 l1l7[6]=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(1.0+gauss_coord[3])*(1.0-gauss_coord[3]);
525
526}
527/*}}}*/
528/*FUNCTION PentaVertexInput::GetNodalFunctionsP1Derivatives {{{1*/
529void PentaVertexInput::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss_coord){
530
531 /*This routine returns the values of the nodal functions derivatives (with respect to the actual coordinate system: */
532 int i;
533 const int NDOF3=3;
534 const int numgrids=6;
535
536 double dh1dh6_ref[NDOF3][numgrids];
537 double Jinv[NDOF3][NDOF3];
538
539 /*Get derivative values with respect to parametric coordinate system: */
540 GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss_coord);
541
542 /*Get Jacobian invert: */
543 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
544
545 /*Build dh1dh3:
546 *
547 * [dhi/dx]= Jinv*[dhi/dr]
548 * [dhi/dy] [dhi/ds]
549 * [dhi/dz] [dhi/dn]
550 */
551
552 for (i=0;i<numgrids;i++){
553 *(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];
554 *(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];
555 *(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];
556 }
557
558}
559/*}}}*/
560/*FUNCTION PentaVertexInput::GetNodalFunctionsMINIDerivatives{{{1*/
561void PentaVertexInput::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss_coord){
562
563 /*This routine returns the values of the nodal functions derivatives (with respect to the
564 * actual coordinate system: */
565
566 int i;
567
568 const int numgrids=7;
569 double dh1dh7_ref[3][numgrids];
570 double Jinv[3][3];
571
572
573 /*Get derivative values with respect to parametric coordinate system: */
574 GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss_coord);
575
576 /*Get Jacobian invert: */
577 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
578
579 /*Build dh1dh6:
580 *
581 * [dhi/dx]= Jinv'*[dhi/dr]
582 * [dhi/dy] [dhi/ds]
583 * [dhi/dz] [dhi/dzeta]
584 */
585
586 for (i=0;i<numgrids;i++){
587 *(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];
588 *(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];
589 *(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];
590 }
591
592}
593/*}}}*/
594/*FUNCTION PentaVertexInput::GetNodalFunctionsP1DerivativesReference {{{1*/
595void PentaVertexInput::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss_coord){
596
597 /*This routine returns the values of the nodal functions derivatives (with respect to the
598 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
599
600 const int numgrids=6;
601 double A1,A2,A3,z;
602
603 A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
604 A2=gauss_coord[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
605 A3=gauss_coord[2]; //third area coordinate value In term of xi and eta: A3=y/SQRT3;
606 z=gauss_coord[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
607
608
609 /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
610 *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
611 *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
612 *(dl1dl6+numgrids*2+0)=-0.5*A1;
613
614 /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
615 *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
616 *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
617 *(dl1dl6+numgrids*2+1)=-0.5*A2;
618
619 /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
620 *(dl1dl6+numgrids*0+2)=0.0;
621 *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
622 *(dl1dl6+numgrids*2+2)=-0.5*A3;
623
624 /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
625 *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
626 *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
627 *(dl1dl6+numgrids*2+3)=0.5*A1;
628
629 /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
630 *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
631 *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
632 *(dl1dl6+numgrids*2+4)=0.5*A2;
633
634 /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
635 *(dl1dl6+numgrids*0+5)=0.0;
636 *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
637 *(dl1dl6+numgrids*2+5)=0.5*A3;
638}
639/*}}}*/
640/*FUNCTION PentaVertexInput::GetNodalFunctionsMINIDerivativesReference{{{1*/
641void PentaVertexInput::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss_coord){
642
643 /*This routine returns the values of the nodal functions derivatives (with respect to the
644 * natural coordinate system) at the gaussian point. */
645
646 int numgrids=7; //six plus bubble grids
647
648 double r=gauss_coord[1]-gauss_coord[0];
649 double s=-3.0/SQRT3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
650 double zeta=gauss_coord[3];
651
652 /*First nodal function: */
653 *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
654 *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
655 *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
656
657 /*Second nodal function: */
658 *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
659 *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
660 *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
661
662 /*Third nodal function: */
663 *(dl1dl7+numgrids*0+2)=0;
664 *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
665 *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
666
667 /*Fourth nodal function: */
668 *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
669 *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
670 *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
671
672 /*Fith nodal function: */
673 *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
674 *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
675 *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
676
677 /*Sixth nodal function: */
678 *(dl1dl7+numgrids*0+5)=0;
679 *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
680 *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
681
682 /*Seventh nodal function: */
683 *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
684 *(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));
685 *(dl1dl7+numgrids*2+6)=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(-2.0*zeta);
686
687}
688/*}}}*/
689/*FUNCTION PentaVertexInput::GetJacobian {{{1*/
690void PentaVertexInput::GetJacobian(double* J, double* xyz_list,double* gauss_coord){
691
692 const int NDOF3=3;
693 int i,j;
694
695 /*The Jacobian is constant over the element, discard the gaussian points.
696 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
697
698 double A1,A2,A3; //area coordinates
699 double xi,eta,zi; //parametric coordinates
700
701 double x1,x2,x3,x4,x5,x6;
702 double y1,y2,y3,y4,y5,y6;
703 double z1,z2,z3,z4,z5,z6;
704
705 /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
706 A1=gauss_coord[0];
707 A2=gauss_coord[1];
708 A3=gauss_coord[2];
709
710 xi=A2-A1;
711 eta=SQRT3*A3;
712 zi=gauss_coord[3];
713
714 x1=*(xyz_list+3*0+0);
715 x2=*(xyz_list+3*1+0);
716 x3=*(xyz_list+3*2+0);
717 x4=*(xyz_list+3*3+0);
718 x5=*(xyz_list+3*4+0);
719 x6=*(xyz_list+3*5+0);
720
721 y1=*(xyz_list+3*0+1);
722 y2=*(xyz_list+3*1+1);
723 y3=*(xyz_list+3*2+1);
724 y4=*(xyz_list+3*3+1);
725 y5=*(xyz_list+3*4+1);
726 y6=*(xyz_list+3*5+1);
727
728 z1=*(xyz_list+3*0+2);
729 z2=*(xyz_list+3*1+2);
730 z3=*(xyz_list+3*2+2);
731 z4=*(xyz_list+3*3+2);
732 z5=*(xyz_list+3*4+2);
733 z6=*(xyz_list+3*5+2);
734
735 *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
736 *(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);
737 *(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);
738
739 *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
740 *(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);
741 *(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);
742
743 *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
744 *(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);
745 *(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);
746
747}
748/*}}}*/
749/*FUNCTION PentaVertexInput::GetJacobianInvert {{{1*/
750void PentaVertexInput::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_coord){
751
752 double Jdet;
753 const int NDOF3=3;
754
755 /*Call Jacobian routine to get the jacobian:*/
756 GetJacobian(Jinv, xyz_list, gauss_coord);
757
758 /*Invert Jacobian matrix: */
759 MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet);
760}
761/*}}}*/
762/*FUNCTION PentaVertexInput::GetBPattyn {{{1*/
763void PentaVertexInput::GetBPattyn(double* B, double* xyz_list, double* gauss_coord){
764 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
765 * For grid i, Bi can be expressed in the actual coordinate system
766 * by:
767 * Bi=[ dh/dx 0 ]
768 * [ 0 dh/dy ]
769 * [ 1/2*dh/dy 1/2*dh/dx ]
770 * [ 1/2*dh/dz 0 ]
771 * [ 0 1/2*dh/dz ]
772 * where h is the interpolation function for grid i.
773 *
774 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
775 */
776
777 int i;
778 const int numgrids=6;
779 const int NDOF3=3;
780 const int NDOF2=2;
781
782 double dh1dh6[NDOF3][numgrids];
783
784 /*Get dh1dh6 in actual coordinate system: */
785 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
786
787 /*Build B: */
788 for (i=0;i<numgrids;i++){
789 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
790 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
791
792 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
793 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
794
795 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
796 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
797
798 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];
799 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
800
801 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
802 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
803 }
804
805}
806/*}}}*/
807/*FUNCTION PentaVertexInput::GetBStokes {{{1*/
808void PentaVertexInput::GetBStokes(double* B, double* xyz_list, double* gauss_coord){
809
810 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
811 * For grid i, Bi can be expressed in the actual coordinate system
812 * by: Bi=[ dh/dx 0 0 0 ]
813 * [ 0 dh/dy 0 0 ]
814 * [ 0 0 dh/dy 0 ]
815 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ]
816 * [ 1/2*dh/dz 0 1/2*dh/dx 0 ]
817 * [ 0 1/2*dh/dz 1/2*dh/dy 0 ]
818 * [ 0 0 0 h ]
819 * [ dh/dx dh/dy dh/dz 0 ]
820 * where h is the interpolation function for grid i.
821 * Same thing for Bb except the last column that does not exist.
822 */
823
824 int i;
825 const int calculationdof=3;
826 const int numgrids=6;
827 int DOFPERGRID=4;
828
829 double dh1dh7[calculationdof][numgrids+1];
830 double l1l6[numgrids];
831
832
833 /*Get dh1dh7 in actual coordinate system: */
834 GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss_coord);
835 GetNodalFunctionsP1(l1l6, gauss_coord);
836
837 /*Build B: */
838 for (i=0;i<numgrids+1;i++){
839 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
840 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
841 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
842 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
843 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
844 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
845 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
846 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
847 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
848 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];
849 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];
850 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
851 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
852 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
853 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
854 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
855 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
856 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
857 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
858 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
859 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
860 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
861 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
862 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
863 }
864
865 for (i=0;i<numgrids;i++){ //last column not for the bubble function
866 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
867 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
868 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
869 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
870 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
871 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
872 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
873 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
874 }
875
876}
877/*}}}*/
878/*FUNCTION PentaVertexInput::SquareMin(double* psquaremin, bool process_units){{{1*/
879void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
880
881 int i;
882 const int numnodes=6;
883 double valuescopy[numnodes];
884 double squaremin;
885
886 /*First, copy values, to process units if requested: */
887 for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];
888
889 /*Process units if requested: */
890 if(process_units)NodalValuesUnitConversion(&valuescopy[0],numnodes,enum_type,parameters);
891
892 /*Now, figure out minimum of valuescopy: */
893 squaremin=pow(valuescopy[0],2);
894 for(i=1;i<numnodes;i++){
895 if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
896 }
897 /*Assign output pointers:*/
898 *psquaremin=squaremin;
899}
900/*}}}*/
901/*FUNCTION PentaVertexInput::Scale(double scale_factor){{{1*/
902void PentaVertexInput::Scale(double scale_factor){
903
904 int i;
905 const int numgrids=6;
906
907 for(i=0;i<numgrids;i++)values[i]=values[i]*scale_factor;
908}
909/*}}}*/
910/*FUNCTION PentaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
911void PentaVertexInput::AXPY(Input* xinput,double scalar){
912
913 int i;
914 const int numgrids=6;
915 PentaVertexInput* xpentavertexinput=NULL;
916
917 /*xinput is of the same type, so cast it: */
918 xpentavertexinput=(PentaVertexInput*)xinput;
919
920 /*Carry out the AXPY operation:*/
921 for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
922
923}
924/*}}}*/
925/*FUNCTION PentaVertexInput::Constrain(double cm_min, double cm_max){{{1*/
926void PentaVertexInput::Constrain(double cm_min, double cm_max){
927
928 int i;
929 const int numgrids=6;
930
931 if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
932 if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
933
934}
935/*}}}*/
936/*FUNCTION PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
937void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
938
939 const int numvertices=6;
940 VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);
941
942
943}
944/*}}}*/
Note: See TracBrowser for help on using the repository browser.