« back to all changes in this revision
Viewing changes to test/regression/references/SubDomain.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/SubDomain.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
//
124
124
125
125
// Reset values.
126
126
*values = 0.00000000;
127
128
// Map degree of freedom to element degree of freedom
129
const unsigned int dof = i;
130
131
// Array of basisvalues.
132
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
133
134
// Declare helper variables.
135
unsigned int rr = 0;
136
unsigned int ss = 0;
137
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
138
139
// Compute basisvalues.
140
basisvalues[0] = 1.00000000;
141
basisvalues[1] = tmp0;
142
for (unsigned int r = 0; r < 1; r++)
143
{
144
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
145
ss = r*(r + 1)*(r + 2)/6;
146
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
147
}// end loop over 'r'
148
for (unsigned int r = 0; r < 1; r++)
149
{
150
for (unsigned int s = 0; s < 1 - r; s++)
151
{
152
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
153
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
154
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
155
}// end loop over 's'
156
}// end loop over 'r'
157
for (unsigned int r = 0; r < 2; r++)
158
{
159
for (unsigned int s = 0; s < 2 - r; s++)
160
{
161
for (unsigned int t = 0; t < 2 - r - s; t++)
162
{
163
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
164
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
165
}// end loop over 't'
166
}// end loop over 's'
167
}// end loop over 'r'
168
169
// Table(s) of coefficients.
170
static const double coefficients0[4][4] = \
171
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
172
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
173
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
174
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
175
176
// Compute value(s).
177
for (unsigned int r = 0; r < 4; r++)
178
{
179
*values += coefficients0[dof][r]*basisvalues[r];
180
}// end loop over 'r'
127
switch (i)
128
{
129
case 0:
130
{
131
132
// Array of basisvalues.
133
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
134
135
// Declare helper variables.
136
unsigned int rr = 0;
137
unsigned int ss = 0;
138
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
139
140
// Compute basisvalues.
141
basisvalues[0] = 1.00000000;
142
basisvalues[1] = tmp0;
143
for (unsigned int r = 0; r < 1; r++)
144
{
145
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
146
ss = r*(r + 1)*(r + 2)/6;
147
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
148
}// end loop over 'r'
149
for (unsigned int r = 0; r < 1; r++)
150
{
151
for (unsigned int s = 0; s < 1 - r; s++)
152
{
153
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
154
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
155
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
156
}// end loop over 's'
157
}// end loop over 'r'
158
for (unsigned int r = 0; r < 2; r++)
159
{
160
for (unsigned int s = 0; s < 2 - r; s++)
161
{
162
for (unsigned int t = 0; t < 2 - r - s; t++)
163
{
164
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
165
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
166
}// end loop over 't'
167
}// end loop over 's'
168
}// end loop over 'r'
169
170
// Table(s) of coefficients.
171
static const double coefficients0[4] = \
172
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
173
174
// Compute value(s).
175
for (unsigned int r = 0; r < 4; r++)
176
{
177
*values += coefficients0[r]*basisvalues[r];
178
}// end loop over 'r'
179
break;
180
}
181
case 1:
182
{
183
184
// Array of basisvalues.
185
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
186
187
// Declare helper variables.
188
unsigned int rr = 0;
189
unsigned int ss = 0;
190
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
191
192
// Compute basisvalues.
193
basisvalues[0] = 1.00000000;
194
basisvalues[1] = tmp0;
195
for (unsigned int r = 0; r < 1; r++)
196
{
197
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
198
ss = r*(r + 1)*(r + 2)/6;
199
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
200
}// end loop over 'r'
201
for (unsigned int r = 0; r < 1; r++)
202
{
203
for (unsigned int s = 0; s < 1 - r; s++)
204
{
205
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
206
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
207
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
208
}// end loop over 's'
209
}// end loop over 'r'
210
for (unsigned int r = 0; r < 2; r++)
211
{
212
for (unsigned int s = 0; s < 2 - r; s++)
213
{
214
for (unsigned int t = 0; t < 2 - r - s; t++)
215
{
216
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
217
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
218
}// end loop over 't'
219
}// end loop over 's'
220
}// end loop over 'r'
221
222
// Table(s) of coefficients.
223
static const double coefficients0[4] = \
224
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
225
226
// Compute value(s).
227
for (unsigned int r = 0; r < 4; r++)
228
{
229
*values += coefficients0[r]*basisvalues[r];
230
}// end loop over 'r'
231
break;
232
}
233
case 2:
234
{
235
236
// Array of basisvalues.
237
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
238
239
// Declare helper variables.
240
unsigned int rr = 0;
241
unsigned int ss = 0;
242
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
243
244
// Compute basisvalues.
245
basisvalues[0] = 1.00000000;
246
basisvalues[1] = tmp0;
247
for (unsigned int r = 0; r < 1; r++)
248
{
249
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
250
ss = r*(r + 1)*(r + 2)/6;
251
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
252
}// end loop over 'r'
253
for (unsigned int r = 0; r < 1; r++)
254
{
255
for (unsigned int s = 0; s < 1 - r; s++)
256
{
257
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
258
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
259
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
260
}// end loop over 's'
261
}// end loop over 'r'
262
for (unsigned int r = 0; r < 2; r++)
263
{
264
for (unsigned int s = 0; s < 2 - r; s++)
265
{
266
for (unsigned int t = 0; t < 2 - r - s; t++)
267
{
268
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
269
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
270
}// end loop over 't'
271
}// end loop over 's'
272
}// end loop over 'r'
273
274
// Table(s) of coefficients.
275
static const double coefficients0[4] = \
276
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
277
278
// Compute value(s).
279
for (unsigned int r = 0; r < 4; r++)
280
{
281
*values += coefficients0[r]*basisvalues[r];
282
}// end loop over 'r'
283
break;
284
}
285
case 3:
286
{
287
288
// Array of basisvalues.
289
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
290
291
// Declare helper variables.
292
unsigned int rr = 0;
293
unsigned int ss = 0;
294
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
295
296
// Compute basisvalues.
297
basisvalues[0] = 1.00000000;
298
basisvalues[1] = tmp0;
299
for (unsigned int r = 0; r < 1; r++)
300
{
301
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
302
ss = r*(r + 1)*(r + 2)/6;
303
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
304
}// end loop over 'r'
305
for (unsigned int r = 0; r < 1; r++)
306
{
307
for (unsigned int s = 0; s < 1 - r; s++)
308
{
309
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
310
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
311
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
312
}// end loop over 's'
313
}// end loop over 'r'
314
for (unsigned int r = 0; r < 2; r++)
315
{
316
for (unsigned int s = 0; s < 2 - r; s++)
317
{
318
for (unsigned int t = 0; t < 2 - r - s; t++)
319
{
320
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
321
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
322
}// end loop over 't'
323
}// end loop over 's'
324
}// end loop over 'r'
325
326
// Table(s) of coefficients.
327
static const double coefficients0[4] = \
328
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
329
330
// Compute value(s).
331
for (unsigned int r = 0; r < 4; r++)
332
{
333
*values += coefficients0[r]*basisvalues[r];
334
}// end loop over 'r'
335
break;
336
}
337
}
338
181
339
}
182
340
183
341
/// Evaluate all basis functions at given point in cell
317
475
values[r] = 0.00000000;
318
476
}// end loop over 'r'
319
477
320
// Map degree of freedom to element degree of freedom
321
const unsigned int dof = i;
322
323
// Array of basisvalues.
324
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
325
326
// Declare helper variables.
327
unsigned int rr = 0;
328
unsigned int ss = 0;
329
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
330
331
// Compute basisvalues.
332
basisvalues[0] = 1.00000000;
333
basisvalues[1] = tmp0;
334
for (unsigned int r = 0; r < 1; r++)
335
{
336
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
337
ss = r*(r + 1)*(r + 2)/6;
338
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
339
}// end loop over 'r'
340
for (unsigned int r = 0; r < 1; r++)
341
{
342
for (unsigned int s = 0; s < 1 - r; s++)
343
{
344
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
345
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
346
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
347
}// end loop over 's'
348
}// end loop over 'r'
349
for (unsigned int r = 0; r < 2; r++)
350
{
351
for (unsigned int s = 0; s < 2 - r; s++)
352
{
353
for (unsigned int t = 0; t < 2 - r - s; t++)
354
{
355
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
356
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
357
}// end loop over 't'
358
}// end loop over 's'
359
}// end loop over 'r'
360
361
// Table(s) of coefficients.
362
static const double coefficients0[4][4] = \
363
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
364
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
365
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
366
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
367
368
// Tables of derivatives of the polynomial base (transpose).
369
static const double dmats0[4][4] = \
370
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
371
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
372
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
373
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
374
375
static const double dmats1[4][4] = \
376
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
377
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
378
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
379
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
380
381
static const double dmats2[4][4] = \
382
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
383
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
384
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
385
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
386
387
// Compute reference derivatives.
388
// Declare pointer to array of derivatives on FIAT element.
389
double *derivatives = new double[num_derivatives];
390
for (unsigned int r = 0; r < num_derivatives; r++)
391
{
392
derivatives[r] = 0.00000000;
393
}// end loop over 'r'
394
395
// Declare derivative matrix (of polynomial basis).
396
double dmats[4][4] = \
397
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
398
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
399
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
400
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
401
402
// Declare (auxiliary) derivative matrix (of polynomial basis).
403
double dmats_old[4][4] = \
404
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
405
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
406
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
407
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
408
409
// Loop possible derivatives.
410
for (unsigned int r = 0; r < num_derivatives; r++)
411
{
412
// Resetting dmats values to compute next derivative.
413
for (unsigned int t = 0; t < 4; t++)
414
{
415
for (unsigned int u = 0; u < 4; u++)
416
{
417
dmats[t][u] = 0.00000000;
418
if (t == u)
419
{
420
dmats[t][u] = 1.00000000;
421
}
422
423
}// end loop over 'u'
424
}// end loop over 't'
425
426
// Looping derivative order to generate dmats.
427
for (unsigned int s = 0; s < n; s++)
428
{
429
// Updating dmats_old with new values and resetting dmats.
430
for (unsigned int t = 0; t < 4; t++)
431
{
432
for (unsigned int u = 0; u < 4; u++)
433
{
434
dmats_old[t][u] = dmats[t][u];
435
dmats[t][u] = 0.00000000;
436
}// end loop over 'u'
437
}// end loop over 't'
438
439
// Update dmats using an inner product.
440
if (combinations[r][s] == 0)
441
{
442
for (unsigned int t = 0; t < 4; t++)
443
{
444
for (unsigned int u = 0; u < 4; u++)
445
{
446
for (unsigned int tu = 0; tu < 4; tu++)
447
{
448
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
449
}// end loop over 'tu'
450
}// end loop over 'u'
451
}// end loop over 't'
452
}
453
454
if (combinations[r][s] == 1)
455
{
456
for (unsigned int t = 0; t < 4; t++)
457
{
458
for (unsigned int u = 0; u < 4; u++)
459
{
460
for (unsigned int tu = 0; tu < 4; tu++)
461
{
462
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
463
}// end loop over 'tu'
464
}// end loop over 'u'
465
}// end loop over 't'
466
}
467
468
if (combinations[r][s] == 2)
469
{
470
for (unsigned int t = 0; t < 4; t++)
471
{
472
for (unsigned int u = 0; u < 4; u++)
473
{
474
for (unsigned int tu = 0; tu < 4; tu++)
475
{
476
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
477
}// end loop over 'tu'
478
}// end loop over 'u'
479
}// end loop over 't'
480
}
481
482
}// end loop over 's'
483
for (unsigned int s = 0; s < 4; s++)
484
{
485
for (unsigned int t = 0; t < 4; t++)
486
{
487
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
488
}// end loop over 't'
489
}// end loop over 's'
490
}// end loop over 'r'
491
492
// Transform derivatives back to physical element
493
for (unsigned int r = 0; r < num_derivatives; r++)
494
{
495
for (unsigned int s = 0; s < num_derivatives; s++)
496
{
497
values[r] += transform[r][s]*derivatives[s];
498
}// end loop over 's'
499
}// end loop over 'r'
500
501
// Delete pointer to array of derivatives on FIAT element
502
delete [] derivatives;
503
504
// Delete pointer to array of combinations of derivatives and transform
505
for (unsigned int r = 0; r < num_derivatives; r++)
506
{
507
delete [] combinations[r];
508
}// end loop over 'r'
509
delete [] combinations;
510
for (unsigned int r = 0; r < num_derivatives; r++)
511
{
512
delete [] transform[r];
513
}// end loop over 'r'
514
delete [] transform;
478
switch (i)
479
{
480
case 0:
481
{
482
483
// Array of basisvalues.
484
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
485
486
// Declare helper variables.
487
unsigned int rr = 0;
488
unsigned int ss = 0;
489
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
490
491
// Compute basisvalues.
492
basisvalues[0] = 1.00000000;
493
basisvalues[1] = tmp0;
494
for (unsigned int r = 0; r < 1; r++)
495
{
496
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
497
ss = r*(r + 1)*(r + 2)/6;
498
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
499
}// end loop over 'r'
500
for (unsigned int r = 0; r < 1; r++)
501
{
502
for (unsigned int s = 0; s < 1 - r; s++)
503
{
504
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
505
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
506
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
507
}// end loop over 's'
508
}// end loop over 'r'
509
for (unsigned int r = 0; r < 2; r++)
510
{
511
for (unsigned int s = 0; s < 2 - r; s++)
512
{
513
for (unsigned int t = 0; t < 2 - r - s; t++)
514
{
515
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
516
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
517
}// end loop over 't'
518
}// end loop over 's'
519
}// end loop over 'r'
520
521
// Table(s) of coefficients.
522
static const double coefficients0[4] = \
523
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
524
525
// Tables of derivatives of the polynomial base (transpose).
526
static const double dmats0[4][4] = \
527
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
528
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
529
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
530
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
531
532
static const double dmats1[4][4] = \
533
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
534
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
535
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
536
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
537
538
static const double dmats2[4][4] = \
539
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
540
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
541
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
542
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
543
544
// Compute reference derivatives.
545
// Declare pointer to array of derivatives on FIAT element.
546
double *derivatives = new double[num_derivatives];
547
for (unsigned int r = 0; r < num_derivatives; r++)
548
{
549
derivatives[r] = 0.00000000;
550
}// end loop over 'r'
551
552
// Declare derivative matrix (of polynomial basis).
553
double dmats[4][4] = \
554
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
555
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
556
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
557
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
558
559
// Declare (auxiliary) derivative matrix (of polynomial basis).
560
double dmats_old[4][4] = \
561
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
562
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
563
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
564
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
565
566
// Loop possible derivatives.
567
for (unsigned int r = 0; r < num_derivatives; r++)
568
{
569
// Resetting dmats values to compute next derivative.
570
for (unsigned int t = 0; t < 4; t++)
571
{
572
for (unsigned int u = 0; u < 4; u++)
573
{
574
dmats[t][u] = 0.00000000;
575
if (t == u)
576
{
577
dmats[t][u] = 1.00000000;
578
}
579
580
}// end loop over 'u'
581
}// end loop over 't'
582
583
// Looping derivative order to generate dmats.
584
for (unsigned int s = 0; s < n; s++)
585
{
586
// Updating dmats_old with new values and resetting dmats.
587
for (unsigned int t = 0; t < 4; t++)
588
{
589
for (unsigned int u = 0; u < 4; u++)
590
{
591
dmats_old[t][u] = dmats[t][u];
592
dmats[t][u] = 0.00000000;
593
}// end loop over 'u'
594
}// end loop over 't'
595
596
// Update dmats using an inner product.
597
if (combinations[r][s] == 0)
598
{
599
for (unsigned int t = 0; t < 4; t++)
600
{
601
for (unsigned int u = 0; u < 4; u++)
602
{
603
for (unsigned int tu = 0; tu < 4; tu++)
604
{
605
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
606
}// end loop over 'tu'
607
}// end loop over 'u'
608
}// end loop over 't'
609
}
610
611
if (combinations[r][s] == 1)
612
{
613
for (unsigned int t = 0; t < 4; t++)
614
{
615
for (unsigned int u = 0; u < 4; u++)
616
{
617
for (unsigned int tu = 0; tu < 4; tu++)
618
{
619
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
620
}// end loop over 'tu'
621
}// end loop over 'u'
622
}// end loop over 't'
623
}
624
625
if (combinations[r][s] == 2)
626
{
627
for (unsigned int t = 0; t < 4; t++)
628
{
629
for (unsigned int u = 0; u < 4; u++)
630
{
631
for (unsigned int tu = 0; tu < 4; tu++)
632
{
633
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
634
}// end loop over 'tu'
635
}// end loop over 'u'
636
}// end loop over 't'
637
}
638
639
}// end loop over 's'
640
for (unsigned int s = 0; s < 4; s++)
641
{
642
for (unsigned int t = 0; t < 4; t++)
643
{
644
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
645
}// end loop over 't'
646
}// end loop over 's'
647
}// end loop over 'r'
648
649
// Transform derivatives back to physical element
650
for (unsigned int r = 0; r < num_derivatives; r++)
651
{
652
for (unsigned int s = 0; s < num_derivatives; s++)
653
{
654
values[r] += transform[r][s]*derivatives[s];
655
}// end loop over 's'
656
}// end loop over 'r'
657
658
// Delete pointer to array of derivatives on FIAT element
659
delete [] derivatives;
660
661
// Delete pointer to array of combinations of derivatives and transform
662
for (unsigned int r = 0; r < num_derivatives; r++)
663
{
664
delete [] combinations[r];
665
}// end loop over 'r'
666
delete [] combinations;
667
for (unsigned int r = 0; r < num_derivatives; r++)
668
{
669
delete [] transform[r];
670
}// end loop over 'r'
671
delete [] transform;
672
break;
673
}
674
case 1:
675
{
676
677
// Array of basisvalues.
678
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
679
680
// Declare helper variables.
681
unsigned int rr = 0;
682
unsigned int ss = 0;
683
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
684
685
// Compute basisvalues.
686
basisvalues[0] = 1.00000000;
687
basisvalues[1] = tmp0;
688
for (unsigned int r = 0; r < 1; r++)
689
{
690
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
691
ss = r*(r + 1)*(r + 2)/6;
692
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
693
}// end loop over 'r'
694
for (unsigned int r = 0; r < 1; r++)
695
{
696
for (unsigned int s = 0; s < 1 - r; s++)
697
{
698
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
699
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
700
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
701
}// end loop over 's'
702
}// end loop over 'r'
703
for (unsigned int r = 0; r < 2; r++)
704
{
705
for (unsigned int s = 0; s < 2 - r; s++)
706
{
707
for (unsigned int t = 0; t < 2 - r - s; t++)
708
{
709
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
710
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
711
}// end loop over 't'
712
}// end loop over 's'
713
}// end loop over 'r'
714
715
// Table(s) of coefficients.
716
static const double coefficients0[4] = \
717
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
718
719
// Tables of derivatives of the polynomial base (transpose).
720
static const double dmats0[4][4] = \
721
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
722
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
723
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
724
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
725
726
static const double dmats1[4][4] = \
727
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
728
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
729
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
730
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
731
732
static const double dmats2[4][4] = \
733
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
734
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
735
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
736
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
737
738
// Compute reference derivatives.
739
// Declare pointer to array of derivatives on FIAT element.
740
double *derivatives = new double[num_derivatives];
741
for (unsigned int r = 0; r < num_derivatives; r++)
742
{
743
derivatives[r] = 0.00000000;
744
}// end loop over 'r'
745
746
// Declare derivative matrix (of polynomial basis).
747
double dmats[4][4] = \
748
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
749
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
750
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
751
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
752
753
// Declare (auxiliary) derivative matrix (of polynomial basis).
754
double dmats_old[4][4] = \
755
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
756
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
757
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
758
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
759
760
// Loop possible derivatives.
761
for (unsigned int r = 0; r < num_derivatives; r++)
762
{
763
// Resetting dmats values to compute next derivative.
764
for (unsigned int t = 0; t < 4; t++)
765
{
766
for (unsigned int u = 0; u < 4; u++)
767
{
768
dmats[t][u] = 0.00000000;
769
if (t == u)
770
{
771
dmats[t][u] = 1.00000000;
772
}
773
774
}// end loop over 'u'
775
}// end loop over 't'
776
777
// Looping derivative order to generate dmats.
778
for (unsigned int s = 0; s < n; s++)
779
{
780
// Updating dmats_old with new values and resetting dmats.
781
for (unsigned int t = 0; t < 4; t++)
782
{
783
for (unsigned int u = 0; u < 4; u++)
784
{
785
dmats_old[t][u] = dmats[t][u];
786
dmats[t][u] = 0.00000000;
787
}// end loop over 'u'
788
}// end loop over 't'
789
790
// Update dmats using an inner product.
791
if (combinations[r][s] == 0)
792
{
793
for (unsigned int t = 0; t < 4; t++)
794
{
795
for (unsigned int u = 0; u < 4; u++)
796
{
797
for (unsigned int tu = 0; tu < 4; tu++)
798
{
799
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
800
}// end loop over 'tu'
801
}// end loop over 'u'
802
}// end loop over 't'
803
}
804
805
if (combinations[r][s] == 1)
806
{
807
for (unsigned int t = 0; t < 4; t++)
808
{
809
for (unsigned int u = 0; u < 4; u++)
810
{
811
for (unsigned int tu = 0; tu < 4; tu++)
812
{
813
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
814
}// end loop over 'tu'
815
}// end loop over 'u'
816
}// end loop over 't'
817
}
818
819
if (combinations[r][s] == 2)
820
{
821
for (unsigned int t = 0; t < 4; t++)
822
{
823
for (unsigned int u = 0; u < 4; u++)
824
{
825
for (unsigned int tu = 0; tu < 4; tu++)
826
{
827
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
828
}// end loop over 'tu'
829
}// end loop over 'u'
830
}// end loop over 't'
831
}
832
833
}// end loop over 's'
834
for (unsigned int s = 0; s < 4; s++)
835
{
836
for (unsigned int t = 0; t < 4; t++)
837
{
838
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
839
}// end loop over 't'
840
}// end loop over 's'
841
}// end loop over 'r'
842
843
// Transform derivatives back to physical element
844
for (unsigned int r = 0; r < num_derivatives; r++)
845
{
846
for (unsigned int s = 0; s < num_derivatives; s++)
847
{
848
values[r] += transform[r][s]*derivatives[s];
849
}// end loop over 's'
850
}// end loop over 'r'
851
852
// Delete pointer to array of derivatives on FIAT element
853
delete [] derivatives;
854
855
// Delete pointer to array of combinations of derivatives and transform
856
for (unsigned int r = 0; r < num_derivatives; r++)
857
{
858
delete [] combinations[r];
859
}// end loop over 'r'
860
delete [] combinations;
861
for (unsigned int r = 0; r < num_derivatives; r++)
862
{
863
delete [] transform[r];
864
}// end loop over 'r'
865
delete [] transform;
866
break;
867
}
868
case 2:
869
{
870
871
// Array of basisvalues.
872
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
873
874
// Declare helper variables.
875
unsigned int rr = 0;
876
unsigned int ss = 0;
877
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
878
879
// Compute basisvalues.
880
basisvalues[0] = 1.00000000;
881
basisvalues[1] = tmp0;
882
for (unsigned int r = 0; r < 1; r++)
883
{
884
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
885
ss = r*(r + 1)*(r + 2)/6;
886
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
887
}// end loop over 'r'
888
for (unsigned int r = 0; r < 1; r++)
889
{
890
for (unsigned int s = 0; s < 1 - r; s++)
891
{
892
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
893
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
894
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
895
}// end loop over 's'
896
}// end loop over 'r'
897
for (unsigned int r = 0; r < 2; r++)
898
{
899
for (unsigned int s = 0; s < 2 - r; s++)
900
{
901
for (unsigned int t = 0; t < 2 - r - s; t++)
902
{
903
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
904
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
905
}// end loop over 't'
906
}// end loop over 's'
907
}// end loop over 'r'
908
909
// Table(s) of coefficients.
910
static const double coefficients0[4] = \
911
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
912
913
// Tables of derivatives of the polynomial base (transpose).
914
static const double dmats0[4][4] = \
915
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
916
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
917
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
918
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
919
920
static const double dmats1[4][4] = \
921
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
922
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
923
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
924
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
925
926
static const double dmats2[4][4] = \
927
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
928
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
929
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
930
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
931
932
// Compute reference derivatives.
933
// Declare pointer to array of derivatives on FIAT element.
934
double *derivatives = new double[num_derivatives];
935
for (unsigned int r = 0; r < num_derivatives; r++)
936
{
937
derivatives[r] = 0.00000000;
938
}// end loop over 'r'
939
940
// Declare derivative matrix (of polynomial basis).
941
double dmats[4][4] = \
942
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
943
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
944
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
945
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
946
947
// Declare (auxiliary) derivative matrix (of polynomial basis).
948
double dmats_old[4][4] = \
949
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
950
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
951
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
952
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
953
954
// Loop possible derivatives.
955
for (unsigned int r = 0; r < num_derivatives; r++)
956
{
957
// Resetting dmats values to compute next derivative.
958
for (unsigned int t = 0; t < 4; t++)
959
{
960
for (unsigned int u = 0; u < 4; u++)
961
{
962
dmats[t][u] = 0.00000000;
963
if (t == u)
964
{
965
dmats[t][u] = 1.00000000;
966
}
967
968
}// end loop over 'u'
969
}// end loop over 't'
970
971
// Looping derivative order to generate dmats.
972
for (unsigned int s = 0; s < n; s++)
973
{
974
// Updating dmats_old with new values and resetting dmats.
975
for (unsigned int t = 0; t < 4; t++)
976
{
977
for (unsigned int u = 0; u < 4; u++)
978
{
979
dmats_old[t][u] = dmats[t][u];
980
dmats[t][u] = 0.00000000;
981
}// end loop over 'u'
982
}// end loop over 't'
983
984
// Update dmats using an inner product.
985
if (combinations[r][s] == 0)
986
{
987
for (unsigned int t = 0; t < 4; t++)
988
{
989
for (unsigned int u = 0; u < 4; u++)
990
{
991
for (unsigned int tu = 0; tu < 4; tu++)
992
{
993
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
994
}// end loop over 'tu'
995
}// end loop over 'u'
996
}// end loop over 't'
997
}
998
999
if (combinations[r][s] == 1)
1000
{
1001
for (unsigned int t = 0; t < 4; t++)
1002
{
1003
for (unsigned int u = 0; u < 4; u++)
1004
{
1005
for (unsigned int tu = 0; tu < 4; tu++)
1006
{
1007
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1008
}// end loop over 'tu'
1009
}// end loop over 'u'
1010
}// end loop over 't'
1011
}
1012
1013
if (combinations[r][s] == 2)
1014
{
1015
for (unsigned int t = 0; t < 4; t++)
1016
{
1017
for (unsigned int u = 0; u < 4; u++)
1018
{
1019
for (unsigned int tu = 0; tu < 4; tu++)
1020
{
1021
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1022
}// end loop over 'tu'
1023
}// end loop over 'u'
1024
}// end loop over 't'
1025
}
1026
1027
}// end loop over 's'
1028
for (unsigned int s = 0; s < 4; s++)
1029
{
1030
for (unsigned int t = 0; t < 4; t++)
1031
{
1032
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1033
}// end loop over 't'
1034
}// end loop over 's'
1035
}// end loop over 'r'
1036
1037
// Transform derivatives back to physical element
1038
for (unsigned int r = 0; r < num_derivatives; r++)
1039
{
1040
for (unsigned int s = 0; s < num_derivatives; s++)
1041
{
1042
values[r] += transform[r][s]*derivatives[s];
1043
}// end loop over 's'
1044
}// end loop over 'r'
1045
1046
// Delete pointer to array of derivatives on FIAT element
1047
delete [] derivatives;
1048
1049
// Delete pointer to array of combinations of derivatives and transform
1050
for (unsigned int r = 0; r < num_derivatives; r++)
1051
{
1052
delete [] combinations[r];
1053
}// end loop over 'r'
1054
delete [] combinations;
1055
for (unsigned int r = 0; r < num_derivatives; r++)
1056
{
1057
delete [] transform[r];
1058
}// end loop over 'r'
1059
delete [] transform;
1060
break;
1061
}
1062
case 3:
1063
{
1064
1065
// Array of basisvalues.
1066
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1067
1068
// Declare helper variables.
1069
unsigned int rr = 0;
1070
unsigned int ss = 0;
1071
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1072
1073
// Compute basisvalues.
1074
basisvalues[0] = 1.00000000;
1075
basisvalues[1] = tmp0;
1076
for (unsigned int r = 0; r < 1; r++)
1077
{
1078
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1079
ss = r*(r + 1)*(r + 2)/6;
1080
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1081
}// end loop over 'r'
1082
for (unsigned int r = 0; r < 1; r++)
1083
{
1084
for (unsigned int s = 0; s < 1 - r; s++)
1085
{
1086
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1087
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1088
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1089
}// end loop over 's'
1090
}// end loop over 'r'
1091
for (unsigned int r = 0; r < 2; r++)
1092
{
1093
for (unsigned int s = 0; s < 2 - r; s++)
1094
{
1095
for (unsigned int t = 0; t < 2 - r - s; t++)
1096
{
1097
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1098
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1099
}// end loop over 't'
1100
}// end loop over 's'
1101
}// end loop over 'r'
1102
1103
// Table(s) of coefficients.
1104
static const double coefficients0[4] = \
1105
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
1106
1107
// Tables of derivatives of the polynomial base (transpose).
1108
static const double dmats0[4][4] = \
1109
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1110
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
1111
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1112
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1113
1114
static const double dmats1[4][4] = \
1115
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1116
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1117
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
1118
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1119
1120
static const double dmats2[4][4] = \
1121
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1122
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1123
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
1124
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
1125
1126
// Compute reference derivatives.
1127
// Declare pointer to array of derivatives on FIAT element.
1128
double *derivatives = new double[num_derivatives];
1129
for (unsigned int r = 0; r < num_derivatives; r++)
1130
{
1131
derivatives[r] = 0.00000000;
1132
}// end loop over 'r'
1133
1134
// Declare derivative matrix (of polynomial basis).
1135
double dmats[4][4] = \
1136
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1137
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1138
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1139
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1140
1141
// Declare (auxiliary) derivative matrix (of polynomial basis).
1142
double dmats_old[4][4] = \
1143
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1144
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1145
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1146
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1147
1148
// Loop possible derivatives.
1149
for (unsigned int r = 0; r < num_derivatives; r++)
1150
{
1151
// Resetting dmats values to compute next derivative.
1152
for (unsigned int t = 0; t < 4; t++)
1153
{
1154
for (unsigned int u = 0; u < 4; u++)
1155
{
1156
dmats[t][u] = 0.00000000;
1157
if (t == u)
1158
{
1159
dmats[t][u] = 1.00000000;
1160
}
1161
1162
}// end loop over 'u'
1163
}// end loop over 't'
1164
1165
// Looping derivative order to generate dmats.
1166
for (unsigned int s = 0; s < n; s++)
1167
{
1168
// Updating dmats_old with new values and resetting dmats.
1169
for (unsigned int t = 0; t < 4; t++)
1170
{
1171
for (unsigned int u = 0; u < 4; u++)
1172
{
1173
dmats_old[t][u] = dmats[t][u];
1174
dmats[t][u] = 0.00000000;
1175
}// end loop over 'u'
1176
}// end loop over 't'
1177
1178
// Update dmats using an inner product.
1179
if (combinations[r][s] == 0)
1180
{
1181
for (unsigned int t = 0; t < 4; t++)
1182
{
1183
for (unsigned int u = 0; u < 4; u++)
1184
{
1185
for (unsigned int tu = 0; tu < 4; tu++)
1186
{
1187
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1188
}// end loop over 'tu'
1189
}// end loop over 'u'
1190
}// end loop over 't'
1191
}
1192
1193
if (combinations[r][s] == 1)
1194
{
1195
for (unsigned int t = 0; t < 4; t++)
1196
{
1197
for (unsigned int u = 0; u < 4; u++)
1198
{
1199
for (unsigned int tu = 0; tu < 4; tu++)
1200
{
1201
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1202
}// end loop over 'tu'
1203
}// end loop over 'u'
1204
}// end loop over 't'
1205
}
1206
1207
if (combinations[r][s] == 2)
1208
{
1209
for (unsigned int t = 0; t < 4; t++)
1210
{
1211
for (unsigned int u = 0; u < 4; u++)
1212
{
1213
for (unsigned int tu = 0; tu < 4; tu++)
1214
{
1215
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1216
}// end loop over 'tu'
1217
}// end loop over 'u'
1218
}// end loop over 't'
1219
}
1220
1221
}// end loop over 's'
1222
for (unsigned int s = 0; s < 4; s++)
1223
{
1224
for (unsigned int t = 0; t < 4; t++)
1225
{
1226
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1227
}// end loop over 't'
1228
}// end loop over 's'
1229
}// end loop over 'r'
1230
1231
// Transform derivatives back to physical element
1232
for (unsigned int r = 0; r < num_derivatives; r++)
1233
{
1234
for (unsigned int s = 0; s < num_derivatives; s++)
1235
{
1236
values[r] += transform[r][s]*derivatives[s];
1237
}// end loop over 's'
1238
}// end loop over 'r'
1239
1240
// Delete pointer to array of derivatives on FIAT element
1241
delete [] derivatives;
1242
1243
// Delete pointer to array of combinations of derivatives and transform
1244
for (unsigned int r = 0; r < num_derivatives; r++)
1245
{
1246
delete [] combinations[r];
1247
}// end loop over 'r'
1248
delete [] combinations;
1249
for (unsigned int r = 0; r < num_derivatives; r++)
1250
{
1251
delete [] transform[r];
1252
}// end loop over 'r'
1253
delete [] transform;
1254
break;
1255
}
1256
}
1257
515
1258
}
516
1259
517
1260
/// 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