~ubuntu-branches/ubuntu/utopic/ffc/utopic
» Revision
17
: test/regression/references/r_quadrature_O/Stokes.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_quadrature_O/Stokes.h
- browse files at revision 17
- compare with another revision
- download diff
- download tarball
- view history from revision 17
- Committer: Package Import Robot
- Author(s): Johannes Ring
- Date: 2013-06-26 14:48:32 UTC
- mfrom: (1.1.13)
- Revision ID: package-import@ubuntu.com-20130626144832-1xd8htax4s3utybz
Tags: 1.2.0-1
* New upstream release.
* debian/control:
- Bump required version for python-ufc, python-fiat, python-instant
and python-ufl in Depends field.
- Bump Standards-Version to 3.9.4.
- Remove DM-Upload-Allowed field.
- Bump required debhelper version in Build-Depends.
- Remove cdbs from Build-Depends.
- Use canonical URIs for Vcs-* fields.
* debian/compat: Bump to compatibility level 9.
* debian/rules: Rewrite for debhelper (drop cdbs).
* debian/control:
- Bump required version for python-ufc, python-fiat, python-instant
and python-ufl in Depends field.
- Bump Standards-Version to 3.9.4.
- Remove DM-Upload-Allowed field.
- Bump required debhelper version in Build-Depends.
- Remove cdbs from Build-Depends.
- Use canonical URIs for Vcs-* fields.
* debian/compat: Bump to compatibility level 9.
* debian/rules: Rewrite for debhelper (drop cdbs).
- files added:
- .pc
- ffc/uflacsrepr
- sandbox/features
- sandbox/multidomain
- files modified:
- ChangeLog
added
removed
test/regression/references/r_quadrature_O/Stokes.h
1
// This code conforms with the UFC specification version 2.0.5
2
// and was automatically generated by FFC version 1.0.0.
1
// This code conforms with the UFC specification version 2.2.0
2
// and was automatically generated by FFC version 1.2.0.
3
3
//
4
4
// This code was generated with the following parameters:
5
5
//
52
52
/// Return a string identifying the finite element
53
53
virtual const char* signature() const
54
54
{
55
return "FiniteElement('Lagrange', Cell('triangle', Space(2)), 2, None)";
55
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
56
56
}
57
57
58
58
/// Return the cell shape
62
62
}
63
63
64
64
/// Return the topological dimension of the cell shape
65
virtual unsigned int topological_dimension() const
65
virtual std::size_t topological_dimension() const
66
66
{
67
67
return 2;
68
68
}
69
69
70
70
/// Return the geometric dimension of the cell shape
71
virtual unsigned int geometric_dimension() const
71
virtual std::size_t geometric_dimension() const
72
72
{
73
73
return 2;
74
74
}
75
75
76
76
/// Return the dimension of the finite element function space
77
virtual unsigned int space_dimension() const
77
virtual std::size_t space_dimension() const
78
78
{
79
79
return 6;
80
80
}
81
81
82
82
/// Return the rank of the value space
83
virtual unsigned int value_rank() const
83
virtual std::size_t value_rank() const
84
84
{
85
85
return 0;
86
86
}
87
87
88
88
/// Return the dimension of the value space for axis i
89
virtual unsigned int value_dimension(unsigned int i) const
89
virtual std::size_t value_dimension(std::size_t i) const
90
90
{
91
91
return 1;
92
92
}
93
93
94
/// Evaluate basis function i at given point in cell
95
virtual void evaluate_basis(unsigned int i,
94
/// Evaluate basis function i at given point x in cell
95
virtual void evaluate_basis(std::size_t i,
96
96
double* values,
97
const double* coordinates,
98
const ufc::cell& c) const
97
const double* x,
98
const double* vertex_coordinates,
99
int cell_orientation) const
99
100
{
100
// Extract vertex coordinates
101
const double * const * x = c.coordinates;
102
103
// Compute Jacobian of affine map from reference cell
104
const double J_00 = x[1][0] - x[0][0];
105
const double J_01 = x[2][0] - x[0][0];
106
const double J_10 = x[1][1] - x[0][1];
107
const double J_11 = x[2][1] - x[0][1];
108
109
// Compute determinant of Jacobian
110
double detJ = J_00*J_11 - J_01*J_10;
111
112
// Compute inverse of Jacobian
101
// Compute Jacobian
102
double J[4];
103
compute_jacobian_triangle_2d(J, vertex_coordinates);
104
105
// Compute Jacobian inverse and determinant
106
double K[4];
107
double detJ;
108
compute_jacobian_inverse_triangle_2d(K, detJ, J);
109
113
110
114
111
// Compute constants
115
const double C0 = x[1][0] + x[2][0];
116
const double C1 = x[1][1] + x[2][1];
112
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
113
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
117
114
118
115
// Get coordinates and map to the reference (FIAT) element
119
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
120
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
116
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
117
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
121
118
122
// Reset values.
119
// Reset values
123
120
*values = 0.0;
124
121
switch (i)
125
122
{
126
123
case 0:
127
124
{
128
125
129
// Array of basisvalues.
126
// Array of basisvalues
130
127
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
131
128
132
// Declare helper variables.
129
// Declare helper variables
133
130
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
134
131
double tmp1 = (1.0 - Y)/2.0;
135
132
double tmp2 = tmp1*tmp1;
136
133
137
// Compute basisvalues.
134
// Compute basisvalues
138
135
basisvalues[0] = 1.0;
139
136
basisvalues[1] = tmp0;
140
137
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
148
145
basisvalues[4] *= std::sqrt(4.5);
149
146
basisvalues[3] *= std::sqrt(7.5);
150
147
151
// Table(s) of coefficients.
148
// Table(s) of coefficients
152
149
static const double coefficients0[6] = \
153
150
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
154
151
155
// Compute value(s).
152
// Compute value(s)
156
153
for (unsigned int r = 0; r < 6; r++)
157
154
{
158
155
*values += coefficients0[r]*basisvalues[r];
162
159
case 1:
163
160
{
164
161
165
// Array of basisvalues.
162
// Array of basisvalues
166
163
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
167
164
168
// Declare helper variables.
165
// Declare helper variables
169
166
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
170
167
double tmp1 = (1.0 - Y)/2.0;
171
168
double tmp2 = tmp1*tmp1;
172
169
173
// Compute basisvalues.
170
// Compute basisvalues
174
171
basisvalues[0] = 1.0;
175
172
basisvalues[1] = tmp0;
176
173
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
184
181
basisvalues[4] *= std::sqrt(4.5);
185
182
basisvalues[3] *= std::sqrt(7.5);
186
183
187
// Table(s) of coefficients.
184
// Table(s) of coefficients
188
185
static const double coefficients0[6] = \
189
186
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
190
187
191
// Compute value(s).
188
// Compute value(s)
192
189
for (unsigned int r = 0; r < 6; r++)
193
190
{
194
191
*values += coefficients0[r]*basisvalues[r];
198
195
case 2:
199
196
{
200
197
201
// Array of basisvalues.
198
// Array of basisvalues
202
199
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
203
200
204
// Declare helper variables.
201
// Declare helper variables
205
202
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
206
203
double tmp1 = (1.0 - Y)/2.0;
207
204
double tmp2 = tmp1*tmp1;
208
205
209
// Compute basisvalues.
206
// Compute basisvalues
210
207
basisvalues[0] = 1.0;
211
208
basisvalues[1] = tmp0;
212
209
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
220
217
basisvalues[4] *= std::sqrt(4.5);
221
218
basisvalues[3] *= std::sqrt(7.5);
222
219
223
// Table(s) of coefficients.
220
// Table(s) of coefficients
224
221
static const double coefficients0[6] = \
225
222
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
226
223
227
// Compute value(s).
224
// Compute value(s)
228
225
for (unsigned int r = 0; r < 6; r++)
229
226
{
230
227
*values += coefficients0[r]*basisvalues[r];
234
231
case 3:
235
232
{
236
233
237
// Array of basisvalues.
234
// Array of basisvalues
238
235
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
239
236
240
// Declare helper variables.
237
// Declare helper variables
241
238
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
242
239
double tmp1 = (1.0 - Y)/2.0;
243
240
double tmp2 = tmp1*tmp1;
244
241
245
// Compute basisvalues.
242
// Compute basisvalues
246
243
basisvalues[0] = 1.0;
247
244
basisvalues[1] = tmp0;
248
245
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
256
253
basisvalues[4] *= std::sqrt(4.5);
257
254
basisvalues[3] *= std::sqrt(7.5);
258
255
259
// Table(s) of coefficients.
256
// Table(s) of coefficients
260
257
static const double coefficients0[6] = \
261
258
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
262
259
263
// Compute value(s).
260
// Compute value(s)
264
261
for (unsigned int r = 0; r < 6; r++)
265
262
{
266
263
*values += coefficients0[r]*basisvalues[r];
270
267
case 4:
271
268
{
272
269
273
// Array of basisvalues.
270
// Array of basisvalues
274
271
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
275
272
276
// Declare helper variables.
273
// Declare helper variables
277
274
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
278
275
double tmp1 = (1.0 - Y)/2.0;
279
276
double tmp2 = tmp1*tmp1;
280
277
281
// Compute basisvalues.
278
// Compute basisvalues
282
279
basisvalues[0] = 1.0;
283
280
basisvalues[1] = tmp0;
284
281
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
292
289
basisvalues[4] *= std::sqrt(4.5);
293
290
basisvalues[3] *= std::sqrt(7.5);
294
291
295
// Table(s) of coefficients.
292
// Table(s) of coefficients
296
293
static const double coefficients0[6] = \
297
294
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
298
295
299
// Compute value(s).
296
// Compute value(s)
300
297
for (unsigned int r = 0; r < 6; r++)
301
298
{
302
299
*values += coefficients0[r]*basisvalues[r];
306
303
case 5:
307
304
{
308
305
309
// Array of basisvalues.
306
// Array of basisvalues
310
307
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
311
308
312
// Declare helper variables.
309
// Declare helper variables
313
310
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
314
311
double tmp1 = (1.0 - Y)/2.0;
315
312
double tmp2 = tmp1*tmp1;
316
313
317
// Compute basisvalues.
314
// Compute basisvalues
318
315
basisvalues[0] = 1.0;
319
316
basisvalues[1] = tmp0;
320
317
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
328
325
basisvalues[4] *= std::sqrt(4.5);
329
326
basisvalues[3] *= std::sqrt(7.5);
330
327
331
// Table(s) of coefficients.
328
// Table(s) of coefficients
332
329
static const double coefficients0[6] = \
333
330
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
334
331
335
// Compute value(s).
332
// Compute value(s)
336
333
for (unsigned int r = 0; r < 6; r++)
337
334
{
338
335
*values += coefficients0[r]*basisvalues[r];
343
340
344
341
}
345
342
346
/// Evaluate all basis functions at given point in cell
343
/// Evaluate all basis functions at given point x in cell
347
344
virtual void evaluate_basis_all(double* values,
348
const double* coordinates,
349
const ufc::cell& c) const
345
const double* x,
346
const double* vertex_coordinates,
347
int cell_orientation) const
350
348
{
351
349
// Helper variable to hold values of a single dof.
352
350
double dof_values = 0.0;
353
351
354
// Loop dofs and call evaluate_basis.
352
// Loop dofs and call evaluate_basis
355
353
for (unsigned int r = 0; r < 6; r++)
356
354
{
357
evaluate_basis(r, &dof_values, coordinates, c);
355
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
358
356
values[r] = dof_values;
359
357
}// end loop over 'r'
360
358
}
361
359
362
/// Evaluate order n derivatives of basis function i at given point in cell
363
virtual void evaluate_basis_derivatives(unsigned int i,
364
unsigned int n,
360
/// Evaluate order n derivatives of basis function i at given point x in cell
361
virtual void evaluate_basis_derivatives(std::size_t i,
362
std::size_t n,
365
363
double* values,
366
const double* coordinates,
367
const ufc::cell& c) const
364
const double* x,
365
const double* vertex_coordinates,
366
int cell_orientation) const
368
367
{
369
// Extract vertex coordinates
370
const double * const * x = c.coordinates;
371
372
// Compute Jacobian of affine map from reference cell
373
const double J_00 = x[1][0] - x[0][0];
374
const double J_01 = x[2][0] - x[0][0];
375
const double J_10 = x[1][1] - x[0][1];
376
const double J_11 = x[2][1] - x[0][1];
377
378
// Compute determinant of Jacobian
379
double detJ = J_00*J_11 - J_01*J_10;
380
381
// Compute inverse of Jacobian
382
const double K_00 = J_11 / detJ;
383
const double K_01 = -J_01 / detJ;
384
const double K_10 = -J_10 / detJ;
385
const double K_11 = J_00 / detJ;
368
// Compute Jacobian
369
double J[4];
370
compute_jacobian_triangle_2d(J, vertex_coordinates);
371
372
// Compute Jacobian inverse and determinant
373
double K[4];
374
double detJ;
375
compute_jacobian_inverse_triangle_2d(K, detJ, J);
376
386
377
387
378
// Compute constants
388
const double C0 = x[1][0] + x[2][0];
389
const double C1 = x[1][1] + x[2][1];
379
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
380
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
390
381
391
382
// Get coordinates and map to the reference (FIAT) element
392
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
393
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
383
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
384
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
394
385
395
386
// Compute number of derivatives.
396
387
unsigned int num_derivatives = 1;
427
418
}
428
419
429
420
// Compute inverse of Jacobian
430
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
421
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
431
422
432
423
// Declare transformation matrix
433
424
// Declare pointer to two dimensional array and initialise
461
452
case 0:
462
453
{
463
454
464
// Array of basisvalues.
455
// Array of basisvalues
465
456
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
466
457
467
// Declare helper variables.
458
// Declare helper variables
468
459
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
469
460
double tmp1 = (1.0 - Y)/2.0;
470
461
double tmp2 = tmp1*tmp1;
471
462
472
// Compute basisvalues.
463
// Compute basisvalues
473
464
basisvalues[0] = 1.0;
474
465
basisvalues[1] = tmp0;
475
466
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
483
474
basisvalues[4] *= std::sqrt(4.5);
484
475
basisvalues[3] *= std::sqrt(7.5);
485
476
486
// Table(s) of coefficients.
477
// Table(s) of coefficients
487
478
static const double coefficients0[6] = \
488
479
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
489
480
627
618
case 1:
628
619
{
629
620
630
// Array of basisvalues.
621
// Array of basisvalues
631
622
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
632
623
633
// Declare helper variables.
624
// Declare helper variables
634
625
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
635
626
double tmp1 = (1.0 - Y)/2.0;
636
627
double tmp2 = tmp1*tmp1;
637
628
638
// Compute basisvalues.
629
// Compute basisvalues
639
630
basisvalues[0] = 1.0;
640
631
basisvalues[1] = tmp0;
641
632
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
649
640
basisvalues[4] *= std::sqrt(4.5);
650
641
basisvalues[3] *= std::sqrt(7.5);
651
642
652
// Table(s) of coefficients.
643
// Table(s) of coefficients
653
644
static const double coefficients0[6] = \
654
645
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
655
646
793
784
case 2:
794
785
{
795
786
796
// Array of basisvalues.
787
// Array of basisvalues
797
788
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
798
789
799
// Declare helper variables.
790
// Declare helper variables
800
791
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
801
792
double tmp1 = (1.0 - Y)/2.0;
802
793
double tmp2 = tmp1*tmp1;
803
794
804
// Compute basisvalues.
795
// Compute basisvalues
805
796
basisvalues[0] = 1.0;
806
797
basisvalues[1] = tmp0;
807
798
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
815
806
basisvalues[4] *= std::sqrt(4.5);
816
807
basisvalues[3] *= std::sqrt(7.5);
817
808
818
// Table(s) of coefficients.
809
// Table(s) of coefficients
819
810
static const double coefficients0[6] = \
820
811
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
821
812
959
950
case 3:
960
951
{
961
952
962
// Array of basisvalues.
953
// Array of basisvalues
963
954
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
964
955
965
// Declare helper variables.
956
// Declare helper variables
966
957
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
967
958
double tmp1 = (1.0 - Y)/2.0;
968
959
double tmp2 = tmp1*tmp1;
969
960
970
// Compute basisvalues.
961
// Compute basisvalues
971
962
basisvalues[0] = 1.0;
972
963
basisvalues[1] = tmp0;
973
964
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
981
972
basisvalues[4] *= std::sqrt(4.5);
982
973
basisvalues[3] *= std::sqrt(7.5);
983
974
984
// Table(s) of coefficients.
975
// Table(s) of coefficients
985
976
static const double coefficients0[6] = \
986
977
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
987
978
1125
1116
case 4:
1126
1117
{
1127
1118
1128
// Array of basisvalues.
1119
// Array of basisvalues
1129
1120
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1130
1121
1131
// Declare helper variables.
1122
// Declare helper variables
1132
1123
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1133
1124
double tmp1 = (1.0 - Y)/2.0;
1134
1125
double tmp2 = tmp1*tmp1;
1135
1126
1136
// Compute basisvalues.
1127
// Compute basisvalues
1137
1128
basisvalues[0] = 1.0;
1138
1129
basisvalues[1] = tmp0;
1139
1130
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1147
1138
basisvalues[4] *= std::sqrt(4.5);
1148
1139
basisvalues[3] *= std::sqrt(7.5);
1149
1140
1150
// Table(s) of coefficients.
1141
// Table(s) of coefficients
1151
1142
static const double coefficients0[6] = \
1152
1143
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
1153
1144
1291
1282
case 5:
1292
1283
{
1293
1284
1294
// Array of basisvalues.
1285
// Array of basisvalues
1295
1286
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1296
1287
1297
// Declare helper variables.
1288
// Declare helper variables
1298
1289
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1299
1290
double tmp1 = (1.0 - Y)/2.0;
1300
1291
double tmp2 = tmp1*tmp1;
1301
1292
1302
// Compute basisvalues.
1293
// Compute basisvalues
1303
1294
basisvalues[0] = 1.0;
1304
1295
basisvalues[1] = tmp0;
1305
1296
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1313
1304
basisvalues[4] *= std::sqrt(4.5);
1314
1305
basisvalues[3] *= std::sqrt(7.5);
1315
1306
1316
// Table(s) of coefficients.
1307
// Table(s) of coefficients
1317
1308
static const double coefficients0[6] = \
1318
1309
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
1319
1310
1458
1449
1459
1450
}
1460
1451
1461
/// Evaluate order n derivatives of all basis functions at given point in cell
1462
virtual void evaluate_basis_derivatives_all(unsigned int n,
1452
/// Evaluate order n derivatives of all basis functions at given point x in cell
1453
virtual void evaluate_basis_derivatives_all(std::size_t n,
1463
1454
double* values,
1464
const double* coordinates,
1465
const ufc::cell& c) const
1455
const double* x,
1456
const double* vertex_coordinates,
1457
int cell_orientation) const
1466
1458
{
1467
1459
// Compute number of derivatives.
1468
1460
unsigned int num_derivatives = 1;
1481
1473
// Loop dofs and call evaluate_basis_derivatives.
1482
1474
for (unsigned int r = 0; r < 6; r++)
1483
1475
{
1484
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
1476
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1485
1477
for (unsigned int s = 0; s < num_derivatives; s++)
1486
1478
{
1487
1479
values[r*num_derivatives + s] = dof_values[s];
1493
1485
}
1494
1486
1495
1487
/// Evaluate linear functional for dof i on the function f
1496
virtual double evaluate_dof(unsigned int i,
1488
virtual double evaluate_dof(std::size_t i,
1497
1489
const ufc::function& f,
1490
const double* vertex_coordinates,
1491
int cell_orientation,
1498
1492
const ufc::cell& c) const
1499
1493
{
1500
// Declare variables for result of evaluation.
1494
// Declare variables for result of evaluation
1501
1495
double vals[1];
1502
1496
1503
// Declare variable for physical coordinates.
1497
// Declare variable for physical coordinates
1504
1498
double y[2];
1505
const double * const * x = c.coordinates;
1506
1499
switch (i)
1507
1500
{
1508
1501
case 0:
1509
1502
{
1510
y[0] = x[0][0];
1511
y[1] = x[0][1];
1503
y[0] = vertex_coordinates[0];
1504
y[1] = vertex_coordinates[1];
1512
1505
f.evaluate(vals, y, c);
1513
1506
return vals[0];
1514
1507
break;
1515
1508
}
1516
1509
case 1:
1517
1510
{
1518
y[0] = x[1][0];
1519
y[1] = x[1][1];
1511
y[0] = vertex_coordinates[2];
1512
y[1] = vertex_coordinates[3];
1520
1513
f.evaluate(vals, y, c);
1521
1514
return vals[0];
1522
1515
break;
1523
1516
}
1524
1517
case 2:
1525
1518
{
1526
y[0] = x[2][0];
1527
y[1] = x[2][1];
1519
y[0] = vertex_coordinates[4];
1520
y[1] = vertex_coordinates[5];
1528
1521
f.evaluate(vals, y, c);
1529
1522
return vals[0];
1530
1523
break;
1531
1524
}
1532
1525
case 3:
1533
1526
{
1534
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
1535
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
1527
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
1528
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1536
1529
f.evaluate(vals, y, c);
1537
1530
return vals[0];
1538
1531
break;
1539
1532
}
1540
1533
case 4:
1541
1534
{
1542
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
1543
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
1535
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
1536
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1544
1537
f.evaluate(vals, y, c);
1545
1538
return vals[0];
1546
1539
break;
1547
1540
}
1548
1541
case 5:
1549
1542
{
1550
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
1551
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
1543
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
1544
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1552
1545
f.evaluate(vals, y, c);
1553
1546
return vals[0];
1554
1547
break;
1561
1554
/// Evaluate linear functionals for all dofs on the function f
1562
1555
virtual void evaluate_dofs(double* values,
1563
1556
const ufc::function& f,
1557
const double* vertex_coordinates,
1558
int cell_orientation,
1564
1559
const ufc::cell& c) const
1565
1560
{
1566
// Declare variables for result of evaluation.
1561
// Declare variables for result of evaluation
1567
1562
double vals[1];
1568
1563
1569
// Declare variable for physical coordinates.
1564
// Declare variable for physical coordinates
1570
1565
double y[2];
1571
const double * const * x = c.coordinates;
1572
y[0] = x[0][0];
1573
y[1] = x[0][1];
1566
y[0] = vertex_coordinates[0];
1567
y[1] = vertex_coordinates[1];
1574
1568
f.evaluate(vals, y, c);
1575
1569
values[0] = vals[0];
1576
y[0] = x[1][0];
1577
y[1] = x[1][1];
1570
y[0] = vertex_coordinates[2];
1571
y[1] = vertex_coordinates[3];
1578
1572
f.evaluate(vals, y, c);
1579
1573
values[1] = vals[0];
1580
y[0] = x[2][0];
1581
y[1] = x[2][1];
1574
y[0] = vertex_coordinates[4];
1575
y[1] = vertex_coordinates[5];
1582
1576
f.evaluate(vals, y, c);
1583
1577
values[2] = vals[0];
1584
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
1585
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
1578
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
1579
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1586
1580
f.evaluate(vals, y, c);
1587
1581
values[3] = vals[0];
1588
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
1589
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
1582
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
1583
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1590
1584
f.evaluate(vals, y, c);
1591
1585
values[4] = vals[0];
1592
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
1593
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
1586
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
1587
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1594
1588
f.evaluate(vals, y, c);
1595
1589
values[5] = vals[0];
1596
1590
}
1598
1592
/// Interpolate vertex values from dof values
1599
1593
virtual void interpolate_vertex_values(double* vertex_values,
1600
1594
const double* dof_values,
1595
const double* vertex_coordinates,
1596
int cell_orientation,
1601
1597
const ufc::cell& c) const
1602
1598
{
1603
1599
// Evaluate function and change variables
1611
1607
const double* xhat,
1612
1608
const ufc::cell& c) const
1613
1609
{
1614
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
1610
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1615
1611
}
1616
1612
1617
1613
/// Map from coordinate x in cell to coordinate xhat in reference cell
1619
1615
const double* x,
1620
1616
const ufc::cell& c) const
1621
1617
{
1622
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
1618
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1623
1619
}
1624
1620
1625
1621
/// Return the number of sub elements (for a mixed element)
1626
virtual unsigned int num_sub_elements() const
1622
virtual std::size_t num_sub_elements() const
1627
1623
{
1628
1624
return 0;
1629
1625
}
1630
1626
1631
1627
/// Create a new finite element for sub element i (for a mixed element)
1632
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1628
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1633
1629
{
1634
1630
return 0;
1635
1631
}
1663
1659
/// Return a string identifying the finite element
1664
1660
virtual const char* signature() const
1665
1661
{
1666
return "VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None)";
1662
return "VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None)";
1667
1663
}
1668
1664
1669
1665
/// Return the cell shape
1673
1669
}
1674
1670
1675
1671
/// Return the topological dimension of the cell shape
1676
virtual unsigned int topological_dimension() const
1672
virtual std::size_t topological_dimension() const
1677
1673
{
1678
1674
return 2;
1679
1675
}
1680
1676
1681
1677
/// Return the geometric dimension of the cell shape
1682
virtual unsigned int geometric_dimension() const
1678
virtual std::size_t geometric_dimension() const
1683
1679
{
1684
1680
return 2;
1685
1681
}
1686
1682
1687
1683
/// Return the dimension of the finite element function space
1688
virtual unsigned int space_dimension() const
1684
virtual std::size_t space_dimension() const
1689
1685
{
1690
1686
return 12;
1691
1687
}
1692
1688
1693
1689
/// Return the rank of the value space
1694
virtual unsigned int value_rank() const
1690
virtual std::size_t value_rank() const
1695
1691
{
1696
1692
return 1;
1697
1693
}
1698
1694
1699
1695
/// Return the dimension of the value space for axis i
1700
virtual unsigned int value_dimension(unsigned int i) const
1696
virtual std::size_t value_dimension(std::size_t i) const
1701
1697
{
1702
1698
switch (i)
1703
1699
{
1711
1707
return 0;
1712
1708
}
1713
1709
1714
/// Evaluate basis function i at given point in cell
1715
virtual void evaluate_basis(unsigned int i,
1710
/// Evaluate basis function i at given point x in cell
1711
virtual void evaluate_basis(std::size_t i,
1716
1712
double* values,
1717
const double* coordinates,
1718
const ufc::cell& c) const
1713
const double* x,
1714
const double* vertex_coordinates,
1715
int cell_orientation) const
1719
1716
{
1720
// Extract vertex coordinates
1721
const double * const * x = c.coordinates;
1722
1723
// Compute Jacobian of affine map from reference cell
1724
const double J_00 = x[1][0] - x[0][0];
1725
const double J_01 = x[2][0] - x[0][0];
1726
const double J_10 = x[1][1] - x[0][1];
1727
const double J_11 = x[2][1] - x[0][1];
1728
1729
// Compute determinant of Jacobian
1730
double detJ = J_00*J_11 - J_01*J_10;
1731
1732
// Compute inverse of Jacobian
1717
// Compute Jacobian
1718
double J[4];
1719
compute_jacobian_triangle_2d(J, vertex_coordinates);
1720
1721
// Compute Jacobian inverse and determinant
1722
double K[4];
1723
double detJ;
1724
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1725
1733
1726
1734
1727
// Compute constants
1735
const double C0 = x[1][0] + x[2][0];
1736
const double C1 = x[1][1] + x[2][1];
1728
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1729
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1737
1730
1738
1731
// Get coordinates and map to the reference (FIAT) element
1739
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1740
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
1732
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1733
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1741
1734
1742
// Reset values.
1735
// Reset values
1743
1736
values[0] = 0.0;
1744
1737
values[1] = 0.0;
1745
1738
switch (i)
1747
1740
case 0:
1748
1741
{
1749
1742
1750
// Array of basisvalues.
1743
// Array of basisvalues
1751
1744
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1752
1745
1753
// Declare helper variables.
1746
// Declare helper variables
1754
1747
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1755
1748
double tmp1 = (1.0 - Y)/2.0;
1756
1749
double tmp2 = tmp1*tmp1;
1757
1750
1758
// Compute basisvalues.
1751
// Compute basisvalues
1759
1752
basisvalues[0] = 1.0;
1760
1753
basisvalues[1] = tmp0;
1761
1754
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1769
1762
basisvalues[4] *= std::sqrt(4.5);
1770
1763
basisvalues[3] *= std::sqrt(7.5);
1771
1764
1772
// Table(s) of coefficients.
1765
// Table(s) of coefficients
1773
1766
static const double coefficients0[6] = \
1774
1767
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
1775
1768
1776
// Compute value(s).
1769
// Compute value(s)
1777
1770
for (unsigned int r = 0; r < 6; r++)
1778
1771
{
1779
1772
values[0] += coefficients0[r]*basisvalues[r];
1783
1776
case 1:
1784
1777
{
1785
1778
1786
// Array of basisvalues.
1779
// Array of basisvalues
1787
1780
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1788
1781
1789
// Declare helper variables.
1782
// Declare helper variables
1790
1783
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1791
1784
double tmp1 = (1.0 - Y)/2.0;
1792
1785
double tmp2 = tmp1*tmp1;
1793
1786
1794
// Compute basisvalues.
1787
// Compute basisvalues
1795
1788
basisvalues[0] = 1.0;
1796
1789
basisvalues[1] = tmp0;
1797
1790
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1805
1798
basisvalues[4] *= std::sqrt(4.5);
1806
1799
basisvalues[3] *= std::sqrt(7.5);
1807
1800
1808
// Table(s) of coefficients.
1801
// Table(s) of coefficients
1809
1802
static const double coefficients0[6] = \
1810
1803
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
1811
1804
1812
// Compute value(s).
1805
// Compute value(s)
1813
1806
for (unsigned int r = 0; r < 6; r++)
1814
1807
{
1815
1808
values[0] += coefficients0[r]*basisvalues[r];
1819
1812
case 2:
1820
1813
{
1821
1814
1822
// Array of basisvalues.
1815
// Array of basisvalues
1823
1816
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1824
1817
1825
// Declare helper variables.
1818
// Declare helper variables
1826
1819
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1827
1820
double tmp1 = (1.0 - Y)/2.0;
1828
1821
double tmp2 = tmp1*tmp1;
1829
1822
1830
// Compute basisvalues.
1823
// Compute basisvalues
1831
1824
basisvalues[0] = 1.0;
1832
1825
basisvalues[1] = tmp0;
1833
1826
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1841
1834
basisvalues[4] *= std::sqrt(4.5);
1842
1835
basisvalues[3] *= std::sqrt(7.5);
1843
1836
1844
// Table(s) of coefficients.
1837
// Table(s) of coefficients
1845
1838
static const double coefficients0[6] = \
1846
1839
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
1847
1840
1848
// Compute value(s).
1841
// Compute value(s)
1849
1842
for (unsigned int r = 0; r < 6; r++)
1850
1843
{
1851
1844
values[0] += coefficients0[r]*basisvalues[r];
1855
1848
case 3:
1856
1849
{
1857
1850
1858
// Array of basisvalues.
1851
// Array of basisvalues
1859
1852
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1860
1853
1861
// Declare helper variables.
1854
// Declare helper variables
1862
1855
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1863
1856
double tmp1 = (1.0 - Y)/2.0;
1864
1857
double tmp2 = tmp1*tmp1;
1865
1858
1866
// Compute basisvalues.
1859
// Compute basisvalues
1867
1860
basisvalues[0] = 1.0;
1868
1861
basisvalues[1] = tmp0;
1869
1862
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1877
1870
basisvalues[4] *= std::sqrt(4.5);
1878
1871
basisvalues[3] *= std::sqrt(7.5);
1879
1872
1880
// Table(s) of coefficients.
1873
// Table(s) of coefficients
1881
1874
static const double coefficients0[6] = \
1882
1875
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
1883
1876
1884
// Compute value(s).
1877
// Compute value(s)
1885
1878
for (unsigned int r = 0; r < 6; r++)
1886
1879
{
1887
1880
values[0] += coefficients0[r]*basisvalues[r];
1891
1884
case 4:
1892
1885
{
1893
1886
1894
// Array of basisvalues.
1887
// Array of basisvalues
1895
1888
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1896
1889
1897
// Declare helper variables.
1890
// Declare helper variables
1898
1891
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1899
1892
double tmp1 = (1.0 - Y)/2.0;
1900
1893
double tmp2 = tmp1*tmp1;
1901
1894
1902
// Compute basisvalues.
1895
// Compute basisvalues
1903
1896
basisvalues[0] = 1.0;
1904
1897
basisvalues[1] = tmp0;
1905
1898
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1913
1906
basisvalues[4] *= std::sqrt(4.5);
1914
1907
basisvalues[3] *= std::sqrt(7.5);
1915
1908
1916
// Table(s) of coefficients.
1909
// Table(s) of coefficients
1917
1910
static const double coefficients0[6] = \
1918
1911
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
1919
1912
1920
// Compute value(s).
1913
// Compute value(s)
1921
1914
for (unsigned int r = 0; r < 6; r++)
1922
1915
{
1923
1916
values[0] += coefficients0[r]*basisvalues[r];
1927
1920
case 5:
1928
1921
{
1929
1922
1930
// Array of basisvalues.
1923
// Array of basisvalues
1931
1924
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1932
1925
1933
// Declare helper variables.
1926
// Declare helper variables
1934
1927
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1935
1928
double tmp1 = (1.0 - Y)/2.0;
1936
1929
double tmp2 = tmp1*tmp1;
1937
1930
1938
// Compute basisvalues.
1931
// Compute basisvalues
1939
1932
basisvalues[0] = 1.0;
1940
1933
basisvalues[1] = tmp0;
1941
1934
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1949
1942
basisvalues[4] *= std::sqrt(4.5);
1950
1943
basisvalues[3] *= std::sqrt(7.5);
1951
1944
1952
// Table(s) of coefficients.
1945
// Table(s) of coefficients
1953
1946
static const double coefficients0[6] = \
1954
1947
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
1955
1948
1956
// Compute value(s).
1949
// Compute value(s)
1957
1950
for (unsigned int r = 0; r < 6; r++)
1958
1951
{
1959
1952
values[0] += coefficients0[r]*basisvalues[r];
1963
1956
case 6:
1964
1957
{
1965
1958
1966
// Array of basisvalues.
1959
// Array of basisvalues
1967
1960
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1968
1961
1969
// Declare helper variables.
1962
// Declare helper variables
1970
1963
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1971
1964
double tmp1 = (1.0 - Y)/2.0;
1972
1965
double tmp2 = tmp1*tmp1;
1973
1966
1974
// Compute basisvalues.
1967
// Compute basisvalues
1975
1968
basisvalues[0] = 1.0;
1976
1969
basisvalues[1] = tmp0;
1977
1970
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1985
1978
basisvalues[4] *= std::sqrt(4.5);
1986
1979
basisvalues[3] *= std::sqrt(7.5);
1987
1980
1988
// Table(s) of coefficients.
1981
// Table(s) of coefficients
1989
1982
static const double coefficients0[6] = \
1990
1983
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
1991
1984
1992
// Compute value(s).
1985
// Compute value(s)
1993
1986
for (unsigned int r = 0; r < 6; r++)
1994
1987
{
1995
1988
values[1] += coefficients0[r]*basisvalues[r];
1999
1992
case 7:
2000
1993
{
2001
1994
2002
// Array of basisvalues.
1995
// Array of basisvalues
2003
1996
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2004
1997
2005
// Declare helper variables.
1998
// Declare helper variables
2006
1999
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2007
2000
double tmp1 = (1.0 - Y)/2.0;
2008
2001
double tmp2 = tmp1*tmp1;
2009
2002
2010
// Compute basisvalues.
2003
// Compute basisvalues
2011
2004
basisvalues[0] = 1.0;
2012
2005
basisvalues[1] = tmp0;
2013
2006
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2021
2014
basisvalues[4] *= std::sqrt(4.5);
2022
2015
basisvalues[3] *= std::sqrt(7.5);
2023
2016
2024
// Table(s) of coefficients.
2017
// Table(s) of coefficients
2025
2018
static const double coefficients0[6] = \
2026
2019
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
2027
2020
2028
// Compute value(s).
2021
// Compute value(s)
2029
2022
for (unsigned int r = 0; r < 6; r++)
2030
2023
{
2031
2024
values[1] += coefficients0[r]*basisvalues[r];
2035
2028
case 8:
2036
2029
{
2037
2030
2038
// Array of basisvalues.
2031
// Array of basisvalues
2039
2032
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2040
2033
2041
// Declare helper variables.
2034
// Declare helper variables
2042
2035
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2043
2036
double tmp1 = (1.0 - Y)/2.0;
2044
2037
double tmp2 = tmp1*tmp1;
2045
2038
2046
// Compute basisvalues.
2039
// Compute basisvalues
2047
2040
basisvalues[0] = 1.0;
2048
2041
basisvalues[1] = tmp0;
2049
2042
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2057
2050
basisvalues[4] *= std::sqrt(4.5);
2058
2051
basisvalues[3] *= std::sqrt(7.5);
2059
2052
2060
// Table(s) of coefficients.
2053
// Table(s) of coefficients
2061
2054
static const double coefficients0[6] = \
2062
2055
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
2063
2056
2064
// Compute value(s).
2057
// Compute value(s)
2065
2058
for (unsigned int r = 0; r < 6; r++)
2066
2059
{
2067
2060
values[1] += coefficients0[r]*basisvalues[r];
2071
2064
case 9:
2072
2065
{
2073
2066
2074
// Array of basisvalues.
2067
// Array of basisvalues
2075
2068
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2076
2069
2077
// Declare helper variables.
2070
// Declare helper variables
2078
2071
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2079
2072
double tmp1 = (1.0 - Y)/2.0;
2080
2073
double tmp2 = tmp1*tmp1;
2081
2074
2082
// Compute basisvalues.
2075
// Compute basisvalues
2083
2076
basisvalues[0] = 1.0;
2084
2077
basisvalues[1] = tmp0;
2085
2078
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2093
2086
basisvalues[4] *= std::sqrt(4.5);
2094
2087
basisvalues[3] *= std::sqrt(7.5);
2095
2088
2096
// Table(s) of coefficients.
2089
// Table(s) of coefficients
2097
2090
static const double coefficients0[6] = \
2098
2091
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
2099
2092
2100
// Compute value(s).
2093
// Compute value(s)
2101
2094
for (unsigned int r = 0; r < 6; r++)
2102
2095
{
2103
2096
values[1] += coefficients0[r]*basisvalues[r];
2107
2100
case 10:
2108
2101
{
2109
2102
2110
// Array of basisvalues.
2103
// Array of basisvalues
2111
2104
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2112
2105
2113
// Declare helper variables.
2106
// Declare helper variables
2114
2107
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2115
2108
double tmp1 = (1.0 - Y)/2.0;
2116
2109
double tmp2 = tmp1*tmp1;
2117
2110
2118
// Compute basisvalues.
2111
// Compute basisvalues
2119
2112
basisvalues[0] = 1.0;
2120
2113
basisvalues[1] = tmp0;
2121
2114
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2129
2122
basisvalues[4] *= std::sqrt(4.5);
2130
2123
basisvalues[3] *= std::sqrt(7.5);
2131
2124
2132
// Table(s) of coefficients.
2125
// Table(s) of coefficients
2133
2126
static const double coefficients0[6] = \
2134
2127
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
2135
2128
2136
// Compute value(s).
2129
// Compute value(s)
2137
2130
for (unsigned int r = 0; r < 6; r++)
2138
2131
{
2139
2132
values[1] += coefficients0[r]*basisvalues[r];
2143
2136
case 11:
2144
2137
{
2145
2138
2146
// Array of basisvalues.
2139
// Array of basisvalues
2147
2140
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2148
2141
2149
// Declare helper variables.
2142
// Declare helper variables
2150
2143
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2151
2144
double tmp1 = (1.0 - Y)/2.0;
2152
2145
double tmp2 = tmp1*tmp1;
2153
2146
2154
// Compute basisvalues.
2147
// Compute basisvalues
2155
2148
basisvalues[0] = 1.0;
2156
2149
basisvalues[1] = tmp0;
2157
2150
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2165
2158
basisvalues[4] *= std::sqrt(4.5);
2166
2159
basisvalues[3] *= std::sqrt(7.5);
2167
2160
2168
// Table(s) of coefficients.
2161
// Table(s) of coefficients
2169
2162
static const double coefficients0[6] = \
2170
2163
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
2171
2164
2172
// Compute value(s).
2165
// Compute value(s)
2173
2166
for (unsigned int r = 0; r < 6; r++)
2174
2167
{
2175
2168
values[1] += coefficients0[r]*basisvalues[r];
2180
2173
2181
2174
}
2182
2175
2183
/// Evaluate all basis functions at given point in cell
2176
/// Evaluate all basis functions at given point x in cell
2184
2177
virtual void evaluate_basis_all(double* values,
2185
const double* coordinates,
2186
const ufc::cell& c) const
2178
const double* x,
2179
const double* vertex_coordinates,
2180
int cell_orientation) const
2187
2181
{
2188
2182
// Helper variable to hold values of a single dof.
2189
2183
double dof_values[2] = {0.0, 0.0};
2190
2184
2191
// Loop dofs and call evaluate_basis.
2185
// Loop dofs and call evaluate_basis
2192
2186
for (unsigned int r = 0; r < 12; r++)
2193
2187
{
2194
evaluate_basis(r, dof_values, coordinates, c);
2188
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
2195
2189
for (unsigned int s = 0; s < 2; s++)
2196
2190
{
2197
2191
values[r*2 + s] = dof_values[s];
2199
2193
}// end loop over 'r'
2200
2194
}
2201
2195
2202
/// Evaluate order n derivatives of basis function i at given point in cell
2203
virtual void evaluate_basis_derivatives(unsigned int i,
2204
unsigned int n,
2196
/// Evaluate order n derivatives of basis function i at given point x in cell
2197
virtual void evaluate_basis_derivatives(std::size_t i,
2198
std::size_t n,
2205
2199
double* values,
2206
const double* coordinates,
2207
const ufc::cell& c) const
2200
const double* x,
2201
const double* vertex_coordinates,
2202
int cell_orientation) const
2208
2203
{
2209
// Extract vertex coordinates
2210
const double * const * x = c.coordinates;
2211
2212
// Compute Jacobian of affine map from reference cell
2213
const double J_00 = x[1][0] - x[0][0];
2214
const double J_01 = x[2][0] - x[0][0];
2215
const double J_10 = x[1][1] - x[0][1];
2216
const double J_11 = x[2][1] - x[0][1];
2217
2218
// Compute determinant of Jacobian
2219
double detJ = J_00*J_11 - J_01*J_10;
2220
2221
// Compute inverse of Jacobian
2222
const double K_00 = J_11 / detJ;
2223
const double K_01 = -J_01 / detJ;
2224
const double K_10 = -J_10 / detJ;
2225
const double K_11 = J_00 / detJ;
2204
// Compute Jacobian
2205
double J[4];
2206
compute_jacobian_triangle_2d(J, vertex_coordinates);
2207
2208
// Compute Jacobian inverse and determinant
2209
double K[4];
2210
double detJ;
2211
compute_jacobian_inverse_triangle_2d(K, detJ, J);
2212
2226
2213
2227
2214
// Compute constants
2228
const double C0 = x[1][0] + x[2][0];
2229
const double C1 = x[1][1] + x[2][1];
2215
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
2216
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
2230
2217
2231
2218
// Get coordinates and map to the reference (FIAT) element
2232
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
2233
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
2219
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
2220
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
2234
2221
2235
2222
// Compute number of derivatives.
2236
2223
unsigned int num_derivatives = 1;
2267
2254
}
2268
2255
2269
2256
// Compute inverse of Jacobian
2270
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
2257
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
2271
2258
2272
2259
// Declare transformation matrix
2273
2260
// Declare pointer to two dimensional array and initialise
2301
2288
case 0:
2302
2289
{
2303
2290
2304
// Array of basisvalues.
2291
// Array of basisvalues
2305
2292
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2306
2293
2307
// Declare helper variables.
2294
// Declare helper variables
2308
2295
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2309
2296
double tmp1 = (1.0 - Y)/2.0;
2310
2297
double tmp2 = tmp1*tmp1;
2311
2298
2312
// Compute basisvalues.
2299
// Compute basisvalues
2313
2300
basisvalues[0] = 1.0;
2314
2301
basisvalues[1] = tmp0;
2315
2302
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2323
2310
basisvalues[4] *= std::sqrt(4.5);
2324
2311
basisvalues[3] *= std::sqrt(7.5);
2325
2312
2326
// Table(s) of coefficients.
2313
// Table(s) of coefficients
2327
2314
static const double coefficients0[6] = \
2328
2315
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
2329
2316
2467
2454
case 1:
2468
2455
{
2469
2456
2470
// Array of basisvalues.
2457
// Array of basisvalues
2471
2458
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2472
2459
2473
// Declare helper variables.
2460
// Declare helper variables
2474
2461
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2475
2462
double tmp1 = (1.0 - Y)/2.0;
2476
2463
double tmp2 = tmp1*tmp1;
2477
2464
2478
// Compute basisvalues.
2465
// Compute basisvalues
2479
2466
basisvalues[0] = 1.0;
2480
2467
basisvalues[1] = tmp0;
2481
2468
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2489
2476
basisvalues[4] *= std::sqrt(4.5);
2490
2477
basisvalues[3] *= std::sqrt(7.5);
2491
2478
2492
// Table(s) of coefficients.
2479
// Table(s) of coefficients
2493
2480
static const double coefficients0[6] = \
2494
2481
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
2495
2482
2633
2620
case 2:
2634
2621
{
2635
2622
2636
// Array of basisvalues.
2623
// Array of basisvalues
2637
2624
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2638
2625
2639
// Declare helper variables.
2626
// Declare helper variables
2640
2627
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2641
2628
double tmp1 = (1.0 - Y)/2.0;
2642
2629
double tmp2 = tmp1*tmp1;
2643
2630
2644
// Compute basisvalues.
2631
// Compute basisvalues
2645
2632
basisvalues[0] = 1.0;
2646
2633
basisvalues[1] = tmp0;
2647
2634
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2655
2642
basisvalues[4] *= std::sqrt(4.5);
2656
2643
basisvalues[3] *= std::sqrt(7.5);
2657
2644
2658
// Table(s) of coefficients.
2645
// Table(s) of coefficients
2659
2646
static const double coefficients0[6] = \
2660
2647
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
2661
2648
2799
2786
case 3:
2800
2787
{
2801
2788
2802
// Array of basisvalues.
2789
// Array of basisvalues
2803
2790
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2804
2791
2805
// Declare helper variables.
2792
// Declare helper variables
2806
2793
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2807
2794
double tmp1 = (1.0 - Y)/2.0;
2808
2795
double tmp2 = tmp1*tmp1;
2809
2796
2810
// Compute basisvalues.
2797
// Compute basisvalues
2811
2798
basisvalues[0] = 1.0;
2812
2799
basisvalues[1] = tmp0;
2813
2800
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2821
2808
basisvalues[4] *= std::sqrt(4.5);
2822
2809
basisvalues[3] *= std::sqrt(7.5);
2823
2810
2824
// Table(s) of coefficients.
2811
// Table(s) of coefficients
2825
2812
static const double coefficients0[6] = \
2826
2813
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
2827
2814
2965
2952
case 4:
2966
2953
{
2967
2954
2968
// Array of basisvalues.
2955
// Array of basisvalues
2969
2956
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2970
2957
2971
// Declare helper variables.
2958
// Declare helper variables
2972
2959
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2973
2960
double tmp1 = (1.0 - Y)/2.0;
2974
2961
double tmp2 = tmp1*tmp1;
2975
2962
2976
// Compute basisvalues.
2963
// Compute basisvalues
2977
2964
basisvalues[0] = 1.0;
2978
2965
basisvalues[1] = tmp0;
2979
2966
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2987
2974
basisvalues[4] *= std::sqrt(4.5);
2988
2975
basisvalues[3] *= std::sqrt(7.5);
2989
2976
2990
// Table(s) of coefficients.
2977
// Table(s) of coefficients
2991
2978
static const double coefficients0[6] = \
2992
2979
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
2993
2980
3131
3118
case 5:
3132
3119
{
3133
3120
3134
// Array of basisvalues.
3121
// Array of basisvalues
3135
3122
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3136
3123
3137
// Declare helper variables.
3124
// Declare helper variables
3138
3125
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3139
3126
double tmp1 = (1.0 - Y)/2.0;
3140
3127
double tmp2 = tmp1*tmp1;
3141
3128
3142
// Compute basisvalues.
3129
// Compute basisvalues
3143
3130
basisvalues[0] = 1.0;
3144
3131
basisvalues[1] = tmp0;
3145
3132
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3153
3140
basisvalues[4] *= std::sqrt(4.5);
3154
3141
basisvalues[3] *= std::sqrt(7.5);
3155
3142
3156
// Table(s) of coefficients.
3143
// Table(s) of coefficients
3157
3144
static const double coefficients0[6] = \
3158
3145
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
3159
3146
3297
3284
case 6:
3298
3285
{
3299
3286
3300
// Array of basisvalues.
3287
// Array of basisvalues
3301
3288
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3302
3289
3303
// Declare helper variables.
3290
// Declare helper variables
3304
3291
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3305
3292
double tmp1 = (1.0 - Y)/2.0;
3306
3293
double tmp2 = tmp1*tmp1;
3307
3294
3308
// Compute basisvalues.
3295
// Compute basisvalues
3309
3296
basisvalues[0] = 1.0;
3310
3297
basisvalues[1] = tmp0;
3311
3298
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3319
3306
basisvalues[4] *= std::sqrt(4.5);
3320
3307
basisvalues[3] *= std::sqrt(7.5);
3321
3308
3322
// Table(s) of coefficients.
3309
// Table(s) of coefficients
3323
3310
static const double coefficients0[6] = \
3324
3311
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
3325
3312
3463
3450
case 7:
3464
3451
{
3465
3452
3466
// Array of basisvalues.
3453
// Array of basisvalues
3467
3454
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3468
3455
3469
// Declare helper variables.
3456
// Declare helper variables
3470
3457
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3471
3458
double tmp1 = (1.0 - Y)/2.0;
3472
3459
double tmp2 = tmp1*tmp1;
3473
3460
3474
// Compute basisvalues.
3461
// Compute basisvalues
3475
3462
basisvalues[0] = 1.0;
3476
3463
basisvalues[1] = tmp0;
3477
3464
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3485
3472
basisvalues[4] *= std::sqrt(4.5);
3486
3473
basisvalues[3] *= std::sqrt(7.5);
3487
3474
3488
// Table(s) of coefficients.
3475
// Table(s) of coefficients
3489
3476
static const double coefficients0[6] = \
3490
3477
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
3491
3478
3629
3616
case 8:
3630
3617
{
3631
3618
3632
// Array of basisvalues.
3619
// Array of basisvalues
3633
3620
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3634
3621
3635
// Declare helper variables.
3622
// Declare helper variables
3636
3623
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3637
3624
double tmp1 = (1.0 - Y)/2.0;
3638
3625
double tmp2 = tmp1*tmp1;
3639
3626
3640
// Compute basisvalues.
3627
// Compute basisvalues
3641
3628
basisvalues[0] = 1.0;
3642
3629
basisvalues[1] = tmp0;
3643
3630
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3651
3638
basisvalues[4] *= std::sqrt(4.5);
3652
3639
basisvalues[3] *= std::sqrt(7.5);
3653
3640
3654
// Table(s) of coefficients.
3641
// Table(s) of coefficients
3655
3642
static const double coefficients0[6] = \
3656
3643
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
3657
3644
3795
3782
case 9:
3796
3783
{
3797
3784
3798
// Array of basisvalues.
3785
// Array of basisvalues
3799
3786
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3800
3787
3801
// Declare helper variables.
3788
// Declare helper variables
3802
3789
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3803
3790
double tmp1 = (1.0 - Y)/2.0;
3804
3791
double tmp2 = tmp1*tmp1;
3805
3792
3806
// Compute basisvalues.
3793
// Compute basisvalues
3807
3794
basisvalues[0] = 1.0;
3808
3795
basisvalues[1] = tmp0;
3809
3796
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3817
3804
basisvalues[4] *= std::sqrt(4.5);
3818
3805
basisvalues[3] *= std::sqrt(7.5);
3819
3806
3820
// Table(s) of coefficients.
3807
// Table(s) of coefficients
3821
3808
static const double coefficients0[6] = \
3822
3809
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
3823
3810
3961
3948
case 10:
3962
3949
{
3963
3950
3964
// Array of basisvalues.
3951
// Array of basisvalues
3965
3952
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3966
3953
3967
// Declare helper variables.
3954
// Declare helper variables
3968
3955
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3969
3956
double tmp1 = (1.0 - Y)/2.0;
3970
3957
double tmp2 = tmp1*tmp1;
3971
3958
3972
// Compute basisvalues.
3959
// Compute basisvalues
3973
3960
basisvalues[0] = 1.0;
3974
3961
basisvalues[1] = tmp0;
3975
3962
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3983
3970
basisvalues[4] *= std::sqrt(4.5);
3984
3971
basisvalues[3] *= std::sqrt(7.5);
3985
3972
3986
// Table(s) of coefficients.
3973
// Table(s) of coefficients
3987
3974
static const double coefficients0[6] = \
3988
3975
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
3989
3976
4127
4114
case 11:
4128
4115
{
4129
4116
4130
// Array of basisvalues.
4117
// Array of basisvalues
4131
4118
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
4132
4119
4133
// Declare helper variables.
4120
// Declare helper variables
4134
4121
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4135
4122
double tmp1 = (1.0 - Y)/2.0;
4136
4123
double tmp2 = tmp1*tmp1;
4137
4124
4138
// Compute basisvalues.
4125
// Compute basisvalues
4139
4126
basisvalues[0] = 1.0;
4140
4127
basisvalues[1] = tmp0;
4141
4128
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
4149
4136
basisvalues[4] *= std::sqrt(4.5);
4150
4137
basisvalues[3] *= std::sqrt(7.5);
4151
4138
4152
// Table(s) of coefficients.
4139
// Table(s) of coefficients
4153
4140
static const double coefficients0[6] = \
4154
4141
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
4155
4142
4294
4281
4295
4282
}
4296
4283
4297
/// Evaluate order n derivatives of all basis functions at given point in cell
4298
virtual void evaluate_basis_derivatives_all(unsigned int n,
4284
/// Evaluate order n derivatives of all basis functions at given point x in cell
4285
virtual void evaluate_basis_derivatives_all(std::size_t n,
4299
4286
double* values,
4300
const double* coordinates,
4301
const ufc::cell& c) const
4287
const double* x,
4288
const double* vertex_coordinates,
4289
int cell_orientation) const
4302
4290
{
4303
4291
// Compute number of derivatives.
4304
4292
unsigned int num_derivatives = 1;
4317
4305
// Loop dofs and call evaluate_basis_derivatives.
4318
4306
for (unsigned int r = 0; r < 12; r++)
4319
4307
{
4320
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
4308
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
4321
4309
for (unsigned int s = 0; s < 2*num_derivatives; s++)
4322
4310
{
4323
4311
values[r*2*num_derivatives + s] = dof_values[s];
4329
4317
}
4330
4318
4331
4319
/// Evaluate linear functional for dof i on the function f
4332
virtual double evaluate_dof(unsigned int i,
4320
virtual double evaluate_dof(std::size_t i,
4333
4321
const ufc::function& f,
4322
const double* vertex_coordinates,
4323
int cell_orientation,
4334
4324
const ufc::cell& c) const
4335
4325
{
4336
// Declare variables for result of evaluation.
4326
// Declare variables for result of evaluation
4337
4327
double vals[2];
4338
4328
4339
// Declare variable for physical coordinates.
4329
// Declare variable for physical coordinates
4340
4330
double y[2];
4341
const double * const * x = c.coordinates;
4342
4331
switch (i)
4343
4332
{
4344
4333
case 0:
4345
4334
{
4346
y[0] = x[0][0];
4347
y[1] = x[0][1];
4335
y[0] = vertex_coordinates[0];
4336
y[1] = vertex_coordinates[1];
4348
4337
f.evaluate(vals, y, c);
4349
4338
return vals[0];
4350
4339
break;
4351
4340
}
4352
4341
case 1:
4353
4342
{
4354
y[0] = x[1][0];
4355
y[1] = x[1][1];
4343
y[0] = vertex_coordinates[2];
4344
y[1] = vertex_coordinates[3];
4356
4345
f.evaluate(vals, y, c);
4357
4346
return vals[0];
4358
4347
break;
4359
4348
}
4360
4349
case 2:
4361
4350
{
4362
y[0] = x[2][0];
4363
y[1] = x[2][1];
4351
y[0] = vertex_coordinates[4];
4352
y[1] = vertex_coordinates[5];
4364
4353
f.evaluate(vals, y, c);
4365
4354
return vals[0];
4366
4355
break;
4367
4356
}
4368
4357
case 3:
4369
4358
{
4370
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
4371
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
4359
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
4360
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
4372
4361
f.evaluate(vals, y, c);
4373
4362
return vals[0];
4374
4363
break;
4375
4364
}
4376
4365
case 4:
4377
4366
{
4378
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
4379
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
4367
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
4368
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
4380
4369
f.evaluate(vals, y, c);
4381
4370
return vals[0];
4382
4371
break;
4383
4372
}
4384
4373
case 5:
4385
4374
{
4386
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
4387
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
4375
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
4376
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
4388
4377
f.evaluate(vals, y, c);
4389
4378
return vals[0];
4390
4379
break;
4391
4380
}
4392
4381
case 6:
4393
4382
{
4394
y[0] = x[0][0];
4395
y[1] = x[0][1];
4383
y[0] = vertex_coordinates[0];
4384
y[1] = vertex_coordinates[1];
4396
4385
f.evaluate(vals, y, c);
4397
4386
return vals[1];
4398
4387
break;
4399
4388
}
4400
4389
case 7:
4401
4390
{
4402
y[0] = x[1][0];
4403
y[1] = x[1][1];
4391
y[0] = vertex_coordinates[2];
4392
y[1] = vertex_coordinates[3];
4404
4393
f.evaluate(vals, y, c);
4405
4394
return vals[1];
4406
4395
break;
4407
4396
}
4408
4397
case 8:
4409
4398
{
4410
y[0] = x[2][0];
4411
y[1] = x[2][1];
4399
y[0] = vertex_coordinates[4];
4400
y[1] = vertex_coordinates[5];
4412
4401
f.evaluate(vals, y, c);
4413
4402
return vals[1];
4414
4403
break;
4415
4404
}
4416
4405
case 9:
4417
4406
{
4418
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
4419
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
4407
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
4408
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
4420
4409
f.evaluate(vals, y, c);
4421
4410
return vals[1];
4422
4411
break;
4423
4412
}
4424
4413
case 10:
4425
4414
{
4426
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
4427
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
4415
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
4416
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
4428
4417
f.evaluate(vals, y, c);
4429
4418
return vals[1];
4430
4419
break;
4431
4420
}
4432
4421
case 11:
4433
4422
{
4434
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
4435
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
4423
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
4424
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
4436
4425
f.evaluate(vals, y, c);
4437
4426
return vals[1];
4438
4427
break;
4445
4434
/// Evaluate linear functionals for all dofs on the function f
4446
4435
virtual void evaluate_dofs(double* values,
4447
4436
const ufc::function& f,
4437
const double* vertex_coordinates,
4438
int cell_orientation,
4448
4439
const ufc::cell& c) const
4449
4440
{
4450
// Declare variables for result of evaluation.
4441
// Declare variables for result of evaluation
4451
4442
double vals[2];
4452
4443
4453
// Declare variable for physical coordinates.
4444
// Declare variable for physical coordinates
4454
4445
double y[2];
4455
const double * const * x = c.coordinates;
4456
y[0] = x[0][0];
4457
y[1] = x[0][1];
4446
y[0] = vertex_coordinates[0];
4447
y[1] = vertex_coordinates[1];
4458
4448
f.evaluate(vals, y, c);
4459
4449
values[0] = vals[0];
4460
y[0] = x[1][0];
4461
y[1] = x[1][1];
4450
y[0] = vertex_coordinates[2];
4451
y[1] = vertex_coordinates[3];
4462
4452
f.evaluate(vals, y, c);
4463
4453
values[1] = vals[0];
4464
y[0] = x[2][0];
4465
y[1] = x[2][1];
4454
y[0] = vertex_coordinates[4];
4455
y[1] = vertex_coordinates[5];
4466
4456
f.evaluate(vals, y, c);
4467
4457
values[2] = vals[0];
4468
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
4469
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
4458
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
4459
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
4470
4460
f.evaluate(vals, y, c);
4471
4461
values[3] = vals[0];
4472
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
4473
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
4462
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
4463
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
4474
4464
f.evaluate(vals, y, c);
4475
4465
values[4] = vals[0];
4476
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
4477
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
4466
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
4467
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
4478
4468
f.evaluate(vals, y, c);
4479
4469
values[5] = vals[0];
4480
y[0] = x[0][0];
4481
y[1] = x[0][1];
4470
y[0] = vertex_coordinates[0];
4471
y[1] = vertex_coordinates[1];
4482
4472
f.evaluate(vals, y, c);
4483
4473
values[6] = vals[1];
4484
y[0] = x[1][0];
4485
y[1] = x[1][1];
4474
y[0] = vertex_coordinates[2];
4475
y[1] = vertex_coordinates[3];
4486
4476
f.evaluate(vals, y, c);
4487
4477
values[7] = vals[1];
4488
y[0] = x[2][0];
4489
y[1] = x[2][1];
4478
y[0] = vertex_coordinates[4];
4479
y[1] = vertex_coordinates[5];
4490
4480
f.evaluate(vals, y, c);
4491
4481
values[8] = vals[1];
4492
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
4493
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
4482
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
4483
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
4494
4484
f.evaluate(vals, y, c);
4495
4485
values[9] = vals[1];
4496
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
4497
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
4486
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
4487
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
4498
4488
f.evaluate(vals, y, c);
4499
4489
values[10] = vals[1];
4500
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
4501
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
4490
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
4491
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
4502
4492
f.evaluate(vals, y, c);
4503
4493
values[11] = vals[1];
4504
4494
}
4506
4496
/// Interpolate vertex values from dof values
4507
4497
virtual void interpolate_vertex_values(double* vertex_values,
4508
4498
const double* dof_values,
4499
const double* vertex_coordinates,
4500
int cell_orientation,
4509
4501
const ufc::cell& c) const
4510
4502
{
4511
4503
// Evaluate function and change variables
4523
4515
const double* xhat,
4524
4516
const ufc::cell& c) const
4525
4517
{
4526
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4518
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
4527
4519
}
4528
4520
4529
4521
/// Map from coordinate x in cell to coordinate xhat in reference cell
4531
4523
const double* x,
4532
4524
const ufc::cell& c) const
4533
4525
{
4534
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4526
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
4535
4527
}
4536
4528
4537
4529
/// Return the number of sub elements (for a mixed element)
4538
virtual unsigned int num_sub_elements() const
4530
virtual std::size_t num_sub_elements() const
4539
4531
{
4540
4532
return 2;
4541
4533
}
4542
4534
4543
4535
/// Create a new finite element for sub element i (for a mixed element)
4544
virtual ufc::finite_element* create_sub_element(unsigned int i) const
4536
virtual ufc::finite_element* create_sub_element(std::size_t i) const
4545
4537
{
4546
4538
switch (i)
4547
4539
{
4589
4581
/// Return a string identifying the finite element
4590
4582
virtual const char* signature() const
4591
4583
{
4592
return "FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)";
4584
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
4593
4585
}
4594
4586
4595
4587
/// Return the cell shape
4599
4591
}
4600
4592
4601
4593
/// Return the topological dimension of the cell shape
4602
virtual unsigned int topological_dimension() const
4594
virtual std::size_t topological_dimension() const
4603
4595
{
4604
4596
return 2;
4605
4597
}
4606
4598
4607
4599
/// Return the geometric dimension of the cell shape
4608
virtual unsigned int geometric_dimension() const
4600
virtual std::size_t geometric_dimension() const
4609
4601
{
4610
4602
return 2;
4611
4603
}
4612
4604
4613
4605
/// Return the dimension of the finite element function space
4614
virtual unsigned int space_dimension() const
4606
virtual std::size_t space_dimension() const
4615
4607
{
4616
4608
return 3;
4617
4609
}
4618
4610
4619
4611
/// Return the rank of the value space
4620
virtual unsigned int value_rank() const
4612
virtual std::size_t value_rank() const
4621
4613
{
4622
4614
return 0;
4623
4615
}
4624
4616
4625
4617
/// Return the dimension of the value space for axis i
4626
virtual unsigned int value_dimension(unsigned int i) const
4618
virtual std::size_t value_dimension(std::size_t i) const
4627
4619
{
4628
4620
return 1;
4629
4621
}
4630
4622
4631
/// Evaluate basis function i at given point in cell
4632
virtual void evaluate_basis(unsigned int i,
4623
/// Evaluate basis function i at given point x in cell
4624
virtual void evaluate_basis(std::size_t i,
4633
4625
double* values,
4634
const double* coordinates,
4635
const ufc::cell& c) const
4626
const double* x,
4627
const double* vertex_coordinates,
4628
int cell_orientation) const
4636
4629
{
4637
// Extract vertex coordinates
4638
const double * const * x = c.coordinates;
4639
4640
// Compute Jacobian of affine map from reference cell
4641
const double J_00 = x[1][0] - x[0][0];
4642
const double J_01 = x[2][0] - x[0][0];
4643
const double J_10 = x[1][1] - x[0][1];
4644
const double J_11 = x[2][1] - x[0][1];
4645
4646
// Compute determinant of Jacobian
4647
double detJ = J_00*J_11 - J_01*J_10;
4648
4649
// Compute inverse of Jacobian
4630
// Compute Jacobian
4631
double J[4];
4632
compute_jacobian_triangle_2d(J, vertex_coordinates);
4633
4634
// Compute Jacobian inverse and determinant
4635
double K[4];
4636
double detJ;
4637
compute_jacobian_inverse_triangle_2d(K, detJ, J);
4638
4650
4639
4651
4640
// Compute constants
4652
const double C0 = x[1][0] + x[2][0];
4653
const double C1 = x[1][1] + x[2][1];
4641
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
4642
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
4654
4643
4655
4644
// Get coordinates and map to the reference (FIAT) element
4656
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
4657
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
4645
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
4646
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
4658
4647
4659
// Reset values.
4648
// Reset values
4660
4649
*values = 0.0;
4661
4650
switch (i)
4662
4651
{
4663
4652
case 0:
4664
4653
{
4665
4654
4666
// Array of basisvalues.
4655
// Array of basisvalues
4667
4656
double basisvalues[3] = {0.0, 0.0, 0.0};
4668
4657
4669
// Declare helper variables.
4658
// Declare helper variables
4670
4659
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4671
4660
4672
// Compute basisvalues.
4661
// Compute basisvalues
4673
4662
basisvalues[0] = 1.0;
4674
4663
basisvalues[1] = tmp0;
4675
4664
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
4677
4666
basisvalues[2] *= std::sqrt(1.0);
4678
4667
basisvalues[1] *= std::sqrt(3.0);
4679
4668
4680
// Table(s) of coefficients.
4669
// Table(s) of coefficients
4681
4670
static const double coefficients0[3] = \
4682
4671
{0.47140452, -0.28867513, -0.16666667};
4683
4672
4684
// Compute value(s).
4673
// Compute value(s)
4685
4674
for (unsigned int r = 0; r < 3; r++)
4686
4675
{
4687
4676
*values += coefficients0[r]*basisvalues[r];
4691
4680
case 1:
4692
4681
{
4693
4682
4694
// Array of basisvalues.
4683
// Array of basisvalues
4695
4684
double basisvalues[3] = {0.0, 0.0, 0.0};
4696
4685
4697
// Declare helper variables.
4686
// Declare helper variables
4698
4687
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4699
4688
4700
// Compute basisvalues.
4689
// Compute basisvalues
4701
4690
basisvalues[0] = 1.0;
4702
4691
basisvalues[1] = tmp0;
4703
4692
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
4705
4694
basisvalues[2] *= std::sqrt(1.0);
4706
4695
basisvalues[1] *= std::sqrt(3.0);
4707
4696
4708
// Table(s) of coefficients.
4697
// Table(s) of coefficients
4709
4698
static const double coefficients0[3] = \
4710
4699
{0.47140452, 0.28867513, -0.16666667};
4711
4700
4712
// Compute value(s).
4701
// Compute value(s)
4713
4702
for (unsigned int r = 0; r < 3; r++)
4714
4703
{
4715
4704
*values += coefficients0[r]*basisvalues[r];
4719
4708
case 2:
4720
4709
{
4721
4710
4722
// Array of basisvalues.
4711
// Array of basisvalues
4723
4712
double basisvalues[3] = {0.0, 0.0, 0.0};
4724
4713
4725
// Declare helper variables.
4714
// Declare helper variables
4726
4715
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4727
4716
4728
// Compute basisvalues.
4717
// Compute basisvalues
4729
4718
basisvalues[0] = 1.0;
4730
4719
basisvalues[1] = tmp0;
4731
4720
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
4733
4722
basisvalues[2] *= std::sqrt(1.0);
4734
4723
basisvalues[1] *= std::sqrt(3.0);
4735
4724
4736
// Table(s) of coefficients.
4725
// Table(s) of coefficients
4737
4726
static const double coefficients0[3] = \
4738
4727
{0.47140452, 0.0, 0.33333333};
4739
4728
4740
// Compute value(s).
4729
// Compute value(s)
4741
4730
for (unsigned int r = 0; r < 3; r++)
4742
4731
{
4743
4732
*values += coefficients0[r]*basisvalues[r];
4748
4737
4749
4738
}
4750
4739
4751
/// Evaluate all basis functions at given point in cell
4740
/// Evaluate all basis functions at given point x in cell
4752
4741
virtual void evaluate_basis_all(double* values,
4753
const double* coordinates,
4754
const ufc::cell& c) const
4742
const double* x,
4743
const double* vertex_coordinates,
4744
int cell_orientation) const
4755
4745
{
4756
4746
// Helper variable to hold values of a single dof.
4757
4747
double dof_values = 0.0;
4758
4748
4759
// Loop dofs and call evaluate_basis.
4749
// Loop dofs and call evaluate_basis
4760
4750
for (unsigned int r = 0; r < 3; r++)
4761
4751
{
4762
evaluate_basis(r, &dof_values, coordinates, c);
4752
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
4763
4753
values[r] = dof_values;
4764
4754
}// end loop over 'r'
4765
4755
}
4766
4756
4767
/// Evaluate order n derivatives of basis function i at given point in cell
4768
virtual void evaluate_basis_derivatives(unsigned int i,
4769
unsigned int n,
4757
/// Evaluate order n derivatives of basis function i at given point x in cell
4758
virtual void evaluate_basis_derivatives(std::size_t i,
4759
std::size_t n,
4770
4760
double* values,
4771
const double* coordinates,
4772
const ufc::cell& c) const
4761
const double* x,
4762
const double* vertex_coordinates,
4763
int cell_orientation) const
4773
4764
{
4774
// Extract vertex coordinates
4775
const double * const * x = c.coordinates;
4776
4777
// Compute Jacobian of affine map from reference cell
4778
const double J_00 = x[1][0] - x[0][0];
4779
const double J_01 = x[2][0] - x[0][0];
4780
const double J_10 = x[1][1] - x[0][1];
4781
const double J_11 = x[2][1] - x[0][1];
4782
4783
// Compute determinant of Jacobian
4784
double detJ = J_00*J_11 - J_01*J_10;
4785
4786
// Compute inverse of Jacobian
4787
const double K_00 = J_11 / detJ;
4788
const double K_01 = -J_01 / detJ;
4789
const double K_10 = -J_10 / detJ;
4790
const double K_11 = J_00 / detJ;
4765
// Compute Jacobian
4766
double J[4];
4767
compute_jacobian_triangle_2d(J, vertex_coordinates);
4768
4769
// Compute Jacobian inverse and determinant
4770
double K[4];
4771
double detJ;
4772
compute_jacobian_inverse_triangle_2d(K, detJ, J);
4773
4791
4774
4792
4775
// Compute constants
4793
const double C0 = x[1][0] + x[2][0];
4794
const double C1 = x[1][1] + x[2][1];
4776
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
4777
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
4795
4778
4796
4779
// Get coordinates and map to the reference (FIAT) element
4797
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
4798
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
4780
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
4781
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
4799
4782
4800
4783
// Compute number of derivatives.
4801
4784
unsigned int num_derivatives = 1;
4832
4815
}
4833
4816
4834
4817
// Compute inverse of Jacobian
4835
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
4818
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
4836
4819
4837
4820
// Declare transformation matrix
4838
4821
// Declare pointer to two dimensional array and initialise
4866
4849
case 0:
4867
4850
{
4868
4851
4869
// Array of basisvalues.
4852
// Array of basisvalues
4870
4853
double basisvalues[3] = {0.0, 0.0, 0.0};
4871
4854
4872
// Declare helper variables.
4855
// Declare helper variables
4873
4856
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4874
4857
4875
// Compute basisvalues.
4858
// Compute basisvalues
4876
4859
basisvalues[0] = 1.0;
4877
4860
basisvalues[1] = tmp0;
4878
4861
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
4880
4863
basisvalues[2] *= std::sqrt(1.0);
4881
4864
basisvalues[1] *= std::sqrt(3.0);
4882
4865
4883
// Table(s) of coefficients.
4866
// Table(s) of coefficients
4884
4867
static const double coefficients0[3] = \
4885
4868
{0.47140452, -0.28867513, -0.16666667};
4886
4869
5012
4995
case 1:
5013
4996
{
5014
4997
5015
// Array of basisvalues.
4998
// Array of basisvalues
5016
4999
double basisvalues[3] = {0.0, 0.0, 0.0};
5017
5000
5018
// Declare helper variables.
5001
// Declare helper variables
5019
5002
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5020
5003
5021
// Compute basisvalues.
5004
// Compute basisvalues
5022
5005
basisvalues[0] = 1.0;
5023
5006
basisvalues[1] = tmp0;
5024
5007
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
5026
5009
basisvalues[2] *= std::sqrt(1.0);
5027
5010
basisvalues[1] *= std::sqrt(3.0);
5028
5011
5029
// Table(s) of coefficients.
5012
// Table(s) of coefficients
5030
5013
static const double coefficients0[3] = \
5031
5014
{0.47140452, 0.28867513, -0.16666667};
5032
5015
5158
5141
case 2:
5159
5142
{
5160
5143
5161
// Array of basisvalues.
5144
// Array of basisvalues
5162
5145
double basisvalues[3] = {0.0, 0.0, 0.0};
5163
5146
5164
// Declare helper variables.
5147
// Declare helper variables
5165
5148
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5166
5149
5167
// Compute basisvalues.
5150
// Compute basisvalues
5168
5151
basisvalues[0] = 1.0;
5169
5152
basisvalues[1] = tmp0;
5170
5153
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
5172
5155
basisvalues[2] *= std::sqrt(1.0);
5173
5156
basisvalues[1] *= std::sqrt(3.0);
5174
5157
5175
// Table(s) of coefficients.
5158
// Table(s) of coefficients
5176
5159
static const double coefficients0[3] = \
5177
5160
{0.47140452, 0.0, 0.33333333};
5178
5161
5305
5288
5306
5289
}
5307
5290
5308
/// Evaluate order n derivatives of all basis functions at given point in cell
5309
virtual void evaluate_basis_derivatives_all(unsigned int n,
5291
/// Evaluate order n derivatives of all basis functions at given point x in cell
5292
virtual void evaluate_basis_derivatives_all(std::size_t n,
5310
5293
double* values,
5311
const double* coordinates,
5312
const ufc::cell& c) const
5294
const double* x,
5295
const double* vertex_coordinates,
5296
int cell_orientation) const
5313
5297
{
5314
5298
// Compute number of derivatives.
5315
5299
unsigned int num_derivatives = 1;
5328
5312
// Loop dofs and call evaluate_basis_derivatives.
5329
5313
for (unsigned int r = 0; r < 3; r++)
5330
5314
{
5331
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
5315
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
5332
5316
for (unsigned int s = 0; s < num_derivatives; s++)
5333
5317
{
5334
5318
values[r*num_derivatives + s] = dof_values[s];
5340
5324
}
5341
5325
5342
5326
/// Evaluate linear functional for dof i on the function f
5343
virtual double evaluate_dof(unsigned int i,
5327
virtual double evaluate_dof(std::size_t i,
5344
5328
const ufc::function& f,
5329
const double* vertex_coordinates,
5330
int cell_orientation,
5345
5331
const ufc::cell& c) const
5346
5332
{
5347
// Declare variables for result of evaluation.
5333
// Declare variables for result of evaluation
5348
5334
double vals[1];
5349
5335
5350
// Declare variable for physical coordinates.
5336
// Declare variable for physical coordinates
5351
5337
double y[2];
5352
const double * const * x = c.coordinates;
5353
5338
switch (i)
5354
5339
{
5355
5340
case 0:
5356
5341
{
5357
y[0] = x[0][0];
5358
y[1] = x[0][1];
5342
y[0] = vertex_coordinates[0];
5343
y[1] = vertex_coordinates[1];
5359
5344
f.evaluate(vals, y, c);
5360
5345
return vals[0];
5361
5346
break;
5362
5347
}
5363
5348
case 1:
5364
5349
{
5365
y[0] = x[1][0];
5366
y[1] = x[1][1];
5350
y[0] = vertex_coordinates[2];
5351
y[1] = vertex_coordinates[3];
5367
5352
f.evaluate(vals, y, c);
5368
5353
return vals[0];
5369
5354
break;
5370
5355
}
5371
5356
case 2:
5372
5357
{
5373
y[0] = x[2][0];
5374
y[1] = x[2][1];
5358
y[0] = vertex_coordinates[4];
5359
y[1] = vertex_coordinates[5];
5375
5360
f.evaluate(vals, y, c);
5376
5361
return vals[0];
5377
5362
break;
5384
5369
/// Evaluate linear functionals for all dofs on the function f
5385
5370
virtual void evaluate_dofs(double* values,
5386
5371
const ufc::function& f,
5372
const double* vertex_coordinates,
5373
int cell_orientation,
5387
5374
const ufc::cell& c) const
5388
5375
{
5389
// Declare variables for result of evaluation.
5376
// Declare variables for result of evaluation
5390
5377
double vals[1];
5391
5378
5392
// Declare variable for physical coordinates.
5379
// Declare variable for physical coordinates
5393
5380
double y[2];
5394
const double * const * x = c.coordinates;
5395
y[0] = x[0][0];
5396
y[1] = x[0][1];
5381
y[0] = vertex_coordinates[0];
5382
y[1] = vertex_coordinates[1];
5397
5383
f.evaluate(vals, y, c);
5398
5384
values[0] = vals[0];
5399
y[0] = x[1][0];
5400
y[1] = x[1][1];
5385
y[0] = vertex_coordinates[2];
5386
y[1] = vertex_coordinates[3];
5401
5387
f.evaluate(vals, y, c);
5402
5388
values[1] = vals[0];
5403
y[0] = x[2][0];
5404
y[1] = x[2][1];
5389
y[0] = vertex_coordinates[4];
5390
y[1] = vertex_coordinates[5];
5405
5391
f.evaluate(vals, y, c);
5406
5392
values[2] = vals[0];
5407
5393
}
5409
5395
/// Interpolate vertex values from dof values
5410
5396
virtual void interpolate_vertex_values(double* vertex_values,
5411
5397
const double* dof_values,
5398
const double* vertex_coordinates,
5399
int cell_orientation,
5412
5400
const ufc::cell& c) const
5413
5401
{
5414
5402
// Evaluate function and change variables
5422
5410
const double* xhat,
5423
5411
const ufc::cell& c) const
5424
5412
{
5425
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
5413
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
5426
5414
}
5427
5415
5428
5416
/// Map from coordinate x in cell to coordinate xhat in reference cell
5430
5418
const double* x,
5431
5419
const ufc::cell& c) const
5432
5420
{
5433
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
5421
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
5434
5422
}
5435
5423
5436
5424
/// Return the number of sub elements (for a mixed element)
5437
virtual unsigned int num_sub_elements() const
5425
virtual std::size_t num_sub_elements() const
5438
5426
{
5439
5427
return 0;
5440
5428
}
5441
5429
5442
5430
/// Create a new finite element for sub element i (for a mixed element)
5443
virtual ufc::finite_element* create_sub_element(unsigned int i) const
5431
virtual ufc::finite_element* create_sub_element(std::size_t i) const
5444
5432
{
5445
5433
return 0;
5446
5434
}
5474
5462
/// Return a string identifying the finite element
5475
5463
virtual const char* signature() const
5476
5464
{
5477
return "MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) })";
5465
return "MixedElement(*[VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (3,) })";
5478
5466
}
5479
5467
5480
5468
/// Return the cell shape
5484
5472
}
5485
5473
5486
5474
/// Return the topological dimension of the cell shape
5487
virtual unsigned int topological_dimension() const
5475
virtual std::size_t topological_dimension() const
5488
5476
{
5489
5477
return 2;
5490
5478
}
5491
5479
5492
5480
/// Return the geometric dimension of the cell shape
5493
virtual unsigned int geometric_dimension() const
5481
virtual std::size_t geometric_dimension() const
5494
5482
{
5495
5483
return 2;
5496
5484
}
5497
5485
5498
5486
/// Return the dimension of the finite element function space
5499
virtual unsigned int space_dimension() const
5487
virtual std::size_t space_dimension() const
5500
5488
{
5501
5489
return 15;
5502
5490
}
5503
5491
5504
5492
/// Return the rank of the value space
5505
virtual unsigned int value_rank() const
5493
virtual std::size_t value_rank() const
5506
5494
{
5507
5495
return 1;
5508
5496
}
5509
5497
5510
5498
/// Return the dimension of the value space for axis i
5511
virtual unsigned int value_dimension(unsigned int i) const
5499
virtual std::size_t value_dimension(std::size_t i) const
5512
5500
{
5513
5501
switch (i)
5514
5502
{
5522
5510
return 0;
5523
5511
}
5524
5512
5525
/// Evaluate basis function i at given point in cell
5526
virtual void evaluate_basis(unsigned int i,
5513
/// Evaluate basis function i at given point x in cell
5514
virtual void evaluate_basis(std::size_t i,
5527
5515
double* values,
5528
const double* coordinates,
5529
const ufc::cell& c) const
5516
const double* x,
5517
const double* vertex_coordinates,
5518
int cell_orientation) const
5530
5519
{
5531
// Extract vertex coordinates
5532
const double * const * x = c.coordinates;
5533
5534
// Compute Jacobian of affine map from reference cell
5535
const double J_00 = x[1][0] - x[0][0];
5536
const double J_01 = x[2][0] - x[0][0];
5537
const double J_10 = x[1][1] - x[0][1];
5538
const double J_11 = x[2][1] - x[0][1];
5539
5540
// Compute determinant of Jacobian
5541
double detJ = J_00*J_11 - J_01*J_10;
5542
5543
// Compute inverse of Jacobian
5520
// Compute Jacobian
5521
double J[4];
5522
compute_jacobian_triangle_2d(J, vertex_coordinates);
5523
5524
// Compute Jacobian inverse and determinant
5525
double K[4];
5526
double detJ;
5527
compute_jacobian_inverse_triangle_2d(K, detJ, J);
5528
5544
5529
5545
5530
// Compute constants
5546
const double C0 = x[1][0] + x[2][0];
5547
const double C1 = x[1][1] + x[2][1];
5531
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
5532
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
5548
5533
5549
5534
// Get coordinates and map to the reference (FIAT) element
5550
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
5551
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
5535
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
5536
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
5552
5537
5553
// Reset values.
5538
// Reset values
5554
5539
values[0] = 0.0;
5555
5540
values[1] = 0.0;
5556
5541
values[2] = 0.0;
5559
5544
case 0:
5560
5545
{
5561
5546
5562
// Array of basisvalues.
5547
// Array of basisvalues
5563
5548
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5564
5549
5565
// Declare helper variables.
5550
// Declare helper variables
5566
5551
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5567
5552
double tmp1 = (1.0 - Y)/2.0;
5568
5553
double tmp2 = tmp1*tmp1;
5569
5554
5570
// Compute basisvalues.
5555
// Compute basisvalues
5571
5556
basisvalues[0] = 1.0;
5572
5557
basisvalues[1] = tmp0;
5573
5558
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5581
5566
basisvalues[4] *= std::sqrt(4.5);
5582
5567
basisvalues[3] *= std::sqrt(7.5);
5583
5568
5584
// Table(s) of coefficients.
5569
// Table(s) of coefficients
5585
5570
static const double coefficients0[6] = \
5586
5571
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
5587
5572
5588
// Compute value(s).
5573
// Compute value(s)
5589
5574
for (unsigned int r = 0; r < 6; r++)
5590
5575
{
5591
5576
values[0] += coefficients0[r]*basisvalues[r];
5595
5580
case 1:
5596
5581
{
5597
5582
5598
// Array of basisvalues.
5583
// Array of basisvalues
5599
5584
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5600
5585
5601
// Declare helper variables.
5586
// Declare helper variables
5602
5587
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5603
5588
double tmp1 = (1.0 - Y)/2.0;
5604
5589
double tmp2 = tmp1*tmp1;
5605
5590
5606
// Compute basisvalues.
5591
// Compute basisvalues
5607
5592
basisvalues[0] = 1.0;
5608
5593
basisvalues[1] = tmp0;
5609
5594
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5617
5602
basisvalues[4] *= std::sqrt(4.5);
5618
5603
basisvalues[3] *= std::sqrt(7.5);
5619
5604
5620
// Table(s) of coefficients.
5605
// Table(s) of coefficients
5621
5606
static const double coefficients0[6] = \
5622
5607
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
5623
5608
5624
// Compute value(s).
5609
// Compute value(s)
5625
5610
for (unsigned int r = 0; r < 6; r++)
5626
5611
{
5627
5612
values[0] += coefficients0[r]*basisvalues[r];
5631
5616
case 2:
5632
5617
{
5633
5618
5634
// Array of basisvalues.
5619
// Array of basisvalues
5635
5620
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5636
5621
5637
// Declare helper variables.
5622
// Declare helper variables
5638
5623
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5639
5624
double tmp1 = (1.0 - Y)/2.0;
5640
5625
double tmp2 = tmp1*tmp1;
5641
5626
5642
// Compute basisvalues.
5627
// Compute basisvalues
5643
5628
basisvalues[0] = 1.0;
5644
5629
basisvalues[1] = tmp0;
5645
5630
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5653
5638
basisvalues[4] *= std::sqrt(4.5);
5654
5639
basisvalues[3] *= std::sqrt(7.5);
5655
5640
5656
// Table(s) of coefficients.
5641
// Table(s) of coefficients
5657
5642
static const double coefficients0[6] = \
5658
5643
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
5659
5644
5660
// Compute value(s).
5645
// Compute value(s)
5661
5646
for (unsigned int r = 0; r < 6; r++)
5662
5647
{
5663
5648
values[0] += coefficients0[r]*basisvalues[r];
5667
5652
case 3:
5668
5653
{
5669
5654
5670
// Array of basisvalues.
5655
// Array of basisvalues
5671
5656
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5672
5657
5673
// Declare helper variables.
5658
// Declare helper variables
5674
5659
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5675
5660
double tmp1 = (1.0 - Y)/2.0;
5676
5661
double tmp2 = tmp1*tmp1;
5677
5662
5678
// Compute basisvalues.
5663
// Compute basisvalues
5679
5664
basisvalues[0] = 1.0;
5680
5665
basisvalues[1] = tmp0;
5681
5666
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5689
5674
basisvalues[4] *= std::sqrt(4.5);
5690
5675
basisvalues[3] *= std::sqrt(7.5);
5691
5676
5692
// Table(s) of coefficients.
5677
// Table(s) of coefficients
5693
5678
static const double coefficients0[6] = \
5694
5679
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
5695
5680
5696
// Compute value(s).
5681
// Compute value(s)
5697
5682
for (unsigned int r = 0; r < 6; r++)
5698
5683
{
5699
5684
values[0] += coefficients0[r]*basisvalues[r];
5703
5688
case 4:
5704
5689
{
5705
5690
5706
// Array of basisvalues.
5691
// Array of basisvalues
5707
5692
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5708
5693
5709
// Declare helper variables.
5694
// Declare helper variables
5710
5695
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5711
5696
double tmp1 = (1.0 - Y)/2.0;
5712
5697
double tmp2 = tmp1*tmp1;
5713
5698
5714
// Compute basisvalues.
5699
// Compute basisvalues
5715
5700
basisvalues[0] = 1.0;
5716
5701
basisvalues[1] = tmp0;
5717
5702
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5725
5710
basisvalues[4] *= std::sqrt(4.5);
5726
5711
basisvalues[3] *= std::sqrt(7.5);
5727
5712
5728
// Table(s) of coefficients.
5713
// Table(s) of coefficients
5729
5714
static const double coefficients0[6] = \
5730
5715
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
5731
5716
5732
// Compute value(s).
5717
// Compute value(s)
5733
5718
for (unsigned int r = 0; r < 6; r++)
5734
5719
{
5735
5720
values[0] += coefficients0[r]*basisvalues[r];
5739
5724
case 5:
5740
5725
{
5741
5726
5742
// Array of basisvalues.
5727
// Array of basisvalues
5743
5728
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5744
5729
5745
// Declare helper variables.
5730
// Declare helper variables
5746
5731
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5747
5732
double tmp1 = (1.0 - Y)/2.0;
5748
5733
double tmp2 = tmp1*tmp1;
5749
5734
5750
// Compute basisvalues.
5735
// Compute basisvalues
5751
5736
basisvalues[0] = 1.0;
5752
5737
basisvalues[1] = tmp0;
5753
5738
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5761
5746
basisvalues[4] *= std::sqrt(4.5);
5762
5747
basisvalues[3] *= std::sqrt(7.5);
5763
5748
5764
// Table(s) of coefficients.
5749
// Table(s) of coefficients
5765
5750
static const double coefficients0[6] = \
5766
5751
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
5767
5752
5768
// Compute value(s).
5753
// Compute value(s)
5769
5754
for (unsigned int r = 0; r < 6; r++)
5770
5755
{
5771
5756
values[0] += coefficients0[r]*basisvalues[r];
5775
5760
case 6:
5776
5761
{
5777
5762
5778
// Array of basisvalues.
5763
// Array of basisvalues
5779
5764
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5780
5765
5781
// Declare helper variables.
5766
// Declare helper variables
5782
5767
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5783
5768
double tmp1 = (1.0 - Y)/2.0;
5784
5769
double tmp2 = tmp1*tmp1;
5785
5770
5786
// Compute basisvalues.
5771
// Compute basisvalues
5787
5772
basisvalues[0] = 1.0;
5788
5773
basisvalues[1] = tmp0;
5789
5774
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5797
5782
basisvalues[4] *= std::sqrt(4.5);
5798
5783
basisvalues[3] *= std::sqrt(7.5);
5799
5784
5800
// Table(s) of coefficients.
5785
// Table(s) of coefficients
5801
5786
static const double coefficients0[6] = \
5802
5787
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
5803
5788
5804
// Compute value(s).
5789
// Compute value(s)
5805
5790
for (unsigned int r = 0; r < 6; r++)
5806
5791
{
5807
5792
values[1] += coefficients0[r]*basisvalues[r];
5811
5796
case 7:
5812
5797
{
5813
5798
5814
// Array of basisvalues.
5799
// Array of basisvalues
5815
5800
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5816
5801
5817
// Declare helper variables.
5802
// Declare helper variables
5818
5803
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5819
5804
double tmp1 = (1.0 - Y)/2.0;
5820
5805
double tmp2 = tmp1*tmp1;
5821
5806
5822
// Compute basisvalues.
5807
// Compute basisvalues
5823
5808
basisvalues[0] = 1.0;
5824
5809
basisvalues[1] = tmp0;
5825
5810
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5833
5818
basisvalues[4] *= std::sqrt(4.5);
5834
5819
basisvalues[3] *= std::sqrt(7.5);
5835
5820
5836
// Table(s) of coefficients.
5821
// Table(s) of coefficients
5837
5822
static const double coefficients0[6] = \
5838
5823
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
5839
5824
5840
// Compute value(s).
5825
// Compute value(s)
5841
5826
for (unsigned int r = 0; r < 6; r++)
5842
5827
{
5843
5828
values[1] += coefficients0[r]*basisvalues[r];
5847
5832
case 8:
5848
5833
{
5849
5834
5850
// Array of basisvalues.
5835
// Array of basisvalues
5851
5836
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5852
5837
5853
// Declare helper variables.
5838
// Declare helper variables
5854
5839
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5855
5840
double tmp1 = (1.0 - Y)/2.0;
5856
5841
double tmp2 = tmp1*tmp1;
5857
5842
5858
// Compute basisvalues.
5843
// Compute basisvalues
5859
5844
basisvalues[0] = 1.0;
5860
5845
basisvalues[1] = tmp0;
5861
5846
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5869
5854
basisvalues[4] *= std::sqrt(4.5);
5870
5855
basisvalues[3] *= std::sqrt(7.5);
5871
5856
5872
// Table(s) of coefficients.
5857
// Table(s) of coefficients
5873
5858
static const double coefficients0[6] = \
5874
5859
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
5875
5860
5876
// Compute value(s).
5861
// Compute value(s)
5877
5862
for (unsigned int r = 0; r < 6; r++)
5878
5863
{
5879
5864
values[1] += coefficients0[r]*basisvalues[r];
5883
5868
case 9:
5884
5869
{
5885
5870
5886
// Array of basisvalues.
5871
// Array of basisvalues
5887
5872
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5888
5873
5889
// Declare helper variables.
5874
// Declare helper variables
5890
5875
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5891
5876
double tmp1 = (1.0 - Y)/2.0;
5892
5877
double tmp2 = tmp1*tmp1;
5893
5878
5894
// Compute basisvalues.
5879
// Compute basisvalues
5895
5880
basisvalues[0] = 1.0;
5896
5881
basisvalues[1] = tmp0;
5897
5882
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5905
5890
basisvalues[4] *= std::sqrt(4.5);
5906
5891
basisvalues[3] *= std::sqrt(7.5);
5907
5892
5908
// Table(s) of coefficients.
5893
// Table(s) of coefficients
5909
5894
static const double coefficients0[6] = \
5910
5895
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
5911
5896
5912
// Compute value(s).
5897
// Compute value(s)
5913
5898
for (unsigned int r = 0; r < 6; r++)
5914
5899
{
5915
5900
values[1] += coefficients0[r]*basisvalues[r];
5919
5904
case 10:
5920
5905
{
5921
5906
5922
// Array of basisvalues.
5907
// Array of basisvalues
5923
5908
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5924
5909
5925
// Declare helper variables.
5910
// Declare helper variables
5926
5911
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5927
5912
double tmp1 = (1.0 - Y)/2.0;
5928
5913
double tmp2 = tmp1*tmp1;
5929
5914
5930
// Compute basisvalues.
5915
// Compute basisvalues
5931
5916
basisvalues[0] = 1.0;
5932
5917
basisvalues[1] = tmp0;
5933
5918
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5941
5926
basisvalues[4] *= std::sqrt(4.5);
5942
5927
basisvalues[3] *= std::sqrt(7.5);
5943
5928
5944
// Table(s) of coefficients.
5929
// Table(s) of coefficients
5945
5930
static const double coefficients0[6] = \
5946
5931
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
5947
5932
5948
// Compute value(s).
5933
// Compute value(s)
5949
5934
for (unsigned int r = 0; r < 6; r++)
5950
5935
{
5951
5936
values[1] += coefficients0[r]*basisvalues[r];
5955
5940
case 11:
5956
5941
{
5957
5942
5958
// Array of basisvalues.
5943
// Array of basisvalues
5959
5944
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5960
5945
5961
// Declare helper variables.
5946
// Declare helper variables
5962
5947
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5963
5948
double tmp1 = (1.0 - Y)/2.0;
5964
5949
double tmp2 = tmp1*tmp1;
5965
5950
5966
// Compute basisvalues.
5951
// Compute basisvalues
5967
5952
basisvalues[0] = 1.0;
5968
5953
basisvalues[1] = tmp0;
5969
5954
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5977
5962
basisvalues[4] *= std::sqrt(4.5);
5978
5963
basisvalues[3] *= std::sqrt(7.5);
5979
5964
5980
// Table(s) of coefficients.
5965
// Table(s) of coefficients
5981
5966
static const double coefficients0[6] = \
5982
5967
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
5983
5968
5984
// Compute value(s).
5969
// Compute value(s)
5985
5970
for (unsigned int r = 0; r < 6; r++)
5986
5971
{
5987
5972
values[1] += coefficients0[r]*basisvalues[r];
5991
5976
case 12:
5992
5977
{
5993
5978
5994
// Array of basisvalues.
5979
// Array of basisvalues
5995
5980
double basisvalues[3] = {0.0, 0.0, 0.0};
5996
5981
5997
// Declare helper variables.
5982
// Declare helper variables
5998
5983
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5999
5984
6000
// Compute basisvalues.
5985
// Compute basisvalues
6001
5986
basisvalues[0] = 1.0;
6002
5987
basisvalues[1] = tmp0;
6003
5988
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
6005
5990
basisvalues[2] *= std::sqrt(1.0);
6006
5991
basisvalues[1] *= std::sqrt(3.0);
6007
5992
6008
// Table(s) of coefficients.
5993
// Table(s) of coefficients
6009
5994
static const double coefficients0[3] = \
6010
5995
{0.47140452, -0.28867513, -0.16666667};
6011
5996
6012
// Compute value(s).
5997
// Compute value(s)
6013
5998
for (unsigned int r = 0; r < 3; r++)
6014
5999
{
6015
6000
values[2] += coefficients0[r]*basisvalues[r];
6019
6004
case 13:
6020
6005
{
6021
6006
6022
// Array of basisvalues.
6007
// Array of basisvalues
6023
6008
double basisvalues[3] = {0.0, 0.0, 0.0};
6024
6009
6025
// Declare helper variables.
6010
// Declare helper variables
6026
6011
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6027
6012
6028
// Compute basisvalues.
6013
// Compute basisvalues
6029
6014
basisvalues[0] = 1.0;
6030
6015
basisvalues[1] = tmp0;
6031
6016
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
6033
6018
basisvalues[2] *= std::sqrt(1.0);
6034
6019
basisvalues[1] *= std::sqrt(3.0);
6035
6020
6036
// Table(s) of coefficients.
6021
// Table(s) of coefficients
6037
6022
static const double coefficients0[3] = \
6038
6023
{0.47140452, 0.28867513, -0.16666667};
6039
6024
6040
// Compute value(s).
6025
// Compute value(s)
6041
6026
for (unsigned int r = 0; r < 3; r++)
6042
6027
{
6043
6028
values[2] += coefficients0[r]*basisvalues[r];
6047
6032
case 14:
6048
6033
{
6049
6034
6050
// Array of basisvalues.
6035
// Array of basisvalues
6051
6036
double basisvalues[3] = {0.0, 0.0, 0.0};
6052
6037
6053
// Declare helper variables.
6038
// Declare helper variables
6054
6039
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6055
6040
6056
// Compute basisvalues.
6041
// Compute basisvalues
6057
6042
basisvalues[0] = 1.0;
6058
6043
basisvalues[1] = tmp0;
6059
6044
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
6061
6046
basisvalues[2] *= std::sqrt(1.0);
6062
6047
basisvalues[1] *= std::sqrt(3.0);
6063
6048
6064
// Table(s) of coefficients.
6049
// Table(s) of coefficients
6065
6050
static const double coefficients0[3] = \
6066
6051
{0.47140452, 0.0, 0.33333333};
6067
6052
6068
// Compute value(s).
6053
// Compute value(s)
6069
6054
for (unsigned int r = 0; r < 3; r++)
6070
6055
{
6071
6056
values[2] += coefficients0[r]*basisvalues[r];
6076
6061
6077
6062
}
6078
6063
6079
/// Evaluate all basis functions at given point in cell
6064
/// Evaluate all basis functions at given point x in cell
6080
6065
virtual void evaluate_basis_all(double* values,
6081
const double* coordinates,
6082
const ufc::cell& c) const
6066
const double* x,
6067
const double* vertex_coordinates,
6068
int cell_orientation) const
6083
6069
{
6084
6070
// Helper variable to hold values of a single dof.
6085
6071
double dof_values[3] = {0.0, 0.0, 0.0};
6086
6072
6087
// Loop dofs and call evaluate_basis.
6073
// Loop dofs and call evaluate_basis
6088
6074
for (unsigned int r = 0; r < 15; r++)
6089
6075
{
6090
evaluate_basis(r, dof_values, coordinates, c);
6076
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
6091
6077
for (unsigned int s = 0; s < 3; s++)
6092
6078
{
6093
6079
values[r*3 + s] = dof_values[s];
6095
6081
}// end loop over 'r'
6096
6082
}
6097
6083
6098
/// Evaluate order n derivatives of basis function i at given point in cell
6099
virtual void evaluate_basis_derivatives(unsigned int i,
6100
unsigned int n,
6084
/// Evaluate order n derivatives of basis function i at given point x in cell
6085
virtual void evaluate_basis_derivatives(std::size_t i,
6086
std::size_t n,
6101
6087
double* values,
6102
const double* coordinates,
6103
const ufc::cell& c) const
6088
const double* x,
6089
const double* vertex_coordinates,
6090
int cell_orientation) const
6104
6091
{
6105
// Extract vertex coordinates
6106
const double * const * x = c.coordinates;
6107
6108
// Compute Jacobian of affine map from reference cell
6109
const double J_00 = x[1][0] - x[0][0];
6110
const double J_01 = x[2][0] - x[0][0];
6111
const double J_10 = x[1][1] - x[0][1];
6112
const double J_11 = x[2][1] - x[0][1];
6113
6114
// Compute determinant of Jacobian
6115
double detJ = J_00*J_11 - J_01*J_10;
6116
6117
// Compute inverse of Jacobian
6118
const double K_00 = J_11 / detJ;
6119
const double K_01 = -J_01 / detJ;
6120
const double K_10 = -J_10 / detJ;
6121
const double K_11 = J_00 / detJ;
6092
// Compute Jacobian
6093
double J[4];
6094
compute_jacobian_triangle_2d(J, vertex_coordinates);
6095
6096
// Compute Jacobian inverse and determinant
6097
double K[4];
6098
double detJ;
6099
compute_jacobian_inverse_triangle_2d(K, detJ, J);
6100
6122
6101
6123
6102
// Compute constants
6124
const double C0 = x[1][0] + x[2][0];
6125
const double C1 = x[1][1] + x[2][1];
6103
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
6104
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
6126
6105
6127
6106
// Get coordinates and map to the reference (FIAT) element
6128
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
6129
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
6107
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
6108
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
6130
6109
6131
6110
// Compute number of derivatives.
6132
6111
unsigned int num_derivatives = 1;
6163
6142
}
6164
6143
6165
6144
// Compute inverse of Jacobian
6166
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
6145
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
6167
6146
6168
6147
// Declare transformation matrix
6169
6148
// Declare pointer to two dimensional array and initialise
6197
6176
case 0:
6198
6177
{
6199
6178
6200
// Array of basisvalues.
6179
// Array of basisvalues
6201
6180
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6202
6181
6203
// Declare helper variables.
6182
// Declare helper variables
6204
6183
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6205
6184
double tmp1 = (1.0 - Y)/2.0;
6206
6185
double tmp2 = tmp1*tmp1;
6207
6186
6208
// Compute basisvalues.
6187
// Compute basisvalues
6209
6188
basisvalues[0] = 1.0;
6210
6189
basisvalues[1] = tmp0;
6211
6190
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6219
6198
basisvalues[4] *= std::sqrt(4.5);
6220
6199
basisvalues[3] *= std::sqrt(7.5);
6221
6200
6222
// Table(s) of coefficients.
6201
// Table(s) of coefficients
6223
6202
static const double coefficients0[6] = \
6224
6203
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
6225
6204
6363
6342
case 1:
6364
6343
{
6365
6344
6366
// Array of basisvalues.
6345
// Array of basisvalues
6367
6346
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6368
6347
6369
// Declare helper variables.
6348
// Declare helper variables
6370
6349
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6371
6350
double tmp1 = (1.0 - Y)/2.0;
6372
6351
double tmp2 = tmp1*tmp1;
6373
6352
6374
// Compute basisvalues.
6353
// Compute basisvalues
6375
6354
basisvalues[0] = 1.0;
6376
6355
basisvalues[1] = tmp0;
6377
6356
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6385
6364
basisvalues[4] *= std::sqrt(4.5);
6386
6365
basisvalues[3] *= std::sqrt(7.5);
6387
6366
6388
// Table(s) of coefficients.
6367
// Table(s) of coefficients
6389
6368
static const double coefficients0[6] = \
6390
6369
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
6391
6370
6529
6508
case 2:
6530
6509
{
6531
6510
6532
// Array of basisvalues.
6511
// Array of basisvalues
6533
6512
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6534
6513
6535
// Declare helper variables.
6514
// Declare helper variables
6536
6515
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6537
6516
double tmp1 = (1.0 - Y)/2.0;
6538
6517
double tmp2 = tmp1*tmp1;
6539
6518
6540
// Compute basisvalues.
6519
// Compute basisvalues
6541
6520
basisvalues[0] = 1.0;
6542
6521
basisvalues[1] = tmp0;
6543
6522
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6551
6530
basisvalues[4] *= std::sqrt(4.5);
6552
6531
basisvalues[3] *= std::sqrt(7.5);
6553
6532
6554
// Table(s) of coefficients.
6533
// Table(s) of coefficients
6555
6534
static const double coefficients0[6] = \
6556
6535
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
6557
6536
6695
6674
case 3:
6696
6675
{
6697
6676
6698
// Array of basisvalues.
6677
// Array of basisvalues
6699
6678
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6700
6679
6701
// Declare helper variables.
6680
// Declare helper variables
6702
6681
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6703
6682
double tmp1 = (1.0 - Y)/2.0;
6704
6683
double tmp2 = tmp1*tmp1;
6705
6684
6706
// Compute basisvalues.
6685
// Compute basisvalues
6707
6686
basisvalues[0] = 1.0;
6708
6687
basisvalues[1] = tmp0;
6709
6688
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6717
6696
basisvalues[4] *= std::sqrt(4.5);
6718
6697
basisvalues[3] *= std::sqrt(7.5);
6719
6698
6720
// Table(s) of coefficients.
6699
// Table(s) of coefficients
6721
6700
static const double coefficients0[6] = \
6722
6701
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
6723
6702
6861
6840
case 4:
6862
6841
{
6863
6842
6864
// Array of basisvalues.
6843
// Array of basisvalues
6865
6844
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6866
6845
6867
// Declare helper variables.
6846
// Declare helper variables
6868
6847
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6869
6848
double tmp1 = (1.0 - Y)/2.0;
6870
6849
double tmp2 = tmp1*tmp1;
6871
6850
6872
// Compute basisvalues.
6851
// Compute basisvalues
6873
6852
basisvalues[0] = 1.0;
6874
6853
basisvalues[1] = tmp0;
6875
6854
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6883
6862
basisvalues[4] *= std::sqrt(4.5);
6884
6863
basisvalues[3] *= std::sqrt(7.5);
6885
6864
6886
// Table(s) of coefficients.
6865
// Table(s) of coefficients
6887
6866
static const double coefficients0[6] = \
6888
6867
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
6889
6868
7027
7006
case 5:
7028
7007
{
7029
7008
7030
// Array of basisvalues.
7009
// Array of basisvalues
7031
7010
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
7032
7011
7033
// Declare helper variables.
7012
// Declare helper variables
7034
7013
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
7035
7014
double tmp1 = (1.0 - Y)/2.0;
7036
7015
double tmp2 = tmp1*tmp1;
7037
7016
7038
// Compute basisvalues.
7017
// Compute basisvalues
7039
7018
basisvalues[0] = 1.0;
7040
7019
basisvalues[1] = tmp0;
7041
7020
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
7049
7028
basisvalues[4] *= std::sqrt(4.5);
7050
7029
basisvalues[3] *= std::sqrt(7.5);
7051
7030
7052
// Table(s) of coefficients.
7031
// Table(s) of coefficients
7053
7032
static const double coefficients0[6] = \
7054
7033
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
7055
7034
7193
7172
case 6:
7194
7173
{
7195
7174
7196
// Array of basisvalues.
7175
// Array of basisvalues
7197
7176
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
7198
7177
7199
// Declare helper variables.
7178
// Declare helper variables
7200
7179
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
7201
7180
double tmp1 = (1.0 - Y)/2.0;
7202
7181
double tmp2 = tmp1*tmp1;
7203
7182
7204
// Compute basisvalues.
7183
// Compute basisvalues
7205
7184
basisvalues[0] = 1.0;
7206
7185
basisvalues[1] = tmp0;
7207
7186
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
7215
7194
basisvalues[4] *= std::sqrt(4.5);
7216
7195
basisvalues[3] *= std::sqrt(7.5);
7217
7196
7218
// Table(s) of coefficients.
7197
// Table(s) of coefficients
7219
7198
static const double coefficients0[6] = \
7220
7199
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
7221
7200
7359
7338
case 7:
7360
7339
{
7361
7340
7362
// Array of basisvalues.
7341
// Array of basisvalues
7363
7342
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
7364
7343
7365
// Declare helper variables.
7344
// Declare helper variables
7366
7345
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
7367
7346
double tmp1 = (1.0 - Y)/2.0;
7368
7347
double tmp2 = tmp1*tmp1;
7369
7348
7370
// Compute basisvalues.
7349
// Compute basisvalues
7371
7350
basisvalues[0] = 1.0;
7372
7351
basisvalues[1] = tmp0;
7373
7352
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
7381
7360
basisvalues[4] *= std::sqrt(4.5);
7382
7361
basisvalues[3] *= std::sqrt(7.5);
7383
7362
7384
// Table(s) of coefficients.
7363
// Table(s) of coefficients
7385
7364
static const double coefficients0[6] = \
7386
7365
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
7387
7366
7525
7504
case 8:
7526
7505
{
7527
7506
7528
// Array of basisvalues.
7507
// Array of basisvalues
7529
7508
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
7530
7509
7531
// Declare helper variables.
7510
// Declare helper variables
7532
7511
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
7533
7512
double tmp1 = (1.0 - Y)/2.0;
7534
7513
double tmp2 = tmp1*tmp1;
7535
7514
7536
// Compute basisvalues.
7515
// Compute basisvalues
7537
7516
basisvalues[0] = 1.0;
7538
7517
basisvalues[1] = tmp0;
7539
7518
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
7547
7526
basisvalues[4] *= std::sqrt(4.5);
7548
7527
basisvalues[3] *= std::sqrt(7.5);
7549
7528
7550
// Table(s) of coefficients.
7529
// Table(s) of coefficients
7551
7530
static const double coefficients0[6] = \
7552
7531
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
7553
7532
7691
7670
case 9:
7692
7671
{
7693
7672
7694
// Array of basisvalues.
7673
// Array of basisvalues
7695
7674
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
7696
7675
7697
// Declare helper variables.
7676
// Declare helper variables
7698
7677
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
7699
7678
double tmp1 = (1.0 - Y)/2.0;
7700
7679
double tmp2 = tmp1*tmp1;
7701
7680
7702
// Compute basisvalues.
7681
// Compute basisvalues
7703
7682
basisvalues[0] = 1.0;
7704
7683
basisvalues[1] = tmp0;
7705
7684
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
7713
7692
basisvalues[4] *= std::sqrt(4.5);
7714
7693
basisvalues[3] *= std::sqrt(7.5);
7715
7694
7716
// Table(s) of coefficients.
7695
// Table(s) of coefficients
7717
7696
static const double coefficients0[6] = \
7718
7697
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
7719
7698
7857
7836
case 10:
7858
7837
{
7859
7838
7860
// Array of basisvalues.
7839
// Array of basisvalues
7861
7840
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
7862
7841
7863
// Declare helper variables.
7842
// Declare helper variables
7864
7843
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
7865
7844
double tmp1 = (1.0 - Y)/2.0;
7866
7845
double tmp2 = tmp1*tmp1;
7867
7846
7868
// Compute basisvalues.
7847
// Compute basisvalues
7869
7848
basisvalues[0] = 1.0;
7870
7849
basisvalues[1] = tmp0;
7871
7850
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
7879
7858
basisvalues[4] *= std::sqrt(4.5);
7880
7859
basisvalues[3] *= std::sqrt(7.5);
7881
7860
7882
// Table(s) of coefficients.
7861
// Table(s) of coefficients
7883
7862
static const double coefficients0[6] = \
7884
7863
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
7885
7864
8023
8002
case 11:
8024
8003
{
8025
8004
8026
// Array of basisvalues.
8005
// Array of basisvalues
8027
8006
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
8028
8007
8029
// Declare helper variables.
8008
// Declare helper variables
8030
8009
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
8031
8010
double tmp1 = (1.0 - Y)/2.0;
8032
8011
double tmp2 = tmp1*tmp1;
8033
8012
8034
// Compute basisvalues.
8013
// Compute basisvalues
8035
8014
basisvalues[0] = 1.0;
8036
8015
basisvalues[1] = tmp0;
8037
8016
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
8045
8024
basisvalues[4] *= std::sqrt(4.5);
8046
8025
basisvalues[3] *= std::sqrt(7.5);
8047
8026
8048
// Table(s) of coefficients.
8027
// Table(s) of coefficients
8049
8028
static const double coefficients0[6] = \
8050
8029
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
8051
8030
8189
8168
case 12:
8190
8169
{
8191
8170
8192
// Array of basisvalues.
8171
// Array of basisvalues
8193
8172
double basisvalues[3] = {0.0, 0.0, 0.0};
8194
8173
8195
// Declare helper variables.
8174
// Declare helper variables
8196
8175
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
8197
8176
8198
// Compute basisvalues.
8177
// Compute basisvalues
8199
8178
basisvalues[0] = 1.0;
8200
8179
basisvalues[1] = tmp0;
8201
8180
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
8203
8182
basisvalues[2] *= std::sqrt(1.0);
8204
8183
basisvalues[1] *= std::sqrt(3.0);
8205
8184
8206
// Table(s) of coefficients.
8185
// Table(s) of coefficients
8207
8186
static const double coefficients0[3] = \
8208
8187
{0.47140452, -0.28867513, -0.16666667};
8209
8188
8335
8314
case 13:
8336
8315
{
8337
8316
8338
// Array of basisvalues.
8317
// Array of basisvalues
8339
8318
double basisvalues[3] = {0.0, 0.0, 0.0};
8340
8319
8341
// Declare helper variables.
8320
// Declare helper variables
8342
8321
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
8343
8322
8344
// Compute basisvalues.
8323
// Compute basisvalues
8345
8324
basisvalues[0] = 1.0;
8346
8325
basisvalues[1] = tmp0;
8347
8326
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
8349
8328
basisvalues[2] *= std::sqrt(1.0);
8350
8329
basisvalues[1] *= std::sqrt(3.0);
8351
8330
8352
// Table(s) of coefficients.
8331
// Table(s) of coefficients
8353
8332
static const double coefficients0[3] = \
8354
8333
{0.47140452, 0.28867513, -0.16666667};
8355
8334
8481
8460
case 14:
8482
8461
{
8483
8462
8484
// Array of basisvalues.
8463
// Array of basisvalues
8485
8464
double basisvalues[3] = {0.0, 0.0, 0.0};
8486
8465
8487
// Declare helper variables.
8466
// Declare helper variables
8488
8467
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
8489
8468
8490
// Compute basisvalues.
8469
// Compute basisvalues
8491
8470
basisvalues[0] = 1.0;
8492
8471
basisvalues[1] = tmp0;
8493
8472
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
8495
8474
basisvalues[2] *= std::sqrt(1.0);
8496
8475
basisvalues[1] *= std::sqrt(3.0);
8497
8476
8498
// Table(s) of coefficients.
8477
// Table(s) of coefficients
8499
8478
static const double coefficients0[3] = \
8500
8479
{0.47140452, 0.0, 0.33333333};
8501
8480
8628
8607
8629
8608
}
8630
8609
8631
/// Evaluate order n derivatives of all basis functions at given point in cell
8632
virtual void evaluate_basis_derivatives_all(unsigned int n,
8610
/// Evaluate order n derivatives of all basis functions at given point x in cell
8611
virtual void evaluate_basis_derivatives_all(std::size_t n,
8633
8612
double* values,
8634
const double* coordinates,
8635
const ufc::cell& c) const
8613
const double* x,
8614
const double* vertex_coordinates,
8615
int cell_orientation) const
8636
8616
{
8637
8617
// Compute number of derivatives.
8638
8618
unsigned int num_derivatives = 1;
8651
8631
// Loop dofs and call evaluate_basis_derivatives.
8652
8632
for (unsigned int r = 0; r < 15; r++)
8653
8633
{
8654
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
8634
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
8655
8635
for (unsigned int s = 0; s < 3*num_derivatives; s++)
8656
8636
{
8657
8637
values[r*3*num_derivatives + s] = dof_values[s];
8663
8643
}
8664
8644
8665
8645
/// Evaluate linear functional for dof i on the function f
8666
virtual double evaluate_dof(unsigned int i,
8646
virtual double evaluate_dof(std::size_t i,
8667
8647
const ufc::function& f,
8648
const double* vertex_coordinates,
8649
int cell_orientation,
8668
8650
const ufc::cell& c) const
8669
8651
{
8670
// Declare variables for result of evaluation.
8652
// Declare variables for result of evaluation
8671
8653
double vals[3];
8672
8654
8673
// Declare variable for physical coordinates.
8655
// Declare variable for physical coordinates
8674
8656
double y[2];
8675
const double * const * x = c.coordinates;
8676
8657
switch (i)
8677
8658
{
8678
8659
case 0:
8679
8660
{
8680
y[0] = x[0][0];
8681
y[1] = x[0][1];
8661
y[0] = vertex_coordinates[0];
8662
y[1] = vertex_coordinates[1];
8682
8663
f.evaluate(vals, y, c);
8683
8664
return vals[0];
8684
8665
break;
8685
8666
}
8686
8667
case 1:
8687
8668
{
8688
y[0] = x[1][0];
8689
y[1] = x[1][1];
8669
y[0] = vertex_coordinates[2];
8670
y[1] = vertex_coordinates[3];
8690
8671
f.evaluate(vals, y, c);
8691
8672
return vals[0];
8692
8673
break;
8693
8674
}
8694
8675
case 2:
8695
8676
{
8696
y[0] = x[2][0];
8697
y[1] = x[2][1];
8677
y[0] = vertex_coordinates[4];
8678
y[1] = vertex_coordinates[5];
8698
8679
f.evaluate(vals, y, c);
8699
8680
return vals[0];
8700
8681
break;
8701
8682
}
8702
8683
case 3:
8703
8684
{
8704
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
8705
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
8685
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
8686
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
8706
8687
f.evaluate(vals, y, c);
8707
8688
return vals[0];
8708
8689
break;
8709
8690
}
8710
8691
case 4:
8711
8692
{
8712
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
8713
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
8693
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
8694
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
8714
8695
f.evaluate(vals, y, c);
8715
8696
return vals[0];
8716
8697
break;
8717
8698
}
8718
8699
case 5:
8719
8700
{
8720
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
8721
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
8701
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
8702
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
8722
8703
f.evaluate(vals, y, c);
8723
8704
return vals[0];
8724
8705
break;
8725
8706
}
8726
8707
case 6:
8727
8708
{
8728
y[0] = x[0][0];
8729
y[1] = x[0][1];
8709
y[0] = vertex_coordinates[0];
8710
y[1] = vertex_coordinates[1];
8730
8711
f.evaluate(vals, y, c);
8731
8712
return vals[1];
8732
8713
break;
8733
8714
}
8734
8715
case 7:
8735
8716
{
8736
y[0] = x[1][0];
8737
y[1] = x[1][1];
8717
y[0] = vertex_coordinates[2];
8718
y[1] = vertex_coordinates[3];
8738
8719
f.evaluate(vals, y, c);
8739
8720
return vals[1];
8740
8721
break;
8741
8722
}
8742
8723
case 8:
8743
8724
{
8744
y[0] = x[2][0];
8745
y[1] = x[2][1];
8725
y[0] = vertex_coordinates[4];
8726
y[1] = vertex_coordinates[5];
8746
8727
f.evaluate(vals, y, c);
8747
8728
return vals[1];
8748
8729
break;
8749
8730
}
8750
8731
case 9:
8751
8732
{
8752
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
8753
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
8733
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
8734
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
8754
8735
f.evaluate(vals, y, c);
8755
8736
return vals[1];
8756
8737
break;
8757
8738
}
8758
8739
case 10:
8759
8740
{
8760
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
8761
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
8741
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
8742
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
8762
8743
f.evaluate(vals, y, c);
8763
8744
return vals[1];
8764
8745
break;
8765
8746
}
8766
8747
case 11:
8767
8748
{
8768
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
8769
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
8749
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
8750
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
8770
8751
f.evaluate(vals, y, c);
8771
8752
return vals[1];
8772
8753
break;
8773
8754
}
8774
8755
case 12:
8775
8756
{
8776
y[0] = x[0][0];
8777
y[1] = x[0][1];
8757
y[0] = vertex_coordinates[0];
8758
y[1] = vertex_coordinates[1];
8778
8759
f.evaluate(vals, y, c);
8779
8760
return vals[2];
8780
8761
break;
8781
8762
}
8782
8763
case 13:
8783
8764
{
8784
y[0] = x[1][0];
8785
y[1] = x[1][1];
8765
y[0] = vertex_coordinates[2];
8766
y[1] = vertex_coordinates[3];
8786
8767
f.evaluate(vals, y, c);
8787
8768
return vals[2];
8788
8769
break;
8789
8770
}
8790
8771
case 14:
8791
8772
{
8792
y[0] = x[2][0];
8793
y[1] = x[2][1];
8773
y[0] = vertex_coordinates[4];
8774
y[1] = vertex_coordinates[5];
8794
8775
f.evaluate(vals, y, c);
8795
8776
return vals[2];
8796
8777
break;
8803
8784
/// Evaluate linear functionals for all dofs on the function f
8804
8785
virtual void evaluate_dofs(double* values,
8805
8786
const ufc::function& f,
8787
const double* vertex_coordinates,
8788
int cell_orientation,
8806
8789
const ufc::cell& c) const
8807
8790
{
8808
// Declare variables for result of evaluation.
8791
// Declare variables for result of evaluation
8809
8792
double vals[3];
8810
8793
8811
// Declare variable for physical coordinates.
8794
// Declare variable for physical coordinates
8812
8795
double y[2];
8813
const double * const * x = c.coordinates;
8814
y[0] = x[0][0];
8815
y[1] = x[0][1];
8796
y[0] = vertex_coordinates[0];
8797
y[1] = vertex_coordinates[1];
8816
8798
f.evaluate(vals, y, c);
8817
8799
values[0] = vals[0];
8818
y[0] = x[1][0];
8819
y[1] = x[1][1];
8800
y[0] = vertex_coordinates[2];
8801
y[1] = vertex_coordinates[3];
8820
8802
f.evaluate(vals, y, c);
8821
8803
values[1] = vals[0];
8822
y[0] = x[2][0];
8823
y[1] = x[2][1];
8804
y[0] = vertex_coordinates[4];
8805
y[1] = vertex_coordinates[5];
8824
8806
f.evaluate(vals, y, c);
8825
8807
values[2] = vals[0];
8826
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
8827
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
8808
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
8809
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
8828
8810
f.evaluate(vals, y, c);
8829
8811
values[3] = vals[0];
8830
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
8831
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
8812
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
8813
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
8832
8814
f.evaluate(vals, y, c);
8833
8815
values[4] = vals[0];
8834
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
8835
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
8816
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
8817
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
8836
8818
f.evaluate(vals, y, c);
8837
8819
values[5] = vals[0];
8838
y[0] = x[0][0];
8839
y[1] = x[0][1];
8820
y[0] = vertex_coordinates[0];
8821
y[1] = vertex_coordinates[1];
8840
8822
f.evaluate(vals, y, c);
8841
8823
values[6] = vals[1];
8842
y[0] = x[1][0];
8843
y[1] = x[1][1];
8824
y[0] = vertex_coordinates[2];
8825
y[1] = vertex_coordinates[3];
8844
8826
f.evaluate(vals, y, c);
8845
8827
values[7] = vals[1];
8846
y[0] = x[2][0];
8847
y[1] = x[2][1];
8828
y[0] = vertex_coordinates[4];
8829
y[1] = vertex_coordinates[5];
8848
8830
f.evaluate(vals, y, c);
8849
8831
values[8] = vals[1];
8850
y[0] = 0.5*x[1][0] + 0.5*x[2][0];
8851
y[1] = 0.5*x[1][1] + 0.5*x[2][1];
8832
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
8833
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
8852
8834
f.evaluate(vals, y, c);
8853
8835
values[9] = vals[1];
8854
y[0] = 0.5*x[0][0] + 0.5*x[2][0];
8855
y[1] = 0.5*x[0][1] + 0.5*x[2][1];
8836
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
8837
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
8856
8838
f.evaluate(vals, y, c);
8857
8839
values[10] = vals[1];
8858
y[0] = 0.5*x[0][0] + 0.5*x[1][0];
8859
y[1] = 0.5*x[0][1] + 0.5*x[1][1];
8840
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
8841
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
8860
8842
f.evaluate(vals, y, c);
8861
8843
values[11] = vals[1];
8862
y[0] = x[0][0];
8863
y[1] = x[0][1];
8844
y[0] = vertex_coordinates[0];
8845
y[1] = vertex_coordinates[1];
8864
8846
f.evaluate(vals, y, c);
8865
8847
values[12] = vals[2];
8866
y[0] = x[1][0];
8867
y[1] = x[1][1];
8848
y[0] = vertex_coordinates[2];
8849
y[1] = vertex_coordinates[3];
8868
8850
f.evaluate(vals, y, c);
8869
8851
values[13] = vals[2];
8870
y[0] = x[2][0];
8871
y[1] = x[2][1];
8852
y[0] = vertex_coordinates[4];
8853
y[1] = vertex_coordinates[5];
8872
8854
f.evaluate(vals, y, c);
8873
8855
values[14] = vals[2];
8874
8856
}
8876
8858
/// Interpolate vertex values from dof values
8877
8859
virtual void interpolate_vertex_values(double* vertex_values,
8878
8860
const double* dof_values,
8861
const double* vertex_coordinates,
8862
int cell_orientation,
8879
8863
const ufc::cell& c) const
8880
8864
{
8881
8865
// Evaluate function and change variables
8897
8881
const double* xhat,
8898
8882
const ufc::cell& c) const
8899
8883
{
8900
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
8884
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
8901
8885
}
8902
8886
8903
8887
/// Map from coordinate x in cell to coordinate xhat in reference cell
8905
8889
const double* x,
8906
8890
const ufc::cell& c) const
8907
8891
{
8908
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
8892
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
8909
8893
}
8910
8894
8911
8895
/// Return the number of sub elements (for a mixed element)
8912
virtual unsigned int num_sub_elements() const
8896
virtual std::size_t num_sub_elements() const
8913
8897
{
8914
8898
return 2;
8915
8899
}
8916
8900
8917
8901
/// Create a new finite element for sub element i (for a mixed element)
8918
virtual ufc::finite_element* create_sub_element(unsigned int i) const
8902
virtual ufc::finite_element* create_sub_element(std::size_t i) const
8919
8903
{
8920
8904
switch (i)
8921
8905
{
8947
8931
8948
8932
class stokes_dofmap_0: public ufc::dofmap
8949
8933
{
8950
private:
8951
8952
unsigned int _global_dimension;
8953
8934
public:
8954
8935
8955
8936
/// Constructor
8956
8937
stokes_dofmap_0() : ufc::dofmap()
8957
8938
{
8958
_global_dimension = 0;
8939
// Do nothing
8959
8940
}
8960
8941
8961
8942
/// Destructor
8967
8948
/// Return a string identifying the dofmap
8968
8949
virtual const char* signature() const
8969
8950
{
8970
return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', Space(2)), 2, None)";
8951
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
8971
8952
}
8972
8953
8973
8954
/// Return true iff mesh entities of topological dimension d are needed
8974
virtual bool needs_mesh_entities(unsigned int d) const
8955
virtual bool needs_mesh_entities(std::size_t d) const
8975
8956
{
8976
8957
switch (d)
8977
8958
{
8995
8976
return false;
8996
8977
}
8997
8978
8998
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
8999
virtual bool init_mesh(const ufc::mesh& m)
9000
{
9001
_global_dimension = m.num_entities[0] + m.num_entities[1];
9002
return false;
9003
}
9004
9005
/// Initialize dofmap for given cell
9006
virtual void init_cell(const ufc::mesh& m,
9007
const ufc::cell& c)
9008
{
9009
// Do nothing
9010
}
9011
9012
/// Finish initialization of dofmap for cells
9013
virtual void init_cell_finalize()
9014
{
9015
// Do nothing
9016
}
9017
9018
8979
/// Return the topological dimension of the associated cell shape
9019
virtual unsigned int topological_dimension() const
8980
virtual std::size_t topological_dimension() const
9020
8981
{
9021
8982
return 2;
9022
8983
}
9023
8984
9024
8985
/// Return the geometric dimension of the associated cell shape
9025
virtual unsigned int geometric_dimension() const
8986
virtual std::size_t geometric_dimension() const
9026
8987
{
9027
8988
return 2;
9028
8989
}
9029
8990
9030
8991
/// Return the dimension of the global finite element function space
9031
virtual unsigned int global_dimension() const
8992
virtual std::size_t global_dimension(const std::vector<std::size_t>&
8993
num_global_entities) const
9032
8994
{
9033
return _global_dimension;
8995
return num_global_entities[0] + num_global_entities[1];
9034
8996
}
9035
8997
9036
8998
/// Return the dimension of the local finite element function space for a cell
9037
virtual unsigned int local_dimension(const ufc::cell& c) const
8999
virtual std::size_t local_dimension(const ufc::cell& c) const
9038
9000
{
9039
9001
return 6;
9040
9002
}
9041
9003
9042
9004
/// Return the maximum dimension of the local finite element function space
9043
virtual unsigned int max_local_dimension() const
9005
virtual std::size_t max_local_dimension() const
9044
9006
{
9045
9007
return 6;
9046
9008
}
9047
9009
9048
9010
/// Return the number of dofs on each cell facet
9049
virtual unsigned int num_facet_dofs() const
9011
virtual std::size_t num_facet_dofs() const
9050
9012
{
9051
9013
return 3;
9052
9014
}
9053
9015
9054
9016
/// Return the number of dofs associated with each cell entity of dimension d
9055
virtual unsigned int num_entity_dofs(unsigned int d) const
9017
virtual std::size_t num_entity_dofs(std::size_t d) const
9056
9018
{
9057
9019
switch (d)
9058
9020
{
9077
9039
}
9078
9040
9079
9041
/// Tabulate the local-to-global mapping of dofs on a cell
9080
virtual void tabulate_dofs(unsigned int* dofs,
9081
const ufc::mesh& m,
9042
virtual void tabulate_dofs(std::size_t* dofs,
9043
const std::vector<std::size_t>& num_global_entities,
9082
9044
const ufc::cell& c) const
9083
9045
{
9084
9046
unsigned int offset = 0;
9085
9047
dofs[0] = offset + c.entity_indices[0][0];
9086
9048
dofs[1] = offset + c.entity_indices[0][1];
9087
9049
dofs[2] = offset + c.entity_indices[0][2];
9088
offset += m.num_entities[0];
9050
offset += num_global_entities[0];
9089
9051
dofs[3] = offset + c.entity_indices[1][0];
9090
9052
dofs[4] = offset + c.entity_indices[1][1];
9091
9053
dofs[5] = offset + c.entity_indices[1][2];
9092
offset += m.num_entities[1];
9054
offset += num_global_entities[1];
9093
9055
}
9094
9056
9095
9057
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
9096
virtual void tabulate_facet_dofs(unsigned int* dofs,
9097
unsigned int facet) const
9058
virtual void tabulate_facet_dofs(std::size_t* dofs,
9059
std::size_t facet) const
9098
9060
{
9099
9061
switch (facet)
9100
9062
{
9124
9086
}
9125
9087
9126
9088
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
9127
virtual void tabulate_entity_dofs(unsigned int* dofs,
9128
unsigned int d, unsigned int i) const
9089
virtual void tabulate_entity_dofs(std::size_t* dofs,
9090
std::size_t d, std::size_t i) const
9129
9091
{
9130
9092
if (d > 2)
9131
9093
{
9200
9162
}
9201
9163
9202
9164
/// Tabulate the coordinates of all dofs on a cell
9203
virtual void tabulate_coordinates(double** coordinates,
9204
const ufc::cell& c) const
9165
virtual void tabulate_coordinates(double** dof_coordinates,
9166
const double* vertex_coordinates) const
9205
9167
{
9206
const double * const * x = c.coordinates;
9207
9208
coordinates[0][0] = x[0][0];
9209
coordinates[0][1] = x[0][1];
9210
coordinates[1][0] = x[1][0];
9211
coordinates[1][1] = x[1][1];
9212
coordinates[2][0] = x[2][0];
9213
coordinates[2][1] = x[2][1];
9214
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
9215
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
9216
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
9217
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
9218
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
9219
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
9168
dof_coordinates[0][0] = vertex_coordinates[0];
9169
dof_coordinates[0][1] = vertex_coordinates[1];
9170
dof_coordinates[1][0] = vertex_coordinates[2];
9171
dof_coordinates[1][1] = vertex_coordinates[3];
9172
dof_coordinates[2][0] = vertex_coordinates[4];
9173
dof_coordinates[2][1] = vertex_coordinates[5];
9174
dof_coordinates[3][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
9175
dof_coordinates[3][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
9176
dof_coordinates[4][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
9177
dof_coordinates[4][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
9178
dof_coordinates[5][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
9179
dof_coordinates[5][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
9220
9180
}
9221
9181
9222
9182
/// Return the number of sub dofmaps (for a mixed element)
9223
virtual unsigned int num_sub_dofmaps() const
9183
virtual std::size_t num_sub_dofmaps() const
9224
9184
{
9225
9185
return 0;
9226
9186
}
9227
9187
9228
9188
/// Create a new dofmap for sub dofmap i (for a mixed element)
9229
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
9189
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
9230
9190
{
9231
9191
return 0;
9232
9192
}
9244
9204
9245
9205
class stokes_dofmap_1: public ufc::dofmap
9246
9206
{
9247
private:
9248
9249
unsigned int _global_dimension;
9250
9207
public:
9251
9208
9252
9209
/// Constructor
9253
9210
stokes_dofmap_1() : ufc::dofmap()
9254
9211
{
9255
_global_dimension = 0;
9212
// Do nothing
9256
9213
}
9257
9214
9258
9215
/// Destructor
9264
9221
/// Return a string identifying the dofmap
9265
9222
virtual const char* signature() const
9266
9223
{
9267
return "FFC dofmap for VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None)";
9224
return "FFC dofmap for VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None)";
9268
9225
}
9269
9226
9270
9227
/// Return true iff mesh entities of topological dimension d are needed
9271
virtual bool needs_mesh_entities(unsigned int d) const
9228
virtual bool needs_mesh_entities(std::size_t d) const
9272
9229
{
9273
9230
switch (d)
9274
9231
{
9292
9249
return false;
9293
9250
}
9294
9251
9295
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
9296
virtual bool init_mesh(const ufc::mesh& m)
9297
{
9298
_global_dimension = 2*m.num_entities[0] + 2*m.num_entities[1];
9299
return false;
9300
}
9301
9302
/// Initialize dofmap for given cell
9303
virtual void init_cell(const ufc::mesh& m,
9304
const ufc::cell& c)
9305
{
9306
// Do nothing
9307
}
9308
9309
/// Finish initialization of dofmap for cells
9310
virtual void init_cell_finalize()
9311
{
9312
// Do nothing
9313
}
9314
9315
9252
/// Return the topological dimension of the associated cell shape
9316
virtual unsigned int topological_dimension() const
9253
virtual std::size_t topological_dimension() const
9317
9254
{
9318
9255
return 2;
9319
9256
}
9320
9257
9321
9258
/// Return the geometric dimension of the associated cell shape
9322
virtual unsigned int geometric_dimension() const
9259
virtual std::size_t geometric_dimension() const
9323
9260
{
9324
9261
return 2;
9325
9262
}
9326
9263
9327
9264
/// Return the dimension of the global finite element function space
9328
virtual unsigned int global_dimension() const
9265
virtual std::size_t global_dimension(const std::vector<std::size_t>&
9266
num_global_entities) const
9329
9267
{
9330
return _global_dimension;
9268
return 2*num_global_entities[0] + 2*num_global_entities[1];
9331
9269
}
9332
9270
9333
9271
/// Return the dimension of the local finite element function space for a cell
9334
virtual unsigned int local_dimension(const ufc::cell& c) const
9272
virtual std::size_t local_dimension(const ufc::cell& c) const
9335
9273
{
9336
9274
return 12;
9337
9275
}
9338
9276
9339
9277
/// Return the maximum dimension of the local finite element function space
9340
virtual unsigned int max_local_dimension() const
9278
virtual std::size_t max_local_dimension() const
9341
9279
{
9342
9280
return 12;
9343
9281
}
9344
9282
9345
9283
/// Return the number of dofs on each cell facet
9346
virtual unsigned int num_facet_dofs() const
9284
virtual std::size_t num_facet_dofs() const
9347
9285
{
9348
9286
return 6;
9349
9287
}
9350
9288
9351
9289
/// Return the number of dofs associated with each cell entity of dimension d
9352
virtual unsigned int num_entity_dofs(unsigned int d) const
9290
virtual std::size_t num_entity_dofs(std::size_t d) const
9353
9291
{
9354
9292
switch (d)
9355
9293
{
9374
9312
}
9375
9313
9376
9314
/// Tabulate the local-to-global mapping of dofs on a cell
9377
virtual void tabulate_dofs(unsigned int* dofs,
9378
const ufc::mesh& m,
9315
virtual void tabulate_dofs(std::size_t* dofs,
9316
const std::vector<std::size_t>& num_global_entities,
9379
9317
const ufc::cell& c) const
9380
9318
{
9381
9319
unsigned int offset = 0;
9382
9320
dofs[0] = offset + c.entity_indices[0][0];
9383
9321
dofs[1] = offset + c.entity_indices[0][1];
9384
9322
dofs[2] = offset + c.entity_indices[0][2];
9385
offset += m.num_entities[0];
9323
offset += num_global_entities[0];
9386
9324
dofs[3] = offset + c.entity_indices[1][0];
9387
9325
dofs[4] = offset + c.entity_indices[1][1];
9388
9326
dofs[5] = offset + c.entity_indices[1][2];
9389
offset += m.num_entities[1];
9327
offset += num_global_entities[1];
9390
9328
dofs[6] = offset + c.entity_indices[0][0];
9391
9329
dofs[7] = offset + c.entity_indices[0][1];
9392
9330
dofs[8] = offset + c.entity_indices[0][2];
9393
offset += m.num_entities[0];
9331
offset += num_global_entities[0];
9394
9332
dofs[9] = offset + c.entity_indices[1][0];
9395
9333
dofs[10] = offset + c.entity_indices[1][1];
9396
9334
dofs[11] = offset + c.entity_indices[1][2];
9397
offset += m.num_entities[1];
9335
offset += num_global_entities[1];
9398
9336
}
9399
9337
9400
9338
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
9401
virtual void tabulate_facet_dofs(unsigned int* dofs,
9402
unsigned int facet) const
9339
virtual void tabulate_facet_dofs(std::size_t* dofs,
9340
std::size_t facet) const
9403
9341
{
9404
9342
switch (facet)
9405
9343
{
9438
9376
}
9439
9377
9440
9378
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
9441
virtual void tabulate_entity_dofs(unsigned int* dofs,
9442
unsigned int d, unsigned int i) const
9379
virtual void tabulate_entity_dofs(std::size_t* dofs,
9380
std::size_t d, std::size_t i) const
9443
9381
{
9444
9382
if (d > 2)
9445
9383
{
9520
9458
}
9521
9459
9522
9460
/// Tabulate the coordinates of all dofs on a cell
9523
virtual void tabulate_coordinates(double** coordinates,
9524
const ufc::cell& c) const
9461
virtual void tabulate_coordinates(double** dof_coordinates,
9462
const double* vertex_coordinates) const
9525
9463
{
9526
const double * const * x = c.coordinates;
9527
9528
coordinates[0][0] = x[0][0];
9529
coordinates[0][1] = x[0][1];
9530
coordinates[1][0] = x[1][0];
9531
coordinates[1][1] = x[1][1];
9532
coordinates[2][0] = x[2][0];
9533
coordinates[2][1] = x[2][1];
9534
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
9535
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
9536
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
9537
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
9538
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
9539
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
9540
coordinates[6][0] = x[0][0];
9541
coordinates[6][1] = x[0][1];
9542
coordinates[7][0] = x[1][0];
9543
coordinates[7][1] = x[1][1];
9544
coordinates[8][0] = x[2][0];
9545
coordinates[8][1] = x[2][1];
9546
coordinates[9][0] = 0.5*x[1][0] + 0.5*x[2][0];
9547
coordinates[9][1] = 0.5*x[1][1] + 0.5*x[2][1];
9548
coordinates[10][0] = 0.5*x[0][0] + 0.5*x[2][0];
9549
coordinates[10][1] = 0.5*x[0][1] + 0.5*x[2][1];
9550
coordinates[11][0] = 0.5*x[0][0] + 0.5*x[1][0];
9551
coordinates[11][1] = 0.5*x[0][1] + 0.5*x[1][1];
9464
dof_coordinates[0][0] = vertex_coordinates[0];
9465
dof_coordinates[0][1] = vertex_coordinates[1];
9466
dof_coordinates[1][0] = vertex_coordinates[2];
9467
dof_coordinates[1][1] = vertex_coordinates[3];
9468
dof_coordinates[2][0] = vertex_coordinates[4];
9469
dof_coordinates[2][1] = vertex_coordinates[5];
9470
dof_coordinates[3][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
9471
dof_coordinates[3][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
9472
dof_coordinates[4][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
9473
dof_coordinates[4][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
9474
dof_coordinates[5][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
9475
dof_coordinates[5][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
9476
dof_coordinates[6][0] = vertex_coordinates[0];
9477
dof_coordinates[6][1] = vertex_coordinates[1];
9478
dof_coordinates[7][0] = vertex_coordinates[2];
9479
dof_coordinates[7][1] = vertex_coordinates[3];
9480
dof_coordinates[8][0] = vertex_coordinates[4];
9481
dof_coordinates[8][1] = vertex_coordinates[5];
9482
dof_coordinates[9][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
9483
dof_coordinates[9][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
9484
dof_coordinates[10][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
9485
dof_coordinates[10][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
9486
dof_coordinates[11][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
9487
dof_coordinates[11][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
9552
9488
}
9553
9489
9554
9490
/// Return the number of sub dofmaps (for a mixed element)
9555
virtual unsigned int num_sub_dofmaps() const
9491
virtual std::size_t num_sub_dofmaps() const
9556
9492
{
9557
9493
return 2;
9558
9494
}
9559
9495
9560
9496
/// Create a new dofmap for sub dofmap i (for a mixed element)
9561
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
9497
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
9562
9498
{
9563
9499
switch (i)
9564
9500
{
9590
9526
9591
9527
class stokes_dofmap_2: public ufc::dofmap
9592
9528
{
9593
private:
9594
9595
unsigned int _global_dimension;
9596
9529
public:
9597
9530
9598
9531
/// Constructor
9599
9532
stokes_dofmap_2() : ufc::dofmap()
9600
9533
{
9601
_global_dimension = 0;
9534
// Do nothing
9602
9535
}
9603
9536
9604
9537
/// Destructor
9610
9543
/// Return a string identifying the dofmap
9611
9544
virtual const char* signature() const
9612
9545
{
9613
return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)";
9546
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
9614
9547
}
9615
9548
9616
9549
/// Return true iff mesh entities of topological dimension d are needed
9617
virtual bool needs_mesh_entities(unsigned int d) const
9550
virtual bool needs_mesh_entities(std::size_t d) const
9618
9551
{
9619
9552
switch (d)
9620
9553
{
9638
9571
return false;
9639
9572
}
9640
9573
9641
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
9642
virtual bool init_mesh(const ufc::mesh& m)
9643
{
9644
_global_dimension = m.num_entities[0];
9645
return false;
9646
}
9647
9648
/// Initialize dofmap for given cell
9649
virtual void init_cell(const ufc::mesh& m,
9650
const ufc::cell& c)
9651
{
9652
// Do nothing
9653
}
9654
9655
/// Finish initialization of dofmap for cells
9656
virtual void init_cell_finalize()
9657
{
9658
// Do nothing
9659
}
9660
9661
9574
/// Return the topological dimension of the associated cell shape
9662
virtual unsigned int topological_dimension() const
9575
virtual std::size_t topological_dimension() const
9663
9576
{
9664
9577
return 2;
9665
9578
}
9666
9579
9667
9580
/// Return the geometric dimension of the associated cell shape
9668
virtual unsigned int geometric_dimension() const
9581
virtual std::size_t geometric_dimension() const
9669
9582
{
9670
9583
return 2;
9671
9584
}
9672
9585
9673
9586
/// Return the dimension of the global finite element function space
9674
virtual unsigned int global_dimension() const
9587
virtual std::size_t global_dimension(const std::vector<std::size_t>&
9588
num_global_entities) const
9675
9589
{
9676
return _global_dimension;
9590
return num_global_entities[0];
9677
9591
}
9678
9592
9679
9593
/// Return the dimension of the local finite element function space for a cell
9680
virtual unsigned int local_dimension(const ufc::cell& c) const
9594
virtual std::size_t local_dimension(const ufc::cell& c) const
9681
9595
{
9682
9596
return 3;
9683
9597
}
9684
9598
9685
9599
/// Return the maximum dimension of the local finite element function space
9686
virtual unsigned int max_local_dimension() const
9600
virtual std::size_t max_local_dimension() const
9687
9601
{
9688
9602
return 3;
9689
9603
}
9690
9604
9691
9605
/// Return the number of dofs on each cell facet
9692
virtual unsigned int num_facet_dofs() const
9606
virtual std::size_t num_facet_dofs() const
9693
9607
{
9694
9608
return 2;
9695
9609
}
9696
9610
9697
9611
/// Return the number of dofs associated with each cell entity of dimension d
9698
virtual unsigned int num_entity_dofs(unsigned int d) const
9612
virtual std::size_t num_entity_dofs(std::size_t d) const
9699
9613
{
9700
9614
switch (d)
9701
9615
{
9720
9634
}
9721
9635
9722
9636
/// Tabulate the local-to-global mapping of dofs on a cell
9723
virtual void tabulate_dofs(unsigned int* dofs,
9724
const ufc::mesh& m,
9637
virtual void tabulate_dofs(std::size_t* dofs,
9638
const std::vector<std::size_t>& num_global_entities,
9725
9639
const ufc::cell& c) const
9726
9640
{
9727
9641
dofs[0] = c.entity_indices[0][0];
9730
9644
}
9731
9645
9732
9646
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
9733
virtual void tabulate_facet_dofs(unsigned int* dofs,
9734
unsigned int facet) const
9647
virtual void tabulate_facet_dofs(std::size_t* dofs,
9648
std::size_t facet) const
9735
9649
{
9736
9650
switch (facet)
9737
9651
{
9758
9672
}
9759
9673
9760
9674
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
9761
virtual void tabulate_entity_dofs(unsigned int* dofs,
9762
unsigned int d, unsigned int i) const
9675
virtual void tabulate_entity_dofs(std::size_t* dofs,
9676
std::size_t d, std::size_t i) const
9763
9677
{
9764
9678
if (d > 2)
9765
9679
{
9811
9725
}
9812
9726
9813
9727
/// Tabulate the coordinates of all dofs on a cell
9814
virtual void tabulate_coordinates(double** coordinates,
9815
const ufc::cell& c) const
9728
virtual void tabulate_coordinates(double** dof_coordinates,
9729
const double* vertex_coordinates) const
9816
9730
{
9817
const double * const * x = c.coordinates;
9818
9819
coordinates[0][0] = x[0][0];
9820
coordinates[0][1] = x[0][1];
9821
coordinates[1][0] = x[1][0];
9822
coordinates[1][1] = x[1][1];
9823
coordinates[2][0] = x[2][0];
9824
coordinates[2][1] = x[2][1];
9731
dof_coordinates[0][0] = vertex_coordinates[0];
9732
dof_coordinates[0][1] = vertex_coordinates[1];
9733
dof_coordinates[1][0] = vertex_coordinates[2];
9734
dof_coordinates[1][1] = vertex_coordinates[3];
9735
dof_coordinates[2][0] = vertex_coordinates[4];
9736
dof_coordinates[2][1] = vertex_coordinates[5];
9825
9737
}
9826
9738
9827
9739
/// Return the number of sub dofmaps (for a mixed element)
9828
virtual unsigned int num_sub_dofmaps() const
9740
virtual std::size_t num_sub_dofmaps() const
9829
9741
{
9830
9742
return 0;
9831
9743
}
9832
9744
9833
9745
/// Create a new dofmap for sub dofmap i (for a mixed element)
9834
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
9746
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
9835
9747
{
9836
9748
return 0;
9837
9749
}
9849
9761
9850
9762
class stokes_dofmap_3: public ufc::dofmap
9851
9763
{
9852
private:
9853
9854
unsigned int _global_dimension;
9855
9764
public:
9856
9765
9857
9766
/// Constructor
9858
9767
stokes_dofmap_3() : ufc::dofmap()
9859
9768
{
9860
_global_dimension = 0;
9769
// Do nothing
9861
9770
}
9862
9771
9863
9772
/// Destructor
9869
9778
/// Return a string identifying the dofmap
9870
9779
virtual const char* signature() const
9871
9780
{
9872
return "FFC dofmap for MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) })";
9781
return "FFC dofmap for MixedElement(*[VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (3,) })";
9873
9782
}
9874
9783
9875
9784
/// Return true iff mesh entities of topological dimension d are needed
9876
virtual bool needs_mesh_entities(unsigned int d) const
9785
virtual bool needs_mesh_entities(std::size_t d) const
9877
9786
{
9878
9787
switch (d)
9879
9788
{
9897
9806
return false;
9898
9807
}
9899
9808
9900
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
9901
virtual bool init_mesh(const ufc::mesh& m)
9902
{
9903
_global_dimension = 3*m.num_entities[0] + 2*m.num_entities[1];
9904
return false;
9905
}
9906
9907
/// Initialize dofmap for given cell
9908
virtual void init_cell(const ufc::mesh& m,
9909
const ufc::cell& c)
9910
{
9911
// Do nothing
9912
}
9913
9914
/// Finish initialization of dofmap for cells
9915
virtual void init_cell_finalize()
9916
{
9917
// Do nothing
9918
}
9919
9920
9809
/// Return the topological dimension of the associated cell shape
9921
virtual unsigned int topological_dimension() const
9810
virtual std::size_t topological_dimension() const
9922
9811
{
9923
9812
return 2;
9924
9813
}
9925
9814
9926
9815
/// Return the geometric dimension of the associated cell shape
9927
virtual unsigned int geometric_dimension() const
9816
virtual std::size_t geometric_dimension() const
9928
9817
{
9929
9818
return 2;
9930
9819
}
9931
9820
9932
9821
/// Return the dimension of the global finite element function space
9933
virtual unsigned int global_dimension() const
9822
virtual std::size_t global_dimension(const std::vector<std::size_t>&
9823
num_global_entities) const
9934
9824
{
9935
return _global_dimension;
9825
return 3*num_global_entities[0] + 2*num_global_entities[1];
9936
9826
}
9937
9827
9938
9828
/// Return the dimension of the local finite element function space for a cell
9939
virtual unsigned int local_dimension(const ufc::cell& c) const
9829
virtual std::size_t local_dimension(const ufc::cell& c) const
9940
9830
{
9941
9831
return 15;
9942
9832
}
9943
9833
9944
9834
/// Return the maximum dimension of the local finite element function space
9945
virtual unsigned int max_local_dimension() const
9835
virtual std::size_t max_local_dimension() const
9946
9836
{
9947
9837
return 15;
9948
9838
}
9949
9839
9950
9840
/// Return the number of dofs on each cell facet
9951
virtual unsigned int num_facet_dofs() const
9841
virtual std::size_t num_facet_dofs() const
9952
9842
{
9953
9843
return 8;
9954
9844
}
9955
9845
9956
9846
/// Return the number of dofs associated with each cell entity of dimension d
9957
virtual unsigned int num_entity_dofs(unsigned int d) const
9847
virtual std::size_t num_entity_dofs(std::size_t d) const
9958
9848
{
9959
9849
switch (d)
9960
9850
{
9979
9869
}
9980
9870
9981
9871
/// Tabulate the local-to-global mapping of dofs on a cell
9982
virtual void tabulate_dofs(unsigned int* dofs,
9983
const ufc::mesh& m,
9872
virtual void tabulate_dofs(std::size_t* dofs,
9873
const std::vector<std::size_t>& num_global_entities,
9984
9874
const ufc::cell& c) const
9985
9875
{
9986
9876
unsigned int offset = 0;
9987
9877
dofs[0] = offset + c.entity_indices[0][0];
9988
9878
dofs[1] = offset + c.entity_indices[0][1];
9989
9879
dofs[2] = offset + c.entity_indices[0][2];
9990
offset += m.num_entities[0];
9880
offset += num_global_entities[0];
9991
9881
dofs[3] = offset + c.entity_indices[1][0];
9992
9882
dofs[4] = offset + c.entity_indices[1][1];
9993
9883
dofs[5] = offset + c.entity_indices[1][2];
9994
offset += m.num_entities[1];
9884
offset += num_global_entities[1];
9995
9885
dofs[6] = offset + c.entity_indices[0][0];
9996
9886
dofs[7] = offset + c.entity_indices[0][1];
9997
9887
dofs[8] = offset + c.entity_indices[0][2];
9998
offset += m.num_entities[0];
9888
offset += num_global_entities[0];
9999
9889
dofs[9] = offset + c.entity_indices[1][0];
10000
9890
dofs[10] = offset + c.entity_indices[1][1];
10001
9891
dofs[11] = offset + c.entity_indices[1][2];
10002
offset += m.num_entities[1];
9892
offset += num_global_entities[1];
10003
9893
dofs[12] = offset + c.entity_indices[0][0];
10004
9894
dofs[13] = offset + c.entity_indices[0][1];
10005
9895
dofs[14] = offset + c.entity_indices[0][2];
10006
offset += m.num_entities[0];
9896
offset += num_global_entities[0];
10007
9897
}
10008
9898
10009
9899
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
10010
virtual void tabulate_facet_dofs(unsigned int* dofs,
10011
unsigned int facet) const
9900
virtual void tabulate_facet_dofs(std::size_t* dofs,
9901
std::size_t facet) const
10012
9902
{
10013
9903
switch (facet)
10014
9904
{
10053
9943
}
10054
9944
10055
9945
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
10056
virtual void tabulate_entity_dofs(unsigned int* dofs,
10057
unsigned int d, unsigned int i) const
9946
virtual void tabulate_entity_dofs(std::size_t* dofs,
9947
std::size_t d, std::size_t i) const
10058
9948
{
10059
9949
if (d > 2)
10060
9950
{
10138
10028
}
10139
10029
10140
10030
/// Tabulate the coordinates of all dofs on a cell
10141
virtual void tabulate_coordinates(double** coordinates,
10142
const ufc::cell& c) const
10031
virtual void tabulate_coordinates(double** dof_coordinates,
10032
const double* vertex_coordinates) const
10143
10033
{
10144
const double * const * x = c.coordinates;
10145
10146
coordinates[0][0] = x[0][0];
10147
coordinates[0][1] = x[0][1];
10148
coordinates[1][0] = x[1][0];
10149
coordinates[1][1] = x[1][1];
10150
coordinates[2][0] = x[2][0];
10151
coordinates[2][1] = x[2][1];
10152
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
10153
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
10154
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
10155
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
10156
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
10157
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
10158
coordinates[6][0] = x[0][0];
10159
coordinates[6][1] = x[0][1];
10160
coordinates[7][0] = x[1][0];
10161
coordinates[7][1] = x[1][1];
10162
coordinates[8][0] = x[2][0];
10163
coordinates[8][1] = x[2][1];
10164
coordinates[9][0] = 0.5*x[1][0] + 0.5*x[2][0];
10165
coordinates[9][1] = 0.5*x[1][1] + 0.5*x[2][1];
10166
coordinates[10][0] = 0.5*x[0][0] + 0.5*x[2][0];
10167
coordinates[10][1] = 0.5*x[0][1] + 0.5*x[2][1];
10168
coordinates[11][0] = 0.5*x[0][0] + 0.5*x[1][0];
10169
coordinates[11][1] = 0.5*x[0][1] + 0.5*x[1][1];
10170
coordinates[12][0] = x[0][0];
10171
coordinates[12][1] = x[0][1];
10172
coordinates[13][0] = x[1][0];
10173
coordinates[13][1] = x[1][1];
10174
coordinates[14][0] = x[2][0];
10175
coordinates[14][1] = x[2][1];
10034
dof_coordinates[0][0] = vertex_coordinates[0];
10035
dof_coordinates[0][1] = vertex_coordinates[1];
10036
dof_coordinates[1][0] = vertex_coordinates[2];
10037
dof_coordinates[1][1] = vertex_coordinates[3];
10038
dof_coordinates[2][0] = vertex_coordinates[4];
10039
dof_coordinates[2][1] = vertex_coordinates[5];
10040
dof_coordinates[3][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
10041
dof_coordinates[3][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
10042
dof_coordinates[4][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
10043
dof_coordinates[4][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
10044
dof_coordinates[5][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
10045
dof_coordinates[5][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
10046
dof_coordinates[6][0] = vertex_coordinates[0];
10047
dof_coordinates[6][1] = vertex_coordinates[1];
10048
dof_coordinates[7][0] = vertex_coordinates[2];
10049
dof_coordinates[7][1] = vertex_coordinates[3];
10050
dof_coordinates[8][0] = vertex_coordinates[4];
10051
dof_coordinates[8][1] = vertex_coordinates[5];
10052
dof_coordinates[9][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
10053
dof_coordinates[9][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
10054
dof_coordinates[10][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
10055
dof_coordinates[10][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
10056
dof_coordinates[11][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
10057
dof_coordinates[11][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
10058
dof_coordinates[12][0] = vertex_coordinates[0];
10059
dof_coordinates[12][1] = vertex_coordinates[1];
10060
dof_coordinates[13][0] = vertex_coordinates[2];
10061
dof_coordinates[13][1] = vertex_coordinates[3];
10062
dof_coordinates[14][0] = vertex_coordinates[4];
10063
dof_coordinates[14][1] = vertex_coordinates[5];
10176
10064
}
10177
10065
10178
10066
/// Return the number of sub dofmaps (for a mixed element)
10179
virtual unsigned int num_sub_dofmaps() const
10067
virtual std::size_t num_sub_dofmaps() const
10180
10068
{
10181
10069
return 2;
10182
10070
}
10183
10071
10184
10072
/// Create a new dofmap for sub dofmap i (for a mixed element)
10185
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
10073
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
10186
10074
{
10187
10075
switch (i)
10188
10076
{
10213
10101
/// tensor corresponding to the local contribution to a form from
10214
10102
/// the integral over a cell.
10215
10103
10216
class stokes_cell_integral_0_0: public ufc::cell_integral
10104
class stokes_cell_integral_0_otherwise: public ufc::cell_integral
10217
10105
{
10218
10106
public:
10219
10107
10220
10108
/// Constructor
10221
stokes_cell_integral_0_0() : ufc::cell_integral()
10109
stokes_cell_integral_0_otherwise() : ufc::cell_integral()
10222
10110
{
10223
10111
// Do nothing
10224
10112
}
10225
10113
10226
10114
/// Destructor
10227
virtual ~stokes_cell_integral_0_0()
10115
virtual ~stokes_cell_integral_0_otherwise()
10228
10116
{
10229
10117
// Do nothing
10230
10118
}
10232
10120
/// Tabulate the tensor for the contribution from a local cell
10233
10121
virtual void tabulate_tensor(double* A,
10234
10122
const double * const * w,
10235
const ufc::cell& c) const
10123
const double* vertex_coordinates,
10124
int cell_orientation) const
10236
10125
{
10237
// Extract vertex coordinates
10238
const double * const * x = c.coordinates;
10239
10240
// Compute Jacobian of affine map from reference cell
10241
const double J_00 = x[1][0] - x[0][0];
10242
const double J_01 = x[2][0] - x[0][0];
10243
const double J_10 = x[1][1] - x[0][1];
10244
const double J_11 = x[2][1] - x[0][1];
10245
10246
// Compute determinant of Jacobian
10247
double detJ = J_00*J_11 - J_01*J_10;
10248
10249
// Compute inverse of Jacobian
10250
const double K_00 = J_11 / detJ;
10251
const double K_01 = -J_01 / detJ;
10252
const double K_10 = -J_10 / detJ;
10253
const double K_11 = J_00 / detJ;
10126
// Compute Jacobian
10127
double J[4];
10128
compute_jacobian_triangle_2d(J, vertex_coordinates);
10129
10130
// Compute Jacobian inverse and determinant
10131
double K[4];
10132
double detJ;
10133
compute_jacobian_inverse_triangle_2d(K, detJ, J);
10254
10134
10255
10135
// Set scale factor
10256
10136
const double det = std::abs(detJ);
10257
10137
10258
// Cell Volume.
10259
10260
// Compute circumradius, assuming triangle is embedded in 2D.
10261
10262
10263
// Facet Area.
10138
// Cell volume
10139
10140
// Compute circumradius of triangle in 2D
10141
10142
10143
// Facet area
10264
10144
10265
10145
// Array of quadrature weights.
10266
10146
static const double W6[6] = {0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333};
10313
10193
}// end loop over 'r'
10314
10194
// Number of operations to compute geometry constants: 20.
10315
10195
double G[11];
10316
G[0] = K_10*det;
10317
G[1] = K_00*det;
10318
G[2] = K_11*det;
10319
G[3] = K_01*det;
10320
G[4] = - K_10*det;
10321
G[5] = - K_00*det;
10322
G[6] = - K_11*det;
10323
G[7] = - K_01*det;
10324
G[8] = det*(K_10*K_10 + K_11*K_11);
10325
G[9] = det*(K_00*K_10 + K_01*K_11);
10326
G[10] = det*(K_00*K_00 + K_01*K_01);
10196
G[0] = K[2]*det;
10197
G[1] = K[0]*det;
10198
G[2] = K[3]*det;
10199
G[3] = K[1]*det;
10200
G[4] = - K[2]*det;
10201
G[5] = - K[0]*det;
10202
G[6] = - K[3]*det;
10203
G[7] = - K[1]*det;
10204
G[8] = det*(K[2]*K[2] + K[3]*K[3]);
10205
G[9] = det*(K[0]*K[2] + K[1]*K[3]);
10206
G[10] = det*(K[0]*K[0] + K[1]*K[1]);
10327
10207
10328
10208
// Compute element tensor using UFL quadrature representation
10329
10209
// Optimisations: ('eliminate zeros', True), ('ignore ones', True), ('ignore zero tables', True), ('optimisation', 'simplify_expressions'), ('remove zero terms', True)
10427
10307
}// end loop over 'ip'
10428
10308
}
10429
10309
10430
/// Tabulate the tensor for the contribution from a local cell
10431
/// using the specified reference cell quadrature points/weights
10432
virtual void tabulate_tensor(double* A,
10433
const double * const * w,
10434
const ufc::cell& c,
10435
unsigned int num_quadrature_points,
10436
const double * const * quadrature_points,
10437
const double* quadrature_weights) const
10438
{
10439
std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
10440
}
10441
10442
10310
};
10443
10311
10444
10312
/// This class defines the interface for the tabulation of the cell
10445
10313
/// tensor corresponding to the local contribution to a form from
10446
10314
/// the integral over a cell.
10447
10315
10448
class stokes_cell_integral_1_0: public ufc::cell_integral
10316
class stokes_cell_integral_1_otherwise: public ufc::cell_integral
10449
10317
{
10450
10318
public:
10451
10319
10452
10320
/// Constructor
10453
stokes_cell_integral_1_0() : ufc::cell_integral()
10321
stokes_cell_integral_1_otherwise() : ufc::cell_integral()
10454
10322
{
10455
10323
// Do nothing
10456
10324
}
10457
10325
10458
10326
/// Destructor
10459
virtual ~stokes_cell_integral_1_0()
10327
virtual ~stokes_cell_integral_1_otherwise()
10460
10328
{
10461
10329
// Do nothing
10462
10330
}
10464
10332
/// Tabulate the tensor for the contribution from a local cell
10465
10333
virtual void tabulate_tensor(double* A,
10466
10334
const double * const * w,
10467
const ufc::cell& c) const
10335
const double* vertex_coordinates,
10336
int cell_orientation) const
10468
10337
{
10469
// Extract vertex coordinates
10470
const double * const * x = c.coordinates;
10471
10472
// Compute Jacobian of affine map from reference cell
10473
const double J_00 = x[1][0] - x[0][0];
10474
const double J_01 = x[2][0] - x[0][0];
10475
const double J_10 = x[1][1] - x[0][1];
10476
const double J_11 = x[2][1] - x[0][1];
10477
10478
// Compute determinant of Jacobian
10479
double detJ = J_00*J_11 - J_01*J_10;
10480
10481
// Compute inverse of Jacobian
10338
// Compute Jacobian
10339
double J[4];
10340
compute_jacobian_triangle_2d(J, vertex_coordinates);
10341
10342
// Compute Jacobian inverse and determinant
10343
double K[4];
10344
double detJ;
10345
compute_jacobian_inverse_triangle_2d(K, detJ, J);
10482
10346
10483
10347
// Set scale factor
10484
10348
const double det = std::abs(detJ);
10485
10349
10486
// Cell Volume.
10487
10488
// Compute circumradius, assuming triangle is embedded in 2D.
10489
10490
10491
// Facet Area.
10350
// Cell volume
10351
10352
// Compute circumradius of triangle in 2D
10353
10354
10355
// Facet area
10492
10356
10493
10357
// Array of quadrature weights.
10494
10358
static const double W6[6] = {0.054975872, 0.054975872, 0.054975872, 0.11169079, 0.11169079, 0.11169079};
10560
10424
}// end loop over 'ip'
10561
10425
}
10562
10426
10563
/// Tabulate the tensor for the contribution from a local cell
10564
/// using the specified reference cell quadrature points/weights
10565
virtual void tabulate_tensor(double* A,
10566
const double * const * w,
10567
const ufc::cell& c,
10568
unsigned int num_quadrature_points,
10569
const double * const * quadrature_points,
10570
const double* quadrature_weights) const
10571
{
10572
std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
10573
}
10574
10575
10427
};
10576
10428
10577
10429
/// This class defines the interface for the assembly of the global
10608
10460
/// Return a string identifying the form
10609
10461
virtual const char* signature() const
10610
10462
{
10611
return "Form([Integral(Sum(Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(2),), {}))), Sum(IndexSum(IndexSum(Product(Indexed(ComponentTensor(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2), Index(1)), {Index(2): 2, Index(1): 2})), MultiIndex((Index(3), Index(4)), {Index(4): 2, Index(3): 2})), Indexed(ComponentTensor(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(5),), {Index(5): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(5),), {Index(5): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(6),), {Index(6): 2})), MultiIndex((Index(6), Index(5)), {Index(5): 2, Index(6): 2})), MultiIndex((Index(3), Index(4)), {Index(4): 2, Index(3): 2}))), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((Index(4),), {Index(4): 2})), Product(IntValue(-1, (), (), {}), Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(7),), {Index(7): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(7),), {Index(7): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(7),), {Index(7): 2})), MultiIndex((Index(7),), {Index(7): 2})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((FixedIndex(2),), {})))))), Measure('cell', 0, None))])";
10463
return "1cb2159c037a79a64f83e9f9cf5fe4a59801e9e8993b5164e21bfbc99d8fca959d238af2f714ec13c9ecd37f4b2396497ae99792db1c82d32a28d1d5b1f94450";
10612
10464
}
10613
10465
10614
10466
/// Return the rank of the global tensor (r)
10615
virtual unsigned int rank() const
10467
virtual std::size_t rank() const
10616
10468
{
10617
10469
return 2;
10618
10470
}
10619
10471
10620
10472
/// Return the number of coefficients (n)
10621
virtual unsigned int num_coefficients() const
10473
virtual std::size_t num_coefficients() const
10622
10474
{
10623
10475
return 0;
10624
10476
}
10625
10477
10626
10478
/// Return the number of cell domains
10627
virtual unsigned int num_cell_domains() const
10479
virtual std::size_t num_cell_domains() const
10628
10480
{
10629
return 1;
10481
return 0;
10630
10482
}
10631
10483
10632
10484
/// Return the number of exterior facet domains
10633
virtual unsigned int num_exterior_facet_domains() const
10485
virtual std::size_t num_exterior_facet_domains() const
10634
10486
{
10635
10487
return 0;
10636
10488
}
10637
10489
10638
10490
/// Return the number of interior facet domains
10639
virtual unsigned int num_interior_facet_domains() const
10640
{
10641
return 0;
10491
virtual std::size_t num_interior_facet_domains() const
10492
{
10493
return 0;
10494
}
10495
10496
/// Return the number of point domains
10497
virtual std::size_t num_point_domains() const
10498
{
10499
return 0;
10500
}
10501
10502
/// Return whether the form has any cell integrals
10503
virtual bool has_cell_integrals() const
10504
{
10505
return true;
10506
}
10507
10508
/// Return whether the form has any exterior facet integrals
10509
virtual bool has_exterior_facet_integrals() const
10510
{
10511
return false;
10512
}
10513
10514
/// Return whether the form has any interior facet integrals
10515
virtual bool has_interior_facet_integrals() const
10516
{
10517
return false;
10518
}
10519
10520
/// Return whether the form has any point integrals
10521
virtual bool has_point_integrals() const
10522
{
10523
return false;
10642
10524
}
10643
10525
10644
10526
/// Create a new finite element for argument function i
10645
virtual ufc::finite_element* create_finite_element(unsigned int i) const
10527
virtual ufc::finite_element* create_finite_element(std::size_t i) const
10646
10528
{
10647
10529
switch (i)
10648
10530
{
10662
10544
}
10663
10545
10664
10546
/// Create a new dofmap for argument function i
10665
virtual ufc::dofmap* create_dofmap(unsigned int i) const
10547
virtual ufc::dofmap* create_dofmap(std::size_t i) const
10666
10548
{
10667
10549
switch (i)
10668
10550
{
10682
10564
}
10683
10565
10684
10566
/// Create a new cell integral on sub domain i
10685
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
10567
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
10686
10568
{
10687
switch (i)
10688
{
10689
case 0:
10690
{
10691
return new stokes_cell_integral_0_0();
10692
break;
10693
}
10694
}
10695
10696
10569
return 0;
10697
10570
}
10698
10571
10699
10572
/// Create a new exterior facet integral on sub domain i
10700
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
10573
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
10701
10574
{
10702
10575
return 0;
10703
10576
}
10704
10577
10705
10578
/// Create a new interior facet integral on sub domain i
10706
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
10579
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
10580
{
10581
return 0;
10582
}
10583
10584
/// Create a new point integral on sub domain i
10585
virtual ufc::point_integral* create_point_integral(std::size_t i) const
10586
{
10587
return 0;
10588
}
10589
10590
/// Create a new cell integral on everywhere else
10591
virtual ufc::cell_integral* create_default_cell_integral() const
10592
{
10593
return new stokes_cell_integral_0_otherwise();
10594
}
10595
10596
/// Create a new exterior facet integral on everywhere else
10597
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
10598
{
10599
return 0;
10600
}
10601
10602
/// Create a new interior facet integral on everywhere else
10603
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
10604
{
10605
return 0;
10606
}
10607
10608
/// Create a new point integral on everywhere else
10609
virtual ufc::point_integral* create_default_point_integral() const
10707
10610
{
10708
10611
return 0;
10709
10612
}
10744
10647
/// Return a string identifying the form
10745
10648
virtual const char* signature() const
10746
10649
{
10747
return "Form([Integral(IndexSum(Product(Indexed(Coefficient(VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), 0), MultiIndex((Index(0),), {Index(0): 2})), Indexed(ListTensor(Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(0),), {})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 2, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(0),), {Index(0): 2}))), MultiIndex((Index(0),), {Index(0): 2})), Measure('cell', 0, None))])";
10650
return "b9394062283c7f39e6ddf66053c2837a34c1885e9c7c6837808355de8e776ebb1ebfdaa298fb13f1bbb2e5eda76fb5ad2bf64c2f236366297a5af53cf6fd0b1b";
10748
10651
}
10749
10652
10750
10653
/// Return the rank of the global tensor (r)
10751
virtual unsigned int rank() const
10654
virtual std::size_t rank() const
10752
10655
{
10753
10656
return 1;
10754
10657
}
10755
10658
10756
10659
/// Return the number of coefficients (n)
10757
virtual unsigned int num_coefficients() const
10660
virtual std::size_t num_coefficients() const
10758
10661
{
10759
10662
return 1;
10760
10663
}
10761
10664
10762
10665
/// Return the number of cell domains
10763
virtual unsigned int num_cell_domains() const
10666
virtual std::size_t num_cell_domains() const
10764
10667
{
10765
return 1;
10668
return 0;
10766
10669
}
10767
10670
10768
10671
/// Return the number of exterior facet domains
10769
virtual unsigned int num_exterior_facet_domains() const
10672
virtual std::size_t num_exterior_facet_domains() const
10770
10673
{
10771
10674
return 0;
10772
10675
}
10773
10676
10774
10677
/// Return the number of interior facet domains
10775
virtual unsigned int num_interior_facet_domains() const
10776
{
10777
return 0;
10678
virtual std::size_t num_interior_facet_domains() const
10679
{
10680
return 0;
10681
}
10682
10683
/// Return the number of point domains
10684
virtual std::size_t num_point_domains() const
10685
{
10686
return 0;
10687
}
10688
10689
/// Return whether the form has any cell integrals
10690
virtual bool has_cell_integrals() const
10691
{
10692
return true;
10693
}
10694
10695
/// Return whether the form has any exterior facet integrals
10696
virtual bool has_exterior_facet_integrals() const
10697
{
10698
return false;
10699
}
10700
10701
/// Return whether the form has any interior facet integrals
10702
virtual bool has_interior_facet_integrals() const
10703
{
10704
return false;
10705
}
10706
10707
/// Return whether the form has any point integrals
10708
virtual bool has_point_integrals() const
10709
{
10710
return false;
10778
10711
}
10779
10712
10780
10713
/// Create a new finite element for argument function i
10781
virtual ufc::finite_element* create_finite_element(unsigned int i) const
10714
virtual ufc::finite_element* create_finite_element(std::size_t i) const
10782
10715
{
10783
10716
switch (i)
10784
10717
{
10798
10731
}
10799
10732
10800
10733
/// Create a new dofmap for argument function i
10801
virtual ufc::dofmap* create_dofmap(unsigned int i) const
10734
virtual ufc::dofmap* create_dofmap(std::size_t i) const
10802
10735
{
10803
10736
switch (i)
10804
10737
{
10818
10751
}
10819
10752
10820
10753
/// Create a new cell integral on sub domain i
10821
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
10754
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
10822
10755
{
10823
switch (i)
10824
{
10825
case 0:
10826
{
10827
return new stokes_cell_integral_1_0();
10828
break;
10829
}
10830
}
10831
10832
10756
return 0;
10833
10757
}
10834
10758
10835
10759
/// Create a new exterior facet integral on sub domain i
10836
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
10760
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
10837
10761
{
10838
10762
return 0;
10839
10763
}
10840
10764
10841
10765
/// Create a new interior facet integral on sub domain i
10842
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
10766
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
10767
{
10768
return 0;
10769
}
10770
10771
/// Create a new point integral on sub domain i
10772
virtual ufc::point_integral* create_point_integral(std::size_t i) const
10773
{
10774
return 0;
10775
}
10776
10777
/// Create a new cell integral on everywhere else
10778
virtual ufc::cell_integral* create_default_cell_integral() const
10779
{
10780
return new stokes_cell_integral_1_otherwise();
10781
}
10782
10783
/// Create a new exterior facet integral on everywhere else
10784
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
10785
{
10786
return 0;
10787
}
10788
10789
/// Create a new interior facet integral on everywhere else
10790
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
10791
{
10792
return 0;
10793
}
10794
10795
/// Create a new point integral on everywhere else
10796
virtual ufc::point_integral* create_default_point_integral() const
10843
10797
{
10844
10798
return 0;
10845
10799
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1