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

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

New results framework. Now, elements have results datasets,
which they fill up as they wish with results found in the inputs.
Then, the results dataset will be processed in OutputResults, to output
something to disk. We ended up putting the results inside the elements,
because they depend on the interpolation, this avoids partitioning of vectors,
and inputs cannot hold different time steps for the same enum!

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