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

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

Serious cleanup of the control solutoin
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*/
202Result* PentaVertexInput::SpawnResult(int step, double time){
203
204 return new PentaVertexResult(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.