« back to all changes in this revision
Viewing changes to test/regression/references/Optimization.h
- browse files at revision 8
- compare with another revision
- download diff
- download tarball
- view history from revision 8
- Committer: Bazaar Package Importer
- Author(s): Johannes Ring
- Date: 2010-07-01 19:54:32 UTC
- mfrom: (1.1.5 upstream)
- Revision ID: james.westby@ubuntu.com-20100701195432-xz3gw5nprdj79jcb
Tags: 0.9.3-1
* New upstream release.
* debian/control:
- Minor fix in Vcs fields.
- Bump Standards-Version to 3.9.0 (no changes needed).
- Update version for python-ufc, python-fiat, and python-ufl in
Depends field.
* Switch to dpkg-source 3.0 (quilt) format.
* Update debian/copyright and debian/copyright_hints.
* debian/control:
- Minor fix in Vcs fields.
- Bump Standards-Version to 3.9.0 (no changes needed).
- Update version for python-ufc, python-fiat, and python-ufl in
Depends field.
* Switch to dpkg-source 3.0 (quilt) format.
* Update debian/copyright and debian/copyright_hints.
- files added:
- bench/MassH1_2D_1.ufl
- bench/results
- debian/source
- files removed:
- bench/NavierStokes.form
- files modified:
- ChangeLog
added
removed
test/regression/references/Optimization.h
1
// This code conforms with the UFC specification version 1.4
2
// and was automatically generated by FFC version 0.9.2.
1
// This code conforms with the UFC specification version 1.4.1
2
// and was automatically generated by FFC version 0.9.3.
3
3
//
4
4
// This code was generated with the following parameters:
5
5
//
105
105
106
106
// Reset values.
107
107
*values = 0.00000000;
108
109
// Map degree of freedom to element degree of freedom
110
const unsigned int dof = i;
111
112
// Array of basisvalues.
113
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
114
115
// Declare helper variables.
116
unsigned int rr = 0;
117
unsigned int ss = 0;
118
unsigned int tt = 0;
119
double tmp5 = 0.00000000;
120
double tmp6 = 0.00000000;
121
double tmp7 = 0.00000000;
122
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
123
double tmp1 = (1.00000000 - Y)/2.00000000;
124
double tmp2 = tmp1*tmp1;
125
126
// Compute basisvalues.
127
basisvalues[0] = 1.00000000;
128
basisvalues[1] = tmp0;
129
for (unsigned int r = 1; r < 3; r++)
130
{
131
rr = (r + 1)*((r + 1) + 1)/2;
132
ss = r*(r + 1)/2;
133
tt = (r - 1)*((r - 1) + 1)/2;
134
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
135
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
136
}// end loop over 'r'
137
for (unsigned int r = 0; r < 3; r++)
138
{
139
rr = (r + 1)*(r + 1 + 1)/2 + 1;
140
ss = r*(r + 1)/2;
141
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
142
}// end loop over 'r'
143
for (unsigned int r = 0; r < 2; r++)
144
{
145
for (unsigned int s = 1; s < 3 - r; s++)
146
{
147
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
148
ss = (r + s)*(r + s + 1)/2 + s;
149
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
150
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
151
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
152
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
153
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
154
}// end loop over 's'
155
}// end loop over 'r'
156
for (unsigned int r = 0; r < 4; r++)
157
{
158
for (unsigned int s = 0; s < 4 - r; s++)
159
{
160
rr = (r + s)*(r + s + 1)/2 + s;
161
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
162
}// end loop over 's'
163
}// end loop over 'r'
164
165
// Table(s) of coefficients.
166
static const double coefficients0[10][10] = \
167
{{0.04714045, -0.02886751, -0.01666667, 0.07824608, 0.06060915, 0.03499271, -0.06013378, -0.05082232, -0.03936680, -0.02272843},
168
{0.04714045, 0.02886751, -0.01666667, 0.07824608, -0.06060915, 0.03499271, 0.06013378, -0.05082232, 0.03936680, -0.02272843},
169
{0.04714045, 0.00000000, 0.03333333, 0.00000000, 0.00000000, 0.10497813, 0.00000000, 0.00000000, 0.00000000, 0.09091373},
170
{0.10606602, 0.25980762, -0.15000000, 0.11736912, 0.06060915, -0.07873360, 0.00000000, 0.10164464, -0.13122266, 0.09091373},
171
{0.10606602, 0.00000000, 0.30000000, 0.00000000, 0.15152288, 0.02624453, 0.00000000, 0.00000000, 0.13122266, -0.13637059},
172
{0.10606602, -0.25980762, -0.15000000, 0.11736912, -0.06060915, -0.07873360, 0.00000000, 0.10164464, 0.13122266, 0.09091373},
173
{0.10606602, 0.00000000, 0.30000000, 0.00000000, -0.15152288, 0.02624453, 0.00000000, 0.00000000, -0.13122266, -0.13637059},
174
{0.10606602, -0.25980762, -0.15000000, -0.07824608, 0.09091373, 0.09622995, 0.18040134, 0.05082232, -0.01312227, -0.02272843},
175
{0.10606602, 0.25980762, -0.15000000, -0.07824608, -0.09091373, 0.09622995, -0.18040134, 0.05082232, 0.01312227, -0.02272843},
176
{0.63639610, 0.00000000, 0.00000000, -0.23473824, 0.00000000, -0.26244533, 0.00000000, -0.20328928, 0.00000000, 0.09091373}};
177
178
// Compute value(s).
179
for (unsigned int r = 0; r < 10; r++)
180
{
181
*values += coefficients0[dof][r]*basisvalues[r];
182
}// end loop over 'r'
108
switch (i)
109
{
110
case 0:
111
{
112
113
// Array of basisvalues.
114
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
115
116
// Declare helper variables.
117
unsigned int rr = 0;
118
unsigned int ss = 0;
119
unsigned int tt = 0;
120
double tmp5 = 0.00000000;
121
double tmp6 = 0.00000000;
122
double tmp7 = 0.00000000;
123
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
124
double tmp1 = (1.00000000 - Y)/2.00000000;
125
double tmp2 = tmp1*tmp1;
126
127
// Compute basisvalues.
128
basisvalues[0] = 1.00000000;
129
basisvalues[1] = tmp0;
130
for (unsigned int r = 1; r < 3; r++)
131
{
132
rr = (r + 1)*((r + 1) + 1)/2;
133
ss = r*(r + 1)/2;
134
tt = (r - 1)*((r - 1) + 1)/2;
135
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
136
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
137
}// end loop over 'r'
138
for (unsigned int r = 0; r < 3; r++)
139
{
140
rr = (r + 1)*(r + 1 + 1)/2 + 1;
141
ss = r*(r + 1)/2;
142
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
143
}// end loop over 'r'
144
for (unsigned int r = 0; r < 2; r++)
145
{
146
for (unsigned int s = 1; s < 3 - r; s++)
147
{
148
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
149
ss = (r + s)*(r + s + 1)/2 + s;
150
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
151
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
152
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
153
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
154
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
155
}// end loop over 's'
156
}// end loop over 'r'
157
for (unsigned int r = 0; r < 4; r++)
158
{
159
for (unsigned int s = 0; s < 4 - r; s++)
160
{
161
rr = (r + s)*(r + s + 1)/2 + s;
162
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
163
}// end loop over 's'
164
}// end loop over 'r'
165
166
// Table(s) of coefficients.
167
static const double coefficients0[10] = \
168
{0.04714045, -0.02886751, -0.01666667, 0.07824608, 0.06060915, 0.03499271, -0.06013378, -0.05082232, -0.03936680, -0.02272843};
169
170
// Compute value(s).
171
for (unsigned int r = 0; r < 10; r++)
172
{
173
*values += coefficients0[r]*basisvalues[r];
174
}// end loop over 'r'
175
break;
176
}
177
case 1:
178
{
179
180
// Array of basisvalues.
181
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
182
183
// Declare helper variables.
184
unsigned int rr = 0;
185
unsigned int ss = 0;
186
unsigned int tt = 0;
187
double tmp5 = 0.00000000;
188
double tmp6 = 0.00000000;
189
double tmp7 = 0.00000000;
190
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
191
double tmp1 = (1.00000000 - Y)/2.00000000;
192
double tmp2 = tmp1*tmp1;
193
194
// Compute basisvalues.
195
basisvalues[0] = 1.00000000;
196
basisvalues[1] = tmp0;
197
for (unsigned int r = 1; r < 3; r++)
198
{
199
rr = (r + 1)*((r + 1) + 1)/2;
200
ss = r*(r + 1)/2;
201
tt = (r - 1)*((r - 1) + 1)/2;
202
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
203
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
204
}// end loop over 'r'
205
for (unsigned int r = 0; r < 3; r++)
206
{
207
rr = (r + 1)*(r + 1 + 1)/2 + 1;
208
ss = r*(r + 1)/2;
209
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
210
}// end loop over 'r'
211
for (unsigned int r = 0; r < 2; r++)
212
{
213
for (unsigned int s = 1; s < 3 - r; s++)
214
{
215
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
216
ss = (r + s)*(r + s + 1)/2 + s;
217
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
218
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
219
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
220
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
221
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
222
}// end loop over 's'
223
}// end loop over 'r'
224
for (unsigned int r = 0; r < 4; r++)
225
{
226
for (unsigned int s = 0; s < 4 - r; s++)
227
{
228
rr = (r + s)*(r + s + 1)/2 + s;
229
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
230
}// end loop over 's'
231
}// end loop over 'r'
232
233
// Table(s) of coefficients.
234
static const double coefficients0[10] = \
235
{0.04714045, 0.02886751, -0.01666667, 0.07824608, -0.06060915, 0.03499271, 0.06013378, -0.05082232, 0.03936680, -0.02272843};
236
237
// Compute value(s).
238
for (unsigned int r = 0; r < 10; r++)
239
{
240
*values += coefficients0[r]*basisvalues[r];
241
}// end loop over 'r'
242
break;
243
}
244
case 2:
245
{
246
247
// Array of basisvalues.
248
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
249
250
// Declare helper variables.
251
unsigned int rr = 0;
252
unsigned int ss = 0;
253
unsigned int tt = 0;
254
double tmp5 = 0.00000000;
255
double tmp6 = 0.00000000;
256
double tmp7 = 0.00000000;
257
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
258
double tmp1 = (1.00000000 - Y)/2.00000000;
259
double tmp2 = tmp1*tmp1;
260
261
// Compute basisvalues.
262
basisvalues[0] = 1.00000000;
263
basisvalues[1] = tmp0;
264
for (unsigned int r = 1; r < 3; r++)
265
{
266
rr = (r + 1)*((r + 1) + 1)/2;
267
ss = r*(r + 1)/2;
268
tt = (r - 1)*((r - 1) + 1)/2;
269
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
270
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
271
}// end loop over 'r'
272
for (unsigned int r = 0; r < 3; r++)
273
{
274
rr = (r + 1)*(r + 1 + 1)/2 + 1;
275
ss = r*(r + 1)/2;
276
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
277
}// end loop over 'r'
278
for (unsigned int r = 0; r < 2; r++)
279
{
280
for (unsigned int s = 1; s < 3 - r; s++)
281
{
282
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
283
ss = (r + s)*(r + s + 1)/2 + s;
284
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
285
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
286
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
287
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
288
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
289
}// end loop over 's'
290
}// end loop over 'r'
291
for (unsigned int r = 0; r < 4; r++)
292
{
293
for (unsigned int s = 0; s < 4 - r; s++)
294
{
295
rr = (r + s)*(r + s + 1)/2 + s;
296
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
297
}// end loop over 's'
298
}// end loop over 'r'
299
300
// Table(s) of coefficients.
301
static const double coefficients0[10] = \
302
{0.04714045, 0.00000000, 0.03333333, 0.00000000, 0.00000000, 0.10497813, 0.00000000, 0.00000000, 0.00000000, 0.09091373};
303
304
// Compute value(s).
305
for (unsigned int r = 0; r < 10; r++)
306
{
307
*values += coefficients0[r]*basisvalues[r];
308
}// end loop over 'r'
309
break;
310
}
311
case 3:
312
{
313
314
// Array of basisvalues.
315
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
316
317
// Declare helper variables.
318
unsigned int rr = 0;
319
unsigned int ss = 0;
320
unsigned int tt = 0;
321
double tmp5 = 0.00000000;
322
double tmp6 = 0.00000000;
323
double tmp7 = 0.00000000;
324
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
325
double tmp1 = (1.00000000 - Y)/2.00000000;
326
double tmp2 = tmp1*tmp1;
327
328
// Compute basisvalues.
329
basisvalues[0] = 1.00000000;
330
basisvalues[1] = tmp0;
331
for (unsigned int r = 1; r < 3; r++)
332
{
333
rr = (r + 1)*((r + 1) + 1)/2;
334
ss = r*(r + 1)/2;
335
tt = (r - 1)*((r - 1) + 1)/2;
336
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
337
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
338
}// end loop over 'r'
339
for (unsigned int r = 0; r < 3; r++)
340
{
341
rr = (r + 1)*(r + 1 + 1)/2 + 1;
342
ss = r*(r + 1)/2;
343
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
344
}// end loop over 'r'
345
for (unsigned int r = 0; r < 2; r++)
346
{
347
for (unsigned int s = 1; s < 3 - r; s++)
348
{
349
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
350
ss = (r + s)*(r + s + 1)/2 + s;
351
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
352
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
353
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
354
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
355
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
356
}// end loop over 's'
357
}// end loop over 'r'
358
for (unsigned int r = 0; r < 4; r++)
359
{
360
for (unsigned int s = 0; s < 4 - r; s++)
361
{
362
rr = (r + s)*(r + s + 1)/2 + s;
363
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
364
}// end loop over 's'
365
}// end loop over 'r'
366
367
// Table(s) of coefficients.
368
static const double coefficients0[10] = \
369
{0.10606602, 0.25980762, -0.15000000, 0.11736912, 0.06060915, -0.07873360, 0.00000000, 0.10164464, -0.13122266, 0.09091373};
370
371
// Compute value(s).
372
for (unsigned int r = 0; r < 10; r++)
373
{
374
*values += coefficients0[r]*basisvalues[r];
375
}// end loop over 'r'
376
break;
377
}
378
case 4:
379
{
380
381
// Array of basisvalues.
382
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
383
384
// Declare helper variables.
385
unsigned int rr = 0;
386
unsigned int ss = 0;
387
unsigned int tt = 0;
388
double tmp5 = 0.00000000;
389
double tmp6 = 0.00000000;
390
double tmp7 = 0.00000000;
391
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
392
double tmp1 = (1.00000000 - Y)/2.00000000;
393
double tmp2 = tmp1*tmp1;
394
395
// Compute basisvalues.
396
basisvalues[0] = 1.00000000;
397
basisvalues[1] = tmp0;
398
for (unsigned int r = 1; r < 3; r++)
399
{
400
rr = (r + 1)*((r + 1) + 1)/2;
401
ss = r*(r + 1)/2;
402
tt = (r - 1)*((r - 1) + 1)/2;
403
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
404
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
405
}// end loop over 'r'
406
for (unsigned int r = 0; r < 3; r++)
407
{
408
rr = (r + 1)*(r + 1 + 1)/2 + 1;
409
ss = r*(r + 1)/2;
410
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
411
}// end loop over 'r'
412
for (unsigned int r = 0; r < 2; r++)
413
{
414
for (unsigned int s = 1; s < 3 - r; s++)
415
{
416
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
417
ss = (r + s)*(r + s + 1)/2 + s;
418
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
419
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
420
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
421
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
422
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
423
}// end loop over 's'
424
}// end loop over 'r'
425
for (unsigned int r = 0; r < 4; r++)
426
{
427
for (unsigned int s = 0; s < 4 - r; s++)
428
{
429
rr = (r + s)*(r + s + 1)/2 + s;
430
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
431
}// end loop over 's'
432
}// end loop over 'r'
433
434
// Table(s) of coefficients.
435
static const double coefficients0[10] = \
436
{0.10606602, 0.00000000, 0.30000000, 0.00000000, 0.15152288, 0.02624453, 0.00000000, 0.00000000, 0.13122266, -0.13637059};
437
438
// Compute value(s).
439
for (unsigned int r = 0; r < 10; r++)
440
{
441
*values += coefficients0[r]*basisvalues[r];
442
}// end loop over 'r'
443
break;
444
}
445
case 5:
446
{
447
448
// Array of basisvalues.
449
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
450
451
// Declare helper variables.
452
unsigned int rr = 0;
453
unsigned int ss = 0;
454
unsigned int tt = 0;
455
double tmp5 = 0.00000000;
456
double tmp6 = 0.00000000;
457
double tmp7 = 0.00000000;
458
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
459
double tmp1 = (1.00000000 - Y)/2.00000000;
460
double tmp2 = tmp1*tmp1;
461
462
// Compute basisvalues.
463
basisvalues[0] = 1.00000000;
464
basisvalues[1] = tmp0;
465
for (unsigned int r = 1; r < 3; r++)
466
{
467
rr = (r + 1)*((r + 1) + 1)/2;
468
ss = r*(r + 1)/2;
469
tt = (r - 1)*((r - 1) + 1)/2;
470
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
471
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
472
}// end loop over 'r'
473
for (unsigned int r = 0; r < 3; r++)
474
{
475
rr = (r + 1)*(r + 1 + 1)/2 + 1;
476
ss = r*(r + 1)/2;
477
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
478
}// end loop over 'r'
479
for (unsigned int r = 0; r < 2; r++)
480
{
481
for (unsigned int s = 1; s < 3 - r; s++)
482
{
483
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
484
ss = (r + s)*(r + s + 1)/2 + s;
485
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
486
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
487
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
488
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
489
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
490
}// end loop over 's'
491
}// end loop over 'r'
492
for (unsigned int r = 0; r < 4; r++)
493
{
494
for (unsigned int s = 0; s < 4 - r; s++)
495
{
496
rr = (r + s)*(r + s + 1)/2 + s;
497
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
498
}// end loop over 's'
499
}// end loop over 'r'
500
501
// Table(s) of coefficients.
502
static const double coefficients0[10] = \
503
{0.10606602, -0.25980762, -0.15000000, 0.11736912, -0.06060915, -0.07873360, 0.00000000, 0.10164464, 0.13122266, 0.09091373};
504
505
// Compute value(s).
506
for (unsigned int r = 0; r < 10; r++)
507
{
508
*values += coefficients0[r]*basisvalues[r];
509
}// end loop over 'r'
510
break;
511
}
512
case 6:
513
{
514
515
// Array of basisvalues.
516
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
517
518
// Declare helper variables.
519
unsigned int rr = 0;
520
unsigned int ss = 0;
521
unsigned int tt = 0;
522
double tmp5 = 0.00000000;
523
double tmp6 = 0.00000000;
524
double tmp7 = 0.00000000;
525
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
526
double tmp1 = (1.00000000 - Y)/2.00000000;
527
double tmp2 = tmp1*tmp1;
528
529
// Compute basisvalues.
530
basisvalues[0] = 1.00000000;
531
basisvalues[1] = tmp0;
532
for (unsigned int r = 1; r < 3; r++)
533
{
534
rr = (r + 1)*((r + 1) + 1)/2;
535
ss = r*(r + 1)/2;
536
tt = (r - 1)*((r - 1) + 1)/2;
537
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
538
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
539
}// end loop over 'r'
540
for (unsigned int r = 0; r < 3; r++)
541
{
542
rr = (r + 1)*(r + 1 + 1)/2 + 1;
543
ss = r*(r + 1)/2;
544
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
545
}// end loop over 'r'
546
for (unsigned int r = 0; r < 2; r++)
547
{
548
for (unsigned int s = 1; s < 3 - r; s++)
549
{
550
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
551
ss = (r + s)*(r + s + 1)/2 + s;
552
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
553
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
554
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
555
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
556
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
557
}// end loop over 's'
558
}// end loop over 'r'
559
for (unsigned int r = 0; r < 4; r++)
560
{
561
for (unsigned int s = 0; s < 4 - r; s++)
562
{
563
rr = (r + s)*(r + s + 1)/2 + s;
564
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
565
}// end loop over 's'
566
}// end loop over 'r'
567
568
// Table(s) of coefficients.
569
static const double coefficients0[10] = \
570
{0.10606602, 0.00000000, 0.30000000, 0.00000000, -0.15152288, 0.02624453, 0.00000000, 0.00000000, -0.13122266, -0.13637059};
571
572
// Compute value(s).
573
for (unsigned int r = 0; r < 10; r++)
574
{
575
*values += coefficients0[r]*basisvalues[r];
576
}// end loop over 'r'
577
break;
578
}
579
case 7:
580
{
581
582
// Array of basisvalues.
583
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
584
585
// Declare helper variables.
586
unsigned int rr = 0;
587
unsigned int ss = 0;
588
unsigned int tt = 0;
589
double tmp5 = 0.00000000;
590
double tmp6 = 0.00000000;
591
double tmp7 = 0.00000000;
592
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
593
double tmp1 = (1.00000000 - Y)/2.00000000;
594
double tmp2 = tmp1*tmp1;
595
596
// Compute basisvalues.
597
basisvalues[0] = 1.00000000;
598
basisvalues[1] = tmp0;
599
for (unsigned int r = 1; r < 3; r++)
600
{
601
rr = (r + 1)*((r + 1) + 1)/2;
602
ss = r*(r + 1)/2;
603
tt = (r - 1)*((r - 1) + 1)/2;
604
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
605
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
606
}// end loop over 'r'
607
for (unsigned int r = 0; r < 3; r++)
608
{
609
rr = (r + 1)*(r + 1 + 1)/2 + 1;
610
ss = r*(r + 1)/2;
611
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
612
}// end loop over 'r'
613
for (unsigned int r = 0; r < 2; r++)
614
{
615
for (unsigned int s = 1; s < 3 - r; s++)
616
{
617
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
618
ss = (r + s)*(r + s + 1)/2 + s;
619
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
620
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
621
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
622
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
623
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
624
}// end loop over 's'
625
}// end loop over 'r'
626
for (unsigned int r = 0; r < 4; r++)
627
{
628
for (unsigned int s = 0; s < 4 - r; s++)
629
{
630
rr = (r + s)*(r + s + 1)/2 + s;
631
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
632
}// end loop over 's'
633
}// end loop over 'r'
634
635
// Table(s) of coefficients.
636
static const double coefficients0[10] = \
637
{0.10606602, -0.25980762, -0.15000000, -0.07824608, 0.09091373, 0.09622995, 0.18040134, 0.05082232, -0.01312227, -0.02272843};
638
639
// Compute value(s).
640
for (unsigned int r = 0; r < 10; r++)
641
{
642
*values += coefficients0[r]*basisvalues[r];
643
}// end loop over 'r'
644
break;
645
}
646
case 8:
647
{
648
649
// Array of basisvalues.
650
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
651
652
// Declare helper variables.
653
unsigned int rr = 0;
654
unsigned int ss = 0;
655
unsigned int tt = 0;
656
double tmp5 = 0.00000000;
657
double tmp6 = 0.00000000;
658
double tmp7 = 0.00000000;
659
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
660
double tmp1 = (1.00000000 - Y)/2.00000000;
661
double tmp2 = tmp1*tmp1;
662
663
// Compute basisvalues.
664
basisvalues[0] = 1.00000000;
665
basisvalues[1] = tmp0;
666
for (unsigned int r = 1; r < 3; r++)
667
{
668
rr = (r + 1)*((r + 1) + 1)/2;
669
ss = r*(r + 1)/2;
670
tt = (r - 1)*((r - 1) + 1)/2;
671
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
672
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
673
}// end loop over 'r'
674
for (unsigned int r = 0; r < 3; r++)
675
{
676
rr = (r + 1)*(r + 1 + 1)/2 + 1;
677
ss = r*(r + 1)/2;
678
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
679
}// end loop over 'r'
680
for (unsigned int r = 0; r < 2; r++)
681
{
682
for (unsigned int s = 1; s < 3 - r; s++)
683
{
684
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
685
ss = (r + s)*(r + s + 1)/2 + s;
686
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
687
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
688
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
689
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
690
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
691
}// end loop over 's'
692
}// end loop over 'r'
693
for (unsigned int r = 0; r < 4; r++)
694
{
695
for (unsigned int s = 0; s < 4 - r; s++)
696
{
697
rr = (r + s)*(r + s + 1)/2 + s;
698
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
699
}// end loop over 's'
700
}// end loop over 'r'
701
702
// Table(s) of coefficients.
703
static const double coefficients0[10] = \
704
{0.10606602, 0.25980762, -0.15000000, -0.07824608, -0.09091373, 0.09622995, -0.18040134, 0.05082232, 0.01312227, -0.02272843};
705
706
// Compute value(s).
707
for (unsigned int r = 0; r < 10; r++)
708
{
709
*values += coefficients0[r]*basisvalues[r];
710
}// end loop over 'r'
711
break;
712
}
713
case 9:
714
{
715
716
// Array of basisvalues.
717
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
718
719
// Declare helper variables.
720
unsigned int rr = 0;
721
unsigned int ss = 0;
722
unsigned int tt = 0;
723
double tmp5 = 0.00000000;
724
double tmp6 = 0.00000000;
725
double tmp7 = 0.00000000;
726
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
727
double tmp1 = (1.00000000 - Y)/2.00000000;
728
double tmp2 = tmp1*tmp1;
729
730
// Compute basisvalues.
731
basisvalues[0] = 1.00000000;
732
basisvalues[1] = tmp0;
733
for (unsigned int r = 1; r < 3; r++)
734
{
735
rr = (r + 1)*((r + 1) + 1)/2;
736
ss = r*(r + 1)/2;
737
tt = (r - 1)*((r - 1) + 1)/2;
738
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
739
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
740
}// end loop over 'r'
741
for (unsigned int r = 0; r < 3; r++)
742
{
743
rr = (r + 1)*(r + 1 + 1)/2 + 1;
744
ss = r*(r + 1)/2;
745
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
746
}// end loop over 'r'
747
for (unsigned int r = 0; r < 2; r++)
748
{
749
for (unsigned int s = 1; s < 3 - r; s++)
750
{
751
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
752
ss = (r + s)*(r + s + 1)/2 + s;
753
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
754
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
755
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
756
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
757
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
758
}// end loop over 's'
759
}// end loop over 'r'
760
for (unsigned int r = 0; r < 4; r++)
761
{
762
for (unsigned int s = 0; s < 4 - r; s++)
763
{
764
rr = (r + s)*(r + s + 1)/2 + s;
765
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
766
}// end loop over 's'
767
}// end loop over 'r'
768
769
// Table(s) of coefficients.
770
static const double coefficients0[10] = \
771
{0.63639610, 0.00000000, 0.00000000, -0.23473824, 0.00000000, -0.26244533, 0.00000000, -0.20328928, 0.00000000, 0.09091373};
772
773
// Compute value(s).
774
for (unsigned int r = 0; r < 10; r++)
775
{
776
*values += coefficients0[r]*basisvalues[r];
777
}// end loop over 'r'
778
break;
779
}
780
}
781
183
782
}
184
783
185
784
/// Evaluate all basis functions at given point in cell
295
894
values[r] = 0.00000000;
296
895
}// end loop over 'r'
297
896
298
// Map degree of freedom to element degree of freedom
299
const unsigned int dof = i;
300
301
// Array of basisvalues.
302
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
303
304
// Declare helper variables.
305
unsigned int rr = 0;
306
unsigned int ss = 0;
307
unsigned int tt = 0;
308
double tmp5 = 0.00000000;
309
double tmp6 = 0.00000000;
310
double tmp7 = 0.00000000;
311
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
312
double tmp1 = (1.00000000 - Y)/2.00000000;
313
double tmp2 = tmp1*tmp1;
314
315
// Compute basisvalues.
316
basisvalues[0] = 1.00000000;
317
basisvalues[1] = tmp0;
318
for (unsigned int r = 1; r < 3; r++)
319
{
320
rr = (r + 1)*((r + 1) + 1)/2;
321
ss = r*(r + 1)/2;
322
tt = (r - 1)*((r - 1) + 1)/2;
323
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
324
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
325
}// end loop over 'r'
326
for (unsigned int r = 0; r < 3; r++)
327
{
328
rr = (r + 1)*(r + 1 + 1)/2 + 1;
329
ss = r*(r + 1)/2;
330
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
331
}// end loop over 'r'
332
for (unsigned int r = 0; r < 2; r++)
333
{
334
for (unsigned int s = 1; s < 3 - r; s++)
335
{
336
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
337
ss = (r + s)*(r + s + 1)/2 + s;
338
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
339
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
340
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
341
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
342
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
343
}// end loop over 's'
344
}// end loop over 'r'
345
for (unsigned int r = 0; r < 4; r++)
346
{
347
for (unsigned int s = 0; s < 4 - r; s++)
348
{
349
rr = (r + s)*(r + s + 1)/2 + s;
350
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
351
}// end loop over 's'
352
}// end loop over 'r'
353
354
// Table(s) of coefficients.
355
static const double coefficients0[10][10] = \
356
{{0.04714045, -0.02886751, -0.01666667, 0.07824608, 0.06060915, 0.03499271, -0.06013378, -0.05082232, -0.03936680, -0.02272843},
357
{0.04714045, 0.02886751, -0.01666667, 0.07824608, -0.06060915, 0.03499271, 0.06013378, -0.05082232, 0.03936680, -0.02272843},
358
{0.04714045, 0.00000000, 0.03333333, 0.00000000, 0.00000000, 0.10497813, 0.00000000, 0.00000000, 0.00000000, 0.09091373},
359
{0.10606602, 0.25980762, -0.15000000, 0.11736912, 0.06060915, -0.07873360, 0.00000000, 0.10164464, -0.13122266, 0.09091373},
360
{0.10606602, 0.00000000, 0.30000000, 0.00000000, 0.15152288, 0.02624453, 0.00000000, 0.00000000, 0.13122266, -0.13637059},
361
{0.10606602, -0.25980762, -0.15000000, 0.11736912, -0.06060915, -0.07873360, 0.00000000, 0.10164464, 0.13122266, 0.09091373},
362
{0.10606602, 0.00000000, 0.30000000, 0.00000000, -0.15152288, 0.02624453, 0.00000000, 0.00000000, -0.13122266, -0.13637059},
363
{0.10606602, -0.25980762, -0.15000000, -0.07824608, 0.09091373, 0.09622995, 0.18040134, 0.05082232, -0.01312227, -0.02272843},
364
{0.10606602, 0.25980762, -0.15000000, -0.07824608, -0.09091373, 0.09622995, -0.18040134, 0.05082232, 0.01312227, -0.02272843},
365
{0.63639610, 0.00000000, 0.00000000, -0.23473824, 0.00000000, -0.26244533, 0.00000000, -0.20328928, 0.00000000, 0.09091373}};
366
367
// Tables of derivatives of the polynomial base (transpose).
368
static const double dmats0[10][10] = \
369
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
370
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
371
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
372
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
373
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
374
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
375
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
376
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
377
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
378
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
379
380
static const double dmats1[10][10] = \
381
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
382
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
383
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
384
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
385
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
386
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
387
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
388
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
389
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
390
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
391
392
// Compute reference derivatives.
393
// Declare pointer to array of derivatives on FIAT element.
394
double *derivatives = new double[num_derivatives];
395
for (unsigned int r = 0; r < num_derivatives; r++)
396
{
397
derivatives[r] = 0.00000000;
398
}// end loop over 'r'
399
400
// Declare derivative matrix (of polynomial basis).
401
double dmats[10][10] = \
402
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
403
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
404
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
405
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
406
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
407
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
408
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
409
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
410
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
411
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
412
413
// Declare (auxiliary) derivative matrix (of polynomial basis).
414
double dmats_old[10][10] = \
415
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
416
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
417
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
418
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
419
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
420
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
421
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
422
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
423
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
424
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
425
426
// Loop possible derivatives.
427
for (unsigned int r = 0; r < num_derivatives; r++)
428
{
429
// Resetting dmats values to compute next derivative.
430
for (unsigned int t = 0; t < 10; t++)
431
{
432
for (unsigned int u = 0; u < 10; u++)
433
{
434
dmats[t][u] = 0.00000000;
435
if (t == u)
436
{
437
dmats[t][u] = 1.00000000;
438
}
439
440
}// end loop over 'u'
441
}// end loop over 't'
442
443
// Looping derivative order to generate dmats.
444
for (unsigned int s = 0; s < n; s++)
445
{
446
// Updating dmats_old with new values and resetting dmats.
447
for (unsigned int t = 0; t < 10; t++)
448
{
449
for (unsigned int u = 0; u < 10; u++)
450
{
451
dmats_old[t][u] = dmats[t][u];
452
dmats[t][u] = 0.00000000;
453
}// end loop over 'u'
454
}// end loop over 't'
455
456
// Update dmats using an inner product.
457
if (combinations[r][s] == 0)
458
{
459
for (unsigned int t = 0; t < 10; t++)
460
{
461
for (unsigned int u = 0; u < 10; u++)
462
{
463
for (unsigned int tu = 0; tu < 10; tu++)
464
{
465
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
466
}// end loop over 'tu'
467
}// end loop over 'u'
468
}// end loop over 't'
469
}
470
471
if (combinations[r][s] == 1)
472
{
473
for (unsigned int t = 0; t < 10; t++)
474
{
475
for (unsigned int u = 0; u < 10; u++)
476
{
477
for (unsigned int tu = 0; tu < 10; tu++)
478
{
479
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
480
}// end loop over 'tu'
481
}// end loop over 'u'
482
}// end loop over 't'
483
}
484
485
}// end loop over 's'
486
for (unsigned int s = 0; s < 10; s++)
487
{
488
for (unsigned int t = 0; t < 10; t++)
489
{
490
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
491
}// end loop over 't'
492
}// end loop over 's'
493
}// end loop over 'r'
494
495
// Transform derivatives back to physical element
496
for (unsigned int r = 0; r < num_derivatives; r++)
497
{
498
for (unsigned int s = 0; s < num_derivatives; s++)
499
{
500
values[r] += transform[r][s]*derivatives[s];
501
}// end loop over 's'
502
}// end loop over 'r'
503
504
// Delete pointer to array of derivatives on FIAT element
505
delete [] derivatives;
506
507
// Delete pointer to array of combinations of derivatives and transform
508
for (unsigned int r = 0; r < num_derivatives; r++)
509
{
510
delete [] combinations[r];
511
}// end loop over 'r'
512
delete [] combinations;
513
for (unsigned int r = 0; r < num_derivatives; r++)
514
{
515
delete [] transform[r];
516
}// end loop over 'r'
517
delete [] transform;
897
switch (i)
898
{
899
case 0:
900
{
901
902
// Array of basisvalues.
903
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
904
905
// Declare helper variables.
906
unsigned int rr = 0;
907
unsigned int ss = 0;
908
unsigned int tt = 0;
909
double tmp5 = 0.00000000;
910
double tmp6 = 0.00000000;
911
double tmp7 = 0.00000000;
912
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
913
double tmp1 = (1.00000000 - Y)/2.00000000;
914
double tmp2 = tmp1*tmp1;
915
916
// Compute basisvalues.
917
basisvalues[0] = 1.00000000;
918
basisvalues[1] = tmp0;
919
for (unsigned int r = 1; r < 3; r++)
920
{
921
rr = (r + 1)*((r + 1) + 1)/2;
922
ss = r*(r + 1)/2;
923
tt = (r - 1)*((r - 1) + 1)/2;
924
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
925
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
926
}// end loop over 'r'
927
for (unsigned int r = 0; r < 3; r++)
928
{
929
rr = (r + 1)*(r + 1 + 1)/2 + 1;
930
ss = r*(r + 1)/2;
931
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
932
}// end loop over 'r'
933
for (unsigned int r = 0; r < 2; r++)
934
{
935
for (unsigned int s = 1; s < 3 - r; s++)
936
{
937
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
938
ss = (r + s)*(r + s + 1)/2 + s;
939
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
940
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
941
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
942
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
943
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
944
}// end loop over 's'
945
}// end loop over 'r'
946
for (unsigned int r = 0; r < 4; r++)
947
{
948
for (unsigned int s = 0; s < 4 - r; s++)
949
{
950
rr = (r + s)*(r + s + 1)/2 + s;
951
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
952
}// end loop over 's'
953
}// end loop over 'r'
954
955
// Table(s) of coefficients.
956
static const double coefficients0[10] = \
957
{0.04714045, -0.02886751, -0.01666667, 0.07824608, 0.06060915, 0.03499271, -0.06013378, -0.05082232, -0.03936680, -0.02272843};
958
959
// Tables of derivatives of the polynomial base (transpose).
960
static const double dmats0[10][10] = \
961
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
962
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
963
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
964
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
965
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
966
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
967
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
968
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
969
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
970
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
971
972
static const double dmats1[10][10] = \
973
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
974
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
975
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
976
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
977
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
978
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
979
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
980
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
981
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
982
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
983
984
// Compute reference derivatives.
985
// Declare pointer to array of derivatives on FIAT element.
986
double *derivatives = new double[num_derivatives];
987
for (unsigned int r = 0; r < num_derivatives; r++)
988
{
989
derivatives[r] = 0.00000000;
990
}// end loop over 'r'
991
992
// Declare derivative matrix (of polynomial basis).
993
double dmats[10][10] = \
994
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
995
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
996
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
997
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
998
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
999
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1000
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1001
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1002
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1003
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1004
1005
// Declare (auxiliary) derivative matrix (of polynomial basis).
1006
double dmats_old[10][10] = \
1007
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1008
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1009
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1010
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1011
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1012
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1013
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1014
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1015
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1016
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1017
1018
// Loop possible derivatives.
1019
for (unsigned int r = 0; r < num_derivatives; r++)
1020
{
1021
// Resetting dmats values to compute next derivative.
1022
for (unsigned int t = 0; t < 10; t++)
1023
{
1024
for (unsigned int u = 0; u < 10; u++)
1025
{
1026
dmats[t][u] = 0.00000000;
1027
if (t == u)
1028
{
1029
dmats[t][u] = 1.00000000;
1030
}
1031
1032
}// end loop over 'u'
1033
}// end loop over 't'
1034
1035
// Looping derivative order to generate dmats.
1036
for (unsigned int s = 0; s < n; s++)
1037
{
1038
// Updating dmats_old with new values and resetting dmats.
1039
for (unsigned int t = 0; t < 10; t++)
1040
{
1041
for (unsigned int u = 0; u < 10; u++)
1042
{
1043
dmats_old[t][u] = dmats[t][u];
1044
dmats[t][u] = 0.00000000;
1045
}// end loop over 'u'
1046
}// end loop over 't'
1047
1048
// Update dmats using an inner product.
1049
if (combinations[r][s] == 0)
1050
{
1051
for (unsigned int t = 0; t < 10; t++)
1052
{
1053
for (unsigned int u = 0; u < 10; u++)
1054
{
1055
for (unsigned int tu = 0; tu < 10; tu++)
1056
{
1057
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1058
}// end loop over 'tu'
1059
}// end loop over 'u'
1060
}// end loop over 't'
1061
}
1062
1063
if (combinations[r][s] == 1)
1064
{
1065
for (unsigned int t = 0; t < 10; t++)
1066
{
1067
for (unsigned int u = 0; u < 10; u++)
1068
{
1069
for (unsigned int tu = 0; tu < 10; tu++)
1070
{
1071
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1072
}// end loop over 'tu'
1073
}// end loop over 'u'
1074
}// end loop over 't'
1075
}
1076
1077
}// end loop over 's'
1078
for (unsigned int s = 0; s < 10; s++)
1079
{
1080
for (unsigned int t = 0; t < 10; t++)
1081
{
1082
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1083
}// end loop over 't'
1084
}// end loop over 's'
1085
}// end loop over 'r'
1086
1087
// Transform derivatives back to physical element
1088
for (unsigned int r = 0; r < num_derivatives; r++)
1089
{
1090
for (unsigned int s = 0; s < num_derivatives; s++)
1091
{
1092
values[r] += transform[r][s]*derivatives[s];
1093
}// end loop over 's'
1094
}// end loop over 'r'
1095
1096
// Delete pointer to array of derivatives on FIAT element
1097
delete [] derivatives;
1098
1099
// Delete pointer to array of combinations of derivatives and transform
1100
for (unsigned int r = 0; r < num_derivatives; r++)
1101
{
1102
delete [] combinations[r];
1103
}// end loop over 'r'
1104
delete [] combinations;
1105
for (unsigned int r = 0; r < num_derivatives; r++)
1106
{
1107
delete [] transform[r];
1108
}// end loop over 'r'
1109
delete [] transform;
1110
break;
1111
}
1112
case 1:
1113
{
1114
1115
// Array of basisvalues.
1116
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1117
1118
// Declare helper variables.
1119
unsigned int rr = 0;
1120
unsigned int ss = 0;
1121
unsigned int tt = 0;
1122
double tmp5 = 0.00000000;
1123
double tmp6 = 0.00000000;
1124
double tmp7 = 0.00000000;
1125
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1126
double tmp1 = (1.00000000 - Y)/2.00000000;
1127
double tmp2 = tmp1*tmp1;
1128
1129
// Compute basisvalues.
1130
basisvalues[0] = 1.00000000;
1131
basisvalues[1] = tmp0;
1132
for (unsigned int r = 1; r < 3; r++)
1133
{
1134
rr = (r + 1)*((r + 1) + 1)/2;
1135
ss = r*(r + 1)/2;
1136
tt = (r - 1)*((r - 1) + 1)/2;
1137
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1138
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1139
}// end loop over 'r'
1140
for (unsigned int r = 0; r < 3; r++)
1141
{
1142
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1143
ss = r*(r + 1)/2;
1144
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1145
}// end loop over 'r'
1146
for (unsigned int r = 0; r < 2; r++)
1147
{
1148
for (unsigned int s = 1; s < 3 - r; s++)
1149
{
1150
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1151
ss = (r + s)*(r + s + 1)/2 + s;
1152
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1153
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1154
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1155
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1156
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1157
}// end loop over 's'
1158
}// end loop over 'r'
1159
for (unsigned int r = 0; r < 4; r++)
1160
{
1161
for (unsigned int s = 0; s < 4 - r; s++)
1162
{
1163
rr = (r + s)*(r + s + 1)/2 + s;
1164
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1165
}// end loop over 's'
1166
}// end loop over 'r'
1167
1168
// Table(s) of coefficients.
1169
static const double coefficients0[10] = \
1170
{0.04714045, 0.02886751, -0.01666667, 0.07824608, -0.06060915, 0.03499271, 0.06013378, -0.05082232, 0.03936680, -0.02272843};
1171
1172
// Tables of derivatives of the polynomial base (transpose).
1173
static const double dmats0[10][10] = \
1174
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1175
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1176
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1177
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1178
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1179
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1180
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1181
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1182
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1183
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1184
1185
static const double dmats1[10][10] = \
1186
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1187
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1188
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1189
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1190
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1191
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1192
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1193
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1194
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1195
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1196
1197
// Compute reference derivatives.
1198
// Declare pointer to array of derivatives on FIAT element.
1199
double *derivatives = new double[num_derivatives];
1200
for (unsigned int r = 0; r < num_derivatives; r++)
1201
{
1202
derivatives[r] = 0.00000000;
1203
}// end loop over 'r'
1204
1205
// Declare derivative matrix (of polynomial basis).
1206
double dmats[10][10] = \
1207
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1208
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1209
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1210
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1211
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1212
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1213
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1214
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1215
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1216
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1217
1218
// Declare (auxiliary) derivative matrix (of polynomial basis).
1219
double dmats_old[10][10] = \
1220
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1221
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1222
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1223
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1224
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1225
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1226
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1227
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1228
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1229
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1230
1231
// Loop possible derivatives.
1232
for (unsigned int r = 0; r < num_derivatives; r++)
1233
{
1234
// Resetting dmats values to compute next derivative.
1235
for (unsigned int t = 0; t < 10; t++)
1236
{
1237
for (unsigned int u = 0; u < 10; u++)
1238
{
1239
dmats[t][u] = 0.00000000;
1240
if (t == u)
1241
{
1242
dmats[t][u] = 1.00000000;
1243
}
1244
1245
}// end loop over 'u'
1246
}// end loop over 't'
1247
1248
// Looping derivative order to generate dmats.
1249
for (unsigned int s = 0; s < n; s++)
1250
{
1251
// Updating dmats_old with new values and resetting dmats.
1252
for (unsigned int t = 0; t < 10; t++)
1253
{
1254
for (unsigned int u = 0; u < 10; u++)
1255
{
1256
dmats_old[t][u] = dmats[t][u];
1257
dmats[t][u] = 0.00000000;
1258
}// end loop over 'u'
1259
}// end loop over 't'
1260
1261
// Update dmats using an inner product.
1262
if (combinations[r][s] == 0)
1263
{
1264
for (unsigned int t = 0; t < 10; t++)
1265
{
1266
for (unsigned int u = 0; u < 10; u++)
1267
{
1268
for (unsigned int tu = 0; tu < 10; tu++)
1269
{
1270
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1271
}// end loop over 'tu'
1272
}// end loop over 'u'
1273
}// end loop over 't'
1274
}
1275
1276
if (combinations[r][s] == 1)
1277
{
1278
for (unsigned int t = 0; t < 10; t++)
1279
{
1280
for (unsigned int u = 0; u < 10; u++)
1281
{
1282
for (unsigned int tu = 0; tu < 10; tu++)
1283
{
1284
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1285
}// end loop over 'tu'
1286
}// end loop over 'u'
1287
}// end loop over 't'
1288
}
1289
1290
}// end loop over 's'
1291
for (unsigned int s = 0; s < 10; s++)
1292
{
1293
for (unsigned int t = 0; t < 10; t++)
1294
{
1295
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1296
}// end loop over 't'
1297
}// end loop over 's'
1298
}// end loop over 'r'
1299
1300
// Transform derivatives back to physical element
1301
for (unsigned int r = 0; r < num_derivatives; r++)
1302
{
1303
for (unsigned int s = 0; s < num_derivatives; s++)
1304
{
1305
values[r] += transform[r][s]*derivatives[s];
1306
}// end loop over 's'
1307
}// end loop over 'r'
1308
1309
// Delete pointer to array of derivatives on FIAT element
1310
delete [] derivatives;
1311
1312
// Delete pointer to array of combinations of derivatives and transform
1313
for (unsigned int r = 0; r < num_derivatives; r++)
1314
{
1315
delete [] combinations[r];
1316
}// end loop over 'r'
1317
delete [] combinations;
1318
for (unsigned int r = 0; r < num_derivatives; r++)
1319
{
1320
delete [] transform[r];
1321
}// end loop over 'r'
1322
delete [] transform;
1323
break;
1324
}
1325
case 2:
1326
{
1327
1328
// Array of basisvalues.
1329
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1330
1331
// Declare helper variables.
1332
unsigned int rr = 0;
1333
unsigned int ss = 0;
1334
unsigned int tt = 0;
1335
double tmp5 = 0.00000000;
1336
double tmp6 = 0.00000000;
1337
double tmp7 = 0.00000000;
1338
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1339
double tmp1 = (1.00000000 - Y)/2.00000000;
1340
double tmp2 = tmp1*tmp1;
1341
1342
// Compute basisvalues.
1343
basisvalues[0] = 1.00000000;
1344
basisvalues[1] = tmp0;
1345
for (unsigned int r = 1; r < 3; r++)
1346
{
1347
rr = (r + 1)*((r + 1) + 1)/2;
1348
ss = r*(r + 1)/2;
1349
tt = (r - 1)*((r - 1) + 1)/2;
1350
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1351
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1352
}// end loop over 'r'
1353
for (unsigned int r = 0; r < 3; r++)
1354
{
1355
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1356
ss = r*(r + 1)/2;
1357
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1358
}// end loop over 'r'
1359
for (unsigned int r = 0; r < 2; r++)
1360
{
1361
for (unsigned int s = 1; s < 3 - r; s++)
1362
{
1363
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1364
ss = (r + s)*(r + s + 1)/2 + s;
1365
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1366
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1367
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1368
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1369
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1370
}// end loop over 's'
1371
}// end loop over 'r'
1372
for (unsigned int r = 0; r < 4; r++)
1373
{
1374
for (unsigned int s = 0; s < 4 - r; s++)
1375
{
1376
rr = (r + s)*(r + s + 1)/2 + s;
1377
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1378
}// end loop over 's'
1379
}// end loop over 'r'
1380
1381
// Table(s) of coefficients.
1382
static const double coefficients0[10] = \
1383
{0.04714045, 0.00000000, 0.03333333, 0.00000000, 0.00000000, 0.10497813, 0.00000000, 0.00000000, 0.00000000, 0.09091373};
1384
1385
// Tables of derivatives of the polynomial base (transpose).
1386
static const double dmats0[10][10] = \
1387
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1388
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1389
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1390
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1391
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1392
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1393
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1394
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1395
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1396
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1397
1398
static const double dmats1[10][10] = \
1399
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1400
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1401
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1402
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1403
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1404
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1405
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1406
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1407
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1408
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1409
1410
// Compute reference derivatives.
1411
// Declare pointer to array of derivatives on FIAT element.
1412
double *derivatives = new double[num_derivatives];
1413
for (unsigned int r = 0; r < num_derivatives; r++)
1414
{
1415
derivatives[r] = 0.00000000;
1416
}// end loop over 'r'
1417
1418
// Declare derivative matrix (of polynomial basis).
1419
double dmats[10][10] = \
1420
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1421
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1422
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1423
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1424
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1425
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1426
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1427
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1428
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1429
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1430
1431
// Declare (auxiliary) derivative matrix (of polynomial basis).
1432
double dmats_old[10][10] = \
1433
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1434
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1435
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1436
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1437
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1438
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1439
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1440
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1441
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1442
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1443
1444
// Loop possible derivatives.
1445
for (unsigned int r = 0; r < num_derivatives; r++)
1446
{
1447
// Resetting dmats values to compute next derivative.
1448
for (unsigned int t = 0; t < 10; t++)
1449
{
1450
for (unsigned int u = 0; u < 10; u++)
1451
{
1452
dmats[t][u] = 0.00000000;
1453
if (t == u)
1454
{
1455
dmats[t][u] = 1.00000000;
1456
}
1457
1458
}// end loop over 'u'
1459
}// end loop over 't'
1460
1461
// Looping derivative order to generate dmats.
1462
for (unsigned int s = 0; s < n; s++)
1463
{
1464
// Updating dmats_old with new values and resetting dmats.
1465
for (unsigned int t = 0; t < 10; t++)
1466
{
1467
for (unsigned int u = 0; u < 10; u++)
1468
{
1469
dmats_old[t][u] = dmats[t][u];
1470
dmats[t][u] = 0.00000000;
1471
}// end loop over 'u'
1472
}// end loop over 't'
1473
1474
// Update dmats using an inner product.
1475
if (combinations[r][s] == 0)
1476
{
1477
for (unsigned int t = 0; t < 10; t++)
1478
{
1479
for (unsigned int u = 0; u < 10; u++)
1480
{
1481
for (unsigned int tu = 0; tu < 10; tu++)
1482
{
1483
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1484
}// end loop over 'tu'
1485
}// end loop over 'u'
1486
}// end loop over 't'
1487
}
1488
1489
if (combinations[r][s] == 1)
1490
{
1491
for (unsigned int t = 0; t < 10; t++)
1492
{
1493
for (unsigned int u = 0; u < 10; u++)
1494
{
1495
for (unsigned int tu = 0; tu < 10; tu++)
1496
{
1497
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1498
}// end loop over 'tu'
1499
}// end loop over 'u'
1500
}// end loop over 't'
1501
}
1502
1503
}// end loop over 's'
1504
for (unsigned int s = 0; s < 10; s++)
1505
{
1506
for (unsigned int t = 0; t < 10; t++)
1507
{
1508
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1509
}// end loop over 't'
1510
}// end loop over 's'
1511
}// end loop over 'r'
1512
1513
// Transform derivatives back to physical element
1514
for (unsigned int r = 0; r < num_derivatives; r++)
1515
{
1516
for (unsigned int s = 0; s < num_derivatives; s++)
1517
{
1518
values[r] += transform[r][s]*derivatives[s];
1519
}// end loop over 's'
1520
}// end loop over 'r'
1521
1522
// Delete pointer to array of derivatives on FIAT element
1523
delete [] derivatives;
1524
1525
// Delete pointer to array of combinations of derivatives and transform
1526
for (unsigned int r = 0; r < num_derivatives; r++)
1527
{
1528
delete [] combinations[r];
1529
}// end loop over 'r'
1530
delete [] combinations;
1531
for (unsigned int r = 0; r < num_derivatives; r++)
1532
{
1533
delete [] transform[r];
1534
}// end loop over 'r'
1535
delete [] transform;
1536
break;
1537
}
1538
case 3:
1539
{
1540
1541
// Array of basisvalues.
1542
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1543
1544
// Declare helper variables.
1545
unsigned int rr = 0;
1546
unsigned int ss = 0;
1547
unsigned int tt = 0;
1548
double tmp5 = 0.00000000;
1549
double tmp6 = 0.00000000;
1550
double tmp7 = 0.00000000;
1551
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1552
double tmp1 = (1.00000000 - Y)/2.00000000;
1553
double tmp2 = tmp1*tmp1;
1554
1555
// Compute basisvalues.
1556
basisvalues[0] = 1.00000000;
1557
basisvalues[1] = tmp0;
1558
for (unsigned int r = 1; r < 3; r++)
1559
{
1560
rr = (r + 1)*((r + 1) + 1)/2;
1561
ss = r*(r + 1)/2;
1562
tt = (r - 1)*((r - 1) + 1)/2;
1563
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1564
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1565
}// end loop over 'r'
1566
for (unsigned int r = 0; r < 3; r++)
1567
{
1568
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1569
ss = r*(r + 1)/2;
1570
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1571
}// end loop over 'r'
1572
for (unsigned int r = 0; r < 2; r++)
1573
{
1574
for (unsigned int s = 1; s < 3 - r; s++)
1575
{
1576
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1577
ss = (r + s)*(r + s + 1)/2 + s;
1578
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1579
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1580
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1581
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1582
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1583
}// end loop over 's'
1584
}// end loop over 'r'
1585
for (unsigned int r = 0; r < 4; r++)
1586
{
1587
for (unsigned int s = 0; s < 4 - r; s++)
1588
{
1589
rr = (r + s)*(r + s + 1)/2 + s;
1590
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1591
}// end loop over 's'
1592
}// end loop over 'r'
1593
1594
// Table(s) of coefficients.
1595
static const double coefficients0[10] = \
1596
{0.10606602, 0.25980762, -0.15000000, 0.11736912, 0.06060915, -0.07873360, 0.00000000, 0.10164464, -0.13122266, 0.09091373};
1597
1598
// Tables of derivatives of the polynomial base (transpose).
1599
static const double dmats0[10][10] = \
1600
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1601
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1602
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1603
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1604
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1605
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1606
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1607
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1608
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1609
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1610
1611
static const double dmats1[10][10] = \
1612
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1613
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1614
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1615
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1616
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1617
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1618
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1619
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1620
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1621
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1622
1623
// Compute reference derivatives.
1624
// Declare pointer to array of derivatives on FIAT element.
1625
double *derivatives = new double[num_derivatives];
1626
for (unsigned int r = 0; r < num_derivatives; r++)
1627
{
1628
derivatives[r] = 0.00000000;
1629
}// end loop over 'r'
1630
1631
// Declare derivative matrix (of polynomial basis).
1632
double dmats[10][10] = \
1633
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1634
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1635
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1636
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1637
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1638
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1639
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1640
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1641
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1642
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1643
1644
// Declare (auxiliary) derivative matrix (of polynomial basis).
1645
double dmats_old[10][10] = \
1646
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1647
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1648
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1649
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1650
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1651
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1652
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1653
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1654
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1655
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1656
1657
// Loop possible derivatives.
1658
for (unsigned int r = 0; r < num_derivatives; r++)
1659
{
1660
// Resetting dmats values to compute next derivative.
1661
for (unsigned int t = 0; t < 10; t++)
1662
{
1663
for (unsigned int u = 0; u < 10; u++)
1664
{
1665
dmats[t][u] = 0.00000000;
1666
if (t == u)
1667
{
1668
dmats[t][u] = 1.00000000;
1669
}
1670
1671
}// end loop over 'u'
1672
}// end loop over 't'
1673
1674
// Looping derivative order to generate dmats.
1675
for (unsigned int s = 0; s < n; s++)
1676
{
1677
// Updating dmats_old with new values and resetting dmats.
1678
for (unsigned int t = 0; t < 10; t++)
1679
{
1680
for (unsigned int u = 0; u < 10; u++)
1681
{
1682
dmats_old[t][u] = dmats[t][u];
1683
dmats[t][u] = 0.00000000;
1684
}// end loop over 'u'
1685
}// end loop over 't'
1686
1687
// Update dmats using an inner product.
1688
if (combinations[r][s] == 0)
1689
{
1690
for (unsigned int t = 0; t < 10; t++)
1691
{
1692
for (unsigned int u = 0; u < 10; u++)
1693
{
1694
for (unsigned int tu = 0; tu < 10; tu++)
1695
{
1696
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1697
}// end loop over 'tu'
1698
}// end loop over 'u'
1699
}// end loop over 't'
1700
}
1701
1702
if (combinations[r][s] == 1)
1703
{
1704
for (unsigned int t = 0; t < 10; t++)
1705
{
1706
for (unsigned int u = 0; u < 10; u++)
1707
{
1708
for (unsigned int tu = 0; tu < 10; tu++)
1709
{
1710
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1711
}// end loop over 'tu'
1712
}// end loop over 'u'
1713
}// end loop over 't'
1714
}
1715
1716
}// end loop over 's'
1717
for (unsigned int s = 0; s < 10; s++)
1718
{
1719
for (unsigned int t = 0; t < 10; t++)
1720
{
1721
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1722
}// end loop over 't'
1723
}// end loop over 's'
1724
}// end loop over 'r'
1725
1726
// Transform derivatives back to physical element
1727
for (unsigned int r = 0; r < num_derivatives; r++)
1728
{
1729
for (unsigned int s = 0; s < num_derivatives; s++)
1730
{
1731
values[r] += transform[r][s]*derivatives[s];
1732
}// end loop over 's'
1733
}// end loop over 'r'
1734
1735
// Delete pointer to array of derivatives on FIAT element
1736
delete [] derivatives;
1737
1738
// Delete pointer to array of combinations of derivatives and transform
1739
for (unsigned int r = 0; r < num_derivatives; r++)
1740
{
1741
delete [] combinations[r];
1742
}// end loop over 'r'
1743
delete [] combinations;
1744
for (unsigned int r = 0; r < num_derivatives; r++)
1745
{
1746
delete [] transform[r];
1747
}// end loop over 'r'
1748
delete [] transform;
1749
break;
1750
}
1751
case 4:
1752
{
1753
1754
// Array of basisvalues.
1755
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1756
1757
// Declare helper variables.
1758
unsigned int rr = 0;
1759
unsigned int ss = 0;
1760
unsigned int tt = 0;
1761
double tmp5 = 0.00000000;
1762
double tmp6 = 0.00000000;
1763
double tmp7 = 0.00000000;
1764
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1765
double tmp1 = (1.00000000 - Y)/2.00000000;
1766
double tmp2 = tmp1*tmp1;
1767
1768
// Compute basisvalues.
1769
basisvalues[0] = 1.00000000;
1770
basisvalues[1] = tmp0;
1771
for (unsigned int r = 1; r < 3; r++)
1772
{
1773
rr = (r + 1)*((r + 1) + 1)/2;
1774
ss = r*(r + 1)/2;
1775
tt = (r - 1)*((r - 1) + 1)/2;
1776
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1777
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1778
}// end loop over 'r'
1779
for (unsigned int r = 0; r < 3; r++)
1780
{
1781
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1782
ss = r*(r + 1)/2;
1783
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1784
}// end loop over 'r'
1785
for (unsigned int r = 0; r < 2; r++)
1786
{
1787
for (unsigned int s = 1; s < 3 - r; s++)
1788
{
1789
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1790
ss = (r + s)*(r + s + 1)/2 + s;
1791
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1792
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1793
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1794
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1795
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1796
}// end loop over 's'
1797
}// end loop over 'r'
1798
for (unsigned int r = 0; r < 4; r++)
1799
{
1800
for (unsigned int s = 0; s < 4 - r; s++)
1801
{
1802
rr = (r + s)*(r + s + 1)/2 + s;
1803
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1804
}// end loop over 's'
1805
}// end loop over 'r'
1806
1807
// Table(s) of coefficients.
1808
static const double coefficients0[10] = \
1809
{0.10606602, 0.00000000, 0.30000000, 0.00000000, 0.15152288, 0.02624453, 0.00000000, 0.00000000, 0.13122266, -0.13637059};
1810
1811
// Tables of derivatives of the polynomial base (transpose).
1812
static const double dmats0[10][10] = \
1813
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1814
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1815
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1816
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1817
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1818
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1819
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1820
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1821
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1822
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1823
1824
static const double dmats1[10][10] = \
1825
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1826
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1827
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1828
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1829
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1830
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1831
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1832
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1833
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1834
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1835
1836
// Compute reference derivatives.
1837
// Declare pointer to array of derivatives on FIAT element.
1838
double *derivatives = new double[num_derivatives];
1839
for (unsigned int r = 0; r < num_derivatives; r++)
1840
{
1841
derivatives[r] = 0.00000000;
1842
}// end loop over 'r'
1843
1844
// Declare derivative matrix (of polynomial basis).
1845
double dmats[10][10] = \
1846
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1847
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1848
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1849
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1850
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1851
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1852
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1853
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1854
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1855
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1856
1857
// Declare (auxiliary) derivative matrix (of polynomial basis).
1858
double dmats_old[10][10] = \
1859
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1860
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1861
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1862
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1863
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1864
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1865
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1866
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1867
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1868
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1869
1870
// Loop possible derivatives.
1871
for (unsigned int r = 0; r < num_derivatives; r++)
1872
{
1873
// Resetting dmats values to compute next derivative.
1874
for (unsigned int t = 0; t < 10; t++)
1875
{
1876
for (unsigned int u = 0; u < 10; u++)
1877
{
1878
dmats[t][u] = 0.00000000;
1879
if (t == u)
1880
{
1881
dmats[t][u] = 1.00000000;
1882
}
1883
1884
}// end loop over 'u'
1885
}// end loop over 't'
1886
1887
// Looping derivative order to generate dmats.
1888
for (unsigned int s = 0; s < n; s++)
1889
{
1890
// Updating dmats_old with new values and resetting dmats.
1891
for (unsigned int t = 0; t < 10; t++)
1892
{
1893
for (unsigned int u = 0; u < 10; u++)
1894
{
1895
dmats_old[t][u] = dmats[t][u];
1896
dmats[t][u] = 0.00000000;
1897
}// end loop over 'u'
1898
}// end loop over 't'
1899
1900
// Update dmats using an inner product.
1901
if (combinations[r][s] == 0)
1902
{
1903
for (unsigned int t = 0; t < 10; t++)
1904
{
1905
for (unsigned int u = 0; u < 10; u++)
1906
{
1907
for (unsigned int tu = 0; tu < 10; tu++)
1908
{
1909
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1910
}// end loop over 'tu'
1911
}// end loop over 'u'
1912
}// end loop over 't'
1913
}
1914
1915
if (combinations[r][s] == 1)
1916
{
1917
for (unsigned int t = 0; t < 10; t++)
1918
{
1919
for (unsigned int u = 0; u < 10; u++)
1920
{
1921
for (unsigned int tu = 0; tu < 10; tu++)
1922
{
1923
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1924
}// end loop over 'tu'
1925
}// end loop over 'u'
1926
}// end loop over 't'
1927
}
1928
1929
}// end loop over 's'
1930
for (unsigned int s = 0; s < 10; s++)
1931
{
1932
for (unsigned int t = 0; t < 10; t++)
1933
{
1934
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1935
}// end loop over 't'
1936
}// end loop over 's'
1937
}// end loop over 'r'
1938
1939
// Transform derivatives back to physical element
1940
for (unsigned int r = 0; r < num_derivatives; r++)
1941
{
1942
for (unsigned int s = 0; s < num_derivatives; s++)
1943
{
1944
values[r] += transform[r][s]*derivatives[s];
1945
}// end loop over 's'
1946
}// end loop over 'r'
1947
1948
// Delete pointer to array of derivatives on FIAT element
1949
delete [] derivatives;
1950
1951
// Delete pointer to array of combinations of derivatives and transform
1952
for (unsigned int r = 0; r < num_derivatives; r++)
1953
{
1954
delete [] combinations[r];
1955
}// end loop over 'r'
1956
delete [] combinations;
1957
for (unsigned int r = 0; r < num_derivatives; r++)
1958
{
1959
delete [] transform[r];
1960
}// end loop over 'r'
1961
delete [] transform;
1962
break;
1963
}
1964
case 5:
1965
{
1966
1967
// Array of basisvalues.
1968
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1969
1970
// Declare helper variables.
1971
unsigned int rr = 0;
1972
unsigned int ss = 0;
1973
unsigned int tt = 0;
1974
double tmp5 = 0.00000000;
1975
double tmp6 = 0.00000000;
1976
double tmp7 = 0.00000000;
1977
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1978
double tmp1 = (1.00000000 - Y)/2.00000000;
1979
double tmp2 = tmp1*tmp1;
1980
1981
// Compute basisvalues.
1982
basisvalues[0] = 1.00000000;
1983
basisvalues[1] = tmp0;
1984
for (unsigned int r = 1; r < 3; r++)
1985
{
1986
rr = (r + 1)*((r + 1) + 1)/2;
1987
ss = r*(r + 1)/2;
1988
tt = (r - 1)*((r - 1) + 1)/2;
1989
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1990
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1991
}// end loop over 'r'
1992
for (unsigned int r = 0; r < 3; r++)
1993
{
1994
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1995
ss = r*(r + 1)/2;
1996
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1997
}// end loop over 'r'
1998
for (unsigned int r = 0; r < 2; r++)
1999
{
2000
for (unsigned int s = 1; s < 3 - r; s++)
2001
{
2002
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
2003
ss = (r + s)*(r + s + 1)/2 + s;
2004
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
2005
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2006
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2007
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2008
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
2009
}// end loop over 's'
2010
}// end loop over 'r'
2011
for (unsigned int r = 0; r < 4; r++)
2012
{
2013
for (unsigned int s = 0; s < 4 - r; s++)
2014
{
2015
rr = (r + s)*(r + s + 1)/2 + s;
2016
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
2017
}// end loop over 's'
2018
}// end loop over 'r'
2019
2020
// Table(s) of coefficients.
2021
static const double coefficients0[10] = \
2022
{0.10606602, -0.25980762, -0.15000000, 0.11736912, -0.06060915, -0.07873360, 0.00000000, 0.10164464, 0.13122266, 0.09091373};
2023
2024
// Tables of derivatives of the polynomial base (transpose).
2025
static const double dmats0[10][10] = \
2026
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2027
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2028
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2029
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2030
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2031
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2032
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2033
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2034
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2035
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2036
2037
static const double dmats1[10][10] = \
2038
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2039
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2040
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2041
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2042
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2043
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2044
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2045
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2046
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2047
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2048
2049
// Compute reference derivatives.
2050
// Declare pointer to array of derivatives on FIAT element.
2051
double *derivatives = new double[num_derivatives];
2052
for (unsigned int r = 0; r < num_derivatives; r++)
2053
{
2054
derivatives[r] = 0.00000000;
2055
}// end loop over 'r'
2056
2057
// Declare derivative matrix (of polynomial basis).
2058
double dmats[10][10] = \
2059
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2060
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2061
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2062
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2063
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2064
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2065
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2066
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2067
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2068
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2069
2070
// Declare (auxiliary) derivative matrix (of polynomial basis).
2071
double dmats_old[10][10] = \
2072
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2073
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2074
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2075
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2076
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2077
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2078
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2079
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2080
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2081
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2082
2083
// Loop possible derivatives.
2084
for (unsigned int r = 0; r < num_derivatives; r++)
2085
{
2086
// Resetting dmats values to compute next derivative.
2087
for (unsigned int t = 0; t < 10; t++)
2088
{
2089
for (unsigned int u = 0; u < 10; u++)
2090
{
2091
dmats[t][u] = 0.00000000;
2092
if (t == u)
2093
{
2094
dmats[t][u] = 1.00000000;
2095
}
2096
2097
}// end loop over 'u'
2098
}// end loop over 't'
2099
2100
// Looping derivative order to generate dmats.
2101
for (unsigned int s = 0; s < n; s++)
2102
{
2103
// Updating dmats_old with new values and resetting dmats.
2104
for (unsigned int t = 0; t < 10; t++)
2105
{
2106
for (unsigned int u = 0; u < 10; u++)
2107
{
2108
dmats_old[t][u] = dmats[t][u];
2109
dmats[t][u] = 0.00000000;
2110
}// end loop over 'u'
2111
}// end loop over 't'
2112
2113
// Update dmats using an inner product.
2114
if (combinations[r][s] == 0)
2115
{
2116
for (unsigned int t = 0; t < 10; t++)
2117
{
2118
for (unsigned int u = 0; u < 10; u++)
2119
{
2120
for (unsigned int tu = 0; tu < 10; tu++)
2121
{
2122
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2123
}// end loop over 'tu'
2124
}// end loop over 'u'
2125
}// end loop over 't'
2126
}
2127
2128
if (combinations[r][s] == 1)
2129
{
2130
for (unsigned int t = 0; t < 10; t++)
2131
{
2132
for (unsigned int u = 0; u < 10; u++)
2133
{
2134
for (unsigned int tu = 0; tu < 10; tu++)
2135
{
2136
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2137
}// end loop over 'tu'
2138
}// end loop over 'u'
2139
}// end loop over 't'
2140
}
2141
2142
}// end loop over 's'
2143
for (unsigned int s = 0; s < 10; s++)
2144
{
2145
for (unsigned int t = 0; t < 10; t++)
2146
{
2147
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2148
}// end loop over 't'
2149
}// end loop over 's'
2150
}// end loop over 'r'
2151
2152
// Transform derivatives back to physical element
2153
for (unsigned int r = 0; r < num_derivatives; r++)
2154
{
2155
for (unsigned int s = 0; s < num_derivatives; s++)
2156
{
2157
values[r] += transform[r][s]*derivatives[s];
2158
}// end loop over 's'
2159
}// end loop over 'r'
2160
2161
// Delete pointer to array of derivatives on FIAT element
2162
delete [] derivatives;
2163
2164
// Delete pointer to array of combinations of derivatives and transform
2165
for (unsigned int r = 0; r < num_derivatives; r++)
2166
{
2167
delete [] combinations[r];
2168
}// end loop over 'r'
2169
delete [] combinations;
2170
for (unsigned int r = 0; r < num_derivatives; r++)
2171
{
2172
delete [] transform[r];
2173
}// end loop over 'r'
2174
delete [] transform;
2175
break;
2176
}
2177
case 6:
2178
{
2179
2180
// Array of basisvalues.
2181
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
2182
2183
// Declare helper variables.
2184
unsigned int rr = 0;
2185
unsigned int ss = 0;
2186
unsigned int tt = 0;
2187
double tmp5 = 0.00000000;
2188
double tmp6 = 0.00000000;
2189
double tmp7 = 0.00000000;
2190
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
2191
double tmp1 = (1.00000000 - Y)/2.00000000;
2192
double tmp2 = tmp1*tmp1;
2193
2194
// Compute basisvalues.
2195
basisvalues[0] = 1.00000000;
2196
basisvalues[1] = tmp0;
2197
for (unsigned int r = 1; r < 3; r++)
2198
{
2199
rr = (r + 1)*((r + 1) + 1)/2;
2200
ss = r*(r + 1)/2;
2201
tt = (r - 1)*((r - 1) + 1)/2;
2202
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
2203
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
2204
}// end loop over 'r'
2205
for (unsigned int r = 0; r < 3; r++)
2206
{
2207
rr = (r + 1)*(r + 1 + 1)/2 + 1;
2208
ss = r*(r + 1)/2;
2209
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
2210
}// end loop over 'r'
2211
for (unsigned int r = 0; r < 2; r++)
2212
{
2213
for (unsigned int s = 1; s < 3 - r; s++)
2214
{
2215
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
2216
ss = (r + s)*(r + s + 1)/2 + s;
2217
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
2218
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2219
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2220
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2221
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
2222
}// end loop over 's'
2223
}// end loop over 'r'
2224
for (unsigned int r = 0; r < 4; r++)
2225
{
2226
for (unsigned int s = 0; s < 4 - r; s++)
2227
{
2228
rr = (r + s)*(r + s + 1)/2 + s;
2229
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
2230
}// end loop over 's'
2231
}// end loop over 'r'
2232
2233
// Table(s) of coefficients.
2234
static const double coefficients0[10] = \
2235
{0.10606602, 0.00000000, 0.30000000, 0.00000000, -0.15152288, 0.02624453, 0.00000000, 0.00000000, -0.13122266, -0.13637059};
2236
2237
// Tables of derivatives of the polynomial base (transpose).
2238
static const double dmats0[10][10] = \
2239
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2240
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2241
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2242
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2243
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2244
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2245
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2246
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2247
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2248
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2249
2250
static const double dmats1[10][10] = \
2251
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2252
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2253
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2254
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2255
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2256
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2257
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2258
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2259
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2260
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2261
2262
// Compute reference derivatives.
2263
// Declare pointer to array of derivatives on FIAT element.
2264
double *derivatives = new double[num_derivatives];
2265
for (unsigned int r = 0; r < num_derivatives; r++)
2266
{
2267
derivatives[r] = 0.00000000;
2268
}// end loop over 'r'
2269
2270
// Declare derivative matrix (of polynomial basis).
2271
double dmats[10][10] = \
2272
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2273
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2274
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2275
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2276
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2277
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2278
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2279
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2280
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2281
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2282
2283
// Declare (auxiliary) derivative matrix (of polynomial basis).
2284
double dmats_old[10][10] = \
2285
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2286
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2287
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2288
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2289
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2290
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2291
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2292
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2293
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2294
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2295
2296
// Loop possible derivatives.
2297
for (unsigned int r = 0; r < num_derivatives; r++)
2298
{
2299
// Resetting dmats values to compute next derivative.
2300
for (unsigned int t = 0; t < 10; t++)
2301
{
2302
for (unsigned int u = 0; u < 10; u++)
2303
{
2304
dmats[t][u] = 0.00000000;
2305
if (t == u)
2306
{
2307
dmats[t][u] = 1.00000000;
2308
}
2309
2310
}// end loop over 'u'
2311
}// end loop over 't'
2312
2313
// Looping derivative order to generate dmats.
2314
for (unsigned int s = 0; s < n; s++)
2315
{
2316
// Updating dmats_old with new values and resetting dmats.
2317
for (unsigned int t = 0; t < 10; t++)
2318
{
2319
for (unsigned int u = 0; u < 10; u++)
2320
{
2321
dmats_old[t][u] = dmats[t][u];
2322
dmats[t][u] = 0.00000000;
2323
}// end loop over 'u'
2324
}// end loop over 't'
2325
2326
// Update dmats using an inner product.
2327
if (combinations[r][s] == 0)
2328
{
2329
for (unsigned int t = 0; t < 10; t++)
2330
{
2331
for (unsigned int u = 0; u < 10; u++)
2332
{
2333
for (unsigned int tu = 0; tu < 10; tu++)
2334
{
2335
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2336
}// end loop over 'tu'
2337
}// end loop over 'u'
2338
}// end loop over 't'
2339
}
2340
2341
if (combinations[r][s] == 1)
2342
{
2343
for (unsigned int t = 0; t < 10; t++)
2344
{
2345
for (unsigned int u = 0; u < 10; u++)
2346
{
2347
for (unsigned int tu = 0; tu < 10; tu++)
2348
{
2349
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2350
}// end loop over 'tu'
2351
}// end loop over 'u'
2352
}// end loop over 't'
2353
}
2354
2355
}// end loop over 's'
2356
for (unsigned int s = 0; s < 10; s++)
2357
{
2358
for (unsigned int t = 0; t < 10; t++)
2359
{
2360
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2361
}// end loop over 't'
2362
}// end loop over 's'
2363
}// end loop over 'r'
2364
2365
// Transform derivatives back to physical element
2366
for (unsigned int r = 0; r < num_derivatives; r++)
2367
{
2368
for (unsigned int s = 0; s < num_derivatives; s++)
2369
{
2370
values[r] += transform[r][s]*derivatives[s];
2371
}// end loop over 's'
2372
}// end loop over 'r'
2373
2374
// Delete pointer to array of derivatives on FIAT element
2375
delete [] derivatives;
2376
2377
// Delete pointer to array of combinations of derivatives and transform
2378
for (unsigned int r = 0; r < num_derivatives; r++)
2379
{
2380
delete [] combinations[r];
2381
}// end loop over 'r'
2382
delete [] combinations;
2383
for (unsigned int r = 0; r < num_derivatives; r++)
2384
{
2385
delete [] transform[r];
2386
}// end loop over 'r'
2387
delete [] transform;
2388
break;
2389
}
2390
case 7:
2391
{
2392
2393
// Array of basisvalues.
2394
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
2395
2396
// Declare helper variables.
2397
unsigned int rr = 0;
2398
unsigned int ss = 0;
2399
unsigned int tt = 0;
2400
double tmp5 = 0.00000000;
2401
double tmp6 = 0.00000000;
2402
double tmp7 = 0.00000000;
2403
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
2404
double tmp1 = (1.00000000 - Y)/2.00000000;
2405
double tmp2 = tmp1*tmp1;
2406
2407
// Compute basisvalues.
2408
basisvalues[0] = 1.00000000;
2409
basisvalues[1] = tmp0;
2410
for (unsigned int r = 1; r < 3; r++)
2411
{
2412
rr = (r + 1)*((r + 1) + 1)/2;
2413
ss = r*(r + 1)/2;
2414
tt = (r - 1)*((r - 1) + 1)/2;
2415
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
2416
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
2417
}// end loop over 'r'
2418
for (unsigned int r = 0; r < 3; r++)
2419
{
2420
rr = (r + 1)*(r + 1 + 1)/2 + 1;
2421
ss = r*(r + 1)/2;
2422
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
2423
}// end loop over 'r'
2424
for (unsigned int r = 0; r < 2; r++)
2425
{
2426
for (unsigned int s = 1; s < 3 - r; s++)
2427
{
2428
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
2429
ss = (r + s)*(r + s + 1)/2 + s;
2430
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
2431
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2432
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2433
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2434
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
2435
}// end loop over 's'
2436
}// end loop over 'r'
2437
for (unsigned int r = 0; r < 4; r++)
2438
{
2439
for (unsigned int s = 0; s < 4 - r; s++)
2440
{
2441
rr = (r + s)*(r + s + 1)/2 + s;
2442
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
2443
}// end loop over 's'
2444
}// end loop over 'r'
2445
2446
// Table(s) of coefficients.
2447
static const double coefficients0[10] = \
2448
{0.10606602, -0.25980762, -0.15000000, -0.07824608, 0.09091373, 0.09622995, 0.18040134, 0.05082232, -0.01312227, -0.02272843};
2449
2450
// Tables of derivatives of the polynomial base (transpose).
2451
static const double dmats0[10][10] = \
2452
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2453
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2454
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2455
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2456
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2457
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2458
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2459
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2460
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2461
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2462
2463
static const double dmats1[10][10] = \
2464
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2465
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2466
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2467
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2468
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2469
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2470
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2471
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2472
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2473
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2474
2475
// Compute reference derivatives.
2476
// Declare pointer to array of derivatives on FIAT element.
2477
double *derivatives = new double[num_derivatives];
2478
for (unsigned int r = 0; r < num_derivatives; r++)
2479
{
2480
derivatives[r] = 0.00000000;
2481
}// end loop over 'r'
2482
2483
// Declare derivative matrix (of polynomial basis).
2484
double dmats[10][10] = \
2485
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2486
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2487
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2488
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2489
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2490
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2491
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2492
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2493
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2494
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2495
2496
// Declare (auxiliary) derivative matrix (of polynomial basis).
2497
double dmats_old[10][10] = \
2498
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2499
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2500
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2501
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2502
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2503
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2504
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2505
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2506
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2507
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2508
2509
// Loop possible derivatives.
2510
for (unsigned int r = 0; r < num_derivatives; r++)
2511
{
2512
// Resetting dmats values to compute next derivative.
2513
for (unsigned int t = 0; t < 10; t++)
2514
{
2515
for (unsigned int u = 0; u < 10; u++)
2516
{
2517
dmats[t][u] = 0.00000000;
2518
if (t == u)
2519
{
2520
dmats[t][u] = 1.00000000;
2521
}
2522
2523
}// end loop over 'u'
2524
}// end loop over 't'
2525
2526
// Looping derivative order to generate dmats.
2527
for (unsigned int s = 0; s < n; s++)
2528
{
2529
// Updating dmats_old with new values and resetting dmats.
2530
for (unsigned int t = 0; t < 10; t++)
2531
{
2532
for (unsigned int u = 0; u < 10; u++)
2533
{
2534
dmats_old[t][u] = dmats[t][u];
2535
dmats[t][u] = 0.00000000;
2536
}// end loop over 'u'
2537
}// end loop over 't'
2538
2539
// Update dmats using an inner product.
2540
if (combinations[r][s] == 0)
2541
{
2542
for (unsigned int t = 0; t < 10; t++)
2543
{
2544
for (unsigned int u = 0; u < 10; u++)
2545
{
2546
for (unsigned int tu = 0; tu < 10; tu++)
2547
{
2548
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2549
}// end loop over 'tu'
2550
}// end loop over 'u'
2551
}// end loop over 't'
2552
}
2553
2554
if (combinations[r][s] == 1)
2555
{
2556
for (unsigned int t = 0; t < 10; t++)
2557
{
2558
for (unsigned int u = 0; u < 10; u++)
2559
{
2560
for (unsigned int tu = 0; tu < 10; tu++)
2561
{
2562
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2563
}// end loop over 'tu'
2564
}// end loop over 'u'
2565
}// end loop over 't'
2566
}
2567
2568
}// end loop over 's'
2569
for (unsigned int s = 0; s < 10; s++)
2570
{
2571
for (unsigned int t = 0; t < 10; t++)
2572
{
2573
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2574
}// end loop over 't'
2575
}// end loop over 's'
2576
}// end loop over 'r'
2577
2578
// Transform derivatives back to physical element
2579
for (unsigned int r = 0; r < num_derivatives; r++)
2580
{
2581
for (unsigned int s = 0; s < num_derivatives; s++)
2582
{
2583
values[r] += transform[r][s]*derivatives[s];
2584
}// end loop over 's'
2585
}// end loop over 'r'
2586
2587
// Delete pointer to array of derivatives on FIAT element
2588
delete [] derivatives;
2589
2590
// Delete pointer to array of combinations of derivatives and transform
2591
for (unsigned int r = 0; r < num_derivatives; r++)
2592
{
2593
delete [] combinations[r];
2594
}// end loop over 'r'
2595
delete [] combinations;
2596
for (unsigned int r = 0; r < num_derivatives; r++)
2597
{
2598
delete [] transform[r];
2599
}// end loop over 'r'
2600
delete [] transform;
2601
break;
2602
}
2603
case 8:
2604
{
2605
2606
// Array of basisvalues.
2607
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
2608
2609
// Declare helper variables.
2610
unsigned int rr = 0;
2611
unsigned int ss = 0;
2612
unsigned int tt = 0;
2613
double tmp5 = 0.00000000;
2614
double tmp6 = 0.00000000;
2615
double tmp7 = 0.00000000;
2616
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
2617
double tmp1 = (1.00000000 - Y)/2.00000000;
2618
double tmp2 = tmp1*tmp1;
2619
2620
// Compute basisvalues.
2621
basisvalues[0] = 1.00000000;
2622
basisvalues[1] = tmp0;
2623
for (unsigned int r = 1; r < 3; r++)
2624
{
2625
rr = (r + 1)*((r + 1) + 1)/2;
2626
ss = r*(r + 1)/2;
2627
tt = (r - 1)*((r - 1) + 1)/2;
2628
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
2629
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
2630
}// end loop over 'r'
2631
for (unsigned int r = 0; r < 3; r++)
2632
{
2633
rr = (r + 1)*(r + 1 + 1)/2 + 1;
2634
ss = r*(r + 1)/2;
2635
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
2636
}// end loop over 'r'
2637
for (unsigned int r = 0; r < 2; r++)
2638
{
2639
for (unsigned int s = 1; s < 3 - r; s++)
2640
{
2641
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
2642
ss = (r + s)*(r + s + 1)/2 + s;
2643
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
2644
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2645
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2646
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2647
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
2648
}// end loop over 's'
2649
}// end loop over 'r'
2650
for (unsigned int r = 0; r < 4; r++)
2651
{
2652
for (unsigned int s = 0; s < 4 - r; s++)
2653
{
2654
rr = (r + s)*(r + s + 1)/2 + s;
2655
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
2656
}// end loop over 's'
2657
}// end loop over 'r'
2658
2659
// Table(s) of coefficients.
2660
static const double coefficients0[10] = \
2661
{0.10606602, 0.25980762, -0.15000000, -0.07824608, -0.09091373, 0.09622995, -0.18040134, 0.05082232, 0.01312227, -0.02272843};
2662
2663
// Tables of derivatives of the polynomial base (transpose).
2664
static const double dmats0[10][10] = \
2665
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2666
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2667
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2668
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2669
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2670
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2671
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2672
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2673
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2674
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2675
2676
static const double dmats1[10][10] = \
2677
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2678
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2679
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2680
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2681
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2682
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2683
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2684
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2685
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2686
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2687
2688
// Compute reference derivatives.
2689
// Declare pointer to array of derivatives on FIAT element.
2690
double *derivatives = new double[num_derivatives];
2691
for (unsigned int r = 0; r < num_derivatives; r++)
2692
{
2693
derivatives[r] = 0.00000000;
2694
}// end loop over 'r'
2695
2696
// Declare derivative matrix (of polynomial basis).
2697
double dmats[10][10] = \
2698
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2699
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2700
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2701
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2702
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2703
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2704
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2705
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2706
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2707
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2708
2709
// Declare (auxiliary) derivative matrix (of polynomial basis).
2710
double dmats_old[10][10] = \
2711
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2712
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2713
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2714
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2715
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2716
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2717
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2718
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2719
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2720
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2721
2722
// Loop possible derivatives.
2723
for (unsigned int r = 0; r < num_derivatives; r++)
2724
{
2725
// Resetting dmats values to compute next derivative.
2726
for (unsigned int t = 0; t < 10; t++)
2727
{
2728
for (unsigned int u = 0; u < 10; u++)
2729
{
2730
dmats[t][u] = 0.00000000;
2731
if (t == u)
2732
{
2733
dmats[t][u] = 1.00000000;
2734
}
2735
2736
}// end loop over 'u'
2737
}// end loop over 't'
2738
2739
// Looping derivative order to generate dmats.
2740
for (unsigned int s = 0; s < n; s++)
2741
{
2742
// Updating dmats_old with new values and resetting dmats.
2743
for (unsigned int t = 0; t < 10; t++)
2744
{
2745
for (unsigned int u = 0; u < 10; u++)
2746
{
2747
dmats_old[t][u] = dmats[t][u];
2748
dmats[t][u] = 0.00000000;
2749
}// end loop over 'u'
2750
}// end loop over 't'
2751
2752
// Update dmats using an inner product.
2753
if (combinations[r][s] == 0)
2754
{
2755
for (unsigned int t = 0; t < 10; t++)
2756
{
2757
for (unsigned int u = 0; u < 10; u++)
2758
{
2759
for (unsigned int tu = 0; tu < 10; tu++)
2760
{
2761
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2762
}// end loop over 'tu'
2763
}// end loop over 'u'
2764
}// end loop over 't'
2765
}
2766
2767
if (combinations[r][s] == 1)
2768
{
2769
for (unsigned int t = 0; t < 10; t++)
2770
{
2771
for (unsigned int u = 0; u < 10; u++)
2772
{
2773
for (unsigned int tu = 0; tu < 10; tu++)
2774
{
2775
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2776
}// end loop over 'tu'
2777
}// end loop over 'u'
2778
}// end loop over 't'
2779
}
2780
2781
}// end loop over 's'
2782
for (unsigned int s = 0; s < 10; s++)
2783
{
2784
for (unsigned int t = 0; t < 10; t++)
2785
{
2786
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2787
}// end loop over 't'
2788
}// end loop over 's'
2789
}// end loop over 'r'
2790
2791
// Transform derivatives back to physical element
2792
for (unsigned int r = 0; r < num_derivatives; r++)
2793
{
2794
for (unsigned int s = 0; s < num_derivatives; s++)
2795
{
2796
values[r] += transform[r][s]*derivatives[s];
2797
}// end loop over 's'
2798
}// end loop over 'r'
2799
2800
// Delete pointer to array of derivatives on FIAT element
2801
delete [] derivatives;
2802
2803
// Delete pointer to array of combinations of derivatives and transform
2804
for (unsigned int r = 0; r < num_derivatives; r++)
2805
{
2806
delete [] combinations[r];
2807
}// end loop over 'r'
2808
delete [] combinations;
2809
for (unsigned int r = 0; r < num_derivatives; r++)
2810
{
2811
delete [] transform[r];
2812
}// end loop over 'r'
2813
delete [] transform;
2814
break;
2815
}
2816
case 9:
2817
{
2818
2819
// Array of basisvalues.
2820
double basisvalues[10] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
2821
2822
// Declare helper variables.
2823
unsigned int rr = 0;
2824
unsigned int ss = 0;
2825
unsigned int tt = 0;
2826
double tmp5 = 0.00000000;
2827
double tmp6 = 0.00000000;
2828
double tmp7 = 0.00000000;
2829
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
2830
double tmp1 = (1.00000000 - Y)/2.00000000;
2831
double tmp2 = tmp1*tmp1;
2832
2833
// Compute basisvalues.
2834
basisvalues[0] = 1.00000000;
2835
basisvalues[1] = tmp0;
2836
for (unsigned int r = 1; r < 3; r++)
2837
{
2838
rr = (r + 1)*((r + 1) + 1)/2;
2839
ss = r*(r + 1)/2;
2840
tt = (r - 1)*((r - 1) + 1)/2;
2841
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
2842
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
2843
}// end loop over 'r'
2844
for (unsigned int r = 0; r < 3; r++)
2845
{
2846
rr = (r + 1)*(r + 1 + 1)/2 + 1;
2847
ss = r*(r + 1)/2;
2848
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
2849
}// end loop over 'r'
2850
for (unsigned int r = 0; r < 2; r++)
2851
{
2852
for (unsigned int s = 1; s < 3 - r; s++)
2853
{
2854
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
2855
ss = (r + s)*(r + s + 1)/2 + s;
2856
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
2857
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2858
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2859
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2860
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
2861
}// end loop over 's'
2862
}// end loop over 'r'
2863
for (unsigned int r = 0; r < 4; r++)
2864
{
2865
for (unsigned int s = 0; s < 4 - r; s++)
2866
{
2867
rr = (r + s)*(r + s + 1)/2 + s;
2868
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
2869
}// end loop over 's'
2870
}// end loop over 'r'
2871
2872
// Table(s) of coefficients.
2873
static const double coefficients0[10] = \
2874
{0.63639610, 0.00000000, 0.00000000, -0.23473824, 0.00000000, -0.26244533, 0.00000000, -0.20328928, 0.00000000, 0.09091373};
2875
2876
// Tables of derivatives of the polynomial base (transpose).
2877
static const double dmats0[10][10] = \
2878
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2879
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2880
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2881
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2882
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2883
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2884
{5.29150262, 0.00000000, -2.99332591, 13.66260102, 0.00000000, 0.61101009, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2885
{0.00000000, 4.38178046, 0.00000000, 0.00000000, 12.52198067, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2886
{3.46410162, 0.00000000, 7.83836718, 0.00000000, 0.00000000, 8.40000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2887
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2888
2889
static const double dmats1[10][10] = \
2890
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2891
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2892
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2893
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2894
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2895
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2896
{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.30550505, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2897
{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2898
{1.73205081, -5.09116882, 3.91918359, 0.00000000, 9.69948452, 4.20000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2899
{5.00000000, 0.00000000, -2.82842712, 0.00000000, 0.00000000, 12.12435565, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2900
2901
// Compute reference derivatives.
2902
// Declare pointer to array of derivatives on FIAT element.
2903
double *derivatives = new double[num_derivatives];
2904
for (unsigned int r = 0; r < num_derivatives; r++)
2905
{
2906
derivatives[r] = 0.00000000;
2907
}// end loop over 'r'
2908
2909
// Declare derivative matrix (of polynomial basis).
2910
double dmats[10][10] = \
2911
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2912
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2913
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2914
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2915
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2916
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2917
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2918
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2919
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2920
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2921
2922
// Declare (auxiliary) derivative matrix (of polynomial basis).
2923
double dmats_old[10][10] = \
2924
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2925
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2926
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2927
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2928
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2929
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2930
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2931
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2932
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2933
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2934
2935
// Loop possible derivatives.
2936
for (unsigned int r = 0; r < num_derivatives; r++)
2937
{
2938
// Resetting dmats values to compute next derivative.
2939
for (unsigned int t = 0; t < 10; t++)
2940
{
2941
for (unsigned int u = 0; u < 10; u++)
2942
{
2943
dmats[t][u] = 0.00000000;
2944
if (t == u)
2945
{
2946
dmats[t][u] = 1.00000000;
2947
}
2948
2949
}// end loop over 'u'
2950
}// end loop over 't'
2951
2952
// Looping derivative order to generate dmats.
2953
for (unsigned int s = 0; s < n; s++)
2954
{
2955
// Updating dmats_old with new values and resetting dmats.
2956
for (unsigned int t = 0; t < 10; t++)
2957
{
2958
for (unsigned int u = 0; u < 10; u++)
2959
{
2960
dmats_old[t][u] = dmats[t][u];
2961
dmats[t][u] = 0.00000000;
2962
}// end loop over 'u'
2963
}// end loop over 't'
2964
2965
// Update dmats using an inner product.
2966
if (combinations[r][s] == 0)
2967
{
2968
for (unsigned int t = 0; t < 10; t++)
2969
{
2970
for (unsigned int u = 0; u < 10; u++)
2971
{
2972
for (unsigned int tu = 0; tu < 10; tu++)
2973
{
2974
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2975
}// end loop over 'tu'
2976
}// end loop over 'u'
2977
}// end loop over 't'
2978
}
2979
2980
if (combinations[r][s] == 1)
2981
{
2982
for (unsigned int t = 0; t < 10; t++)
2983
{
2984
for (unsigned int u = 0; u < 10; u++)
2985
{
2986
for (unsigned int tu = 0; tu < 10; tu++)
2987
{
2988
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2989
}// end loop over 'tu'
2990
}// end loop over 'u'
2991
}// end loop over 't'
2992
}
2993
2994
}// end loop over 's'
2995
for (unsigned int s = 0; s < 10; s++)
2996
{
2997
for (unsigned int t = 0; t < 10; t++)
2998
{
2999
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3000
}// end loop over 't'
3001
}// end loop over 's'
3002
}// end loop over 'r'
3003
3004
// Transform derivatives back to physical element
3005
for (unsigned int r = 0; r < num_derivatives; r++)
3006
{
3007
for (unsigned int s = 0; s < num_derivatives; s++)
3008
{
3009
values[r] += transform[r][s]*derivatives[s];
3010
}// end loop over 's'
3011
}// end loop over 'r'
3012
3013
// Delete pointer to array of derivatives on FIAT element
3014
delete [] derivatives;
3015
3016
// Delete pointer to array of combinations of derivatives and transform
3017
for (unsigned int r = 0; r < num_derivatives; r++)
3018
{
3019
delete [] combinations[r];
3020
}// end loop over 'r'
3021
delete [] combinations;
3022
for (unsigned int r = 0; r < num_derivatives; r++)
3023
{
3024
delete [] transform[r];
3025
}// end loop over 'r'
3026
delete [] transform;
3027
break;
3028
}
3029
}
3030
518
3031
}
519
3032
520
3033
/// Evaluate order n derivatives of all basis functions at given point in cell
Loggerhead is a web-based interface for Breezy
Version: 2.0.1