~ubuntu-branches/ubuntu/raring/ffc/raring
» Revision
15
: test/regression/references/r_quadrature_O/NavierStokes.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_quadrature_O/NavierStokes.h
- browse files at revision 15
- compare with another revision
- download diff
- download tarball
- view history from revision 15
- Committer: Package Import Robot
- Author(s): Johannes Ring
- Date: 2011-11-29 11:38:38 UTC
- mfrom: (1.1.11)
- Revision ID: package-import@ubuntu.com-20111129113838-vx815cxjdf50zuq3
Tags: 1.0-rc1-1
New upstream release.
- files added:
- doc/sphinx
- doc/sphinx/scripts
- doc/sphinx/source
- doc/sphinx/source/_static
- doc/sphinx/source/_templates
- test/regression/references/output
- test/regression/references/r_quadrature_O
- files removed:
- test/regression/references/r_auto/AdaptivePoisson.out
- files modified:
- ChangeLog
added
removed
test/regression/references/r_quadrature_O/NavierStokes.h
1
// This code conforms with the UFC specification version 2.0.3
2
// and was automatically generated by FFC version 1.0-beta2+.
3
//
4
// This code was generated with the following parameters:
5
//
6
// cache_dir: ''
7
// convert_exceptions_to_warnings: True
8
// cpp_optimize: False
9
// cpp_optimize_flags: '-O2'
10
// epsilon: 1e-14
11
// error_control: False
12
// form_postfix: True
13
// format: 'ufc'
14
// log_level: 20
15
// log_prefix: ''
16
// optimize: True
17
// output_dir: '.'
18
// precision: '8'
19
// quadrature_degree: 'auto'
20
// quadrature_rule: 'auto'
21
// representation: 'quadrature'
22
// split: False
23
// swig_binary: 'swig'
24
// swig_path: ''
25
26
#ifndef __NAVIERSTOKES_H
27
#define __NAVIERSTOKES_H
28
29
#include <cmath>
30
#include <stdexcept>
31
#include <fstream>
32
#include <ufc.h>
33
34
/// This class defines the interface for a finite element.
35
36
class navierstokes_finite_element_0: public ufc::finite_element
37
{
38
public:
39
40
/// Constructor
41
navierstokes_finite_element_0() : ufc::finite_element()
42
{
43
// Do nothing
44
}
45
46
/// Destructor
47
virtual ~navierstokes_finite_element_0()
48
{
49
// Do nothing
50
}
51
52
/// Return a string identifying the finite element
53
virtual const char* signature() const
54
{
55
return "FiniteElement('Lagrange', Cell('tetrahedron', Space(3)), 1, None)";
56
}
57
58
/// Return the cell shape
59
virtual ufc::shape cell_shape() const
60
{
61
return ufc::tetrahedron;
62
}
63
64
/// Return the topological dimension of the cell shape
65
virtual unsigned int topological_dimension() const
66
{
67
return 3;
68
}
69
70
/// Return the geometric dimension of the cell shape
71
virtual unsigned int geometric_dimension() const
72
{
73
return 3;
74
}
75
76
/// Return the dimension of the finite element function space
77
virtual unsigned int space_dimension() const
78
{
79
return 4;
80
}
81
82
/// Return the rank of the value space
83
virtual unsigned int value_rank() const
84
{
85
return 0;
86
}
87
88
/// Return the dimension of the value space for axis i
89
virtual unsigned int value_dimension(unsigned int i) const
90
{
91
return 1;
92
}
93
94
/// Evaluate basis function i at given point in cell
95
virtual void evaluate_basis(unsigned int i,
96
double* values,
97
const double* coordinates,
98
const ufc::cell& c) const
99
{
100
// Extract vertex coordinates
101
const double * const * x = c.coordinates;
102
103
// Compute Jacobian of affine map from reference cell
104
const double J_00 = x[1][0] - x[0][0];
105
const double J_01 = x[2][0] - x[0][0];
106
const double J_02 = x[3][0] - x[0][0];
107
const double J_10 = x[1][1] - x[0][1];
108
const double J_11 = x[2][1] - x[0][1];
109
const double J_12 = x[3][1] - x[0][1];
110
const double J_20 = x[1][2] - x[0][2];
111
const double J_21 = x[2][2] - x[0][2];
112
const double J_22 = x[3][2] - x[0][2];
113
114
// Compute sub determinants
115
const double d_00 = J_11*J_22 - J_12*J_21;
116
const double d_01 = J_12*J_20 - J_10*J_22;
117
const double d_02 = J_10*J_21 - J_11*J_20;
118
const double d_10 = J_02*J_21 - J_01*J_22;
119
const double d_11 = J_00*J_22 - J_02*J_20;
120
const double d_12 = J_01*J_20 - J_00*J_21;
121
const double d_20 = J_01*J_12 - J_02*J_11;
122
const double d_21 = J_02*J_10 - J_00*J_12;
123
const double d_22 = J_00*J_11 - J_01*J_10;
124
125
// Compute determinant of Jacobian
126
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
127
128
// Compute inverse of Jacobian
129
130
// Compute constants
131
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
132
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
133
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
134
135
// Get coordinates and map to the reference (FIAT) element
136
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
137
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
138
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
139
140
141
// Reset values.
142
*values = 0.0;
143
switch (i)
144
{
145
case 0:
146
{
147
148
// Array of basisvalues.
149
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
150
151
// Declare helper variables.
152
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
153
154
// Compute basisvalues.
155
basisvalues[0] = 1.0;
156
basisvalues[1] = tmp0;
157
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
158
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
159
basisvalues[0] *= std::sqrt(0.75);
160
basisvalues[3] *= std::sqrt(1.25);
161
basisvalues[2] *= std::sqrt(2.5);
162
basisvalues[1] *= std::sqrt(7.5);
163
164
// Table(s) of coefficients.
165
static const double coefficients0[4] = \
166
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
167
168
// Compute value(s).
169
for (unsigned int r = 0; r < 4; r++)
170
{
171
*values += coefficients0[r]*basisvalues[r];
172
}// end loop over 'r'
173
break;
174
}
175
case 1:
176
{
177
178
// Array of basisvalues.
179
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
180
181
// Declare helper variables.
182
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
183
184
// Compute basisvalues.
185
basisvalues[0] = 1.0;
186
basisvalues[1] = tmp0;
187
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
188
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
189
basisvalues[0] *= std::sqrt(0.75);
190
basisvalues[3] *= std::sqrt(1.25);
191
basisvalues[2] *= std::sqrt(2.5);
192
basisvalues[1] *= std::sqrt(7.5);
193
194
// Table(s) of coefficients.
195
static const double coefficients0[4] = \
196
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
197
198
// Compute value(s).
199
for (unsigned int r = 0; r < 4; r++)
200
{
201
*values += coefficients0[r]*basisvalues[r];
202
}// end loop over 'r'
203
break;
204
}
205
case 2:
206
{
207
208
// Array of basisvalues.
209
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
210
211
// Declare helper variables.
212
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
213
214
// Compute basisvalues.
215
basisvalues[0] = 1.0;
216
basisvalues[1] = tmp0;
217
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
218
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
219
basisvalues[0] *= std::sqrt(0.75);
220
basisvalues[3] *= std::sqrt(1.25);
221
basisvalues[2] *= std::sqrt(2.5);
222
basisvalues[1] *= std::sqrt(7.5);
223
224
// Table(s) of coefficients.
225
static const double coefficients0[4] = \
226
{0.28867513, 0.0, 0.21081851, -0.074535599};
227
228
// Compute value(s).
229
for (unsigned int r = 0; r < 4; r++)
230
{
231
*values += coefficients0[r]*basisvalues[r];
232
}// end loop over 'r'
233
break;
234
}
235
case 3:
236
{
237
238
// Array of basisvalues.
239
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
240
241
// Declare helper variables.
242
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
243
244
// Compute basisvalues.
245
basisvalues[0] = 1.0;
246
basisvalues[1] = tmp0;
247
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
248
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
249
basisvalues[0] *= std::sqrt(0.75);
250
basisvalues[3] *= std::sqrt(1.25);
251
basisvalues[2] *= std::sqrt(2.5);
252
basisvalues[1] *= std::sqrt(7.5);
253
254
// Table(s) of coefficients.
255
static const double coefficients0[4] = \
256
{0.28867513, 0.0, 0.0, 0.2236068};
257
258
// Compute value(s).
259
for (unsigned int r = 0; r < 4; r++)
260
{
261
*values += coefficients0[r]*basisvalues[r];
262
}// end loop over 'r'
263
break;
264
}
265
}
266
267
}
268
269
/// Evaluate all basis functions at given point in cell
270
virtual void evaluate_basis_all(double* values,
271
const double* coordinates,
272
const ufc::cell& c) const
273
{
274
// Helper variable to hold values of a single dof.
275
double dof_values = 0.0;
276
277
// Loop dofs and call evaluate_basis.
278
for (unsigned int r = 0; r < 4; r++)
279
{
280
evaluate_basis(r, &dof_values, coordinates, c);
281
values[r] = dof_values;
282
}// end loop over 'r'
283
}
284
285
/// Evaluate order n derivatives of basis function i at given point in cell
286
virtual void evaluate_basis_derivatives(unsigned int i,
287
unsigned int n,
288
double* values,
289
const double* coordinates,
290
const ufc::cell& c) const
291
{
292
// Extract vertex coordinates
293
const double * const * x = c.coordinates;
294
295
// Compute Jacobian of affine map from reference cell
296
const double J_00 = x[1][0] - x[0][0];
297
const double J_01 = x[2][0] - x[0][0];
298
const double J_02 = x[3][0] - x[0][0];
299
const double J_10 = x[1][1] - x[0][1];
300
const double J_11 = x[2][1] - x[0][1];
301
const double J_12 = x[3][1] - x[0][1];
302
const double J_20 = x[1][2] - x[0][2];
303
const double J_21 = x[2][2] - x[0][2];
304
const double J_22 = x[3][2] - x[0][2];
305
306
// Compute sub determinants
307
const double d_00 = J_11*J_22 - J_12*J_21;
308
const double d_01 = J_12*J_20 - J_10*J_22;
309
const double d_02 = J_10*J_21 - J_11*J_20;
310
const double d_10 = J_02*J_21 - J_01*J_22;
311
const double d_11 = J_00*J_22 - J_02*J_20;
312
const double d_12 = J_01*J_20 - J_00*J_21;
313
const double d_20 = J_01*J_12 - J_02*J_11;
314
const double d_21 = J_02*J_10 - J_00*J_12;
315
const double d_22 = J_00*J_11 - J_01*J_10;
316
317
// Compute determinant of Jacobian
318
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
319
320
// Compute inverse of Jacobian
321
const double K_00 = d_00 / detJ;
322
const double K_01 = d_10 / detJ;
323
const double K_02 = d_20 / detJ;
324
const double K_10 = d_01 / detJ;
325
const double K_11 = d_11 / detJ;
326
const double K_12 = d_21 / detJ;
327
const double K_20 = d_02 / detJ;
328
const double K_21 = d_12 / detJ;
329
const double K_22 = d_22 / detJ;
330
331
// Compute constants
332
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
333
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
334
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
335
336
// Get coordinates and map to the reference (FIAT) element
337
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
338
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
339
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
340
341
342
// Compute number of derivatives.
343
unsigned int num_derivatives = 1;
344
for (unsigned int r = 0; r < n; r++)
345
{
346
num_derivatives *= 3;
347
}// end loop over 'r'
348
349
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
350
unsigned int **combinations = new unsigned int *[num_derivatives];
351
for (unsigned int row = 0; row < num_derivatives; row++)
352
{
353
combinations[row] = new unsigned int [n];
354
for (unsigned int col = 0; col < n; col++)
355
combinations[row][col] = 0;
356
}
357
358
// Generate combinations of derivatives
359
for (unsigned int row = 1; row < num_derivatives; row++)
360
{
361
for (unsigned int num = 0; num < row; num++)
362
{
363
for (unsigned int col = n-1; col+1 > 0; col--)
364
{
365
if (combinations[row][col] + 1 > 2)
366
combinations[row][col] = 0;
367
else
368
{
369
combinations[row][col] += 1;
370
break;
371
}
372
}
373
}
374
}
375
376
// Compute inverse of Jacobian
377
const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}};
378
379
// Declare transformation matrix
380
// Declare pointer to two dimensional array and initialise
381
double **transform = new double *[num_derivatives];
382
383
for (unsigned int j = 0; j < num_derivatives; j++)
384
{
385
transform[j] = new double [num_derivatives];
386
for (unsigned int k = 0; k < num_derivatives; k++)
387
transform[j][k] = 1;
388
}
389
390
// Construct transformation matrix
391
for (unsigned int row = 0; row < num_derivatives; row++)
392
{
393
for (unsigned int col = 0; col < num_derivatives; col++)
394
{
395
for (unsigned int k = 0; k < n; k++)
396
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
397
}
398
}
399
400
// Reset values. Assuming that values is always an array.
401
for (unsigned int r = 0; r < num_derivatives; r++)
402
{
403
values[r] = 0.0;
404
}// end loop over 'r'
405
406
switch (i)
407
{
408
case 0:
409
{
410
411
// Array of basisvalues.
412
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
413
414
// Declare helper variables.
415
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
416
417
// Compute basisvalues.
418
basisvalues[0] = 1.0;
419
basisvalues[1] = tmp0;
420
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
421
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
422
basisvalues[0] *= std::sqrt(0.75);
423
basisvalues[3] *= std::sqrt(1.25);
424
basisvalues[2] *= std::sqrt(2.5);
425
basisvalues[1] *= std::sqrt(7.5);
426
427
// Table(s) of coefficients.
428
static const double coefficients0[4] = \
429
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
430
431
// Tables of derivatives of the polynomial base (transpose).
432
static const double dmats0[4][4] = \
433
{{0.0, 0.0, 0.0, 0.0},
434
{6.3245553, 0.0, 0.0, 0.0},
435
{0.0, 0.0, 0.0, 0.0},
436
{0.0, 0.0, 0.0, 0.0}};
437
438
static const double dmats1[4][4] = \
439
{{0.0, 0.0, 0.0, 0.0},
440
{3.1622777, 0.0, 0.0, 0.0},
441
{5.4772256, 0.0, 0.0, 0.0},
442
{0.0, 0.0, 0.0, 0.0}};
443
444
static const double dmats2[4][4] = \
445
{{0.0, 0.0, 0.0, 0.0},
446
{3.1622777, 0.0, 0.0, 0.0},
447
{1.8257419, 0.0, 0.0, 0.0},
448
{5.1639778, 0.0, 0.0, 0.0}};
449
450
// Compute reference derivatives.
451
// Declare pointer to array of derivatives on FIAT element.
452
double *derivatives = new double[num_derivatives];
453
for (unsigned int r = 0; r < num_derivatives; r++)
454
{
455
derivatives[r] = 0.0;
456
}// end loop over 'r'
457
458
// Declare derivative matrix (of polynomial basis).
459
double dmats[4][4] = \
460
{{1.0, 0.0, 0.0, 0.0},
461
{0.0, 1.0, 0.0, 0.0},
462
{0.0, 0.0, 1.0, 0.0},
463
{0.0, 0.0, 0.0, 1.0}};
464
465
// Declare (auxiliary) derivative matrix (of polynomial basis).
466
double dmats_old[4][4] = \
467
{{1.0, 0.0, 0.0, 0.0},
468
{0.0, 1.0, 0.0, 0.0},
469
{0.0, 0.0, 1.0, 0.0},
470
{0.0, 0.0, 0.0, 1.0}};
471
472
// Loop possible derivatives.
473
for (unsigned int r = 0; r < num_derivatives; r++)
474
{
475
// Resetting dmats values to compute next derivative.
476
for (unsigned int t = 0; t < 4; t++)
477
{
478
for (unsigned int u = 0; u < 4; u++)
479
{
480
dmats[t][u] = 0.0;
481
if (t == u)
482
{
483
dmats[t][u] = 1.0;
484
}
485
486
}// end loop over 'u'
487
}// end loop over 't'
488
489
// Looping derivative order to generate dmats.
490
for (unsigned int s = 0; s < n; s++)
491
{
492
// Updating dmats_old with new values and resetting dmats.
493
for (unsigned int t = 0; t < 4; t++)
494
{
495
for (unsigned int u = 0; u < 4; u++)
496
{
497
dmats_old[t][u] = dmats[t][u];
498
dmats[t][u] = 0.0;
499
}// end loop over 'u'
500
}// end loop over 't'
501
502
// Update dmats using an inner product.
503
if (combinations[r][s] == 0)
504
{
505
for (unsigned int t = 0; t < 4; t++)
506
{
507
for (unsigned int u = 0; u < 4; u++)
508
{
509
for (unsigned int tu = 0; tu < 4; tu++)
510
{
511
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
512
}// end loop over 'tu'
513
}// end loop over 'u'
514
}// end loop over 't'
515
}
516
517
if (combinations[r][s] == 1)
518
{
519
for (unsigned int t = 0; t < 4; t++)
520
{
521
for (unsigned int u = 0; u < 4; u++)
522
{
523
for (unsigned int tu = 0; tu < 4; tu++)
524
{
525
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
526
}// end loop over 'tu'
527
}// end loop over 'u'
528
}// end loop over 't'
529
}
530
531
if (combinations[r][s] == 2)
532
{
533
for (unsigned int t = 0; t < 4; t++)
534
{
535
for (unsigned int u = 0; u < 4; u++)
536
{
537
for (unsigned int tu = 0; tu < 4; tu++)
538
{
539
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
540
}// end loop over 'tu'
541
}// end loop over 'u'
542
}// end loop over 't'
543
}
544
545
}// end loop over 's'
546
for (unsigned int s = 0; s < 4; s++)
547
{
548
for (unsigned int t = 0; t < 4; t++)
549
{
550
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
551
}// end loop over 't'
552
}// end loop over 's'
553
}// end loop over 'r'
554
555
// Transform derivatives back to physical element
556
for (unsigned int r = 0; r < num_derivatives; r++)
557
{
558
for (unsigned int s = 0; s < num_derivatives; s++)
559
{
560
values[r] += transform[r][s]*derivatives[s];
561
}// end loop over 's'
562
}// end loop over 'r'
563
564
// Delete pointer to array of derivatives on FIAT element
565
delete [] derivatives;
566
567
// Delete pointer to array of combinations of derivatives and transform
568
for (unsigned int r = 0; r < num_derivatives; r++)
569
{
570
delete [] combinations[r];
571
}// end loop over 'r'
572
delete [] combinations;
573
for (unsigned int r = 0; r < num_derivatives; r++)
574
{
575
delete [] transform[r];
576
}// end loop over 'r'
577
delete [] transform;
578
break;
579
}
580
case 1:
581
{
582
583
// Array of basisvalues.
584
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
585
586
// Declare helper variables.
587
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
588
589
// Compute basisvalues.
590
basisvalues[0] = 1.0;
591
basisvalues[1] = tmp0;
592
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
593
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
594
basisvalues[0] *= std::sqrt(0.75);
595
basisvalues[3] *= std::sqrt(1.25);
596
basisvalues[2] *= std::sqrt(2.5);
597
basisvalues[1] *= std::sqrt(7.5);
598
599
// Table(s) of coefficients.
600
static const double coefficients0[4] = \
601
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
602
603
// Tables of derivatives of the polynomial base (transpose).
604
static const double dmats0[4][4] = \
605
{{0.0, 0.0, 0.0, 0.0},
606
{6.3245553, 0.0, 0.0, 0.0},
607
{0.0, 0.0, 0.0, 0.0},
608
{0.0, 0.0, 0.0, 0.0}};
609
610
static const double dmats1[4][4] = \
611
{{0.0, 0.0, 0.0, 0.0},
612
{3.1622777, 0.0, 0.0, 0.0},
613
{5.4772256, 0.0, 0.0, 0.0},
614
{0.0, 0.0, 0.0, 0.0}};
615
616
static const double dmats2[4][4] = \
617
{{0.0, 0.0, 0.0, 0.0},
618
{3.1622777, 0.0, 0.0, 0.0},
619
{1.8257419, 0.0, 0.0, 0.0},
620
{5.1639778, 0.0, 0.0, 0.0}};
621
622
// Compute reference derivatives.
623
// Declare pointer to array of derivatives on FIAT element.
624
double *derivatives = new double[num_derivatives];
625
for (unsigned int r = 0; r < num_derivatives; r++)
626
{
627
derivatives[r] = 0.0;
628
}// end loop over 'r'
629
630
// Declare derivative matrix (of polynomial basis).
631
double dmats[4][4] = \
632
{{1.0, 0.0, 0.0, 0.0},
633
{0.0, 1.0, 0.0, 0.0},
634
{0.0, 0.0, 1.0, 0.0},
635
{0.0, 0.0, 0.0, 1.0}};
636
637
// Declare (auxiliary) derivative matrix (of polynomial basis).
638
double dmats_old[4][4] = \
639
{{1.0, 0.0, 0.0, 0.0},
640
{0.0, 1.0, 0.0, 0.0},
641
{0.0, 0.0, 1.0, 0.0},
642
{0.0, 0.0, 0.0, 1.0}};
643
644
// Loop possible derivatives.
645
for (unsigned int r = 0; r < num_derivatives; r++)
646
{
647
// Resetting dmats values to compute next derivative.
648
for (unsigned int t = 0; t < 4; t++)
649
{
650
for (unsigned int u = 0; u < 4; u++)
651
{
652
dmats[t][u] = 0.0;
653
if (t == u)
654
{
655
dmats[t][u] = 1.0;
656
}
657
658
}// end loop over 'u'
659
}// end loop over 't'
660
661
// Looping derivative order to generate dmats.
662
for (unsigned int s = 0; s < n; s++)
663
{
664
// Updating dmats_old with new values and resetting dmats.
665
for (unsigned int t = 0; t < 4; t++)
666
{
667
for (unsigned int u = 0; u < 4; u++)
668
{
669
dmats_old[t][u] = dmats[t][u];
670
dmats[t][u] = 0.0;
671
}// end loop over 'u'
672
}// end loop over 't'
673
674
// Update dmats using an inner product.
675
if (combinations[r][s] == 0)
676
{
677
for (unsigned int t = 0; t < 4; t++)
678
{
679
for (unsigned int u = 0; u < 4; u++)
680
{
681
for (unsigned int tu = 0; tu < 4; tu++)
682
{
683
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
684
}// end loop over 'tu'
685
}// end loop over 'u'
686
}// end loop over 't'
687
}
688
689
if (combinations[r][s] == 1)
690
{
691
for (unsigned int t = 0; t < 4; t++)
692
{
693
for (unsigned int u = 0; u < 4; u++)
694
{
695
for (unsigned int tu = 0; tu < 4; tu++)
696
{
697
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
698
}// end loop over 'tu'
699
}// end loop over 'u'
700
}// end loop over 't'
701
}
702
703
if (combinations[r][s] == 2)
704
{
705
for (unsigned int t = 0; t < 4; t++)
706
{
707
for (unsigned int u = 0; u < 4; u++)
708
{
709
for (unsigned int tu = 0; tu < 4; tu++)
710
{
711
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
712
}// end loop over 'tu'
713
}// end loop over 'u'
714
}// end loop over 't'
715
}
716
717
}// end loop over 's'
718
for (unsigned int s = 0; s < 4; s++)
719
{
720
for (unsigned int t = 0; t < 4; t++)
721
{
722
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
723
}// end loop over 't'
724
}// end loop over 's'
725
}// end loop over 'r'
726
727
// Transform derivatives back to physical element
728
for (unsigned int r = 0; r < num_derivatives; r++)
729
{
730
for (unsigned int s = 0; s < num_derivatives; s++)
731
{
732
values[r] += transform[r][s]*derivatives[s];
733
}// end loop over 's'
734
}// end loop over 'r'
735
736
// Delete pointer to array of derivatives on FIAT element
737
delete [] derivatives;
738
739
// Delete pointer to array of combinations of derivatives and transform
740
for (unsigned int r = 0; r < num_derivatives; r++)
741
{
742
delete [] combinations[r];
743
}// end loop over 'r'
744
delete [] combinations;
745
for (unsigned int r = 0; r < num_derivatives; r++)
746
{
747
delete [] transform[r];
748
}// end loop over 'r'
749
delete [] transform;
750
break;
751
}
752
case 2:
753
{
754
755
// Array of basisvalues.
756
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
757
758
// Declare helper variables.
759
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
760
761
// Compute basisvalues.
762
basisvalues[0] = 1.0;
763
basisvalues[1] = tmp0;
764
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
765
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
766
basisvalues[0] *= std::sqrt(0.75);
767
basisvalues[3] *= std::sqrt(1.25);
768
basisvalues[2] *= std::sqrt(2.5);
769
basisvalues[1] *= std::sqrt(7.5);
770
771
// Table(s) of coefficients.
772
static const double coefficients0[4] = \
773
{0.28867513, 0.0, 0.21081851, -0.074535599};
774
775
// Tables of derivatives of the polynomial base (transpose).
776
static const double dmats0[4][4] = \
777
{{0.0, 0.0, 0.0, 0.0},
778
{6.3245553, 0.0, 0.0, 0.0},
779
{0.0, 0.0, 0.0, 0.0},
780
{0.0, 0.0, 0.0, 0.0}};
781
782
static const double dmats1[4][4] = \
783
{{0.0, 0.0, 0.0, 0.0},
784
{3.1622777, 0.0, 0.0, 0.0},
785
{5.4772256, 0.0, 0.0, 0.0},
786
{0.0, 0.0, 0.0, 0.0}};
787
788
static const double dmats2[4][4] = \
789
{{0.0, 0.0, 0.0, 0.0},
790
{3.1622777, 0.0, 0.0, 0.0},
791
{1.8257419, 0.0, 0.0, 0.0},
792
{5.1639778, 0.0, 0.0, 0.0}};
793
794
// Compute reference derivatives.
795
// Declare pointer to array of derivatives on FIAT element.
796
double *derivatives = new double[num_derivatives];
797
for (unsigned int r = 0; r < num_derivatives; r++)
798
{
799
derivatives[r] = 0.0;
800
}// end loop over 'r'
801
802
// Declare derivative matrix (of polynomial basis).
803
double dmats[4][4] = \
804
{{1.0, 0.0, 0.0, 0.0},
805
{0.0, 1.0, 0.0, 0.0},
806
{0.0, 0.0, 1.0, 0.0},
807
{0.0, 0.0, 0.0, 1.0}};
808
809
// Declare (auxiliary) derivative matrix (of polynomial basis).
810
double dmats_old[4][4] = \
811
{{1.0, 0.0, 0.0, 0.0},
812
{0.0, 1.0, 0.0, 0.0},
813
{0.0, 0.0, 1.0, 0.0},
814
{0.0, 0.0, 0.0, 1.0}};
815
816
// Loop possible derivatives.
817
for (unsigned int r = 0; r < num_derivatives; r++)
818
{
819
// Resetting dmats values to compute next derivative.
820
for (unsigned int t = 0; t < 4; t++)
821
{
822
for (unsigned int u = 0; u < 4; u++)
823
{
824
dmats[t][u] = 0.0;
825
if (t == u)
826
{
827
dmats[t][u] = 1.0;
828
}
829
830
}// end loop over 'u'
831
}// end loop over 't'
832
833
// Looping derivative order to generate dmats.
834
for (unsigned int s = 0; s < n; s++)
835
{
836
// Updating dmats_old with new values and resetting dmats.
837
for (unsigned int t = 0; t < 4; t++)
838
{
839
for (unsigned int u = 0; u < 4; u++)
840
{
841
dmats_old[t][u] = dmats[t][u];
842
dmats[t][u] = 0.0;
843
}// end loop over 'u'
844
}// end loop over 't'
845
846
// Update dmats using an inner product.
847
if (combinations[r][s] == 0)
848
{
849
for (unsigned int t = 0; t < 4; t++)
850
{
851
for (unsigned int u = 0; u < 4; u++)
852
{
853
for (unsigned int tu = 0; tu < 4; tu++)
854
{
855
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
856
}// end loop over 'tu'
857
}// end loop over 'u'
858
}// end loop over 't'
859
}
860
861
if (combinations[r][s] == 1)
862
{
863
for (unsigned int t = 0; t < 4; t++)
864
{
865
for (unsigned int u = 0; u < 4; u++)
866
{
867
for (unsigned int tu = 0; tu < 4; tu++)
868
{
869
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
870
}// end loop over 'tu'
871
}// end loop over 'u'
872
}// end loop over 't'
873
}
874
875
if (combinations[r][s] == 2)
876
{
877
for (unsigned int t = 0; t < 4; t++)
878
{
879
for (unsigned int u = 0; u < 4; u++)
880
{
881
for (unsigned int tu = 0; tu < 4; tu++)
882
{
883
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
884
}// end loop over 'tu'
885
}// end loop over 'u'
886
}// end loop over 't'
887
}
888
889
}// end loop over 's'
890
for (unsigned int s = 0; s < 4; s++)
891
{
892
for (unsigned int t = 0; t < 4; t++)
893
{
894
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
895
}// end loop over 't'
896
}// end loop over 's'
897
}// end loop over 'r'
898
899
// Transform derivatives back to physical element
900
for (unsigned int r = 0; r < num_derivatives; r++)
901
{
902
for (unsigned int s = 0; s < num_derivatives; s++)
903
{
904
values[r] += transform[r][s]*derivatives[s];
905
}// end loop over 's'
906
}// end loop over 'r'
907
908
// Delete pointer to array of derivatives on FIAT element
909
delete [] derivatives;
910
911
// Delete pointer to array of combinations of derivatives and transform
912
for (unsigned int r = 0; r < num_derivatives; r++)
913
{
914
delete [] combinations[r];
915
}// end loop over 'r'
916
delete [] combinations;
917
for (unsigned int r = 0; r < num_derivatives; r++)
918
{
919
delete [] transform[r];
920
}// end loop over 'r'
921
delete [] transform;
922
break;
923
}
924
case 3:
925
{
926
927
// Array of basisvalues.
928
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
929
930
// Declare helper variables.
931
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
932
933
// Compute basisvalues.
934
basisvalues[0] = 1.0;
935
basisvalues[1] = tmp0;
936
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
937
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
938
basisvalues[0] *= std::sqrt(0.75);
939
basisvalues[3] *= std::sqrt(1.25);
940
basisvalues[2] *= std::sqrt(2.5);
941
basisvalues[1] *= std::sqrt(7.5);
942
943
// Table(s) of coefficients.
944
static const double coefficients0[4] = \
945
{0.28867513, 0.0, 0.0, 0.2236068};
946
947
// Tables of derivatives of the polynomial base (transpose).
948
static const double dmats0[4][4] = \
949
{{0.0, 0.0, 0.0, 0.0},
950
{6.3245553, 0.0, 0.0, 0.0},
951
{0.0, 0.0, 0.0, 0.0},
952
{0.0, 0.0, 0.0, 0.0}};
953
954
static const double dmats1[4][4] = \
955
{{0.0, 0.0, 0.0, 0.0},
956
{3.1622777, 0.0, 0.0, 0.0},
957
{5.4772256, 0.0, 0.0, 0.0},
958
{0.0, 0.0, 0.0, 0.0}};
959
960
static const double dmats2[4][4] = \
961
{{0.0, 0.0, 0.0, 0.0},
962
{3.1622777, 0.0, 0.0, 0.0},
963
{1.8257419, 0.0, 0.0, 0.0},
964
{5.1639778, 0.0, 0.0, 0.0}};
965
966
// Compute reference derivatives.
967
// Declare pointer to array of derivatives on FIAT element.
968
double *derivatives = new double[num_derivatives];
969
for (unsigned int r = 0; r < num_derivatives; r++)
970
{
971
derivatives[r] = 0.0;
972
}// end loop over 'r'
973
974
// Declare derivative matrix (of polynomial basis).
975
double dmats[4][4] = \
976
{{1.0, 0.0, 0.0, 0.0},
977
{0.0, 1.0, 0.0, 0.0},
978
{0.0, 0.0, 1.0, 0.0},
979
{0.0, 0.0, 0.0, 1.0}};
980
981
// Declare (auxiliary) derivative matrix (of polynomial basis).
982
double dmats_old[4][4] = \
983
{{1.0, 0.0, 0.0, 0.0},
984
{0.0, 1.0, 0.0, 0.0},
985
{0.0, 0.0, 1.0, 0.0},
986
{0.0, 0.0, 0.0, 1.0}};
987
988
// Loop possible derivatives.
989
for (unsigned int r = 0; r < num_derivatives; r++)
990
{
991
// Resetting dmats values to compute next derivative.
992
for (unsigned int t = 0; t < 4; t++)
993
{
994
for (unsigned int u = 0; u < 4; u++)
995
{
996
dmats[t][u] = 0.0;
997
if (t == u)
998
{
999
dmats[t][u] = 1.0;
1000
}
1001
1002
}// end loop over 'u'
1003
}// end loop over 't'
1004
1005
// Looping derivative order to generate dmats.
1006
for (unsigned int s = 0; s < n; s++)
1007
{
1008
// Updating dmats_old with new values and resetting dmats.
1009
for (unsigned int t = 0; t < 4; t++)
1010
{
1011
for (unsigned int u = 0; u < 4; u++)
1012
{
1013
dmats_old[t][u] = dmats[t][u];
1014
dmats[t][u] = 0.0;
1015
}// end loop over 'u'
1016
}// end loop over 't'
1017
1018
// Update dmats using an inner product.
1019
if (combinations[r][s] == 0)
1020
{
1021
for (unsigned int t = 0; t < 4; t++)
1022
{
1023
for (unsigned int u = 0; u < 4; u++)
1024
{
1025
for (unsigned int tu = 0; tu < 4; tu++)
1026
{
1027
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1028
}// end loop over 'tu'
1029
}// end loop over 'u'
1030
}// end loop over 't'
1031
}
1032
1033
if (combinations[r][s] == 1)
1034
{
1035
for (unsigned int t = 0; t < 4; t++)
1036
{
1037
for (unsigned int u = 0; u < 4; u++)
1038
{
1039
for (unsigned int tu = 0; tu < 4; tu++)
1040
{
1041
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1042
}// end loop over 'tu'
1043
}// end loop over 'u'
1044
}// end loop over 't'
1045
}
1046
1047
if (combinations[r][s] == 2)
1048
{
1049
for (unsigned int t = 0; t < 4; t++)
1050
{
1051
for (unsigned int u = 0; u < 4; u++)
1052
{
1053
for (unsigned int tu = 0; tu < 4; tu++)
1054
{
1055
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1056
}// end loop over 'tu'
1057
}// end loop over 'u'
1058
}// end loop over 't'
1059
}
1060
1061
}// end loop over 's'
1062
for (unsigned int s = 0; s < 4; s++)
1063
{
1064
for (unsigned int t = 0; t < 4; t++)
1065
{
1066
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1067
}// end loop over 't'
1068
}// end loop over 's'
1069
}// end loop over 'r'
1070
1071
// Transform derivatives back to physical element
1072
for (unsigned int r = 0; r < num_derivatives; r++)
1073
{
1074
for (unsigned int s = 0; s < num_derivatives; s++)
1075
{
1076
values[r] += transform[r][s]*derivatives[s];
1077
}// end loop over 's'
1078
}// end loop over 'r'
1079
1080
// Delete pointer to array of derivatives on FIAT element
1081
delete [] derivatives;
1082
1083
// Delete pointer to array of combinations of derivatives and transform
1084
for (unsigned int r = 0; r < num_derivatives; r++)
1085
{
1086
delete [] combinations[r];
1087
}// end loop over 'r'
1088
delete [] combinations;
1089
for (unsigned int r = 0; r < num_derivatives; r++)
1090
{
1091
delete [] transform[r];
1092
}// end loop over 'r'
1093
delete [] transform;
1094
break;
1095
}
1096
}
1097
1098
}
1099
1100
/// Evaluate order n derivatives of all basis functions at given point in cell
1101
virtual void evaluate_basis_derivatives_all(unsigned int n,
1102
double* values,
1103
const double* coordinates,
1104
const ufc::cell& c) const
1105
{
1106
// Compute number of derivatives.
1107
unsigned int num_derivatives = 1;
1108
for (unsigned int r = 0; r < n; r++)
1109
{
1110
num_derivatives *= 3;
1111
}// end loop over 'r'
1112
1113
// Helper variable to hold values of a single dof.
1114
double *dof_values = new double[num_derivatives];
1115
for (unsigned int r = 0; r < num_derivatives; r++)
1116
{
1117
dof_values[r] = 0.0;
1118
}// end loop over 'r'
1119
1120
// Loop dofs and call evaluate_basis_derivatives.
1121
for (unsigned int r = 0; r < 4; r++)
1122
{
1123
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
1124
for (unsigned int s = 0; s < num_derivatives; s++)
1125
{
1126
values[r*num_derivatives + s] = dof_values[s];
1127
}// end loop over 's'
1128
}// end loop over 'r'
1129
1130
// Delete pointer.
1131
delete [] dof_values;
1132
}
1133
1134
/// Evaluate linear functional for dof i on the function f
1135
virtual double evaluate_dof(unsigned int i,
1136
const ufc::function& f,
1137
const ufc::cell& c) const
1138
{
1139
// Declare variables for result of evaluation.
1140
double vals[1];
1141
1142
// Declare variable for physical coordinates.
1143
double y[3];
1144
const double * const * x = c.coordinates;
1145
switch (i)
1146
{
1147
case 0:
1148
{
1149
y[0] = x[0][0];
1150
y[1] = x[0][1];
1151
y[2] = x[0][2];
1152
f.evaluate(vals, y, c);
1153
return vals[0];
1154
break;
1155
}
1156
case 1:
1157
{
1158
y[0] = x[1][0];
1159
y[1] = x[1][1];
1160
y[2] = x[1][2];
1161
f.evaluate(vals, y, c);
1162
return vals[0];
1163
break;
1164
}
1165
case 2:
1166
{
1167
y[0] = x[2][0];
1168
y[1] = x[2][1];
1169
y[2] = x[2][2];
1170
f.evaluate(vals, y, c);
1171
return vals[0];
1172
break;
1173
}
1174
case 3:
1175
{
1176
y[0] = x[3][0];
1177
y[1] = x[3][1];
1178
y[2] = x[3][2];
1179
f.evaluate(vals, y, c);
1180
return vals[0];
1181
break;
1182
}
1183
}
1184
1185
return 0.0;
1186
}
1187
1188
/// Evaluate linear functionals for all dofs on the function f
1189
virtual void evaluate_dofs(double* values,
1190
const ufc::function& f,
1191
const ufc::cell& c) const
1192
{
1193
// Declare variables for result of evaluation.
1194
double vals[1];
1195
1196
// Declare variable for physical coordinates.
1197
double y[3];
1198
const double * const * x = c.coordinates;
1199
y[0] = x[0][0];
1200
y[1] = x[0][1];
1201
y[2] = x[0][2];
1202
f.evaluate(vals, y, c);
1203
values[0] = vals[0];
1204
y[0] = x[1][0];
1205
y[1] = x[1][1];
1206
y[2] = x[1][2];
1207
f.evaluate(vals, y, c);
1208
values[1] = vals[0];
1209
y[0] = x[2][0];
1210
y[1] = x[2][1];
1211
y[2] = x[2][2];
1212
f.evaluate(vals, y, c);
1213
values[2] = vals[0];
1214
y[0] = x[3][0];
1215
y[1] = x[3][1];
1216
y[2] = x[3][2];
1217
f.evaluate(vals, y, c);
1218
values[3] = vals[0];
1219
}
1220
1221
/// Interpolate vertex values from dof values
1222
virtual void interpolate_vertex_values(double* vertex_values,
1223
const double* dof_values,
1224
const ufc::cell& c) const
1225
{
1226
// Evaluate function and change variables
1227
vertex_values[0] = dof_values[0];
1228
vertex_values[1] = dof_values[1];
1229
vertex_values[2] = dof_values[2];
1230
vertex_values[3] = dof_values[3];
1231
}
1232
1233
/// Map coordinate xhat from reference cell to coordinate x in cell
1234
virtual void map_from_reference_cell(double* x,
1235
const double* xhat,
1236
const ufc::cell& c) const
1237
{
1238
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
1239
}
1240
1241
/// Map from coordinate x in cell to coordinate xhat in reference cell
1242
virtual void map_to_reference_cell(double* xhat,
1243
const double* x,
1244
const ufc::cell& c) const
1245
{
1246
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
1247
}
1248
1249
/// Return the number of sub elements (for a mixed element)
1250
virtual unsigned int num_sub_elements() const
1251
{
1252
return 0;
1253
}
1254
1255
/// Create a new finite element for sub element i (for a mixed element)
1256
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1257
{
1258
return 0;
1259
}
1260
1261
/// Create a new class instance
1262
virtual ufc::finite_element* create() const
1263
{
1264
return new navierstokes_finite_element_0();
1265
}
1266
1267
};
1268
1269
/// This class defines the interface for a finite element.
1270
1271
class navierstokes_finite_element_1: public ufc::finite_element
1272
{
1273
public:
1274
1275
/// Constructor
1276
navierstokes_finite_element_1() : ufc::finite_element()
1277
{
1278
// Do nothing
1279
}
1280
1281
/// Destructor
1282
virtual ~navierstokes_finite_element_1()
1283
{
1284
// Do nothing
1285
}
1286
1287
/// Return a string identifying the finite element
1288
virtual const char* signature() const
1289
{
1290
return "VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None)";
1291
}
1292
1293
/// Return the cell shape
1294
virtual ufc::shape cell_shape() const
1295
{
1296
return ufc::tetrahedron;
1297
}
1298
1299
/// Return the topological dimension of the cell shape
1300
virtual unsigned int topological_dimension() const
1301
{
1302
return 3;
1303
}
1304
1305
/// Return the geometric dimension of the cell shape
1306
virtual unsigned int geometric_dimension() const
1307
{
1308
return 3;
1309
}
1310
1311
/// Return the dimension of the finite element function space
1312
virtual unsigned int space_dimension() const
1313
{
1314
return 12;
1315
}
1316
1317
/// Return the rank of the value space
1318
virtual unsigned int value_rank() const
1319
{
1320
return 1;
1321
}
1322
1323
/// Return the dimension of the value space for axis i
1324
virtual unsigned int value_dimension(unsigned int i) const
1325
{
1326
switch (i)
1327
{
1328
case 0:
1329
{
1330
return 3;
1331
break;
1332
}
1333
}
1334
1335
return 0;
1336
}
1337
1338
/// Evaluate basis function i at given point in cell
1339
virtual void evaluate_basis(unsigned int i,
1340
double* values,
1341
const double* coordinates,
1342
const ufc::cell& c) const
1343
{
1344
// Extract vertex coordinates
1345
const double * const * x = c.coordinates;
1346
1347
// Compute Jacobian of affine map from reference cell
1348
const double J_00 = x[1][0] - x[0][0];
1349
const double J_01 = x[2][0] - x[0][0];
1350
const double J_02 = x[3][0] - x[0][0];
1351
const double J_10 = x[1][1] - x[0][1];
1352
const double J_11 = x[2][1] - x[0][1];
1353
const double J_12 = x[3][1] - x[0][1];
1354
const double J_20 = x[1][2] - x[0][2];
1355
const double J_21 = x[2][2] - x[0][2];
1356
const double J_22 = x[3][2] - x[0][2];
1357
1358
// Compute sub determinants
1359
const double d_00 = J_11*J_22 - J_12*J_21;
1360
const double d_01 = J_12*J_20 - J_10*J_22;
1361
const double d_02 = J_10*J_21 - J_11*J_20;
1362
const double d_10 = J_02*J_21 - J_01*J_22;
1363
const double d_11 = J_00*J_22 - J_02*J_20;
1364
const double d_12 = J_01*J_20 - J_00*J_21;
1365
const double d_20 = J_01*J_12 - J_02*J_11;
1366
const double d_21 = J_02*J_10 - J_00*J_12;
1367
const double d_22 = J_00*J_11 - J_01*J_10;
1368
1369
// Compute determinant of Jacobian
1370
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
1371
1372
// Compute inverse of Jacobian
1373
1374
// Compute constants
1375
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
1376
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
1377
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
1378
1379
// Get coordinates and map to the reference (FIAT) element
1380
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
1381
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
1382
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
1383
1384
1385
// Reset values.
1386
values[0] = 0.0;
1387
values[1] = 0.0;
1388
values[2] = 0.0;
1389
switch (i)
1390
{
1391
case 0:
1392
{
1393
1394
// Array of basisvalues.
1395
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1396
1397
// Declare helper variables.
1398
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1399
1400
// Compute basisvalues.
1401
basisvalues[0] = 1.0;
1402
basisvalues[1] = tmp0;
1403
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1404
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1405
basisvalues[0] *= std::sqrt(0.75);
1406
basisvalues[3] *= std::sqrt(1.25);
1407
basisvalues[2] *= std::sqrt(2.5);
1408
basisvalues[1] *= std::sqrt(7.5);
1409
1410
// Table(s) of coefficients.
1411
static const double coefficients0[4] = \
1412
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1413
1414
// Compute value(s).
1415
for (unsigned int r = 0; r < 4; r++)
1416
{
1417
values[0] += coefficients0[r]*basisvalues[r];
1418
}// end loop over 'r'
1419
break;
1420
}
1421
case 1:
1422
{
1423
1424
// Array of basisvalues.
1425
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1426
1427
// Declare helper variables.
1428
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1429
1430
// Compute basisvalues.
1431
basisvalues[0] = 1.0;
1432
basisvalues[1] = tmp0;
1433
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1434
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1435
basisvalues[0] *= std::sqrt(0.75);
1436
basisvalues[3] *= std::sqrt(1.25);
1437
basisvalues[2] *= std::sqrt(2.5);
1438
basisvalues[1] *= std::sqrt(7.5);
1439
1440
// Table(s) of coefficients.
1441
static const double coefficients0[4] = \
1442
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
1443
1444
// Compute value(s).
1445
for (unsigned int r = 0; r < 4; r++)
1446
{
1447
values[0] += coefficients0[r]*basisvalues[r];
1448
}// end loop over 'r'
1449
break;
1450
}
1451
case 2:
1452
{
1453
1454
// Array of basisvalues.
1455
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1456
1457
// Declare helper variables.
1458
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1459
1460
// Compute basisvalues.
1461
basisvalues[0] = 1.0;
1462
basisvalues[1] = tmp0;
1463
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1464
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1465
basisvalues[0] *= std::sqrt(0.75);
1466
basisvalues[3] *= std::sqrt(1.25);
1467
basisvalues[2] *= std::sqrt(2.5);
1468
basisvalues[1] *= std::sqrt(7.5);
1469
1470
// Table(s) of coefficients.
1471
static const double coefficients0[4] = \
1472
{0.28867513, 0.0, 0.21081851, -0.074535599};
1473
1474
// Compute value(s).
1475
for (unsigned int r = 0; r < 4; r++)
1476
{
1477
values[0] += coefficients0[r]*basisvalues[r];
1478
}// end loop over 'r'
1479
break;
1480
}
1481
case 3:
1482
{
1483
1484
// Array of basisvalues.
1485
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1486
1487
// Declare helper variables.
1488
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1489
1490
// Compute basisvalues.
1491
basisvalues[0] = 1.0;
1492
basisvalues[1] = tmp0;
1493
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1494
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1495
basisvalues[0] *= std::sqrt(0.75);
1496
basisvalues[3] *= std::sqrt(1.25);
1497
basisvalues[2] *= std::sqrt(2.5);
1498
basisvalues[1] *= std::sqrt(7.5);
1499
1500
// Table(s) of coefficients.
1501
static const double coefficients0[4] = \
1502
{0.28867513, 0.0, 0.0, 0.2236068};
1503
1504
// Compute value(s).
1505
for (unsigned int r = 0; r < 4; r++)
1506
{
1507
values[0] += coefficients0[r]*basisvalues[r];
1508
}// end loop over 'r'
1509
break;
1510
}
1511
case 4:
1512
{
1513
1514
// Array of basisvalues.
1515
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1516
1517
// Declare helper variables.
1518
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1519
1520
// Compute basisvalues.
1521
basisvalues[0] = 1.0;
1522
basisvalues[1] = tmp0;
1523
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1524
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1525
basisvalues[0] *= std::sqrt(0.75);
1526
basisvalues[3] *= std::sqrt(1.25);
1527
basisvalues[2] *= std::sqrt(2.5);
1528
basisvalues[1] *= std::sqrt(7.5);
1529
1530
// Table(s) of coefficients.
1531
static const double coefficients0[4] = \
1532
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1533
1534
// Compute value(s).
1535
for (unsigned int r = 0; r < 4; r++)
1536
{
1537
values[1] += coefficients0[r]*basisvalues[r];
1538
}// end loop over 'r'
1539
break;
1540
}
1541
case 5:
1542
{
1543
1544
// Array of basisvalues.
1545
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1546
1547
// Declare helper variables.
1548
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1549
1550
// Compute basisvalues.
1551
basisvalues[0] = 1.0;
1552
basisvalues[1] = tmp0;
1553
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1554
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1555
basisvalues[0] *= std::sqrt(0.75);
1556
basisvalues[3] *= std::sqrt(1.25);
1557
basisvalues[2] *= std::sqrt(2.5);
1558
basisvalues[1] *= std::sqrt(7.5);
1559
1560
// Table(s) of coefficients.
1561
static const double coefficients0[4] = \
1562
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
1563
1564
// Compute value(s).
1565
for (unsigned int r = 0; r < 4; r++)
1566
{
1567
values[1] += coefficients0[r]*basisvalues[r];
1568
}// end loop over 'r'
1569
break;
1570
}
1571
case 6:
1572
{
1573
1574
// Array of basisvalues.
1575
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1576
1577
// Declare helper variables.
1578
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1579
1580
// Compute basisvalues.
1581
basisvalues[0] = 1.0;
1582
basisvalues[1] = tmp0;
1583
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1584
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1585
basisvalues[0] *= std::sqrt(0.75);
1586
basisvalues[3] *= std::sqrt(1.25);
1587
basisvalues[2] *= std::sqrt(2.5);
1588
basisvalues[1] *= std::sqrt(7.5);
1589
1590
// Table(s) of coefficients.
1591
static const double coefficients0[4] = \
1592
{0.28867513, 0.0, 0.21081851, -0.074535599};
1593
1594
// Compute value(s).
1595
for (unsigned int r = 0; r < 4; r++)
1596
{
1597
values[1] += coefficients0[r]*basisvalues[r];
1598
}// end loop over 'r'
1599
break;
1600
}
1601
case 7:
1602
{
1603
1604
// Array of basisvalues.
1605
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1606
1607
// Declare helper variables.
1608
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1609
1610
// Compute basisvalues.
1611
basisvalues[0] = 1.0;
1612
basisvalues[1] = tmp0;
1613
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1614
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1615
basisvalues[0] *= std::sqrt(0.75);
1616
basisvalues[3] *= std::sqrt(1.25);
1617
basisvalues[2] *= std::sqrt(2.5);
1618
basisvalues[1] *= std::sqrt(7.5);
1619
1620
// Table(s) of coefficients.
1621
static const double coefficients0[4] = \
1622
{0.28867513, 0.0, 0.0, 0.2236068};
1623
1624
// Compute value(s).
1625
for (unsigned int r = 0; r < 4; r++)
1626
{
1627
values[1] += coefficients0[r]*basisvalues[r];
1628
}// end loop over 'r'
1629
break;
1630
}
1631
case 8:
1632
{
1633
1634
// Array of basisvalues.
1635
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1636
1637
// Declare helper variables.
1638
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1639
1640
// Compute basisvalues.
1641
basisvalues[0] = 1.0;
1642
basisvalues[1] = tmp0;
1643
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1644
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1645
basisvalues[0] *= std::sqrt(0.75);
1646
basisvalues[3] *= std::sqrt(1.25);
1647
basisvalues[2] *= std::sqrt(2.5);
1648
basisvalues[1] *= std::sqrt(7.5);
1649
1650
// Table(s) of coefficients.
1651
static const double coefficients0[4] = \
1652
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1653
1654
// Compute value(s).
1655
for (unsigned int r = 0; r < 4; r++)
1656
{
1657
values[2] += coefficients0[r]*basisvalues[r];
1658
}// end loop over 'r'
1659
break;
1660
}
1661
case 9:
1662
{
1663
1664
// Array of basisvalues.
1665
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1666
1667
// Declare helper variables.
1668
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1669
1670
// Compute basisvalues.
1671
basisvalues[0] = 1.0;
1672
basisvalues[1] = tmp0;
1673
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1674
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1675
basisvalues[0] *= std::sqrt(0.75);
1676
basisvalues[3] *= std::sqrt(1.25);
1677
basisvalues[2] *= std::sqrt(2.5);
1678
basisvalues[1] *= std::sqrt(7.5);
1679
1680
// Table(s) of coefficients.
1681
static const double coefficients0[4] = \
1682
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
1683
1684
// Compute value(s).
1685
for (unsigned int r = 0; r < 4; r++)
1686
{
1687
values[2] += coefficients0[r]*basisvalues[r];
1688
}// end loop over 'r'
1689
break;
1690
}
1691
case 10:
1692
{
1693
1694
// Array of basisvalues.
1695
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1696
1697
// Declare helper variables.
1698
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1699
1700
// Compute basisvalues.
1701
basisvalues[0] = 1.0;
1702
basisvalues[1] = tmp0;
1703
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1704
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1705
basisvalues[0] *= std::sqrt(0.75);
1706
basisvalues[3] *= std::sqrt(1.25);
1707
basisvalues[2] *= std::sqrt(2.5);
1708
basisvalues[1] *= std::sqrt(7.5);
1709
1710
// Table(s) of coefficients.
1711
static const double coefficients0[4] = \
1712
{0.28867513, 0.0, 0.21081851, -0.074535599};
1713
1714
// Compute value(s).
1715
for (unsigned int r = 0; r < 4; r++)
1716
{
1717
values[2] += coefficients0[r]*basisvalues[r];
1718
}// end loop over 'r'
1719
break;
1720
}
1721
case 11:
1722
{
1723
1724
// Array of basisvalues.
1725
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1726
1727
// Declare helper variables.
1728
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1729
1730
// Compute basisvalues.
1731
basisvalues[0] = 1.0;
1732
basisvalues[1] = tmp0;
1733
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1734
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1735
basisvalues[0] *= std::sqrt(0.75);
1736
basisvalues[3] *= std::sqrt(1.25);
1737
basisvalues[2] *= std::sqrt(2.5);
1738
basisvalues[1] *= std::sqrt(7.5);
1739
1740
// Table(s) of coefficients.
1741
static const double coefficients0[4] = \
1742
{0.28867513, 0.0, 0.0, 0.2236068};
1743
1744
// Compute value(s).
1745
for (unsigned int r = 0; r < 4; r++)
1746
{
1747
values[2] += coefficients0[r]*basisvalues[r];
1748
}// end loop over 'r'
1749
break;
1750
}
1751
}
1752
1753
}
1754
1755
/// Evaluate all basis functions at given point in cell
1756
virtual void evaluate_basis_all(double* values,
1757
const double* coordinates,
1758
const ufc::cell& c) const
1759
{
1760
// Helper variable to hold values of a single dof.
1761
double dof_values[3] = {0.0, 0.0, 0.0};
1762
1763
// Loop dofs and call evaluate_basis.
1764
for (unsigned int r = 0; r < 12; r++)
1765
{
1766
evaluate_basis(r, dof_values, coordinates, c);
1767
for (unsigned int s = 0; s < 3; s++)
1768
{
1769
values[r*3 + s] = dof_values[s];
1770
}// end loop over 's'
1771
}// end loop over 'r'
1772
}
1773
1774
/// Evaluate order n derivatives of basis function i at given point in cell
1775
virtual void evaluate_basis_derivatives(unsigned int i,
1776
unsigned int n,
1777
double* values,
1778
const double* coordinates,
1779
const ufc::cell& c) const
1780
{
1781
// Extract vertex coordinates
1782
const double * const * x = c.coordinates;
1783
1784
// Compute Jacobian of affine map from reference cell
1785
const double J_00 = x[1][0] - x[0][0];
1786
const double J_01 = x[2][0] - x[0][0];
1787
const double J_02 = x[3][0] - x[0][0];
1788
const double J_10 = x[1][1] - x[0][1];
1789
const double J_11 = x[2][1] - x[0][1];
1790
const double J_12 = x[3][1] - x[0][1];
1791
const double J_20 = x[1][2] - x[0][2];
1792
const double J_21 = x[2][2] - x[0][2];
1793
const double J_22 = x[3][2] - x[0][2];
1794
1795
// Compute sub determinants
1796
const double d_00 = J_11*J_22 - J_12*J_21;
1797
const double d_01 = J_12*J_20 - J_10*J_22;
1798
const double d_02 = J_10*J_21 - J_11*J_20;
1799
const double d_10 = J_02*J_21 - J_01*J_22;
1800
const double d_11 = J_00*J_22 - J_02*J_20;
1801
const double d_12 = J_01*J_20 - J_00*J_21;
1802
const double d_20 = J_01*J_12 - J_02*J_11;
1803
const double d_21 = J_02*J_10 - J_00*J_12;
1804
const double d_22 = J_00*J_11 - J_01*J_10;
1805
1806
// Compute determinant of Jacobian
1807
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
1808
1809
// Compute inverse of Jacobian
1810
const double K_00 = d_00 / detJ;
1811
const double K_01 = d_10 / detJ;
1812
const double K_02 = d_20 / detJ;
1813
const double K_10 = d_01 / detJ;
1814
const double K_11 = d_11 / detJ;
1815
const double K_12 = d_21 / detJ;
1816
const double K_20 = d_02 / detJ;
1817
const double K_21 = d_12 / detJ;
1818
const double K_22 = d_22 / detJ;
1819
1820
// Compute constants
1821
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
1822
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
1823
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
1824
1825
// Get coordinates and map to the reference (FIAT) element
1826
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
1827
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
1828
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
1829
1830
1831
// Compute number of derivatives.
1832
unsigned int num_derivatives = 1;
1833
for (unsigned int r = 0; r < n; r++)
1834
{
1835
num_derivatives *= 3;
1836
}// end loop over 'r'
1837
1838
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1839
unsigned int **combinations = new unsigned int *[num_derivatives];
1840
for (unsigned int row = 0; row < num_derivatives; row++)
1841
{
1842
combinations[row] = new unsigned int [n];
1843
for (unsigned int col = 0; col < n; col++)
1844
combinations[row][col] = 0;
1845
}
1846
1847
// Generate combinations of derivatives
1848
for (unsigned int row = 1; row < num_derivatives; row++)
1849
{
1850
for (unsigned int num = 0; num < row; num++)
1851
{
1852
for (unsigned int col = n-1; col+1 > 0; col--)
1853
{
1854
if (combinations[row][col] + 1 > 2)
1855
combinations[row][col] = 0;
1856
else
1857
{
1858
combinations[row][col] += 1;
1859
break;
1860
}
1861
}
1862
}
1863
}
1864
1865
// Compute inverse of Jacobian
1866
const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}};
1867
1868
// Declare transformation matrix
1869
// Declare pointer to two dimensional array and initialise
1870
double **transform = new double *[num_derivatives];
1871
1872
for (unsigned int j = 0; j < num_derivatives; j++)
1873
{
1874
transform[j] = new double [num_derivatives];
1875
for (unsigned int k = 0; k < num_derivatives; k++)
1876
transform[j][k] = 1;
1877
}
1878
1879
// Construct transformation matrix
1880
for (unsigned int row = 0; row < num_derivatives; row++)
1881
{
1882
for (unsigned int col = 0; col < num_derivatives; col++)
1883
{
1884
for (unsigned int k = 0; k < n; k++)
1885
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1886
}
1887
}
1888
1889
// Reset values. Assuming that values is always an array.
1890
for (unsigned int r = 0; r < 3*num_derivatives; r++)
1891
{
1892
values[r] = 0.0;
1893
}// end loop over 'r'
1894
1895
switch (i)
1896
{
1897
case 0:
1898
{
1899
1900
// Array of basisvalues.
1901
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1902
1903
// Declare helper variables.
1904
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1905
1906
// Compute basisvalues.
1907
basisvalues[0] = 1.0;
1908
basisvalues[1] = tmp0;
1909
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1910
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1911
basisvalues[0] *= std::sqrt(0.75);
1912
basisvalues[3] *= std::sqrt(1.25);
1913
basisvalues[2] *= std::sqrt(2.5);
1914
basisvalues[1] *= std::sqrt(7.5);
1915
1916
// Table(s) of coefficients.
1917
static const double coefficients0[4] = \
1918
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1919
1920
// Tables of derivatives of the polynomial base (transpose).
1921
static const double dmats0[4][4] = \
1922
{{0.0, 0.0, 0.0, 0.0},
1923
{6.3245553, 0.0, 0.0, 0.0},
1924
{0.0, 0.0, 0.0, 0.0},
1925
{0.0, 0.0, 0.0, 0.0}};
1926
1927
static const double dmats1[4][4] = \
1928
{{0.0, 0.0, 0.0, 0.0},
1929
{3.1622777, 0.0, 0.0, 0.0},
1930
{5.4772256, 0.0, 0.0, 0.0},
1931
{0.0, 0.0, 0.0, 0.0}};
1932
1933
static const double dmats2[4][4] = \
1934
{{0.0, 0.0, 0.0, 0.0},
1935
{3.1622777, 0.0, 0.0, 0.0},
1936
{1.8257419, 0.0, 0.0, 0.0},
1937
{5.1639778, 0.0, 0.0, 0.0}};
1938
1939
// Compute reference derivatives.
1940
// Declare pointer to array of derivatives on FIAT element.
1941
double *derivatives = new double[num_derivatives];
1942
for (unsigned int r = 0; r < num_derivatives; r++)
1943
{
1944
derivatives[r] = 0.0;
1945
}// end loop over 'r'
1946
1947
// Declare derivative matrix (of polynomial basis).
1948
double dmats[4][4] = \
1949
{{1.0, 0.0, 0.0, 0.0},
1950
{0.0, 1.0, 0.0, 0.0},
1951
{0.0, 0.0, 1.0, 0.0},
1952
{0.0, 0.0, 0.0, 1.0}};
1953
1954
// Declare (auxiliary) derivative matrix (of polynomial basis).
1955
double dmats_old[4][4] = \
1956
{{1.0, 0.0, 0.0, 0.0},
1957
{0.0, 1.0, 0.0, 0.0},
1958
{0.0, 0.0, 1.0, 0.0},
1959
{0.0, 0.0, 0.0, 1.0}};
1960
1961
// Loop possible derivatives.
1962
for (unsigned int r = 0; r < num_derivatives; r++)
1963
{
1964
// Resetting dmats values to compute next derivative.
1965
for (unsigned int t = 0; t < 4; t++)
1966
{
1967
for (unsigned int u = 0; u < 4; u++)
1968
{
1969
dmats[t][u] = 0.0;
1970
if (t == u)
1971
{
1972
dmats[t][u] = 1.0;
1973
}
1974
1975
}// end loop over 'u'
1976
}// end loop over 't'
1977
1978
// Looping derivative order to generate dmats.
1979
for (unsigned int s = 0; s < n; s++)
1980
{
1981
// Updating dmats_old with new values and resetting dmats.
1982
for (unsigned int t = 0; t < 4; t++)
1983
{
1984
for (unsigned int u = 0; u < 4; u++)
1985
{
1986
dmats_old[t][u] = dmats[t][u];
1987
dmats[t][u] = 0.0;
1988
}// end loop over 'u'
1989
}// end loop over 't'
1990
1991
// Update dmats using an inner product.
1992
if (combinations[r][s] == 0)
1993
{
1994
for (unsigned int t = 0; t < 4; t++)
1995
{
1996
for (unsigned int u = 0; u < 4; u++)
1997
{
1998
for (unsigned int tu = 0; tu < 4; tu++)
1999
{
2000
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2001
}// end loop over 'tu'
2002
}// end loop over 'u'
2003
}// end loop over 't'
2004
}
2005
2006
if (combinations[r][s] == 1)
2007
{
2008
for (unsigned int t = 0; t < 4; t++)
2009
{
2010
for (unsigned int u = 0; u < 4; u++)
2011
{
2012
for (unsigned int tu = 0; tu < 4; tu++)
2013
{
2014
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2015
}// end loop over 'tu'
2016
}// end loop over 'u'
2017
}// end loop over 't'
2018
}
2019
2020
if (combinations[r][s] == 2)
2021
{
2022
for (unsigned int t = 0; t < 4; t++)
2023
{
2024
for (unsigned int u = 0; u < 4; u++)
2025
{
2026
for (unsigned int tu = 0; tu < 4; tu++)
2027
{
2028
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2029
}// end loop over 'tu'
2030
}// end loop over 'u'
2031
}// end loop over 't'
2032
}
2033
2034
}// end loop over 's'
2035
for (unsigned int s = 0; s < 4; s++)
2036
{
2037
for (unsigned int t = 0; t < 4; t++)
2038
{
2039
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2040
}// end loop over 't'
2041
}// end loop over 's'
2042
}// end loop over 'r'
2043
2044
// Transform derivatives back to physical element
2045
for (unsigned int r = 0; r < num_derivatives; r++)
2046
{
2047
for (unsigned int s = 0; s < num_derivatives; s++)
2048
{
2049
values[r] += transform[r][s]*derivatives[s];
2050
}// end loop over 's'
2051
}// end loop over 'r'
2052
2053
// Delete pointer to array of derivatives on FIAT element
2054
delete [] derivatives;
2055
2056
// Delete pointer to array of combinations of derivatives and transform
2057
for (unsigned int r = 0; r < num_derivatives; r++)
2058
{
2059
delete [] combinations[r];
2060
}// end loop over 'r'
2061
delete [] combinations;
2062
for (unsigned int r = 0; r < num_derivatives; r++)
2063
{
2064
delete [] transform[r];
2065
}// end loop over 'r'
2066
delete [] transform;
2067
break;
2068
}
2069
case 1:
2070
{
2071
2072
// Array of basisvalues.
2073
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2074
2075
// Declare helper variables.
2076
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2077
2078
// Compute basisvalues.
2079
basisvalues[0] = 1.0;
2080
basisvalues[1] = tmp0;
2081
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2082
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2083
basisvalues[0] *= std::sqrt(0.75);
2084
basisvalues[3] *= std::sqrt(1.25);
2085
basisvalues[2] *= std::sqrt(2.5);
2086
basisvalues[1] *= std::sqrt(7.5);
2087
2088
// Table(s) of coefficients.
2089
static const double coefficients0[4] = \
2090
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
2091
2092
// Tables of derivatives of the polynomial base (transpose).
2093
static const double dmats0[4][4] = \
2094
{{0.0, 0.0, 0.0, 0.0},
2095
{6.3245553, 0.0, 0.0, 0.0},
2096
{0.0, 0.0, 0.0, 0.0},
2097
{0.0, 0.0, 0.0, 0.0}};
2098
2099
static const double dmats1[4][4] = \
2100
{{0.0, 0.0, 0.0, 0.0},
2101
{3.1622777, 0.0, 0.0, 0.0},
2102
{5.4772256, 0.0, 0.0, 0.0},
2103
{0.0, 0.0, 0.0, 0.0}};
2104
2105
static const double dmats2[4][4] = \
2106
{{0.0, 0.0, 0.0, 0.0},
2107
{3.1622777, 0.0, 0.0, 0.0},
2108
{1.8257419, 0.0, 0.0, 0.0},
2109
{5.1639778, 0.0, 0.0, 0.0}};
2110
2111
// Compute reference derivatives.
2112
// Declare pointer to array of derivatives on FIAT element.
2113
double *derivatives = new double[num_derivatives];
2114
for (unsigned int r = 0; r < num_derivatives; r++)
2115
{
2116
derivatives[r] = 0.0;
2117
}// end loop over 'r'
2118
2119
// Declare derivative matrix (of polynomial basis).
2120
double dmats[4][4] = \
2121
{{1.0, 0.0, 0.0, 0.0},
2122
{0.0, 1.0, 0.0, 0.0},
2123
{0.0, 0.0, 1.0, 0.0},
2124
{0.0, 0.0, 0.0, 1.0}};
2125
2126
// Declare (auxiliary) derivative matrix (of polynomial basis).
2127
double dmats_old[4][4] = \
2128
{{1.0, 0.0, 0.0, 0.0},
2129
{0.0, 1.0, 0.0, 0.0},
2130
{0.0, 0.0, 1.0, 0.0},
2131
{0.0, 0.0, 0.0, 1.0}};
2132
2133
// Loop possible derivatives.
2134
for (unsigned int r = 0; r < num_derivatives; r++)
2135
{
2136
// Resetting dmats values to compute next derivative.
2137
for (unsigned int t = 0; t < 4; t++)
2138
{
2139
for (unsigned int u = 0; u < 4; u++)
2140
{
2141
dmats[t][u] = 0.0;
2142
if (t == u)
2143
{
2144
dmats[t][u] = 1.0;
2145
}
2146
2147
}// end loop over 'u'
2148
}// end loop over 't'
2149
2150
// Looping derivative order to generate dmats.
2151
for (unsigned int s = 0; s < n; s++)
2152
{
2153
// Updating dmats_old with new values and resetting dmats.
2154
for (unsigned int t = 0; t < 4; t++)
2155
{
2156
for (unsigned int u = 0; u < 4; u++)
2157
{
2158
dmats_old[t][u] = dmats[t][u];
2159
dmats[t][u] = 0.0;
2160
}// end loop over 'u'
2161
}// end loop over 't'
2162
2163
// Update dmats using an inner product.
2164
if (combinations[r][s] == 0)
2165
{
2166
for (unsigned int t = 0; t < 4; t++)
2167
{
2168
for (unsigned int u = 0; u < 4; u++)
2169
{
2170
for (unsigned int tu = 0; tu < 4; tu++)
2171
{
2172
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2173
}// end loop over 'tu'
2174
}// end loop over 'u'
2175
}// end loop over 't'
2176
}
2177
2178
if (combinations[r][s] == 1)
2179
{
2180
for (unsigned int t = 0; t < 4; t++)
2181
{
2182
for (unsigned int u = 0; u < 4; u++)
2183
{
2184
for (unsigned int tu = 0; tu < 4; tu++)
2185
{
2186
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2187
}// end loop over 'tu'
2188
}// end loop over 'u'
2189
}// end loop over 't'
2190
}
2191
2192
if (combinations[r][s] == 2)
2193
{
2194
for (unsigned int t = 0; t < 4; t++)
2195
{
2196
for (unsigned int u = 0; u < 4; u++)
2197
{
2198
for (unsigned int tu = 0; tu < 4; tu++)
2199
{
2200
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2201
}// end loop over 'tu'
2202
}// end loop over 'u'
2203
}// end loop over 't'
2204
}
2205
2206
}// end loop over 's'
2207
for (unsigned int s = 0; s < 4; s++)
2208
{
2209
for (unsigned int t = 0; t < 4; t++)
2210
{
2211
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2212
}// end loop over 't'
2213
}// end loop over 's'
2214
}// end loop over 'r'
2215
2216
// Transform derivatives back to physical element
2217
for (unsigned int r = 0; r < num_derivatives; r++)
2218
{
2219
for (unsigned int s = 0; s < num_derivatives; s++)
2220
{
2221
values[r] += transform[r][s]*derivatives[s];
2222
}// end loop over 's'
2223
}// end loop over 'r'
2224
2225
// Delete pointer to array of derivatives on FIAT element
2226
delete [] derivatives;
2227
2228
// Delete pointer to array of combinations of derivatives and transform
2229
for (unsigned int r = 0; r < num_derivatives; r++)
2230
{
2231
delete [] combinations[r];
2232
}// end loop over 'r'
2233
delete [] combinations;
2234
for (unsigned int r = 0; r < num_derivatives; r++)
2235
{
2236
delete [] transform[r];
2237
}// end loop over 'r'
2238
delete [] transform;
2239
break;
2240
}
2241
case 2:
2242
{
2243
2244
// Array of basisvalues.
2245
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2246
2247
// Declare helper variables.
2248
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2249
2250
// Compute basisvalues.
2251
basisvalues[0] = 1.0;
2252
basisvalues[1] = tmp0;
2253
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2254
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2255
basisvalues[0] *= std::sqrt(0.75);
2256
basisvalues[3] *= std::sqrt(1.25);
2257
basisvalues[2] *= std::sqrt(2.5);
2258
basisvalues[1] *= std::sqrt(7.5);
2259
2260
// Table(s) of coefficients.
2261
static const double coefficients0[4] = \
2262
{0.28867513, 0.0, 0.21081851, -0.074535599};
2263
2264
// Tables of derivatives of the polynomial base (transpose).
2265
static const double dmats0[4][4] = \
2266
{{0.0, 0.0, 0.0, 0.0},
2267
{6.3245553, 0.0, 0.0, 0.0},
2268
{0.0, 0.0, 0.0, 0.0},
2269
{0.0, 0.0, 0.0, 0.0}};
2270
2271
static const double dmats1[4][4] = \
2272
{{0.0, 0.0, 0.0, 0.0},
2273
{3.1622777, 0.0, 0.0, 0.0},
2274
{5.4772256, 0.0, 0.0, 0.0},
2275
{0.0, 0.0, 0.0, 0.0}};
2276
2277
static const double dmats2[4][4] = \
2278
{{0.0, 0.0, 0.0, 0.0},
2279
{3.1622777, 0.0, 0.0, 0.0},
2280
{1.8257419, 0.0, 0.0, 0.0},
2281
{5.1639778, 0.0, 0.0, 0.0}};
2282
2283
// Compute reference derivatives.
2284
// Declare pointer to array of derivatives on FIAT element.
2285
double *derivatives = new double[num_derivatives];
2286
for (unsigned int r = 0; r < num_derivatives; r++)
2287
{
2288
derivatives[r] = 0.0;
2289
}// end loop over 'r'
2290
2291
// Declare derivative matrix (of polynomial basis).
2292
double dmats[4][4] = \
2293
{{1.0, 0.0, 0.0, 0.0},
2294
{0.0, 1.0, 0.0, 0.0},
2295
{0.0, 0.0, 1.0, 0.0},
2296
{0.0, 0.0, 0.0, 1.0}};
2297
2298
// Declare (auxiliary) derivative matrix (of polynomial basis).
2299
double dmats_old[4][4] = \
2300
{{1.0, 0.0, 0.0, 0.0},
2301
{0.0, 1.0, 0.0, 0.0},
2302
{0.0, 0.0, 1.0, 0.0},
2303
{0.0, 0.0, 0.0, 1.0}};
2304
2305
// Loop possible derivatives.
2306
for (unsigned int r = 0; r < num_derivatives; r++)
2307
{
2308
// Resetting dmats values to compute next derivative.
2309
for (unsigned int t = 0; t < 4; t++)
2310
{
2311
for (unsigned int u = 0; u < 4; u++)
2312
{
2313
dmats[t][u] = 0.0;
2314
if (t == u)
2315
{
2316
dmats[t][u] = 1.0;
2317
}
2318
2319
}// end loop over 'u'
2320
}// end loop over 't'
2321
2322
// Looping derivative order to generate dmats.
2323
for (unsigned int s = 0; s < n; s++)
2324
{
2325
// Updating dmats_old with new values and resetting dmats.
2326
for (unsigned int t = 0; t < 4; t++)
2327
{
2328
for (unsigned int u = 0; u < 4; u++)
2329
{
2330
dmats_old[t][u] = dmats[t][u];
2331
dmats[t][u] = 0.0;
2332
}// end loop over 'u'
2333
}// end loop over 't'
2334
2335
// Update dmats using an inner product.
2336
if (combinations[r][s] == 0)
2337
{
2338
for (unsigned int t = 0; t < 4; t++)
2339
{
2340
for (unsigned int u = 0; u < 4; u++)
2341
{
2342
for (unsigned int tu = 0; tu < 4; tu++)
2343
{
2344
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2345
}// end loop over 'tu'
2346
}// end loop over 'u'
2347
}// end loop over 't'
2348
}
2349
2350
if (combinations[r][s] == 1)
2351
{
2352
for (unsigned int t = 0; t < 4; t++)
2353
{
2354
for (unsigned int u = 0; u < 4; u++)
2355
{
2356
for (unsigned int tu = 0; tu < 4; tu++)
2357
{
2358
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2359
}// end loop over 'tu'
2360
}// end loop over 'u'
2361
}// end loop over 't'
2362
}
2363
2364
if (combinations[r][s] == 2)
2365
{
2366
for (unsigned int t = 0; t < 4; t++)
2367
{
2368
for (unsigned int u = 0; u < 4; u++)
2369
{
2370
for (unsigned int tu = 0; tu < 4; tu++)
2371
{
2372
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2373
}// end loop over 'tu'
2374
}// end loop over 'u'
2375
}// end loop over 't'
2376
}
2377
2378
}// end loop over 's'
2379
for (unsigned int s = 0; s < 4; s++)
2380
{
2381
for (unsigned int t = 0; t < 4; t++)
2382
{
2383
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2384
}// end loop over 't'
2385
}// end loop over 's'
2386
}// end loop over 'r'
2387
2388
// Transform derivatives back to physical element
2389
for (unsigned int r = 0; r < num_derivatives; r++)
2390
{
2391
for (unsigned int s = 0; s < num_derivatives; s++)
2392
{
2393
values[r] += transform[r][s]*derivatives[s];
2394
}// end loop over 's'
2395
}// end loop over 'r'
2396
2397
// Delete pointer to array of derivatives on FIAT element
2398
delete [] derivatives;
2399
2400
// Delete pointer to array of combinations of derivatives and transform
2401
for (unsigned int r = 0; r < num_derivatives; r++)
2402
{
2403
delete [] combinations[r];
2404
}// end loop over 'r'
2405
delete [] combinations;
2406
for (unsigned int r = 0; r < num_derivatives; r++)
2407
{
2408
delete [] transform[r];
2409
}// end loop over 'r'
2410
delete [] transform;
2411
break;
2412
}
2413
case 3:
2414
{
2415
2416
// Array of basisvalues.
2417
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2418
2419
// Declare helper variables.
2420
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2421
2422
// Compute basisvalues.
2423
basisvalues[0] = 1.0;
2424
basisvalues[1] = tmp0;
2425
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2426
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2427
basisvalues[0] *= std::sqrt(0.75);
2428
basisvalues[3] *= std::sqrt(1.25);
2429
basisvalues[2] *= std::sqrt(2.5);
2430
basisvalues[1] *= std::sqrt(7.5);
2431
2432
// Table(s) of coefficients.
2433
static const double coefficients0[4] = \
2434
{0.28867513, 0.0, 0.0, 0.2236068};
2435
2436
// Tables of derivatives of the polynomial base (transpose).
2437
static const double dmats0[4][4] = \
2438
{{0.0, 0.0, 0.0, 0.0},
2439
{6.3245553, 0.0, 0.0, 0.0},
2440
{0.0, 0.0, 0.0, 0.0},
2441
{0.0, 0.0, 0.0, 0.0}};
2442
2443
static const double dmats1[4][4] = \
2444
{{0.0, 0.0, 0.0, 0.0},
2445
{3.1622777, 0.0, 0.0, 0.0},
2446
{5.4772256, 0.0, 0.0, 0.0},
2447
{0.0, 0.0, 0.0, 0.0}};
2448
2449
static const double dmats2[4][4] = \
2450
{{0.0, 0.0, 0.0, 0.0},
2451
{3.1622777, 0.0, 0.0, 0.0},
2452
{1.8257419, 0.0, 0.0, 0.0},
2453
{5.1639778, 0.0, 0.0, 0.0}};
2454
2455
// Compute reference derivatives.
2456
// Declare pointer to array of derivatives on FIAT element.
2457
double *derivatives = new double[num_derivatives];
2458
for (unsigned int r = 0; r < num_derivatives; r++)
2459
{
2460
derivatives[r] = 0.0;
2461
}// end loop over 'r'
2462
2463
// Declare derivative matrix (of polynomial basis).
2464
double dmats[4][4] = \
2465
{{1.0, 0.0, 0.0, 0.0},
2466
{0.0, 1.0, 0.0, 0.0},
2467
{0.0, 0.0, 1.0, 0.0},
2468
{0.0, 0.0, 0.0, 1.0}};
2469
2470
// Declare (auxiliary) derivative matrix (of polynomial basis).
2471
double dmats_old[4][4] = \
2472
{{1.0, 0.0, 0.0, 0.0},
2473
{0.0, 1.0, 0.0, 0.0},
2474
{0.0, 0.0, 1.0, 0.0},
2475
{0.0, 0.0, 0.0, 1.0}};
2476
2477
// Loop possible derivatives.
2478
for (unsigned int r = 0; r < num_derivatives; r++)
2479
{
2480
// Resetting dmats values to compute next derivative.
2481
for (unsigned int t = 0; t < 4; t++)
2482
{
2483
for (unsigned int u = 0; u < 4; u++)
2484
{
2485
dmats[t][u] = 0.0;
2486
if (t == u)
2487
{
2488
dmats[t][u] = 1.0;
2489
}
2490
2491
}// end loop over 'u'
2492
}// end loop over 't'
2493
2494
// Looping derivative order to generate dmats.
2495
for (unsigned int s = 0; s < n; s++)
2496
{
2497
// Updating dmats_old with new values and resetting dmats.
2498
for (unsigned int t = 0; t < 4; t++)
2499
{
2500
for (unsigned int u = 0; u < 4; u++)
2501
{
2502
dmats_old[t][u] = dmats[t][u];
2503
dmats[t][u] = 0.0;
2504
}// end loop over 'u'
2505
}// end loop over 't'
2506
2507
// Update dmats using an inner product.
2508
if (combinations[r][s] == 0)
2509
{
2510
for (unsigned int t = 0; t < 4; t++)
2511
{
2512
for (unsigned int u = 0; u < 4; u++)
2513
{
2514
for (unsigned int tu = 0; tu < 4; tu++)
2515
{
2516
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2517
}// end loop over 'tu'
2518
}// end loop over 'u'
2519
}// end loop over 't'
2520
}
2521
2522
if (combinations[r][s] == 1)
2523
{
2524
for (unsigned int t = 0; t < 4; t++)
2525
{
2526
for (unsigned int u = 0; u < 4; u++)
2527
{
2528
for (unsigned int tu = 0; tu < 4; tu++)
2529
{
2530
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2531
}// end loop over 'tu'
2532
}// end loop over 'u'
2533
}// end loop over 't'
2534
}
2535
2536
if (combinations[r][s] == 2)
2537
{
2538
for (unsigned int t = 0; t < 4; t++)
2539
{
2540
for (unsigned int u = 0; u < 4; u++)
2541
{
2542
for (unsigned int tu = 0; tu < 4; tu++)
2543
{
2544
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2545
}// end loop over 'tu'
2546
}// end loop over 'u'
2547
}// end loop over 't'
2548
}
2549
2550
}// end loop over 's'
2551
for (unsigned int s = 0; s < 4; s++)
2552
{
2553
for (unsigned int t = 0; t < 4; t++)
2554
{
2555
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2556
}// end loop over 't'
2557
}// end loop over 's'
2558
}// end loop over 'r'
2559
2560
// Transform derivatives back to physical element
2561
for (unsigned int r = 0; r < num_derivatives; r++)
2562
{
2563
for (unsigned int s = 0; s < num_derivatives; s++)
2564
{
2565
values[r] += transform[r][s]*derivatives[s];
2566
}// end loop over 's'
2567
}// end loop over 'r'
2568
2569
// Delete pointer to array of derivatives on FIAT element
2570
delete [] derivatives;
2571
2572
// Delete pointer to array of combinations of derivatives and transform
2573
for (unsigned int r = 0; r < num_derivatives; r++)
2574
{
2575
delete [] combinations[r];
2576
}// end loop over 'r'
2577
delete [] combinations;
2578
for (unsigned int r = 0; r < num_derivatives; r++)
2579
{
2580
delete [] transform[r];
2581
}// end loop over 'r'
2582
delete [] transform;
2583
break;
2584
}
2585
case 4:
2586
{
2587
2588
// Array of basisvalues.
2589
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2590
2591
// Declare helper variables.
2592
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2593
2594
// Compute basisvalues.
2595
basisvalues[0] = 1.0;
2596
basisvalues[1] = tmp0;
2597
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2598
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2599
basisvalues[0] *= std::sqrt(0.75);
2600
basisvalues[3] *= std::sqrt(1.25);
2601
basisvalues[2] *= std::sqrt(2.5);
2602
basisvalues[1] *= std::sqrt(7.5);
2603
2604
// Table(s) of coefficients.
2605
static const double coefficients0[4] = \
2606
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
2607
2608
// Tables of derivatives of the polynomial base (transpose).
2609
static const double dmats0[4][4] = \
2610
{{0.0, 0.0, 0.0, 0.0},
2611
{6.3245553, 0.0, 0.0, 0.0},
2612
{0.0, 0.0, 0.0, 0.0},
2613
{0.0, 0.0, 0.0, 0.0}};
2614
2615
static const double dmats1[4][4] = \
2616
{{0.0, 0.0, 0.0, 0.0},
2617
{3.1622777, 0.0, 0.0, 0.0},
2618
{5.4772256, 0.0, 0.0, 0.0},
2619
{0.0, 0.0, 0.0, 0.0}};
2620
2621
static const double dmats2[4][4] = \
2622
{{0.0, 0.0, 0.0, 0.0},
2623
{3.1622777, 0.0, 0.0, 0.0},
2624
{1.8257419, 0.0, 0.0, 0.0},
2625
{5.1639778, 0.0, 0.0, 0.0}};
2626
2627
// Compute reference derivatives.
2628
// Declare pointer to array of derivatives on FIAT element.
2629
double *derivatives = new double[num_derivatives];
2630
for (unsigned int r = 0; r < num_derivatives; r++)
2631
{
2632
derivatives[r] = 0.0;
2633
}// end loop over 'r'
2634
2635
// Declare derivative matrix (of polynomial basis).
2636
double dmats[4][4] = \
2637
{{1.0, 0.0, 0.0, 0.0},
2638
{0.0, 1.0, 0.0, 0.0},
2639
{0.0, 0.0, 1.0, 0.0},
2640
{0.0, 0.0, 0.0, 1.0}};
2641
2642
// Declare (auxiliary) derivative matrix (of polynomial basis).
2643
double dmats_old[4][4] = \
2644
{{1.0, 0.0, 0.0, 0.0},
2645
{0.0, 1.0, 0.0, 0.0},
2646
{0.0, 0.0, 1.0, 0.0},
2647
{0.0, 0.0, 0.0, 1.0}};
2648
2649
// Loop possible derivatives.
2650
for (unsigned int r = 0; r < num_derivatives; r++)
2651
{
2652
// Resetting dmats values to compute next derivative.
2653
for (unsigned int t = 0; t < 4; t++)
2654
{
2655
for (unsigned int u = 0; u < 4; u++)
2656
{
2657
dmats[t][u] = 0.0;
2658
if (t == u)
2659
{
2660
dmats[t][u] = 1.0;
2661
}
2662
2663
}// end loop over 'u'
2664
}// end loop over 't'
2665
2666
// Looping derivative order to generate dmats.
2667
for (unsigned int s = 0; s < n; s++)
2668
{
2669
// Updating dmats_old with new values and resetting dmats.
2670
for (unsigned int t = 0; t < 4; t++)
2671
{
2672
for (unsigned int u = 0; u < 4; u++)
2673
{
2674
dmats_old[t][u] = dmats[t][u];
2675
dmats[t][u] = 0.0;
2676
}// end loop over 'u'
2677
}// end loop over 't'
2678
2679
// Update dmats using an inner product.
2680
if (combinations[r][s] == 0)
2681
{
2682
for (unsigned int t = 0; t < 4; t++)
2683
{
2684
for (unsigned int u = 0; u < 4; u++)
2685
{
2686
for (unsigned int tu = 0; tu < 4; tu++)
2687
{
2688
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2689
}// end loop over 'tu'
2690
}// end loop over 'u'
2691
}// end loop over 't'
2692
}
2693
2694
if (combinations[r][s] == 1)
2695
{
2696
for (unsigned int t = 0; t < 4; t++)
2697
{
2698
for (unsigned int u = 0; u < 4; u++)
2699
{
2700
for (unsigned int tu = 0; tu < 4; tu++)
2701
{
2702
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2703
}// end loop over 'tu'
2704
}// end loop over 'u'
2705
}// end loop over 't'
2706
}
2707
2708
if (combinations[r][s] == 2)
2709
{
2710
for (unsigned int t = 0; t < 4; t++)
2711
{
2712
for (unsigned int u = 0; u < 4; u++)
2713
{
2714
for (unsigned int tu = 0; tu < 4; tu++)
2715
{
2716
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2717
}// end loop over 'tu'
2718
}// end loop over 'u'
2719
}// end loop over 't'
2720
}
2721
2722
}// end loop over 's'
2723
for (unsigned int s = 0; s < 4; s++)
2724
{
2725
for (unsigned int t = 0; t < 4; t++)
2726
{
2727
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2728
}// end loop over 't'
2729
}// end loop over 's'
2730
}// end loop over 'r'
2731
2732
// Transform derivatives back to physical element
2733
for (unsigned int r = 0; r < num_derivatives; r++)
2734
{
2735
for (unsigned int s = 0; s < num_derivatives; s++)
2736
{
2737
values[num_derivatives + r] += transform[r][s]*derivatives[s];
2738
}// end loop over 's'
2739
}// end loop over 'r'
2740
2741
// Delete pointer to array of derivatives on FIAT element
2742
delete [] derivatives;
2743
2744
// Delete pointer to array of combinations of derivatives and transform
2745
for (unsigned int r = 0; r < num_derivatives; r++)
2746
{
2747
delete [] combinations[r];
2748
}// end loop over 'r'
2749
delete [] combinations;
2750
for (unsigned int r = 0; r < num_derivatives; r++)
2751
{
2752
delete [] transform[r];
2753
}// end loop over 'r'
2754
delete [] transform;
2755
break;
2756
}
2757
case 5:
2758
{
2759
2760
// Array of basisvalues.
2761
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2762
2763
// Declare helper variables.
2764
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2765
2766
// Compute basisvalues.
2767
basisvalues[0] = 1.0;
2768
basisvalues[1] = tmp0;
2769
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2770
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2771
basisvalues[0] *= std::sqrt(0.75);
2772
basisvalues[3] *= std::sqrt(1.25);
2773
basisvalues[2] *= std::sqrt(2.5);
2774
basisvalues[1] *= std::sqrt(7.5);
2775
2776
// Table(s) of coefficients.
2777
static const double coefficients0[4] = \
2778
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
2779
2780
// Tables of derivatives of the polynomial base (transpose).
2781
static const double dmats0[4][4] = \
2782
{{0.0, 0.0, 0.0, 0.0},
2783
{6.3245553, 0.0, 0.0, 0.0},
2784
{0.0, 0.0, 0.0, 0.0},
2785
{0.0, 0.0, 0.0, 0.0}};
2786
2787
static const double dmats1[4][4] = \
2788
{{0.0, 0.0, 0.0, 0.0},
2789
{3.1622777, 0.0, 0.0, 0.0},
2790
{5.4772256, 0.0, 0.0, 0.0},
2791
{0.0, 0.0, 0.0, 0.0}};
2792
2793
static const double dmats2[4][4] = \
2794
{{0.0, 0.0, 0.0, 0.0},
2795
{3.1622777, 0.0, 0.0, 0.0},
2796
{1.8257419, 0.0, 0.0, 0.0},
2797
{5.1639778, 0.0, 0.0, 0.0}};
2798
2799
// Compute reference derivatives.
2800
// Declare pointer to array of derivatives on FIAT element.
2801
double *derivatives = new double[num_derivatives];
2802
for (unsigned int r = 0; r < num_derivatives; r++)
2803
{
2804
derivatives[r] = 0.0;
2805
}// end loop over 'r'
2806
2807
// Declare derivative matrix (of polynomial basis).
2808
double dmats[4][4] = \
2809
{{1.0, 0.0, 0.0, 0.0},
2810
{0.0, 1.0, 0.0, 0.0},
2811
{0.0, 0.0, 1.0, 0.0},
2812
{0.0, 0.0, 0.0, 1.0}};
2813
2814
// Declare (auxiliary) derivative matrix (of polynomial basis).
2815
double dmats_old[4][4] = \
2816
{{1.0, 0.0, 0.0, 0.0},
2817
{0.0, 1.0, 0.0, 0.0},
2818
{0.0, 0.0, 1.0, 0.0},
2819
{0.0, 0.0, 0.0, 1.0}};
2820
2821
// Loop possible derivatives.
2822
for (unsigned int r = 0; r < num_derivatives; r++)
2823
{
2824
// Resetting dmats values to compute next derivative.
2825
for (unsigned int t = 0; t < 4; t++)
2826
{
2827
for (unsigned int u = 0; u < 4; u++)
2828
{
2829
dmats[t][u] = 0.0;
2830
if (t == u)
2831
{
2832
dmats[t][u] = 1.0;
2833
}
2834
2835
}// end loop over 'u'
2836
}// end loop over 't'
2837
2838
// Looping derivative order to generate dmats.
2839
for (unsigned int s = 0; s < n; s++)
2840
{
2841
// Updating dmats_old with new values and resetting dmats.
2842
for (unsigned int t = 0; t < 4; t++)
2843
{
2844
for (unsigned int u = 0; u < 4; u++)
2845
{
2846
dmats_old[t][u] = dmats[t][u];
2847
dmats[t][u] = 0.0;
2848
}// end loop over 'u'
2849
}// end loop over 't'
2850
2851
// Update dmats using an inner product.
2852
if (combinations[r][s] == 0)
2853
{
2854
for (unsigned int t = 0; t < 4; t++)
2855
{
2856
for (unsigned int u = 0; u < 4; u++)
2857
{
2858
for (unsigned int tu = 0; tu < 4; tu++)
2859
{
2860
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2861
}// end loop over 'tu'
2862
}// end loop over 'u'
2863
}// end loop over 't'
2864
}
2865
2866
if (combinations[r][s] == 1)
2867
{
2868
for (unsigned int t = 0; t < 4; t++)
2869
{
2870
for (unsigned int u = 0; u < 4; u++)
2871
{
2872
for (unsigned int tu = 0; tu < 4; tu++)
2873
{
2874
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2875
}// end loop over 'tu'
2876
}// end loop over 'u'
2877
}// end loop over 't'
2878
}
2879
2880
if (combinations[r][s] == 2)
2881
{
2882
for (unsigned int t = 0; t < 4; t++)
2883
{
2884
for (unsigned int u = 0; u < 4; u++)
2885
{
2886
for (unsigned int tu = 0; tu < 4; tu++)
2887
{
2888
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2889
}// end loop over 'tu'
2890
}// end loop over 'u'
2891
}// end loop over 't'
2892
}
2893
2894
}// end loop over 's'
2895
for (unsigned int s = 0; s < 4; s++)
2896
{
2897
for (unsigned int t = 0; t < 4; t++)
2898
{
2899
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2900
}// end loop over 't'
2901
}// end loop over 's'
2902
}// end loop over 'r'
2903
2904
// Transform derivatives back to physical element
2905
for (unsigned int r = 0; r < num_derivatives; r++)
2906
{
2907
for (unsigned int s = 0; s < num_derivatives; s++)
2908
{
2909
values[num_derivatives + r] += transform[r][s]*derivatives[s];
2910
}// end loop over 's'
2911
}// end loop over 'r'
2912
2913
// Delete pointer to array of derivatives on FIAT element
2914
delete [] derivatives;
2915
2916
// Delete pointer to array of combinations of derivatives and transform
2917
for (unsigned int r = 0; r < num_derivatives; r++)
2918
{
2919
delete [] combinations[r];
2920
}// end loop over 'r'
2921
delete [] combinations;
2922
for (unsigned int r = 0; r < num_derivatives; r++)
2923
{
2924
delete [] transform[r];
2925
}// end loop over 'r'
2926
delete [] transform;
2927
break;
2928
}
2929
case 6:
2930
{
2931
2932
// Array of basisvalues.
2933
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2934
2935
// Declare helper variables.
2936
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2937
2938
// Compute basisvalues.
2939
basisvalues[0] = 1.0;
2940
basisvalues[1] = tmp0;
2941
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2942
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2943
basisvalues[0] *= std::sqrt(0.75);
2944
basisvalues[3] *= std::sqrt(1.25);
2945
basisvalues[2] *= std::sqrt(2.5);
2946
basisvalues[1] *= std::sqrt(7.5);
2947
2948
// Table(s) of coefficients.
2949
static const double coefficients0[4] = \
2950
{0.28867513, 0.0, 0.21081851, -0.074535599};
2951
2952
// Tables of derivatives of the polynomial base (transpose).
2953
static const double dmats0[4][4] = \
2954
{{0.0, 0.0, 0.0, 0.0},
2955
{6.3245553, 0.0, 0.0, 0.0},
2956
{0.0, 0.0, 0.0, 0.0},
2957
{0.0, 0.0, 0.0, 0.0}};
2958
2959
static const double dmats1[4][4] = \
2960
{{0.0, 0.0, 0.0, 0.0},
2961
{3.1622777, 0.0, 0.0, 0.0},
2962
{5.4772256, 0.0, 0.0, 0.0},
2963
{0.0, 0.0, 0.0, 0.0}};
2964
2965
static const double dmats2[4][4] = \
2966
{{0.0, 0.0, 0.0, 0.0},
2967
{3.1622777, 0.0, 0.0, 0.0},
2968
{1.8257419, 0.0, 0.0, 0.0},
2969
{5.1639778, 0.0, 0.0, 0.0}};
2970
2971
// Compute reference derivatives.
2972
// Declare pointer to array of derivatives on FIAT element.
2973
double *derivatives = new double[num_derivatives];
2974
for (unsigned int r = 0; r < num_derivatives; r++)
2975
{
2976
derivatives[r] = 0.0;
2977
}// end loop over 'r'
2978
2979
// Declare derivative matrix (of polynomial basis).
2980
double dmats[4][4] = \
2981
{{1.0, 0.0, 0.0, 0.0},
2982
{0.0, 1.0, 0.0, 0.0},
2983
{0.0, 0.0, 1.0, 0.0},
2984
{0.0, 0.0, 0.0, 1.0}};
2985
2986
// Declare (auxiliary) derivative matrix (of polynomial basis).
2987
double dmats_old[4][4] = \
2988
{{1.0, 0.0, 0.0, 0.0},
2989
{0.0, 1.0, 0.0, 0.0},
2990
{0.0, 0.0, 1.0, 0.0},
2991
{0.0, 0.0, 0.0, 1.0}};
2992
2993
// Loop possible derivatives.
2994
for (unsigned int r = 0; r < num_derivatives; r++)
2995
{
2996
// Resetting dmats values to compute next derivative.
2997
for (unsigned int t = 0; t < 4; t++)
2998
{
2999
for (unsigned int u = 0; u < 4; u++)
3000
{
3001
dmats[t][u] = 0.0;
3002
if (t == u)
3003
{
3004
dmats[t][u] = 1.0;
3005
}
3006
3007
}// end loop over 'u'
3008
}// end loop over 't'
3009
3010
// Looping derivative order to generate dmats.
3011
for (unsigned int s = 0; s < n; s++)
3012
{
3013
// Updating dmats_old with new values and resetting dmats.
3014
for (unsigned int t = 0; t < 4; t++)
3015
{
3016
for (unsigned int u = 0; u < 4; u++)
3017
{
3018
dmats_old[t][u] = dmats[t][u];
3019
dmats[t][u] = 0.0;
3020
}// end loop over 'u'
3021
}// end loop over 't'
3022
3023
// Update dmats using an inner product.
3024
if (combinations[r][s] == 0)
3025
{
3026
for (unsigned int t = 0; t < 4; t++)
3027
{
3028
for (unsigned int u = 0; u < 4; u++)
3029
{
3030
for (unsigned int tu = 0; tu < 4; tu++)
3031
{
3032
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3033
}// end loop over 'tu'
3034
}// end loop over 'u'
3035
}// end loop over 't'
3036
}
3037
3038
if (combinations[r][s] == 1)
3039
{
3040
for (unsigned int t = 0; t < 4; t++)
3041
{
3042
for (unsigned int u = 0; u < 4; u++)
3043
{
3044
for (unsigned int tu = 0; tu < 4; tu++)
3045
{
3046
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3047
}// end loop over 'tu'
3048
}// end loop over 'u'
3049
}// end loop over 't'
3050
}
3051
3052
if (combinations[r][s] == 2)
3053
{
3054
for (unsigned int t = 0; t < 4; t++)
3055
{
3056
for (unsigned int u = 0; u < 4; u++)
3057
{
3058
for (unsigned int tu = 0; tu < 4; tu++)
3059
{
3060
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3061
}// end loop over 'tu'
3062
}// end loop over 'u'
3063
}// end loop over 't'
3064
}
3065
3066
}// end loop over 's'
3067
for (unsigned int s = 0; s < 4; s++)
3068
{
3069
for (unsigned int t = 0; t < 4; t++)
3070
{
3071
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3072
}// end loop over 't'
3073
}// end loop over 's'
3074
}// end loop over 'r'
3075
3076
// Transform derivatives back to physical element
3077
for (unsigned int r = 0; r < num_derivatives; r++)
3078
{
3079
for (unsigned int s = 0; s < num_derivatives; s++)
3080
{
3081
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3082
}// end loop over 's'
3083
}// end loop over 'r'
3084
3085
// Delete pointer to array of derivatives on FIAT element
3086
delete [] derivatives;
3087
3088
// Delete pointer to array of combinations of derivatives and transform
3089
for (unsigned int r = 0; r < num_derivatives; r++)
3090
{
3091
delete [] combinations[r];
3092
}// end loop over 'r'
3093
delete [] combinations;
3094
for (unsigned int r = 0; r < num_derivatives; r++)
3095
{
3096
delete [] transform[r];
3097
}// end loop over 'r'
3098
delete [] transform;
3099
break;
3100
}
3101
case 7:
3102
{
3103
3104
// Array of basisvalues.
3105
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3106
3107
// Declare helper variables.
3108
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3109
3110
// Compute basisvalues.
3111
basisvalues[0] = 1.0;
3112
basisvalues[1] = tmp0;
3113
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3114
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3115
basisvalues[0] *= std::sqrt(0.75);
3116
basisvalues[3] *= std::sqrt(1.25);
3117
basisvalues[2] *= std::sqrt(2.5);
3118
basisvalues[1] *= std::sqrt(7.5);
3119
3120
// Table(s) of coefficients.
3121
static const double coefficients0[4] = \
3122
{0.28867513, 0.0, 0.0, 0.2236068};
3123
3124
// Tables of derivatives of the polynomial base (transpose).
3125
static const double dmats0[4][4] = \
3126
{{0.0, 0.0, 0.0, 0.0},
3127
{6.3245553, 0.0, 0.0, 0.0},
3128
{0.0, 0.0, 0.0, 0.0},
3129
{0.0, 0.0, 0.0, 0.0}};
3130
3131
static const double dmats1[4][4] = \
3132
{{0.0, 0.0, 0.0, 0.0},
3133
{3.1622777, 0.0, 0.0, 0.0},
3134
{5.4772256, 0.0, 0.0, 0.0},
3135
{0.0, 0.0, 0.0, 0.0}};
3136
3137
static const double dmats2[4][4] = \
3138
{{0.0, 0.0, 0.0, 0.0},
3139
{3.1622777, 0.0, 0.0, 0.0},
3140
{1.8257419, 0.0, 0.0, 0.0},
3141
{5.1639778, 0.0, 0.0, 0.0}};
3142
3143
// Compute reference derivatives.
3144
// Declare pointer to array of derivatives on FIAT element.
3145
double *derivatives = new double[num_derivatives];
3146
for (unsigned int r = 0; r < num_derivatives; r++)
3147
{
3148
derivatives[r] = 0.0;
3149
}// end loop over 'r'
3150
3151
// Declare derivative matrix (of polynomial basis).
3152
double dmats[4][4] = \
3153
{{1.0, 0.0, 0.0, 0.0},
3154
{0.0, 1.0, 0.0, 0.0},
3155
{0.0, 0.0, 1.0, 0.0},
3156
{0.0, 0.0, 0.0, 1.0}};
3157
3158
// Declare (auxiliary) derivative matrix (of polynomial basis).
3159
double dmats_old[4][4] = \
3160
{{1.0, 0.0, 0.0, 0.0},
3161
{0.0, 1.0, 0.0, 0.0},
3162
{0.0, 0.0, 1.0, 0.0},
3163
{0.0, 0.0, 0.0, 1.0}};
3164
3165
// Loop possible derivatives.
3166
for (unsigned int r = 0; r < num_derivatives; r++)
3167
{
3168
// Resetting dmats values to compute next derivative.
3169
for (unsigned int t = 0; t < 4; t++)
3170
{
3171
for (unsigned int u = 0; u < 4; u++)
3172
{
3173
dmats[t][u] = 0.0;
3174
if (t == u)
3175
{
3176
dmats[t][u] = 1.0;
3177
}
3178
3179
}// end loop over 'u'
3180
}// end loop over 't'
3181
3182
// Looping derivative order to generate dmats.
3183
for (unsigned int s = 0; s < n; s++)
3184
{
3185
// Updating dmats_old with new values and resetting dmats.
3186
for (unsigned int t = 0; t < 4; t++)
3187
{
3188
for (unsigned int u = 0; u < 4; u++)
3189
{
3190
dmats_old[t][u] = dmats[t][u];
3191
dmats[t][u] = 0.0;
3192
}// end loop over 'u'
3193
}// end loop over 't'
3194
3195
// Update dmats using an inner product.
3196
if (combinations[r][s] == 0)
3197
{
3198
for (unsigned int t = 0; t < 4; t++)
3199
{
3200
for (unsigned int u = 0; u < 4; u++)
3201
{
3202
for (unsigned int tu = 0; tu < 4; tu++)
3203
{
3204
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3205
}// end loop over 'tu'
3206
}// end loop over 'u'
3207
}// end loop over 't'
3208
}
3209
3210
if (combinations[r][s] == 1)
3211
{
3212
for (unsigned int t = 0; t < 4; t++)
3213
{
3214
for (unsigned int u = 0; u < 4; u++)
3215
{
3216
for (unsigned int tu = 0; tu < 4; tu++)
3217
{
3218
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3219
}// end loop over 'tu'
3220
}// end loop over 'u'
3221
}// end loop over 't'
3222
}
3223
3224
if (combinations[r][s] == 2)
3225
{
3226
for (unsigned int t = 0; t < 4; t++)
3227
{
3228
for (unsigned int u = 0; u < 4; u++)
3229
{
3230
for (unsigned int tu = 0; tu < 4; tu++)
3231
{
3232
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3233
}// end loop over 'tu'
3234
}// end loop over 'u'
3235
}// end loop over 't'
3236
}
3237
3238
}// end loop over 's'
3239
for (unsigned int s = 0; s < 4; s++)
3240
{
3241
for (unsigned int t = 0; t < 4; t++)
3242
{
3243
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3244
}// end loop over 't'
3245
}// end loop over 's'
3246
}// end loop over 'r'
3247
3248
// Transform derivatives back to physical element
3249
for (unsigned int r = 0; r < num_derivatives; r++)
3250
{
3251
for (unsigned int s = 0; s < num_derivatives; s++)
3252
{
3253
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3254
}// end loop over 's'
3255
}// end loop over 'r'
3256
3257
// Delete pointer to array of derivatives on FIAT element
3258
delete [] derivatives;
3259
3260
// Delete pointer to array of combinations of derivatives and transform
3261
for (unsigned int r = 0; r < num_derivatives; r++)
3262
{
3263
delete [] combinations[r];
3264
}// end loop over 'r'
3265
delete [] combinations;
3266
for (unsigned int r = 0; r < num_derivatives; r++)
3267
{
3268
delete [] transform[r];
3269
}// end loop over 'r'
3270
delete [] transform;
3271
break;
3272
}
3273
case 8:
3274
{
3275
3276
// Array of basisvalues.
3277
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3278
3279
// Declare helper variables.
3280
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3281
3282
// Compute basisvalues.
3283
basisvalues[0] = 1.0;
3284
basisvalues[1] = tmp0;
3285
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3286
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3287
basisvalues[0] *= std::sqrt(0.75);
3288
basisvalues[3] *= std::sqrt(1.25);
3289
basisvalues[2] *= std::sqrt(2.5);
3290
basisvalues[1] *= std::sqrt(7.5);
3291
3292
// Table(s) of coefficients.
3293
static const double coefficients0[4] = \
3294
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
3295
3296
// Tables of derivatives of the polynomial base (transpose).
3297
static const double dmats0[4][4] = \
3298
{{0.0, 0.0, 0.0, 0.0},
3299
{6.3245553, 0.0, 0.0, 0.0},
3300
{0.0, 0.0, 0.0, 0.0},
3301
{0.0, 0.0, 0.0, 0.0}};
3302
3303
static const double dmats1[4][4] = \
3304
{{0.0, 0.0, 0.0, 0.0},
3305
{3.1622777, 0.0, 0.0, 0.0},
3306
{5.4772256, 0.0, 0.0, 0.0},
3307
{0.0, 0.0, 0.0, 0.0}};
3308
3309
static const double dmats2[4][4] = \
3310
{{0.0, 0.0, 0.0, 0.0},
3311
{3.1622777, 0.0, 0.0, 0.0},
3312
{1.8257419, 0.0, 0.0, 0.0},
3313
{5.1639778, 0.0, 0.0, 0.0}};
3314
3315
// Compute reference derivatives.
3316
// Declare pointer to array of derivatives on FIAT element.
3317
double *derivatives = new double[num_derivatives];
3318
for (unsigned int r = 0; r < num_derivatives; r++)
3319
{
3320
derivatives[r] = 0.0;
3321
}// end loop over 'r'
3322
3323
// Declare derivative matrix (of polynomial basis).
3324
double dmats[4][4] = \
3325
{{1.0, 0.0, 0.0, 0.0},
3326
{0.0, 1.0, 0.0, 0.0},
3327
{0.0, 0.0, 1.0, 0.0},
3328
{0.0, 0.0, 0.0, 1.0}};
3329
3330
// Declare (auxiliary) derivative matrix (of polynomial basis).
3331
double dmats_old[4][4] = \
3332
{{1.0, 0.0, 0.0, 0.0},
3333
{0.0, 1.0, 0.0, 0.0},
3334
{0.0, 0.0, 1.0, 0.0},
3335
{0.0, 0.0, 0.0, 1.0}};
3336
3337
// Loop possible derivatives.
3338
for (unsigned int r = 0; r < num_derivatives; r++)
3339
{
3340
// Resetting dmats values to compute next derivative.
3341
for (unsigned int t = 0; t < 4; t++)
3342
{
3343
for (unsigned int u = 0; u < 4; u++)
3344
{
3345
dmats[t][u] = 0.0;
3346
if (t == u)
3347
{
3348
dmats[t][u] = 1.0;
3349
}
3350
3351
}// end loop over 'u'
3352
}// end loop over 't'
3353
3354
// Looping derivative order to generate dmats.
3355
for (unsigned int s = 0; s < n; s++)
3356
{
3357
// Updating dmats_old with new values and resetting dmats.
3358
for (unsigned int t = 0; t < 4; t++)
3359
{
3360
for (unsigned int u = 0; u < 4; u++)
3361
{
3362
dmats_old[t][u] = dmats[t][u];
3363
dmats[t][u] = 0.0;
3364
}// end loop over 'u'
3365
}// end loop over 't'
3366
3367
// Update dmats using an inner product.
3368
if (combinations[r][s] == 0)
3369
{
3370
for (unsigned int t = 0; t < 4; t++)
3371
{
3372
for (unsigned int u = 0; u < 4; u++)
3373
{
3374
for (unsigned int tu = 0; tu < 4; tu++)
3375
{
3376
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3377
}// end loop over 'tu'
3378
}// end loop over 'u'
3379
}// end loop over 't'
3380
}
3381
3382
if (combinations[r][s] == 1)
3383
{
3384
for (unsigned int t = 0; t < 4; t++)
3385
{
3386
for (unsigned int u = 0; u < 4; u++)
3387
{
3388
for (unsigned int tu = 0; tu < 4; tu++)
3389
{
3390
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3391
}// end loop over 'tu'
3392
}// end loop over 'u'
3393
}// end loop over 't'
3394
}
3395
3396
if (combinations[r][s] == 2)
3397
{
3398
for (unsigned int t = 0; t < 4; t++)
3399
{
3400
for (unsigned int u = 0; u < 4; u++)
3401
{
3402
for (unsigned int tu = 0; tu < 4; tu++)
3403
{
3404
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3405
}// end loop over 'tu'
3406
}// end loop over 'u'
3407
}// end loop over 't'
3408
}
3409
3410
}// end loop over 's'
3411
for (unsigned int s = 0; s < 4; s++)
3412
{
3413
for (unsigned int t = 0; t < 4; t++)
3414
{
3415
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3416
}// end loop over 't'
3417
}// end loop over 's'
3418
}// end loop over 'r'
3419
3420
// Transform derivatives back to physical element
3421
for (unsigned int r = 0; r < num_derivatives; r++)
3422
{
3423
for (unsigned int s = 0; s < num_derivatives; s++)
3424
{
3425
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
3426
}// end loop over 's'
3427
}// end loop over 'r'
3428
3429
// Delete pointer to array of derivatives on FIAT element
3430
delete [] derivatives;
3431
3432
// Delete pointer to array of combinations of derivatives and transform
3433
for (unsigned int r = 0; r < num_derivatives; r++)
3434
{
3435
delete [] combinations[r];
3436
}// end loop over 'r'
3437
delete [] combinations;
3438
for (unsigned int r = 0; r < num_derivatives; r++)
3439
{
3440
delete [] transform[r];
3441
}// end loop over 'r'
3442
delete [] transform;
3443
break;
3444
}
3445
case 9:
3446
{
3447
3448
// Array of basisvalues.
3449
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3450
3451
// Declare helper variables.
3452
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3453
3454
// Compute basisvalues.
3455
basisvalues[0] = 1.0;
3456
basisvalues[1] = tmp0;
3457
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3458
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3459
basisvalues[0] *= std::sqrt(0.75);
3460
basisvalues[3] *= std::sqrt(1.25);
3461
basisvalues[2] *= std::sqrt(2.5);
3462
basisvalues[1] *= std::sqrt(7.5);
3463
3464
// Table(s) of coefficients.
3465
static const double coefficients0[4] = \
3466
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
3467
3468
// Tables of derivatives of the polynomial base (transpose).
3469
static const double dmats0[4][4] = \
3470
{{0.0, 0.0, 0.0, 0.0},
3471
{6.3245553, 0.0, 0.0, 0.0},
3472
{0.0, 0.0, 0.0, 0.0},
3473
{0.0, 0.0, 0.0, 0.0}};
3474
3475
static const double dmats1[4][4] = \
3476
{{0.0, 0.0, 0.0, 0.0},
3477
{3.1622777, 0.0, 0.0, 0.0},
3478
{5.4772256, 0.0, 0.0, 0.0},
3479
{0.0, 0.0, 0.0, 0.0}};
3480
3481
static const double dmats2[4][4] = \
3482
{{0.0, 0.0, 0.0, 0.0},
3483
{3.1622777, 0.0, 0.0, 0.0},
3484
{1.8257419, 0.0, 0.0, 0.0},
3485
{5.1639778, 0.0, 0.0, 0.0}};
3486
3487
// Compute reference derivatives.
3488
// Declare pointer to array of derivatives on FIAT element.
3489
double *derivatives = new double[num_derivatives];
3490
for (unsigned int r = 0; r < num_derivatives; r++)
3491
{
3492
derivatives[r] = 0.0;
3493
}// end loop over 'r'
3494
3495
// Declare derivative matrix (of polynomial basis).
3496
double dmats[4][4] = \
3497
{{1.0, 0.0, 0.0, 0.0},
3498
{0.0, 1.0, 0.0, 0.0},
3499
{0.0, 0.0, 1.0, 0.0},
3500
{0.0, 0.0, 0.0, 1.0}};
3501
3502
// Declare (auxiliary) derivative matrix (of polynomial basis).
3503
double dmats_old[4][4] = \
3504
{{1.0, 0.0, 0.0, 0.0},
3505
{0.0, 1.0, 0.0, 0.0},
3506
{0.0, 0.0, 1.0, 0.0},
3507
{0.0, 0.0, 0.0, 1.0}};
3508
3509
// Loop possible derivatives.
3510
for (unsigned int r = 0; r < num_derivatives; r++)
3511
{
3512
// Resetting dmats values to compute next derivative.
3513
for (unsigned int t = 0; t < 4; t++)
3514
{
3515
for (unsigned int u = 0; u < 4; u++)
3516
{
3517
dmats[t][u] = 0.0;
3518
if (t == u)
3519
{
3520
dmats[t][u] = 1.0;
3521
}
3522
3523
}// end loop over 'u'
3524
}// end loop over 't'
3525
3526
// Looping derivative order to generate dmats.
3527
for (unsigned int s = 0; s < n; s++)
3528
{
3529
// Updating dmats_old with new values and resetting dmats.
3530
for (unsigned int t = 0; t < 4; t++)
3531
{
3532
for (unsigned int u = 0; u < 4; u++)
3533
{
3534
dmats_old[t][u] = dmats[t][u];
3535
dmats[t][u] = 0.0;
3536
}// end loop over 'u'
3537
}// end loop over 't'
3538
3539
// Update dmats using an inner product.
3540
if (combinations[r][s] == 0)
3541
{
3542
for (unsigned int t = 0; t < 4; t++)
3543
{
3544
for (unsigned int u = 0; u < 4; u++)
3545
{
3546
for (unsigned int tu = 0; tu < 4; tu++)
3547
{
3548
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3549
}// end loop over 'tu'
3550
}// end loop over 'u'
3551
}// end loop over 't'
3552
}
3553
3554
if (combinations[r][s] == 1)
3555
{
3556
for (unsigned int t = 0; t < 4; t++)
3557
{
3558
for (unsigned int u = 0; u < 4; u++)
3559
{
3560
for (unsigned int tu = 0; tu < 4; tu++)
3561
{
3562
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3563
}// end loop over 'tu'
3564
}// end loop over 'u'
3565
}// end loop over 't'
3566
}
3567
3568
if (combinations[r][s] == 2)
3569
{
3570
for (unsigned int t = 0; t < 4; t++)
3571
{
3572
for (unsigned int u = 0; u < 4; u++)
3573
{
3574
for (unsigned int tu = 0; tu < 4; tu++)
3575
{
3576
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3577
}// end loop over 'tu'
3578
}// end loop over 'u'
3579
}// end loop over 't'
3580
}
3581
3582
}// end loop over 's'
3583
for (unsigned int s = 0; s < 4; s++)
3584
{
3585
for (unsigned int t = 0; t < 4; t++)
3586
{
3587
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3588
}// end loop over 't'
3589
}// end loop over 's'
3590
}// end loop over 'r'
3591
3592
// Transform derivatives back to physical element
3593
for (unsigned int r = 0; r < num_derivatives; r++)
3594
{
3595
for (unsigned int s = 0; s < num_derivatives; s++)
3596
{
3597
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
3598
}// end loop over 's'
3599
}// end loop over 'r'
3600
3601
// Delete pointer to array of derivatives on FIAT element
3602
delete [] derivatives;
3603
3604
// Delete pointer to array of combinations of derivatives and transform
3605
for (unsigned int r = 0; r < num_derivatives; r++)
3606
{
3607
delete [] combinations[r];
3608
}// end loop over 'r'
3609
delete [] combinations;
3610
for (unsigned int r = 0; r < num_derivatives; r++)
3611
{
3612
delete [] transform[r];
3613
}// end loop over 'r'
3614
delete [] transform;
3615
break;
3616
}
3617
case 10:
3618
{
3619
3620
// Array of basisvalues.
3621
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3622
3623
// Declare helper variables.
3624
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3625
3626
// Compute basisvalues.
3627
basisvalues[0] = 1.0;
3628
basisvalues[1] = tmp0;
3629
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3630
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3631
basisvalues[0] *= std::sqrt(0.75);
3632
basisvalues[3] *= std::sqrt(1.25);
3633
basisvalues[2] *= std::sqrt(2.5);
3634
basisvalues[1] *= std::sqrt(7.5);
3635
3636
// Table(s) of coefficients.
3637
static const double coefficients0[4] = \
3638
{0.28867513, 0.0, 0.21081851, -0.074535599};
3639
3640
// Tables of derivatives of the polynomial base (transpose).
3641
static const double dmats0[4][4] = \
3642
{{0.0, 0.0, 0.0, 0.0},
3643
{6.3245553, 0.0, 0.0, 0.0},
3644
{0.0, 0.0, 0.0, 0.0},
3645
{0.0, 0.0, 0.0, 0.0}};
3646
3647
static const double dmats1[4][4] = \
3648
{{0.0, 0.0, 0.0, 0.0},
3649
{3.1622777, 0.0, 0.0, 0.0},
3650
{5.4772256, 0.0, 0.0, 0.0},
3651
{0.0, 0.0, 0.0, 0.0}};
3652
3653
static const double dmats2[4][4] = \
3654
{{0.0, 0.0, 0.0, 0.0},
3655
{3.1622777, 0.0, 0.0, 0.0},
3656
{1.8257419, 0.0, 0.0, 0.0},
3657
{5.1639778, 0.0, 0.0, 0.0}};
3658
3659
// Compute reference derivatives.
3660
// Declare pointer to array of derivatives on FIAT element.
3661
double *derivatives = new double[num_derivatives];
3662
for (unsigned int r = 0; r < num_derivatives; r++)
3663
{
3664
derivatives[r] = 0.0;
3665
}// end loop over 'r'
3666
3667
// Declare derivative matrix (of polynomial basis).
3668
double dmats[4][4] = \
3669
{{1.0, 0.0, 0.0, 0.0},
3670
{0.0, 1.0, 0.0, 0.0},
3671
{0.0, 0.0, 1.0, 0.0},
3672
{0.0, 0.0, 0.0, 1.0}};
3673
3674
// Declare (auxiliary) derivative matrix (of polynomial basis).
3675
double dmats_old[4][4] = \
3676
{{1.0, 0.0, 0.0, 0.0},
3677
{0.0, 1.0, 0.0, 0.0},
3678
{0.0, 0.0, 1.0, 0.0},
3679
{0.0, 0.0, 0.0, 1.0}};
3680
3681
// Loop possible derivatives.
3682
for (unsigned int r = 0; r < num_derivatives; r++)
3683
{
3684
// Resetting dmats values to compute next derivative.
3685
for (unsigned int t = 0; t < 4; t++)
3686
{
3687
for (unsigned int u = 0; u < 4; u++)
3688
{
3689
dmats[t][u] = 0.0;
3690
if (t == u)
3691
{
3692
dmats[t][u] = 1.0;
3693
}
3694
3695
}// end loop over 'u'
3696
}// end loop over 't'
3697
3698
// Looping derivative order to generate dmats.
3699
for (unsigned int s = 0; s < n; s++)
3700
{
3701
// Updating dmats_old with new values and resetting dmats.
3702
for (unsigned int t = 0; t < 4; t++)
3703
{
3704
for (unsigned int u = 0; u < 4; u++)
3705
{
3706
dmats_old[t][u] = dmats[t][u];
3707
dmats[t][u] = 0.0;
3708
}// end loop over 'u'
3709
}// end loop over 't'
3710
3711
// Update dmats using an inner product.
3712
if (combinations[r][s] == 0)
3713
{
3714
for (unsigned int t = 0; t < 4; t++)
3715
{
3716
for (unsigned int u = 0; u < 4; u++)
3717
{
3718
for (unsigned int tu = 0; tu < 4; tu++)
3719
{
3720
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3721
}// end loop over 'tu'
3722
}// end loop over 'u'
3723
}// end loop over 't'
3724
}
3725
3726
if (combinations[r][s] == 1)
3727
{
3728
for (unsigned int t = 0; t < 4; t++)
3729
{
3730
for (unsigned int u = 0; u < 4; u++)
3731
{
3732
for (unsigned int tu = 0; tu < 4; tu++)
3733
{
3734
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3735
}// end loop over 'tu'
3736
}// end loop over 'u'
3737
}// end loop over 't'
3738
}
3739
3740
if (combinations[r][s] == 2)
3741
{
3742
for (unsigned int t = 0; t < 4; t++)
3743
{
3744
for (unsigned int u = 0; u < 4; u++)
3745
{
3746
for (unsigned int tu = 0; tu < 4; tu++)
3747
{
3748
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3749
}// end loop over 'tu'
3750
}// end loop over 'u'
3751
}// end loop over 't'
3752
}
3753
3754
}// end loop over 's'
3755
for (unsigned int s = 0; s < 4; s++)
3756
{
3757
for (unsigned int t = 0; t < 4; t++)
3758
{
3759
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3760
}// end loop over 't'
3761
}// end loop over 's'
3762
}// end loop over 'r'
3763
3764
// Transform derivatives back to physical element
3765
for (unsigned int r = 0; r < num_derivatives; r++)
3766
{
3767
for (unsigned int s = 0; s < num_derivatives; s++)
3768
{
3769
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
3770
}// end loop over 's'
3771
}// end loop over 'r'
3772
3773
// Delete pointer to array of derivatives on FIAT element
3774
delete [] derivatives;
3775
3776
// Delete pointer to array of combinations of derivatives and transform
3777
for (unsigned int r = 0; r < num_derivatives; r++)
3778
{
3779
delete [] combinations[r];
3780
}// end loop over 'r'
3781
delete [] combinations;
3782
for (unsigned int r = 0; r < num_derivatives; r++)
3783
{
3784
delete [] transform[r];
3785
}// end loop over 'r'
3786
delete [] transform;
3787
break;
3788
}
3789
case 11:
3790
{
3791
3792
// Array of basisvalues.
3793
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3794
3795
// Declare helper variables.
3796
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3797
3798
// Compute basisvalues.
3799
basisvalues[0] = 1.0;
3800
basisvalues[1] = tmp0;
3801
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3802
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3803
basisvalues[0] *= std::sqrt(0.75);
3804
basisvalues[3] *= std::sqrt(1.25);
3805
basisvalues[2] *= std::sqrt(2.5);
3806
basisvalues[1] *= std::sqrt(7.5);
3807
3808
// Table(s) of coefficients.
3809
static const double coefficients0[4] = \
3810
{0.28867513, 0.0, 0.0, 0.2236068};
3811
3812
// Tables of derivatives of the polynomial base (transpose).
3813
static const double dmats0[4][4] = \
3814
{{0.0, 0.0, 0.0, 0.0},
3815
{6.3245553, 0.0, 0.0, 0.0},
3816
{0.0, 0.0, 0.0, 0.0},
3817
{0.0, 0.0, 0.0, 0.0}};
3818
3819
static const double dmats1[4][4] = \
3820
{{0.0, 0.0, 0.0, 0.0},
3821
{3.1622777, 0.0, 0.0, 0.0},
3822
{5.4772256, 0.0, 0.0, 0.0},
3823
{0.0, 0.0, 0.0, 0.0}};
3824
3825
static const double dmats2[4][4] = \
3826
{{0.0, 0.0, 0.0, 0.0},
3827
{3.1622777, 0.0, 0.0, 0.0},
3828
{1.8257419, 0.0, 0.0, 0.0},
3829
{5.1639778, 0.0, 0.0, 0.0}};
3830
3831
// Compute reference derivatives.
3832
// Declare pointer to array of derivatives on FIAT element.
3833
double *derivatives = new double[num_derivatives];
3834
for (unsigned int r = 0; r < num_derivatives; r++)
3835
{
3836
derivatives[r] = 0.0;
3837
}// end loop over 'r'
3838
3839
// Declare derivative matrix (of polynomial basis).
3840
double dmats[4][4] = \
3841
{{1.0, 0.0, 0.0, 0.0},
3842
{0.0, 1.0, 0.0, 0.0},
3843
{0.0, 0.0, 1.0, 0.0},
3844
{0.0, 0.0, 0.0, 1.0}};
3845
3846
// Declare (auxiliary) derivative matrix (of polynomial basis).
3847
double dmats_old[4][4] = \
3848
{{1.0, 0.0, 0.0, 0.0},
3849
{0.0, 1.0, 0.0, 0.0},
3850
{0.0, 0.0, 1.0, 0.0},
3851
{0.0, 0.0, 0.0, 1.0}};
3852
3853
// Loop possible derivatives.
3854
for (unsigned int r = 0; r < num_derivatives; r++)
3855
{
3856
// Resetting dmats values to compute next derivative.
3857
for (unsigned int t = 0; t < 4; t++)
3858
{
3859
for (unsigned int u = 0; u < 4; u++)
3860
{
3861
dmats[t][u] = 0.0;
3862
if (t == u)
3863
{
3864
dmats[t][u] = 1.0;
3865
}
3866
3867
}// end loop over 'u'
3868
}// end loop over 't'
3869
3870
// Looping derivative order to generate dmats.
3871
for (unsigned int s = 0; s < n; s++)
3872
{
3873
// Updating dmats_old with new values and resetting dmats.
3874
for (unsigned int t = 0; t < 4; t++)
3875
{
3876
for (unsigned int u = 0; u < 4; u++)
3877
{
3878
dmats_old[t][u] = dmats[t][u];
3879
dmats[t][u] = 0.0;
3880
}// end loop over 'u'
3881
}// end loop over 't'
3882
3883
// Update dmats using an inner product.
3884
if (combinations[r][s] == 0)
3885
{
3886
for (unsigned int t = 0; t < 4; t++)
3887
{
3888
for (unsigned int u = 0; u < 4; u++)
3889
{
3890
for (unsigned int tu = 0; tu < 4; tu++)
3891
{
3892
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3893
}// end loop over 'tu'
3894
}// end loop over 'u'
3895
}// end loop over 't'
3896
}
3897
3898
if (combinations[r][s] == 1)
3899
{
3900
for (unsigned int t = 0; t < 4; t++)
3901
{
3902
for (unsigned int u = 0; u < 4; u++)
3903
{
3904
for (unsigned int tu = 0; tu < 4; tu++)
3905
{
3906
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3907
}// end loop over 'tu'
3908
}// end loop over 'u'
3909
}// end loop over 't'
3910
}
3911
3912
if (combinations[r][s] == 2)
3913
{
3914
for (unsigned int t = 0; t < 4; t++)
3915
{
3916
for (unsigned int u = 0; u < 4; u++)
3917
{
3918
for (unsigned int tu = 0; tu < 4; tu++)
3919
{
3920
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3921
}// end loop over 'tu'
3922
}// end loop over 'u'
3923
}// end loop over 't'
3924
}
3925
3926
}// end loop over 's'
3927
for (unsigned int s = 0; s < 4; s++)
3928
{
3929
for (unsigned int t = 0; t < 4; t++)
3930
{
3931
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3932
}// end loop over 't'
3933
}// end loop over 's'
3934
}// end loop over 'r'
3935
3936
// Transform derivatives back to physical element
3937
for (unsigned int r = 0; r < num_derivatives; r++)
3938
{
3939
for (unsigned int s = 0; s < num_derivatives; s++)
3940
{
3941
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
3942
}// end loop over 's'
3943
}// end loop over 'r'
3944
3945
// Delete pointer to array of derivatives on FIAT element
3946
delete [] derivatives;
3947
3948
// Delete pointer to array of combinations of derivatives and transform
3949
for (unsigned int r = 0; r < num_derivatives; r++)
3950
{
3951
delete [] combinations[r];
3952
}// end loop over 'r'
3953
delete [] combinations;
3954
for (unsigned int r = 0; r < num_derivatives; r++)
3955
{
3956
delete [] transform[r];
3957
}// end loop over 'r'
3958
delete [] transform;
3959
break;
3960
}
3961
}
3962
3963
}
3964
3965
/// Evaluate order n derivatives of all basis functions at given point in cell
3966
virtual void evaluate_basis_derivatives_all(unsigned int n,
3967
double* values,
3968
const double* coordinates,
3969
const ufc::cell& c) const
3970
{
3971
// Compute number of derivatives.
3972
unsigned int num_derivatives = 1;
3973
for (unsigned int r = 0; r < n; r++)
3974
{
3975
num_derivatives *= 3;
3976
}// end loop over 'r'
3977
3978
// Helper variable to hold values of a single dof.
3979
double *dof_values = new double[3*num_derivatives];
3980
for (unsigned int r = 0; r < 3*num_derivatives; r++)
3981
{
3982
dof_values[r] = 0.0;
3983
}// end loop over 'r'
3984
3985
// Loop dofs and call evaluate_basis_derivatives.
3986
for (unsigned int r = 0; r < 12; r++)
3987
{
3988
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
3989
for (unsigned int s = 0; s < 3*num_derivatives; s++)
3990
{
3991
values[r*3*num_derivatives + s] = dof_values[s];
3992
}// end loop over 's'
3993
}// end loop over 'r'
3994
3995
// Delete pointer.
3996
delete [] dof_values;
3997
}
3998
3999
/// Evaluate linear functional for dof i on the function f
4000
virtual double evaluate_dof(unsigned int i,
4001
const ufc::function& f,
4002
const ufc::cell& c) const
4003
{
4004
// Declare variables for result of evaluation.
4005
double vals[3];
4006
4007
// Declare variable for physical coordinates.
4008
double y[3];
4009
const double * const * x = c.coordinates;
4010
switch (i)
4011
{
4012
case 0:
4013
{
4014
y[0] = x[0][0];
4015
y[1] = x[0][1];
4016
y[2] = x[0][2];
4017
f.evaluate(vals, y, c);
4018
return vals[0];
4019
break;
4020
}
4021
case 1:
4022
{
4023
y[0] = x[1][0];
4024
y[1] = x[1][1];
4025
y[2] = x[1][2];
4026
f.evaluate(vals, y, c);
4027
return vals[0];
4028
break;
4029
}
4030
case 2:
4031
{
4032
y[0] = x[2][0];
4033
y[1] = x[2][1];
4034
y[2] = x[2][2];
4035
f.evaluate(vals, y, c);
4036
return vals[0];
4037
break;
4038
}
4039
case 3:
4040
{
4041
y[0] = x[3][0];
4042
y[1] = x[3][1];
4043
y[2] = x[3][2];
4044
f.evaluate(vals, y, c);
4045
return vals[0];
4046
break;
4047
}
4048
case 4:
4049
{
4050
y[0] = x[0][0];
4051
y[1] = x[0][1];
4052
y[2] = x[0][2];
4053
f.evaluate(vals, y, c);
4054
return vals[1];
4055
break;
4056
}
4057
case 5:
4058
{
4059
y[0] = x[1][0];
4060
y[1] = x[1][1];
4061
y[2] = x[1][2];
4062
f.evaluate(vals, y, c);
4063
return vals[1];
4064
break;
4065
}
4066
case 6:
4067
{
4068
y[0] = x[2][0];
4069
y[1] = x[2][1];
4070
y[2] = x[2][2];
4071
f.evaluate(vals, y, c);
4072
return vals[1];
4073
break;
4074
}
4075
case 7:
4076
{
4077
y[0] = x[3][0];
4078
y[1] = x[3][1];
4079
y[2] = x[3][2];
4080
f.evaluate(vals, y, c);
4081
return vals[1];
4082
break;
4083
}
4084
case 8:
4085
{
4086
y[0] = x[0][0];
4087
y[1] = x[0][1];
4088
y[2] = x[0][2];
4089
f.evaluate(vals, y, c);
4090
return vals[2];
4091
break;
4092
}
4093
case 9:
4094
{
4095
y[0] = x[1][0];
4096
y[1] = x[1][1];
4097
y[2] = x[1][2];
4098
f.evaluate(vals, y, c);
4099
return vals[2];
4100
break;
4101
}
4102
case 10:
4103
{
4104
y[0] = x[2][0];
4105
y[1] = x[2][1];
4106
y[2] = x[2][2];
4107
f.evaluate(vals, y, c);
4108
return vals[2];
4109
break;
4110
}
4111
case 11:
4112
{
4113
y[0] = x[3][0];
4114
y[1] = x[3][1];
4115
y[2] = x[3][2];
4116
f.evaluate(vals, y, c);
4117
return vals[2];
4118
break;
4119
}
4120
}
4121
4122
return 0.0;
4123
}
4124
4125
/// Evaluate linear functionals for all dofs on the function f
4126
virtual void evaluate_dofs(double* values,
4127
const ufc::function& f,
4128
const ufc::cell& c) const
4129
{
4130
// Declare variables for result of evaluation.
4131
double vals[3];
4132
4133
// Declare variable for physical coordinates.
4134
double y[3];
4135
const double * const * x = c.coordinates;
4136
y[0] = x[0][0];
4137
y[1] = x[0][1];
4138
y[2] = x[0][2];
4139
f.evaluate(vals, y, c);
4140
values[0] = vals[0];
4141
y[0] = x[1][0];
4142
y[1] = x[1][1];
4143
y[2] = x[1][2];
4144
f.evaluate(vals, y, c);
4145
values[1] = vals[0];
4146
y[0] = x[2][0];
4147
y[1] = x[2][1];
4148
y[2] = x[2][2];
4149
f.evaluate(vals, y, c);
4150
values[2] = vals[0];
4151
y[0] = x[3][0];
4152
y[1] = x[3][1];
4153
y[2] = x[3][2];
4154
f.evaluate(vals, y, c);
4155
values[3] = vals[0];
4156
y[0] = x[0][0];
4157
y[1] = x[0][1];
4158
y[2] = x[0][2];
4159
f.evaluate(vals, y, c);
4160
values[4] = vals[1];
4161
y[0] = x[1][0];
4162
y[1] = x[1][1];
4163
y[2] = x[1][2];
4164
f.evaluate(vals, y, c);
4165
values[5] = vals[1];
4166
y[0] = x[2][0];
4167
y[1] = x[2][1];
4168
y[2] = x[2][2];
4169
f.evaluate(vals, y, c);
4170
values[6] = vals[1];
4171
y[0] = x[3][0];
4172
y[1] = x[3][1];
4173
y[2] = x[3][2];
4174
f.evaluate(vals, y, c);
4175
values[7] = vals[1];
4176
y[0] = x[0][0];
4177
y[1] = x[0][1];
4178
y[2] = x[0][2];
4179
f.evaluate(vals, y, c);
4180
values[8] = vals[2];
4181
y[0] = x[1][0];
4182
y[1] = x[1][1];
4183
y[2] = x[1][2];
4184
f.evaluate(vals, y, c);
4185
values[9] = vals[2];
4186
y[0] = x[2][0];
4187
y[1] = x[2][1];
4188
y[2] = x[2][2];
4189
f.evaluate(vals, y, c);
4190
values[10] = vals[2];
4191
y[0] = x[3][0];
4192
y[1] = x[3][1];
4193
y[2] = x[3][2];
4194
f.evaluate(vals, y, c);
4195
values[11] = vals[2];
4196
}
4197
4198
/// Interpolate vertex values from dof values
4199
virtual void interpolate_vertex_values(double* vertex_values,
4200
const double* dof_values,
4201
const ufc::cell& c) const
4202
{
4203
// Evaluate function and change variables
4204
vertex_values[0] = dof_values[0];
4205
vertex_values[3] = dof_values[1];
4206
vertex_values[6] = dof_values[2];
4207
vertex_values[9] = dof_values[3];
4208
// Evaluate function and change variables
4209
vertex_values[1] = dof_values[4];
4210
vertex_values[4] = dof_values[5];
4211
vertex_values[7] = dof_values[6];
4212
vertex_values[10] = dof_values[7];
4213
// Evaluate function and change variables
4214
vertex_values[2] = dof_values[8];
4215
vertex_values[5] = dof_values[9];
4216
vertex_values[8] = dof_values[10];
4217
vertex_values[11] = dof_values[11];
4218
}
4219
4220
/// Map coordinate xhat from reference cell to coordinate x in cell
4221
virtual void map_from_reference_cell(double* x,
4222
const double* xhat,
4223
const ufc::cell& c) const
4224
{
4225
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4226
}
4227
4228
/// Map from coordinate x in cell to coordinate xhat in reference cell
4229
virtual void map_to_reference_cell(double* xhat,
4230
const double* x,
4231
const ufc::cell& c) const
4232
{
4233
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4234
}
4235
4236
/// Return the number of sub elements (for a mixed element)
4237
virtual unsigned int num_sub_elements() const
4238
{
4239
return 3;
4240
}
4241
4242
/// Create a new finite element for sub element i (for a mixed element)
4243
virtual ufc::finite_element* create_sub_element(unsigned int i) const
4244
{
4245
switch (i)
4246
{
4247
case 0:
4248
{
4249
return new navierstokes_finite_element_0();
4250
break;
4251
}
4252
case 1:
4253
{
4254
return new navierstokes_finite_element_0();
4255
break;
4256
}
4257
case 2:
4258
{
4259
return new navierstokes_finite_element_0();
4260
break;
4261
}
4262
}
4263
4264
return 0;
4265
}
4266
4267
/// Create a new class instance
4268
virtual ufc::finite_element* create() const
4269
{
4270
return new navierstokes_finite_element_1();
4271
}
4272
4273
};
4274
4275
/// This class defines the interface for a local-to-global mapping of
4276
/// degrees of freedom (dofs).
4277
4278
class navierstokes_dofmap_0: public ufc::dofmap
4279
{
4280
private:
4281
4282
unsigned int _global_dimension;
4283
public:
4284
4285
/// Constructor
4286
navierstokes_dofmap_0() : ufc::dofmap()
4287
{
4288
_global_dimension = 0;
4289
}
4290
4291
/// Destructor
4292
virtual ~navierstokes_dofmap_0()
4293
{
4294
// Do nothing
4295
}
4296
4297
/// Return a string identifying the dofmap
4298
virtual const char* signature() const
4299
{
4300
return "FFC dofmap for FiniteElement('Lagrange', Cell('tetrahedron', Space(3)), 1, None)";
4301
}
4302
4303
/// Return true iff mesh entities of topological dimension d are needed
4304
virtual bool needs_mesh_entities(unsigned int d) const
4305
{
4306
switch (d)
4307
{
4308
case 0:
4309
{
4310
return true;
4311
break;
4312
}
4313
case 1:
4314
{
4315
return false;
4316
break;
4317
}
4318
case 2:
4319
{
4320
return false;
4321
break;
4322
}
4323
case 3:
4324
{
4325
return false;
4326
break;
4327
}
4328
}
4329
4330
return false;
4331
}
4332
4333
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
4334
virtual bool init_mesh(const ufc::mesh& m)
4335
{
4336
_global_dimension = m.num_entities[0];
4337
return false;
4338
}
4339
4340
/// Initialize dofmap for given cell
4341
virtual void init_cell(const ufc::mesh& m,
4342
const ufc::cell& c)
4343
{
4344
// Do nothing
4345
}
4346
4347
/// Finish initialization of dofmap for cells
4348
virtual void init_cell_finalize()
4349
{
4350
// Do nothing
4351
}
4352
4353
/// Return the topological dimension of the associated cell shape
4354
virtual unsigned int topological_dimension() const
4355
{
4356
return 3;
4357
}
4358
4359
/// Return the geometric dimension of the associated cell shape
4360
virtual unsigned int geometric_dimension() const
4361
{
4362
return 3;
4363
}
4364
4365
/// Return the dimension of the global finite element function space
4366
virtual unsigned int global_dimension() const
4367
{
4368
return _global_dimension;
4369
}
4370
4371
/// Return the dimension of the local finite element function space for a cell
4372
virtual unsigned int local_dimension(const ufc::cell& c) const
4373
{
4374
return 4;
4375
}
4376
4377
/// Return the maximum dimension of the local finite element function space
4378
virtual unsigned int max_local_dimension() const
4379
{
4380
return 4;
4381
}
4382
4383
/// Return the number of dofs on each cell facet
4384
virtual unsigned int num_facet_dofs() const
4385
{
4386
return 3;
4387
}
4388
4389
/// Return the number of dofs associated with each cell entity of dimension d
4390
virtual unsigned int num_entity_dofs(unsigned int d) const
4391
{
4392
switch (d)
4393
{
4394
case 0:
4395
{
4396
return 1;
4397
break;
4398
}
4399
case 1:
4400
{
4401
return 0;
4402
break;
4403
}
4404
case 2:
4405
{
4406
return 0;
4407
break;
4408
}
4409
case 3:
4410
{
4411
return 0;
4412
break;
4413
}
4414
}
4415
4416
return 0;
4417
}
4418
4419
/// Tabulate the local-to-global mapping of dofs on a cell
4420
virtual void tabulate_dofs(unsigned int* dofs,
4421
const ufc::mesh& m,
4422
const ufc::cell& c) const
4423
{
4424
dofs[0] = c.entity_indices[0][0];
4425
dofs[1] = c.entity_indices[0][1];
4426
dofs[2] = c.entity_indices[0][2];
4427
dofs[3] = c.entity_indices[0][3];
4428
}
4429
4430
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4431
virtual void tabulate_facet_dofs(unsigned int* dofs,
4432
unsigned int facet) const
4433
{
4434
switch (facet)
4435
{
4436
case 0:
4437
{
4438
dofs[0] = 1;
4439
dofs[1] = 2;
4440
dofs[2] = 3;
4441
break;
4442
}
4443
case 1:
4444
{
4445
dofs[0] = 0;
4446
dofs[1] = 2;
4447
dofs[2] = 3;
4448
break;
4449
}
4450
case 2:
4451
{
4452
dofs[0] = 0;
4453
dofs[1] = 1;
4454
dofs[2] = 3;
4455
break;
4456
}
4457
case 3:
4458
{
4459
dofs[0] = 0;
4460
dofs[1] = 1;
4461
dofs[2] = 2;
4462
break;
4463
}
4464
}
4465
4466
}
4467
4468
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4469
virtual void tabulate_entity_dofs(unsigned int* dofs,
4470
unsigned int d, unsigned int i) const
4471
{
4472
if (d > 3)
4473
{
4474
std::cerr << "*** FFC warning: " << "d is larger than dimension (3)" << std::endl;
4475
}
4476
4477
switch (d)
4478
{
4479
case 0:
4480
{
4481
if (i > 3)
4482
{
4483
std::cerr << "*** FFC warning: " << "i is larger than number of entities (3)" << std::endl;
4484
}
4485
4486
switch (i)
4487
{
4488
case 0:
4489
{
4490
dofs[0] = 0;
4491
break;
4492
}
4493
case 1:
4494
{
4495
dofs[0] = 1;
4496
break;
4497
}
4498
case 2:
4499
{
4500
dofs[0] = 2;
4501
break;
4502
}
4503
case 3:
4504
{
4505
dofs[0] = 3;
4506
break;
4507
}
4508
}
4509
4510
break;
4511
}
4512
case 1:
4513
{
4514
4515
break;
4516
}
4517
case 2:
4518
{
4519
4520
break;
4521
}
4522
case 3:
4523
{
4524
4525
break;
4526
}
4527
}
4528
4529
}
4530
4531
/// Tabulate the coordinates of all dofs on a cell
4532
virtual void tabulate_coordinates(double** coordinates,
4533
const ufc::cell& c) const
4534
{
4535
const double * const * x = c.coordinates;
4536
4537
coordinates[0][0] = x[0][0];
4538
coordinates[0][1] = x[0][1];
4539
coordinates[0][2] = x[0][2];
4540
coordinates[1][0] = x[1][0];
4541
coordinates[1][1] = x[1][1];
4542
coordinates[1][2] = x[1][2];
4543
coordinates[2][0] = x[2][0];
4544
coordinates[2][1] = x[2][1];
4545
coordinates[2][2] = x[2][2];
4546
coordinates[3][0] = x[3][0];
4547
coordinates[3][1] = x[3][1];
4548
coordinates[3][2] = x[3][2];
4549
}
4550
4551
/// Return the number of sub dofmaps (for a mixed element)
4552
virtual unsigned int num_sub_dofmaps() const
4553
{
4554
return 0;
4555
}
4556
4557
/// Create a new dofmap for sub dofmap i (for a mixed element)
4558
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
4559
{
4560
return 0;
4561
}
4562
4563
/// Create a new class instance
4564
virtual ufc::dofmap* create() const
4565
{
4566
return new navierstokes_dofmap_0();
4567
}
4568
4569
};
4570
4571
/// This class defines the interface for a local-to-global mapping of
4572
/// degrees of freedom (dofs).
4573
4574
class navierstokes_dofmap_1: public ufc::dofmap
4575
{
4576
private:
4577
4578
unsigned int _global_dimension;
4579
public:
4580
4581
/// Constructor
4582
navierstokes_dofmap_1() : ufc::dofmap()
4583
{
4584
_global_dimension = 0;
4585
}
4586
4587
/// Destructor
4588
virtual ~navierstokes_dofmap_1()
4589
{
4590
// Do nothing
4591
}
4592
4593
/// Return a string identifying the dofmap
4594
virtual const char* signature() const
4595
{
4596
return "FFC dofmap for VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None)";
4597
}
4598
4599
/// Return true iff mesh entities of topological dimension d are needed
4600
virtual bool needs_mesh_entities(unsigned int d) const
4601
{
4602
switch (d)
4603
{
4604
case 0:
4605
{
4606
return true;
4607
break;
4608
}
4609
case 1:
4610
{
4611
return false;
4612
break;
4613
}
4614
case 2:
4615
{
4616
return false;
4617
break;
4618
}
4619
case 3:
4620
{
4621
return false;
4622
break;
4623
}
4624
}
4625
4626
return false;
4627
}
4628
4629
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
4630
virtual bool init_mesh(const ufc::mesh& m)
4631
{
4632
_global_dimension = 3*m.num_entities[0];
4633
return false;
4634
}
4635
4636
/// Initialize dofmap for given cell
4637
virtual void init_cell(const ufc::mesh& m,
4638
const ufc::cell& c)
4639
{
4640
// Do nothing
4641
}
4642
4643
/// Finish initialization of dofmap for cells
4644
virtual void init_cell_finalize()
4645
{
4646
// Do nothing
4647
}
4648
4649
/// Return the topological dimension of the associated cell shape
4650
virtual unsigned int topological_dimension() const
4651
{
4652
return 3;
4653
}
4654
4655
/// Return the geometric dimension of the associated cell shape
4656
virtual unsigned int geometric_dimension() const
4657
{
4658
return 3;
4659
}
4660
4661
/// Return the dimension of the global finite element function space
4662
virtual unsigned int global_dimension() const
4663
{
4664
return _global_dimension;
4665
}
4666
4667
/// Return the dimension of the local finite element function space for a cell
4668
virtual unsigned int local_dimension(const ufc::cell& c) const
4669
{
4670
return 12;
4671
}
4672
4673
/// Return the maximum dimension of the local finite element function space
4674
virtual unsigned int max_local_dimension() const
4675
{
4676
return 12;
4677
}
4678
4679
/// Return the number of dofs on each cell facet
4680
virtual unsigned int num_facet_dofs() const
4681
{
4682
return 9;
4683
}
4684
4685
/// Return the number of dofs associated with each cell entity of dimension d
4686
virtual unsigned int num_entity_dofs(unsigned int d) const
4687
{
4688
switch (d)
4689
{
4690
case 0:
4691
{
4692
return 3;
4693
break;
4694
}
4695
case 1:
4696
{
4697
return 0;
4698
break;
4699
}
4700
case 2:
4701
{
4702
return 0;
4703
break;
4704
}
4705
case 3:
4706
{
4707
return 0;
4708
break;
4709
}
4710
}
4711
4712
return 0;
4713
}
4714
4715
/// Tabulate the local-to-global mapping of dofs on a cell
4716
virtual void tabulate_dofs(unsigned int* dofs,
4717
const ufc::mesh& m,
4718
const ufc::cell& c) const
4719
{
4720
unsigned int offset = 0;
4721
dofs[0] = offset + c.entity_indices[0][0];
4722
dofs[1] = offset + c.entity_indices[0][1];
4723
dofs[2] = offset + c.entity_indices[0][2];
4724
dofs[3] = offset + c.entity_indices[0][3];
4725
offset += m.num_entities[0];
4726
dofs[4] = offset + c.entity_indices[0][0];
4727
dofs[5] = offset + c.entity_indices[0][1];
4728
dofs[6] = offset + c.entity_indices[0][2];
4729
dofs[7] = offset + c.entity_indices[0][3];
4730
offset += m.num_entities[0];
4731
dofs[8] = offset + c.entity_indices[0][0];
4732
dofs[9] = offset + c.entity_indices[0][1];
4733
dofs[10] = offset + c.entity_indices[0][2];
4734
dofs[11] = offset + c.entity_indices[0][3];
4735
offset += m.num_entities[0];
4736
}
4737
4738
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4739
virtual void tabulate_facet_dofs(unsigned int* dofs,
4740
unsigned int facet) const
4741
{
4742
switch (facet)
4743
{
4744
case 0:
4745
{
4746
dofs[0] = 1;
4747
dofs[1] = 2;
4748
dofs[2] = 3;
4749
dofs[3] = 5;
4750
dofs[4] = 6;
4751
dofs[5] = 7;
4752
dofs[6] = 9;
4753
dofs[7] = 10;
4754
dofs[8] = 11;
4755
break;
4756
}
4757
case 1:
4758
{
4759
dofs[0] = 0;
4760
dofs[1] = 2;
4761
dofs[2] = 3;
4762
dofs[3] = 4;
4763
dofs[4] = 6;
4764
dofs[5] = 7;
4765
dofs[6] = 8;
4766
dofs[7] = 10;
4767
dofs[8] = 11;
4768
break;
4769
}
4770
case 2:
4771
{
4772
dofs[0] = 0;
4773
dofs[1] = 1;
4774
dofs[2] = 3;
4775
dofs[3] = 4;
4776
dofs[4] = 5;
4777
dofs[5] = 7;
4778
dofs[6] = 8;
4779
dofs[7] = 9;
4780
dofs[8] = 11;
4781
break;
4782
}
4783
case 3:
4784
{
4785
dofs[0] = 0;
4786
dofs[1] = 1;
4787
dofs[2] = 2;
4788
dofs[3] = 4;
4789
dofs[4] = 5;
4790
dofs[5] = 6;
4791
dofs[6] = 8;
4792
dofs[7] = 9;
4793
dofs[8] = 10;
4794
break;
4795
}
4796
}
4797
4798
}
4799
4800
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4801
virtual void tabulate_entity_dofs(unsigned int* dofs,
4802
unsigned int d, unsigned int i) const
4803
{
4804
if (d > 3)
4805
{
4806
std::cerr << "*** FFC warning: " << "d is larger than dimension (3)" << std::endl;
4807
}
4808
4809
switch (d)
4810
{
4811
case 0:
4812
{
4813
if (i > 3)
4814
{
4815
std::cerr << "*** FFC warning: " << "i is larger than number of entities (3)" << std::endl;
4816
}
4817
4818
switch (i)
4819
{
4820
case 0:
4821
{
4822
dofs[0] = 0;
4823
dofs[1] = 4;
4824
dofs[2] = 8;
4825
break;
4826
}
4827
case 1:
4828
{
4829
dofs[0] = 1;
4830
dofs[1] = 5;
4831
dofs[2] = 9;
4832
break;
4833
}
4834
case 2:
4835
{
4836
dofs[0] = 2;
4837
dofs[1] = 6;
4838
dofs[2] = 10;
4839
break;
4840
}
4841
case 3:
4842
{
4843
dofs[0] = 3;
4844
dofs[1] = 7;
4845
dofs[2] = 11;
4846
break;
4847
}
4848
}
4849
4850
break;
4851
}
4852
case 1:
4853
{
4854
4855
break;
4856
}
4857
case 2:
4858
{
4859
4860
break;
4861
}
4862
case 3:
4863
{
4864
4865
break;
4866
}
4867
}
4868
4869
}
4870
4871
/// Tabulate the coordinates of all dofs on a cell
4872
virtual void tabulate_coordinates(double** coordinates,
4873
const ufc::cell& c) const
4874
{
4875
const double * const * x = c.coordinates;
4876
4877
coordinates[0][0] = x[0][0];
4878
coordinates[0][1] = x[0][1];
4879
coordinates[0][2] = x[0][2];
4880
coordinates[1][0] = x[1][0];
4881
coordinates[1][1] = x[1][1];
4882
coordinates[1][2] = x[1][2];
4883
coordinates[2][0] = x[2][0];
4884
coordinates[2][1] = x[2][1];
4885
coordinates[2][2] = x[2][2];
4886
coordinates[3][0] = x[3][0];
4887
coordinates[3][1] = x[3][1];
4888
coordinates[3][2] = x[3][2];
4889
coordinates[4][0] = x[0][0];
4890
coordinates[4][1] = x[0][1];
4891
coordinates[4][2] = x[0][2];
4892
coordinates[5][0] = x[1][0];
4893
coordinates[5][1] = x[1][1];
4894
coordinates[5][2] = x[1][2];
4895
coordinates[6][0] = x[2][0];
4896
coordinates[6][1] = x[2][1];
4897
coordinates[6][2] = x[2][2];
4898
coordinates[7][0] = x[3][0];
4899
coordinates[7][1] = x[3][1];
4900
coordinates[7][2] = x[3][2];
4901
coordinates[8][0] = x[0][0];
4902
coordinates[8][1] = x[0][1];
4903
coordinates[8][2] = x[0][2];
4904
coordinates[9][0] = x[1][0];
4905
coordinates[9][1] = x[1][1];
4906
coordinates[9][2] = x[1][2];
4907
coordinates[10][0] = x[2][0];
4908
coordinates[10][1] = x[2][1];
4909
coordinates[10][2] = x[2][2];
4910
coordinates[11][0] = x[3][0];
4911
coordinates[11][1] = x[3][1];
4912
coordinates[11][2] = x[3][2];
4913
}
4914
4915
/// Return the number of sub dofmaps (for a mixed element)
4916
virtual unsigned int num_sub_dofmaps() const
4917
{
4918
return 3;
4919
}
4920
4921
/// Create a new dofmap for sub dofmap i (for a mixed element)
4922
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
4923
{
4924
switch (i)
4925
{
4926
case 0:
4927
{
4928
return new navierstokes_dofmap_0();
4929
break;
4930
}
4931
case 1:
4932
{
4933
return new navierstokes_dofmap_0();
4934
break;
4935
}
4936
case 2:
4937
{
4938
return new navierstokes_dofmap_0();
4939
break;
4940
}
4941
}
4942
4943
return 0;
4944
}
4945
4946
/// Create a new class instance
4947
virtual ufc::dofmap* create() const
4948
{
4949
return new navierstokes_dofmap_1();
4950
}
4951
4952
};
4953
4954
/// This class defines the interface for the tabulation of the cell
4955
/// tensor corresponding to the local contribution to a form from
4956
/// the integral over a cell.
4957
4958
class navierstokes_cell_integral_0_0: public ufc::cell_integral
4959
{
4960
public:
4961
4962
/// Constructor
4963
navierstokes_cell_integral_0_0() : ufc::cell_integral()
4964
{
4965
// Do nothing
4966
}
4967
4968
/// Destructor
4969
virtual ~navierstokes_cell_integral_0_0()
4970
{
4971
// Do nothing
4972
}
4973
4974
/// Tabulate the tensor for the contribution from a local cell
4975
virtual void tabulate_tensor(double* A,
4976
const double * const * w,
4977
const ufc::cell& c) const
4978
{
4979
// Extract vertex coordinates
4980
const double * const * x = c.coordinates;
4981
4982
// Compute Jacobian of affine map from reference cell
4983
const double J_00 = x[1][0] - x[0][0];
4984
const double J_01 = x[2][0] - x[0][0];
4985
const double J_02 = x[3][0] - x[0][0];
4986
const double J_10 = x[1][1] - x[0][1];
4987
const double J_11 = x[2][1] - x[0][1];
4988
const double J_12 = x[3][1] - x[0][1];
4989
const double J_20 = x[1][2] - x[0][2];
4990
const double J_21 = x[2][2] - x[0][2];
4991
const double J_22 = x[3][2] - x[0][2];
4992
4993
// Compute sub determinants
4994
const double d_00 = J_11*J_22 - J_12*J_21;
4995
const double d_01 = J_12*J_20 - J_10*J_22;
4996
const double d_02 = J_10*J_21 - J_11*J_20;
4997
const double d_10 = J_02*J_21 - J_01*J_22;
4998
const double d_11 = J_00*J_22 - J_02*J_20;
4999
const double d_12 = J_01*J_20 - J_00*J_21;
5000
const double d_20 = J_01*J_12 - J_02*J_11;
5001
const double d_21 = J_02*J_10 - J_00*J_12;
5002
const double d_22 = J_00*J_11 - J_01*J_10;
5003
5004
// Compute determinant of Jacobian
5005
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
5006
5007
// Compute inverse of Jacobian
5008
const double K_00 = d_00 / detJ;
5009
const double K_01 = d_10 / detJ;
5010
const double K_02 = d_20 / detJ;
5011
const double K_10 = d_01 / detJ;
5012
const double K_11 = d_11 / detJ;
5013
const double K_12 = d_21 / detJ;
5014
const double K_20 = d_02 / detJ;
5015
const double K_21 = d_12 / detJ;
5016
const double K_22 = d_22 / detJ;
5017
5018
// Set scale factor
5019
const double det = std::abs(detJ);
5020
5021
// Cell Volume.
5022
5023
// Compute circumradius.
5024
5025
5026
// Facet Area (divide by two because 'det' is scaled by area of reference triangle).
5027
5028
// Array of quadrature weights.
5029
static const double W4[4] = {0.041666667, 0.041666667, 0.041666667, 0.041666667};
5030
// Quadrature points on the UFC reference element: (0.5854102, 0.1381966, 0.1381966), (0.1381966, 0.5854102, 0.1381966), (0.1381966, 0.1381966, 0.5854102), (0.1381966, 0.1381966, 0.1381966)
5031
5032
// Value of basis functions at quadrature points.
5033
static const double FE0_C0[4][4] = \
5034
{{0.1381966, 0.5854102, 0.1381966, 0.1381966},
5035
{0.1381966, 0.1381966, 0.5854102, 0.1381966},
5036
{0.1381966, 0.1381966, 0.1381966, 0.5854102},
5037
{0.5854102, 0.1381966, 0.1381966, 0.1381966}};
5038
5039
// Array of non-zero columns
5040
static const unsigned int nzc8[4] = {8, 9, 10, 11};
5041
5042
// Array of non-zero columns
5043
static const unsigned int nzc4[4] = {4, 5, 6, 7};
5044
5045
// Array of non-zero columns
5046
static const unsigned int nzc0[4] = {0, 1, 2, 3};
5047
5048
static const double FE0_C0_D001[4][2] = \
5049
{{-1.0, 1.0},
5050
{-1.0, 1.0},
5051
{-1.0, 1.0},
5052
{-1.0, 1.0}};
5053
5054
// Array of non-zero columns
5055
static const unsigned int nzc9[2] = {8, 11};
5056
5057
// Array of non-zero columns
5058
static const unsigned int nzc5[2] = {4, 7};
5059
5060
// Array of non-zero columns
5061
static const unsigned int nzc6[2] = {4, 6};
5062
5063
// Array of non-zero columns
5064
static const unsigned int nzc10[2] = {8, 10};
5065
5066
// Array of non-zero columns
5067
static const unsigned int nzc7[2] = {4, 5};
5068
5069
// Array of non-zero columns
5070
static const unsigned int nzc3[2] = {0, 1};
5071
5072
// Array of non-zero columns
5073
static const unsigned int nzc2[2] = {0, 2};
5074
5075
// Array of non-zero columns
5076
static const unsigned int nzc1[2] = {0, 3};
5077
5078
// Array of non-zero columns
5079
static const unsigned int nzc11[2] = {8, 9};
5080
5081
// Reset values in the element tensor.
5082
for (unsigned int r = 0; r < 144; r++)
5083
{
5084
A[r] = 0.0;
5085
}// end loop over 'r'
5086
// Number of operations to compute geometry constants: 9.
5087
double G[9];
5088
G[0] = K_20*det;
5089
G[1] = K_21*det;
5090
G[2] = K_22*det;
5091
G[3] = K_10*det;
5092
G[4] = K_11*det;
5093
G[5] = K_12*det;
5094
G[6] = K_00*det;
5095
G[7] = K_01*det;
5096
G[8] = K_02*det;
5097
5098
// Compute element tensor using UFL quadrature representation
5099
// Optimisations: ('eliminate zeros', True), ('ignore ones', True), ('ignore zero tables', True), ('optimisation', 'simplify_expressions'), ('remove zero terms', True)
5100
5101
// Loop quadrature points for integral.
5102
// Number of operations to compute element tensor for following IP loop = 1032
5103
for (unsigned int ip = 0; ip < 4; ip++)
5104
{
5105
5106
// Coefficient declarations.
5107
double F0 = 0.0;
5108
double F1 = 0.0;
5109
double F2 = 0.0;
5110
5111
// Total number of operations to compute function values = 24
5112
for (unsigned int r = 0; r < 4; r++)
5113
{
5114
F0 += FE0_C0[ip][r]*w[0][nzc0[r]];
5115
F1 += FE0_C0[ip][r]*w[0][nzc4[r]];
5116
F2 += FE0_C0[ip][r]*w[0][nzc8[r]];
5117
}// end loop over 'r'
5118
5119
// Number of operations to compute ip constants: 18
5120
double I[3];
5121
// Number of operations: 6
5122
I[0] = W4[ip]*(F0*G[0] + F1*G[1] + F2*G[2]);
5123
5124
// Number of operations: 6
5125
I[1] = W4[ip]*(F0*G[3] + F1*G[4] + F2*G[5]);
5126
5127
// Number of operations: 6
5128
I[2] = W4[ip]*(F0*G[6] + F1*G[7] + F2*G[8]);
5129
5130
5131
// Number of operations for primary indices: 216
5132
for (unsigned int j = 0; j < 4; j++)
5133
{
5134
for (unsigned int k = 0; k < 2; k++)
5135
{
5136
// Number of operations to compute entry: 3
5137
A[nzc0[j]*12 + nzc1[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[0];
5138
// Number of operations to compute entry: 3
5139
A[nzc0[j]*12 + nzc2[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[1];
5140
// Number of operations to compute entry: 3
5141
A[nzc0[j]*12 + nzc3[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[2];
5142
// Number of operations to compute entry: 3
5143
A[nzc4[j]*12 + nzc5[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[0];
5144
// Number of operations to compute entry: 3
5145
A[nzc4[j]*12 + nzc6[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[1];
5146
// Number of operations to compute entry: 3
5147
A[nzc4[j]*12 + nzc7[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[2];
5148
// Number of operations to compute entry: 3
5149
A[nzc8[j]*12 + nzc10[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[1];
5150
// Number of operations to compute entry: 3
5151
A[nzc8[j]*12 + nzc11[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[2];
5152
// Number of operations to compute entry: 3
5153
A[nzc8[j]*12 + nzc9[k]] += FE0_C0[ip][j]*FE0_C0_D001[ip][k]*I[0];
5154
}// end loop over 'k'
5155
}// end loop over 'j'
5156
}// end loop over 'ip'
5157
}
5158
5159
/// Tabulate the tensor for the contribution from a local cell
5160
/// using the specified reference cell quadrature points/weights
5161
virtual void tabulate_tensor(double* A,
5162
const double * const * w,
5163
const ufc::cell& c,
5164
unsigned int num_quadrature_points,
5165
const double * const * quadrature_points,
5166
const double* quadrature_weights) const
5167
{
5168
std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
5169
}
5170
5171
};
5172
5173
/// This class defines the interface for the assembly of the global
5174
/// tensor corresponding to a form with r + n arguments, that is, a
5175
/// mapping
5176
///
5177
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
5178
///
5179
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
5180
/// global tensor A is defined by
5181
///
5182
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
5183
///
5184
/// where each argument Vj represents the application to the
5185
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
5186
/// fixed functions (coefficients).
5187
5188
class navierstokes_form_0: public ufc::form
5189
{
5190
public:
5191
5192
/// Constructor
5193
navierstokes_form_0() : ufc::form()
5194
{
5195
// Do nothing
5196
}
5197
5198
/// Destructor
5199
virtual ~navierstokes_form_0()
5200
{
5201
// Do nothing
5202
}
5203
5204
/// Return a string identifying the form
5205
virtual const char* signature() const
5206
{
5207
return "Form([Integral(IndexSum(Product(IndexSum(Product(Indexed(Coefficient(VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None), 0), MultiIndex((Index(0),), {Index(0): 3})), Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None), 1), MultiIndex((Index(0),), {Index(0): 3})), MultiIndex((Index(1),), {Index(1): 3}))), MultiIndex((Index(0),), {Index(0): 3})), Indexed(Argument(VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None), 0), MultiIndex((Index(1),), {Index(1): 3}))), MultiIndex((Index(1),), {Index(1): 3})), Measure('cell', 0, None))])";
5208
}
5209
5210
/// Return the rank of the global tensor (r)
5211
virtual unsigned int rank() const
5212
{
5213
return 2;
5214
}
5215
5216
/// Return the number of coefficients (n)
5217
virtual unsigned int num_coefficients() const
5218
{
5219
return 1;
5220
}
5221
5222
/// Return the number of cell domains
5223
virtual unsigned int num_cell_domains() const
5224
{
5225
return 1;
5226
}
5227
5228
/// Return the number of exterior facet domains
5229
virtual unsigned int num_exterior_facet_domains() const
5230
{
5231
return 0;
5232
}
5233
5234
/// Return the number of interior facet domains
5235
virtual unsigned int num_interior_facet_domains() const
5236
{
5237
return 0;
5238
}
5239
5240
/// Create a new finite element for argument function i
5241
virtual ufc::finite_element* create_finite_element(unsigned int i) const
5242
{
5243
switch (i)
5244
{
5245
case 0:
5246
{
5247
return new navierstokes_finite_element_1();
5248
break;
5249
}
5250
case 1:
5251
{
5252
return new navierstokes_finite_element_1();
5253
break;
5254
}
5255
case 2:
5256
{
5257
return new navierstokes_finite_element_1();
5258
break;
5259
}
5260
}
5261
5262
return 0;
5263
}
5264
5265
/// Create a new dofmap for argument function i
5266
virtual ufc::dofmap* create_dofmap(unsigned int i) const
5267
{
5268
switch (i)
5269
{
5270
case 0:
5271
{
5272
return new navierstokes_dofmap_1();
5273
break;
5274
}
5275
case 1:
5276
{
5277
return new navierstokes_dofmap_1();
5278
break;
5279
}
5280
case 2:
5281
{
5282
return new navierstokes_dofmap_1();
5283
break;
5284
}
5285
}
5286
5287
return 0;
5288
}
5289
5290
/// Create a new cell integral on sub domain i
5291
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
5292
{
5293
switch (i)
5294
{
5295
case 0:
5296
{
5297
return new navierstokes_cell_integral_0_0();
5298
break;
5299
}
5300
}
5301
5302
return 0;
5303
}
5304
5305
/// Create a new exterior facet integral on sub domain i
5306
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
5307
{
5308
return 0;
5309
}
5310
5311
/// Create a new interior facet integral on sub domain i
5312
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
5313
{
5314
return 0;
5315
}
5316
5317
};
5318
5319
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1