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