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