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