~ubuntu-branches/ubuntu/utopic/ffc/utopic
» Revision
17
: test/regression/references/r_quadrature/NavierStokes.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_quadrature/NavierStokes.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/NavierStokes.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('tetrahedron', Space(3)), 1, None)";
55
return "FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, 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 3;
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 3;
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 4;
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_02 = x[3][0] - x[0][0];
107
const double J_10 = x[1][1] - x[0][1];
108
const double J_11 = x[2][1] - x[0][1];
109
const double J_12 = x[3][1] - x[0][1];
110
const double J_20 = x[1][2] - x[0][2];
111
const double J_21 = x[2][2] - x[0][2];
112
const double J_22 = x[3][2] - x[0][2];
113
114
// Compute sub determinants
115
const double d_00 = J_11*J_22 - J_12*J_21;
116
const double d_01 = J_12*J_20 - J_10*J_22;
117
const double d_02 = J_10*J_21 - J_11*J_20;
118
const double d_10 = J_02*J_21 - J_01*J_22;
119
const double d_11 = J_00*J_22 - J_02*J_20;
120
const double d_12 = J_01*J_20 - J_00*J_21;
121
const double d_20 = J_01*J_12 - J_02*J_11;
122
const double d_21 = J_02*J_10 - J_00*J_12;
123
const double d_22 = J_00*J_11 - J_01*J_10;
124
125
// Compute determinant of Jacobian
126
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
127
128
// Compute inverse of Jacobian
101
// Compute Jacobian
102
double J[9];
103
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
104
105
// Compute Jacobian inverse and determinant
106
double K[9];
107
double detJ;
108
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
109
129
110
130
111
// Compute constants
131
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
132
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
133
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
112
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
113
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
114
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
115
116
// Compute subdeterminants
117
const double d_00 = J[4]*J[8] - J[5]*J[7];
118
const double d_01 = J[5]*J[6] - J[3]*J[8];
119
const double d_02 = J[3]*J[7] - J[4]*J[6];
120
const double d_10 = J[2]*J[7] - J[1]*J[8];
121
const double d_11 = J[0]*J[8] - J[2]*J[6];
122
const double d_12 = J[1]*J[6] - J[0]*J[7];
123
const double d_20 = J[1]*J[5] - J[2]*J[4];
124
const double d_21 = J[2]*J[3] - J[0]*J[5];
125
const double d_22 = J[0]*J[4] - J[1]*J[3];
134
126
135
127
// Get coordinates and map to the reference (FIAT) element
136
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
137
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
138
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
139
140
141
// Reset values.
128
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
129
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
130
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
131
132
133
// Reset values
142
134
*values = 0.0;
143
135
switch (i)
144
136
{
145
137
case 0:
146
138
{
147
139
148
// Array of basisvalues.
140
// Array of basisvalues
149
141
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
150
142
151
// Declare helper variables.
143
// Declare helper variables
152
144
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
153
145
154
// Compute basisvalues.
146
// Compute basisvalues
155
147
basisvalues[0] = 1.0;
156
148
basisvalues[1] = tmp0;
157
149
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
161
153
basisvalues[2] *= std::sqrt(2.5);
162
154
basisvalues[1] *= std::sqrt(7.5);
163
155
164
// Table(s) of coefficients.
156
// Table(s) of coefficients
165
157
static const double coefficients0[4] = \
166
158
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
167
159
168
// Compute value(s).
160
// Compute value(s)
169
161
for (unsigned int r = 0; r < 4; r++)
170
162
{
171
163
*values += coefficients0[r]*basisvalues[r];
175
167
case 1:
176
168
{
177
169
178
// Array of basisvalues.
170
// Array of basisvalues
179
171
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
180
172
181
// Declare helper variables.
173
// Declare helper variables
182
174
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
183
175
184
// Compute basisvalues.
176
// Compute basisvalues
185
177
basisvalues[0] = 1.0;
186
178
basisvalues[1] = tmp0;
187
179
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
191
183
basisvalues[2] *= std::sqrt(2.5);
192
184
basisvalues[1] *= std::sqrt(7.5);
193
185
194
// Table(s) of coefficients.
186
// Table(s) of coefficients
195
187
static const double coefficients0[4] = \
196
188
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
197
189
198
// Compute value(s).
190
// Compute value(s)
199
191
for (unsigned int r = 0; r < 4; r++)
200
192
{
201
193
*values += coefficients0[r]*basisvalues[r];
205
197
case 2:
206
198
{
207
199
208
// Array of basisvalues.
200
// Array of basisvalues
209
201
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
210
202
211
// Declare helper variables.
203
// Declare helper variables
212
204
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
213
205
214
// Compute basisvalues.
206
// Compute basisvalues
215
207
basisvalues[0] = 1.0;
216
208
basisvalues[1] = tmp0;
217
209
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
221
213
basisvalues[2] *= std::sqrt(2.5);
222
214
basisvalues[1] *= std::sqrt(7.5);
223
215
224
// Table(s) of coefficients.
216
// Table(s) of coefficients
225
217
static const double coefficients0[4] = \
226
218
{0.28867513, 0.0, 0.21081851, -0.074535599};
227
219
228
// Compute value(s).
220
// Compute value(s)
229
221
for (unsigned int r = 0; r < 4; r++)
230
222
{
231
223
*values += coefficients0[r]*basisvalues[r];
235
227
case 3:
236
228
{
237
229
238
// Array of basisvalues.
230
// Array of basisvalues
239
231
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
240
232
241
// Declare helper variables.
233
// Declare helper variables
242
234
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
243
235
244
// Compute basisvalues.
236
// Compute basisvalues
245
237
basisvalues[0] = 1.0;
246
238
basisvalues[1] = tmp0;
247
239
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
251
243
basisvalues[2] *= std::sqrt(2.5);
252
244
basisvalues[1] *= std::sqrt(7.5);
253
245
254
// Table(s) of coefficients.
246
// Table(s) of coefficients
255
247
static const double coefficients0[4] = \
256
248
{0.28867513, 0.0, 0.0, 0.2236068};
257
249
258
// Compute value(s).
250
// Compute value(s)
259
251
for (unsigned int r = 0; r < 4; r++)
260
252
{
261
253
*values += coefficients0[r]*basisvalues[r];
266
258
267
259
}
268
260
269
/// Evaluate all basis functions at given point in cell
261
/// Evaluate all basis functions at given point x in cell
270
262
virtual void evaluate_basis_all(double* values,
271
const double* coordinates,
272
const ufc::cell& c) const
263
const double* x,
264
const double* vertex_coordinates,
265
int cell_orientation) const
273
266
{
274
267
// Helper variable to hold values of a single dof.
275
268
double dof_values = 0.0;
276
269
277
// Loop dofs and call evaluate_basis.
270
// Loop dofs and call evaluate_basis
278
271
for (unsigned int r = 0; r < 4; r++)
279
272
{
280
evaluate_basis(r, &dof_values, coordinates, c);
273
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
281
274
values[r] = dof_values;
282
275
}// end loop over 'r'
283
276
}
284
277
285
/// Evaluate order n derivatives of basis function i at given point in cell
286
virtual void evaluate_basis_derivatives(unsigned int i,
287
unsigned int n,
278
/// Evaluate order n derivatives of basis function i at given point x in cell
279
virtual void evaluate_basis_derivatives(std::size_t i,
280
std::size_t n,
288
281
double* values,
289
const double* coordinates,
290
const ufc::cell& c) const
282
const double* x,
283
const double* vertex_coordinates,
284
int cell_orientation) const
291
285
{
292
// Extract vertex coordinates
293
const double * const * x = c.coordinates;
294
295
// Compute Jacobian of affine map from reference cell
296
const double J_00 = x[1][0] - x[0][0];
297
const double J_01 = x[2][0] - x[0][0];
298
const double J_02 = x[3][0] - x[0][0];
299
const double J_10 = x[1][1] - x[0][1];
300
const double J_11 = x[2][1] - x[0][1];
301
const double J_12 = x[3][1] - x[0][1];
302
const double J_20 = x[1][2] - x[0][2];
303
const double J_21 = x[2][2] - x[0][2];
304
const double J_22 = x[3][2] - x[0][2];
305
306
// Compute sub determinants
307
const double d_00 = J_11*J_22 - J_12*J_21;
308
const double d_01 = J_12*J_20 - J_10*J_22;
309
const double d_02 = J_10*J_21 - J_11*J_20;
310
const double d_10 = J_02*J_21 - J_01*J_22;
311
const double d_11 = J_00*J_22 - J_02*J_20;
312
const double d_12 = J_01*J_20 - J_00*J_21;
313
const double d_20 = J_01*J_12 - J_02*J_11;
314
const double d_21 = J_02*J_10 - J_00*J_12;
315
const double d_22 = J_00*J_11 - J_01*J_10;
316
317
// Compute determinant of Jacobian
318
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
319
320
// Compute inverse of Jacobian
321
const double K_00 = d_00 / detJ;
322
const double K_01 = d_10 / detJ;
323
const double K_02 = d_20 / detJ;
324
const double K_10 = d_01 / detJ;
325
const double K_11 = d_11 / detJ;
326
const double K_12 = d_21 / detJ;
327
const double K_20 = d_02 / detJ;
328
const double K_21 = d_12 / detJ;
329
const double K_22 = d_22 / detJ;
286
// Compute Jacobian
287
double J[9];
288
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
289
290
// Compute Jacobian inverse and determinant
291
double K[9];
292
double detJ;
293
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
294
330
295
331
296
// Compute constants
332
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
333
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
334
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
297
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
298
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
299
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
300
301
// Compute subdeterminants
302
const double d_00 = J[4]*J[8] - J[5]*J[7];
303
const double d_01 = J[5]*J[6] - J[3]*J[8];
304
const double d_02 = J[3]*J[7] - J[4]*J[6];
305
const double d_10 = J[2]*J[7] - J[1]*J[8];
306
const double d_11 = J[0]*J[8] - J[2]*J[6];
307
const double d_12 = J[1]*J[6] - J[0]*J[7];
308
const double d_20 = J[1]*J[5] - J[2]*J[4];
309
const double d_21 = J[2]*J[3] - J[0]*J[5];
310
const double d_22 = J[0]*J[4] - J[1]*J[3];
335
311
336
312
// Get coordinates and map to the reference (FIAT) element
337
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
338
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
339
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
313
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
314
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
315
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
340
316
341
317
342
318
// Compute number of derivatives.
374
350
}
375
351
376
352
// Compute inverse of Jacobian
377
const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}};
353
const double Jinv[3][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}, {K[6], K[7], K[8]}};
378
354
379
355
// Declare transformation matrix
380
356
// Declare pointer to two dimensional array and initialise
408
384
case 0:
409
385
{
410
386
411
// Array of basisvalues.
387
// Array of basisvalues
412
388
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
413
389
414
// Declare helper variables.
390
// Declare helper variables
415
391
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
416
392
417
// Compute basisvalues.
393
// Compute basisvalues
418
394
basisvalues[0] = 1.0;
419
395
basisvalues[1] = tmp0;
420
396
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
424
400
basisvalues[2] *= std::sqrt(2.5);
425
401
basisvalues[1] *= std::sqrt(7.5);
426
402
427
// Table(s) of coefficients.
403
// Table(s) of coefficients
428
404
static const double coefficients0[4] = \
429
405
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
430
406
580
556
case 1:
581
557
{
582
558
583
// Array of basisvalues.
559
// Array of basisvalues
584
560
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
585
561
586
// Declare helper variables.
562
// Declare helper variables
587
563
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
588
564
589
// Compute basisvalues.
565
// Compute basisvalues
590
566
basisvalues[0] = 1.0;
591
567
basisvalues[1] = tmp0;
592
568
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
596
572
basisvalues[2] *= std::sqrt(2.5);
597
573
basisvalues[1] *= std::sqrt(7.5);
598
574
599
// Table(s) of coefficients.
575
// Table(s) of coefficients
600
576
static const double coefficients0[4] = \
601
577
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
602
578
752
728
case 2:
753
729
{
754
730
755
// Array of basisvalues.
731
// Array of basisvalues
756
732
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
757
733
758
// Declare helper variables.
734
// Declare helper variables
759
735
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
760
736
761
// Compute basisvalues.
737
// Compute basisvalues
762
738
basisvalues[0] = 1.0;
763
739
basisvalues[1] = tmp0;
764
740
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
768
744
basisvalues[2] *= std::sqrt(2.5);
769
745
basisvalues[1] *= std::sqrt(7.5);
770
746
771
// Table(s) of coefficients.
747
// Table(s) of coefficients
772
748
static const double coefficients0[4] = \
773
749
{0.28867513, 0.0, 0.21081851, -0.074535599};
774
750
924
900
case 3:
925
901
{
926
902
927
// Array of basisvalues.
903
// Array of basisvalues
928
904
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
929
905
930
// Declare helper variables.
906
// Declare helper variables
931
907
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
932
908
933
// Compute basisvalues.
909
// Compute basisvalues
934
910
basisvalues[0] = 1.0;
935
911
basisvalues[1] = tmp0;
936
912
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
940
916
basisvalues[2] *= std::sqrt(2.5);
941
917
basisvalues[1] *= std::sqrt(7.5);
942
918
943
// Table(s) of coefficients.
919
// Table(s) of coefficients
944
920
static const double coefficients0[4] = \
945
921
{0.28867513, 0.0, 0.0, 0.2236068};
946
922
1097
1073
1098
1074
}
1099
1075
1100
/// Evaluate order n derivatives of all basis functions at given point in cell
1101
virtual void evaluate_basis_derivatives_all(unsigned int n,
1076
/// Evaluate order n derivatives of all basis functions at given point x in cell
1077
virtual void evaluate_basis_derivatives_all(std::size_t n,
1102
1078
double* values,
1103
const double* coordinates,
1104
const ufc::cell& c) const
1079
const double* x,
1080
const double* vertex_coordinates,
1081
int cell_orientation) const
1105
1082
{
1106
1083
// Compute number of derivatives.
1107
1084
unsigned int num_derivatives = 1;
1120
1097
// Loop dofs and call evaluate_basis_derivatives.
1121
1098
for (unsigned int r = 0; r < 4; r++)
1122
1099
{
1123
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
1100
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1124
1101
for (unsigned int s = 0; s < num_derivatives; s++)
1125
1102
{
1126
1103
values[r*num_derivatives + s] = dof_values[s];
1132
1109
}
1133
1110
1134
1111
/// Evaluate linear functional for dof i on the function f
1135
virtual double evaluate_dof(unsigned int i,
1112
virtual double evaluate_dof(std::size_t i,
1136
1113
const ufc::function& f,
1114
const double* vertex_coordinates,
1115
int cell_orientation,
1137
1116
const ufc::cell& c) const
1138
1117
{
1139
// Declare variables for result of evaluation.
1118
// Declare variables for result of evaluation
1140
1119
double vals[1];
1141
1120
1142
// Declare variable for physical coordinates.
1121
// Declare variable for physical coordinates
1143
1122
double y[3];
1144
const double * const * x = c.coordinates;
1145
1123
switch (i)
1146
1124
{
1147
1125
case 0:
1148
1126
{
1149
y[0] = x[0][0];
1150
y[1] = x[0][1];
1151
y[2] = x[0][2];
1127
y[0] = vertex_coordinates[0];
1128
y[1] = vertex_coordinates[1];
1129
y[2] = vertex_coordinates[2];
1152
1130
f.evaluate(vals, y, c);
1153
1131
return vals[0];
1154
1132
break;
1155
1133
}
1156
1134
case 1:
1157
1135
{
1158
y[0] = x[1][0];
1159
y[1] = x[1][1];
1160
y[2] = x[1][2];
1136
y[0] = vertex_coordinates[3];
1137
y[1] = vertex_coordinates[4];
1138
y[2] = vertex_coordinates[5];
1161
1139
f.evaluate(vals, y, c);
1162
1140
return vals[0];
1163
1141
break;
1164
1142
}
1165
1143
case 2:
1166
1144
{
1167
y[0] = x[2][0];
1168
y[1] = x[2][1];
1169
y[2] = x[2][2];
1145
y[0] = vertex_coordinates[6];
1146
y[1] = vertex_coordinates[7];
1147
y[2] = vertex_coordinates[8];
1170
1148
f.evaluate(vals, y, c);
1171
1149
return vals[0];
1172
1150
break;
1173
1151
}
1174
1152
case 3:
1175
1153
{
1176
y[0] = x[3][0];
1177
y[1] = x[3][1];
1178
y[2] = x[3][2];
1154
y[0] = vertex_coordinates[9];
1155
y[1] = vertex_coordinates[10];
1156
y[2] = vertex_coordinates[11];
1179
1157
f.evaluate(vals, y, c);
1180
1158
return vals[0];
1181
1159
break;
1188
1166
/// Evaluate linear functionals for all dofs on the function f
1189
1167
virtual void evaluate_dofs(double* values,
1190
1168
const ufc::function& f,
1169
const double* vertex_coordinates,
1170
int cell_orientation,
1191
1171
const ufc::cell& c) const
1192
1172
{
1193
// Declare variables for result of evaluation.
1173
// Declare variables for result of evaluation
1194
1174
double vals[1];
1195
1175
1196
// Declare variable for physical coordinates.
1176
// Declare variable for physical coordinates
1197
1177
double y[3];
1198
const double * const * x = c.coordinates;
1199
y[0] = x[0][0];
1200
y[1] = x[0][1];
1201
y[2] = x[0][2];
1178
y[0] = vertex_coordinates[0];
1179
y[1] = vertex_coordinates[1];
1180
y[2] = vertex_coordinates[2];
1202
1181
f.evaluate(vals, y, c);
1203
1182
values[0] = vals[0];
1204
y[0] = x[1][0];
1205
y[1] = x[1][1];
1206
y[2] = x[1][2];
1183
y[0] = vertex_coordinates[3];
1184
y[1] = vertex_coordinates[4];
1185
y[2] = vertex_coordinates[5];
1207
1186
f.evaluate(vals, y, c);
1208
1187
values[1] = vals[0];
1209
y[0] = x[2][0];
1210
y[1] = x[2][1];
1211
y[2] = x[2][2];
1188
y[0] = vertex_coordinates[6];
1189
y[1] = vertex_coordinates[7];
1190
y[2] = vertex_coordinates[8];
1212
1191
f.evaluate(vals, y, c);
1213
1192
values[2] = vals[0];
1214
y[0] = x[3][0];
1215
y[1] = x[3][1];
1216
y[2] = x[3][2];
1193
y[0] = vertex_coordinates[9];
1194
y[1] = vertex_coordinates[10];
1195
y[2] = vertex_coordinates[11];
1217
1196
f.evaluate(vals, y, c);
1218
1197
values[3] = vals[0];
1219
1198
}
1221
1200
/// Interpolate vertex values from dof values
1222
1201
virtual void interpolate_vertex_values(double* vertex_values,
1223
1202
const double* dof_values,
1203
const double* vertex_coordinates,
1204
int cell_orientation,
1224
1205
const ufc::cell& c) const
1225
1206
{
1226
1207
// Evaluate function and change variables
1235
1216
const double* xhat,
1236
1217
const ufc::cell& c) const
1237
1218
{
1238
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
1219
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1239
1220
}
1240
1221
1241
1222
/// Map from coordinate x in cell to coordinate xhat in reference cell
1243
1224
const double* x,
1244
1225
const ufc::cell& c) const
1245
1226
{
1246
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
1227
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1247
1228
}
1248
1229
1249
1230
/// Return the number of sub elements (for a mixed element)
1250
virtual unsigned int num_sub_elements() const
1231
virtual std::size_t num_sub_elements() const
1251
1232
{
1252
1233
return 0;
1253
1234
}
1254
1235
1255
1236
/// Create a new finite element for sub element i (for a mixed element)
1256
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1237
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1257
1238
{
1258
1239
return 0;
1259
1240
}
1287
1268
/// Return a string identifying the finite element
1288
1269
virtual const char* signature() const
1289
1270
{
1290
return "VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None)";
1271
return "VectorElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, 3, None)";
1291
1272
}
1292
1273
1293
1274
/// Return the cell shape
1297
1278
}
1298
1279
1299
1280
/// Return the topological dimension of the cell shape
1300
virtual unsigned int topological_dimension() const
1281
virtual std::size_t topological_dimension() const
1301
1282
{
1302
1283
return 3;
1303
1284
}
1304
1285
1305
1286
/// Return the geometric dimension of the cell shape
1306
virtual unsigned int geometric_dimension() const
1287
virtual std::size_t geometric_dimension() const
1307
1288
{
1308
1289
return 3;
1309
1290
}
1310
1291
1311
1292
/// Return the dimension of the finite element function space
1312
virtual unsigned int space_dimension() const
1293
virtual std::size_t space_dimension() const
1313
1294
{
1314
1295
return 12;
1315
1296
}
1316
1297
1317
1298
/// Return the rank of the value space
1318
virtual unsigned int value_rank() const
1299
virtual std::size_t value_rank() const
1319
1300
{
1320
1301
return 1;
1321
1302
}
1322
1303
1323
1304
/// Return the dimension of the value space for axis i
1324
virtual unsigned int value_dimension(unsigned int i) const
1305
virtual std::size_t value_dimension(std::size_t i) const
1325
1306
{
1326
1307
switch (i)
1327
1308
{
1335
1316
return 0;
1336
1317
}
1337
1318
1338
/// Evaluate basis function i at given point in cell
1339
virtual void evaluate_basis(unsigned int i,
1319
/// Evaluate basis function i at given point x in cell
1320
virtual void evaluate_basis(std::size_t i,
1340
1321
double* values,
1341
const double* coordinates,
1342
const ufc::cell& c) const
1322
const double* x,
1323
const double* vertex_coordinates,
1324
int cell_orientation) const
1343
1325
{
1344
// Extract vertex coordinates
1345
const double * const * x = c.coordinates;
1346
1347
// Compute Jacobian of affine map from reference cell
1348
const double J_00 = x[1][0] - x[0][0];
1349
const double J_01 = x[2][0] - x[0][0];
1350
const double J_02 = x[3][0] - x[0][0];
1351
const double J_10 = x[1][1] - x[0][1];
1352
const double J_11 = x[2][1] - x[0][1];
1353
const double J_12 = x[3][1] - x[0][1];
1354
const double J_20 = x[1][2] - x[0][2];
1355
const double J_21 = x[2][2] - x[0][2];
1356
const double J_22 = x[3][2] - x[0][2];
1357
1358
// Compute sub determinants
1359
const double d_00 = J_11*J_22 - J_12*J_21;
1360
const double d_01 = J_12*J_20 - J_10*J_22;
1361
const double d_02 = J_10*J_21 - J_11*J_20;
1362
const double d_10 = J_02*J_21 - J_01*J_22;
1363
const double d_11 = J_00*J_22 - J_02*J_20;
1364
const double d_12 = J_01*J_20 - J_00*J_21;
1365
const double d_20 = J_01*J_12 - J_02*J_11;
1366
const double d_21 = J_02*J_10 - J_00*J_12;
1367
const double d_22 = J_00*J_11 - J_01*J_10;
1368
1369
// Compute determinant of Jacobian
1370
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
1371
1372
// Compute inverse of Jacobian
1326
// Compute Jacobian
1327
double J[9];
1328
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
1329
1330
// Compute Jacobian inverse and determinant
1331
double K[9];
1332
double detJ;
1333
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
1334
1373
1335
1374
1336
// Compute constants
1375
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
1376
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
1377
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
1337
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
1338
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
1339
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
1340
1341
// Compute subdeterminants
1342
const double d_00 = J[4]*J[8] - J[5]*J[7];
1343
const double d_01 = J[5]*J[6] - J[3]*J[8];
1344
const double d_02 = J[3]*J[7] - J[4]*J[6];
1345
const double d_10 = J[2]*J[7] - J[1]*J[8];
1346
const double d_11 = J[0]*J[8] - J[2]*J[6];
1347
const double d_12 = J[1]*J[6] - J[0]*J[7];
1348
const double d_20 = J[1]*J[5] - J[2]*J[4];
1349
const double d_21 = J[2]*J[3] - J[0]*J[5];
1350
const double d_22 = J[0]*J[4] - J[1]*J[3];
1378
1351
1379
1352
// Get coordinates and map to the reference (FIAT) element
1380
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
1381
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
1382
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
1383
1384
1385
// Reset values.
1353
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
1354
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
1355
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
1356
1357
1358
// Reset values
1386
1359
values[0] = 0.0;
1387
1360
values[1] = 0.0;
1388
1361
values[2] = 0.0;
1391
1364
case 0:
1392
1365
{
1393
1366
1394
// Array of basisvalues.
1367
// Array of basisvalues
1395
1368
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1396
1369
1397
// Declare helper variables.
1370
// Declare helper variables
1398
1371
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1399
1372
1400
// Compute basisvalues.
1373
// Compute basisvalues
1401
1374
basisvalues[0] = 1.0;
1402
1375
basisvalues[1] = tmp0;
1403
1376
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1407
1380
basisvalues[2] *= std::sqrt(2.5);
1408
1381
basisvalues[1] *= std::sqrt(7.5);
1409
1382
1410
// Table(s) of coefficients.
1383
// Table(s) of coefficients
1411
1384
static const double coefficients0[4] = \
1412
1385
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1413
1386
1414
// Compute value(s).
1387
// Compute value(s)
1415
1388
for (unsigned int r = 0; r < 4; r++)
1416
1389
{
1417
1390
values[0] += coefficients0[r]*basisvalues[r];
1421
1394
case 1:
1422
1395
{
1423
1396
1424
// Array of basisvalues.
1397
// Array of basisvalues
1425
1398
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1426
1399
1427
// Declare helper variables.
1400
// Declare helper variables
1428
1401
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1429
1402
1430
// Compute basisvalues.
1403
// Compute basisvalues
1431
1404
basisvalues[0] = 1.0;
1432
1405
basisvalues[1] = tmp0;
1433
1406
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1437
1410
basisvalues[2] *= std::sqrt(2.5);
1438
1411
basisvalues[1] *= std::sqrt(7.5);
1439
1412
1440
// Table(s) of coefficients.
1413
// Table(s) of coefficients
1441
1414
static const double coefficients0[4] = \
1442
1415
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
1443
1416
1444
// Compute value(s).
1417
// Compute value(s)
1445
1418
for (unsigned int r = 0; r < 4; r++)
1446
1419
{
1447
1420
values[0] += coefficients0[r]*basisvalues[r];
1451
1424
case 2:
1452
1425
{
1453
1426
1454
// Array of basisvalues.
1427
// Array of basisvalues
1455
1428
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1456
1429
1457
// Declare helper variables.
1430
// Declare helper variables
1458
1431
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1459
1432
1460
// Compute basisvalues.
1433
// Compute basisvalues
1461
1434
basisvalues[0] = 1.0;
1462
1435
basisvalues[1] = tmp0;
1463
1436
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1467
1440
basisvalues[2] *= std::sqrt(2.5);
1468
1441
basisvalues[1] *= std::sqrt(7.5);
1469
1442
1470
// Table(s) of coefficients.
1443
// Table(s) of coefficients
1471
1444
static const double coefficients0[4] = \
1472
1445
{0.28867513, 0.0, 0.21081851, -0.074535599};
1473
1446
1474
// Compute value(s).
1447
// Compute value(s)
1475
1448
for (unsigned int r = 0; r < 4; r++)
1476
1449
{
1477
1450
values[0] += coefficients0[r]*basisvalues[r];
1481
1454
case 3:
1482
1455
{
1483
1456
1484
// Array of basisvalues.
1457
// Array of basisvalues
1485
1458
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1486
1459
1487
// Declare helper variables.
1460
// Declare helper variables
1488
1461
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1489
1462
1490
// Compute basisvalues.
1463
// Compute basisvalues
1491
1464
basisvalues[0] = 1.0;
1492
1465
basisvalues[1] = tmp0;
1493
1466
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1497
1470
basisvalues[2] *= std::sqrt(2.5);
1498
1471
basisvalues[1] *= std::sqrt(7.5);
1499
1472
1500
// Table(s) of coefficients.
1473
// Table(s) of coefficients
1501
1474
static const double coefficients0[4] = \
1502
1475
{0.28867513, 0.0, 0.0, 0.2236068};
1503
1476
1504
// Compute value(s).
1477
// Compute value(s)
1505
1478
for (unsigned int r = 0; r < 4; r++)
1506
1479
{
1507
1480
values[0] += coefficients0[r]*basisvalues[r];
1511
1484
case 4:
1512
1485
{
1513
1486
1514
// Array of basisvalues.
1487
// Array of basisvalues
1515
1488
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1516
1489
1517
// Declare helper variables.
1490
// Declare helper variables
1518
1491
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1519
1492
1520
// Compute basisvalues.
1493
// Compute basisvalues
1521
1494
basisvalues[0] = 1.0;
1522
1495
basisvalues[1] = tmp0;
1523
1496
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1527
1500
basisvalues[2] *= std::sqrt(2.5);
1528
1501
basisvalues[1] *= std::sqrt(7.5);
1529
1502
1530
// Table(s) of coefficients.
1503
// Table(s) of coefficients
1531
1504
static const double coefficients0[4] = \
1532
1505
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1533
1506
1534
// Compute value(s).
1507
// Compute value(s)
1535
1508
for (unsigned int r = 0; r < 4; r++)
1536
1509
{
1537
1510
values[1] += coefficients0[r]*basisvalues[r];
1541
1514
case 5:
1542
1515
{
1543
1516
1544
// Array of basisvalues.
1517
// Array of basisvalues
1545
1518
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1546
1519
1547
// Declare helper variables.
1520
// Declare helper variables
1548
1521
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1549
1522
1550
// Compute basisvalues.
1523
// Compute basisvalues
1551
1524
basisvalues[0] = 1.0;
1552
1525
basisvalues[1] = tmp0;
1553
1526
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1557
1530
basisvalues[2] *= std::sqrt(2.5);
1558
1531
basisvalues[1] *= std::sqrt(7.5);
1559
1532
1560
// Table(s) of coefficients.
1533
// Table(s) of coefficients
1561
1534
static const double coefficients0[4] = \
1562
1535
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
1563
1536
1564
// Compute value(s).
1537
// Compute value(s)
1565
1538
for (unsigned int r = 0; r < 4; r++)
1566
1539
{
1567
1540
values[1] += coefficients0[r]*basisvalues[r];
1571
1544
case 6:
1572
1545
{
1573
1546
1574
// Array of basisvalues.
1547
// Array of basisvalues
1575
1548
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1576
1549
1577
// Declare helper variables.
1550
// Declare helper variables
1578
1551
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1579
1552
1580
// Compute basisvalues.
1553
// Compute basisvalues
1581
1554
basisvalues[0] = 1.0;
1582
1555
basisvalues[1] = tmp0;
1583
1556
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1587
1560
basisvalues[2] *= std::sqrt(2.5);
1588
1561
basisvalues[1] *= std::sqrt(7.5);
1589
1562
1590
// Table(s) of coefficients.
1563
// Table(s) of coefficients
1591
1564
static const double coefficients0[4] = \
1592
1565
{0.28867513, 0.0, 0.21081851, -0.074535599};
1593
1566
1594
// Compute value(s).
1567
// Compute value(s)
1595
1568
for (unsigned int r = 0; r < 4; r++)
1596
1569
{
1597
1570
values[1] += coefficients0[r]*basisvalues[r];
1601
1574
case 7:
1602
1575
{
1603
1576
1604
// Array of basisvalues.
1577
// Array of basisvalues
1605
1578
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1606
1579
1607
// Declare helper variables.
1580
// Declare helper variables
1608
1581
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1609
1582
1610
// Compute basisvalues.
1583
// Compute basisvalues
1611
1584
basisvalues[0] = 1.0;
1612
1585
basisvalues[1] = tmp0;
1613
1586
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1617
1590
basisvalues[2] *= std::sqrt(2.5);
1618
1591
basisvalues[1] *= std::sqrt(7.5);
1619
1592
1620
// Table(s) of coefficients.
1593
// Table(s) of coefficients
1621
1594
static const double coefficients0[4] = \
1622
1595
{0.28867513, 0.0, 0.0, 0.2236068};
1623
1596
1624
// Compute value(s).
1597
// Compute value(s)
1625
1598
for (unsigned int r = 0; r < 4; r++)
1626
1599
{
1627
1600
values[1] += coefficients0[r]*basisvalues[r];
1631
1604
case 8:
1632
1605
{
1633
1606
1634
// Array of basisvalues.
1607
// Array of basisvalues
1635
1608
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1636
1609
1637
// Declare helper variables.
1610
// Declare helper variables
1638
1611
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1639
1612
1640
// Compute basisvalues.
1613
// Compute basisvalues
1641
1614
basisvalues[0] = 1.0;
1642
1615
basisvalues[1] = tmp0;
1643
1616
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1647
1620
basisvalues[2] *= std::sqrt(2.5);
1648
1621
basisvalues[1] *= std::sqrt(7.5);
1649
1622
1650
// Table(s) of coefficients.
1623
// Table(s) of coefficients
1651
1624
static const double coefficients0[4] = \
1652
1625
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1653
1626
1654
// Compute value(s).
1627
// Compute value(s)
1655
1628
for (unsigned int r = 0; r < 4; r++)
1656
1629
{
1657
1630
values[2] += coefficients0[r]*basisvalues[r];
1661
1634
case 9:
1662
1635
{
1663
1636
1664
// Array of basisvalues.
1637
// Array of basisvalues
1665
1638
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1666
1639
1667
// Declare helper variables.
1640
// Declare helper variables
1668
1641
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1669
1642
1670
// Compute basisvalues.
1643
// Compute basisvalues
1671
1644
basisvalues[0] = 1.0;
1672
1645
basisvalues[1] = tmp0;
1673
1646
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1677
1650
basisvalues[2] *= std::sqrt(2.5);
1678
1651
basisvalues[1] *= std::sqrt(7.5);
1679
1652
1680
// Table(s) of coefficients.
1653
// Table(s) of coefficients
1681
1654
static const double coefficients0[4] = \
1682
1655
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
1683
1656
1684
// Compute value(s).
1657
// Compute value(s)
1685
1658
for (unsigned int r = 0; r < 4; r++)
1686
1659
{
1687
1660
values[2] += coefficients0[r]*basisvalues[r];
1691
1664
case 10:
1692
1665
{
1693
1666
1694
// Array of basisvalues.
1667
// Array of basisvalues
1695
1668
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1696
1669
1697
// Declare helper variables.
1670
// Declare helper variables
1698
1671
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1699
1672
1700
// Compute basisvalues.
1673
// Compute basisvalues
1701
1674
basisvalues[0] = 1.0;
1702
1675
basisvalues[1] = tmp0;
1703
1676
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1707
1680
basisvalues[2] *= std::sqrt(2.5);
1708
1681
basisvalues[1] *= std::sqrt(7.5);
1709
1682
1710
// Table(s) of coefficients.
1683
// Table(s) of coefficients
1711
1684
static const double coefficients0[4] = \
1712
1685
{0.28867513, 0.0, 0.21081851, -0.074535599};
1713
1686
1714
// Compute value(s).
1687
// Compute value(s)
1715
1688
for (unsigned int r = 0; r < 4; r++)
1716
1689
{
1717
1690
values[2] += coefficients0[r]*basisvalues[r];
1721
1694
case 11:
1722
1695
{
1723
1696
1724
// Array of basisvalues.
1697
// Array of basisvalues
1725
1698
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1726
1699
1727
// Declare helper variables.
1700
// Declare helper variables
1728
1701
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1729
1702
1730
// Compute basisvalues.
1703
// Compute basisvalues
1731
1704
basisvalues[0] = 1.0;
1732
1705
basisvalues[1] = tmp0;
1733
1706
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1737
1710
basisvalues[2] *= std::sqrt(2.5);
1738
1711
basisvalues[1] *= std::sqrt(7.5);
1739
1712
1740
// Table(s) of coefficients.
1713
// Table(s) of coefficients
1741
1714
static const double coefficients0[4] = \
1742
1715
{0.28867513, 0.0, 0.0, 0.2236068};
1743
1716
1744
// Compute value(s).
1717
// Compute value(s)
1745
1718
for (unsigned int r = 0; r < 4; r++)
1746
1719
{
1747
1720
values[2] += coefficients0[r]*basisvalues[r];
1752
1725
1753
1726
}
1754
1727
1755
/// Evaluate all basis functions at given point in cell
1728
/// Evaluate all basis functions at given point x in cell
1756
1729
virtual void evaluate_basis_all(double* values,
1757
const double* coordinates,
1758
const ufc::cell& c) const
1730
const double* x,
1731
const double* vertex_coordinates,
1732
int cell_orientation) const
1759
1733
{
1760
1734
// Helper variable to hold values of a single dof.
1761
1735
double dof_values[3] = {0.0, 0.0, 0.0};
1762
1736
1763
// Loop dofs and call evaluate_basis.
1737
// Loop dofs and call evaluate_basis
1764
1738
for (unsigned int r = 0; r < 12; r++)
1765
1739
{
1766
evaluate_basis(r, dof_values, coordinates, c);
1740
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
1767
1741
for (unsigned int s = 0; s < 3; s++)
1768
1742
{
1769
1743
values[r*3 + s] = dof_values[s];
1771
1745
}// end loop over 'r'
1772
1746
}
1773
1747
1774
/// Evaluate order n derivatives of basis function i at given point in cell
1775
virtual void evaluate_basis_derivatives(unsigned int i,
1776
unsigned int n,
1748
/// Evaluate order n derivatives of basis function i at given point x in cell
1749
virtual void evaluate_basis_derivatives(std::size_t i,
1750
std::size_t n,
1777
1751
double* values,
1778
const double* coordinates,
1779
const ufc::cell& c) const
1752
const double* x,
1753
const double* vertex_coordinates,
1754
int cell_orientation) const
1780
1755
{
1781
// Extract vertex coordinates
1782
const double * const * x = c.coordinates;
1783
1784
// Compute Jacobian of affine map from reference cell
1785
const double J_00 = x[1][0] - x[0][0];
1786
const double J_01 = x[2][0] - x[0][0];
1787
const double J_02 = x[3][0] - x[0][0];
1788
const double J_10 = x[1][1] - x[0][1];
1789
const double J_11 = x[2][1] - x[0][1];
1790
const double J_12 = x[3][1] - x[0][1];
1791
const double J_20 = x[1][2] - x[0][2];
1792
const double J_21 = x[2][2] - x[0][2];
1793
const double J_22 = x[3][2] - x[0][2];
1794
1795
// Compute sub determinants
1796
const double d_00 = J_11*J_22 - J_12*J_21;
1797
const double d_01 = J_12*J_20 - J_10*J_22;
1798
const double d_02 = J_10*J_21 - J_11*J_20;
1799
const double d_10 = J_02*J_21 - J_01*J_22;
1800
const double d_11 = J_00*J_22 - J_02*J_20;
1801
const double d_12 = J_01*J_20 - J_00*J_21;
1802
const double d_20 = J_01*J_12 - J_02*J_11;
1803
const double d_21 = J_02*J_10 - J_00*J_12;
1804
const double d_22 = J_00*J_11 - J_01*J_10;
1805
1806
// Compute determinant of Jacobian
1807
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
1808
1809
// Compute inverse of Jacobian
1810
const double K_00 = d_00 / detJ;
1811
const double K_01 = d_10 / detJ;
1812
const double K_02 = d_20 / detJ;
1813
const double K_10 = d_01 / detJ;
1814
const double K_11 = d_11 / detJ;
1815
const double K_12 = d_21 / detJ;
1816
const double K_20 = d_02 / detJ;
1817
const double K_21 = d_12 / detJ;
1818
const double K_22 = d_22 / detJ;
1756
// Compute Jacobian
1757
double J[9];
1758
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
1759
1760
// Compute Jacobian inverse and determinant
1761
double K[9];
1762
double detJ;
1763
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
1764
1819
1765
1820
1766
// Compute constants
1821
const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
1822
const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
1823
const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
1767
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
1768
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
1769
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
1770
1771
// Compute subdeterminants
1772
const double d_00 = J[4]*J[8] - J[5]*J[7];
1773
const double d_01 = J[5]*J[6] - J[3]*J[8];
1774
const double d_02 = J[3]*J[7] - J[4]*J[6];
1775
const double d_10 = J[2]*J[7] - J[1]*J[8];
1776
const double d_11 = J[0]*J[8] - J[2]*J[6];
1777
const double d_12 = J[1]*J[6] - J[0]*J[7];
1778
const double d_20 = J[1]*J[5] - J[2]*J[4];
1779
const double d_21 = J[2]*J[3] - J[0]*J[5];
1780
const double d_22 = J[0]*J[4] - J[1]*J[3];
1824
1781
1825
1782
// Get coordinates and map to the reference (FIAT) element
1826
double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
1827
double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
1828
double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
1783
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
1784
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
1785
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
1829
1786
1830
1787
1831
1788
// Compute number of derivatives.
1863
1820
}
1864
1821
1865
1822
// Compute inverse of Jacobian
1866
const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}};
1823
const double Jinv[3][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}, {K[6], K[7], K[8]}};
1867
1824
1868
1825
// Declare transformation matrix
1869
1826
// Declare pointer to two dimensional array and initialise
1897
1854
case 0:
1898
1855
{
1899
1856
1900
// Array of basisvalues.
1857
// Array of basisvalues
1901
1858
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1902
1859
1903
// Declare helper variables.
1860
// Declare helper variables
1904
1861
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1905
1862
1906
// Compute basisvalues.
1863
// Compute basisvalues
1907
1864
basisvalues[0] = 1.0;
1908
1865
basisvalues[1] = tmp0;
1909
1866
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1913
1870
basisvalues[2] *= std::sqrt(2.5);
1914
1871
basisvalues[1] *= std::sqrt(7.5);
1915
1872
1916
// Table(s) of coefficients.
1873
// Table(s) of coefficients
1917
1874
static const double coefficients0[4] = \
1918
1875
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
1919
1876
2069
2026
case 1:
2070
2027
{
2071
2028
2072
// Array of basisvalues.
2029
// Array of basisvalues
2073
2030
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2074
2031
2075
// Declare helper variables.
2032
// Declare helper variables
2076
2033
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2077
2034
2078
// Compute basisvalues.
2035
// Compute basisvalues
2079
2036
basisvalues[0] = 1.0;
2080
2037
basisvalues[1] = tmp0;
2081
2038
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2085
2042
basisvalues[2] *= std::sqrt(2.5);
2086
2043
basisvalues[1] *= std::sqrt(7.5);
2087
2044
2088
// Table(s) of coefficients.
2045
// Table(s) of coefficients
2089
2046
static const double coefficients0[4] = \
2090
2047
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
2091
2048
2241
2198
case 2:
2242
2199
{
2243
2200
2244
// Array of basisvalues.
2201
// Array of basisvalues
2245
2202
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2246
2203
2247
// Declare helper variables.
2204
// Declare helper variables
2248
2205
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2249
2206
2250
// Compute basisvalues.
2207
// Compute basisvalues
2251
2208
basisvalues[0] = 1.0;
2252
2209
basisvalues[1] = tmp0;
2253
2210
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2257
2214
basisvalues[2] *= std::sqrt(2.5);
2258
2215
basisvalues[1] *= std::sqrt(7.5);
2259
2216
2260
// Table(s) of coefficients.
2217
// Table(s) of coefficients
2261
2218
static const double coefficients0[4] = \
2262
2219
{0.28867513, 0.0, 0.21081851, -0.074535599};
2263
2220
2413
2370
case 3:
2414
2371
{
2415
2372
2416
// Array of basisvalues.
2373
// Array of basisvalues
2417
2374
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2418
2375
2419
// Declare helper variables.
2376
// Declare helper variables
2420
2377
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2421
2378
2422
// Compute basisvalues.
2379
// Compute basisvalues
2423
2380
basisvalues[0] = 1.0;
2424
2381
basisvalues[1] = tmp0;
2425
2382
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2429
2386
basisvalues[2] *= std::sqrt(2.5);
2430
2387
basisvalues[1] *= std::sqrt(7.5);
2431
2388
2432
// Table(s) of coefficients.
2389
// Table(s) of coefficients
2433
2390
static const double coefficients0[4] = \
2434
2391
{0.28867513, 0.0, 0.0, 0.2236068};
2435
2392
2585
2542
case 4:
2586
2543
{
2587
2544
2588
// Array of basisvalues.
2545
// Array of basisvalues
2589
2546
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2590
2547
2591
// Declare helper variables.
2548
// Declare helper variables
2592
2549
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2593
2550
2594
// Compute basisvalues.
2551
// Compute basisvalues
2595
2552
basisvalues[0] = 1.0;
2596
2553
basisvalues[1] = tmp0;
2597
2554
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2601
2558
basisvalues[2] *= std::sqrt(2.5);
2602
2559
basisvalues[1] *= std::sqrt(7.5);
2603
2560
2604
// Table(s) of coefficients.
2561
// Table(s) of coefficients
2605
2562
static const double coefficients0[4] = \
2606
2563
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
2607
2564
2757
2714
case 5:
2758
2715
{
2759
2716
2760
// Array of basisvalues.
2717
// Array of basisvalues
2761
2718
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2762
2719
2763
// Declare helper variables.
2720
// Declare helper variables
2764
2721
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2765
2722
2766
// Compute basisvalues.
2723
// Compute basisvalues
2767
2724
basisvalues[0] = 1.0;
2768
2725
basisvalues[1] = tmp0;
2769
2726
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2773
2730
basisvalues[2] *= std::sqrt(2.5);
2774
2731
basisvalues[1] *= std::sqrt(7.5);
2775
2732
2776
// Table(s) of coefficients.
2733
// Table(s) of coefficients
2777
2734
static const double coefficients0[4] = \
2778
2735
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
2779
2736
2929
2886
case 6:
2930
2887
{
2931
2888
2932
// Array of basisvalues.
2889
// Array of basisvalues
2933
2890
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2934
2891
2935
// Declare helper variables.
2892
// Declare helper variables
2936
2893
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2937
2894
2938
// Compute basisvalues.
2895
// Compute basisvalues
2939
2896
basisvalues[0] = 1.0;
2940
2897
basisvalues[1] = tmp0;
2941
2898
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2945
2902
basisvalues[2] *= std::sqrt(2.5);
2946
2903
basisvalues[1] *= std::sqrt(7.5);
2947
2904
2948
// Table(s) of coefficients.
2905
// Table(s) of coefficients
2949
2906
static const double coefficients0[4] = \
2950
2907
{0.28867513, 0.0, 0.21081851, -0.074535599};
2951
2908
3101
3058
case 7:
3102
3059
{
3103
3060
3104
// Array of basisvalues.
3061
// Array of basisvalues
3105
3062
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3106
3063
3107
// Declare helper variables.
3064
// Declare helper variables
3108
3065
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3109
3066
3110
// Compute basisvalues.
3067
// Compute basisvalues
3111
3068
basisvalues[0] = 1.0;
3112
3069
basisvalues[1] = tmp0;
3113
3070
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3117
3074
basisvalues[2] *= std::sqrt(2.5);
3118
3075
basisvalues[1] *= std::sqrt(7.5);
3119
3076
3120
// Table(s) of coefficients.
3077
// Table(s) of coefficients
3121
3078
static const double coefficients0[4] = \
3122
3079
{0.28867513, 0.0, 0.0, 0.2236068};
3123
3080
3273
3230
case 8:
3274
3231
{
3275
3232
3276
// Array of basisvalues.
3233
// Array of basisvalues
3277
3234
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3278
3235
3279
// Declare helper variables.
3236
// Declare helper variables
3280
3237
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3281
3238
3282
// Compute basisvalues.
3239
// Compute basisvalues
3283
3240
basisvalues[0] = 1.0;
3284
3241
basisvalues[1] = tmp0;
3285
3242
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3289
3246
basisvalues[2] *= std::sqrt(2.5);
3290
3247
basisvalues[1] *= std::sqrt(7.5);
3291
3248
3292
// Table(s) of coefficients.
3249
// Table(s) of coefficients
3293
3250
static const double coefficients0[4] = \
3294
3251
{0.28867513, -0.18257419, -0.10540926, -0.074535599};
3295
3252
3445
3402
case 9:
3446
3403
{
3447
3404
3448
// Array of basisvalues.
3405
// Array of basisvalues
3449
3406
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3450
3407
3451
// Declare helper variables.
3408
// Declare helper variables
3452
3409
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3453
3410
3454
// Compute basisvalues.
3411
// Compute basisvalues
3455
3412
basisvalues[0] = 1.0;
3456
3413
basisvalues[1] = tmp0;
3457
3414
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3461
3418
basisvalues[2] *= std::sqrt(2.5);
3462
3419
basisvalues[1] *= std::sqrt(7.5);
3463
3420
3464
// Table(s) of coefficients.
3421
// Table(s) of coefficients
3465
3422
static const double coefficients0[4] = \
3466
3423
{0.28867513, 0.18257419, -0.10540926, -0.074535599};
3467
3424
3617
3574
case 10:
3618
3575
{
3619
3576
3620
// Array of basisvalues.
3577
// Array of basisvalues
3621
3578
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3622
3579
3623
// Declare helper variables.
3580
// Declare helper variables
3624
3581
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3625
3582
3626
// Compute basisvalues.
3583
// Compute basisvalues
3627
3584
basisvalues[0] = 1.0;
3628
3585
basisvalues[1] = tmp0;
3629
3586
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3633
3590
basisvalues[2] *= std::sqrt(2.5);
3634
3591
basisvalues[1] *= std::sqrt(7.5);
3635
3592
3636
// Table(s) of coefficients.
3593
// Table(s) of coefficients
3637
3594
static const double coefficients0[4] = \
3638
3595
{0.28867513, 0.0, 0.21081851, -0.074535599};
3639
3596
3789
3746
case 11:
3790
3747
{
3791
3748
3792
// Array of basisvalues.
3749
// Array of basisvalues
3793
3750
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3794
3751
3795
// Declare helper variables.
3752
// Declare helper variables
3796
3753
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3797
3754
3798
// Compute basisvalues.
3755
// Compute basisvalues
3799
3756
basisvalues[0] = 1.0;
3800
3757
basisvalues[1] = tmp0;
3801
3758
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3805
3762
basisvalues[2] *= std::sqrt(2.5);
3806
3763
basisvalues[1] *= std::sqrt(7.5);
3807
3764
3808
// Table(s) of coefficients.
3765
// Table(s) of coefficients
3809
3766
static const double coefficients0[4] = \
3810
3767
{0.28867513, 0.0, 0.0, 0.2236068};
3811
3768
3962
3919
3963
3920
}
3964
3921
3965
/// Evaluate order n derivatives of all basis functions at given point in cell
3966
virtual void evaluate_basis_derivatives_all(unsigned int n,
3922
/// Evaluate order n derivatives of all basis functions at given point x in cell
3923
virtual void evaluate_basis_derivatives_all(std::size_t n,
3967
3924
double* values,
3968
const double* coordinates,
3969
const ufc::cell& c) const
3925
const double* x,
3926
const double* vertex_coordinates,
3927
int cell_orientation) const
3970
3928
{
3971
3929
// Compute number of derivatives.
3972
3930
unsigned int num_derivatives = 1;
3985
3943
// Loop dofs and call evaluate_basis_derivatives.
3986
3944
for (unsigned int r = 0; r < 12; r++)
3987
3945
{
3988
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
3946
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
3989
3947
for (unsigned int s = 0; s < 3*num_derivatives; s++)
3990
3948
{
3991
3949
values[r*3*num_derivatives + s] = dof_values[s];
3997
3955
}
3998
3956
3999
3957
/// Evaluate linear functional for dof i on the function f
4000
virtual double evaluate_dof(unsigned int i,
3958
virtual double evaluate_dof(std::size_t i,
4001
3959
const ufc::function& f,
3960
const double* vertex_coordinates,
3961
int cell_orientation,
4002
3962
const ufc::cell& c) const
4003
3963
{
4004
// Declare variables for result of evaluation.
3964
// Declare variables for result of evaluation
4005
3965
double vals[3];
4006
3966
4007
// Declare variable for physical coordinates.
3967
// Declare variable for physical coordinates
4008
3968
double y[3];
4009
const double * const * x = c.coordinates;
4010
3969
switch (i)
4011
3970
{
4012
3971
case 0:
4013
3972
{
4014
y[0] = x[0][0];
4015
y[1] = x[0][1];
4016
y[2] = x[0][2];
3973
y[0] = vertex_coordinates[0];
3974
y[1] = vertex_coordinates[1];
3975
y[2] = vertex_coordinates[2];
4017
3976
f.evaluate(vals, y, c);
4018
3977
return vals[0];
4019
3978
break;
4020
3979
}
4021
3980
case 1:
4022
3981
{
4023
y[0] = x[1][0];
4024
y[1] = x[1][1];
4025
y[2] = x[1][2];
3982
y[0] = vertex_coordinates[3];
3983
y[1] = vertex_coordinates[4];
3984
y[2] = vertex_coordinates[5];
4026
3985
f.evaluate(vals, y, c);
4027
3986
return vals[0];
4028
3987
break;
4029
3988
}
4030
3989
case 2:
4031
3990
{
4032
y[0] = x[2][0];
4033
y[1] = x[2][1];
4034
y[2] = x[2][2];
3991
y[0] = vertex_coordinates[6];
3992
y[1] = vertex_coordinates[7];
3993
y[2] = vertex_coordinates[8];
4035
3994
f.evaluate(vals, y, c);
4036
3995
return vals[0];
4037
3996
break;
4038
3997
}
4039
3998
case 3:
4040
3999
{
4041
y[0] = x[3][0];
4042
y[1] = x[3][1];
4043
y[2] = x[3][2];
4000
y[0] = vertex_coordinates[9];
4001
y[1] = vertex_coordinates[10];
4002
y[2] = vertex_coordinates[11];
4044
4003
f.evaluate(vals, y, c);
4045
4004
return vals[0];
4046
4005
break;
4047
4006
}
4048
4007
case 4:
4049
4008
{
4050
y[0] = x[0][0];
4051
y[1] = x[0][1];
4052
y[2] = x[0][2];
4009
y[0] = vertex_coordinates[0];
4010
y[1] = vertex_coordinates[1];
4011
y[2] = vertex_coordinates[2];
4053
4012
f.evaluate(vals, y, c);
4054
4013
return vals[1];
4055
4014
break;
4056
4015
}
4057
4016
case 5:
4058
4017
{
4059
y[0] = x[1][0];
4060
y[1] = x[1][1];
4061
y[2] = x[1][2];
4018
y[0] = vertex_coordinates[3];
4019
y[1] = vertex_coordinates[4];
4020
y[2] = vertex_coordinates[5];
4062
4021
f.evaluate(vals, y, c);
4063
4022
return vals[1];
4064
4023
break;
4065
4024
}
4066
4025
case 6:
4067
4026
{
4068
y[0] = x[2][0];
4069
y[1] = x[2][1];
4070
y[2] = x[2][2];
4027
y[0] = vertex_coordinates[6];
4028
y[1] = vertex_coordinates[7];
4029
y[2] = vertex_coordinates[8];
4071
4030
f.evaluate(vals, y, c);
4072
4031
return vals[1];
4073
4032
break;
4074
4033
}
4075
4034
case 7:
4076
4035
{
4077
y[0] = x[3][0];
4078
y[1] = x[3][1];
4079
y[2] = x[3][2];
4036
y[0] = vertex_coordinates[9];
4037
y[1] = vertex_coordinates[10];
4038
y[2] = vertex_coordinates[11];
4080
4039
f.evaluate(vals, y, c);
4081
4040
return vals[1];
4082
4041
break;
4083
4042
}
4084
4043
case 8:
4085
4044
{
4086
y[0] = x[0][0];
4087
y[1] = x[0][1];
4088
y[2] = x[0][2];
4045
y[0] = vertex_coordinates[0];
4046
y[1] = vertex_coordinates[1];
4047
y[2] = vertex_coordinates[2];
4089
4048
f.evaluate(vals, y, c);
4090
4049
return vals[2];
4091
4050
break;
4092
4051
}
4093
4052
case 9:
4094
4053
{
4095
y[0] = x[1][0];
4096
y[1] = x[1][1];
4097
y[2] = x[1][2];
4054
y[0] = vertex_coordinates[3];
4055
y[1] = vertex_coordinates[4];
4056
y[2] = vertex_coordinates[5];
4098
4057
f.evaluate(vals, y, c);
4099
4058
return vals[2];
4100
4059
break;
4101
4060
}
4102
4061
case 10:
4103
4062
{
4104
y[0] = x[2][0];
4105
y[1] = x[2][1];
4106
y[2] = x[2][2];
4063
y[0] = vertex_coordinates[6];
4064
y[1] = vertex_coordinates[7];
4065
y[2] = vertex_coordinates[8];
4107
4066
f.evaluate(vals, y, c);
4108
4067
return vals[2];
4109
4068
break;
4110
4069
}
4111
4070
case 11:
4112
4071
{
4113
y[0] = x[3][0];
4114
y[1] = x[3][1];
4115
y[2] = x[3][2];
4072
y[0] = vertex_coordinates[9];
4073
y[1] = vertex_coordinates[10];
4074
y[2] = vertex_coordinates[11];
4116
4075
f.evaluate(vals, y, c);
4117
4076
return vals[2];
4118
4077
break;
4125
4084
/// Evaluate linear functionals for all dofs on the function f
4126
4085
virtual void evaluate_dofs(double* values,
4127
4086
const ufc::function& f,
4087
const double* vertex_coordinates,
4088
int cell_orientation,
4128
4089
const ufc::cell& c) const
4129
4090
{
4130
// Declare variables for result of evaluation.
4091
// Declare variables for result of evaluation
4131
4092
double vals[3];
4132
4093
4133
// Declare variable for physical coordinates.
4094
// Declare variable for physical coordinates
4134
4095
double y[3];
4135
const double * const * x = c.coordinates;
4136
y[0] = x[0][0];
4137
y[1] = x[0][1];
4138
y[2] = x[0][2];
4096
y[0] = vertex_coordinates[0];
4097
y[1] = vertex_coordinates[1];
4098
y[2] = vertex_coordinates[2];
4139
4099
f.evaluate(vals, y, c);
4140
4100
values[0] = vals[0];
4141
y[0] = x[1][0];
4142
y[1] = x[1][1];
4143
y[2] = x[1][2];
4101
y[0] = vertex_coordinates[3];
4102
y[1] = vertex_coordinates[4];
4103
y[2] = vertex_coordinates[5];
4144
4104
f.evaluate(vals, y, c);
4145
4105
values[1] = vals[0];
4146
y[0] = x[2][0];
4147
y[1] = x[2][1];
4148
y[2] = x[2][2];
4106
y[0] = vertex_coordinates[6];
4107
y[1] = vertex_coordinates[7];
4108
y[2] = vertex_coordinates[8];
4149
4109
f.evaluate(vals, y, c);
4150
4110
values[2] = vals[0];
4151
y[0] = x[3][0];
4152
y[1] = x[3][1];
4153
y[2] = x[3][2];
4111
y[0] = vertex_coordinates[9];
4112
y[1] = vertex_coordinates[10];
4113
y[2] = vertex_coordinates[11];
4154
4114
f.evaluate(vals, y, c);
4155
4115
values[3] = vals[0];
4156
y[0] = x[0][0];
4157
y[1] = x[0][1];
4158
y[2] = x[0][2];
4116
y[0] = vertex_coordinates[0];
4117
y[1] = vertex_coordinates[1];
4118
y[2] = vertex_coordinates[2];
4159
4119
f.evaluate(vals, y, c);
4160
4120
values[4] = vals[1];
4161
y[0] = x[1][0];
4162
y[1] = x[1][1];
4163
y[2] = x[1][2];
4121
y[0] = vertex_coordinates[3];
4122
y[1] = vertex_coordinates[4];
4123
y[2] = vertex_coordinates[5];
4164
4124
f.evaluate(vals, y, c);
4165
4125
values[5] = vals[1];
4166
y[0] = x[2][0];
4167
y[1] = x[2][1];
4168
y[2] = x[2][2];
4126
y[0] = vertex_coordinates[6];
4127
y[1] = vertex_coordinates[7];
4128
y[2] = vertex_coordinates[8];
4169
4129
f.evaluate(vals, y, c);
4170
4130
values[6] = vals[1];
4171
y[0] = x[3][0];
4172
y[1] = x[3][1];
4173
y[2] = x[3][2];
4131
y[0] = vertex_coordinates[9];
4132
y[1] = vertex_coordinates[10];
4133
y[2] = vertex_coordinates[11];
4174
4134
f.evaluate(vals, y, c);
4175
4135
values[7] = vals[1];
4176
y[0] = x[0][0];
4177
y[1] = x[0][1];
4178
y[2] = x[0][2];
4136
y[0] = vertex_coordinates[0];
4137
y[1] = vertex_coordinates[1];
4138
y[2] = vertex_coordinates[2];
4179
4139
f.evaluate(vals, y, c);
4180
4140
values[8] = vals[2];
4181
y[0] = x[1][0];
4182
y[1] = x[1][1];
4183
y[2] = x[1][2];
4141
y[0] = vertex_coordinates[3];
4142
y[1] = vertex_coordinates[4];
4143
y[2] = vertex_coordinates[5];
4184
4144
f.evaluate(vals, y, c);
4185
4145
values[9] = vals[2];
4186
y[0] = x[2][0];
4187
y[1] = x[2][1];
4188
y[2] = x[2][2];
4146
y[0] = vertex_coordinates[6];
4147
y[1] = vertex_coordinates[7];
4148
y[2] = vertex_coordinates[8];
4189
4149
f.evaluate(vals, y, c);
4190
4150
values[10] = vals[2];
4191
y[0] = x[3][0];
4192
y[1] = x[3][1];
4193
y[2] = x[3][2];
4151
y[0] = vertex_coordinates[9];
4152
y[1] = vertex_coordinates[10];
4153
y[2] = vertex_coordinates[11];
4194
4154
f.evaluate(vals, y, c);
4195
4155
values[11] = vals[2];
4196
4156
}
4198
4158
/// Interpolate vertex values from dof values
4199
4159
virtual void interpolate_vertex_values(double* vertex_values,
4200
4160
const double* dof_values,
4161
const double* vertex_coordinates,
4162
int cell_orientation,
4201
4163
const ufc::cell& c) const
4202
4164
{
4203
4165
// Evaluate function and change variables
4222
4184
const double* xhat,
4223
4185
const ufc::cell& c) const
4224
4186
{
4225
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4187
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
4226
4188
}
4227
4189
4228
4190
/// Map from coordinate x in cell to coordinate xhat in reference cell
4230
4192
const double* x,
4231
4193
const ufc::cell& c) const
4232
4194
{
4233
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4195
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
4234
4196
}
4235
4197
4236
4198
/// Return the number of sub elements (for a mixed element)
4237
virtual unsigned int num_sub_elements() const
4199
virtual std::size_t num_sub_elements() const
4238
4200
{
4239
4201
return 3;
4240
4202
}
4241
4203
4242
4204
/// Create a new finite element for sub element i (for a mixed element)
4243
virtual ufc::finite_element* create_sub_element(unsigned int i) const
4205
virtual ufc::finite_element* create_sub_element(std::size_t i) const
4244
4206
{
4245
4207
switch (i)
4246
4208
{
4277
4239
4278
4240
class navierstokes_dofmap_0: public ufc::dofmap
4279
4241
{
4280
private:
4281
4282
unsigned int _global_dimension;
4283
4242
public:
4284
4243
4285
4244
/// Constructor
4286
4245
navierstokes_dofmap_0() : ufc::dofmap()
4287
4246
{
4288
_global_dimension = 0;
4247
// Do nothing
4289
4248
}
4290
4249
4291
4250
/// Destructor
4297
4256
/// Return a string identifying the dofmap
4298
4257
virtual const char* signature() const
4299
4258
{
4300
return "FFC dofmap for FiniteElement('Lagrange', Cell('tetrahedron', Space(3)), 1, None)";
4259
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, None)";
4301
4260
}
4302
4261
4303
4262
/// Return true iff mesh entities of topological dimension d are needed
4304
virtual bool needs_mesh_entities(unsigned int d) const
4263
virtual bool needs_mesh_entities(std::size_t d) const
4305
4264
{
4306
4265
switch (d)
4307
4266
{
4330
4289
return false;
4331
4290
}
4332
4291
4333
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
4334
virtual bool init_mesh(const ufc::mesh& m)
4335
{
4336
_global_dimension = m.num_entities[0];
4337
return false;
4338
}
4339
4340
/// Initialize dofmap for given cell
4341
virtual void init_cell(const ufc::mesh& m,
4342
const ufc::cell& c)
4343
{
4344
// Do nothing
4345
}
4346
4347
/// Finish initialization of dofmap for cells
4348
virtual void init_cell_finalize()
4349
{
4350
// Do nothing
4351
}
4352
4353
4292
/// Return the topological dimension of the associated cell shape
4354
virtual unsigned int topological_dimension() const
4293
virtual std::size_t topological_dimension() const
4355
4294
{
4356
4295
return 3;
4357
4296
}
4358
4297
4359
4298
/// Return the geometric dimension of the associated cell shape
4360
virtual unsigned int geometric_dimension() const
4299
virtual std::size_t geometric_dimension() const
4361
4300
{
4362
4301
return 3;
4363
4302
}
4364
4303
4365
4304
/// Return the dimension of the global finite element function space
4366
virtual unsigned int global_dimension() const
4305
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4306
num_global_entities) const
4367
4307
{
4368
return _global_dimension;
4308
return num_global_entities[0];
4369
4309
}
4370
4310
4371
4311
/// Return the dimension of the local finite element function space for a cell
4372
virtual unsigned int local_dimension(const ufc::cell& c) const
4312
virtual std::size_t local_dimension(const ufc::cell& c) const
4373
4313
{
4374
4314
return 4;
4375
4315
}
4376
4316
4377
4317
/// Return the maximum dimension of the local finite element function space
4378
virtual unsigned int max_local_dimension() const
4318
virtual std::size_t max_local_dimension() const
4379
4319
{
4380
4320
return 4;
4381
4321
}
4382
4322
4383
4323
/// Return the number of dofs on each cell facet
4384
virtual unsigned int num_facet_dofs() const
4324
virtual std::size_t num_facet_dofs() const
4385
4325
{
4386
4326
return 3;
4387
4327
}
4388
4328
4389
4329
/// Return the number of dofs associated with each cell entity of dimension d
4390
virtual unsigned int num_entity_dofs(unsigned int d) const
4330
virtual std::size_t num_entity_dofs(std::size_t d) const
4391
4331
{
4392
4332
switch (d)
4393
4333
{
4417
4357
}
4418
4358
4419
4359
/// Tabulate the local-to-global mapping of dofs on a cell
4420
virtual void tabulate_dofs(unsigned int* dofs,
4421
const ufc::mesh& m,
4360
virtual void tabulate_dofs(std::size_t* dofs,
4361
const std::vector<std::size_t>& num_global_entities,
4422
4362
const ufc::cell& c) const
4423
4363
{
4424
4364
dofs[0] = c.entity_indices[0][0];
4428
4368
}
4429
4369
4430
4370
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4431
virtual void tabulate_facet_dofs(unsigned int* dofs,
4432
unsigned int facet) const
4371
virtual void tabulate_facet_dofs(std::size_t* dofs,
4372
std::size_t facet) const
4433
4373
{
4434
4374
switch (facet)
4435
4375
{
4466
4406
}
4467
4407
4468
4408
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4469
virtual void tabulate_entity_dofs(unsigned int* dofs,
4470
unsigned int d, unsigned int i) const
4409
virtual void tabulate_entity_dofs(std::size_t* dofs,
4410
std::size_t d, std::size_t i) const
4471
4411
{
4472
4412
if (d > 3)
4473
4413
{
4529
4469
}
4530
4470
4531
4471
/// Tabulate the coordinates of all dofs on a cell
4532
virtual void tabulate_coordinates(double** coordinates,
4533
const ufc::cell& c) const
4472
virtual void tabulate_coordinates(double** dof_coordinates,
4473
const double* vertex_coordinates) const
4534
4474
{
4535
const double * const * x = c.coordinates;
4536
4537
coordinates[0][0] = x[0][0];
4538
coordinates[0][1] = x[0][1];
4539
coordinates[0][2] = x[0][2];
4540
coordinates[1][0] = x[1][0];
4541
coordinates[1][1] = x[1][1];
4542
coordinates[1][2] = x[1][2];
4543
coordinates[2][0] = x[2][0];
4544
coordinates[2][1] = x[2][1];
4545
coordinates[2][2] = x[2][2];
4546
coordinates[3][0] = x[3][0];
4547
coordinates[3][1] = x[3][1];
4548
coordinates[3][2] = x[3][2];
4475
dof_coordinates[0][0] = vertex_coordinates[0];
4476
dof_coordinates[0][1] = vertex_coordinates[1];
4477
dof_coordinates[0][2] = vertex_coordinates[2];
4478
dof_coordinates[1][0] = vertex_coordinates[3];
4479
dof_coordinates[1][1] = vertex_coordinates[4];
4480
dof_coordinates[1][2] = vertex_coordinates[5];
4481
dof_coordinates[2][0] = vertex_coordinates[6];
4482
dof_coordinates[2][1] = vertex_coordinates[7];
4483
dof_coordinates[2][2] = vertex_coordinates[8];
4484
dof_coordinates[3][0] = vertex_coordinates[9];
4485
dof_coordinates[3][1] = vertex_coordinates[10];
4486
dof_coordinates[3][2] = vertex_coordinates[11];
4549
4487
}
4550
4488
4551
4489
/// Return the number of sub dofmaps (for a mixed element)
4552
virtual unsigned int num_sub_dofmaps() const
4490
virtual std::size_t num_sub_dofmaps() const
4553
4491
{
4554
4492
return 0;
4555
4493
}
4556
4494
4557
4495
/// Create a new dofmap for sub dofmap i (for a mixed element)
4558
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
4496
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4559
4497
{
4560
4498
return 0;
4561
4499
}
4573
4511
4574
4512
class navierstokes_dofmap_1: public ufc::dofmap
4575
4513
{
4576
private:
4577
4578
unsigned int _global_dimension;
4579
4514
public:
4580
4515
4581
4516
/// Constructor
4582
4517
navierstokes_dofmap_1() : ufc::dofmap()
4583
4518
{
4584
_global_dimension = 0;
4519
// Do nothing
4585
4520
}
4586
4521
4587
4522
/// Destructor
4593
4528
/// Return a string identifying the dofmap
4594
4529
virtual const char* signature() const
4595
4530
{
4596
return "FFC dofmap for VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None)";
4531
return "FFC dofmap for VectorElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, 3, None)";
4597
4532
}
4598
4533
4599
4534
/// Return true iff mesh entities of topological dimension d are needed
4600
virtual bool needs_mesh_entities(unsigned int d) const
4535
virtual bool needs_mesh_entities(std::size_t d) const
4601
4536
{
4602
4537
switch (d)
4603
4538
{
4626
4561
return false;
4627
4562
}
4628
4563
4629
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
4630
virtual bool init_mesh(const ufc::mesh& m)
4631
{
4632
_global_dimension = 3*m.num_entities[0];
4633
return false;
4634
}
4635
4636
/// Initialize dofmap for given cell
4637
virtual void init_cell(const ufc::mesh& m,
4638
const ufc::cell& c)
4639
{
4640
// Do nothing
4641
}
4642
4643
/// Finish initialization of dofmap for cells
4644
virtual void init_cell_finalize()
4645
{
4646
// Do nothing
4647
}
4648
4649
4564
/// Return the topological dimension of the associated cell shape
4650
virtual unsigned int topological_dimension() const
4565
virtual std::size_t topological_dimension() const
4651
4566
{
4652
4567
return 3;
4653
4568
}
4654
4569
4655
4570
/// Return the geometric dimension of the associated cell shape
4656
virtual unsigned int geometric_dimension() const
4571
virtual std::size_t geometric_dimension() const
4657
4572
{
4658
4573
return 3;
4659
4574
}
4660
4575
4661
4576
/// Return the dimension of the global finite element function space
4662
virtual unsigned int global_dimension() const
4577
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4578
num_global_entities) const
4663
4579
{
4664
return _global_dimension;
4580
return 3*num_global_entities[0];
4665
4581
}
4666
4582
4667
4583
/// Return the dimension of the local finite element function space for a cell
4668
virtual unsigned int local_dimension(const ufc::cell& c) const
4584
virtual std::size_t local_dimension(const ufc::cell& c) const
4669
4585
{
4670
4586
return 12;
4671
4587
}
4672
4588
4673
4589
/// Return the maximum dimension of the local finite element function space
4674
virtual unsigned int max_local_dimension() const
4590
virtual std::size_t max_local_dimension() const
4675
4591
{
4676
4592
return 12;
4677
4593
}
4678
4594
4679
4595
/// Return the number of dofs on each cell facet
4680
virtual unsigned int num_facet_dofs() const
4596
virtual std::size_t num_facet_dofs() const
4681
4597
{
4682
4598
return 9;
4683
4599
}
4684
4600
4685
4601
/// Return the number of dofs associated with each cell entity of dimension d
4686
virtual unsigned int num_entity_dofs(unsigned int d) const
4602
virtual std::size_t num_entity_dofs(std::size_t d) const
4687
4603
{
4688
4604
switch (d)
4689
4605
{
4713
4629
}
4714
4630
4715
4631
/// Tabulate the local-to-global mapping of dofs on a cell
4716
virtual void tabulate_dofs(unsigned int* dofs,
4717
const ufc::mesh& m,
4632
virtual void tabulate_dofs(std::size_t* dofs,
4633
const std::vector<std::size_t>& num_global_entities,
4718
4634
const ufc::cell& c) const
4719
4635
{
4720
4636
unsigned int offset = 0;
4722
4638
dofs[1] = offset + c.entity_indices[0][1];
4723
4639
dofs[2] = offset + c.entity_indices[0][2];
4724
4640
dofs[3] = offset + c.entity_indices[0][3];
4725
offset += m.num_entities[0];
4641
offset += num_global_entities[0];
4726
4642
dofs[4] = offset + c.entity_indices[0][0];
4727
4643
dofs[5] = offset + c.entity_indices[0][1];
4728
4644
dofs[6] = offset + c.entity_indices[0][2];
4729
4645
dofs[7] = offset + c.entity_indices[0][3];
4730
offset += m.num_entities[0];
4646
offset += num_global_entities[0];
4731
4647
dofs[8] = offset + c.entity_indices[0][0];
4732
4648
dofs[9] = offset + c.entity_indices[0][1];
4733
4649
dofs[10] = offset + c.entity_indices[0][2];
4734
4650
dofs[11] = offset + c.entity_indices[0][3];
4735
offset += m.num_entities[0];
4651
offset += num_global_entities[0];
4736
4652
}
4737
4653
4738
4654
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4739
virtual void tabulate_facet_dofs(unsigned int* dofs,
4740
unsigned int facet) const
4655
virtual void tabulate_facet_dofs(std::size_t* dofs,
4656
std::size_t facet) const
4741
4657
{
4742
4658
switch (facet)
4743
4659
{
4798
4714
}
4799
4715
4800
4716
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4801
virtual void tabulate_entity_dofs(unsigned int* dofs,
4802
unsigned int d, unsigned int i) const
4717
virtual void tabulate_entity_dofs(std::size_t* dofs,
4718
std::size_t d, std::size_t i) const
4803
4719
{
4804
4720
if (d > 3)
4805
4721
{
4869
4785
}
4870
4786
4871
4787
/// Tabulate the coordinates of all dofs on a cell
4872
virtual void tabulate_coordinates(double** coordinates,
4873
const ufc::cell& c) const
4788
virtual void tabulate_coordinates(double** dof_coordinates,
4789
const double* vertex_coordinates) const
4874
4790
{
4875
const double * const * x = c.coordinates;
4876
4877
coordinates[0][0] = x[0][0];
4878
coordinates[0][1] = x[0][1];
4879
coordinates[0][2] = x[0][2];
4880
coordinates[1][0] = x[1][0];
4881
coordinates[1][1] = x[1][1];
4882
coordinates[1][2] = x[1][2];
4883
coordinates[2][0] = x[2][0];
4884
coordinates[2][1] = x[2][1];
4885
coordinates[2][2] = x[2][2];
4886
coordinates[3][0] = x[3][0];
4887
coordinates[3][1] = x[3][1];
4888
coordinates[3][2] = x[3][2];
4889
coordinates[4][0] = x[0][0];
4890
coordinates[4][1] = x[0][1];
4891
coordinates[4][2] = x[0][2];
4892
coordinates[5][0] = x[1][0];
4893
coordinates[5][1] = x[1][1];
4894
coordinates[5][2] = x[1][2];
4895
coordinates[6][0] = x[2][0];
4896
coordinates[6][1] = x[2][1];
4897
coordinates[6][2] = x[2][2];
4898
coordinates[7][0] = x[3][0];
4899
coordinates[7][1] = x[3][1];
4900
coordinates[7][2] = x[3][2];
4901
coordinates[8][0] = x[0][0];
4902
coordinates[8][1] = x[0][1];
4903
coordinates[8][2] = x[0][2];
4904
coordinates[9][0] = x[1][0];
4905
coordinates[9][1] = x[1][1];
4906
coordinates[9][2] = x[1][2];
4907
coordinates[10][0] = x[2][0];
4908
coordinates[10][1] = x[2][1];
4909
coordinates[10][2] = x[2][2];
4910
coordinates[11][0] = x[3][0];
4911
coordinates[11][1] = x[3][1];
4912
coordinates[11][2] = x[3][2];
4791
dof_coordinates[0][0] = vertex_coordinates[0];
4792
dof_coordinates[0][1] = vertex_coordinates[1];
4793
dof_coordinates[0][2] = vertex_coordinates[2];
4794
dof_coordinates[1][0] = vertex_coordinates[3];
4795
dof_coordinates[1][1] = vertex_coordinates[4];
4796
dof_coordinates[1][2] = vertex_coordinates[5];
4797
dof_coordinates[2][0] = vertex_coordinates[6];
4798
dof_coordinates[2][1] = vertex_coordinates[7];
4799
dof_coordinates[2][2] = vertex_coordinates[8];
4800
dof_coordinates[3][0] = vertex_coordinates[9];
4801
dof_coordinates[3][1] = vertex_coordinates[10];
4802
dof_coordinates[3][2] = vertex_coordinates[11];
4803
dof_coordinates[4][0] = vertex_coordinates[0];
4804
dof_coordinates[4][1] = vertex_coordinates[1];
4805
dof_coordinates[4][2] = vertex_coordinates[2];
4806
dof_coordinates[5][0] = vertex_coordinates[3];
4807
dof_coordinates[5][1] = vertex_coordinates[4];
4808
dof_coordinates[5][2] = vertex_coordinates[5];
4809
dof_coordinates[6][0] = vertex_coordinates[6];
4810
dof_coordinates[6][1] = vertex_coordinates[7];
4811
dof_coordinates[6][2] = vertex_coordinates[8];
4812
dof_coordinates[7][0] = vertex_coordinates[9];
4813
dof_coordinates[7][1] = vertex_coordinates[10];
4814
dof_coordinates[7][2] = vertex_coordinates[11];
4815
dof_coordinates[8][0] = vertex_coordinates[0];
4816
dof_coordinates[8][1] = vertex_coordinates[1];
4817
dof_coordinates[8][2] = vertex_coordinates[2];
4818
dof_coordinates[9][0] = vertex_coordinates[3];
4819
dof_coordinates[9][1] = vertex_coordinates[4];
4820
dof_coordinates[9][2] = vertex_coordinates[5];
4821
dof_coordinates[10][0] = vertex_coordinates[6];
4822
dof_coordinates[10][1] = vertex_coordinates[7];
4823
dof_coordinates[10][2] = vertex_coordinates[8];
4824
dof_coordinates[11][0] = vertex_coordinates[9];
4825
dof_coordinates[11][1] = vertex_coordinates[10];
4826
dof_coordinates[11][2] = vertex_coordinates[11];
4913
4827
}
4914
4828
4915
4829
/// Return the number of sub dofmaps (for a mixed element)
4916
virtual unsigned int num_sub_dofmaps() const
4830
virtual std::size_t num_sub_dofmaps() const
4917
4831
{
4918
4832
return 3;
4919
4833
}
4920
4834
4921
4835
/// Create a new dofmap for sub dofmap i (for a mixed element)
4922
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
4836
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4923
4837
{
4924
4838
switch (i)
4925
4839
{
4955
4869
/// tensor corresponding to the local contribution to a form from
4956
4870
/// the integral over a cell.
4957
4871
4958
class navierstokes_cell_integral_0_0: public ufc::cell_integral
4872
class navierstokes_cell_integral_0_otherwise: public ufc::cell_integral
4959
4873
{
4960
4874
public:
4961
4875
4962
4876
/// Constructor
4963
navierstokes_cell_integral_0_0() : ufc::cell_integral()
4877
navierstokes_cell_integral_0_otherwise() : ufc::cell_integral()
4964
4878
{
4965
4879
// Do nothing
4966
4880
}
4967
4881
4968
4882
/// Destructor
4969
virtual ~navierstokes_cell_integral_0_0()
4883
virtual ~navierstokes_cell_integral_0_otherwise()
4970
4884
{
4971
4885
// Do nothing
4972
4886
}
4974
4888
/// Tabulate the tensor for the contribution from a local cell
4975
4889
virtual void tabulate_tensor(double* A,
4976
4890
const double * const * w,
4977
const ufc::cell& c) const
4891
const double* vertex_coordinates,
4892
int cell_orientation) const
4978
4893
{
4979
// Extract vertex coordinates
4980
const double * const * x = c.coordinates;
4981
4982
// Compute Jacobian of affine map from reference cell
4983
const double J_00 = x[1][0] - x[0][0];
4984
const double J_01 = x[2][0] - x[0][0];
4985
const double J_02 = x[3][0] - x[0][0];
4986
const double J_10 = x[1][1] - x[0][1];
4987
const double J_11 = x[2][1] - x[0][1];
4988
const double J_12 = x[3][1] - x[0][1];
4989
const double J_20 = x[1][2] - x[0][2];
4990
const double J_21 = x[2][2] - x[0][2];
4991
const double J_22 = x[3][2] - x[0][2];
4992
4993
// Compute sub determinants
4994
const double d_00 = J_11*J_22 - J_12*J_21;
4995
const double d_01 = J_12*J_20 - J_10*J_22;
4996
const double d_02 = J_10*J_21 - J_11*J_20;
4997
const double d_10 = J_02*J_21 - J_01*J_22;
4998
const double d_11 = J_00*J_22 - J_02*J_20;
4999
const double d_12 = J_01*J_20 - J_00*J_21;
5000
const double d_20 = J_01*J_12 - J_02*J_11;
5001
const double d_21 = J_02*J_10 - J_00*J_12;
5002
const double d_22 = J_00*J_11 - J_01*J_10;
5003
5004
// Compute determinant of Jacobian
5005
double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
5006
5007
// Compute inverse of Jacobian
5008
const double K_00 = d_00 / detJ;
5009
const double K_01 = d_10 / detJ;
5010
const double K_02 = d_20 / detJ;
5011
const double K_10 = d_01 / detJ;
5012
const double K_11 = d_11 / detJ;
5013
const double K_12 = d_21 / detJ;
5014
const double K_20 = d_02 / detJ;
5015
const double K_21 = d_12 / detJ;
5016
const double K_22 = d_22 / detJ;
4894
// Compute Jacobian
4895
double J[9];
4896
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
4897
4898
// Compute Jacobian inverse and determinant
4899
double K[9];
4900
double detJ;
4901
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
5017
4902
5018
4903
// Set scale factor
5019
4904
const double det = std::abs(detJ);
5020
4905
5021
// Cell Volume.
5022
5023
// Compute circumradius.
5024
5025
5026
// Facet Area (divide by two because 'det' is scaled by area of reference triangle).
4906
// Cell volume
4907
4908
// Compute circumradius
4909
4910
4911
// Facet area (divide by two because 'det' is scaled by area of reference triangle)
5027
4912
5028
4913
// Array of quadrature weights.
5029
4914
static const double W4[4] = {0.041666667, 0.041666667, 0.041666667, 0.041666667};
5135
5020
for (unsigned int k = 0; k < 12; k++)
5136
5021
{
5137
5022
// Number of operations to compute entry: 68
5138
A[j*12 + k] += (((((K_01*FE0_C2_D100[ip][k] + K_11*FE0_C2_D010[ip][k] + K_21*FE0_C2_D001[ip][k]))*F1 + ((K_02*FE0_C2_D100[ip][k] + K_12*FE0_C2_D010[ip][k] + K_22*FE0_C2_D001[ip][k]))*F2 + ((K_00*FE0_C2_D100[ip][k] + K_10*FE0_C2_D010[ip][k] + K_20*FE0_C2_D001[ip][k]))*F0))*FE0_C2[ip][j] + ((((K_01*FE0_C0_D100[ip][k] + K_11*FE0_C0_D010[ip][k] + K_21*FE0_C0_D001[ip][k]))*F1 + ((K_02*FE0_C0_D100[ip][k] + K_12*FE0_C0_D010[ip][k] + K_22*FE0_C0_D001[ip][k]))*F2 + ((K_00*FE0_C0_D100[ip][k] + K_10*FE0_C0_D010[ip][k] + K_20*FE0_C0_D001[ip][k]))*F0))*FE0_C0[ip][j] + ((((K_00*FE0_C1_D100[ip][k] + K_10*FE0_C1_D010[ip][k] + K_20*FE0_C1_D001[ip][k]))*F0 + ((K_01*FE0_C1_D100[ip][k] + K_11*FE0_C1_D010[ip][k] + K_21*FE0_C1_D001[ip][k]))*F1 + ((K_02*FE0_C1_D100[ip][k] + K_12*FE0_C1_D010[ip][k] + K_22*FE0_C1_D001[ip][k]))*F2))*FE0_C1[ip][j])*W4[ip]*det;
5023
A[j*12 + k] += (((((K[0]*FE0_C0_D100[ip][k] + K[3]*FE0_C0_D010[ip][k] + K[6]*FE0_C0_D001[ip][k]))*F0 + ((K[1]*FE0_C0_D100[ip][k] + K[4]*FE0_C0_D010[ip][k] + K[7]*FE0_C0_D001[ip][k]))*F1 + ((K[2]*FE0_C0_D100[ip][k] + K[5]*FE0_C0_D010[ip][k] + K[8]*FE0_C0_D001[ip][k]))*F2))*FE0_C0[ip][j] + ((((K[1]*FE0_C2_D100[ip][k] + K[4]*FE0_C2_D010[ip][k] + K[7]*FE0_C2_D001[ip][k]))*F1 + ((K[0]*FE0_C2_D100[ip][k] + K[3]*FE0_C2_D010[ip][k] + K[6]*FE0_C2_D001[ip][k]))*F0 + ((K[2]*FE0_C2_D100[ip][k] + K[5]*FE0_C2_D010[ip][k] + K[8]*FE0_C2_D001[ip][k]))*F2))*FE0_C2[ip][j] + ((((K[0]*FE0_C1_D100[ip][k] + K[3]*FE0_C1_D010[ip][k] + K[6]*FE0_C1_D001[ip][k]))*F0 + ((K[1]*FE0_C1_D100[ip][k] + K[4]*FE0_C1_D010[ip][k] + K[7]*FE0_C1_D001[ip][k]))*F1 + ((K[2]*FE0_C1_D100[ip][k] + K[5]*FE0_C1_D010[ip][k] + K[8]*FE0_C1_D001[ip][k]))*F2))*FE0_C1[ip][j])*W4[ip]*det;
5139
5024
}// end loop over 'k'
5140
5025
}// end loop over 'j'
5141
5026
}// end loop over 'ip'
5142
5027
}
5143
5028
5144
/// Tabulate the tensor for the contribution from a local cell
5145
/// using the specified reference cell quadrature points/weights
5146
virtual void tabulate_tensor(double* A,
5147
const double * const * w,
5148
const ufc::cell& c,
5149
unsigned int num_quadrature_points,
5150
const double * const * quadrature_points,
5151
const double* quadrature_weights) const
5152
{
5153
std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
5154
}
5155
5156
5029
};
5157
5030
5158
5031
/// This class defines the interface for the assembly of the global
5189
5062
/// Return a string identifying the form
5190
5063
virtual const char* signature() const
5191
5064
{
5192
return "Form([Integral(IndexSum(Product(IndexSum(Product(Indexed(Coefficient(VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None), 0), MultiIndex((Index(0),), {Index(0): 3})), Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None), 1), MultiIndex((Index(0),), {Index(0): 3})), MultiIndex((Index(1),), {Index(1): 3}))), MultiIndex((Index(0),), {Index(0): 3})), Indexed(Argument(VectorElement('Lagrange', Cell('tetrahedron', Space(3)), 1, 3, None), 0), MultiIndex((Index(1),), {Index(1): 3}))), MultiIndex((Index(1),), {Index(1): 3})), Measure('cell', 0, None))])";
5065
return "562aab673b326bf1528829b2d9fd338b3aa2c324573441a5499c132135ac2cf3a28713f77266e32ebf3700a5ffdb8d4502f13c1ac5e7213db4dc525ac9602e4d";
5193
5066
}
5194
5067
5195
5068
/// Return the rank of the global tensor (r)
5196
virtual unsigned int rank() const
5069
virtual std::size_t rank() const
5197
5070
{
5198
5071
return 2;
5199
5072
}
5200
5073
5201
5074
/// Return the number of coefficients (n)
5202
virtual unsigned int num_coefficients() const
5075
virtual std::size_t num_coefficients() const
5203
5076
{
5204
5077
return 1;
5205
5078
}
5206
5079
5207
5080
/// Return the number of cell domains
5208
virtual unsigned int num_cell_domains() const
5081
virtual std::size_t num_cell_domains() const
5209
5082
{
5210
return 1;
5083
return 0;
5211
5084
}
5212
5085
5213
5086
/// Return the number of exterior facet domains
5214
virtual unsigned int num_exterior_facet_domains() const
5087
virtual std::size_t num_exterior_facet_domains() const
5215
5088
{
5216
5089
return 0;
5217
5090
}
5218
5091
5219
5092
/// Return the number of interior facet domains
5220
virtual unsigned int num_interior_facet_domains() const
5221
{
5222
return 0;
5093
virtual std::size_t num_interior_facet_domains() const
5094
{
5095
return 0;
5096
}
5097
5098
/// Return the number of point domains
5099
virtual std::size_t num_point_domains() const
5100
{
5101
return 0;
5102
}
5103
5104
/// Return whether the form has any cell integrals
5105
virtual bool has_cell_integrals() const
5106
{
5107
return true;
5108
}
5109
5110
/// Return whether the form has any exterior facet integrals
5111
virtual bool has_exterior_facet_integrals() const
5112
{
5113
return false;
5114
}
5115
5116
/// Return whether the form has any interior facet integrals
5117
virtual bool has_interior_facet_integrals() const
5118
{
5119
return false;
5120
}
5121
5122
/// Return whether the form has any point integrals
5123
virtual bool has_point_integrals() const
5124
{
5125
return false;
5223
5126
}
5224
5127
5225
5128
/// Create a new finite element for argument function i
5226
virtual ufc::finite_element* create_finite_element(unsigned int i) const
5129
virtual ufc::finite_element* create_finite_element(std::size_t i) const
5227
5130
{
5228
5131
switch (i)
5229
5132
{
5248
5151
}
5249
5152
5250
5153
/// Create a new dofmap for argument function i
5251
virtual ufc::dofmap* create_dofmap(unsigned int i) const
5154
virtual ufc::dofmap* create_dofmap(std::size_t i) const
5252
5155
{
5253
5156
switch (i)
5254
5157
{
5273
5176
}
5274
5177
5275
5178
/// Create a new cell integral on sub domain i
5276
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
5179
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
5277
5180
{
5278
switch (i)
5279
{
5280
case 0:
5281
{
5282
return new navierstokes_cell_integral_0_0();
5283
break;
5284
}
5285
}
5286
5287
5181
return 0;
5288
5182
}
5289
5183
5290
5184
/// Create a new exterior facet integral on sub domain i
5291
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
5185
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
5292
5186
{
5293
5187
return 0;
5294
5188
}
5295
5189
5296
5190
/// Create a new interior facet integral on sub domain i
5297
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
5191
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
5192
{
5193
return 0;
5194
}
5195
5196
/// Create a new point integral on sub domain i
5197
virtual ufc::point_integral* create_point_integral(std::size_t i) const
5198
{
5199
return 0;
5200
}
5201
5202
/// Create a new cell integral on everywhere else
5203
virtual ufc::cell_integral* create_default_cell_integral() const
5204
{
5205
return new navierstokes_cell_integral_0_otherwise();
5206
}
5207
5208
/// Create a new exterior facet integral on everywhere else
5209
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
5210
{
5211
return 0;
5212
}
5213
5214
/// Create a new interior facet integral on everywhere else
5215
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
5216
{
5217
return 0;
5218
}
5219
5220
/// Create a new point integral on everywhere else
5221
virtual ufc::point_integral* create_default_point_integral() const
5298
5222
{
5299
5223
return 0;
5300
5224
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1