~ubuntu-branches/ubuntu/raring/ffc/raring
» Revision
4
: test/regression/references/StabilisedStokes.h
« back to all changes in this revision
Viewing changes to test/regression/references/StabilisedStokes.h
- browse files at revision 4
- compare with another revision
- download diff
- download tarball
- view history from revision 4
- Committer: Bazaar Package Importer
- Author(s): Johannes Ring
- Date: 2010-02-03 20:22:35 UTC
- mfrom: (1.1.2 upstream)
- Revision ID: james.westby@ubuntu.com-20100203202235-fe8d0kajuvgy2sqn
Tags: 0.9.0-1
* New upstream release.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.
- files added:
- demo/VectorPoisson.ufl
- ffc/quadrature
- ffc/tensor
- test/evaluate_basis
- test/evaluate_basis_derivatives
- test/regression/references
- test/unit/evaluate_basis
- test/unit/misc
- test/unit/symbolics
- files removed:
- ffc/common
- ffc/compiler
- ffc/compiler/quadrature
- ffc/compiler/tensor
- ffc/fem
- ffc/jit
- test/benchmarks
- test/regression/reference
- test/regression/reference/quadrature
- test/regression/reference/tensor
- test/simple_verify_tensors
- test/verify_tensors
- test/verify_tensors/references
- test/verify_tensors/test_cache
- files modified:
- AUTHORS
added
removed
test/regression/references/StabilisedStokes.h
1
// This code conforms with the UFC specification version 1.4
2
// and was automatically generated by FFC version 0.9.0.
3
4
#ifndef __STABILISEDSTOKES_H
5
#define __STABILISEDSTOKES_H
6
7
#include <cmath>
8
#include <stdexcept>
9
#include <fstream>
10
#include <ufc.h>
11
12
/// This class defines the interface for a finite element.
13
14
class stabilisedstokes_finite_element_0: public ufc::finite_element
15
{
16
public:
17
18
/// Constructor
19
stabilisedstokes_finite_element_0() : ufc::finite_element()
20
{
21
// Do nothing
22
}
23
24
/// Destructor
25
virtual ~stabilisedstokes_finite_element_0()
26
{
27
// Do nothing
28
}
29
30
/// Return a string identifying the finite element
31
virtual const char* signature() const
32
{
33
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
34
}
35
36
/// Return the cell shape
37
virtual ufc::shape cell_shape() const
38
{
39
return ufc::triangle;
40
}
41
42
/// Return the dimension of the finite element function space
43
virtual unsigned int space_dimension() const
44
{
45
return 3;
46
}
47
48
/// Return the rank of the value space
49
virtual unsigned int value_rank() const
50
{
51
return 0;
52
}
53
54
/// Return the dimension of the value space for axis i
55
virtual unsigned int value_dimension(unsigned int i) const
56
{
57
return 1;
58
}
59
60
/// Evaluate basis function i at given point in cell
61
virtual void evaluate_basis(unsigned int i,
62
double* values,
63
const double* coordinates,
64
const ufc::cell& c) const
65
{
66
// Extract vertex coordinates
67
const double * const * x = c.coordinates;
68
69
// Compute Jacobian of affine map from reference cell
70
const double J_00 = x[1][0] - x[0][0];
71
const double J_01 = x[2][0] - x[0][0];
72
const double J_10 = x[1][1] - x[0][1];
73
const double J_11 = x[2][1] - x[0][1];
74
75
// Compute determinant of Jacobian
76
double detJ = J_00*J_11 - J_01*J_10;
77
78
// Compute inverse of Jacobian
79
80
// Compute constants
81
const double C0 = x[1][0] + x[2][0];
82
const double C1 = x[1][1] + x[2][1];
83
84
// Get coordinates and map to the reference (FIAT) element
85
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
86
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
87
88
// Reset values
89
*values = 0.00000000;
90
91
// Map degree of freedom to element degree of freedom
92
const unsigned int dof = i;
93
94
// Array of basisvalues
95
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
96
97
// Declare helper variables
98
unsigned int rr = 0;
99
unsigned int ss = 0;
100
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
101
102
// Compute basisvalues
103
basisvalues[0] = 1.00000000;
104
basisvalues[1] = tmp0;
105
for (unsigned int r = 0; r < 1; r++)
106
{
107
rr = (r + 1)*(r + 1 + 1)/2 + 1;
108
ss = r*(r + 1)/2;
109
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
110
}// end loop over 'r'
111
for (unsigned int r = 0; r < 2; r++)
112
{
113
for (unsigned int s = 0; s < 2 - r; s++)
114
{
115
rr = (r + s)*(r + s + 1)/2 + s;
116
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
117
}// end loop over 's'
118
}// end loop over 'r'
119
120
// Table(s) of coefficients
121
static const double coefficients0[3][3] = \
122
{{0.47140452, -0.28867513, -0.16666667},
123
{0.47140452, 0.28867513, -0.16666667},
124
{0.47140452, 0.00000000, 0.33333333}};
125
126
// Compute value(s).
127
for (unsigned int r = 0; r < 3; r++)
128
{
129
*values += coefficients0[dof][r]*basisvalues[r];
130
}// end loop over 'r'
131
}
132
133
/// Evaluate all basis functions at given point in cell
134
virtual void evaluate_basis_all(double* values,
135
const double* coordinates,
136
const ufc::cell& c) const
137
{
138
// Helper variable to hold values of a single dof.
139
double dof_values = 0.00000000;
140
141
// Loop dofs and call evaluate_basis.
142
for (unsigned int r = 0; r < 3; r++)
143
{
144
evaluate_basis(r, &dof_values, coordinates, c);
145
values[r] = dof_values;
146
}// end loop over 'r'
147
}
148
149
/// Evaluate order n derivatives of basis function i at given point in cell
150
virtual void evaluate_basis_derivatives(unsigned int i,
151
unsigned int n,
152
double* values,
153
const double* coordinates,
154
const ufc::cell& c) const
155
{
156
// Extract vertex coordinates
157
const double * const * x = c.coordinates;
158
159
// Compute Jacobian of affine map from reference cell
160
const double J_00 = x[1][0] - x[0][0];
161
const double J_01 = x[2][0] - x[0][0];
162
const double J_10 = x[1][1] - x[0][1];
163
const double J_11 = x[2][1] - x[0][1];
164
165
// Compute determinant of Jacobian
166
double detJ = J_00*J_11 - J_01*J_10;
167
168
// Compute inverse of Jacobian
169
const double K_00 = J_11 / detJ;
170
const double K_01 = -J_01 / detJ;
171
const double K_10 = -J_10 / detJ;
172
const double K_11 = J_00 / detJ;
173
174
// Compute constants
175
const double C0 = x[1][0] + x[2][0];
176
const double C1 = x[1][1] + x[2][1];
177
178
// Get coordinates and map to the reference (FIAT) element
179
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
180
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
181
182
// Compute number of derivatives.
183
unsigned int num_derivatives = 1;
184
185
for (unsigned int r = 0; r < n; r++)
186
{
187
num_derivatives *= 2;
188
}// end loop over 'r'
189
190
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
191
unsigned int **combinations = new unsigned int *[num_derivatives];
192
for (unsigned int row = 0; row < num_derivatives; row++)
193
{
194
combinations[row] = new unsigned int [n];
195
for (unsigned int col = 0; col < n; col++)
196
combinations[row][col] = 0;
197
}
198
199
// Generate combinations of derivatives
200
for (unsigned int row = 1; row < num_derivatives; row++)
201
{
202
for (unsigned int num = 0; num < row; num++)
203
{
204
for (unsigned int col = n-1; col+1 > 0; col--)
205
{
206
if (combinations[row][col] + 1 > 1)
207
combinations[row][col] = 0;
208
else
209
{
210
combinations[row][col] += 1;
211
break;
212
}
213
}
214
}
215
}
216
217
// Compute inverse of Jacobian
218
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
219
220
// Declare transformation matrix
221
// Declare pointer to two dimensional array and initialise
222
double **transform = new double *[num_derivatives];
223
224
for (unsigned int j = 0; j < num_derivatives; j++)
225
{
226
transform[j] = new double [num_derivatives];
227
for (unsigned int k = 0; k < num_derivatives; k++)
228
transform[j][k] = 1;
229
}
230
231
// Construct transformation matrix
232
for (unsigned int row = 0; row < num_derivatives; row++)
233
{
234
for (unsigned int col = 0; col < num_derivatives; col++)
235
{
236
for (unsigned int k = 0; k < n; k++)
237
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
238
}
239
}
240
241
// Reset values. Assuming that values is always an array.
242
for (unsigned int r = 0; r < num_derivatives; r++)
243
{
244
values[r] = 0.00000000;
245
}// end loop over 'r'
246
247
// Map degree of freedom to element degree of freedom
248
const unsigned int dof = i;
249
250
// Array of basisvalues
251
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
252
253
// Declare helper variables
254
unsigned int rr = 0;
255
unsigned int ss = 0;
256
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
257
258
// Compute basisvalues
259
basisvalues[0] = 1.00000000;
260
basisvalues[1] = tmp0;
261
for (unsigned int r = 0; r < 1; r++)
262
{
263
rr = (r + 1)*(r + 1 + 1)/2 + 1;
264
ss = r*(r + 1)/2;
265
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
266
}// end loop over 'r'
267
for (unsigned int r = 0; r < 2; r++)
268
{
269
for (unsigned int s = 0; s < 2 - r; s++)
270
{
271
rr = (r + s)*(r + s + 1)/2 + s;
272
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
273
}// end loop over 's'
274
}// end loop over 'r'
275
276
// Table(s) of coefficients
277
static const double coefficients0[3][3] = \
278
{{0.47140452, -0.28867513, -0.16666667},
279
{0.47140452, 0.28867513, -0.16666667},
280
{0.47140452, 0.00000000, 0.33333333}};
281
282
// Tables of derivatives of the polynomial base (transpose).
283
static const double dmats0[3][3] = \
284
{{0.00000000, 0.00000000, 0.00000000},
285
{4.89897949, 0.00000000, 0.00000000},
286
{0.00000000, 0.00000000, 0.00000000}};
287
288
static const double dmats1[3][3] = \
289
{{0.00000000, 0.00000000, 0.00000000},
290
{2.44948974, 0.00000000, 0.00000000},
291
{4.24264069, 0.00000000, 0.00000000}};
292
293
// Compute reference derivatives
294
// Declare pointer to array of derivatives on FIAT element
295
double *derivatives = new double [num_derivatives];
296
for (unsigned int r = 0; r < num_derivatives; r++)
297
{
298
derivatives[r] = 0.00000000;
299
}// end loop over 'r'
300
301
// Declare derivative matrix (of polynomial basis).
302
double dmats[3][3] = \
303
{{1.00000000, 0.00000000, 0.00000000},
304
{0.00000000, 1.00000000, 0.00000000},
305
{0.00000000, 0.00000000, 1.00000000}};
306
307
// Declare (auxiliary) derivative matrix (of polynomial basis).
308
double dmats_old[3][3] = \
309
{{1.00000000, 0.00000000, 0.00000000},
310
{0.00000000, 1.00000000, 0.00000000},
311
{0.00000000, 0.00000000, 1.00000000}};
312
313
// Loop possible derivatives.
314
for (unsigned int r = 0; r < num_derivatives; r++)
315
{
316
// Resetting dmats values to compute next derivative.
317
for (unsigned int t = 0; t < 3; t++)
318
{
319
for (unsigned int u = 0; u < 3; u++)
320
{
321
dmats[t][u] = 0.00000000;
322
if (t == u)
323
{
324
dmats[t][u] = 1.00000000;
325
}
326
327
}// end loop over 'u'
328
}// end loop over 't'
329
330
// Looping derivative order to generate dmats.
331
for (unsigned int s = 0; s < n; s++)
332
{
333
// Updating dmats_old with new values and resetting dmats.
334
for (unsigned int t = 0; t < 3; t++)
335
{
336
for (unsigned int u = 0; u < 3; u++)
337
{
338
dmats_old[t][u] = dmats[t][u];
339
dmats[t][u] = 0.00000000;
340
}// end loop over 'u'
341
}// end loop over 't'
342
343
// Update dmats using an inner product.
344
if (combinations[r][s] == 0)
345
{
346
for (unsigned int t = 0; t < 3; t++)
347
{
348
for (unsigned int u = 0; u < 3; u++)
349
{
350
for (unsigned int tu = 0; tu < 3; tu++)
351
{
352
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
353
}// end loop over 'tu'
354
}// end loop over 'u'
355
}// end loop over 't'
356
}
357
358
if (combinations[r][s] == 1)
359
{
360
for (unsigned int t = 0; t < 3; t++)
361
{
362
for (unsigned int u = 0; u < 3; u++)
363
{
364
for (unsigned int tu = 0; tu < 3; tu++)
365
{
366
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
367
}// end loop over 'tu'
368
}// end loop over 'u'
369
}// end loop over 't'
370
}
371
372
}// end loop over 's'
373
for (unsigned int s = 0; s < 3; s++)
374
{
375
for (unsigned int t = 0; t < 3; t++)
376
{
377
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
378
}// end loop over 't'
379
}// end loop over 's'
380
}// end loop over 'r'
381
382
// Transform derivatives back to physical element
383
for (unsigned int row = 0; row < num_derivatives; row++)
384
{
385
for (unsigned int col = 0; col < num_derivatives; col++)
386
{
387
values[row] += transform[row][col]*derivatives[col];
388
}
389
}
390
391
// Delete pointer to array of derivatives on FIAT element
392
delete [] derivatives;
393
394
// Delete pointer to array of combinations of derivatives and transform
395
for (unsigned int r = 0; r < num_derivatives; r++)
396
{
397
delete [] combinations[r];
398
delete [] transform[r];
399
}// end loop over 'r'
400
delete [] combinations;
401
delete [] transform;
402
}
403
404
/// Evaluate order n derivatives of all basis functions at given point in cell
405
virtual void evaluate_basis_derivatives_all(unsigned int n,
406
double* values,
407
const double* coordinates,
408
const ufc::cell& c) const
409
{
410
// Compute number of derivatives.
411
unsigned int num_derivatives = 1;
412
413
for (unsigned int r = 0; r < n; r++)
414
{
415
num_derivatives *= 2;
416
}// end loop over 'r'
417
418
// Helper variable to hold values of a single dof.
419
double *dof_values = new double [num_derivatives];
420
for (unsigned int r = 0; r < num_derivatives; r++)
421
{
422
dof_values[r] = 0.00000000;
423
}// end loop over 'r'
424
425
// Loop dofs and call evaluate_basis_derivatives.
426
for (unsigned int r = 0; r < 3; r++)
427
{
428
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
429
for (unsigned int s = 0; s < num_derivatives; s++)
430
{
431
values[r*num_derivatives + s] = dof_values[s];
432
}// end loop over 's'
433
}// end loop over 'r'
434
435
// Delete pointer.
436
delete [] dof_values;
437
}
438
439
/// Evaluate linear functional for dof i on the function f
440
virtual double evaluate_dof(unsigned int i,
441
const ufc::function& f,
442
const ufc::cell& c) const
443
{
444
// Declare variables for result of evaluation
445
double vals[1];
446
447
// Declare variable for physical coordinates
448
double y[2];
449
450
const double * const * x = c.coordinates;
451
switch (i)
452
{
453
case 0:
454
{
455
y[0] = x[0][0];
456
y[1] = x[0][1];
457
f.evaluate(vals, y, c);
458
return vals[0];
459
break;
460
}
461
case 1:
462
{
463
y[0] = x[1][0];
464
y[1] = x[1][1];
465
f.evaluate(vals, y, c);
466
return vals[0];
467
break;
468
}
469
case 2:
470
{
471
y[0] = x[2][0];
472
y[1] = x[2][1];
473
f.evaluate(vals, y, c);
474
return vals[0];
475
break;
476
}
477
}
478
479
return 0.0;
480
}
481
482
/// Evaluate linear functionals for all dofs on the function f
483
virtual void evaluate_dofs(double* values,
484
const ufc::function& f,
485
const ufc::cell& c) const
486
{
487
// Declare variables for result of evaluation
488
double vals[1];
489
490
// Declare variable for physical coordinates
491
double y[2];
492
493
const double * const * x = c.coordinates;
494
y[0] = x[0][0];
495
y[1] = x[0][1];
496
f.evaluate(vals, y, c);
497
values[0] = vals[0];
498
y[0] = x[1][0];
499
y[1] = x[1][1];
500
f.evaluate(vals, y, c);
501
values[1] = vals[0];
502
y[0] = x[2][0];
503
y[1] = x[2][1];
504
f.evaluate(vals, y, c);
505
values[2] = vals[0];
506
}
507
508
/// Interpolate vertex values from dof values
509
virtual void interpolate_vertex_values(double* vertex_values,
510
const double* dof_values,
511
const ufc::cell& c) const
512
{
513
// Evaluate function and change variables
514
vertex_values[0] = dof_values[0];
515
vertex_values[1] = dof_values[1];
516
vertex_values[2] = dof_values[2];
517
}
518
519
/// Return the number of sub elements (for a mixed element)
520
virtual unsigned int num_sub_elements() const
521
{
522
return 0;
523
}
524
525
/// Create a new finite element for sub element i (for a mixed element)
526
virtual ufc::finite_element* create_sub_element(unsigned int i) const
527
{
528
return 0;
529
}
530
531
};
532
533
/// This class defines the interface for a finite element.
534
535
class stabilisedstokes_finite_element_1: public ufc::finite_element
536
{
537
public:
538
539
/// Constructor
540
stabilisedstokes_finite_element_1() : ufc::finite_element()
541
{
542
// Do nothing
543
}
544
545
/// Destructor
546
virtual ~stabilisedstokes_finite_element_1()
547
{
548
// Do nothing
549
}
550
551
/// Return a string identifying the finite element
552
virtual const char* signature() const
553
{
554
return "VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2)";
555
}
556
557
/// Return the cell shape
558
virtual ufc::shape cell_shape() const
559
{
560
return ufc::triangle;
561
}
562
563
/// Return the dimension of the finite element function space
564
virtual unsigned int space_dimension() const
565
{
566
return 6;
567
}
568
569
/// Return the rank of the value space
570
virtual unsigned int value_rank() const
571
{
572
return 1;
573
}
574
575
/// Return the dimension of the value space for axis i
576
virtual unsigned int value_dimension(unsigned int i) const
577
{
578
switch (i)
579
{
580
case 0:
581
{
582
return 2;
583
break;
584
}
585
}
586
587
return 0;
588
}
589
590
/// Evaluate basis function i at given point in cell
591
virtual void evaluate_basis(unsigned int i,
592
double* values,
593
const double* coordinates,
594
const ufc::cell& c) const
595
{
596
// Extract vertex coordinates
597
const double * const * x = c.coordinates;
598
599
// Compute Jacobian of affine map from reference cell
600
const double J_00 = x[1][0] - x[0][0];
601
const double J_01 = x[2][0] - x[0][0];
602
const double J_10 = x[1][1] - x[0][1];
603
const double J_11 = x[2][1] - x[0][1];
604
605
// Compute determinant of Jacobian
606
double detJ = J_00*J_11 - J_01*J_10;
607
608
// Compute inverse of Jacobian
609
610
// Compute constants
611
const double C0 = x[1][0] + x[2][0];
612
const double C1 = x[1][1] + x[2][1];
613
614
// Get coordinates and map to the reference (FIAT) element
615
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
616
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
617
618
// Reset values
619
values[0] = 0.00000000;
620
values[1] = 0.00000000;
621
if (0 <= i && i <= 2)
622
{
623
// Map degree of freedom to element degree of freedom
624
const unsigned int dof = i;
625
626
// Array of basisvalues
627
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
628
629
// Declare helper variables
630
unsigned int rr = 0;
631
unsigned int ss = 0;
632
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
633
634
// Compute basisvalues
635
basisvalues[0] = 1.00000000;
636
basisvalues[1] = tmp0;
637
for (unsigned int r = 0; r < 1; r++)
638
{
639
rr = (r + 1)*(r + 1 + 1)/2 + 1;
640
ss = r*(r + 1)/2;
641
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
642
}// end loop over 'r'
643
for (unsigned int r = 0; r < 2; r++)
644
{
645
for (unsigned int s = 0; s < 2 - r; s++)
646
{
647
rr = (r + s)*(r + s + 1)/2 + s;
648
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
649
}// end loop over 's'
650
}// end loop over 'r'
651
652
// Table(s) of coefficients
653
static const double coefficients0[3][3] = \
654
{{0.47140452, -0.28867513, -0.16666667},
655
{0.47140452, 0.28867513, -0.16666667},
656
{0.47140452, 0.00000000, 0.33333333}};
657
658
// Compute value(s).
659
for (unsigned int r = 0; r < 3; r++)
660
{
661
values[0] += coefficients0[dof][r]*basisvalues[r];
662
}// end loop over 'r'
663
}
664
665
if (3 <= i && i <= 5)
666
{
667
// Map degree of freedom to element degree of freedom
668
const unsigned int dof = i - 3;
669
670
// Array of basisvalues
671
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
672
673
// Declare helper variables
674
unsigned int rr = 0;
675
unsigned int ss = 0;
676
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
677
678
// Compute basisvalues
679
basisvalues[0] = 1.00000000;
680
basisvalues[1] = tmp0;
681
for (unsigned int r = 0; r < 1; r++)
682
{
683
rr = (r + 1)*(r + 1 + 1)/2 + 1;
684
ss = r*(r + 1)/2;
685
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
686
}// end loop over 'r'
687
for (unsigned int r = 0; r < 2; r++)
688
{
689
for (unsigned int s = 0; s < 2 - r; s++)
690
{
691
rr = (r + s)*(r + s + 1)/2 + s;
692
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
693
}// end loop over 's'
694
}// end loop over 'r'
695
696
// Table(s) of coefficients
697
static const double coefficients0[3][3] = \
698
{{0.47140452, -0.28867513, -0.16666667},
699
{0.47140452, 0.28867513, -0.16666667},
700
{0.47140452, 0.00000000, 0.33333333}};
701
702
// Compute value(s).
703
for (unsigned int r = 0; r < 3; r++)
704
{
705
values[1] += coefficients0[dof][r]*basisvalues[r];
706
}// end loop over 'r'
707
}
708
709
}
710
711
/// Evaluate all basis functions at given point in cell
712
virtual void evaluate_basis_all(double* values,
713
const double* coordinates,
714
const ufc::cell& c) const
715
{
716
// Helper variable to hold values of a single dof.
717
double dof_values[2] = {0.00000000, 0.00000000};
718
719
// Loop dofs and call evaluate_basis.
720
for (unsigned int r = 0; r < 6; r++)
721
{
722
evaluate_basis(r, dof_values, coordinates, c);
723
for (unsigned int s = 0; s < 2; s++)
724
{
725
values[r*2 + s] = dof_values[s];
726
}// end loop over 's'
727
}// end loop over 'r'
728
}
729
730
/// Evaluate order n derivatives of basis function i at given point in cell
731
virtual void evaluate_basis_derivatives(unsigned int i,
732
unsigned int n,
733
double* values,
734
const double* coordinates,
735
const ufc::cell& c) const
736
{
737
// Extract vertex coordinates
738
const double * const * x = c.coordinates;
739
740
// Compute Jacobian of affine map from reference cell
741
const double J_00 = x[1][0] - x[0][0];
742
const double J_01 = x[2][0] - x[0][0];
743
const double J_10 = x[1][1] - x[0][1];
744
const double J_11 = x[2][1] - x[0][1];
745
746
// Compute determinant of Jacobian
747
double detJ = J_00*J_11 - J_01*J_10;
748
749
// Compute inverse of Jacobian
750
const double K_00 = J_11 / detJ;
751
const double K_01 = -J_01 / detJ;
752
const double K_10 = -J_10 / detJ;
753
const double K_11 = J_00 / detJ;
754
755
// Compute constants
756
const double C0 = x[1][0] + x[2][0];
757
const double C1 = x[1][1] + x[2][1];
758
759
// Get coordinates and map to the reference (FIAT) element
760
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
761
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
762
763
// Compute number of derivatives.
764
unsigned int num_derivatives = 1;
765
766
for (unsigned int r = 0; r < n; r++)
767
{
768
num_derivatives *= 2;
769
}// end loop over 'r'
770
771
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
772
unsigned int **combinations = new unsigned int *[num_derivatives];
773
for (unsigned int row = 0; row < num_derivatives; row++)
774
{
775
combinations[row] = new unsigned int [n];
776
for (unsigned int col = 0; col < n; col++)
777
combinations[row][col] = 0;
778
}
779
780
// Generate combinations of derivatives
781
for (unsigned int row = 1; row < num_derivatives; row++)
782
{
783
for (unsigned int num = 0; num < row; num++)
784
{
785
for (unsigned int col = n-1; col+1 > 0; col--)
786
{
787
if (combinations[row][col] + 1 > 1)
788
combinations[row][col] = 0;
789
else
790
{
791
combinations[row][col] += 1;
792
break;
793
}
794
}
795
}
796
}
797
798
// Compute inverse of Jacobian
799
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
800
801
// Declare transformation matrix
802
// Declare pointer to two dimensional array and initialise
803
double **transform = new double *[num_derivatives];
804
805
for (unsigned int j = 0; j < num_derivatives; j++)
806
{
807
transform[j] = new double [num_derivatives];
808
for (unsigned int k = 0; k < num_derivatives; k++)
809
transform[j][k] = 1;
810
}
811
812
// Construct transformation matrix
813
for (unsigned int row = 0; row < num_derivatives; row++)
814
{
815
for (unsigned int col = 0; col < num_derivatives; col++)
816
{
817
for (unsigned int k = 0; k < n; k++)
818
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
819
}
820
}
821
822
// Reset values. Assuming that values is always an array.
823
for (unsigned int r = 0; r < 2*num_derivatives; r++)
824
{
825
values[r] = 0.00000000;
826
}// end loop over 'r'
827
828
if (0 <= i && i <= 2)
829
{
830
// Map degree of freedom to element degree of freedom
831
const unsigned int dof = i;
832
833
// Array of basisvalues
834
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
835
836
// Declare helper variables
837
unsigned int rr = 0;
838
unsigned int ss = 0;
839
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
840
841
// Compute basisvalues
842
basisvalues[0] = 1.00000000;
843
basisvalues[1] = tmp0;
844
for (unsigned int r = 0; r < 1; r++)
845
{
846
rr = (r + 1)*(r + 1 + 1)/2 + 1;
847
ss = r*(r + 1)/2;
848
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
849
}// end loop over 'r'
850
for (unsigned int r = 0; r < 2; r++)
851
{
852
for (unsigned int s = 0; s < 2 - r; s++)
853
{
854
rr = (r + s)*(r + s + 1)/2 + s;
855
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
856
}// end loop over 's'
857
}// end loop over 'r'
858
859
// Table(s) of coefficients
860
static const double coefficients0[3][3] = \
861
{{0.47140452, -0.28867513, -0.16666667},
862
{0.47140452, 0.28867513, -0.16666667},
863
{0.47140452, 0.00000000, 0.33333333}};
864
865
// Tables of derivatives of the polynomial base (transpose).
866
static const double dmats0[3][3] = \
867
{{0.00000000, 0.00000000, 0.00000000},
868
{4.89897949, 0.00000000, 0.00000000},
869
{0.00000000, 0.00000000, 0.00000000}};
870
871
static const double dmats1[3][3] = \
872
{{0.00000000, 0.00000000, 0.00000000},
873
{2.44948974, 0.00000000, 0.00000000},
874
{4.24264069, 0.00000000, 0.00000000}};
875
876
// Compute reference derivatives
877
// Declare pointer to array of derivatives on FIAT element
878
double *derivatives = new double [num_derivatives];
879
for (unsigned int r = 0; r < num_derivatives; r++)
880
{
881
derivatives[r] = 0.00000000;
882
}// end loop over 'r'
883
884
// Declare derivative matrix (of polynomial basis).
885
double dmats[3][3] = \
886
{{1.00000000, 0.00000000, 0.00000000},
887
{0.00000000, 1.00000000, 0.00000000},
888
{0.00000000, 0.00000000, 1.00000000}};
889
890
// Declare (auxiliary) derivative matrix (of polynomial basis).
891
double dmats_old[3][3] = \
892
{{1.00000000, 0.00000000, 0.00000000},
893
{0.00000000, 1.00000000, 0.00000000},
894
{0.00000000, 0.00000000, 1.00000000}};
895
896
// Loop possible derivatives.
897
for (unsigned int r = 0; r < num_derivatives; r++)
898
{
899
// Resetting dmats values to compute next derivative.
900
for (unsigned int t = 0; t < 3; t++)
901
{
902
for (unsigned int u = 0; u < 3; u++)
903
{
904
dmats[t][u] = 0.00000000;
905
if (t == u)
906
{
907
dmats[t][u] = 1.00000000;
908
}
909
910
}// end loop over 'u'
911
}// end loop over 't'
912
913
// Looping derivative order to generate dmats.
914
for (unsigned int s = 0; s < n; s++)
915
{
916
// Updating dmats_old with new values and resetting dmats.
917
for (unsigned int t = 0; t < 3; t++)
918
{
919
for (unsigned int u = 0; u < 3; u++)
920
{
921
dmats_old[t][u] = dmats[t][u];
922
dmats[t][u] = 0.00000000;
923
}// end loop over 'u'
924
}// end loop over 't'
925
926
// Update dmats using an inner product.
927
if (combinations[r][s] == 0)
928
{
929
for (unsigned int t = 0; t < 3; t++)
930
{
931
for (unsigned int u = 0; u < 3; u++)
932
{
933
for (unsigned int tu = 0; tu < 3; tu++)
934
{
935
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
936
}// end loop over 'tu'
937
}// end loop over 'u'
938
}// end loop over 't'
939
}
940
941
if (combinations[r][s] == 1)
942
{
943
for (unsigned int t = 0; t < 3; t++)
944
{
945
for (unsigned int u = 0; u < 3; u++)
946
{
947
for (unsigned int tu = 0; tu < 3; tu++)
948
{
949
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
950
}// end loop over 'tu'
951
}// end loop over 'u'
952
}// end loop over 't'
953
}
954
955
}// end loop over 's'
956
for (unsigned int s = 0; s < 3; s++)
957
{
958
for (unsigned int t = 0; t < 3; t++)
959
{
960
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
961
}// end loop over 't'
962
}// end loop over 's'
963
}// end loop over 'r'
964
965
// Transform derivatives back to physical element
966
for (unsigned int row = 0; row < num_derivatives; row++)
967
{
968
for (unsigned int col = 0; col < num_derivatives; col++)
969
{
970
values[row] += transform[row][col]*derivatives[col];
971
}
972
}
973
974
// Delete pointer to array of derivatives on FIAT element
975
delete [] derivatives;
976
977
// Delete pointer to array of combinations of derivatives and transform
978
for (unsigned int r = 0; r < num_derivatives; r++)
979
{
980
delete [] combinations[r];
981
delete [] transform[r];
982
}// end loop over 'r'
983
delete [] combinations;
984
delete [] transform;
985
}
986
987
if (3 <= i && i <= 5)
988
{
989
// Map degree of freedom to element degree of freedom
990
const unsigned int dof = i - 3;
991
992
// Array of basisvalues
993
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
994
995
// Declare helper variables
996
unsigned int rr = 0;
997
unsigned int ss = 0;
998
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
999
1000
// Compute basisvalues
1001
basisvalues[0] = 1.00000000;
1002
basisvalues[1] = tmp0;
1003
for (unsigned int r = 0; r < 1; r++)
1004
{
1005
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1006
ss = r*(r + 1)/2;
1007
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1008
}// end loop over 'r'
1009
for (unsigned int r = 0; r < 2; r++)
1010
{
1011
for (unsigned int s = 0; s < 2 - r; s++)
1012
{
1013
rr = (r + s)*(r + s + 1)/2 + s;
1014
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1015
}// end loop over 's'
1016
}// end loop over 'r'
1017
1018
// Table(s) of coefficients
1019
static const double coefficients0[3][3] = \
1020
{{0.47140452, -0.28867513, -0.16666667},
1021
{0.47140452, 0.28867513, -0.16666667},
1022
{0.47140452, 0.00000000, 0.33333333}};
1023
1024
// Tables of derivatives of the polynomial base (transpose).
1025
static const double dmats0[3][3] = \
1026
{{0.00000000, 0.00000000, 0.00000000},
1027
{4.89897949, 0.00000000, 0.00000000},
1028
{0.00000000, 0.00000000, 0.00000000}};
1029
1030
static const double dmats1[3][3] = \
1031
{{0.00000000, 0.00000000, 0.00000000},
1032
{2.44948974, 0.00000000, 0.00000000},
1033
{4.24264069, 0.00000000, 0.00000000}};
1034
1035
// Compute reference derivatives
1036
// Declare pointer to array of derivatives on FIAT element
1037
double *derivatives = new double [num_derivatives];
1038
for (unsigned int r = 0; r < num_derivatives; r++)
1039
{
1040
derivatives[r] = 0.00000000;
1041
}// end loop over 'r'
1042
1043
// Declare derivative matrix (of polynomial basis).
1044
double dmats[3][3] = \
1045
{{1.00000000, 0.00000000, 0.00000000},
1046
{0.00000000, 1.00000000, 0.00000000},
1047
{0.00000000, 0.00000000, 1.00000000}};
1048
1049
// Declare (auxiliary) derivative matrix (of polynomial basis).
1050
double dmats_old[3][3] = \
1051
{{1.00000000, 0.00000000, 0.00000000},
1052
{0.00000000, 1.00000000, 0.00000000},
1053
{0.00000000, 0.00000000, 1.00000000}};
1054
1055
// Loop possible derivatives.
1056
for (unsigned int r = 0; r < num_derivatives; r++)
1057
{
1058
// Resetting dmats values to compute next derivative.
1059
for (unsigned int t = 0; t < 3; t++)
1060
{
1061
for (unsigned int u = 0; u < 3; u++)
1062
{
1063
dmats[t][u] = 0.00000000;
1064
if (t == u)
1065
{
1066
dmats[t][u] = 1.00000000;
1067
}
1068
1069
}// end loop over 'u'
1070
}// end loop over 't'
1071
1072
// Looping derivative order to generate dmats.
1073
for (unsigned int s = 0; s < n; s++)
1074
{
1075
// Updating dmats_old with new values and resetting dmats.
1076
for (unsigned int t = 0; t < 3; t++)
1077
{
1078
for (unsigned int u = 0; u < 3; u++)
1079
{
1080
dmats_old[t][u] = dmats[t][u];
1081
dmats[t][u] = 0.00000000;
1082
}// end loop over 'u'
1083
}// end loop over 't'
1084
1085
// Update dmats using an inner product.
1086
if (combinations[r][s] == 0)
1087
{
1088
for (unsigned int t = 0; t < 3; t++)
1089
{
1090
for (unsigned int u = 0; u < 3; u++)
1091
{
1092
for (unsigned int tu = 0; tu < 3; tu++)
1093
{
1094
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1095
}// end loop over 'tu'
1096
}// end loop over 'u'
1097
}// end loop over 't'
1098
}
1099
1100
if (combinations[r][s] == 1)
1101
{
1102
for (unsigned int t = 0; t < 3; t++)
1103
{
1104
for (unsigned int u = 0; u < 3; u++)
1105
{
1106
for (unsigned int tu = 0; tu < 3; tu++)
1107
{
1108
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1109
}// end loop over 'tu'
1110
}// end loop over 'u'
1111
}// end loop over 't'
1112
}
1113
1114
}// end loop over 's'
1115
for (unsigned int s = 0; s < 3; s++)
1116
{
1117
for (unsigned int t = 0; t < 3; t++)
1118
{
1119
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1120
}// end loop over 't'
1121
}// end loop over 's'
1122
}// end loop over 'r'
1123
1124
// Transform derivatives back to physical element
1125
for (unsigned int row = 0; row < num_derivatives; row++)
1126
{
1127
for (unsigned int col = 0; col < num_derivatives; col++)
1128
{
1129
values[num_derivatives + row] += transform[row][col]*derivatives[col];
1130
}
1131
}
1132
1133
// Delete pointer to array of derivatives on FIAT element
1134
delete [] derivatives;
1135
1136
// Delete pointer to array of combinations of derivatives and transform
1137
for (unsigned int r = 0; r < num_derivatives; r++)
1138
{
1139
delete [] combinations[r];
1140
delete [] transform[r];
1141
}// end loop over 'r'
1142
delete [] combinations;
1143
delete [] transform;
1144
}
1145
1146
}
1147
1148
/// Evaluate order n derivatives of all basis functions at given point in cell
1149
virtual void evaluate_basis_derivatives_all(unsigned int n,
1150
double* values,
1151
const double* coordinates,
1152
const ufc::cell& c) const
1153
{
1154
// Compute number of derivatives.
1155
unsigned int num_derivatives = 1;
1156
1157
for (unsigned int r = 0; r < n; r++)
1158
{
1159
num_derivatives *= 2;
1160
}// end loop over 'r'
1161
1162
// Helper variable to hold values of a single dof.
1163
double *dof_values = new double [2*num_derivatives];
1164
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1165
{
1166
dof_values[r] = 0.00000000;
1167
}// end loop over 'r'
1168
1169
// Loop dofs and call evaluate_basis_derivatives.
1170
for (unsigned int r = 0; r < 6; r++)
1171
{
1172
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
1173
for (unsigned int s = 0; s < 2*num_derivatives; s++)
1174
{
1175
values[r*2*num_derivatives + s] = dof_values[s];
1176
}// end loop over 's'
1177
}// end loop over 'r'
1178
1179
// Delete pointer.
1180
delete [] dof_values;
1181
}
1182
1183
/// Evaluate linear functional for dof i on the function f
1184
virtual double evaluate_dof(unsigned int i,
1185
const ufc::function& f,
1186
const ufc::cell& c) const
1187
{
1188
// Declare variables for result of evaluation
1189
double vals[2];
1190
1191
// Declare variable for physical coordinates
1192
double y[2];
1193
1194
const double * const * x = c.coordinates;
1195
switch (i)
1196
{
1197
case 0:
1198
{
1199
y[0] = x[0][0];
1200
y[1] = x[0][1];
1201
f.evaluate(vals, y, c);
1202
return vals[0];
1203
break;
1204
}
1205
case 1:
1206
{
1207
y[0] = x[1][0];
1208
y[1] = x[1][1];
1209
f.evaluate(vals, y, c);
1210
return vals[0];
1211
break;
1212
}
1213
case 2:
1214
{
1215
y[0] = x[2][0];
1216
y[1] = x[2][1];
1217
f.evaluate(vals, y, c);
1218
return vals[0];
1219
break;
1220
}
1221
case 3:
1222
{
1223
y[0] = x[0][0];
1224
y[1] = x[0][1];
1225
f.evaluate(vals, y, c);
1226
return vals[1];
1227
break;
1228
}
1229
case 4:
1230
{
1231
y[0] = x[1][0];
1232
y[1] = x[1][1];
1233
f.evaluate(vals, y, c);
1234
return vals[1];
1235
break;
1236
}
1237
case 5:
1238
{
1239
y[0] = x[2][0];
1240
y[1] = x[2][1];
1241
f.evaluate(vals, y, c);
1242
return vals[1];
1243
break;
1244
}
1245
}
1246
1247
return 0.0;
1248
}
1249
1250
/// Evaluate linear functionals for all dofs on the function f
1251
virtual void evaluate_dofs(double* values,
1252
const ufc::function& f,
1253
const ufc::cell& c) const
1254
{
1255
// Declare variables for result of evaluation
1256
double vals[2];
1257
1258
// Declare variable for physical coordinates
1259
double y[2];
1260
1261
const double * const * x = c.coordinates;
1262
y[0] = x[0][0];
1263
y[1] = x[0][1];
1264
f.evaluate(vals, y, c);
1265
values[0] = vals[0];
1266
y[0] = x[1][0];
1267
y[1] = x[1][1];
1268
f.evaluate(vals, y, c);
1269
values[1] = vals[0];
1270
y[0] = x[2][0];
1271
y[1] = x[2][1];
1272
f.evaluate(vals, y, c);
1273
values[2] = vals[0];
1274
y[0] = x[0][0];
1275
y[1] = x[0][1];
1276
f.evaluate(vals, y, c);
1277
values[3] = vals[1];
1278
y[0] = x[1][0];
1279
y[1] = x[1][1];
1280
f.evaluate(vals, y, c);
1281
values[4] = vals[1];
1282
y[0] = x[2][0];
1283
y[1] = x[2][1];
1284
f.evaluate(vals, y, c);
1285
values[5] = vals[1];
1286
}
1287
1288
/// Interpolate vertex values from dof values
1289
virtual void interpolate_vertex_values(double* vertex_values,
1290
const double* dof_values,
1291
const ufc::cell& c) const
1292
{
1293
// Evaluate function and change variables
1294
vertex_values[0] = dof_values[0];
1295
vertex_values[2] = dof_values[1];
1296
vertex_values[4] = dof_values[2];
1297
// Evaluate function and change variables
1298
vertex_values[1] = dof_values[3];
1299
vertex_values[3] = dof_values[4];
1300
vertex_values[5] = dof_values[5];
1301
}
1302
1303
/// Return the number of sub elements (for a mixed element)
1304
virtual unsigned int num_sub_elements() const
1305
{
1306
return 2;
1307
}
1308
1309
/// Create a new finite element for sub element i (for a mixed element)
1310
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1311
{
1312
switch (i)
1313
{
1314
case 0:
1315
{
1316
return new stabilisedstokes_finite_element_0();
1317
break;
1318
}
1319
case 1:
1320
{
1321
return new stabilisedstokes_finite_element_0();
1322
break;
1323
}
1324
}
1325
1326
return 0;
1327
}
1328
1329
};
1330
1331
/// This class defines the interface for a finite element.
1332
1333
class stabilisedstokes_finite_element_2: public ufc::finite_element
1334
{
1335
public:
1336
1337
/// Constructor
1338
stabilisedstokes_finite_element_2() : ufc::finite_element()
1339
{
1340
// Do nothing
1341
}
1342
1343
/// Destructor
1344
virtual ~stabilisedstokes_finite_element_2()
1345
{
1346
// Do nothing
1347
}
1348
1349
/// Return a string identifying the finite element
1350
virtual const char* signature() const
1351
{
1352
return "MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) })";
1353
}
1354
1355
/// Return the cell shape
1356
virtual ufc::shape cell_shape() const
1357
{
1358
return ufc::triangle;
1359
}
1360
1361
/// Return the dimension of the finite element function space
1362
virtual unsigned int space_dimension() const
1363
{
1364
return 9;
1365
}
1366
1367
/// Return the rank of the value space
1368
virtual unsigned int value_rank() const
1369
{
1370
return 1;
1371
}
1372
1373
/// Return the dimension of the value space for axis i
1374
virtual unsigned int value_dimension(unsigned int i) const
1375
{
1376
switch (i)
1377
{
1378
case 0:
1379
{
1380
return 3;
1381
break;
1382
}
1383
}
1384
1385
return 0;
1386
}
1387
1388
/// Evaluate basis function i at given point in cell
1389
virtual void evaluate_basis(unsigned int i,
1390
double* values,
1391
const double* coordinates,
1392
const ufc::cell& c) const
1393
{
1394
// Extract vertex coordinates
1395
const double * const * x = c.coordinates;
1396
1397
// Compute Jacobian of affine map from reference cell
1398
const double J_00 = x[1][0] - x[0][0];
1399
const double J_01 = x[2][0] - x[0][0];
1400
const double J_10 = x[1][1] - x[0][1];
1401
const double J_11 = x[2][1] - x[0][1];
1402
1403
// Compute determinant of Jacobian
1404
double detJ = J_00*J_11 - J_01*J_10;
1405
1406
// Compute inverse of Jacobian
1407
1408
// Compute constants
1409
const double C0 = x[1][0] + x[2][0];
1410
const double C1 = x[1][1] + x[2][1];
1411
1412
// Get coordinates and map to the reference (FIAT) element
1413
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1414
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
1415
1416
// Reset values
1417
values[0] = 0.00000000;
1418
values[1] = 0.00000000;
1419
values[2] = 0.00000000;
1420
if (0 <= i && i <= 2)
1421
{
1422
// Map degree of freedom to element degree of freedom
1423
const unsigned int dof = i;
1424
1425
// Array of basisvalues
1426
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1427
1428
// Declare helper variables
1429
unsigned int rr = 0;
1430
unsigned int ss = 0;
1431
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1432
1433
// Compute basisvalues
1434
basisvalues[0] = 1.00000000;
1435
basisvalues[1] = tmp0;
1436
for (unsigned int r = 0; r < 1; r++)
1437
{
1438
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1439
ss = r*(r + 1)/2;
1440
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1441
}// end loop over 'r'
1442
for (unsigned int r = 0; r < 2; r++)
1443
{
1444
for (unsigned int s = 0; s < 2 - r; s++)
1445
{
1446
rr = (r + s)*(r + s + 1)/2 + s;
1447
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1448
}// end loop over 's'
1449
}// end loop over 'r'
1450
1451
// Table(s) of coefficients
1452
static const double coefficients0[3][3] = \
1453
{{0.47140452, -0.28867513, -0.16666667},
1454
{0.47140452, 0.28867513, -0.16666667},
1455
{0.47140452, 0.00000000, 0.33333333}};
1456
1457
// Compute value(s).
1458
for (unsigned int r = 0; r < 3; r++)
1459
{
1460
values[0] += coefficients0[dof][r]*basisvalues[r];
1461
}// end loop over 'r'
1462
}
1463
1464
if (3 <= i && i <= 5)
1465
{
1466
// Map degree of freedom to element degree of freedom
1467
const unsigned int dof = i - 3;
1468
1469
// Array of basisvalues
1470
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1471
1472
// Declare helper variables
1473
unsigned int rr = 0;
1474
unsigned int ss = 0;
1475
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1476
1477
// Compute basisvalues
1478
basisvalues[0] = 1.00000000;
1479
basisvalues[1] = tmp0;
1480
for (unsigned int r = 0; r < 1; r++)
1481
{
1482
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1483
ss = r*(r + 1)/2;
1484
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1485
}// end loop over 'r'
1486
for (unsigned int r = 0; r < 2; r++)
1487
{
1488
for (unsigned int s = 0; s < 2 - r; s++)
1489
{
1490
rr = (r + s)*(r + s + 1)/2 + s;
1491
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1492
}// end loop over 's'
1493
}// end loop over 'r'
1494
1495
// Table(s) of coefficients
1496
static const double coefficients0[3][3] = \
1497
{{0.47140452, -0.28867513, -0.16666667},
1498
{0.47140452, 0.28867513, -0.16666667},
1499
{0.47140452, 0.00000000, 0.33333333}};
1500
1501
// Compute value(s).
1502
for (unsigned int r = 0; r < 3; r++)
1503
{
1504
values[1] += coefficients0[dof][r]*basisvalues[r];
1505
}// end loop over 'r'
1506
}
1507
1508
if (6 <= i && i <= 8)
1509
{
1510
// Map degree of freedom to element degree of freedom
1511
const unsigned int dof = i - 6;
1512
1513
// Array of basisvalues
1514
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1515
1516
// Declare helper variables
1517
unsigned int rr = 0;
1518
unsigned int ss = 0;
1519
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1520
1521
// Compute basisvalues
1522
basisvalues[0] = 1.00000000;
1523
basisvalues[1] = tmp0;
1524
for (unsigned int r = 0; r < 1; r++)
1525
{
1526
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1527
ss = r*(r + 1)/2;
1528
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1529
}// end loop over 'r'
1530
for (unsigned int r = 0; r < 2; r++)
1531
{
1532
for (unsigned int s = 0; s < 2 - r; s++)
1533
{
1534
rr = (r + s)*(r + s + 1)/2 + s;
1535
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1536
}// end loop over 's'
1537
}// end loop over 'r'
1538
1539
// Table(s) of coefficients
1540
static const double coefficients0[3][3] = \
1541
{{0.47140452, -0.28867513, -0.16666667},
1542
{0.47140452, 0.28867513, -0.16666667},
1543
{0.47140452, 0.00000000, 0.33333333}};
1544
1545
// Compute value(s).
1546
for (unsigned int r = 0; r < 3; r++)
1547
{
1548
values[2] += coefficients0[dof][r]*basisvalues[r];
1549
}// end loop over 'r'
1550
}
1551
1552
}
1553
1554
/// Evaluate all basis functions at given point in cell
1555
virtual void evaluate_basis_all(double* values,
1556
const double* coordinates,
1557
const ufc::cell& c) const
1558
{
1559
// Helper variable to hold values of a single dof.
1560
double dof_values[3] = {0.00000000, 0.00000000, 0.00000000};
1561
1562
// Loop dofs and call evaluate_basis.
1563
for (unsigned int r = 0; r < 9; r++)
1564
{
1565
evaluate_basis(r, dof_values, coordinates, c);
1566
for (unsigned int s = 0; s < 3; s++)
1567
{
1568
values[r*3 + s] = dof_values[s];
1569
}// end loop over 's'
1570
}// end loop over 'r'
1571
}
1572
1573
/// Evaluate order n derivatives of basis function i at given point in cell
1574
virtual void evaluate_basis_derivatives(unsigned int i,
1575
unsigned int n,
1576
double* values,
1577
const double* coordinates,
1578
const ufc::cell& c) const
1579
{
1580
// Extract vertex coordinates
1581
const double * const * x = c.coordinates;
1582
1583
// Compute Jacobian of affine map from reference cell
1584
const double J_00 = x[1][0] - x[0][0];
1585
const double J_01 = x[2][0] - x[0][0];
1586
const double J_10 = x[1][1] - x[0][1];
1587
const double J_11 = x[2][1] - x[0][1];
1588
1589
// Compute determinant of Jacobian
1590
double detJ = J_00*J_11 - J_01*J_10;
1591
1592
// Compute inverse of Jacobian
1593
const double K_00 = J_11 / detJ;
1594
const double K_01 = -J_01 / detJ;
1595
const double K_10 = -J_10 / detJ;
1596
const double K_11 = J_00 / detJ;
1597
1598
// Compute constants
1599
const double C0 = x[1][0] + x[2][0];
1600
const double C1 = x[1][1] + x[2][1];
1601
1602
// Get coordinates and map to the reference (FIAT) element
1603
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1604
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
1605
1606
// Compute number of derivatives.
1607
unsigned int num_derivatives = 1;
1608
1609
for (unsigned int r = 0; r < n; r++)
1610
{
1611
num_derivatives *= 2;
1612
}// end loop over 'r'
1613
1614
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1615
unsigned int **combinations = new unsigned int *[num_derivatives];
1616
for (unsigned int row = 0; row < num_derivatives; row++)
1617
{
1618
combinations[row] = new unsigned int [n];
1619
for (unsigned int col = 0; col < n; col++)
1620
combinations[row][col] = 0;
1621
}
1622
1623
// Generate combinations of derivatives
1624
for (unsigned int row = 1; row < num_derivatives; row++)
1625
{
1626
for (unsigned int num = 0; num < row; num++)
1627
{
1628
for (unsigned int col = n-1; col+1 > 0; col--)
1629
{
1630
if (combinations[row][col] + 1 > 1)
1631
combinations[row][col] = 0;
1632
else
1633
{
1634
combinations[row][col] += 1;
1635
break;
1636
}
1637
}
1638
}
1639
}
1640
1641
// Compute inverse of Jacobian
1642
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
1643
1644
// Declare transformation matrix
1645
// Declare pointer to two dimensional array and initialise
1646
double **transform = new double *[num_derivatives];
1647
1648
for (unsigned int j = 0; j < num_derivatives; j++)
1649
{
1650
transform[j] = new double [num_derivatives];
1651
for (unsigned int k = 0; k < num_derivatives; k++)
1652
transform[j][k] = 1;
1653
}
1654
1655
// Construct transformation matrix
1656
for (unsigned int row = 0; row < num_derivatives; row++)
1657
{
1658
for (unsigned int col = 0; col < num_derivatives; col++)
1659
{
1660
for (unsigned int k = 0; k < n; k++)
1661
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1662
}
1663
}
1664
1665
// Reset values. Assuming that values is always an array.
1666
for (unsigned int r = 0; r < 3*num_derivatives; r++)
1667
{
1668
values[r] = 0.00000000;
1669
}// end loop over 'r'
1670
1671
if (0 <= i && i <= 2)
1672
{
1673
// Map degree of freedom to element degree of freedom
1674
const unsigned int dof = i;
1675
1676
// Array of basisvalues
1677
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1678
1679
// Declare helper variables
1680
unsigned int rr = 0;
1681
unsigned int ss = 0;
1682
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1683
1684
// Compute basisvalues
1685
basisvalues[0] = 1.00000000;
1686
basisvalues[1] = tmp0;
1687
for (unsigned int r = 0; r < 1; r++)
1688
{
1689
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1690
ss = r*(r + 1)/2;
1691
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1692
}// end loop over 'r'
1693
for (unsigned int r = 0; r < 2; r++)
1694
{
1695
for (unsigned int s = 0; s < 2 - r; s++)
1696
{
1697
rr = (r + s)*(r + s + 1)/2 + s;
1698
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1699
}// end loop over 's'
1700
}// end loop over 'r'
1701
1702
// Table(s) of coefficients
1703
static const double coefficients0[3][3] = \
1704
{{0.47140452, -0.28867513, -0.16666667},
1705
{0.47140452, 0.28867513, -0.16666667},
1706
{0.47140452, 0.00000000, 0.33333333}};
1707
1708
// Tables of derivatives of the polynomial base (transpose).
1709
static const double dmats0[3][3] = \
1710
{{0.00000000, 0.00000000, 0.00000000},
1711
{4.89897949, 0.00000000, 0.00000000},
1712
{0.00000000, 0.00000000, 0.00000000}};
1713
1714
static const double dmats1[3][3] = \
1715
{{0.00000000, 0.00000000, 0.00000000},
1716
{2.44948974, 0.00000000, 0.00000000},
1717
{4.24264069, 0.00000000, 0.00000000}};
1718
1719
// Compute reference derivatives
1720
// Declare pointer to array of derivatives on FIAT element
1721
double *derivatives = new double [num_derivatives];
1722
for (unsigned int r = 0; r < num_derivatives; r++)
1723
{
1724
derivatives[r] = 0.00000000;
1725
}// end loop over 'r'
1726
1727
// Declare derivative matrix (of polynomial basis).
1728
double dmats[3][3] = \
1729
{{1.00000000, 0.00000000, 0.00000000},
1730
{0.00000000, 1.00000000, 0.00000000},
1731
{0.00000000, 0.00000000, 1.00000000}};
1732
1733
// Declare (auxiliary) derivative matrix (of polynomial basis).
1734
double dmats_old[3][3] = \
1735
{{1.00000000, 0.00000000, 0.00000000},
1736
{0.00000000, 1.00000000, 0.00000000},
1737
{0.00000000, 0.00000000, 1.00000000}};
1738
1739
// Loop possible derivatives.
1740
for (unsigned int r = 0; r < num_derivatives; r++)
1741
{
1742
// Resetting dmats values to compute next derivative.
1743
for (unsigned int t = 0; t < 3; t++)
1744
{
1745
for (unsigned int u = 0; u < 3; u++)
1746
{
1747
dmats[t][u] = 0.00000000;
1748
if (t == u)
1749
{
1750
dmats[t][u] = 1.00000000;
1751
}
1752
1753
}// end loop over 'u'
1754
}// end loop over 't'
1755
1756
// Looping derivative order to generate dmats.
1757
for (unsigned int s = 0; s < n; s++)
1758
{
1759
// Updating dmats_old with new values and resetting dmats.
1760
for (unsigned int t = 0; t < 3; t++)
1761
{
1762
for (unsigned int u = 0; u < 3; u++)
1763
{
1764
dmats_old[t][u] = dmats[t][u];
1765
dmats[t][u] = 0.00000000;
1766
}// end loop over 'u'
1767
}// end loop over 't'
1768
1769
// Update dmats using an inner product.
1770
if (combinations[r][s] == 0)
1771
{
1772
for (unsigned int t = 0; t < 3; t++)
1773
{
1774
for (unsigned int u = 0; u < 3; u++)
1775
{
1776
for (unsigned int tu = 0; tu < 3; tu++)
1777
{
1778
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1779
}// end loop over 'tu'
1780
}// end loop over 'u'
1781
}// end loop over 't'
1782
}
1783
1784
if (combinations[r][s] == 1)
1785
{
1786
for (unsigned int t = 0; t < 3; t++)
1787
{
1788
for (unsigned int u = 0; u < 3; u++)
1789
{
1790
for (unsigned int tu = 0; tu < 3; tu++)
1791
{
1792
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1793
}// end loop over 'tu'
1794
}// end loop over 'u'
1795
}// end loop over 't'
1796
}
1797
1798
}// end loop over 's'
1799
for (unsigned int s = 0; s < 3; s++)
1800
{
1801
for (unsigned int t = 0; t < 3; t++)
1802
{
1803
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1804
}// end loop over 't'
1805
}// end loop over 's'
1806
}// end loop over 'r'
1807
1808
// Transform derivatives back to physical element
1809
for (unsigned int row = 0; row < num_derivatives; row++)
1810
{
1811
for (unsigned int col = 0; col < num_derivatives; col++)
1812
{
1813
values[row] += transform[row][col]*derivatives[col];
1814
}
1815
}
1816
1817
// Delete pointer to array of derivatives on FIAT element
1818
delete [] derivatives;
1819
1820
// Delete pointer to array of combinations of derivatives and transform
1821
for (unsigned int r = 0; r < num_derivatives; r++)
1822
{
1823
delete [] combinations[r];
1824
delete [] transform[r];
1825
}// end loop over 'r'
1826
delete [] combinations;
1827
delete [] transform;
1828
}
1829
1830
if (3 <= i && i <= 5)
1831
{
1832
// Map degree of freedom to element degree of freedom
1833
const unsigned int dof = i - 3;
1834
1835
// Array of basisvalues
1836
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1837
1838
// Declare helper variables
1839
unsigned int rr = 0;
1840
unsigned int ss = 0;
1841
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1842
1843
// Compute basisvalues
1844
basisvalues[0] = 1.00000000;
1845
basisvalues[1] = tmp0;
1846
for (unsigned int r = 0; r < 1; r++)
1847
{
1848
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1849
ss = r*(r + 1)/2;
1850
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1851
}// end loop over 'r'
1852
for (unsigned int r = 0; r < 2; r++)
1853
{
1854
for (unsigned int s = 0; s < 2 - r; s++)
1855
{
1856
rr = (r + s)*(r + s + 1)/2 + s;
1857
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1858
}// end loop over 's'
1859
}// end loop over 'r'
1860
1861
// Table(s) of coefficients
1862
static const double coefficients0[3][3] = \
1863
{{0.47140452, -0.28867513, -0.16666667},
1864
{0.47140452, 0.28867513, -0.16666667},
1865
{0.47140452, 0.00000000, 0.33333333}};
1866
1867
// Tables of derivatives of the polynomial base (transpose).
1868
static const double dmats0[3][3] = \
1869
{{0.00000000, 0.00000000, 0.00000000},
1870
{4.89897949, 0.00000000, 0.00000000},
1871
{0.00000000, 0.00000000, 0.00000000}};
1872
1873
static const double dmats1[3][3] = \
1874
{{0.00000000, 0.00000000, 0.00000000},
1875
{2.44948974, 0.00000000, 0.00000000},
1876
{4.24264069, 0.00000000, 0.00000000}};
1877
1878
// Compute reference derivatives
1879
// Declare pointer to array of derivatives on FIAT element
1880
double *derivatives = new double [num_derivatives];
1881
for (unsigned int r = 0; r < num_derivatives; r++)
1882
{
1883
derivatives[r] = 0.00000000;
1884
}// end loop over 'r'
1885
1886
// Declare derivative matrix (of polynomial basis).
1887
double dmats[3][3] = \
1888
{{1.00000000, 0.00000000, 0.00000000},
1889
{0.00000000, 1.00000000, 0.00000000},
1890
{0.00000000, 0.00000000, 1.00000000}};
1891
1892
// Declare (auxiliary) derivative matrix (of polynomial basis).
1893
double dmats_old[3][3] = \
1894
{{1.00000000, 0.00000000, 0.00000000},
1895
{0.00000000, 1.00000000, 0.00000000},
1896
{0.00000000, 0.00000000, 1.00000000}};
1897
1898
// Loop possible derivatives.
1899
for (unsigned int r = 0; r < num_derivatives; r++)
1900
{
1901
// Resetting dmats values to compute next derivative.
1902
for (unsigned int t = 0; t < 3; t++)
1903
{
1904
for (unsigned int u = 0; u < 3; u++)
1905
{
1906
dmats[t][u] = 0.00000000;
1907
if (t == u)
1908
{
1909
dmats[t][u] = 1.00000000;
1910
}
1911
1912
}// end loop over 'u'
1913
}// end loop over 't'
1914
1915
// Looping derivative order to generate dmats.
1916
for (unsigned int s = 0; s < n; s++)
1917
{
1918
// Updating dmats_old with new values and resetting dmats.
1919
for (unsigned int t = 0; t < 3; t++)
1920
{
1921
for (unsigned int u = 0; u < 3; u++)
1922
{
1923
dmats_old[t][u] = dmats[t][u];
1924
dmats[t][u] = 0.00000000;
1925
}// end loop over 'u'
1926
}// end loop over 't'
1927
1928
// Update dmats using an inner product.
1929
if (combinations[r][s] == 0)
1930
{
1931
for (unsigned int t = 0; t < 3; t++)
1932
{
1933
for (unsigned int u = 0; u < 3; u++)
1934
{
1935
for (unsigned int tu = 0; tu < 3; tu++)
1936
{
1937
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1938
}// end loop over 'tu'
1939
}// end loop over 'u'
1940
}// end loop over 't'
1941
}
1942
1943
if (combinations[r][s] == 1)
1944
{
1945
for (unsigned int t = 0; t < 3; t++)
1946
{
1947
for (unsigned int u = 0; u < 3; u++)
1948
{
1949
for (unsigned int tu = 0; tu < 3; tu++)
1950
{
1951
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1952
}// end loop over 'tu'
1953
}// end loop over 'u'
1954
}// end loop over 't'
1955
}
1956
1957
}// end loop over 's'
1958
for (unsigned int s = 0; s < 3; s++)
1959
{
1960
for (unsigned int t = 0; t < 3; t++)
1961
{
1962
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1963
}// end loop over 't'
1964
}// end loop over 's'
1965
}// end loop over 'r'
1966
1967
// Transform derivatives back to physical element
1968
for (unsigned int row = 0; row < num_derivatives; row++)
1969
{
1970
for (unsigned int col = 0; col < num_derivatives; col++)
1971
{
1972
values[num_derivatives + row] += transform[row][col]*derivatives[col];
1973
}
1974
}
1975
1976
// Delete pointer to array of derivatives on FIAT element
1977
delete [] derivatives;
1978
1979
// Delete pointer to array of combinations of derivatives and transform
1980
for (unsigned int r = 0; r < num_derivatives; r++)
1981
{
1982
delete [] combinations[r];
1983
delete [] transform[r];
1984
}// end loop over 'r'
1985
delete [] combinations;
1986
delete [] transform;
1987
}
1988
1989
if (6 <= i && i <= 8)
1990
{
1991
// Map degree of freedom to element degree of freedom
1992
const unsigned int dof = i - 6;
1993
1994
// Array of basisvalues
1995
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1996
1997
// Declare helper variables
1998
unsigned int rr = 0;
1999
unsigned int ss = 0;
2000
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
2001
2002
// Compute basisvalues
2003
basisvalues[0] = 1.00000000;
2004
basisvalues[1] = tmp0;
2005
for (unsigned int r = 0; r < 1; r++)
2006
{
2007
rr = (r + 1)*(r + 1 + 1)/2 + 1;
2008
ss = r*(r + 1)/2;
2009
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
2010
}// end loop over 'r'
2011
for (unsigned int r = 0; r < 2; r++)
2012
{
2013
for (unsigned int s = 0; s < 2 - r; s++)
2014
{
2015
rr = (r + s)*(r + s + 1)/2 + s;
2016
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
2017
}// end loop over 's'
2018
}// end loop over 'r'
2019
2020
// Table(s) of coefficients
2021
static const double coefficients0[3][3] = \
2022
{{0.47140452, -0.28867513, -0.16666667},
2023
{0.47140452, 0.28867513, -0.16666667},
2024
{0.47140452, 0.00000000, 0.33333333}};
2025
2026
// Tables of derivatives of the polynomial base (transpose).
2027
static const double dmats0[3][3] = \
2028
{{0.00000000, 0.00000000, 0.00000000},
2029
{4.89897949, 0.00000000, 0.00000000},
2030
{0.00000000, 0.00000000, 0.00000000}};
2031
2032
static const double dmats1[3][3] = \
2033
{{0.00000000, 0.00000000, 0.00000000},
2034
{2.44948974, 0.00000000, 0.00000000},
2035
{4.24264069, 0.00000000, 0.00000000}};
2036
2037
// Compute reference derivatives
2038
// Declare pointer to array of derivatives on FIAT element
2039
double *derivatives = new double [num_derivatives];
2040
for (unsigned int r = 0; r < num_derivatives; r++)
2041
{
2042
derivatives[r] = 0.00000000;
2043
}// end loop over 'r'
2044
2045
// Declare derivative matrix (of polynomial basis).
2046
double dmats[3][3] = \
2047
{{1.00000000, 0.00000000, 0.00000000},
2048
{0.00000000, 1.00000000, 0.00000000},
2049
{0.00000000, 0.00000000, 1.00000000}};
2050
2051
// Declare (auxiliary) derivative matrix (of polynomial basis).
2052
double dmats_old[3][3] = \
2053
{{1.00000000, 0.00000000, 0.00000000},
2054
{0.00000000, 1.00000000, 0.00000000},
2055
{0.00000000, 0.00000000, 1.00000000}};
2056
2057
// Loop possible derivatives.
2058
for (unsigned int r = 0; r < num_derivatives; r++)
2059
{
2060
// Resetting dmats values to compute next derivative.
2061
for (unsigned int t = 0; t < 3; t++)
2062
{
2063
for (unsigned int u = 0; u < 3; u++)
2064
{
2065
dmats[t][u] = 0.00000000;
2066
if (t == u)
2067
{
2068
dmats[t][u] = 1.00000000;
2069
}
2070
2071
}// end loop over 'u'
2072
}// end loop over 't'
2073
2074
// Looping derivative order to generate dmats.
2075
for (unsigned int s = 0; s < n; s++)
2076
{
2077
// Updating dmats_old with new values and resetting dmats.
2078
for (unsigned int t = 0; t < 3; t++)
2079
{
2080
for (unsigned int u = 0; u < 3; u++)
2081
{
2082
dmats_old[t][u] = dmats[t][u];
2083
dmats[t][u] = 0.00000000;
2084
}// end loop over 'u'
2085
}// end loop over 't'
2086
2087
// Update dmats using an inner product.
2088
if (combinations[r][s] == 0)
2089
{
2090
for (unsigned int t = 0; t < 3; t++)
2091
{
2092
for (unsigned int u = 0; u < 3; u++)
2093
{
2094
for (unsigned int tu = 0; tu < 3; tu++)
2095
{
2096
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2097
}// end loop over 'tu'
2098
}// end loop over 'u'
2099
}// end loop over 't'
2100
}
2101
2102
if (combinations[r][s] == 1)
2103
{
2104
for (unsigned int t = 0; t < 3; t++)
2105
{
2106
for (unsigned int u = 0; u < 3; u++)
2107
{
2108
for (unsigned int tu = 0; tu < 3; tu++)
2109
{
2110
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2111
}// end loop over 'tu'
2112
}// end loop over 'u'
2113
}// end loop over 't'
2114
}
2115
2116
}// end loop over 's'
2117
for (unsigned int s = 0; s < 3; s++)
2118
{
2119
for (unsigned int t = 0; t < 3; t++)
2120
{
2121
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
2122
}// end loop over 't'
2123
}// end loop over 's'
2124
}// end loop over 'r'
2125
2126
// Transform derivatives back to physical element
2127
for (unsigned int row = 0; row < num_derivatives; row++)
2128
{
2129
for (unsigned int col = 0; col < num_derivatives; col++)
2130
{
2131
values[2*num_derivatives + row] += transform[row][col]*derivatives[col];
2132
}
2133
}
2134
2135
// Delete pointer to array of derivatives on FIAT element
2136
delete [] derivatives;
2137
2138
// Delete pointer to array of combinations of derivatives and transform
2139
for (unsigned int r = 0; r < num_derivatives; r++)
2140
{
2141
delete [] combinations[r];
2142
delete [] transform[r];
2143
}// end loop over 'r'
2144
delete [] combinations;
2145
delete [] transform;
2146
}
2147
2148
}
2149
2150
/// Evaluate order n derivatives of all basis functions at given point in cell
2151
virtual void evaluate_basis_derivatives_all(unsigned int n,
2152
double* values,
2153
const double* coordinates,
2154
const ufc::cell& c) const
2155
{
2156
// Compute number of derivatives.
2157
unsigned int num_derivatives = 1;
2158
2159
for (unsigned int r = 0; r < n; r++)
2160
{
2161
num_derivatives *= 2;
2162
}// end loop over 'r'
2163
2164
// Helper variable to hold values of a single dof.
2165
double *dof_values = new double [3*num_derivatives];
2166
for (unsigned int r = 0; r < 3*num_derivatives; r++)
2167
{
2168
dof_values[r] = 0.00000000;
2169
}// end loop over 'r'
2170
2171
// Loop dofs and call evaluate_basis_derivatives.
2172
for (unsigned int r = 0; r < 9; r++)
2173
{
2174
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
2175
for (unsigned int s = 0; s < 3*num_derivatives; s++)
2176
{
2177
values[r*3*num_derivatives + s] = dof_values[s];
2178
}// end loop over 's'
2179
}// end loop over 'r'
2180
2181
// Delete pointer.
2182
delete [] dof_values;
2183
}
2184
2185
/// Evaluate linear functional for dof i on the function f
2186
virtual double evaluate_dof(unsigned int i,
2187
const ufc::function& f,
2188
const ufc::cell& c) const
2189
{
2190
// Declare variables for result of evaluation
2191
double vals[3];
2192
2193
// Declare variable for physical coordinates
2194
double y[2];
2195
2196
const double * const * x = c.coordinates;
2197
switch (i)
2198
{
2199
case 0:
2200
{
2201
y[0] = x[0][0];
2202
y[1] = x[0][1];
2203
f.evaluate(vals, y, c);
2204
return vals[0];
2205
break;
2206
}
2207
case 1:
2208
{
2209
y[0] = x[1][0];
2210
y[1] = x[1][1];
2211
f.evaluate(vals, y, c);
2212
return vals[0];
2213
break;
2214
}
2215
case 2:
2216
{
2217
y[0] = x[2][0];
2218
y[1] = x[2][1];
2219
f.evaluate(vals, y, c);
2220
return vals[0];
2221
break;
2222
}
2223
case 3:
2224
{
2225
y[0] = x[0][0];
2226
y[1] = x[0][1];
2227
f.evaluate(vals, y, c);
2228
return vals[1];
2229
break;
2230
}
2231
case 4:
2232
{
2233
y[0] = x[1][0];
2234
y[1] = x[1][1];
2235
f.evaluate(vals, y, c);
2236
return vals[1];
2237
break;
2238
}
2239
case 5:
2240
{
2241
y[0] = x[2][0];
2242
y[1] = x[2][1];
2243
f.evaluate(vals, y, c);
2244
return vals[1];
2245
break;
2246
}
2247
case 6:
2248
{
2249
y[0] = x[0][0];
2250
y[1] = x[0][1];
2251
f.evaluate(vals, y, c);
2252
return vals[2];
2253
break;
2254
}
2255
case 7:
2256
{
2257
y[0] = x[1][0];
2258
y[1] = x[1][1];
2259
f.evaluate(vals, y, c);
2260
return vals[2];
2261
break;
2262
}
2263
case 8:
2264
{
2265
y[0] = x[2][0];
2266
y[1] = x[2][1];
2267
f.evaluate(vals, y, c);
2268
return vals[2];
2269
break;
2270
}
2271
}
2272
2273
return 0.0;
2274
}
2275
2276
/// Evaluate linear functionals for all dofs on the function f
2277
virtual void evaluate_dofs(double* values,
2278
const ufc::function& f,
2279
const ufc::cell& c) const
2280
{
2281
// Declare variables for result of evaluation
2282
double vals[3];
2283
2284
// Declare variable for physical coordinates
2285
double y[2];
2286
2287
const double * const * x = c.coordinates;
2288
y[0] = x[0][0];
2289
y[1] = x[0][1];
2290
f.evaluate(vals, y, c);
2291
values[0] = vals[0];
2292
y[0] = x[1][0];
2293
y[1] = x[1][1];
2294
f.evaluate(vals, y, c);
2295
values[1] = vals[0];
2296
y[0] = x[2][0];
2297
y[1] = x[2][1];
2298
f.evaluate(vals, y, c);
2299
values[2] = vals[0];
2300
y[0] = x[0][0];
2301
y[1] = x[0][1];
2302
f.evaluate(vals, y, c);
2303
values[3] = vals[1];
2304
y[0] = x[1][0];
2305
y[1] = x[1][1];
2306
f.evaluate(vals, y, c);
2307
values[4] = vals[1];
2308
y[0] = x[2][0];
2309
y[1] = x[2][1];
2310
f.evaluate(vals, y, c);
2311
values[5] = vals[1];
2312
y[0] = x[0][0];
2313
y[1] = x[0][1];
2314
f.evaluate(vals, y, c);
2315
values[6] = vals[2];
2316
y[0] = x[1][0];
2317
y[1] = x[1][1];
2318
f.evaluate(vals, y, c);
2319
values[7] = vals[2];
2320
y[0] = x[2][0];
2321
y[1] = x[2][1];
2322
f.evaluate(vals, y, c);
2323
values[8] = vals[2];
2324
}
2325
2326
/// Interpolate vertex values from dof values
2327
virtual void interpolate_vertex_values(double* vertex_values,
2328
const double* dof_values,
2329
const ufc::cell& c) const
2330
{
2331
// Evaluate function and change variables
2332
vertex_values[0] = dof_values[0];
2333
vertex_values[3] = dof_values[1];
2334
vertex_values[6] = dof_values[2];
2335
// Evaluate function and change variables
2336
vertex_values[1] = dof_values[3];
2337
vertex_values[4] = dof_values[4];
2338
vertex_values[7] = dof_values[5];
2339
// Evaluate function and change variables
2340
vertex_values[2] = dof_values[6];
2341
vertex_values[5] = dof_values[7];
2342
vertex_values[8] = dof_values[8];
2343
}
2344
2345
/// Return the number of sub elements (for a mixed element)
2346
virtual unsigned int num_sub_elements() const
2347
{
2348
return 2;
2349
}
2350
2351
/// Create a new finite element for sub element i (for a mixed element)
2352
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2353
{
2354
switch (i)
2355
{
2356
case 0:
2357
{
2358
return new stabilisedstokes_finite_element_1();
2359
break;
2360
}
2361
case 1:
2362
{
2363
return new stabilisedstokes_finite_element_0();
2364
break;
2365
}
2366
}
2367
2368
return 0;
2369
}
2370
2371
};
2372
2373
/// This class defines the interface for a local-to-global mapping of
2374
/// degrees of freedom (dofs).
2375
2376
class stabilisedstokes_dof_map_0: public ufc::dof_map
2377
{
2378
private:
2379
2380
unsigned int _global_dimension;
2381
2382
public:
2383
2384
/// Constructor
2385
stabilisedstokes_dof_map_0() : ufc::dof_map()
2386
{
2387
_global_dimension = 0;
2388
}
2389
2390
/// Destructor
2391
virtual ~stabilisedstokes_dof_map_0()
2392
{
2393
// Do nothing
2394
}
2395
2396
/// Return a string identifying the dof map
2397
virtual const char* signature() const
2398
{
2399
return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
2400
}
2401
2402
/// Return true iff mesh entities of topological dimension d are needed
2403
virtual bool needs_mesh_entities(unsigned int d) const
2404
{
2405
switch (d)
2406
{
2407
case 0:
2408
{
2409
return true;
2410
break;
2411
}
2412
case 1:
2413
{
2414
return false;
2415
break;
2416
}
2417
case 2:
2418
{
2419
return false;
2420
break;
2421
}
2422
}
2423
2424
return false;
2425
}
2426
2427
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2428
virtual bool init_mesh(const ufc::mesh& m)
2429
{
2430
_global_dimension = m.num_entities[0];
2431
return false;
2432
}
2433
2434
/// Initialize dof map for given cell
2435
virtual void init_cell(const ufc::mesh& m,
2436
const ufc::cell& c)
2437
{
2438
// Do nothing
2439
}
2440
2441
/// Finish initialization of dof map for cells
2442
virtual void init_cell_finalize()
2443
{
2444
// Do nothing
2445
}
2446
2447
/// Return the dimension of the global finite element function space
2448
virtual unsigned int global_dimension() const
2449
{
2450
return _global_dimension;
2451
}
2452
2453
/// Return the dimension of the local finite element function space for a cell
2454
virtual unsigned int local_dimension(const ufc::cell& c) const
2455
{
2456
return 3;
2457
}
2458
2459
/// Return the maximum dimension of the local finite element function space
2460
virtual unsigned int max_local_dimension() const
2461
{
2462
return 3;
2463
}
2464
2465
// Return the geometric dimension of the coordinates this dof map provides
2466
virtual unsigned int geometric_dimension() const
2467
{
2468
return 2;
2469
}
2470
2471
/// Return the number of dofs on each cell facet
2472
virtual unsigned int num_facet_dofs() const
2473
{
2474
return 2;
2475
}
2476
2477
/// Return the number of dofs associated with each cell entity of dimension d
2478
virtual unsigned int num_entity_dofs(unsigned int d) const
2479
{
2480
switch (d)
2481
{
2482
case 0:
2483
{
2484
return 1;
2485
break;
2486
}
2487
case 1:
2488
{
2489
return 0;
2490
break;
2491
}
2492
case 2:
2493
{
2494
return 0;
2495
break;
2496
}
2497
}
2498
2499
return 0;
2500
}
2501
2502
/// Tabulate the local-to-global mapping of dofs on a cell
2503
virtual void tabulate_dofs(unsigned int* dofs,
2504
const ufc::mesh& m,
2505
const ufc::cell& c) const
2506
{
2507
dofs[0] = c.entity_indices[0][0];
2508
dofs[1] = c.entity_indices[0][1];
2509
dofs[2] = c.entity_indices[0][2];
2510
}
2511
2512
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2513
virtual void tabulate_facet_dofs(unsigned int* dofs,
2514
unsigned int facet) const
2515
{
2516
switch (facet)
2517
{
2518
case 0:
2519
{
2520
dofs[0] = 1;
2521
dofs[1] = 2;
2522
break;
2523
}
2524
case 1:
2525
{
2526
dofs[0] = 0;
2527
dofs[1] = 2;
2528
break;
2529
}
2530
case 2:
2531
{
2532
dofs[0] = 0;
2533
dofs[1] = 1;
2534
break;
2535
}
2536
}
2537
2538
}
2539
2540
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2541
virtual void tabulate_entity_dofs(unsigned int* dofs,
2542
unsigned int d, unsigned int i) const
2543
{
2544
if (d > 2)
2545
{
2546
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2547
}
2548
2549
switch (d)
2550
{
2551
case 0:
2552
{
2553
if (i > 2)
2554
{
2555
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2556
}
2557
2558
switch (i)
2559
{
2560
case 0:
2561
{
2562
dofs[0] = 0;
2563
break;
2564
}
2565
case 1:
2566
{
2567
dofs[0] = 1;
2568
break;
2569
}
2570
case 2:
2571
{
2572
dofs[0] = 2;
2573
break;
2574
}
2575
}
2576
2577
break;
2578
}
2579
case 1:
2580
{
2581
2582
break;
2583
}
2584
case 2:
2585
{
2586
2587
break;
2588
}
2589
}
2590
2591
}
2592
2593
/// Tabulate the coordinates of all dofs on a cell
2594
virtual void tabulate_coordinates(double** coordinates,
2595
const ufc::cell& c) const
2596
{
2597
const double * const * x = c.coordinates;
2598
2599
coordinates[0][0] = x[0][0];
2600
coordinates[0][1] = x[0][1];
2601
coordinates[1][0] = x[1][0];
2602
coordinates[1][1] = x[1][1];
2603
coordinates[2][0] = x[2][0];
2604
coordinates[2][1] = x[2][1];
2605
}
2606
2607
/// Return the number of sub dof maps (for a mixed element)
2608
virtual unsigned int num_sub_dof_maps() const
2609
{
2610
return 0;
2611
}
2612
2613
/// Create a new dof_map for sub dof map i (for a mixed element)
2614
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2615
{
2616
return 0;
2617
}
2618
2619
};
2620
2621
/// This class defines the interface for a local-to-global mapping of
2622
/// degrees of freedom (dofs).
2623
2624
class stabilisedstokes_dof_map_1: public ufc::dof_map
2625
{
2626
private:
2627
2628
unsigned int _global_dimension;
2629
2630
public:
2631
2632
/// Constructor
2633
stabilisedstokes_dof_map_1() : ufc::dof_map()
2634
{
2635
_global_dimension = 0;
2636
}
2637
2638
/// Destructor
2639
virtual ~stabilisedstokes_dof_map_1()
2640
{
2641
// Do nothing
2642
}
2643
2644
/// Return a string identifying the dof map
2645
virtual const char* signature() const
2646
{
2647
return "FFC dofmap for VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2)";
2648
}
2649
2650
/// Return true iff mesh entities of topological dimension d are needed
2651
virtual bool needs_mesh_entities(unsigned int d) const
2652
{
2653
switch (d)
2654
{
2655
case 0:
2656
{
2657
return true;
2658
break;
2659
}
2660
case 1:
2661
{
2662
return false;
2663
break;
2664
}
2665
case 2:
2666
{
2667
return false;
2668
break;
2669
}
2670
}
2671
2672
return false;
2673
}
2674
2675
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2676
virtual bool init_mesh(const ufc::mesh& m)
2677
{
2678
_global_dimension = 2*m.num_entities[0];
2679
return false;
2680
}
2681
2682
/// Initialize dof map for given cell
2683
virtual void init_cell(const ufc::mesh& m,
2684
const ufc::cell& c)
2685
{
2686
// Do nothing
2687
}
2688
2689
/// Finish initialization of dof map for cells
2690
virtual void init_cell_finalize()
2691
{
2692
// Do nothing
2693
}
2694
2695
/// Return the dimension of the global finite element function space
2696
virtual unsigned int global_dimension() const
2697
{
2698
return _global_dimension;
2699
}
2700
2701
/// Return the dimension of the local finite element function space for a cell
2702
virtual unsigned int local_dimension(const ufc::cell& c) const
2703
{
2704
return 6;
2705
}
2706
2707
/// Return the maximum dimension of the local finite element function space
2708
virtual unsigned int max_local_dimension() const
2709
{
2710
return 6;
2711
}
2712
2713
// Return the geometric dimension of the coordinates this dof map provides
2714
virtual unsigned int geometric_dimension() const
2715
{
2716
return 2;
2717
}
2718
2719
/// Return the number of dofs on each cell facet
2720
virtual unsigned int num_facet_dofs() const
2721
{
2722
return 4;
2723
}
2724
2725
/// Return the number of dofs associated with each cell entity of dimension d
2726
virtual unsigned int num_entity_dofs(unsigned int d) const
2727
{
2728
switch (d)
2729
{
2730
case 0:
2731
{
2732
return 2;
2733
break;
2734
}
2735
case 1:
2736
{
2737
return 0;
2738
break;
2739
}
2740
case 2:
2741
{
2742
return 0;
2743
break;
2744
}
2745
}
2746
2747
return 0;
2748
}
2749
2750
/// Tabulate the local-to-global mapping of dofs on a cell
2751
virtual void tabulate_dofs(unsigned int* dofs,
2752
const ufc::mesh& m,
2753
const ufc::cell& c) const
2754
{
2755
unsigned int offset = 0;
2756
2757
dofs[0] = offset + c.entity_indices[0][0];
2758
dofs[1] = offset + c.entity_indices[0][1];
2759
dofs[2] = offset + c.entity_indices[0][2];
2760
offset += m.num_entities[0];
2761
dofs[3] = offset + c.entity_indices[0][0];
2762
dofs[4] = offset + c.entity_indices[0][1];
2763
dofs[5] = offset + c.entity_indices[0][2];
2764
offset += m.num_entities[0];
2765
}
2766
2767
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2768
virtual void tabulate_facet_dofs(unsigned int* dofs,
2769
unsigned int facet) const
2770
{
2771
switch (facet)
2772
{
2773
case 0:
2774
{
2775
dofs[0] = 1;
2776
dofs[1] = 2;
2777
dofs[2] = 4;
2778
dofs[3] = 5;
2779
break;
2780
}
2781
case 1:
2782
{
2783
dofs[0] = 0;
2784
dofs[1] = 2;
2785
dofs[2] = 3;
2786
dofs[3] = 5;
2787
break;
2788
}
2789
case 2:
2790
{
2791
dofs[0] = 0;
2792
dofs[1] = 1;
2793
dofs[2] = 3;
2794
dofs[3] = 4;
2795
break;
2796
}
2797
}
2798
2799
}
2800
2801
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2802
virtual void tabulate_entity_dofs(unsigned int* dofs,
2803
unsigned int d, unsigned int i) const
2804
{
2805
if (d > 2)
2806
{
2807
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2808
}
2809
2810
switch (d)
2811
{
2812
case 0:
2813
{
2814
if (i > 2)
2815
{
2816
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2817
}
2818
2819
switch (i)
2820
{
2821
case 0:
2822
{
2823
dofs[0] = 0;
2824
dofs[1] = 3;
2825
break;
2826
}
2827
case 1:
2828
{
2829
dofs[0] = 1;
2830
dofs[1] = 4;
2831
break;
2832
}
2833
case 2:
2834
{
2835
dofs[0] = 2;
2836
dofs[1] = 5;
2837
break;
2838
}
2839
}
2840
2841
break;
2842
}
2843
case 1:
2844
{
2845
2846
break;
2847
}
2848
case 2:
2849
{
2850
2851
break;
2852
}
2853
}
2854
2855
}
2856
2857
/// Tabulate the coordinates of all dofs on a cell
2858
virtual void tabulate_coordinates(double** coordinates,
2859
const ufc::cell& c) const
2860
{
2861
const double * const * x = c.coordinates;
2862
2863
coordinates[0][0] = x[0][0];
2864
coordinates[0][1] = x[0][1];
2865
coordinates[1][0] = x[1][0];
2866
coordinates[1][1] = x[1][1];
2867
coordinates[2][0] = x[2][0];
2868
coordinates[2][1] = x[2][1];
2869
coordinates[3][0] = x[0][0];
2870
coordinates[3][1] = x[0][1];
2871
coordinates[4][0] = x[1][0];
2872
coordinates[4][1] = x[1][1];
2873
coordinates[5][0] = x[2][0];
2874
coordinates[5][1] = x[2][1];
2875
}
2876
2877
/// Return the number of sub dof maps (for a mixed element)
2878
virtual unsigned int num_sub_dof_maps() const
2879
{
2880
return 2;
2881
}
2882
2883
/// Create a new dof_map for sub dof map i (for a mixed element)
2884
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2885
{
2886
switch (i)
2887
{
2888
case 0:
2889
{
2890
return new stabilisedstokes_dof_map_0();
2891
break;
2892
}
2893
case 1:
2894
{
2895
return new stabilisedstokes_dof_map_0();
2896
break;
2897
}
2898
}
2899
2900
return 0;
2901
}
2902
2903
};
2904
2905
/// This class defines the interface for a local-to-global mapping of
2906
/// degrees of freedom (dofs).
2907
2908
class stabilisedstokes_dof_map_2: public ufc::dof_map
2909
{
2910
private:
2911
2912
unsigned int _global_dimension;
2913
2914
public:
2915
2916
/// Constructor
2917
stabilisedstokes_dof_map_2() : ufc::dof_map()
2918
{
2919
_global_dimension = 0;
2920
}
2921
2922
/// Destructor
2923
virtual ~stabilisedstokes_dof_map_2()
2924
{
2925
// Do nothing
2926
}
2927
2928
/// Return a string identifying the dof map
2929
virtual const char* signature() const
2930
{
2931
return "FFC dofmap for MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) })";
2932
}
2933
2934
/// Return true iff mesh entities of topological dimension d are needed
2935
virtual bool needs_mesh_entities(unsigned int d) const
2936
{
2937
switch (d)
2938
{
2939
case 0:
2940
{
2941
return true;
2942
break;
2943
}
2944
case 1:
2945
{
2946
return false;
2947
break;
2948
}
2949
case 2:
2950
{
2951
return false;
2952
break;
2953
}
2954
}
2955
2956
return false;
2957
}
2958
2959
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2960
virtual bool init_mesh(const ufc::mesh& m)
2961
{
2962
_global_dimension = 3*m.num_entities[0];
2963
return false;
2964
}
2965
2966
/// Initialize dof map for given cell
2967
virtual void init_cell(const ufc::mesh& m,
2968
const ufc::cell& c)
2969
{
2970
// Do nothing
2971
}
2972
2973
/// Finish initialization of dof map for cells
2974
virtual void init_cell_finalize()
2975
{
2976
// Do nothing
2977
}
2978
2979
/// Return the dimension of the global finite element function space
2980
virtual unsigned int global_dimension() const
2981
{
2982
return _global_dimension;
2983
}
2984
2985
/// Return the dimension of the local finite element function space for a cell
2986
virtual unsigned int local_dimension(const ufc::cell& c) const
2987
{
2988
return 9;
2989
}
2990
2991
/// Return the maximum dimension of the local finite element function space
2992
virtual unsigned int max_local_dimension() const
2993
{
2994
return 9;
2995
}
2996
2997
// Return the geometric dimension of the coordinates this dof map provides
2998
virtual unsigned int geometric_dimension() const
2999
{
3000
return 2;
3001
}
3002
3003
/// Return the number of dofs on each cell facet
3004
virtual unsigned int num_facet_dofs() const
3005
{
3006
return 6;
3007
}
3008
3009
/// Return the number of dofs associated with each cell entity of dimension d
3010
virtual unsigned int num_entity_dofs(unsigned int d) const
3011
{
3012
switch (d)
3013
{
3014
case 0:
3015
{
3016
return 3;
3017
break;
3018
}
3019
case 1:
3020
{
3021
return 0;
3022
break;
3023
}
3024
case 2:
3025
{
3026
return 0;
3027
break;
3028
}
3029
}
3030
3031
return 0;
3032
}
3033
3034
/// Tabulate the local-to-global mapping of dofs on a cell
3035
virtual void tabulate_dofs(unsigned int* dofs,
3036
const ufc::mesh& m,
3037
const ufc::cell& c) const
3038
{
3039
unsigned int offset = 0;
3040
3041
dofs[0] = offset + c.entity_indices[0][0];
3042
dofs[1] = offset + c.entity_indices[0][1];
3043
dofs[2] = offset + c.entity_indices[0][2];
3044
offset += m.num_entities[0];
3045
dofs[3] = offset + c.entity_indices[0][0];
3046
dofs[4] = offset + c.entity_indices[0][1];
3047
dofs[5] = offset + c.entity_indices[0][2];
3048
offset += m.num_entities[0];
3049
dofs[6] = offset + c.entity_indices[0][0];
3050
dofs[7] = offset + c.entity_indices[0][1];
3051
dofs[8] = offset + c.entity_indices[0][2];
3052
offset += m.num_entities[0];
3053
}
3054
3055
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3056
virtual void tabulate_facet_dofs(unsigned int* dofs,
3057
unsigned int facet) const
3058
{
3059
switch (facet)
3060
{
3061
case 0:
3062
{
3063
dofs[0] = 1;
3064
dofs[1] = 2;
3065
dofs[2] = 4;
3066
dofs[3] = 5;
3067
dofs[4] = 7;
3068
dofs[5] = 8;
3069
break;
3070
}
3071
case 1:
3072
{
3073
dofs[0] = 0;
3074
dofs[1] = 2;
3075
dofs[2] = 3;
3076
dofs[3] = 5;
3077
dofs[4] = 6;
3078
dofs[5] = 8;
3079
break;
3080
}
3081
case 2:
3082
{
3083
dofs[0] = 0;
3084
dofs[1] = 1;
3085
dofs[2] = 3;
3086
dofs[3] = 4;
3087
dofs[4] = 6;
3088
dofs[5] = 7;
3089
break;
3090
}
3091
}
3092
3093
}
3094
3095
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3096
virtual void tabulate_entity_dofs(unsigned int* dofs,
3097
unsigned int d, unsigned int i) const
3098
{
3099
if (d > 2)
3100
{
3101
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3102
}
3103
3104
switch (d)
3105
{
3106
case 0:
3107
{
3108
if (i > 2)
3109
{
3110
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
3111
}
3112
3113
switch (i)
3114
{
3115
case 0:
3116
{
3117
dofs[0] = 0;
3118
dofs[1] = 3;
3119
dofs[2] = 6;
3120
break;
3121
}
3122
case 1:
3123
{
3124
dofs[0] = 1;
3125
dofs[1] = 4;
3126
dofs[2] = 7;
3127
break;
3128
}
3129
case 2:
3130
{
3131
dofs[0] = 2;
3132
dofs[1] = 5;
3133
dofs[2] = 8;
3134
break;
3135
}
3136
}
3137
3138
break;
3139
}
3140
case 1:
3141
{
3142
3143
break;
3144
}
3145
case 2:
3146
{
3147
3148
break;
3149
}
3150
}
3151
3152
}
3153
3154
/// Tabulate the coordinates of all dofs on a cell
3155
virtual void tabulate_coordinates(double** coordinates,
3156
const ufc::cell& c) const
3157
{
3158
const double * const * x = c.coordinates;
3159
3160
coordinates[0][0] = x[0][0];
3161
coordinates[0][1] = x[0][1];
3162
coordinates[1][0] = x[1][0];
3163
coordinates[1][1] = x[1][1];
3164
coordinates[2][0] = x[2][0];
3165
coordinates[2][1] = x[2][1];
3166
coordinates[3][0] = x[0][0];
3167
coordinates[3][1] = x[0][1];
3168
coordinates[4][0] = x[1][0];
3169
coordinates[4][1] = x[1][1];
3170
coordinates[5][0] = x[2][0];
3171
coordinates[5][1] = x[2][1];
3172
coordinates[6][0] = x[0][0];
3173
coordinates[6][1] = x[0][1];
3174
coordinates[7][0] = x[1][0];
3175
coordinates[7][1] = x[1][1];
3176
coordinates[8][0] = x[2][0];
3177
coordinates[8][1] = x[2][1];
3178
}
3179
3180
/// Return the number of sub dof maps (for a mixed element)
3181
virtual unsigned int num_sub_dof_maps() const
3182
{
3183
return 2;
3184
}
3185
3186
/// Create a new dof_map for sub dof map i (for a mixed element)
3187
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
3188
{
3189
switch (i)
3190
{
3191
case 0:
3192
{
3193
return new stabilisedstokes_dof_map_1();
3194
break;
3195
}
3196
case 1:
3197
{
3198
return new stabilisedstokes_dof_map_0();
3199
break;
3200
}
3201
}
3202
3203
return 0;
3204
}
3205
3206
};
3207
3208
/// This class defines the interface for the tabulation of the cell
3209
/// tensor corresponding to the local contribution to a form from
3210
/// the integral over a cell.
3211
3212
class stabilisedstokes_cell_integral_0_0: public ufc::cell_integral
3213
{
3214
public:
3215
3216
/// Constructor
3217
stabilisedstokes_cell_integral_0_0() : ufc::cell_integral()
3218
{
3219
// Do nothing
3220
}
3221
3222
/// Destructor
3223
virtual ~stabilisedstokes_cell_integral_0_0()
3224
{
3225
// Do nothing
3226
}
3227
3228
/// Tabulate the tensor for the contribution from a local cell
3229
virtual void tabulate_tensor(double* A,
3230
const double * const * w,
3231
const ufc::cell& c) const
3232
{
3233
// Extract vertex coordinates
3234
const double * const * x = c.coordinates;
3235
3236
// Compute Jacobian of affine map from reference cell
3237
const double J_00 = x[1][0] - x[0][0];
3238
const double J_01 = x[2][0] - x[0][0];
3239
const double J_10 = x[1][1] - x[0][1];
3240
const double J_11 = x[2][1] - x[0][1];
3241
3242
// Compute determinant of Jacobian
3243
double detJ = J_00*J_11 - J_01*J_10;
3244
3245
// Compute inverse of Jacobian
3246
const double K_00 = J_11 / detJ;
3247
const double K_01 = -J_01 / detJ;
3248
const double K_10 = -J_10 / detJ;
3249
const double K_11 = J_00 / detJ;
3250
3251
// Set scale factor
3252
const double det = std::abs(detJ);
3253
3254
// Array of quadrature weights
3255
static const double W4[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
3256
// Quadrature points on the UFC reference element: (0.17855873, 0.15505103), (0.07503111, 0.64494897), (0.66639025, 0.15505103), (0.28001992, 0.64494897)
3257
3258
// Value of basis functions at quadrature points.
3259
static const double FE0[4][3] = \
3260
{{0.66639025, 0.17855873, 0.15505103},
3261
{0.28001992, 0.07503111, 0.64494897},
3262
{0.17855873, 0.66639025, 0.15505103},
3263
{0.07503111, 0.28001992, 0.64494897}};
3264
3265
static const double FE1_C0_D01[4][9] = \
3266
{{-1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3267
{-1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3268
{-1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3269
{-1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3270
3271
static const double FE1_C0_D10[4][9] = \
3272
{{-1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3273
{-1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3274
{-1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3275
{-1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3276
3277
static const double FE1_C1_D01[4][9] = \
3278
{{0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
3279
{0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
3280
{0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
3281
{0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000}};
3282
3283
static const double FE1_C1_D10[4][9] = \
3284
{{0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3285
{0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3286
{0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3287
{0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3288
3289
static const double FE1_C2[4][9] = \
3290
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.66639025, 0.17855873, 0.15505103},
3291
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.28001992, 0.07503111, 0.64494897},
3292
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.17855873, 0.66639025, 0.15505103},
3293
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.07503111, 0.28001992, 0.64494897}};
3294
3295
static const double FE1_C2_D01[4][9] = \
3296
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000},
3297
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000},
3298
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000},
3299
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000}};
3300
3301
static const double FE1_C2_D10[4][9] = \
3302
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000},
3303
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000},
3304
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000},
3305
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000}};
3306
3307
for (unsigned int r = 0; r < 81; r++)
3308
{
3309
A[r] = 0.00000000;
3310
}// end loop over 'r'
3311
3312
// Compute element tensor using UFL quadrature representation
3313
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
3314
3315
// Loop quadrature points for integral
3316
// Number of operations to compute element tensor for following IP loop = 23352
3317
for (unsigned int ip = 0; ip < 4; ip++)
3318
{
3319
3320
// Coefficient declarations
3321
double F0 = 0.00000000;
3322
3323
// Total number of operations to compute function values = 6
3324
for (unsigned int r = 0; r < 3; r++)
3325
{
3326
F0 += FE0[ip][r]*w[0][r];
3327
}// end loop over 'r'
3328
3329
// Number of operations for primary indices: 5832
3330
for (unsigned int j = 0; j < 9; j++)
3331
{
3332
for (unsigned int k = 0; k < 9; k++)
3333
{
3334
// Number of operations to compute entry: 72
3335
A[j*9 + k] += ((((((K_01*FE1_C2_D10[ip][j] + K_11*FE1_C2_D01[ip][j]))*((K_01*FE1_C2_D10[ip][k] + K_11*FE1_C2_D01[ip][k])) + ((K_00*FE1_C2_D10[ip][j] + K_10*FE1_C2_D01[ip][j]))*((K_00*FE1_C2_D10[ip][k] + K_10*FE1_C2_D01[ip][k]))))*F0*0.20000000*F0 + ((((K_00*FE1_C0_D10[ip][k] + K_10*FE1_C0_D01[ip][k]) + (K_01*FE1_C1_D10[ip][k] + K_11*FE1_C1_D01[ip][k])))*FE1_C2[ip][j] + (((((K_00*FE1_C0_D10[ip][j] + K_10*FE1_C0_D01[ip][j]) + (K_01*FE1_C1_D10[ip][j] + K_11*FE1_C1_D01[ip][j])))*FE1_C2[ip][k])*(-1.00000000) + ((((K_01*FE1_C1_D10[ip][j] + K_11*FE1_C1_D01[ip][j]))*((K_01*FE1_C1_D10[ip][k] + K_11*FE1_C1_D01[ip][k])) + ((K_01*FE1_C0_D10[ip][j] + K_11*FE1_C0_D01[ip][j]))*((K_01*FE1_C0_D10[ip][k] + K_11*FE1_C0_D01[ip][k]))) + (((K_00*FE1_C1_D10[ip][j] + K_10*FE1_C1_D01[ip][j]))*((K_00*FE1_C1_D10[ip][k] + K_10*FE1_C1_D01[ip][k])) + ((K_00*FE1_C0_D10[ip][j] + K_10*FE1_C0_D01[ip][j]))*((K_00*FE1_C0_D10[ip][k] + K_10*FE1_C0_D01[ip][k]))))))))*W4[ip]*det;
3336
}// end loop over 'k'
3337
}// end loop over 'j'
3338
}// end loop over 'ip'
3339
}
3340
3341
};
3342
3343
/// This class defines the interface for the tabulation of the cell
3344
/// tensor corresponding to the local contribution to a form from
3345
/// the integral over a cell.
3346
3347
class stabilisedstokes_cell_integral_1_0: public ufc::cell_integral
3348
{
3349
public:
3350
3351
/// Constructor
3352
stabilisedstokes_cell_integral_1_0() : ufc::cell_integral()
3353
{
3354
// Do nothing
3355
}
3356
3357
/// Destructor
3358
virtual ~stabilisedstokes_cell_integral_1_0()
3359
{
3360
// Do nothing
3361
}
3362
3363
/// Tabulate the tensor for the contribution from a local cell
3364
virtual void tabulate_tensor(double* A,
3365
const double * const * w,
3366
const ufc::cell& c) const
3367
{
3368
// Extract vertex coordinates
3369
const double * const * x = c.coordinates;
3370
3371
// Compute Jacobian of affine map from reference cell
3372
const double J_00 = x[1][0] - x[0][0];
3373
const double J_01 = x[2][0] - x[0][0];
3374
const double J_10 = x[1][1] - x[0][1];
3375
const double J_11 = x[2][1] - x[0][1];
3376
3377
// Compute determinant of Jacobian
3378
double detJ = J_00*J_11 - J_01*J_10;
3379
3380
// Compute inverse of Jacobian
3381
const double K_00 = J_11 / detJ;
3382
const double K_01 = -J_01 / detJ;
3383
const double K_10 = -J_10 / detJ;
3384
const double K_11 = J_00 / detJ;
3385
3386
// Set scale factor
3387
const double det = std::abs(detJ);
3388
3389
// Array of quadrature weights
3390
static const double W4[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
3391
// Quadrature points on the UFC reference element: (0.17855873, 0.15505103), (0.07503111, 0.64494897), (0.66639025, 0.15505103), (0.28001992, 0.64494897)
3392
3393
// Value of basis functions at quadrature points.
3394
static const double FE0[4][3] = \
3395
{{0.66639025, 0.17855873, 0.15505103},
3396
{0.28001992, 0.07503111, 0.64494897},
3397
{0.17855873, 0.66639025, 0.15505103},
3398
{0.07503111, 0.28001992, 0.64494897}};
3399
3400
static const double FE1_C0[4][6] = \
3401
{{0.66639025, 0.17855873, 0.15505103, 0.00000000, 0.00000000, 0.00000000},
3402
{0.28001992, 0.07503111, 0.64494897, 0.00000000, 0.00000000, 0.00000000},
3403
{0.17855873, 0.66639025, 0.15505103, 0.00000000, 0.00000000, 0.00000000},
3404
{0.07503111, 0.28001992, 0.64494897, 0.00000000, 0.00000000, 0.00000000}};
3405
3406
static const double FE1_C1[4][6] = \
3407
{{0.00000000, 0.00000000, 0.00000000, 0.66639025, 0.17855873, 0.15505103},
3408
{0.00000000, 0.00000000, 0.00000000, 0.28001992, 0.07503111, 0.64494897},
3409
{0.00000000, 0.00000000, 0.00000000, 0.17855873, 0.66639025, 0.15505103},
3410
{0.00000000, 0.00000000, 0.00000000, 0.07503111, 0.28001992, 0.64494897}};
3411
3412
static const double FE2_C0[4][9] = \
3413
{{0.66639025, 0.17855873, 0.15505103, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3414
{0.28001992, 0.07503111, 0.64494897, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3415
{0.17855873, 0.66639025, 0.15505103, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
3416
{0.07503111, 0.28001992, 0.64494897, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3417
3418
static const double FE2_C1[4][9] = \
3419
{{0.00000000, 0.00000000, 0.00000000, 0.66639025, 0.17855873, 0.15505103, 0.00000000, 0.00000000, 0.00000000},
3420
{0.00000000, 0.00000000, 0.00000000, 0.28001992, 0.07503111, 0.64494897, 0.00000000, 0.00000000, 0.00000000},
3421
{0.00000000, 0.00000000, 0.00000000, 0.17855873, 0.66639025, 0.15505103, 0.00000000, 0.00000000, 0.00000000},
3422
{0.00000000, 0.00000000, 0.00000000, 0.07503111, 0.28001992, 0.64494897, 0.00000000, 0.00000000, 0.00000000}};
3423
3424
static const double FE2_C2_D01[4][9] = \
3425
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000},
3426
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000},
3427
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000},
3428
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 0.00000000, 1.00000000}};
3429
3430
static const double FE2_C2_D10[4][9] = \
3431
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000},
3432
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000},
3433
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000},
3434
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, -1.00000000, 1.00000000, 0.00000000}};
3435
3436
for (unsigned int r = 0; r < 9; r++)
3437
{
3438
A[r] = 0.00000000;
3439
}// end loop over 'r'
3440
3441
// Compute element tensor using UFL quadrature representation
3442
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
3443
3444
// Loop quadrature points for integral
3445
// Number of operations to compute element tensor for following IP loop = 840
3446
for (unsigned int ip = 0; ip < 4; ip++)
3447
{
3448
3449
// Coefficient declarations
3450
double F0 = 0.00000000;
3451
double F1 = 0.00000000;
3452
double F2 = 0.00000000;
3453
3454
// Total number of operations to compute function values = 6
3455
for (unsigned int r = 0; r < 3; r++)
3456
{
3457
F1 += FE0[ip][r]*w[1][r];
3458
}// end loop over 'r'
3459
3460
// Total number of operations to compute function values = 24
3461
for (unsigned int r = 0; r < 6; r++)
3462
{
3463
F0 += FE1_C0[ip][r]*w[0][r];
3464
F2 += FE1_C1[ip][r]*w[0][r];
3465
}// end loop over 'r'
3466
3467
// Number of operations for primary indices: 180
3468
for (unsigned int j = 0; j < 9; j++)
3469
{
3470
// Number of operations to compute entry: 20
3471
A[j] += ((((((K_00*FE2_C2_D10[ip][j] + K_10*FE2_C2_D01[ip][j]))*F1*0.20000000*F1 + FE2_C0[ip][j]))*F0 + ((FE2_C1[ip][j] + ((K_01*FE2_C2_D10[ip][j] + K_11*FE2_C2_D01[ip][j]))*F1*0.20000000*F1))*F2))*W4[ip]*det;
3472
}// end loop over 'j'
3473
}// end loop over 'ip'
3474
}
3475
3476
};
3477
3478
/// This class defines the interface for the assembly of the global
3479
/// tensor corresponding to a form with r + n arguments, that is, a
3480
/// mapping
3481
///
3482
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3483
///
3484
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3485
/// global tensor A is defined by
3486
///
3487
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3488
///
3489
/// where each argument Vj represents the application to the
3490
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3491
/// fixed functions (coefficients).
3492
3493
class stabilisedstokes_form_0: public ufc::form
3494
{
3495
public:
3496
3497
/// Constructor
3498
stabilisedstokes_form_0() : ufc::form()
3499
{
3500
// Do nothing
3501
}
3502
3503
/// Destructor
3504
virtual ~stabilisedstokes_form_0()
3505
{
3506
// Do nothing
3507
}
3508
3509
/// Return a string identifying the form
3510
virtual const char* signature() const
3511
{
3512
return "Form([Integral(Sum(Product(IndexSum(Product(Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((FixedIndex(2),), {})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((FixedIndex(2),), {})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Product(FloatValue(0.20000000000000001, (), (), {}), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0)))), Sum(Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 1), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 1), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((Index(3),), {Index(3): 2})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(2),), {FixedIndex(2): 3}))), Sum(IndexSum(IndexSum(Product(Indexed(ComponentTensor(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((Index(4),), {Index(4): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((Index(4),), {Index(4): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(5),), {Index(5): 2})), MultiIndex((Index(5), Index(4)), {Index(4): 2, Index(5): 2})), MultiIndex((Index(6), Index(7)), {Index(7): 2, Index(6): 2})), Indexed(ComponentTensor(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 1), MultiIndex((Index(8),), {Index(8): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 1), MultiIndex((Index(8),), {Index(8): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(9),), {Index(9): 2})), MultiIndex((Index(9), Index(8)), {Index(8): 2, Index(9): 2})), MultiIndex((Index(6), Index(7)), {Index(7): 2, Index(6): 2}))), MultiIndex((Index(6),), {Index(6): 2})), MultiIndex((Index(7),), {Index(7): 2})), Product(IntValue(-1, (), (), {}), Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((Index(10),), {Index(10): 2})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 1), MultiIndex((FixedIndex(2),), {FixedIndex(2): 3}))))))), Measure('cell', 0, None))])";
3513
}
3514
3515
/// Return the rank of the global tensor (r)
3516
virtual unsigned int rank() const
3517
{
3518
return 2;
3519
}
3520
3521
/// Return the number of coefficients (n)
3522
virtual unsigned int num_coefficients() const
3523
{
3524
return 1;
3525
}
3526
3527
/// Return the number of cell integrals
3528
virtual unsigned int num_cell_integrals() const
3529
{
3530
return 1;
3531
}
3532
3533
/// Return the number of exterior facet integrals
3534
virtual unsigned int num_exterior_facet_integrals() const
3535
{
3536
return 0;
3537
}
3538
3539
/// Return the number of interior facet integrals
3540
virtual unsigned int num_interior_facet_integrals() const
3541
{
3542
return 0;
3543
}
3544
3545
/// Create a new finite element for argument function i
3546
virtual ufc::finite_element* create_finite_element(unsigned int i) const
3547
{
3548
switch (i)
3549
{
3550
case 0:
3551
{
3552
return new stabilisedstokes_finite_element_2();
3553
break;
3554
}
3555
case 1:
3556
{
3557
return new stabilisedstokes_finite_element_2();
3558
break;
3559
}
3560
case 2:
3561
{
3562
return new stabilisedstokes_finite_element_0();
3563
break;
3564
}
3565
}
3566
3567
return 0;
3568
}
3569
3570
/// Create a new dof map for argument function i
3571
virtual ufc::dof_map* create_dof_map(unsigned int i) const
3572
{
3573
switch (i)
3574
{
3575
case 0:
3576
{
3577
return new stabilisedstokes_dof_map_2();
3578
break;
3579
}
3580
case 1:
3581
{
3582
return new stabilisedstokes_dof_map_2();
3583
break;
3584
}
3585
case 2:
3586
{
3587
return new stabilisedstokes_dof_map_0();
3588
break;
3589
}
3590
}
3591
3592
return 0;
3593
}
3594
3595
/// Create a new cell integral on sub domain i
3596
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
3597
{
3598
switch (i)
3599
{
3600
case 0:
3601
{
3602
return new stabilisedstokes_cell_integral_0_0();
3603
break;
3604
}
3605
}
3606
3607
return 0;
3608
}
3609
3610
/// Create a new exterior facet integral on sub domain i
3611
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
3612
{
3613
return 0;
3614
}
3615
3616
/// Create a new interior facet integral on sub domain i
3617
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
3618
{
3619
return 0;
3620
}
3621
3622
};
3623
3624
/// This class defines the interface for the assembly of the global
3625
/// tensor corresponding to a form with r + n arguments, that is, a
3626
/// mapping
3627
///
3628
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3629
///
3630
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3631
/// global tensor A is defined by
3632
///
3633
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3634
///
3635
/// where each argument Vj represents the application to the
3636
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3637
/// fixed functions (coefficients).
3638
3639
class stabilisedstokes_form_1: public ufc::form
3640
{
3641
public:
3642
3643
/// Constructor
3644
stabilisedstokes_form_1() : ufc::form()
3645
{
3646
// Do nothing
3647
}
3648
3649
/// Destructor
3650
virtual ~stabilisedstokes_form_1()
3651
{
3652
// Do nothing
3653
}
3654
3655
/// Return a string identifying the form
3656
virtual const char* signature() const
3657
{
3658
return "Form([Integral(IndexSum(Product(Indexed(Coefficient(VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), 0), MultiIndex((Index(0),), {Index(0): 2})), Indexed(Sum(ComponentTensor(Product(Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((FixedIndex(2),), {})), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((Index(2),), {Index(2): 2})), Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1), Product(FloatValue(0.20000000000000001, (), (), {}), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1)))), MultiIndex((Index(2),), {Index(2): 2})), ListTensor(Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(0),), {FixedIndex(0): 3})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 1, 2), FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(1),), {FixedIndex(1): 3})))), MultiIndex((Index(0),), {Index(0): 2}))), MultiIndex((Index(0),), {Index(0): 2})), Measure('cell', 0, None))])";
3659
}
3660
3661
/// Return the rank of the global tensor (r)
3662
virtual unsigned int rank() const
3663
{
3664
return 1;
3665
}
3666
3667
/// Return the number of coefficients (n)
3668
virtual unsigned int num_coefficients() const
3669
{
3670
return 2;
3671
}
3672
3673
/// Return the number of cell integrals
3674
virtual unsigned int num_cell_integrals() const
3675
{
3676
return 1;
3677
}
3678
3679
/// Return the number of exterior facet integrals
3680
virtual unsigned int num_exterior_facet_integrals() const
3681
{
3682
return 0;
3683
}
3684
3685
/// Return the number of interior facet integrals
3686
virtual unsigned int num_interior_facet_integrals() const
3687
{
3688
return 0;
3689
}
3690
3691
/// Create a new finite element for argument function i
3692
virtual ufc::finite_element* create_finite_element(unsigned int i) const
3693
{
3694
switch (i)
3695
{
3696
case 0:
3697
{
3698
return new stabilisedstokes_finite_element_2();
3699
break;
3700
}
3701
case 1:
3702
{
3703
return new stabilisedstokes_finite_element_1();
3704
break;
3705
}
3706
case 2:
3707
{
3708
return new stabilisedstokes_finite_element_0();
3709
break;
3710
}
3711
}
3712
3713
return 0;
3714
}
3715
3716
/// Create a new dof map for argument function i
3717
virtual ufc::dof_map* create_dof_map(unsigned int i) const
3718
{
3719
switch (i)
3720
{
3721
case 0:
3722
{
3723
return new stabilisedstokes_dof_map_2();
3724
break;
3725
}
3726
case 1:
3727
{
3728
return new stabilisedstokes_dof_map_1();
3729
break;
3730
}
3731
case 2:
3732
{
3733
return new stabilisedstokes_dof_map_0();
3734
break;
3735
}
3736
}
3737
3738
return 0;
3739
}
3740
3741
/// Create a new cell integral on sub domain i
3742
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
3743
{
3744
switch (i)
3745
{
3746
case 0:
3747
{
3748
return new stabilisedstokes_cell_integral_1_0();
3749
break;
3750
}
3751
}
3752
3753
return 0;
3754
}
3755
3756
/// Create a new exterior facet integral on sub domain i
3757
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
3758
{
3759
return 0;
3760
}
3761
3762
/// Create a new interior facet integral on sub domain i
3763
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
3764
{
3765
return 0;
3766
}
3767
3768
};
3769
3770
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1