Ice Sheet System Model  4.18
Code documentation
PentaRef.cpp
Go to the documentation of this file.
1 
5 /*Headers:*/
6 /*{{{*/
7 #ifdef HAVE_CONFIG_H
8  #include <config.h>
9 #else
10 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11 #endif
12 
13 #include "../classes.h"
14 #include "../../shared/shared.h"
15 /*}}}*/
16 
17 /*Element macros*/
18 #define NUMNODESP0 1
19 #define NUMNODESP1 6
20 #define NUMNODESP1_2d 3
21 #define NUMNODESP1b 7
22 #define NUMNODESP1xP2 9
23 #define NUMNODESP1xP3 12
24 #define NUMNODESP1xP4 15
25 #define NUMNODESP2xP1 12
26 #define NUMNODESP2 18
27 #define NUMNODESP2b 19
28 #define NUMNODESP2xP4 30
29 #define NUMNODESMAX 30
30 
31 /*Object constructors and destructor*/
32 PentaRef::PentaRef(){/*{{{*/
33 }
34 /*}}}*/
35 PentaRef::~PentaRef(){/*{{{*/
36 }
37 /*}}}*/
38 
39 /*Reference Element numerics*/
40 void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
41 
42  /*Output*/
43  int numindices;
44  int* indices = NULL;
45 
46  switch(finiteelement){
47  case P1Enum: case P1DGEnum:
48  numindices = 3;
49  indices = xNew<int>(numindices);
50  indices[0] = 0;
51  indices[1] = 1;
52  indices[2] = 2;
53  break;
54  case P1bubbleEnum: case P1bubblecondensedEnum:
55  numindices = 3;
56  indices = xNew<int>(numindices);
57  indices[0] = 0;
58  indices[1] = 1;
59  indices[2] = 2;
60  break;
61  case P2xP1Enum:
62  numindices = 6;
63  indices = xNew<int>(numindices);
64  indices[0] = 0;
65  indices[1] = 1;
66  indices[2] = 2;
67  indices[3] = 6;
68  indices[4] = 7;
69  indices[5] = 8;
70  break;
71  case P1xP2Enum:
72  numindices = 3;
73  indices = xNew<int>(numindices);
74  indices[0] = 0;
75  indices[1] = 1;
76  indices[2] = 2;
77  break;
78  case P1xP3Enum:
79  numindices = 3;
80  indices = xNew<int>(numindices);
81  indices[0] = 0;
82  indices[1] = 1;
83  indices[2] = 2;
84  break;
85  case P1xP4Enum:
86  numindices = 3;
87  indices = xNew<int>(numindices);
88  indices[0] = 0;
89  indices[1] = 1;
90  indices[2] = 2;
91  break;
92  case P2Enum:
93  numindices = 6;
94  indices = xNew<int>(numindices);
95  indices[0] = 0;
96  indices[1] = 1;
97  indices[2] = 2;
98  indices[3] = 9;
99  indices[4] = 10;
100  indices[5] = 11;
101  break;
102  case P2bubbleEnum:
103  numindices = 6;
104  indices = xNew<int>(numindices);
105  indices[0] = 0;
106  indices[1] = 1;
107  indices[2] = 2;
108  indices[3] = 9;
109  indices[4] = 10;
110  indices[5] = 11;
111  break;
112  case P2xP4Enum:
113  numindices = 6;
114  indices = xNew<int>(numindices);
115  indices[0] = 0;
116  indices[1] = 1;
117  indices[2] = 2;
118  indices[3] = 9;
119  indices[4] = 10;
120  indices[5] = 11;
121  break;
122  default:
123  _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
124  }
125 
126  /*Assign output pointer*/
127  *pnumindices = numindices;
128  *pindices = indices;
129 }
130 /*}}}*/
131 void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
132  /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
133  * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
134  * gaussian point specified by gauss_coord:
135  * 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;
136  * 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;
137  * 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;
138  *
139  * p is a vector of size 3x1 already allocated.
140  *
141  * WARNING: For a significant gain in performance, it is better to use
142  * static memory allocation instead of dynamic.
143  */
144 
145  /*Allocate derivatives of basis functions*/
146  IssmDouble dbasis[3*NUMNODESMAX];
147 
148  /*Fetch number of nodes for this finite element*/
149  int numnodes = this->NumberofNodes(finiteelement);
150  _assert_(numnodes<=NUMNODESMAX);
151 
152  /*Get basis functions derivatives at this point*/
153  GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
154 
155  /*Calculate parameter for this Gauss point*/
156  IssmDouble dpx=0.;
157  IssmDouble dpy=0.;
158  IssmDouble dpz=0.;
159  for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
160  for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
161  for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
162 
163  /*Assign values*/
164  p[0]=dpx;
165  p[1]=dpy;
166  p[2]=dpz;
167 }
168 /*}}}*/
169 void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){/*{{{*/
170  /* WARNING: For a significant gain in performance, it is better to use
171  * static memory allocation instead of dynamic.*/
172 
173  /*Allocate basis functions*/
174  IssmDouble basis[NUMNODESMAX];
175 
176  /*Fetch number of nodes for this finite element*/
177  int numnodes = this->NumberofNodes(finiteelement);
178  _assert_(numnodes<=NUMNODESMAX);
179 
180  /*Get basis functions at this point*/
181  GetNodalFunctions(&basis[0],gauss,finiteelement);
182 
183  /*Calculate parameter for this Gauss point*/
184  IssmDouble value =0.;
185  for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
186 
187  /*Assign output pointer*/
188  *pvalue = value;
189 }
190 /*}}}*/
191 void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
192  /*The Jacobian is constant over the element, discard the gaussian points.
193  * J is assumed to have been allocated of size 2x2.*/
194 
195  IssmDouble A1,A2,A3; // area coordinates
196  IssmDouble xi,eta,zi; // parametric coordinates
197  IssmDouble x1,x2,x3,x4,x5,x6;
198  IssmDouble y1,y2,y3,y4,y5,y6;
199  IssmDouble z1,z2,z3,z4,z5,z6;
200  IssmDouble j_const_reciprocal; // SQRT3/12.0
201 
202  /*Cast gauss to GaussPenta*/
203  _assert_(gauss_in->Enum()==GaussPentaEnum);
204  GaussPenta* gauss = xDynamicCast<GaussPenta*>(gauss_in);
205 
206  /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
207  A1 = gauss->coord1;
208  A2 = gauss->coord2;
209  A3 = gauss->coord3;
210  xi = A2-A1;
211  eta = SQRT3*A3;
212  zi = gauss->coord4;
213 
214  x1=xyz_list[3*0+0];
215  x2=xyz_list[3*1+0];
216  x3=xyz_list[3*2+0];
217  x4=xyz_list[3*3+0];
218  x5=xyz_list[3*4+0];
219  x6=xyz_list[3*5+0];
220 
221  y1=xyz_list[3*0+1];
222  y2=xyz_list[3*1+1];
223  y3=xyz_list[3*2+1];
224  y4=xyz_list[3*3+1];
225  y5=xyz_list[3*4+1];
226  y6=xyz_list[3*5+1];
227 
228  z1=xyz_list[3*0+2];
229  z2=xyz_list[3*1+2];
230  z3=xyz_list[3*2+2];
231  z4=xyz_list[3*3+2];
232  z5=xyz_list[3*4+2];
233  z6=xyz_list[3*5+2];
234 
235  j_const_reciprocal=SQRT3/12;
236 
237  J[3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
238  J[3*1+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*zi+j_const_reciprocal*(-x1-x2+2*x3-x4-x5+2*x6);
239  J[3*2+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*eta+0.25*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
240 
241  J[3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
242  J[3*1+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*zi+j_const_reciprocal*(-y1-y2+2*y3-y4-y5+2*y6);
243  J[3*2+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
244 
245  J[3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
246  J[3*1+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*zi+j_const_reciprocal*(-z1-z2+2*z3-z4-z5+2*z6);
247  J[3*2+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
248 }
249 /*}}}*/
250 void PentaRef::GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
251  /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take
252  * the determinant of it: */
253  IssmDouble J[3][3];
254 
255  /*Get Jacobian*/
256  GetJacobian(&J[0][0],xyz_list,gauss);
257 
258  /*Get Determinant*/
259  Matrix3x3Determinant(Jdet,&J[0][0]);
260  if(*Jdet<0) _error_("negative jacobian determinant!");
261 
262 }
263 /*}}}*/
264 void PentaRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
265 
266  /*Jacobian*/
267  IssmDouble J[3][3];
268 
269  /*Call Jacobian routine to get the jacobian:*/
270  GetJacobian(&J[0][0], xyz_list, gauss);
271 
272  /*Invert Jacobian matrix: */
273  Matrix3x3Invert(Jinv,&J[0][0]);
274 }
275 /*}}}*/
276 void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement){/*{{{*/
277  /*This routine returns the values of the nodal functions at the gaussian point.*/
278 
279  _assert_(basis);
280 
281  /*Cast gauss to GaussPenta*/
282  _assert_(gauss_in->Enum()==GaussPentaEnum);
283  GaussPenta* gauss = xDynamicCast<GaussPenta*>(gauss_in);
284 
285  /*Get current coordinates in reference element*/
286  IssmDouble zeta=gauss->coord4;
287 
288  switch(finiteelement){
289  case P0Enum:
290  basis[0]=1.;
291  return;
292  case P1Enum: case P1DGEnum:
293  basis[0]=gauss->coord1*(1.-zeta)/2.;
294  basis[1]=gauss->coord2*(1.-zeta)/2.;
295  basis[2]=gauss->coord3*(1.-zeta)/2.;
296  basis[3]=gauss->coord1*(1.+zeta)/2.;
297  basis[4]=gauss->coord2*(1.+zeta)/2.;
298  basis[5]=gauss->coord3*(1.+zeta)/2.;
299  return;
300  case P1bubbleEnum: case P1bubblecondensedEnum:
301  basis[0]=gauss->coord1*(1.-zeta)/2.;
302  basis[1]=gauss->coord2*(1.-zeta)/2.;
303  basis[2]=gauss->coord3*(1.-zeta)/2.;
304  basis[3]=gauss->coord1*(1.+zeta)/2.;
305  basis[4]=gauss->coord2*(1.+zeta)/2.;
306  basis[5]=gauss->coord3*(1.+zeta)/2.;
307  basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta);
308  return;
309  case P2xP1Enum:
310  /*Corner nodes*/
311  basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta)/2.;
312  basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta)/2.;
313  basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta)/2.;
314  basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*(1.+zeta)/2.;
315  basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*(1.+zeta)/2.;
316  basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*(1.+zeta)/2.;
317  /*mid-sides of triangles*/
318  basis[ 6]=4.*gauss->coord3*gauss->coord2*(1.-zeta)/2.;
319  basis[ 7]=4.*gauss->coord3*gauss->coord1*(1.-zeta)/2.;
320  basis[ 8]=4.*gauss->coord1*gauss->coord2*(1.-zeta)/2.;
321  basis[ 9]=4.*gauss->coord3*gauss->coord2*(1.+zeta)/2.;
322  basis[10]=4.*gauss->coord3*gauss->coord1*(1.+zeta)/2.;
323  basis[11]=4.*gauss->coord1*gauss->coord2*(1.+zeta)/2.;
324  return;
325  case P1xP2Enum:
326  /*Corner nodes*/
327  basis[ 0]=gauss->coord1*zeta*(zeta-1.)/2.;
328  basis[ 1]=gauss->coord2*zeta*(zeta-1.)/2.;
329  basis[ 2]=gauss->coord3*zeta*(zeta-1.)/2.;
330  basis[ 3]=gauss->coord1*zeta*(zeta+1.)/2.;
331  basis[ 4]=gauss->coord2*zeta*(zeta+1.)/2.;
332  basis[ 5]=gauss->coord3*zeta*(zeta+1.)/2.;
333  /*mid-sides of quads*/
334  basis[ 6]=gauss->coord1*(1.-zeta*zeta);
335  basis[ 7]=gauss->coord2*(1.-zeta*zeta);
336  basis[ 8]=gauss->coord3*(1.-zeta*zeta);
337  return;
338  case P2Enum:
339  /*Corner nodes*/
340  basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta-1.)/2.;
341  basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta-1.)/2.;
342  basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta-1.)/2.;
343  basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta+1.)/2.;
344  basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta+1.)/2.;
345  basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.;
346  /*mid-sides of quads*/
347  basis[ 6]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta*zeta);
348  basis[ 7]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta*zeta);
349  basis[ 8]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta*zeta);
350  /*mid-sides of triangles*/
351  basis[ 9]=4.*gauss->coord3*gauss->coord2*zeta*(zeta-1.)/2.;
352  basis[10]=4.*gauss->coord3*gauss->coord1*zeta*(zeta-1.)/2.;
353  basis[11]=4.*gauss->coord1*gauss->coord2*zeta*(zeta-1.)/2.;
354  basis[12]=4.*gauss->coord3*gauss->coord2*zeta*(zeta+1.)/2.;
355  basis[13]=4.*gauss->coord3*gauss->coord1*zeta*(zeta+1.)/2.;
356  basis[14]=4.*gauss->coord1*gauss->coord2*zeta*(zeta+1.)/2.;
357  /*quad faces*/
358  basis[15]=4.*gauss->coord3*gauss->coord2*(1.-zeta*zeta);
359  basis[16]=4.*gauss->coord3*gauss->coord1*(1.-zeta*zeta);
360  basis[17]=4.*gauss->coord1*gauss->coord2*(1.-zeta*zeta);
361  return;
362  case P2bubbleEnum: case P2bubblecondensedEnum:
363  /*Corner nodes*/
364  basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta-1.)/2.;
365  basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta-1.)/2.;
366  basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta-1.)/2.;
367  basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta+1.)/2.;
368  basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta+1.)/2.;
369  basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.;
370  /*mid-sides of quads*/
371  basis[ 6]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta*zeta);
372  basis[ 7]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta*zeta);
373  basis[ 8]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta*zeta);
374  /*mid-sides of triangles*/
375  basis[ 9]=4.*gauss->coord3*gauss->coord2*zeta*(zeta-1.)/2.;
376  basis[10]=4.*gauss->coord3*gauss->coord1*zeta*(zeta-1.)/2.;
377  basis[11]=4.*gauss->coord1*gauss->coord2*zeta*(zeta-1.)/2.;
378  basis[12]=4.*gauss->coord3*gauss->coord2*zeta*(zeta+1.)/2.;
379  basis[13]=4.*gauss->coord3*gauss->coord1*zeta*(zeta+1.)/2.;
380  basis[14]=4.*gauss->coord1*gauss->coord2*zeta*(zeta+1.)/2.;
381  /*quad faces*/
382  basis[15]=4.*gauss->coord3*gauss->coord2*(1.-zeta*zeta);
383  basis[16]=4.*gauss->coord3*gauss->coord1*(1.-zeta*zeta);
384  basis[17]=4.*gauss->coord1*gauss->coord2*(1.-zeta*zeta);
385  /*bubble*/
386  basis[18]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta);
387  return;
388  case P2xP4Enum :
389  /*Corner nodes*/
390  basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
391  basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
392  basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
393  basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
394  basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
395  basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
396  /*mid-sides of quads*/
397  basis[ 6]=gauss->coord1*(2.*gauss->coord1-1)*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
398  basis[ 7]=gauss->coord2*(2.*gauss->coord2-1)*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
399  basis[ 8]=gauss->coord3*(2.*gauss->coord3-1)*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
400  /*mid-sides of triangles*/
401  basis[ 9]=4.*gauss->coord2*gauss->coord3*(2./3.)*(zeta-1.)*(zeta-0.5)*zeta*(zeta+0.5);
402  basis[10]=4.*gauss->coord1*gauss->coord3*(2./3.)*(zeta-1.)*(zeta-0.5)*zeta*(zeta+0.5);
403  basis[11]=4.*gauss->coord1*gauss->coord2*(2./3.)*(zeta-1.)*(zeta-0.5)*zeta*(zeta+0.5);
404  basis[12]=4.*gauss->coord2*gauss->coord3*(2./3.)*(zeta-0.5)*zeta*(zeta+0.5)*(zeta+1.);
405  basis[13]=4.*gauss->coord1*gauss->coord3*(2./3.)*(zeta-0.5)*zeta*(zeta+0.5)*(zeta+1.);
406  basis[14]=4.*gauss->coord1*gauss->coord2*(2./3.)*(zeta-0.5)*zeta*(zeta+0.5)*(zeta+1.);
407  /*quarter-sides of quads*/
408  basis[15]=gauss->coord1*(2.*gauss->coord1-1.)*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
409  basis[16]=gauss->coord2*(2.*gauss->coord2-1.)*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
410  basis[17]=gauss->coord3*(2.*gauss->coord3-1.)*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
411  basis[18]=gauss->coord1*(2.*gauss->coord1-1.)*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
412  basis[19]=gauss->coord2*(2.*gauss->coord2-1.)*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
413  basis[20]=gauss->coord3*(2.*gauss->coord3-1.)*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
414  /* mid-sides of interior triangles*/
415  basis[21]=4.*gauss->coord2*gauss->coord3*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
416  basis[22]=4.*gauss->coord1*gauss->coord3*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
417  basis[23]=4.*gauss->coord1*gauss->coord2*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
418  basis[24]=4.*gauss->coord2*gauss->coord3*4.*(zeta-1.0)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
419  basis[25]=4.*gauss->coord1*gauss->coord3*4.*(zeta-1.0)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
420  basis[26]=4.*gauss->coord1*gauss->coord2*4.*(zeta-1.0)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
421  basis[27]=4.*gauss->coord2*gauss->coord3*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
422  basis[28]=4.*gauss->coord1*gauss->coord3*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
423  basis[29]=4.*gauss->coord1*gauss->coord2*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
424  return;
425  case P1xP3Enum :
426  /*Corner nodes*/
427  basis[ 0]=-(9.)/(16.)*gauss->coord1*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
428  basis[ 1]=-(9.)/(16.)*gauss->coord2*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
429  basis[ 2]=-(9.)/(16.)*gauss->coord3*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
430  basis[ 3]=(9.)/(16.)*gauss->coord1*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
431  basis[ 4]=(9.)/(16.)*gauss->coord2*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
432  basis[ 5]=(9.)/(16.)*gauss->coord3*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
433  /*third-sides of quads*/
434  basis[ 6]=(27.)/(16.)*gauss->coord1*(zeta-1)*(zeta-1./3.)*(zeta+1.);
435  basis[ 7]=(27.)/(16.)*gauss->coord2*(zeta-1)*(zeta-1./3.)*(zeta+1.);
436  basis[ 8]=(27.)/(16.)*gauss->coord3*(zeta-1)*(zeta-1./3.)*(zeta+1.);
437  basis[ 9]=-(27.)/(16.)*gauss->coord1*(zeta-1)*(zeta+1./3.)*(zeta+1.);
438  basis[10]=-(27.)/(16.)*gauss->coord2*(zeta-1)*(zeta+1./3.)*(zeta+1.);
439  basis[11]=-(27.)/(16.)*gauss->coord3*(zeta-1)*(zeta+1./3.)*(zeta+1.);
440  return;
441  case P1xP4Enum :
442  /*Corner nodes*/
443  basis[ 0]=gauss->coord1*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
444  basis[ 1]=gauss->coord2*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
445  basis[ 2]=gauss->coord3*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
446  basis[ 3]=gauss->coord1*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
447  basis[ 4]=gauss->coord2*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
448  basis[ 5]=gauss->coord3*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
449  /*mid-sides of quads (center of vertical edges)*/
450  basis[ 6]=gauss->coord1*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
451  basis[ 7]=gauss->coord2*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
452  basis[ 8]=gauss->coord3*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
453  /*quarter-sides of quads (-0.5 and +0.5 of vertical edges)*/
454  basis[ 9]=gauss->coord1*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
455  basis[10]=gauss->coord2*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
456  basis[11]=gauss->coord3*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
457  basis[12]=gauss->coord1*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
458  basis[13]=gauss->coord2*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
459  basis[14]=gauss->coord3*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
460  return;
461  default:
462  _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
463  }
464 }
465 /*}}}*/
466 void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
467 
468  /*This routine returns the values of the nodal functions derivatives (with respect to the
469  * actual coordinate system): */
470  IssmDouble Jinv[3][3];
471 
472  /*Fetch number of nodes for this finite element*/
473  int numnodes = this->NumberofNodes(finiteelement);
474 
475  /*Get nodal functions derivatives in reference triangle*/
476  IssmDouble dbasis_ref[3*NUMNODESMAX];
477  GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
478 
479  /*Get Jacobian invert: */
480  GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
481 
482  /*Build dbasis:
483  *
484  * [dhi/dx]= Jinv'*[dhi/dr]
485  * [dhi/dy] [dhi/ds]
486  * [dhi/dz] [dhi/dzeta]
487  */
488 
489  for(int i=0;i<numnodes;i++){
490  dbasis[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i]+Jinv[0][2]*dbasis_ref[2*numnodes+i];
491  dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i]+Jinv[1][2]*dbasis_ref[2*numnodes+i];
492  dbasis[numnodes*2+i]=Jinv[2][0]*dbasis_ref[0*numnodes+i]+Jinv[2][1]*dbasis_ref[1*numnodes+i]+Jinv[2][2]*dbasis_ref[2*numnodes+i];
493  }
494 }
495 /*}}}*/
496 void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss_in,int finiteelement){/*{{{*/
497 
498  /*This routine returns the values of the nodal functions derivatives (with respect to the
499  * natural coordinate system) at the gaussian point. */
500 
501  _assert_(dbasis && gauss_in);
502 
503  /*Cast gauss to GaussPenta*/
504  _assert_(gauss_in->Enum()==GaussPentaEnum);
505  GaussPenta* gauss = xDynamicCast<GaussPenta*>(gauss_in);
506 
507  /*Get current coordinates in reference element*/
508  IssmDouble zeta=gauss->coord4;
509 
510  switch(finiteelement){
511  case P0Enum:
512  /*Zero derivative*/
513  dbasis[NUMNODESP0*0+0] = 0.;
514  dbasis[NUMNODESP0*1+0] = 0.;
515  dbasis[NUMNODESP0*2+0] = 0.;
516  return;
517  case P1Enum: case P1DGEnum:
518  /*Nodal function 1*/
519  dbasis[NUMNODESP1*0+0] = (zeta-1.)/4.;
520  dbasis[NUMNODESP1*1+0] = SQRT3/12.*(zeta-1.);
521  dbasis[NUMNODESP1*2+0] = -.5*gauss->coord1;
522  /*Nodal function 2*/
523  dbasis[NUMNODESP1*0+1] = (1.-zeta)/4.;
524  dbasis[NUMNODESP1*1+1] = SQRT3/12.*(zeta-1);
525  dbasis[NUMNODESP1*2+1] = -.5*gauss->coord2;
526  /*Nodal function 3*/
527  dbasis[NUMNODESP1*0+2] = 0.;
528  dbasis[NUMNODESP1*1+2] = SQRT3/6.*(1.-zeta);
529  dbasis[NUMNODESP1*2+2] = -.5*gauss->coord3;
530  /*Nodal function 4*/
531  dbasis[NUMNODESP1*0+3] = -(1.+zeta)/4.;
532  dbasis[NUMNODESP1*1+3] = -SQRT3/12.*(1.+zeta);
533  dbasis[NUMNODESP1*2+3] = .5*gauss->coord1;
534  /*Nodal function 5*/
535  dbasis[NUMNODESP1*0+4] = (1.+zeta)/4.;
536  dbasis[NUMNODESP1*1+4] = -SQRT3/12.*(1.+zeta);
537  dbasis[NUMNODESP1*2+4] = .5*gauss->coord2;
538  /*Nodal function 6*/
539  dbasis[NUMNODESP1*0+5] = 0.;
540  dbasis[NUMNODESP1*1+5] = SQRT3/6.*(1.+zeta);
541  dbasis[NUMNODESP1*2+5] = .5*gauss->coord3;
542  return;
543  case P1bubbleEnum: case P1bubblecondensedEnum:
544  /*Nodal function 1*/
545  dbasis[NUMNODESP1b*0+0] = (zeta-1.)/4.;
546  dbasis[NUMNODESP1b*1+0] = SQRT3/12.*(zeta-1.);
547  dbasis[NUMNODESP1b*2+0] = -.5*gauss->coord1;
548  /*Nodal function 2*/
549  dbasis[NUMNODESP1b*0+1] = (1.-zeta)/4.;
550  dbasis[NUMNODESP1b*1+1] = SQRT3/12.*(zeta-1);
551  dbasis[NUMNODESP1b*2+1] = -.5*gauss->coord2;
552  /*Nodal function 3*/
553  dbasis[NUMNODESP1b*0+2] = 0.;
554  dbasis[NUMNODESP1b*1+2] = SQRT3/6.*(1.-zeta);
555  dbasis[NUMNODESP1b*2+2] = -.5*gauss->coord3;
556  /*Nodal function 4*/
557  dbasis[NUMNODESP1b*0+3] = -(1.+zeta)/4.;
558  dbasis[NUMNODESP1b*1+3] = -SQRT3/12.*(1.+zeta);
559  dbasis[NUMNODESP1b*2+3] = .5*gauss->coord1;
560  /*Nodal function 5*/
561  dbasis[NUMNODESP1b*0+4] = (1.+zeta)/4.;
562  dbasis[NUMNODESP1b*1+4] = -SQRT3/12.*(1.+zeta);
563  dbasis[NUMNODESP1b*2+4] = .5*gauss->coord2;
564  /*Nodal function 6*/
565  dbasis[NUMNODESP1b*0+5] = 0.;
566  dbasis[NUMNODESP1b*1+5] = SQRT3/6.*(1.+zeta);
567  dbasis[NUMNODESP1b*2+5] = .5*gauss->coord3;
568  /*Nodal function 7*/
569  dbasis[NUMNODESP1b*0+6] = 27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
570  dbasis[NUMNODESP1b*1+6] = 27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
571  dbasis[NUMNODESP1b*2+6] = -54*gauss->coord1*gauss->coord2*gauss->coord3*zeta;
572  return;
573  case P2xP1Enum:
574  /*Nodal function 1*/
575  dbasis[NUMNODESP2xP1*0+0 ] = .5*(1.-zeta)*(-2.*gauss->coord1 + 0.5);
576  dbasis[NUMNODESP2xP1*1+0 ] = .5*(1.-zeta)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
577  dbasis[NUMNODESP2xP1*2+0 ] = -.5*gauss->coord1*(2.*gauss->coord1-1.);
578  /*Nodal function 2*/
579  dbasis[NUMNODESP2xP1*0+1 ] = .5*(1.-zeta)*(+2.*gauss->coord2 - 0.5);
580  dbasis[NUMNODESP2xP1*1+1 ] = .5*(1.-zeta)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
581  dbasis[NUMNODESP2xP1*2+1 ] = -.5*gauss->coord2*(2.*gauss->coord2-1.);
582  /*Nodal function 3*/
583  dbasis[NUMNODESP2xP1*0+2 ] = 0.;
584  dbasis[NUMNODESP2xP1*1+2 ] = .5*(1.-zeta)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
585  dbasis[NUMNODESP2xP1*2+2 ] = -.5*gauss->coord3*(2.*gauss->coord3-1.);
586  /*Nodal function 4*/
587  dbasis[NUMNODESP2xP1*0+3 ] = .5*(1.+zeta)*(-2.*gauss->coord1 + 0.5);
588  dbasis[NUMNODESP2xP1*1+3 ] = .5*(1.+zeta)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
589  dbasis[NUMNODESP2xP1*2+3 ] = .5*gauss->coord1*(2.*gauss->coord1-1.);
590  /*Nodal function 5*/
591  dbasis[NUMNODESP2xP1*0+4 ] = .5*(1.+zeta)*(+2.*gauss->coord2 - 0.5);
592  dbasis[NUMNODESP2xP1*1+4 ] = .5*(1.+zeta)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
593  dbasis[NUMNODESP2xP1*2+4 ] = .5*gauss->coord2*(2.*gauss->coord2-1.);
594  /*Nodal function 6*/
595  dbasis[NUMNODESP2xP1*0+5 ] = 0.;
596  dbasis[NUMNODESP2xP1*1+5 ] = .5*(1.+zeta)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
597  dbasis[NUMNODESP2xP1*2+5 ] = .5*gauss->coord3*(2.*gauss->coord3-1.);
598 
599  /*Nodal function 7*/
600  dbasis[NUMNODESP2xP1*0+6 ] = (1.-zeta)*gauss->coord3;
601  dbasis[NUMNODESP2xP1*1+6 ] = .5*(1.-zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
602  dbasis[NUMNODESP2xP1*2+6 ] = -2.*gauss->coord3*gauss->coord2;
603  /*Nodal function 8*/
604  dbasis[NUMNODESP2xP1*0+7 ] = -(1.-zeta)*gauss->coord3;
605  dbasis[NUMNODESP2xP1*1+7 ] = .5*(1.-zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
606  dbasis[NUMNODESP2xP1*2+7 ] = -2.*gauss->coord3*gauss->coord1;
607  /*Nodal function 9*/
608  dbasis[NUMNODESP2xP1*0+8 ] = (1.-zeta)*(gauss->coord1-gauss->coord2);
609  dbasis[NUMNODESP2xP1*1+8 ] = .5*(1.-zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
610  dbasis[NUMNODESP2xP1*2+8 ] = -2.*gauss->coord1*gauss->coord2;
611  /*Nodal function 10*/
612  dbasis[NUMNODESP2xP1*0+9 ] = (1.+zeta)*gauss->coord3;
613  dbasis[NUMNODESP2xP1*1+9 ] = .5*(1.+zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
614  dbasis[NUMNODESP2xP1*2+9 ] = 2.*gauss->coord3*gauss->coord2;
615  /*Nodal function 11*/
616  dbasis[NUMNODESP2xP1*0+10] = -(1.+zeta)*gauss->coord3;
617  dbasis[NUMNODESP2xP1*1+10] = .5*(1.+zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
618  dbasis[NUMNODESP2xP1*2+10] = 2.*gauss->coord3*gauss->coord1;
619  /*Nodal function 12*/
620  dbasis[NUMNODESP2xP1*0+11] = (1.+zeta)*(gauss->coord1-gauss->coord2);
621  dbasis[NUMNODESP2xP1*1+11] = .5*(1.+zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
622  dbasis[NUMNODESP2xP1*2+11] = 2.*gauss->coord1*gauss->coord2;
623  return;
624  case P1xP2Enum:
625  /*Nodal function 1*/
626  dbasis[NUMNODESP1xP2*0+0] = -zeta*(zeta-1.)/4.;
627  dbasis[NUMNODESP1xP2*1+0] = -SQRT3/12.*zeta*(zeta-1.);
628  dbasis[NUMNODESP1xP2*2+0] = .5*(2.*zeta-1.)*gauss->coord1;
629  /*Nodal function 2*/
630  dbasis[NUMNODESP1xP2*0+1] = zeta*(zeta-1.)/4.;
631  dbasis[NUMNODESP1xP2*1+1] = -SQRT3/12.*zeta*(zeta-1);
632  dbasis[NUMNODESP1xP2*2+1] = .5*(2.*zeta-1.)*gauss->coord2;
633  /*Nodal function 3*/
634  dbasis[NUMNODESP1xP2*0+2] = 0.;
635  dbasis[NUMNODESP1xP2*1+2] = SQRT3/6.*zeta*(zeta-1.);
636  dbasis[NUMNODESP1xP2*2+2] = .5*(2.*zeta-1.)*gauss->coord3;
637  /*Nodal function 4*/
638  dbasis[NUMNODESP1xP2*0+3] = -zeta*(zeta+1)/4.;
639  dbasis[NUMNODESP1xP2*1+3] = -SQRT3/12.*zeta*(zeta+1.);
640  dbasis[NUMNODESP1xP2*2+3] = .5*(2.*zeta+1.)*gauss->coord1;
641  /*Nodal function 5*/
642  dbasis[NUMNODESP1xP2*0+4] = zeta*(zeta+1.)/4.;
643  dbasis[NUMNODESP1xP2*1+4] = -SQRT3/12.*zeta*(zeta+1.);
644  dbasis[NUMNODESP1xP2*2+4] = .5*(2.*zeta+1.)*gauss->coord2;
645  /*Nodal function 6*/
646  dbasis[NUMNODESP1xP2*0+5] = 0.;
647  dbasis[NUMNODESP1xP2*1+5] = SQRT3/6.*zeta*(zeta+1.);
648  dbasis[NUMNODESP1xP2*2+5] = .5*(2.*zeta+1.)*gauss->coord3;
649 
650  /*Nodal function 7*/
651  dbasis[NUMNODESP1xP2*0+6 ] = -0.5*(1.-zeta*zeta);
652  dbasis[NUMNODESP1xP2*1+6 ] = -SQRT3/6.*(1.-zeta*zeta);
653  dbasis[NUMNODESP1xP2*2+6 ] = -2.*zeta*gauss->coord1;
654  /*Nodal function 8*/
655  dbasis[NUMNODESP1xP2*0+7 ] = 0.5*(1.-zeta*zeta);
656  dbasis[NUMNODESP1xP2*1+7 ] = -SQRT3/6.*(1.-zeta*zeta);
657  dbasis[NUMNODESP1xP2*2+7 ] = -2.*zeta*gauss->coord2;
658  /*Nodal function 9*/
659  dbasis[NUMNODESP1xP2*0+8 ] = 0.;
660  dbasis[NUMNODESP1xP2*1+8 ] = SQRT3/3.*(1.-zeta*zeta);
661  dbasis[NUMNODESP1xP2*2+8 ] = -2.*zeta*gauss->coord3;
662  return;
663  case P2Enum:
664  /*Nodal function 1*/
665  dbasis[NUMNODESP2*0+0 ] = .5*zeta*(zeta-1.)*(-2.*gauss->coord1 + 0.5);
666  dbasis[NUMNODESP2*1+0 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
667  dbasis[NUMNODESP2*2+0 ] = .5*(2.*zeta-1.)*gauss->coord1*(2.*gauss->coord1-1.);
668  /*Nodal function 2*/
669  dbasis[NUMNODESP2*0+1 ] = .5*zeta*(zeta-1.)*(+2.*gauss->coord2 - 0.5);
670  dbasis[NUMNODESP2*1+1 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
671  dbasis[NUMNODESP2*2+1 ] = .5*(2.*zeta-1.)*gauss->coord2*(2.*gauss->coord2-1.);
672  /*Nodal function 3*/
673  dbasis[NUMNODESP2*0+2 ] = 0.;
674  dbasis[NUMNODESP2*1+2 ] = .5*zeta*(zeta-1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
675  dbasis[NUMNODESP2*2+2 ] = .5*(2.*zeta-1.)*gauss->coord3*(2.*gauss->coord3-1.);
676  /*Nodal function 4*/
677  dbasis[NUMNODESP2*0+3 ] = .5*zeta*(zeta+1.)*(-2.*gauss->coord1 + 0.5);
678  dbasis[NUMNODESP2*1+3 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
679  dbasis[NUMNODESP2*2+3 ] = .5*(2.*zeta+1.)*gauss->coord1*(2.*gauss->coord1-1.);
680  /*Nodal function 5*/
681  dbasis[NUMNODESP2*0+4 ] = .5*zeta*(zeta+1.)*(+2.*gauss->coord2 - 0.5);
682  dbasis[NUMNODESP2*1+4 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
683  dbasis[NUMNODESP2*2+4 ] = .5*(2.*zeta+1.)*gauss->coord2*(2.*gauss->coord2-1.);
684  /*Nodal function 6*/
685  dbasis[NUMNODESP2*0+5 ] = 0.;
686  dbasis[NUMNODESP2*1+5 ] = .5*zeta*(zeta+1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
687  dbasis[NUMNODESP2*2+5 ] = .5*(2.*zeta+1.)*gauss->coord3*(2.*gauss->coord3-1.);
688 
689  /*Nodal function 7*/
690  dbasis[NUMNODESP2*0+6 ] = (-2.*gauss->coord1 + 0.5)*(1.-zeta*zeta);
691  dbasis[NUMNODESP2*1+6 ] = (-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.)*(1.-zeta*zeta);
692  dbasis[NUMNODESP2*2+6 ] = -2.*zeta*gauss->coord1*(2.*gauss->coord1-1.);
693  /*Nodal function 8*/
694  dbasis[NUMNODESP2*0+7 ] = (+2.*gauss->coord2 - 0.5)*(1.-zeta*zeta);
695  dbasis[NUMNODESP2*1+7 ] = (-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.)*(1.-zeta*zeta);
696  dbasis[NUMNODESP2*2+7 ] = -2.*zeta*gauss->coord2*(2.*gauss->coord2-1.);
697  /*Nodal function 9*/
698  dbasis[NUMNODESP2*0+8 ] = 0.;
699  dbasis[NUMNODESP2*1+8 ] = (+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.)*(1.-zeta*zeta);
700  dbasis[NUMNODESP2*2+8 ] = -2.*zeta*gauss->coord3*(2.*gauss->coord3-1.);
701 
702  /*Nodal function 10*/
703  dbasis[NUMNODESP2*0+9 ] = zeta*(zeta-1.)*gauss->coord3;
704  dbasis[NUMNODESP2*1+9 ] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
705  dbasis[NUMNODESP2*2+9 ] = 2.*gauss->coord3*gauss->coord2*(2.*zeta-1.);
706  /*Nodal function 11*/
707  dbasis[NUMNODESP2*0+10] = -zeta*(zeta-1.)*gauss->coord3;
708  dbasis[NUMNODESP2*1+10] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
709  dbasis[NUMNODESP2*2+10] = 2.*gauss->coord3*gauss->coord1*(2.*zeta-1.);
710  /*Nodal function 12*/
711  dbasis[NUMNODESP2*0+11] = zeta*(zeta-1.)*(gauss->coord1-gauss->coord2);
712  dbasis[NUMNODESP2*1+11] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
713  dbasis[NUMNODESP2*2+11] = 2.*gauss->coord1*gauss->coord2*(2.*zeta-1.);
714  /*Nodal function 13*/
715  dbasis[NUMNODESP2*0+12] = zeta*(zeta+1.)*gauss->coord3;
716  dbasis[NUMNODESP2*1+12] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
717  dbasis[NUMNODESP2*2+12] = 2.*gauss->coord3*gauss->coord2*(2.*zeta+1.);
718  /*Nodal function 14*/
719  dbasis[NUMNODESP2*0+13] = -zeta*(zeta+1.)*gauss->coord3;
720  dbasis[NUMNODESP2*1+13] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
721  dbasis[NUMNODESP2*2+13] = 2.*gauss->coord3*gauss->coord1*(2.*zeta+1.);
722  /*Nodal function 15*/
723  dbasis[NUMNODESP2*0+14] = zeta*(zeta+1.)*(gauss->coord1-gauss->coord2);
724  dbasis[NUMNODESP2*1+14] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
725  dbasis[NUMNODESP2*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.);
726 
727  /*Nodal function 16*/
728  dbasis[NUMNODESP2*0+15] = 2.*gauss->coord3*(1.-zeta*zeta);
729  dbasis[NUMNODESP2*1+15] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
730  dbasis[NUMNODESP2*2+15] = -2.*zeta*4.*gauss->coord3*gauss->coord2;
731  /*Nodal function 17*/
732  dbasis[NUMNODESP2*0+16] = -2.*gauss->coord3*(1.-zeta*zeta);
733  dbasis[NUMNODESP2*1+16] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
734  dbasis[NUMNODESP2*2+16] = -2.*zeta*4.*gauss->coord3*gauss->coord1;
735  /*Nodal function 18*/
736  dbasis[NUMNODESP2*0+17] = 2.*(gauss->coord1-gauss->coord2)*(1.-zeta*zeta);
737  dbasis[NUMNODESP2*1+17] = (1.-zeta*zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
738  dbasis[NUMNODESP2*2+17] = -2.*zeta*4.*gauss->coord1*gauss->coord2;
739  return;
740  case P2bubbleEnum: case P2bubblecondensedEnum:
741  /*Nodal function 1*/
742  dbasis[NUMNODESP2b*0+0 ] = .5*zeta*(zeta-1.)*(-2.*gauss->coord1 + 0.5);
743  dbasis[NUMNODESP2b*1+0 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
744  dbasis[NUMNODESP2b*2+0 ] = .5*(2.*zeta-1.)*gauss->coord1*(2.*gauss->coord1-1.);
745  /*Nodal function 2*/
746  dbasis[NUMNODESP2b*0+1 ] = .5*zeta*(zeta-1.)*(+2.*gauss->coord2 - 0.5);
747  dbasis[NUMNODESP2b*1+1 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
748  dbasis[NUMNODESP2b*2+1 ] = .5*(2.*zeta-1.)*gauss->coord2*(2.*gauss->coord2-1.);
749  /*Nodal function 3*/
750  dbasis[NUMNODESP2b*0+2 ] = 0.;
751  dbasis[NUMNODESP2b*1+2 ] = .5*zeta*(zeta-1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
752  dbasis[NUMNODESP2b*2+2 ] = .5*(2.*zeta-1.)*gauss->coord3*(2.*gauss->coord3-1.);
753  /*Nodal function 4*/
754  dbasis[NUMNODESP2b*0+3 ] = .5*zeta*(zeta+1.)*(-2.*gauss->coord1 + 0.5);
755  dbasis[NUMNODESP2b*1+3 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
756  dbasis[NUMNODESP2b*2+3 ] = .5*(2.*zeta+1.)*gauss->coord1*(2.*gauss->coord1-1.);
757  /*Nodal function 5*/
758  dbasis[NUMNODESP2b*0+4 ] = .5*zeta*(zeta+1.)*(+2.*gauss->coord2 - 0.5);
759  dbasis[NUMNODESP2b*1+4 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
760  dbasis[NUMNODESP2b*2+4 ] = .5*(2.*zeta+1.)*gauss->coord2*(2.*gauss->coord2-1.);
761  /*Nodal function 6*/
762  dbasis[NUMNODESP2b*0+5 ] = 0.;
763  dbasis[NUMNODESP2b*1+5 ] = .5*zeta*(zeta+1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
764  dbasis[NUMNODESP2b*2+5 ] = .5*(2.*zeta+1.)*gauss->coord3*(2.*gauss->coord3-1.);
765 
766  /*Nodal function 7*/
767  dbasis[NUMNODESP2b*0+6 ] = (-2.*gauss->coord1 + 0.5)*(1.-zeta*zeta);
768  dbasis[NUMNODESP2b*1+6 ] = (-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.)*(1.-zeta*zeta);
769  dbasis[NUMNODESP2b*2+6 ] = -2.*zeta*gauss->coord1*(2.*gauss->coord1-1.);
770  /*Nodal function 8*/
771  dbasis[NUMNODESP2b*0+7 ] = (+2.*gauss->coord2 - 0.5)*(1.-zeta*zeta);
772  dbasis[NUMNODESP2b*1+7 ] = (-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.)*(1.-zeta*zeta);
773  dbasis[NUMNODESP2b*2+7 ] = -2.*zeta*gauss->coord2*(2.*gauss->coord2-1.);
774  /*Nodal function 9*/
775  dbasis[NUMNODESP2b*0+8 ] = 0.;
776  dbasis[NUMNODESP2b*1+8 ] = (+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.)*(1.-zeta*zeta);
777  dbasis[NUMNODESP2b*2+8 ] = -2.*zeta*gauss->coord3*(2.*gauss->coord3-1.);
778 
779  /*Nodal function 10*/
780  dbasis[NUMNODESP2b*0+9 ] = zeta*(zeta-1.)*gauss->coord3;
781  dbasis[NUMNODESP2b*1+9 ] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
782  dbasis[NUMNODESP2b*2+9 ] = 2.*gauss->coord3*gauss->coord2*(2.*zeta-1.);
783  /*Nodal function 11*/
784  dbasis[NUMNODESP2b*0+10] = -zeta*(zeta-1.)*gauss->coord3;
785  dbasis[NUMNODESP2b*1+10] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
786  dbasis[NUMNODESP2b*2+10] = 2.*gauss->coord3*gauss->coord1*(2.*zeta-1.);
787  /*Nodal function 12*/
788  dbasis[NUMNODESP2b*0+11] = zeta*(zeta-1.)*(gauss->coord1-gauss->coord2);
789  dbasis[NUMNODESP2b*1+11] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
790  dbasis[NUMNODESP2b*2+11] = 2.*gauss->coord1*gauss->coord2*(2.*zeta-1.);
791  /*Nodal function 13*/
792  dbasis[NUMNODESP2b*0+12] = zeta*(zeta+1.)*gauss->coord3;
793  dbasis[NUMNODESP2b*1+12] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
794  dbasis[NUMNODESP2b*2+12] = 2.*gauss->coord3*gauss->coord2*(2.*zeta+1.);
795  /*Nodal function 14*/
796  dbasis[NUMNODESP2b*0+13] = -zeta*(zeta+1.)*gauss->coord3;
797  dbasis[NUMNODESP2b*1+13] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
798  dbasis[NUMNODESP2b*2+13] = 2.*gauss->coord3*gauss->coord1*(2.*zeta+1.);
799  /*Nodal function 15*/
800  dbasis[NUMNODESP2b*0+14] = zeta*(zeta+1.)*(gauss->coord1-gauss->coord2);
801  dbasis[NUMNODESP2b*1+14] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
802  dbasis[NUMNODESP2b*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.);
803 
804  /*Nodal function 16*/
805  dbasis[NUMNODESP2b*0+15] = 2.*gauss->coord3*(1.-zeta*zeta);
806  dbasis[NUMNODESP2b*1+15] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
807  dbasis[NUMNODESP2b*2+15] = -2.*zeta*4.*gauss->coord3*gauss->coord2;
808  /*Nodal function 17*/
809  dbasis[NUMNODESP2b*0+16] = -2.*gauss->coord3*(1.-zeta*zeta);
810  dbasis[NUMNODESP2b*1+16] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
811  dbasis[NUMNODESP2b*2+16] = -2.*zeta*4.*gauss->coord3*gauss->coord1;
812  /*Nodal function 18*/
813  dbasis[NUMNODESP2b*0+17] = 2.*(gauss->coord1-gauss->coord2)*(1.-zeta*zeta);
814  dbasis[NUMNODESP2b*1+17] = (1.-zeta*zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
815  dbasis[NUMNODESP2b*2+17] = -2.*zeta*4.*gauss->coord1*gauss->coord2;
816 
817  /*Nodal function 19*/
818  dbasis[NUMNODESP2b*0+18] = 27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
819  dbasis[NUMNODESP2b*1+18] = 27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
820  dbasis[NUMNODESP2b*2+18] = -54*gauss->coord1*gauss->coord2*gauss->coord3*zeta;
821  return;
822  case P2xP4Enum :
823  /*Nodal function 1*/
824  dbasis[NUMNODESP2xP4*0+0 ] = (-2* gauss->coord1 + 0.5 ) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
825  dbasis[NUMNODESP2xP4*1+0 ] = (-((2.*SQRT3)/(3.))*gauss->coord1 + (SQRT3/6.) )*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
826  dbasis[NUMNODESP2xP4*2+0 ] = gauss->coord1 *(2.* gauss->coord1 -1)* 2./3.*( (2.*zeta-1)*(zeta -0.5)*(zeta +0.5) + 2.* zeta *zeta *(zeta -1.));
827  /*Nodal function 2*/
828  dbasis[NUMNODESP2xP4*0+1 ] = (2.*gauss->coord2 - 0.5 ) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
829  dbasis[NUMNODESP2xP4*1+1 ] = (-((2.*SQRT3)/(3.))*gauss->coord2 + (SQRT3/6.) )*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
830  dbasis[NUMNODESP2xP4*2+1 ] = gauss->coord2*(2.*gauss->coord2 -1.)* 2./3.* ((2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2. * zeta *zeta*(zeta -1.));
831  /*Nodal function 3*/
832  dbasis[NUMNODESP2xP4*0+2 ] = 0. ;
833  dbasis[NUMNODESP2xP4*1+2 ] = (((4.*SQRT3)/(3.))*gauss->coord3 - (SQRT3)/(3.))*(2./3.)*(zeta -1.)*(zeta-0.5)*(zeta)*(zeta+0.5);
834  dbasis[NUMNODESP2xP4*2+2 ] = gauss->coord3*(2.* gauss->coord3 -1.)* 2./3.*( (2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta -1.));
835  /*Nodal function 4*/
836  dbasis[NUMNODESP2xP4*0+3 ] = (-2.* gauss->coord1 + 0.5 ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
837  dbasis[NUMNODESP2xP4*1+3 ] = (-((2.*SQRT3)/(3.)) *gauss->coord1 + (SQRT3)/(6.) ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
838  dbasis[NUMNODESP2xP4*2+3 ] = gauss->coord1*(2.*gauss->coord1 -1.)* 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.));
839  /*Nodal function 5*/
840  dbasis[NUMNODESP2xP4*0+4 ] = (2*gauss->coord2 - 0.5 ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
841  dbasis[NUMNODESP2xP4*1+4 ] = (-((2.*SQRT3)/(3.)) *gauss->coord2 + (SQRT3/6.))*(2./3.)*(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
842  dbasis[NUMNODESP2xP4*2+4 ] = gauss->coord2 *(2.*gauss->coord2 -1.)* 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.));
843  /*Nodal function 6*/
844  dbasis[NUMNODESP2xP4*0+5 ] = 0. ;
845  dbasis[NUMNODESP2xP4*1+5 ] = (((4.*SQRT3)/(3.))*gauss->coord3 - SQRT3/3. ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
846  dbasis[NUMNODESP2xP4*2+5 ] = gauss->coord3 *(2.*gauss->coord3 -1.)* 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1));
847  /*Nodal function 7*/
848  dbasis[NUMNODESP2xP4*0+6 ] = (-2.* gauss->coord1 + 0.5 ) * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
849  dbasis[NUMNODESP2xP4*1+6 ] = (-((2.*SQRT3)/(3.)) *gauss->coord1 + (SQRT3)/(6.) )* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1.) ;
850  dbasis[NUMNODESP2xP4*2+6 ] = gauss->coord1*(2.*gauss->coord1-1)* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta );
851  /*Nodal function 8*/
852  dbasis[NUMNODESP2xP4*0+7 ] = (2*gauss->coord2 - 0.5 )* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
853  dbasis[NUMNODESP2xP4*1+7 ] = (-((2.*SQRT3)/(3.)) *gauss->coord2 + (SQRT3)/(6.)) * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
854  dbasis[NUMNODESP2xP4*2+7 ] = gauss->coord2*(2.*gauss->coord2-1)* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta );
855  /*Nodal function 9*/
856  dbasis[NUMNODESP2xP4*0+8 ] = 0. ;
857  dbasis[NUMNODESP2xP4*1+8 ] = (((4.*SQRT3)/(3.))*gauss->coord3 - SQRT3/3. ) * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
858  dbasis[NUMNODESP2xP4*2+8 ] = gauss->coord3*(2.*gauss->coord3-1)* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta );
859  /*Nodal function 10*/
860  dbasis[NUMNODESP2xP4*0+9 ] = 2.*gauss->coord3 * 2./3.*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
861  dbasis[NUMNODESP2xP4*1+9 ] = (4.* SQRT3/3.* gauss->coord2- 2.*SQRT3/3. *gauss->coord3) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
862  dbasis[NUMNODESP2xP4*2+9 ] = 4.* gauss->coord2 * gauss->coord3 *(2./3.)*((2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2.* zeta*zeta*(zeta -1.));
863  /*Nodal function 11*/
864  dbasis[NUMNODESP2xP4*0+10] = -2.* gauss->coord3* 2./3.*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
865  dbasis[NUMNODESP2xP4*1+10] = (4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3) * (2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
866  dbasis[NUMNODESP2xP4*2+10] = 4.* gauss->coord3*gauss->coord1 *(2./3.)*( (2*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2* zeta * zeta*(zeta -1));
867  /*Nodal function 12*/
868  dbasis[NUMNODESP2xP4*0+11] = 2.* (gauss->coord1- gauss->coord2)* (2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
869  dbasis[NUMNODESP2xP4*1+11] = -2.* SQRT3/3.*(gauss->coord2 +gauss->coord1) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
870  dbasis[NUMNODESP2xP4*2+11] = 4.* gauss->coord1*gauss->coord2 *(2./3.) *( (2.*zeta-1)*(zeta -0.5)*(zeta +0.5) + 2* zeta* zeta*(zeta -1));
871  /*Nodal function 13*/
872  dbasis[NUMNODESP2xP4*0+12] = 2.* gauss->coord3 * 2./3.*(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
873  dbasis[NUMNODESP2xP4*1+12] = (4.*SQRT3/3.* gauss->coord2 - 2.*SQRT3/3. *gauss->coord3) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
874  dbasis[NUMNODESP2xP4*2+12] = 4.*gauss->coord2*gauss->coord3 *(2./3.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2*zeta*zeta*(zeta +1.));
875  /*Nodal function 14*/
876  dbasis[NUMNODESP2xP4*0+13] = - 2.*gauss->coord3* 2./3.*(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
877  dbasis[NUMNODESP2xP4*1+13] = (4.*SQRT3/3.*gauss->coord1- 2.*SQRT3/3.*gauss->coord3) * (2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
878  dbasis[NUMNODESP2xP4*2+13] = 4.*gauss->coord3*gauss->coord1 *(2./3.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.));
879  /*Nodal function 15*/
880  dbasis[NUMNODESP2xP4*0+14] = 2.* (gauss->coord1- gauss->coord2)* (2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
881  dbasis[NUMNODESP2xP4*1+14] = -2.* SQRT3/3.*(gauss->coord2 +gauss->coord1) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
882  dbasis[NUMNODESP2xP4*2+14] = 4.*gauss->coord1*gauss->coord2 *(2./3.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2* zeta *zeta*(zeta +1.));
883  /*Nodal function 16*/
884  dbasis[NUMNODESP2xP4*0+15] = (-2.* gauss->coord1 + 0.5 )* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
885  dbasis[NUMNODESP2xP4*1+15] = (-2.*SQRT3/3. *gauss->coord1 + SQRT3/6.) * (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
886  dbasis[NUMNODESP2xP4*2+15] = gauss->coord1*(2.*gauss->coord1-1) * (-8./3.)*((2.*zeta -1.)*(zeta-0.5)*(zeta +1.) +zeta*(zeta -1.)*( 2.*zeta + 0.5));
887  /*Nodal function 17*/
888  dbasis[NUMNODESP2xP4*0+16] = (2*gauss->coord2 - 0.5 )* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
889  dbasis[NUMNODESP2xP4*1+16] = (-2.*SQRT3/3. *gauss->coord2 + SQRT3/6.)* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
890  dbasis[NUMNODESP2xP4*2+16] = gauss->coord2*(2.*gauss->coord2-1) * (-8./3.)*((2.*zeta -1.)*(zeta-0.5)*(zeta +1.) +zeta *(zeta -1.)*( 2.*zeta + 0.5));
891  /*Nodal function 18*/
892  dbasis[NUMNODESP2xP4*0+17] = 0. ;
893  dbasis[NUMNODESP2xP4*1+17] = (4.*SQRT3/3.*gauss->coord3 - SQRT3/3. )* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
894  dbasis[NUMNODESP2xP4*2+17] = gauss->coord3*(2*gauss->coord3-1) * (-8./3.)*((2.*zeta-1.)*(zeta-0.5)*(zeta +1.) +zeta *(zeta -1.)*( 2.*zeta + 0.5));
895  /*Nodal function 19*/
896  dbasis[NUMNODESP2xP4*0+18] = (-2.* gauss->coord1 + 0.5 ) * (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
897  dbasis[NUMNODESP2xP4*1+18] = (-2.*SQRT3/3. *gauss->coord1 + SQRT3/6. )* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
898  dbasis[NUMNODESP2xP4*2+18] = gauss->coord1*(2.*gauss->coord1-1) * (-8./3.)*((2.*zeta -1. ) *(zeta+0.5)* (zeta +1.) + zeta* (zeta -1.)*( 2.*zeta + 3./2.));
899  /*Nodal function 20*/
900  dbasis[NUMNODESP2xP4*0+19] = (2*gauss->coord2 - 0.5 )* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
901  dbasis[NUMNODESP2xP4*1+19] = (-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.) * (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
902  dbasis[NUMNODESP2xP4*2+19] = gauss->coord2*(2.*gauss->coord2-1)* (-8./3.)*((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) + zeta*(zeta -1.)*( 2.*zeta + 3./2.));
903  /*Nodal function 21*/
904  dbasis[NUMNODESP2xP4*0+20] = 0 ;
905  dbasis[NUMNODESP2xP4*1+20] = (4.*SQRT3/3.*gauss->coord3 - SQRT3/3.)* (-8./3.)*(zeta - 1.)*(zeta + 0.5)*(zeta)*(zeta +1. ) ;
906  dbasis[NUMNODESP2xP4*2+20] = gauss->coord3*(2*gauss->coord3-1) * (-8./3.)*((2. *zeta -1. )*(zeta+0.5)*(zeta +1.) + zeta*(zeta -1.)*( 2.*zeta + 3./2.));
907  /*Nodal function 22*/
908  dbasis[NUMNODESP2xP4*0+21] = 2. *gauss->coord3 * (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. );
909  dbasis[NUMNODESP2xP4*1+21] = (4.* SQRT3/3.*gauss->coord2- 2.*SQRT3/3.*gauss->coord3) * (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
910  dbasis[NUMNODESP2xP4*2+21] = 4.*gauss->coord2 *gauss->coord3* (-8./3.)*((2.*zeta -1. )*(zeta-0.5)*(zeta +1.) + zeta*(zeta -1.)*( 2.*zeta + 0.5));
911  /*Nodal function 23*/
912  dbasis[NUMNODESP2xP4*0+22] = -2. *gauss->coord3 *( -8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. );
913  dbasis[NUMNODESP2xP4*1+22] = (4.* SQRT3/3.*gauss->coord1- 2.*SQRT3/3.*gauss->coord3 )*(-8./3.)*(zeta -1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
914  dbasis[NUMNODESP2xP4*2+22] = 4.*gauss->coord1*gauss->coord3* (-8./3.)*((2.*zeta -1. )*(zeta-0.5)*(zeta +1.) + zeta*(zeta -1.)*( 2.*zeta + 0.5));
915  /*Nodal function 24*/
916  dbasis[NUMNODESP2xP4*0+23] = 2.*(gauss->coord1- gauss->coord2) * (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. );
917  dbasis[NUMNODESP2xP4*1+23] = -2.*SQRT3/3.*(gauss->coord2+ gauss->coord1) * (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
918  dbasis[NUMNODESP2xP4*2+23] = 4.*gauss->coord1*gauss->coord2* (-8./3.)*((2.*zeta -1. )* (zeta-0.5) *(zeta +1.) + zeta* (zeta -1.)*( 2.*zeta + 0.5));
919  /*Nodal function 25*/
920  dbasis[NUMNODESP2xP4*0+24] = 2. *gauss->coord3 *4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. );
921  dbasis[NUMNODESP2xP4*1+24] = (4.*SQRT3/3.*gauss->coord2 - 2.* SQRT3/3. *gauss->coord3) *4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
922  dbasis[NUMNODESP2xP4*2+24] = 4.* gauss->coord2 * gauss->coord3* 4.*( 4.* zeta*zeta*zeta - 5./2. *zeta );
923  /*Nodal function 26*/
924  dbasis[NUMNODESP2xP4*0+25] = -2. *gauss->coord3*4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. );
925  dbasis[NUMNODESP2xP4*1+25] = (4.*SQRT3/3.*gauss->coord1- 2.*SQRT3/3.*gauss->coord3 )*4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
926  dbasis[NUMNODESP2xP4*2+25] = 4. * gauss->coord1 * gauss->coord3 *4.*( 4. *zeta*zeta*zeta - 5./2.* zeta );
927  /*Nodal function 27*/
928  dbasis[NUMNODESP2xP4*0+26] = 2.*( gauss->coord1-gauss->coord2) * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. );
929  dbasis[NUMNODESP2xP4*1+26] = -2.* SQRT3/3.*(gauss->coord1+ gauss->coord2 )*4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
930  dbasis[NUMNODESP2xP4*2+26] = 4. *gauss->coord1 *gauss->coord2 *4.*( 4.* zeta*zeta*zeta - 5./2.*zeta );
931  /*Nodal function 28*/
932  dbasis[NUMNODESP2xP4*0+27] = 2.* gauss->coord3* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. );
933  dbasis[NUMNODESP2xP4*1+27] = (4.*SQRT3/3.*gauss->coord2- 2.*SQRT3/3.*gauss->coord3) * (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
934  dbasis[NUMNODESP2xP4*2+27] = 4.* gauss->coord2*gauss->coord3* (-8./3.)*((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) +zeta*(zeta -1.)*( 2.*zeta + 3./2.));
935  /*Nodal function 29*/
936  dbasis[NUMNODESP2xP4*0+28] = -2.*gauss->coord3 *(-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. );
937  dbasis[NUMNODESP2xP4*1+28] = (4.*SQRT3/3.*gauss->coord1- 2.*SQRT3/3.*gauss->coord3) * (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
938  dbasis[NUMNODESP2xP4*2+28] = 4.* gauss->coord1*gauss->coord3* (-8./3.)*((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) +zeta*(zeta -1.)*( 2.*zeta + 3./2.));
939  /*Nodal function 30*/
940  dbasis[NUMNODESP2xP4*0+29] = 2.*(gauss->coord1-gauss->coord2)* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. );
941  dbasis[NUMNODESP2xP4*1+29] = -2.*SQRT3/3.*(gauss->coord1+gauss->coord2) * (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
942  dbasis[NUMNODESP2xP4*2+29] = 4.* gauss->coord1*gauss->coord2 * (-8./3.)*((2.*zeta -1. )*(zeta+0.5)*(zeta +1) +zeta*(zeta -1.)*( 2.*zeta + 3./2.));
943  return;
944  case P1xP3Enum :
945  /*Nodal function 1*/
946  dbasis[NUMNODESP1xP3*0+0 ] = (9./32.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
947  dbasis[NUMNODESP1xP3*1+0 ] = ((3.*SQRT3)/32.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
948  dbasis[NUMNODESP1xP3*2+0 ] =- (9./16.)* gauss->coord1 *( 2. *zeta *( zeta -1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
949  /*Nodal function 2*/
950  dbasis[NUMNODESP1xP3*0+1 ] = - (9./32.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
951  dbasis[NUMNODESP1xP3*1+1 ] = ((3.*SQRT3)/32.) *(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
952  dbasis[NUMNODESP1xP3*2+1 ] =- (9./16.)*gauss->coord2 *( 2.* zeta* ( zeta -1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
953  /*Nodal function 3*/
954  dbasis[NUMNODESP1xP3*0+2 ] = 0.;
955  dbasis[NUMNODESP1xP3*1+2 ] = - ((3.*SQRT3)/16.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
956  dbasis[NUMNODESP1xP3*2+2 ] = - (9./16.)* gauss->coord3* ( 2. *zeta *( zeta -1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
957  /*Nodal function 4*/
958  dbasis[NUMNODESP1xP3*0+3 ] = - (9./32.)*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
959  dbasis[NUMNODESP1xP3*1+3 ] = -((3.*SQRT3)/32.) *(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
960  dbasis[NUMNODESP1xP3*2+3 ] = (9./16.)* gauss->coord1*( 2.* zeta *( zeta +1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
961  /*Nodal function 5*/
962  dbasis[NUMNODESP1xP3*0+4 ] = (9./32.)* (zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
963  dbasis[NUMNODESP1xP3*1+4 ] = - ((3.*SQRT3)/32.) *(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
964  dbasis[NUMNODESP1xP3*2+4 ] = (9./16.)* gauss->coord2* ( 2.* zeta *( zeta +1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
965  /*Nodal function 6*/
966  dbasis[NUMNODESP1xP3*0+5 ] = 0.;
967  dbasis[NUMNODESP1xP3*1+5 ] = ((3.*SQRT3)/16.) *(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
968  dbasis[NUMNODESP1xP3*2+5 ] = (9./16.)* gauss->coord3 *( 2.* zeta * ( zeta + 1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
969  /*Nodal function 7*/
970  dbasis[NUMNODESP1xP3*0+6 ] = - (27./32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
971  dbasis[NUMNODESP1xP3*1+6 ] = - (9.*SQRT3/32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
972  dbasis[NUMNODESP1xP3*2+6 ] = gauss->coord1*(27./16.)*( 2.* zeta *( zeta - (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
973  /*Nodal function 8*/
974  dbasis[NUMNODESP1xP3*0+7 ] = (27./32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
975  dbasis[NUMNODESP1xP3*1+7 ] = -((9.*SQRT3)/32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
976  dbasis[NUMNODESP1xP3*2+7 ] = gauss->coord2*(27./16.)*( 2.* zeta *( zeta - (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
977  /*Nodal function 9*/
978  dbasis[NUMNODESP1xP3*0+8 ] = 0.;
979  dbasis[NUMNODESP1xP3*1+8 ] = ((9.*SQRT3)/16.) *(zeta-1.)*(zeta-1./3.)*(zeta+1.);
980  dbasis[NUMNODESP1xP3*2+8 ] = gauss->coord3*(27./16.)*( 2. *zeta *( zeta - (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
981  /*Nodal function 10*/
982  dbasis[NUMNODESP1xP3*0+9 ] = (27./32.) *(zeta-1.)*(zeta+1./3.)*(zeta+1.);
983  dbasis[NUMNODESP1xP3*1+9 ] = ((9.*SQRT3)/32.) *(zeta-1.)*(zeta+1./3.)*(zeta+1.);
984  dbasis[NUMNODESP1xP3*2+9 ] = -gauss->coord1 *(27./16.)*( 2* zeta *( zeta + (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
985  /*Nodal function 11*/
986  dbasis[NUMNODESP1xP3*0+10] = - (27./32.) *(zeta-1)*(zeta+1./3.)*(zeta+1);
987  dbasis[NUMNODESP1xP3*1+10] = ((9.*SQRT3)/32.) *(zeta-1.)*(zeta+1./3.)*(zeta+1);
988  dbasis[NUMNODESP1xP3*2+10] = -gauss->coord2 *(27./16.) *( 2.* zeta *( zeta + (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
989  /*Nodal function 12*/
990  dbasis[NUMNODESP1xP3*0+11] = 0.;
991  dbasis[NUMNODESP1xP3*1+11] = -((9.*SQRT3)/16.) *(zeta-1.)*(zeta+1./3.)*(zeta+1);
992  dbasis[NUMNODESP1xP3*2+11] = -gauss->coord3 *(27./16.)*( 2.* zeta *( zeta + (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
993  return;
994  case P1xP4Enum :
995  /*Nodal function 1*/
996  dbasis[NUMNODESP1xP4*0+0 ] = -0.5*(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
997  dbasis[NUMNODESP1xP4*1+0 ] = -SQRT3/6.*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
998  dbasis[NUMNODESP1xP4*2+0 ] = gauss->coord1 * 2./3.*( (2.*zeta-1)*(zeta -0.5)*(zeta +0.5) + 2.* zeta *zeta *(zeta -1.));
999  /*Nodal function 2*/
1000  dbasis[NUMNODESP1xP4*0+1 ] = +0.5*(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
1001  dbasis[NUMNODESP1xP4*1+1 ] = -SQRT3/6.*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
1002  dbasis[NUMNODESP1xP4*2+1 ] = gauss->coord2* 2./3.* ((2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2. * zeta *zeta*(zeta -1.));
1003  /*Nodal function 3*/
1004  dbasis[NUMNODESP1xP4*0+2 ] = 0. ;
1005  dbasis[NUMNODESP1xP4*1+2 ] = SQRT3/3.*(2./3.)*(zeta -1.)*(zeta-0.5)*(zeta)*(zeta+0.5);
1006  dbasis[NUMNODESP1xP4*2+2 ] = gauss->coord3* 2./3.*( (2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta -1.));
1007  /*Nodal function 4*/
1008  dbasis[NUMNODESP1xP4*0+3 ] = -0.5 *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
1009  dbasis[NUMNODESP1xP4*1+3 ] = -SQRT3/6.*(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
1010  dbasis[NUMNODESP1xP4*2+3 ] = gauss->coord1* 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.));
1011  /*Nodal function 5*/
1012  dbasis[NUMNODESP1xP4*0+4 ] = +0.5 * (2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
1013  dbasis[NUMNODESP1xP4*1+4 ] = -SQRT3/6.*(2./3.)*(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
1014  dbasis[NUMNODESP1xP4*2+4 ] = gauss->coord2 * 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.));
1015  /*Nodal function 6*/
1016  dbasis[NUMNODESP1xP4*0+5 ] = 0. ;
1017  dbasis[NUMNODESP1xP4*1+5 ] = SQRT3/3.*(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
1018  dbasis[NUMNODESP1xP4*2+5 ] = gauss->coord3 * 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1));
1019 
1020  /*Nodal function 7*/
1021  dbasis[NUMNODESP1xP4*0+6 ] = -0.5 * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
1022  dbasis[NUMNODESP1xP4*1+6 ] = -SQRT3/6.* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1.) ;
1023  dbasis[NUMNODESP1xP4*2+6 ] = gauss->coord1* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta );
1024  /*Nodal function 8*/
1025  dbasis[NUMNODESP1xP4*0+7 ] = +0.5* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
1026  dbasis[NUMNODESP1xP4*1+7 ] = -SQRT3/6.* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
1027  dbasis[NUMNODESP1xP4*2+7 ] = gauss->coord2* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta );
1028  /*Nodal function 9*/
1029  dbasis[NUMNODESP1xP4*0+8 ] = 0. ;
1030  dbasis[NUMNODESP1xP4*1+8 ] = SQRT3/3. * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ;
1031  dbasis[NUMNODESP1xP4*2+8 ] = gauss->coord3* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta );
1032 
1033  /*Nodal function 10*/
1034  dbasis[NUMNODESP1xP4*0+9 ] = -0.5* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
1035  dbasis[NUMNODESP1xP4*1+9 ] = -SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
1036  dbasis[NUMNODESP1xP4*2+9 ] = gauss->coord1* (-8./3.)*((2.*zeta -1.)*(zeta-0.5)*(zeta +1.) +zeta*(zeta -1.)*( 2.*zeta + 0.5));
1037  /*Nodal function 11*/
1038  dbasis[NUMNODESP1xP4*0+10] = +0.5* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
1039  dbasis[NUMNODESP1xP4*1+10] = -SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
1040  dbasis[NUMNODESP1xP4*2+10] = gauss->coord2* (-8./3.)*((2.*zeta -1.)*(zeta-0.5)*(zeta +1.) +zeta *(zeta -1.)*( 2.*zeta + 0.5));
1041  /*Nodal function 12*/
1042  dbasis[NUMNODESP1xP4*0+11] = 0. ;
1043  dbasis[NUMNODESP1xP4*1+11] = SQRT3/3.* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ;
1044  dbasis[NUMNODESP1xP4*2+11] = gauss->coord3* (-8./3.)*((2.*zeta-1.)*(zeta-0.5)*(zeta +1.) +zeta *(zeta -1.)*( 2.*zeta + 0.5));
1045  /*Nodal function 13*/
1046  dbasis[NUMNODESP1xP4*0+12] = -0.5* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
1047  dbasis[NUMNODESP1xP4*1+12] = -SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
1048  dbasis[NUMNODESP1xP4*2+12] = gauss->coord1* (-8./3.)*((2.*zeta -1. ) *(zeta+0.5)* (zeta +1.) + zeta* (zeta -1.)*( 2.*zeta + 3./2.));
1049  /*Nodal function 14*/
1050  dbasis[NUMNODESP1xP4*0+13] = +0.5* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
1051  dbasis[NUMNODESP1xP4*1+13] = -SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ;
1052  dbasis[NUMNODESP1xP4*2+13] = gauss->coord2* (-8./3.)*((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) + zeta*(zeta -1.)*( 2.*zeta + 3./2.));
1053  /*Nodal function 15*/
1054  dbasis[NUMNODESP1xP4*0+14] = 0 ;
1055  dbasis[NUMNODESP1xP4*1+14] = SQRT3/3.* (-8./3.)*(zeta - 1.)*(zeta + 0.5)*(zeta)*(zeta +1. ) ;
1056  dbasis[NUMNODESP1xP4*2+14] = gauss->coord3* (-8./3.)*((2. *zeta -1. )*(zeta+0.5)*(zeta +1.) + zeta*(zeta -1.)*( 2.*zeta + 3./2.));
1057 
1058  return;
1059  default:
1060  _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
1061  }
1062 
1063 }
1064 /*}}}*/
1065 void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1066  /*This routine returns the values of the nodal functions at the gaussian point.*/
1067 
1068  IssmDouble x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
1069 
1070  x1=xyz_list[0*3+0];
1071  y1=xyz_list[0*3+1];
1072  z1=xyz_list[0*3+2];
1073  x2=xyz_list[1*3+0];
1074  y2=xyz_list[1*3+1];
1075  z2=xyz_list[1*3+2];
1076  x3=xyz_list[2*3+0];
1077  y3=xyz_list[2*3+1];
1078  z3=xyz_list[2*3+2];
1079  x4=xyz_list[3*3+0];
1080  y4=xyz_list[3*3+1];
1081  z4=xyz_list[3*3+2];
1082 
1083  /*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/
1084  /*Area of a trabezoid = altitude * (base1 + base2)/2 */
1085  *Jdet= sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) * (z4-z1 + z3-z2)/8.;
1086  if(*Jdet<0.) _error_("negative jacobian determinant!");
1087 
1088 }
1089 /*}}}*/
1090 void PentaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1091  /*The Jacobian determinant is constant over the element, discard the gaussian points.
1092  * J is assumed to have been allocated of size 2x2.*/
1093 
1094  IssmDouble x1=xyz_list[3*0+0];
1095  IssmDouble y1=xyz_list[3*0+1];
1096  IssmDouble z1=xyz_list[3*0+2];
1097  IssmDouble x2=xyz_list[3*1+0];
1098  IssmDouble y2=xyz_list[3*1+1];
1099  IssmDouble z2=xyz_list[3*1+2];
1100 
1101  *Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
1102  if(*Jdet<0) _error_("negative jacobian determinant!");
1103 
1104 }
1105 /*}}}*/
1106 void PentaRef::GetTriaJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1107  /*The Jacobian determinant is constant over the element, discard the gaussian points.
1108  * J is assumed to have been allocated of size 2x2.*/
1109 
1110  IssmDouble x1=xyz_list[3*0+0];
1111  IssmDouble y1=xyz_list[3*0+1];
1112  IssmDouble z1=xyz_list[3*0+2];
1113  IssmDouble x2=xyz_list[3*1+0];
1114  IssmDouble y2=xyz_list[3*1+1];
1115  IssmDouble z2=xyz_list[3*1+2];
1116  IssmDouble x3=xyz_list[3*2+0];
1117  IssmDouble y3=xyz_list[3*2+1];
1118  IssmDouble z3=xyz_list[3*2+2];
1119 
1120  /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
1121  *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
1122  if(*Jdet<0) _error_("negative jacobian determinant!");
1123 }
1124 /*}}}*/
1125 void PentaRef::VerticalSegmentIndicesBase(int** pindices,int* pnumseg,int finiteelement){/*{{{*/
1126 
1127  /*Output*/
1128  int numseg;
1129  int* indices = NULL;
1130 
1131  switch(finiteelement){
1132  case P1Enum: case P1DGEnum:
1133  numseg = 3;
1134  indices = xNew<int>(numseg*2);
1135  indices[0*2 + 0] = 0; indices[0*2 + 1] = 3;
1136  indices[1*2 + 0] = 1; indices[1*2 + 1] = 4;
1137  indices[2*2 + 0] = 2; indices[2*2 + 1] = 5;
1138  break;
1139  case P1bubbleEnum: case P1bubblecondensedEnum:
1140  numseg = 3;
1141  indices = xNew<int>(numseg*2);
1142  indices[0*2 + 0] = 0; indices[0*2 + 1] = 3;
1143  indices[1*2 + 0] = 1; indices[1*2 + 1] = 4;
1144  indices[2*2 + 0] = 2; indices[2*2 + 1] = 5;
1145  break;
1146  case P2xP1Enum:
1147  numseg = 6;
1148  indices = xNew<int>(numseg*2);
1149  indices[0*2 + 0] = 0; indices[0*2 + 1] = 3;
1150  indices[1*2 + 0] = 1; indices[1*2 + 1] = 4;
1151  indices[2*2 + 0] = 2; indices[2*2 + 1] = 5;
1152  indices[3*2 + 0] = 6; indices[3*2 + 1] = 9;
1153  indices[4*2 + 0] = 7; indices[4*2 + 1] = 10;
1154  indices[5*2 + 0] = 8; indices[5*2 + 1] = 11;
1155  break;
1156  case P1xP2Enum:
1157  numseg = 3;
1158  indices = xNew<int>(numseg*2);
1159  indices[0*2 + 0] = 0; indices[0*2 + 1] = 6;
1160  indices[1*2 + 0] = 1; indices[1*2 + 1] = 7;
1161  indices[2*2 + 0] = 2; indices[2*2 + 1] = 8;
1162  break;
1163  case P1xP3Enum:
1164  numseg = 3;
1165  indices = xNew<int>(numseg*2);
1166  indices[0*2 + 0] = 0; indices[0*2 + 1] = 6;
1167  indices[1*2 + 0] = 1; indices[1*2 + 1] = 7;
1168  indices[2*2 + 0] = 2; indices[2*2 + 1] = 8;
1169  break;
1170  case P1xP4Enum:
1171  numseg = 3;
1172  indices = xNew<int>(numseg*2);
1173  indices[0*2 + 0] = 0; indices[0*2 + 1] = 6;
1174  indices[1*2 + 0] = 1; indices[1*2 + 1] = 7;
1175  indices[2*2 + 0] = 2; indices[2*2 + 1] = 8;
1176  break;
1177  case P2Enum:
1178  numseg = 6;
1179  indices = xNew<int>(numseg*2);
1180  indices[0*2 + 0] = 0; indices[0*2 + 1] = 6;
1181  indices[1*2 + 0] = 1; indices[1*2 + 1] = 7;
1182  indices[2*2 + 0] = 2; indices[2*2 + 1] = 8;
1183  indices[3*2 + 0] = 9; indices[3*2 + 1] = 15;
1184  indices[4*2 + 0] = 10; indices[4*2 + 1] = 16;
1185  indices[5*2 + 0] = 11; indices[5*2 + 1] = 17;
1186  break;
1187  case P2bubbleEnum:
1188  numseg = 6;
1189  indices = xNew<int>(numseg*2);
1190  indices[0*2 + 0] = 0; indices[0*2 + 1] = 6;
1191  indices[1*2 + 0] = 1; indices[1*2 + 1] = 7;
1192  indices[2*2 + 0] = 2; indices[2*2 + 1] = 8;
1193  indices[3*2 + 0] = 9; indices[3*2 + 1] = 15;
1194  indices[4*2 + 0] = 10; indices[4*2 + 1] = 16;
1195  indices[5*2 + 0] = 11; indices[5*2 + 1] = 17;
1196  break;
1197  case P2xP4Enum:
1198  numseg = 6;
1199  indices = xNew<int>(numseg*2);
1200  indices[0*2 + 0] = 0; indices[0*2 + 1] = 6;
1201  indices[1*2 + 0] = 1; indices[1*2 + 1] = 7;
1202  indices[2*2 + 0] = 2; indices[2*2 + 1] = 8;
1203  indices[3*2 + 0] = 9; indices[3*2 + 1] = 15;
1204  indices[4*2 + 0] = 10; indices[4*2 + 1] = 16;
1205  indices[5*2 + 0] = 11; indices[5*2 + 1] = 17;
1206  break;
1207  default:
1208  _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
1209  }
1210 
1211  /*Assign output pointer*/
1212  *pnumseg = numseg;
1213  *pindices = indices;
1214 }
1215 /*}}}*/
1216 int PentaRef::NumberofNodes(int finiteelement){/*{{{*/
1217 
1218  switch(finiteelement){
1219  case NoneEnum: return 0;
1220  case P0Enum: return NUMNODESP0;
1221  case P1Enum: return NUMNODESP1;
1222  case P1DGEnum: return NUMNODESP1;
1223  case P1bubbleEnum: return NUMNODESP1b;
1224  case P1bubblecondensedEnum: return NUMNODESP1b;
1225  case P2Enum: return NUMNODESP2;
1226  case P2bubbleEnum: return NUMNODESP2b;
1227  case P2bubblecondensedEnum: return NUMNODESP2b;
1228  case P2xP1Enum: return NUMNODESP2xP1;
1229  case P1xP2Enum: return NUMNODESP1xP2;
1230  case P2xP4Enum: return NUMNODESP2xP4;
1231  case P1xP3Enum: return NUMNODESP1xP3;
1232  case P1xP4Enum: return NUMNODESP1xP4;
1233  case P1P1Enum: return NUMNODESP1*2;
1234  case P1P1GLSEnum: return NUMNODESP1*2;
1235  case MINIcondensedEnum: return NUMNODESP1b+NUMNODESP1;
1236  case MINIEnum: return NUMNODESP1b+NUMNODESP1;
1237  case TaylorHoodEnum: return NUMNODESP2+NUMNODESP1;
1238  case LATaylorHoodEnum: return NUMNODESP2;
1239  case OneLayerP4zEnum: return NUMNODESP2xP4+NUMNODESP1;
1240  case CrouzeixRaviartEnum: return NUMNODESP2b+NUMNODESP1;
1241  case LACrouzeixRaviartEnum: return NUMNODESP2b;
1242  default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
1243  }
1244 
1245  return -1;
1246 }
1247 /*}}}*/
1248 int PentaRef::PressureInterpolation(int fe_stokes){/*{{{*/
1249 
1250  switch(fe_stokes){
1251  case P1P1Enum: return P1Enum;
1252  case P1P1GLSEnum: return P1Enum;
1253  case MINIcondensedEnum: return P1Enum;
1254  case MINIEnum: return P1Enum;
1255  case TaylorHoodEnum: return P1Enum;
1256  case LATaylorHoodEnum: return NoneEnum;
1257  case OneLayerP4zEnum: return P1Enum;
1258  case CrouzeixRaviartEnum: return P1DGEnum;
1259  case LACrouzeixRaviartEnum: return NoneEnum;
1260  default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
1261  }
1262 
1263  return -1;
1264 }
1265 /*}}}*/
1266 void PentaRef::SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
1267 
1268  /*Output*/
1269  int numindices;
1270  int* indices = NULL;
1271 
1272  switch(finiteelement){
1273  case P1Enum: case P1DGEnum:
1274  numindices = 3;
1275  indices = xNew<int>(numindices);
1276  indices[0] = 3;
1277  indices[1] = 4;
1278  indices[2] = 5;
1279  break;
1280  case P1bubbleEnum: case P1bubblecondensedEnum:
1281  numindices = 3;
1282  indices = xNew<int>(numindices);
1283  indices[0] = 3;
1284  indices[1] = 4;
1285  indices[2] = 5;
1286  break;
1287  case P2xP1Enum:
1288  numindices = 6;
1289  indices = xNew<int>(numindices);
1290  indices[0] = 3;
1291  indices[1] = 4;
1292  indices[2] = 5;
1293  indices[3] = 9;
1294  indices[4] = 10;
1295  indices[5] = 11;
1296  break;
1297  case P1xP2Enum:
1298  case P1xP3Enum:
1299  case P1xP4Enum:
1300  numindices = 3;
1301  indices = xNew<int>(numindices);
1302  indices[0] = 3;
1303  indices[1] = 4;
1304  indices[2] = 5;
1305  break;
1306  case P2Enum:
1307  numindices = 6;
1308  indices = xNew<int>(numindices);
1309  indices[0] = 3;
1310  indices[1] = 4;
1311  indices[2] = 5;
1312  indices[3] = 12;
1313  indices[4] = 13;
1314  indices[5] = 14;
1315  break;
1316  default:
1317  _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
1318  }
1319 
1320  /*Assign output pointer*/
1321  *pnumindices = numindices;
1322  *pindices = indices;
1323 }
1324 /*}}}*/
1325 int PentaRef::TensorInterpolation(int fe_stokes){/*{{{*/
1326 
1327  switch(fe_stokes){
1328  case XTaylorHoodEnum: return P1DGEnum;
1329  default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
1330  }
1331 
1332  return -1;
1333 }
1334 /*}}}*/
1335 int PentaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
1336 
1337  switch(fe_stokes){
1338  case P1P1Enum: return P1Enum;
1339  case P1P1GLSEnum: return P1Enum;
1340  case MINIcondensedEnum: return P1bubbleEnum;
1341  case MINIEnum: return P1bubbleEnum;
1342  case TaylorHoodEnum: return P2Enum;
1343  case LATaylorHoodEnum: return P2Enum;
1344  case OneLayerP4zEnum: return P2xP4Enum;
1345  case CrouzeixRaviartEnum: return P2bubbleEnum;
1346  case LACrouzeixRaviartEnum: return P2bubbleEnum;
1347  default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
1348  }
1349 
1350  return -1;
1351 }
1352 /*}}}*/
NUMNODESMAX
#define NUMNODESMAX
Definition: PentaRef.cpp:29
PentaRef::GetJacobianInvert
void GetJacobianInvert(IssmDouble *Jinv, IssmDouble *xyz_list, Gauss *gauss)
Definition: PentaRef.cpp:264
CrouzeixRaviartEnum
@ CrouzeixRaviartEnum
Definition: EnumDefinitions.h:1023
PentaRef::TensorInterpolation
int TensorInterpolation(int fe_stokes)
Definition: PentaRef.cpp:1325
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
PentaRef::GetNodalFunctionsDerivativesReference
void GetNodalFunctionsDerivativesReference(IssmDouble *dbasis, Gauss *gauss, int finiteelement)
Definition: PentaRef.cpp:496
GaussPenta::coord1
IssmDouble coord1
Definition: GaussPenta.h:24
MINIEnum
@ MINIEnum
Definition: EnumDefinitions.h:1156
P0Enum
@ P0Enum
Definition: EnumDefinitions.h:661
NUMNODESP1xP4
#define NUMNODESP1xP4
Definition: PentaRef.cpp:24
GaussPenta::coord3
IssmDouble coord3
Definition: GaussPenta.h:26
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
NUMNODESP1b
#define NUMNODESP1b
Definition: PentaRef.cpp:21
P2xP1Enum
@ P2xP1Enum
Definition: EnumDefinitions.h:1226
PentaRef::VerticalSegmentIndicesBase
void VerticalSegmentIndicesBase(int **pindices, int *pnumseg, int finiteelement)
Definition: PentaRef.cpp:1125
P1bubblecondensedEnum
@ P1bubblecondensedEnum
Definition: EnumDefinitions.h:1219
PentaRef::VelocityInterpolation
int VelocityInterpolation(int fe_stokes)
Definition: PentaRef.cpp:1335
NUMNODESP1
#define NUMNODESP1
Definition: PentaRef.cpp:19
NUMNODESP2xP4
#define NUMNODESP2xP4
Definition: PentaRef.cpp:28
NUMNODESP2b
#define NUMNODESP2b
Definition: PentaRef.cpp:27
PentaRef::GetNodalFunctionsDerivatives
void GetNodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss, int finiteelement)
Definition: PentaRef.cpp:466
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
TaylorHoodEnum
@ TaylorHoodEnum
Definition: EnumDefinitions.h:1299
Matrix3x3Invert
void Matrix3x3Invert(IssmDouble *Ainv, IssmDouble *A)
Definition: MatrixUtils.cpp:448
LATaylorHoodEnum
@ LATaylorHoodEnum
Definition: EnumDefinitions.h:1139
PentaRef::NumberofNodes
int NumberofNodes(int finiteelement)
Definition: PentaRef.cpp:1216
PentaRef::SurfaceNodeIndices
void SurfaceNodeIndices(int *pnumindices, int **pindices, int finiteelement)
Definition: PentaRef.cpp:1266
NUMNODESP1xP3
#define NUMNODESP1xP3
Definition: PentaRef.cpp:23
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
OneLayerP4zEnum
@ OneLayerP4zEnum
Definition: EnumDefinitions.h:1208
GaussPenta
Definition: GaussPenta.h:13
NUMNODESP2
#define NUMNODESP2
Definition: PentaRef.cpp:26
GaussPenta::coord2
IssmDouble coord2
Definition: GaussPenta.h:25
P2bubbleEnum
@ P2bubbleEnum
Definition: EnumDefinitions.h:1224
P2xP4Enum
@ P2xP4Enum
Definition: EnumDefinitions.h:1227
P1P1Enum
@ P1P1Enum
Definition: EnumDefinitions.h:1216
MINIcondensedEnum
@ MINIcondensedEnum
Definition: EnumDefinitions.h:1157
Gauss::Enum
virtual int Enum(void)=0
PentaRef::GetInputDerivativeValue
void GetInputDerivativeValue(IssmDouble *pvalues, IssmDouble *plist, IssmDouble *xyz_list, Gauss *gauss, int finiteelement)
Definition: PentaRef.cpp:131
P1bubbleEnum
@ P1bubbleEnum
Definition: EnumDefinitions.h:1218
PentaRef::PressureInterpolation
int PressureInterpolation(int fe_stokes)
Definition: PentaRef.cpp:1248
GaussPenta::coord4
IssmDouble coord4
Definition: GaussPenta.h:27
PentaRef::GetQuadJacobianDeterminant
void GetQuadJacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)
Definition: PentaRef.cpp:1065
PentaRef::BasalNodeIndices
void BasalNodeIndices(int *pnumindices, int **pindices, int finiteelement)
Definition: PentaRef.cpp:40
PentaRef::GetTriaJacobianDeterminant
void GetTriaJacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)
Definition: PentaRef.cpp:1106
LACrouzeixRaviartEnum
@ LACrouzeixRaviartEnum
Definition: EnumDefinitions.h:1138
NoneEnum
@ NoneEnum
Definition: EnumDefinitions.h:1202
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Matrix3x3Determinant
void Matrix3x3Determinant(IssmDouble *Adet, IssmDouble *A)
Definition: MatrixUtils.cpp:431
PentaRef::GetJacobian
void GetJacobian(IssmDouble *J, IssmDouble *xyz_list, Gauss *gauss)
Definition: PentaRef.cpp:191
PentaRef::PentaRef
PentaRef()
Definition: PentaRef.cpp:32
P2bubblecondensedEnum
@ P2bubblecondensedEnum
Definition: EnumDefinitions.h:1225
PentaRef::GetJacobianDeterminant
void GetJacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)
Definition: PentaRef.cpp:250
PentaRef::~PentaRef
~PentaRef()
Definition: PentaRef.cpp:35
P2Enum
@ P2Enum
Definition: EnumDefinitions.h:1223
P1P1GLSEnum
@ P1P1GLSEnum
Definition: EnumDefinitions.h:1217
P1xP4Enum
@ P1xP4Enum
Definition: EnumDefinitions.h:1222
PentaRef::GetInputValue
void GetInputValue(IssmDouble *pvalue, IssmDouble *plist, Gauss *gauss, int finiteelement)
Definition: PentaRef.cpp:169
P1xP3Enum
@ P1xP3Enum
Definition: EnumDefinitions.h:1221
SQRT3
#define SQRT3
Definition: constants.h:10
NUMNODESP2xP1
#define NUMNODESP2xP1
Definition: PentaRef.cpp:25
NUMNODESP1xP2
#define NUMNODESP1xP2
Definition: PentaRef.cpp:22
NUMNODESP0
#define NUMNODESP0
Definition: PentaRef.cpp:18
Gauss
Definition: Gauss.h:8
XTaylorHoodEnum
@ XTaylorHoodEnum
Definition: EnumDefinitions.h:1329
P1xP2Enum
@ P1xP2Enum
Definition: EnumDefinitions.h:1220
PentaRef::GetNodalFunctions
void GetNodalFunctions(IssmDouble *basis, Gauss *gauss, int finiteelement)
Definition: PentaRef.cpp:276
PentaRef::GetSegmentJacobianDeterminant
void GetSegmentJacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)
Definition: PentaRef.cpp:1090
GaussPentaEnum
@ GaussPentaEnum
Definition: EnumDefinitions.h:1078
Generated on Thu Jul 2 2020 08:09:17 for Ice Sheet System Model by   doxygen 1.8.19