source: issm/oecreview/Archive/15392-16133/ISSM-15469-15470.diff@ 16134

Last change on this file since 16134 was 16134, checked in by Mathieu Morlighem, 12 years ago

Added Archive/15392-16133
File size: 66.5 KB
RevLine 
[16134]1Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15469)
4+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15470)
5@@ -70,15 +70,15 @@
6 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
7
8 /*Build B: */
9- for (int i=0;i<NUMNODESP1;i++){
10- *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
11- *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
12+ for(int i=0;i<NUMNODESP1;i++){
13+ B[NDOF2*NUMNODESP1*0+NDOF2*i+0] = dbasis[0][i];
14+ B[NDOF2*NUMNODESP1*0+NDOF2*i+1] = 0.;
15
16- *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
17- *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
18+ B[NDOF2*NUMNODESP1*1+NDOF2*i+0] = 0.;
19+ B[NDOF2*NUMNODESP1*1+NDOF2*i+1] = dbasis[1][i];
20
21- *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i];
22- *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i];
23+ B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
24+ B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
25 }
26 }
27 /*}}}*/
28@@ -96,35 +96,35 @@
29 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
30 */
31
32- int i;
33- IssmDouble dh1dh7[3][NUMNODESMINI];
34- IssmDouble l1l6[NUMNODESP1];
35+ int i;
36+ IssmDouble dbasismini[3][NUMNODESMINI];
37+ IssmDouble basis[NUMNODESP1];
38
39- /*Get dh1dh6 in actual coordinate system: */
40- GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
41- GetNodalFunctionsP1(l1l6, gauss);
42+ /*Get dbasis in actual coordinate system: */
43+ GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
44+ GetNodalFunctionsP1(basis,gauss);
45
46 /*Build B: */
47- for (i=0;i<NUMNODESMINI;i++){
48- *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
49- *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
50- *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
51- *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
52- *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
53- *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
54- *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.5*dh1dh7[1][i];
55- *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.5*dh1dh7[0][i];
56- *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
57- *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;
58- *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;
59- *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
60+ for(i=0;i<NUMNODESMINI;i++){
61+ B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
62+ B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
63+ B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
64+ B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
65+ B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
66+ B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
67+ B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.5*dbasismini[1][i];
68+ B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.5*dbasismini[0][i];
69+ B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
70+ B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = 0.;
71+ B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = 0.;
72+ B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
73 }
74
75- for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
76- *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
77- *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
78- *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
79- *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i];
80+ for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
81+ B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0;
82+ B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0;
83+ B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0;
84+ B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = basis[i];
85 }
86 }
87 /*}}}*/
88@@ -150,20 +150,20 @@
89
90 /*Build B: */
91 for (int i=0;i<NUMNODESP1;i++){
92- *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
93- *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
94+ B[NDOF2*NUMNODESP1*0+NDOF2*i+0] = dbasis[0][i];
95+ B[NDOF2*NUMNODESP1*0+NDOF2*i+1] = 0.;
96
97- *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
98- *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
99+ B[NDOF2*NUMNODESP1*1+NDOF2*i+0] = 0.;
100+ B[NDOF2*NUMNODESP1*1+NDOF2*i+1] = dbasis[1][i];
101
102- *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i];
103- *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i];
104+ B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
105+ B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
106
107- *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dbasis[2][i];
108- *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
109+ B[NDOF2*NUMNODESP1*3+NDOF2*i+0] = .5*dbasis[2][i];
110+ B[NDOF2*NUMNODESP1*3+NDOF2*i+1] = 0.;
111
112- *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
113- *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dbasis[2][i];
114+ B[NDOF2*NUMNODESP1*4+NDOF2*i+0] = 0.;
115+ B[NDOF2*NUMNODESP1*4+NDOF2*i+1] = .5*dbasis[2][i];
116 }
117
118 }
119@@ -188,21 +188,21 @@
120 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss_coord);
121
122 /*Build BPrime: */
123- for (int i=0;i<NUMNODESP1;i++){
124- *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dbasis[0][i];
125- *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i];
126+ for(int i=0;i<NUMNODESP1;i++){
127+ B[NDOF2*NUMNODESP1*0+NDOF2*i+0]=2.*dbasis[0][i];
128+ B[NDOF2*NUMNODESP1*0+NDOF2*i+1]=dbasis[1][i];
129
130- *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
131- *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dbasis[1][i];
132+ B[NDOF2*NUMNODESP1*1+NDOF2*i+0]=dbasis[0][i];
133+ B[NDOF2*NUMNODESP1*1+NDOF2*i+1]=2.*dbasis[1][i];
134
135- *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i];
136- *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i];
137+ B[NDOF2*NUMNODESP1*2+NDOF2*i+0]=dbasis[1][i];
138+ B[NDOF2*NUMNODESP1*2+NDOF2*i+1]=dbasis[0][i];
139
140- *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[2][i];
141- *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
142+ B[NDOF2*NUMNODESP1*3+NDOF2*i+0]=dbasis[2][i];
143+ B[NDOF2*NUMNODESP1*3+NDOF2*i+1]=0.;
144
145- *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
146- *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dbasis[2][i];
147+ B[NDOF2*NUMNODESP1*4+NDOF2*i+0]=0.;
148+ B[NDOF2*NUMNODESP1*4+NDOF2*i+1]=dbasis[2][i];
149 }
150 }
151 /*}}}*/
152@@ -220,28 +220,28 @@
153 */
154
155 int i;
156- IssmDouble dh1dh7[3][NUMNODESMINI];
157+ IssmDouble dbasismini[3][NUMNODESMINI];
158
159- /*Get dh1dh6 in actual coordinate system: */
160- GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
161+ /*Get dbasis in actual coordinate system: */
162+ GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
163
164 /*Build Bprime: */
165- for (i=0;i<NUMNODESMINI;i++){
166- *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=2*dh1dh7[0][i]; //Bprime[0][NDOF4*i]=dh1dh6[0][i];
167- *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=dh1dh7[1][i];
168- *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
169- *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=dh1dh7[0][i];
170- *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=2*dh1dh7[1][i];
171- *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
172- *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=dh1dh7[1][i];
173- *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i];
174- *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
175+ for(i=0;i<NUMNODESMINI;i++){
176+ Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = 2.*dbasismini[0][i];
177+ Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = dbasismini[1][i];
178+ Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
179+ Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = dbasismini[0][i];
180+ Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = 2.*dbasismini[1][i];
181+ Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
182+ Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = dbasismini[1][i];
183+ Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = dbasismini[0][i];
184+ Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
185 }
186
187- for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
188- *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
189- *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
190- *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
191+ for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
192+ Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
193+ Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
194+ Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
195 }
196
197 }
198@@ -265,50 +265,50 @@
199
200 int i;
201
202- IssmDouble dh1dh7[3][NUMNODESMINI];
203- IssmDouble l1l6[NUMNODESP1];
204+ IssmDouble dbasismini[3][NUMNODESMINI];
205+ IssmDouble basis[NUMNODESP1];
206
207- /*Get dh1dh7 in actual coordinate system: */
208- GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
209- GetNodalFunctionsP1(l1l6, gauss);
210+ /*Get dbasismini in actual coordinate system: */
211+ GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
212+ GetNodalFunctionsP1(basis, gauss);
213
214 /*Build B: */
215- for (i=0;i<NUMNODESMINI;i++){
216- B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dh1dh7[0][i+0]; //B[0][NDOF4*i+0] = dh1dh6[0][i+0];
217+ for(i=0;i<NUMNODESMINI;i++){
218+ B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i+0]; //B[0][NDOF4*i+0] = dbasis[0][i+0];
219 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
220 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
221 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
222- B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dh1dh7[1][i+0];
223+ B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i+0];
224 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
225 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
226 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
227- B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dh1dh7[2][i+0];
228- B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = (float).5*dh1dh7[1][i+0];
229- B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = (float).5*dh1dh7[0][i+0];
230+ B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i+0];
231+ B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = .5*dbasismini[1][i+0];
232+ B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = .5*dbasismini[0][i+0];
233 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
234- B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = (float).5*dh1dh7[2][i+0];
235+ B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = .5*dbasismini[2][i+0];
236 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
237- B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = (float).5*dh1dh7[0][i+0];
238+ B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = .5*dbasismini[0][i+0];
239 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
240- B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = (float).5*dh1dh7[2][i+0];
241- B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = (float).5*dh1dh7[1][i+0];
242+ B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = .5*dbasismini[2][i+0];
243+ B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = .5*dbasismini[1][i+0];
244 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = 0.;
245 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = 0.;
246 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = 0.;
247- B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dh1dh7[0][i+0];
248- B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dh1dh7[1][i+0];
249- B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dh1dh7[2][i+0];
250+ B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dbasismini[0][i+0];
251+ B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dbasismini[1][i+0];
252+ B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dbasismini[2][i+0];
253 }
254
255- for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
256- *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
257- *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
258- *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
259- *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
260- *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
261- *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
262- *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i];
263- *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0;
264+ for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
265+ B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
266+ B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
267+ B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
268+ B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
269+ B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
270+ B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
271+ B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = basis[i];
272+ B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = 0.;
273 }
274
275 }
276@@ -331,53 +331,50 @@
277 */
278
279 int i;
280+ IssmDouble dbasis[3][NUMNODESP1];
281+ IssmDouble basis[NUMNODESP1];
282
283- IssmDouble dh1dh6[3][NUMNODESP1];
284- IssmDouble l1l6[NUMNODESP1];
285+ /*Get dbasismini in actual coordinate system: */
286+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
287+ GetNodalFunctionsP1(&basis[0], gauss);
288
289- /*Get dh1dh7 in actual coordinate system: */
290- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
291- GetNodalFunctionsP1(&l1l6[0], gauss);
292-
293 /*Build B: */
294- for (i=0;i<NUMNODESP1;i++){
295- *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
296- *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.;
297- *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.;
298- *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.;
299- *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i];
300- *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.;
301- *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.;
302- *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.;
303- *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i];
304- *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i)=.5*dh1dh6[1][i];
305- *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=.5*dh1dh6[0][i];
306- *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.;
307- *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i)=.5*dh1dh6[2][i];
308- *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.;
309- *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=.5*dh1dh6[0][i];
310- *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.;
311- *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=.5*dh1dh6[2][i];
312- *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=.5*dh1dh6[1][i];
313- *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i)=0.;
314- *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=0.;
315- *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=0.;
316- *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i)=dh1dh6[0][i];
317- *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=dh1dh6[1][i];
318- *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=dh1dh6[2][i];
319- _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
320+ for(i=0;i<NUMNODESP1;i++){
321+ B[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
322+ B[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
323+ B[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
324+ B[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
325+ B[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
326+ B[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
327+ B[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
328+ B[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
329+ B[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
330+ B[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = .5*dbasis[1][i];
331+ B[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = .5*dbasis[0][i];
332+ B[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
333+ B[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = .5*dbasis[2][i];
334+ B[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
335+ B[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = .5*dbasis[0][i];
336+ B[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
337+ B[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = .5*dbasis[2][i];
338+ B[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = .5*dbasis[1][i];
339+ B[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = 0.;
340+ B[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = 0.;
341+ B[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = 0.;
342+ B[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = dbasis[0][i];
343+ B[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = dbasis[1][i];
344+ B[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = dbasis[2][i];
345 }
346
347- for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
348- B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3]=0.;
349- B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3]=0.;
350- B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3]=0.;
351- B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3]=0.;
352- B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3]=0.;
353- B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3]=0.;
354- B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3]=l1l6[i];
355- B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3]=0.;
356- _assert_(((NDOF4*NUMNODESP1)*7+NDOF4*i+3)<8*24);
357+ for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
358+ B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
359+ B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
360+ B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
361+ B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
362+ B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
363+ B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
364+ B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = basis[i];
365+ B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = 0.;
366 }
367
368 }
369@@ -401,50 +398,50 @@
370 */
371
372 int i;
373- IssmDouble dh1dh7[3][NUMNODESMINI];
374- IssmDouble l1l6[NUMNODESP1];
375+ IssmDouble dbasismini[3][NUMNODESMINI];
376+ IssmDouble basis[NUMNODESP1];
377
378- /*Get dh1dh7 in actual coordinate system: */
379- GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
380- GetNodalFunctionsP1(l1l6, gauss);
381+ /*Get dbasismini in actual coordinate system: */
382+ GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
383+ GetNodalFunctionsP1(basis, gauss);
384
385 /*B_primeuild B_prime: */
386- for (i=0;i<NUMNODESMINI;i++){
387- *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
388- *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
389- *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
390- *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
391- *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
392- *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
393- *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;
394- *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;
395- *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];
396- *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh7[1][i];
397- *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh7[0][i];
398- *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
399- *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh7[2][i];
400- *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;
401- *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh7[0][i];
402- *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;
403- *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh7[2][i];
404- *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh7[1][i];
405- *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh7[0][i];
406- *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh7[1][i];
407- *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh7[2][i];
408- *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0;
409- *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0;
410- *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0;
411+ for(i=0;i<NUMNODESMINI;i++){
412+ B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
413+ B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
414+ B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
415+ B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
416+ B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
417+ B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
418+ B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
419+ B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
420+ B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i];
421+ B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = dbasismini[1][i];
422+ B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = dbasismini[0][i];
423+ B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
424+ B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = dbasismini[2][i];
425+ B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
426+ B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = dbasismini[0][i];
427+ B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
428+ B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = dbasismini[2][i];
429+ B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = dbasismini[1][i];
430+ B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = dbasismini[0][i];
431+ B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = dbasismini[1][i];
432+ B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = dbasismini[2][i];
433+ B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = 0.;
434+ B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = 0.;
435+ B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = 0.;
436 }
437
438- for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
439- *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
440- *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
441- *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
442- *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
443- *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
444- *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
445- *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0;
446- *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i];
447+ for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
448+ B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
449+ B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
450+ B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
451+ B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
452+ B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
453+ B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
454+ B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = 0.;
455+ B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = basis[i];
456 }
457
458 }
459@@ -468,52 +465,50 @@
460 */
461
462 int i;
463- IssmDouble dh1dh6[3][NUMNODESP1];
464- IssmDouble l1l6[NUMNODESP1];
465+ IssmDouble dbasis[3][NUMNODESP1];
466+ IssmDouble basis[NUMNODESP1];
467
468- /*Get dh1dh7 in actual coordinate system: */
469- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
470- GetNodalFunctionsP1(l1l6, gauss);
471+ /*Get dbasismini in actual coordinate system: */
472+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
473+ GetNodalFunctionsP1(basis, gauss);
474
475 /*B_primeuild B_prime: */
476- for (i=0;i<NUMNODESP1;i++){
477- *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
478- *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.;
479- *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.;
480- *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.;
481- *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i];
482- *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.;
483- *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.;
484- *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.;
485- *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i];
486- *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i)=dh1dh6[1][i];
487- *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=dh1dh6[0][i];
488- *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.;
489- *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i)=dh1dh6[2][i];
490- *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.;
491- *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=dh1dh6[0][i];
492- *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.;
493- *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=dh1dh6[2][i];
494- *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=dh1dh6[1][i];
495- *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i)=dh1dh6[0][i];
496- *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=dh1dh6[1][i];
497- *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=dh1dh6[2][i];
498- *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i)=0.;
499- *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=0.;
500- *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=0.;
501- _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
502+ for(i=0;i<NUMNODESP1;i++){
503+ B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
504+ B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
505+ B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
506+ B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
507+ B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
508+ B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
509+ B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
510+ B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
511+ B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
512+ B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = dbasis[1][i];
513+ B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = dbasis[0][i];
514+ B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
515+ B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = dbasis[2][i];
516+ B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
517+ B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = dbasis[0][i];
518+ B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
519+ B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = dbasis[2][i];
520+ B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = dbasis[1][i];
521+ B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = dbasis[0][i];
522+ B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = dbasis[1][i];
523+ B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = dbasis[2][i];
524+ B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = 0.;
525+ B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = 0.;
526+ B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = 0.;
527 }
528
529- for (i=0;i<NUMNODESP1;i++){ //last column
530- *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+3)=0.;
531- *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+3)=0.;
532- *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+3)=0.;
533- *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+3)=0.;
534- *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+3)=0.;
535- *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+3)=0.;
536- *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+3)=0.;
537- *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+3)=l1l6[i];
538- _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
539+ for(i=0;i<NUMNODESP1;i++){ //last column
540+ B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
541+ B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
542+ B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
543+ B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
544+ B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
545+ B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
546+ B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = 0.;
547+ B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = basis[i];
548 }
549
550 }
551@@ -532,16 +527,16 @@
552 */
553
554 /*Same thing in the actual coordinate system: */
555- IssmDouble l1l6[6];
556+ IssmDouble basis[6];
557
558 /*Get dh1dh2dh3 in actual coordinates system : */
559- GetNodalFunctionsP1(l1l6, gauss);
560+ GetNodalFunctionsP1(basis, gauss);
561
562 /*Build B': */
563- for (int i=0;i<NUMNODESP1;i++){
564- *(B_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=l1l6[i];
565- *(B_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=l1l6[i];
566- *(B_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=l1l6[i];
567+ for(int i=0;i<NUMNODESP1;i++){
568+ B_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = basis[i];
569+ B_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = basis[i];
570+ B_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = basis[i];
571 }
572 }
573 /*}}}*/
574@@ -559,16 +554,16 @@
575 */
576
577 /*Same thing in the actual coordinate system: */
578- IssmDouble dh1dh6[3][NUMNODESP1];
579+ IssmDouble dbasis[3][NUMNODESP1];
580
581 /*Get dh1dh2dh3 in actual coordinates system : */
582- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
583+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
584
585 /*Build B': */
586- for (int i=0;i<NUMNODESP1;i++){
587- *(B_conduct+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
588- *(B_conduct+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
589- *(B_conduct+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i];
590+ for(int i=0;i<NUMNODESP1;i++){
591+ B_conduct[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
592+ B_conduct[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
593+ B_conduct[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
594 }
595 }
596 /*}}}*/
597@@ -577,15 +572,13 @@
598 /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
599 where hi is the interpolation function for node i.*/
600
601- int i;
602- IssmDouble dh1dh6[3][NUMNODESP1];
603+ /*Get dbasis in actual coordinate system: */
604+ IssmDouble dbasis[3][NUMNODESP1];
605+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
606
607- /*Get dh1dh6 in actual coordinate system: */
608- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
609-
610 /*Build B: */
611- for (i=0;i<NUMNODESP1;i++){
612- B[i]=dh1dh6[2][i];
613+ for(int i=0;i<NUMNODESP1;i++){
614+ B[i] = dbasis[2][i];
615 }
616
617 }
618@@ -603,23 +596,20 @@
619 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
620 */
621
622- /*Same thing in the actual coordinate system: */
623- IssmDouble dh1dh6[3][NUMNODESP1];
624+ /*Get nodal function derivatives in actual coordinates system : */
625+ IssmDouble dbasis[3][NUMNODESP1];
626+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
627
628- /*Get dh1dh2dh3 in actual coordinates system : */
629- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
630-
631 /*Build B': */
632- for (int i=0;i<NUMNODESP1;i++){
633- *(Bprime_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
634- *(Bprime_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
635- *(Bprime_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i];
636+ for(int i=0;i<NUMNODESP1;i++){
637+ Bprime_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
638+ Bprime_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
639+ Bprime_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
640 }
641 }
642 /*}}}*/
643 /*FUNCTION PentaRef::GetBprimeVert{{{*/
644 void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
645- /* Compute Bprime matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for node i*/
646
647 GetNodalFunctionsP1(B, gauss);
648
649@@ -637,23 +627,21 @@
650 ** We assume B has been allocated already, of size: 2 (2 x numnodes)
651 **/
652
653+ /*Get basis in actual coordinate system: */
654 IssmDouble basis[6];
655-
656- /*Get l1l6 in actual coordinate system: */
657 GetNodalFunctionsP1(&basis[0],gauss);
658
659 for(int i=0;i<NUMNODESP1;i++){
660- B[2*NUMNODESP1*0+2*i+0]=basis[i];
661- B[2*NUMNODESP1*0+2*i+1]=0.;
662- B[2*NUMNODESP1*1+2*i+0]=0.;
663- B[2*NUMNODESP1*1+2*i+1]=basis[i];
664+ B[2*NUMNODESP1*0+2*i+0] = basis[i];
665+ B[2*NUMNODESP1*0+2*i+1] = 0.;
666+ B[2*NUMNODESP1*1+2*i+0] = 0.;
667+ B[2*NUMNODESP1*1+2*i+1] = basis[i];
668 }
669 }
670 /*}}}*/
671 /*FUNCTION PentaRef::GetLStokes{{{*/
672 void PentaRef::GetLStokes(IssmDouble* LStokes, GaussPenta* gauss){
673- /*
674- * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
675+ /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
676 * For node i, Li can be expressed in the actual coordinate system
677 * by:
678 * Li=[ h 0 ]
679@@ -672,24 +660,22 @@
680 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
681
682 /*Build LStokes: */
683- for (int i=0;i<3;i++){
684- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i];
685- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
686- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
687- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
688+ for(int i=0;i<3;i++){
689+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
690+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
691+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
692+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
693
694- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
695- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
696- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
697- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
698+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
699+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
700+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
701+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
702 }
703 }
704 /*}}}*/
705 /*FUNCTION PentaRef::GetLprimeStokes {{{*/
706 void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
707-
708- /*
709- * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
710+ /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
711 * For node i, Lpi can be expressed in the actual coordinate system
712 * by:
713 * Lpi=[ h 0 0 0]1
714@@ -727,73 +713,73 @@
715 int num_dof=4;
716
717 IssmDouble l1l2l3[NUMNODESP1_2d];
718- IssmDouble dh1dh6[3][NUMNODESP1];
719+ IssmDouble dbasis[3][NUMNODESP1];
720
721 /*Get l1l2l3 in actual coordinate system: */
722 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
723 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
724 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
725
726- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
727+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
728
729 /*Build LprimeStokes: */
730- for (i=0;i<3;i++){
731- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
732- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
733- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
734- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
735- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
736- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
737- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
738- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
739- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
740- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
741- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=0;
742- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
743- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
744- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
745- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0;
746- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
747- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
748- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
749- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=l1l2l3[i];
750- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
751- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
752- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
753- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=l1l2l3[i];
754- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
755- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
756- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
757- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=dh1dh6[2][i];
758- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
759- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
760- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
761- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=dh1dh6[2][i];
762- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
763- *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
764- *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
765- *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=dh1dh6[2][i];
766- *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
767- *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=dh1dh6[2][i];
768- *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
769- *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=dh1dh6[0][i];
770- *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
771- *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
772- *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=dh1dh6[2][i];
773- *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=dh1dh6[1][i];
774- *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
775- *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=0;
776- *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
777- *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
778- *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=l1l2l3[i];
779- *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
780- *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=0;
781- *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
782- *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=l1l2l3[i];
783- *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
784- *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
785- *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=0;
786- *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=l1l2l3[i];
787+ for(int i=0;i<3;i++){
788+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
789+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
790+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
791+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
792+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
793+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
794+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
795+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
796+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
797+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
798+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = 0.;
799+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
800+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
801+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
802+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = 0.;
803+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
804+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
805+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
806+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = l1l2l3[i];
807+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
808+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
809+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
810+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = l1l2l3[i];
811+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
812+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
813+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
814+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = dbasis[2][i];
815+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = 0.;
816+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
817+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
818+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = dbasis[2][i];
819+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = 0.;
820+ LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+0] = 0.;
821+ LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+1] = 0.;
822+ LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+2] = dbasis[2][i];
823+ LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+3] = 0.;
824+ LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+0] = dbasis[2][i];
825+ LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+1] = 0.;
826+ LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+2] = dbasis[0][i];
827+ LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+3] = 0.;
828+ LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
829+ LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
830+ LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
831+ LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
832+ LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
833+ LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
834+ LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
835+ LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = l1l2l3[i];
836+ LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
837+ LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
838+ LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
839+ LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = l1l2l3[i];
840+ LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
841+ LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
842+ LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
843+ LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = l1l2l3[i];
844 }
845 }
846 /*}}}*/
847@@ -814,9 +800,7 @@
848 * where h is the interpolation function for node i.
849 */
850
851- int i;
852 int num_dof=2;
853-
854 IssmDouble l1l2l3[NUMNODESP1_2d];
855
856 /*Get l1l2l3 in actual coordinate system: */
857@@ -825,32 +809,29 @@
858 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
859
860 /*Build LStokes: */
861- for (i=0;i<3;i++){
862- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
863- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
864- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
865- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
866- *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
867- *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
868- *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
869- *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
870- *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i];
871- *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
872- *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
873- *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i];
874- *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i];
875- *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
876- *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
877- *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i];
878-
879+ for(int i=0;i<3;i++){
880+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
881+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
882+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
883+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
884+ LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
885+ LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
886+ LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
887+ LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
888+ LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = l1l2l3[i];
889+ LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
890+ LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
891+ LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = l1l2l3[i];
892+ LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = l1l2l3[i];
893+ LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
894+ LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
895+ LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = l1l2l3[i];
896 }
897 }
898 /*}}}*/
899 /*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/
900 void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
901-
902- /*
903- * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
904+ /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
905 * For node i, Lpi can be expressed in the actual coordinate system
906 * by:
907 * Lpi=[ h 0 0 0]
908@@ -863,60 +844,57 @@
909 * [ 0 0 0 h]
910 * where h is the interpolation function for node i.
911 */
912- int i;
913 int num_dof=4;
914-
915 IssmDouble l1l2l3[NUMNODESP1_2d];
916- IssmDouble dh1dh6[3][NUMNODESP1];
917+ IssmDouble dbasis[3][NUMNODESP1];
918
919 /*Get l1l2l3 in actual coordinate system: */
920 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
921 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
922 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
923
924- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
925+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
926
927 /*Build LprimeStokes: */
928- for (i=0;i<3;i++){
929- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
930- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
931- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
932- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
933- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
934- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
935- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
936- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
937- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
938- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
939- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
940- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
941- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
942- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
943- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
944- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
945- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
946- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
947- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=dh1dh6[2][i];
948- *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
949- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
950- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
951- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=dh1dh6[2][i];
952- *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
953- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
954- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
955- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0;
956- *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=l1l2l3[i];
957- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
958- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
959- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0;
960- *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=l1l2l3[i];
961+ for(int i=0;i<3;i++){
962+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
963+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
964+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
965+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
966+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
967+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
968+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
969+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
970+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
971+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
972+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
973+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
974+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
975+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
976+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
977+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
978+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
979+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
980+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
981+ LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
982+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
983+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
984+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
985+ LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
986+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
987+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
988+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
989+ LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = l1l2l3[i];
990+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
991+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
992+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
993+ LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = l1l2l3[i];
994 }
995 }
996 /*}}}*/
997 /*FUNCTION PentaRef::GetLStokesMacAyeal {{{*/
998 void PentaRef::GetLStokesMacAyeal(IssmDouble* LStokes, GaussPenta* gauss){
999- /*
1000- * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
1001+ /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
1002 * For node i, Li can be expressed in the actual coordinate system
1003 * by:
1004 * Li=[ h 0 0 0]
1005@@ -926,9 +904,7 @@
1006 * where h is the interpolation function for node i.
1007 */
1008
1009- int i;
1010 int num_dof=4;
1011-
1012 IssmDouble l1l2l3[NUMNODESP1_2d];
1013
1014 /*Get l1l2l3 in actual coordinate system: */
1015@@ -937,32 +913,29 @@
1016 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
1017
1018 /*Build LStokes: */
1019- for (i=0;i<3;i++){
1020- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
1021- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
1022- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
1023- *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
1024- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
1025- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
1026- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
1027- *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
1028- *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
1029- *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
1030- *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
1031- *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
1032- *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
1033- *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
1034- *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
1035- *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
1036-
1037+ for(int i=0;i<3;i++){
1038+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
1039+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
1040+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
1041+ LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
1042+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
1043+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
1044+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
1045+ LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
1046+ LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
1047+ LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
1048+ LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
1049+ LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
1050+ LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
1051+ LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
1052+ LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
1053+ LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
1054 }
1055 }
1056 /*}}}*/
1057 /*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/
1058 void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
1059-
1060- /*
1061- * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
1062+ /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
1063 * For node i, Lpi can be expressed in the actual coordinate system
1064 * by:
1065 * Lpi=[ h 0 ]
1066@@ -971,87 +944,80 @@
1067 * [ 0 h ]
1068 * where h is the interpolation function for node i.
1069 */
1070- int i;
1071 int num_dof=2;
1072-
1073 IssmDouble l1l2l3[NUMNODESP1_2d];
1074- IssmDouble dh1dh6[3][NUMNODESP1];
1075+ IssmDouble dbasis[3][NUMNODESP1];
1076
1077 /*Get l1l2l3 in actual coordinate system: */
1078 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
1079 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
1080 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
1081+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
1082
1083- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
1084-
1085 /*Build LprimeStokes: */
1086- for (i=0;i<3;i++){
1087- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
1088- *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
1089- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
1090- *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
1091- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
1092- *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
1093- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
1094- *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
1095+ for(int i=0;i<3;i++){
1096+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
1097+ LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
1098+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
1099+ LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
1100+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
1101+ LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
1102+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
1103+ LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
1104 }
1105 }
1106 /*}}}*/
1107 /*FUNCTION PentaRef::GetJacobian {{{*/
1108 void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussPenta* gauss){
1109-
1110 /*The Jacobian is constant over the element, discard the gaussian points.
1111 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
1112
1113- IssmDouble A1,A2,A3; //area coordinates
1114- IssmDouble xi,eta,zi; //parametric coordinates
1115-
1116+ IssmDouble A1,A2,A3; // area coordinates
1117+ IssmDouble xi,eta,zi; // parametric coordinates
1118 IssmDouble x1,x2,x3,x4,x5,x6;
1119 IssmDouble y1,y2,y3,y4,y5,y6;
1120 IssmDouble z1,z2,z3,z4,z5,z6;
1121
1122 /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
1123- A1=gauss->coord1;
1124- A2=gauss->coord2;
1125- A3=gauss->coord3;
1126-
1127+ A1 = gauss->coord1;
1128+ A2 = gauss->coord2;
1129+ A3 = gauss->coord3;
1130 xi = A2-A1;
1131 eta = SQRT3*A3;
1132 zi = gauss->coord4;
1133
1134- x1=*(xyz_list+3*0+0);
1135- x2=*(xyz_list+3*1+0);
1136- x3=*(xyz_list+3*2+0);
1137- x4=*(xyz_list+3*3+0);
1138- x5=*(xyz_list+3*4+0);
1139- x6=*(xyz_list+3*5+0);
1140+ x1=xyz_list[3*0+0];
1141+ x2=xyz_list[3*1+0];
1142+ x3=xyz_list[3*2+0];
1143+ x4=xyz_list[3*3+0];
1144+ x5=xyz_list[3*4+0];
1145+ x6=xyz_list[3*5+0];
1146
1147- y1=*(xyz_list+3*0+1);
1148- y2=*(xyz_list+3*1+1);
1149- y3=*(xyz_list+3*2+1);
1150- y4=*(xyz_list+3*3+1);
1151- y5=*(xyz_list+3*4+1);
1152- y6=*(xyz_list+3*5+1);
1153+ y1=xyz_list[3*0+1];
1154+ y2=xyz_list[3*1+1];
1155+ y3=xyz_list[3*2+1];
1156+ y4=xyz_list[3*3+1];
1157+ y5=xyz_list[3*4+1];
1158+ y6=xyz_list[3*5+1];
1159
1160- z1=*(xyz_list+3*0+2);
1161- z2=*(xyz_list+3*1+2);
1162- z3=*(xyz_list+3*2+2);
1163- z4=*(xyz_list+3*3+2);
1164- z5=*(xyz_list+3*4+2);
1165- z6=*(xyz_list+3*5+2);
1166+ z1=xyz_list[3*0+2];
1167+ z2=xyz_list[3*1+2];
1168+ z3=xyz_list[3*2+2];
1169+ z4=xyz_list[3*3+2];
1170+ z5=xyz_list[3*4+2];
1171+ z6=xyz_list[3*5+2];
1172
1173- *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
1174- *(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);
1175- *(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);
1176+ J[NDOF3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
1177+ 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);
1178+ 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);
1179
1180- *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
1181- *(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);
1182- *(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);
1183+ J[NDOF3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
1184+ 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);
1185+ 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);
1186
1187- *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
1188- *(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);
1189- *(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);
1190-
1191+ J[NDOF3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
1192+ 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);
1193+ 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);
1194 }
1195 /*}}}*/
1196 /*FUNCTION PentaRef::GetJacobianDeterminant {{{*/
1197@@ -1074,20 +1040,18 @@
1198 /*The Jacobian determinant is constant over the element, discard the gaussian points.
1199 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
1200
1201- IssmDouble x1,x2,x3,y1,y2,y3,z1,z2,z3;
1202+ IssmDouble x1=xyz_list[3*0+0];
1203+ IssmDouble y1=xyz_list[3*0+1];
1204+ IssmDouble z1=xyz_list[3*0+2];
1205+ IssmDouble x2=xyz_list[3*1+0];
1206+ IssmDouble y2=xyz_list[3*1+1];
1207+ IssmDouble z2=xyz_list[3*1+2];
1208+ IssmDouble x3=xyz_list[3*2+0];
1209+ IssmDouble y3=xyz_list[3*2+1];
1210+ IssmDouble z3=xyz_list[3*2+2];
1211
1212- x1=*(xyz_list+3*0+0);
1213- y1=*(xyz_list+3*0+1);
1214- z1=*(xyz_list+3*0+2);
1215- x2=*(xyz_list+3*1+0);
1216- y2=*(xyz_list+3*1+1);
1217- z2=*(xyz_list+3*1+2);
1218- x3=*(xyz_list+3*2+0);
1219- y3=*(xyz_list+3*2+1);
1220- z3=*(xyz_list+3*2+2);
1221-
1222 /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
1223- *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
1224+ *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);
1225 if(*Jdet<0) _error_("negative jacobian determinant!");
1226 }
1227 /*}}}*
1228@@ -1096,16 +1060,14 @@
1229 /*The Jacobian determinant is constant over the element, discard the gaussian points.
1230 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
1231
1232- IssmDouble x1,x2,y1,y2,z1,z2;
1233+ IssmDouble x1=xyz_list[3*0+0];
1234+ IssmDouble y1=xyz_list[3*0+1];
1235+ IssmDouble z1=xyz_list[3*0+2];
1236+ IssmDouble x2=xyz_list[3*1+0];
1237+ IssmDouble y2=xyz_list[3*1+1];
1238+ IssmDouble z2=xyz_list[3*1+2];
1239
1240- x1=*(xyz_list+3*0+0);
1241- y1=*(xyz_list+3*0+1);
1242- z1=*(xyz_list+3*0+2);
1243- x2=*(xyz_list+3*1+0);
1244- y2=*(xyz_list+3*1+1);
1245- z2=*(xyz_list+3*1+2);
1246-
1247- *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.) + pow(z2-z1,2.));
1248+ *Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
1249 if(*Jdet<0) _error_("negative jacobian determinant!");
1250
1251 }
1252@@ -1274,21 +1236,21 @@
1253 }
1254 /*}}}*/
1255 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{*/
1256-void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dh1dh7,IssmDouble* xyz_list, GaussPenta* gauss){
1257+void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dbasismini,IssmDouble* xyz_list, GaussPenta* gauss){
1258
1259 /*This routine returns the values of the nodal functions derivatives (with respect to the
1260 * actual coordinate system): */
1261
1262- IssmDouble dh1dh7_ref[3][NUMNODESMINI];
1263+ IssmDouble dbasismini_ref[3][NUMNODESMINI];
1264 IssmDouble Jinv[3][3];
1265
1266 /*Get derivative values with respect to parametric coordinate system: */
1267- GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss);
1268+ GetNodalFunctionsMINIDerivativesReference(&dbasismini_ref[0][0], gauss);
1269
1270 /*Get Jacobian invert: */
1271 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
1272
1273- /*Build dh1dh6:
1274+ /*Build dbasis:
1275 *
1276 * [dhi/dx]= Jinv'*[dhi/dr]
1277 * [dhi/dy] [dhi/ds]
1278@@ -1296,9 +1258,9 @@
1279 */
1280
1281 for(int i=0;i<NUMNODESMINI;i++){
1282- *(dh1dh7+NUMNODESMINI*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
1283- *(dh1dh7+NUMNODESMINI*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
1284- *(dh1dh7+NUMNODESMINI*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
1285+ *(dbasismini+NUMNODESMINI*0+i)=Jinv[0][0]*dbasismini_ref[0][i]+Jinv[0][1]*dbasismini_ref[1][i]+Jinv[0][2]*dbasismini_ref[2][i];
1286+ *(dbasismini+NUMNODESMINI*1+i)=Jinv[1][0]*dbasismini_ref[0][i]+Jinv[1][1]*dbasismini_ref[1][i]+Jinv[1][2]*dbasismini_ref[2][i];
1287+ *(dbasismini+NUMNODESMINI*2+i)=Jinv[2][0]*dbasismini_ref[0][i]+Jinv[2][1]*dbasismini_ref[1][i]+Jinv[2][2]*dbasismini_ref[2][i];
1288 }
1289
1290 }
1291@@ -1341,28 +1303,28 @@
1292 }
1293 /*}}}*/
1294 /*FUNCTION PentaRef::GetNodalFunctionsP1 {{{*/
1295-void PentaRef::GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss){
1296+void PentaRef::GetNodalFunctionsP1(IssmDouble* basis, GaussPenta* gauss){
1297 /*This routine returns the values of the nodal functions at the gaussian point.*/
1298
1299- l1l6[0]=gauss->coord1*(1-gauss->coord4)/2.0;
1300- l1l6[1]=gauss->coord2*(1-gauss->coord4)/2.0;
1301- l1l6[2]=gauss->coord3*(1-gauss->coord4)/2.0;
1302- l1l6[3]=gauss->coord1*(1+gauss->coord4)/2.0;
1303- l1l6[4]=gauss->coord2*(1+gauss->coord4)/2.0;
1304- l1l6[5]=gauss->coord3*(1+gauss->coord4)/2.0;
1305+ basis[0]=gauss->coord1*(1-gauss->coord4)/2.0;
1306+ basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
1307+ basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
1308+ basis[3]=gauss->coord1*(1+gauss->coord4)/2.0;
1309+ basis[4]=gauss->coord2*(1+gauss->coord4)/2.0;
1310+ basis[5]=gauss->coord3*(1+gauss->coord4)/2.0;
1311
1312 }
1313 /*}}}*/
1314 /*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{*/
1315-void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dh1dh6,IssmDouble* xyz_list, GaussPenta* gauss){
1316+void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
1317
1318 /*This routine returns the values of the nodal functions derivatives (with respect to the
1319 * actual coordinate system): */
1320- IssmDouble dh1dh6_ref[NDOF3][NUMNODESP1];
1321+ IssmDouble dbasis_ref[NDOF3][NUMNODESP1];
1322 IssmDouble Jinv[NDOF3][NDOF3];
1323
1324 /*Get derivative values with respect to parametric coordinate system: */
1325- GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss);
1326+ GetNodalFunctionsP1DerivativesReference(&dbasis_ref[0][0], gauss);
1327
1328 /*Get Jacobian invert: */
1329 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
1330@@ -1375,9 +1337,9 @@
1331 */
1332
1333 for (int i=0;i<NUMNODESP1;i++){
1334- *(dh1dh6+NUMNODESP1*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
1335- *(dh1dh6+NUMNODESP1*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
1336- *(dh1dh6+NUMNODESP1*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
1337+ *(dbasis+NUMNODESP1*0+i)=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i]+Jinv[0][2]*dbasis_ref[2][i];
1338+ *(dbasis+NUMNODESP1*1+i)=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i]+Jinv[1][2]*dbasis_ref[2][i];
1339+ *(dbasis+NUMNODESP1*2+i)=Jinv[2][0]*dbasis_ref[0][i]+Jinv[2][1]*dbasis_ref[1][i]+Jinv[2][2]*dbasis_ref[2][i];
1340 }
1341
1342 }
1343@@ -1467,34 +1429,36 @@
1344 /*P1 interpolation on Gauss point*/
1345
1346 /*intermediary*/
1347- IssmDouble l1l6[6];
1348+ IssmDouble basis[6];
1349
1350 /*nodal functions: */
1351- GetNodalFunctionsP1(&l1l6[0],gauss);
1352+ GetNodalFunctionsP1(&basis[0],gauss);
1353
1354 /*Assign output pointers:*/
1355- *pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5];
1356+ *pvalue=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2]+basis[3]*plist[3]+basis[4]*plist[4]+basis[5]*plist[5];
1357
1358 }
1359 /*}}}*/
1360 /*FUNCTION PentaRef::GetInputDerivativeValue{{{*/
1361 void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussPenta* gauss){
1362- /*From node 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:
1363+ /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
1364+ * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
1365+ * gaussian point specified by gauss_coord:
1366 * 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;
1367 * 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;
1368 * 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;
1369 *
1370 * p is a vector of size 3x1 already allocated.
1371 */
1372- IssmDouble dh1dh6[3][NUMNODESP1];
1373+ IssmDouble dbasis[3][NUMNODESP1];
1374
1375 /*Get nodal funnctions derivatives in actual coordinate system: */
1376- GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
1377+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
1378
1379 /*Assign output*/
1380- p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5];
1381- p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5];
1382- p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5];
1383+ p[0]=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2]+plist[3]*dbasis[0][3]+plist[4]*dbasis[0][4]+plist[5]*dbasis[0][5];
1384+ p[1]=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2]+plist[3]*dbasis[1][3]+plist[4]*dbasis[1][4]+plist[5]*dbasis[1][5];
1385+ p[2]=plist[0]*dbasis[2][0]+plist[1]*dbasis[2][1]+plist[2]*dbasis[2][2]+plist[3]*dbasis[2][3]+plist[4]*dbasis[2][4]+plist[5]*dbasis[2][5];
1386
1387 }
1388 /*}}}*/
1389Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp
1390===================================================================
1391--- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 15469)
1392+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 15470)
1393@@ -71,8 +71,8 @@
1394
1395 /*Build B: */
1396 for(int i=0;i<numnodes;i++){
1397- B[numnodes*0+i]=dbasis[0*numnodes+i];
1398- B[numnodes*1+i]=dbasis[1*numnodes+i];
1399+ B[numnodes*0+i] = dbasis[0*numnodes+i];
1400+ B[numnodes*1+i] = dbasis[1*numnodes+i];
1401 }
1402
1403 /*Clean-up*/
1404@@ -101,12 +101,12 @@
1405
1406 /*Build B: */
1407 for(int i=0;i<numnodes;i++){
1408- B[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i];
1409- B[NDOF2*numnodes*0+NDOF2*i+1]=0.;
1410- B[NDOF2*numnodes*1+NDOF2*i+0]=0.;
1411- B[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i];
1412- B[NDOF2*numnodes*2+NDOF2*i+0]=.5*dbasis[1*numnodes+i];
1413- B[NDOF2*numnodes*2+NDOF2*i+1]=.5*dbasis[0*numnodes+i];
1414+ B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
1415+ B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
1416+ B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
1417+ B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
1418+ B[NDOF2*numnodes*2+NDOF2*i+0] = .5*dbasis[1*numnodes+i];
1419+ B[NDOF2*numnodes*2+NDOF2*i+1] = .5*dbasis[0*numnodes+i];
1420 }
1421
1422 /*Clean-up*/
1423@@ -135,13 +135,13 @@
1424 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
1425
1426 /*Build B': */
1427- for (int i=0;i<numnodes;i++){
1428- B[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i];
1429- B[NDOF2*numnodes*0+NDOF2*i+1]=0.;
1430- B[NDOF2*numnodes*1+NDOF2*i+0]=0.;
1431- B[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i];
1432- B[NDOF2*numnodes*2+NDOF2*i+0]=0.5*dbasis[1*numnodes+i];
1433- B[NDOF2*numnodes*2+NDOF2*i+1]=0.5*dbasis[0*numnodes+i];
1434+ for(int i=0;i<numnodes;i++){
1435+ B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
1436+ B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
1437+ B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
1438+ B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
1439+ B[NDOF2*numnodes*2+NDOF2*i+0] = 0.5*dbasis[1*numnodes+i];
1440+ B[NDOF2*numnodes*2+NDOF2*i+1] = 0.5*dbasis[0*numnodes+i];
1441 }
1442
1443 /*Clean-up*/
1444@@ -220,9 +220,9 @@
1445 GetNodalFunctions(basis,gauss);
1446
1447 /*Build B: */
1448- for (int i=0;i<numnodes;i++){
1449- B[numnodes*0+i]=basis[i];
1450- B[numnodes*1+i]=basis[i];
1451+ for(int i=0;i<numnodes;i++){
1452+ B[numnodes*0+i] = basis[i];
1453+ B[numnodes*1+i] = basis[i];
1454 }
1455
1456 /*Clean-up*/
1457@@ -251,7 +251,7 @@
1458 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
1459
1460 /*Build B': */
1461- for (int i=0;i<numnodes;i++){
1462+ for(int i=0;i<numnodes;i++){
1463 Bprime[NDOF2*numnodes*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
1464 Bprime[NDOF2*numnodes*0+NDOF2*i+1] = dbasis[1*numnodes+i];
1465 Bprime[NDOF2*numnodes*1+NDOF2*i+0] = dbasis[0*numnodes+i];
1466@@ -287,14 +287,14 @@
1467
1468 /*Build Bprime: */
1469 for(int i=0;i<numnodes;i++){
1470- Bprime[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i];
1471- Bprime[NDOF2*numnodes*0+NDOF2*i+1]=0.;
1472- Bprime[NDOF2*numnodes*1+NDOF2*i+0]=0.;
1473- Bprime[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i];
1474- Bprime[NDOF2*numnodes*2+NDOF2*i+0]=dbasis[1*numnodes+i];
1475- Bprime[NDOF2*numnodes*2+NDOF2*i+1]=dbasis[0*numnodes+i];
1476- Bprime[NDOF2*numnodes*3+NDOF2*i+0]=dbasis[0*numnodes+i];
1477- Bprime[NDOF2*numnodes*3+NDOF2*i+1]=dbasis[1*numnodes+i];
1478+ Bprime[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
1479+ Bprime[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
1480+ Bprime[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
1481+ Bprime[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
1482+ Bprime[NDOF2*numnodes*2+NDOF2*i+0] = dbasis[1*numnodes+i];
1483+ Bprime[NDOF2*numnodes*2+NDOF2*i+1] = dbasis[0*numnodes+i];
1484+ Bprime[NDOF2*numnodes*3+NDOF2*i+0] = dbasis[0*numnodes+i];
1485+ Bprime[NDOF2*numnodes*3+NDOF2*i+1] = dbasis[1*numnodes+i];
1486 }
1487
1488 /*Clean-up*/
1489@@ -322,8 +322,8 @@
1490
1491 /*Build B': */
1492 for(int i=0;i<numnodes;i++){
1493- Bprime[numnodes*0+i]=dbasis[0*numnodes+i];
1494- Bprime[numnodes*1+i]=dbasis[1*numnodes+i];
1495+ Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
1496+ Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
1497 }
1498
1499 /*Clean-up*/
1500@@ -351,10 +351,10 @@
1501
1502 /*Build L: */
1503 for(int i=0;i<numnodes;i++){
1504- B[2*numnodes*0+2*i+0]=basis[i];
1505- B[2*numnodes*0+2*i+1]=0.;
1506- B[2*numnodes*1+2*i+0]=0.;
1507- B[2*numnodes*1+2*i+1]=basis[i];
1508+ B[2*numnodes*0+2*i+0] = basis[i];
1509+ B[2*numnodes*0+2*i+1] = 0.;
1510+ B[2*numnodes*1+2*i+0] = 0.;
1511+ B[2*numnodes*1+2*i+1] = basis[i];
1512 }
1513
1514 /*Clean-up*/
1515@@ -507,8 +507,8 @@
1516 * [dhi/dy] [dhi/ds]
1517 */
1518 for(int i=0;i<numnodes;i++){
1519- dbasis[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
1520- dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
1521+ dbasis[numnodes*0+i] = Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
1522+ dbasis[numnodes*1+i] = Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
1523 }
1524
1525 /*Clean up*/
1526@@ -526,34 +526,34 @@
1527 switch(this->element_type){
1528 case P1Enum: case P1DGEnum:
1529 /*Nodal function 1*/
1530- dbasis[NUMNODESP1*0+0]=-0.5;
1531- dbasis[NUMNODESP1*1+0]=-1.0/(2.0*SQRT3);
1532+ dbasis[NUMNODESP1*0+0] = -0.5;
1533+ dbasis[NUMNODESP1*1+0] = -1.0/(2.0*SQRT3);
1534 /*Nodal function 2*/
1535- dbasis[NUMNODESP1*0+1]=0.5;
1536- dbasis[NUMNODESP1*1+1]=-1.0/(2.0*SQRT3);
1537+ dbasis[NUMNODESP1*0+1] = 0.5;
1538+ dbasis[NUMNODESP1*1+1] = -1.0/(2.0*SQRT3);
1539 /*Nodal function 3*/
1540- dbasis[NUMNODESP1*0+2]=0;
1541- dbasis[NUMNODESP1*1+2]=1.0/SQRT3;
1542+ dbasis[NUMNODESP1*0+2] = 0;
1543+ dbasis[NUMNODESP1*1+2] = 1.0/SQRT3;
1544 return;
1545 case P2Enum:
1546 /*Nodal function 1*/
1547- dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5;
1548- dbasis[NUMNODESP2*1+0]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
1549+ dbasis[NUMNODESP2*0+0] = -2.*gauss->coord1 + 0.5;
1550+ dbasis[NUMNODESP2*1+0] = -2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
1551 /*Nodal function 2*/
1552- dbasis[NUMNODESP2*0+1]=+2.*gauss->coord2 - 0.5;
1553- dbasis[NUMNODESP2*1+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
1554+ dbasis[NUMNODESP2*0+1] = +2.*gauss->coord2 - 0.5;
1555+ dbasis[NUMNODESP2*1+1] = -2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
1556 /*Nodal function 3*/
1557- dbasis[NUMNODESP2*0+2]=0.;
1558- dbasis[NUMNODESP2*1+2]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
1559+ dbasis[NUMNODESP2*0+2] = 0.;
1560+ dbasis[NUMNODESP2*1+2] = +4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
1561 /*Nodal function 4*/
1562- dbasis[NUMNODESP2*0+3]=+2.*gauss->coord3;
1563- dbasis[NUMNODESP2*1+3]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
1564+ dbasis[NUMNODESP2*0+3] = +2.*gauss->coord3;
1565+ dbasis[NUMNODESP2*1+3] = +4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
1566 /*Nodal function 5*/
1567- dbasis[NUMNODESP2*0+4]=-2.*gauss->coord3;
1568- dbasis[NUMNODESP2*1+4]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
1569+ dbasis[NUMNODESP2*0+4] = -2.*gauss->coord3;
1570+ dbasis[NUMNODESP2*1+4] = +4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
1571 /*Nodal function 6*/
1572- dbasis[NUMNODESP2*0+5]=2.*(gauss->coord1-gauss->coord2);
1573- dbasis[NUMNODESP2*1+5]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
1574+ dbasis[NUMNODESP2*0+5] = 2.*(gauss->coord1-gauss->coord2);
1575+ dbasis[NUMNODESP2*1+5] = -2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
1576 return;
1577 default:
1578 _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
1579@@ -589,8 +589,8 @@
1580
1581 /*Assign values*/
1582 xDelete<IssmDouble>(dbasis);
1583- *(p+0)=dpx;
1584- *(p+1)=dpy;
1585+ p[0]=dpx;
1586+ p[1]=dpy;
1587
1588 }
1589 /*}}}*/
Note: See TracBrowser for help on using the repository browser.