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

Last change on this file since 4502 was 4502, checked in by seroussi, 15 years ago

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