1
// This code conforms with the UFC specification version 2.1.0+
2
// and was automatically generated by FFC version 1.1.0+.
4
// This code was generated with the option '-l dolfin' and
5
// contains DOLFIN-specific wrappers that depend on DOLFIN.
7
// This code was generated with the following parameters:
10
// convert_exceptions_to_warnings: False
11
// cpp_optimize: False
12
// cpp_optimize_flags: '-O2'
14
// error_control: False
23
// quadrature_degree: 'auto'
24
// quadrature_rule: 'auto'
25
// representation: 'auto'
27
// swig_binary: 'swig'
30
#ifndef __NAVIERSTOKES_H
31
#define __NAVIERSTOKES_H
38
/// This class defines the interface for a finite element.
40
class navierstokes_finite_element_0: public ufc::finite_element
45
navierstokes_finite_element_0() : ufc::finite_element()
51
virtual ~navierstokes_finite_element_0()
56
/// Return a string identifying the finite element
57
virtual const char* signature() const
59
return "FiniteElement('Discontinuous Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 0, None)";
62
/// Return the cell shape
63
virtual ufc::shape cell_shape() const
65
return ufc::tetrahedron;
68
/// Return the topological dimension of the cell shape
69
virtual std::size_t topological_dimension() const
74
/// Return the geometric dimension of the cell shape
75
virtual std::size_t geometric_dimension() const
80
/// Return the dimension of the finite element function space
81
virtual std::size_t space_dimension() const
86
/// Return the rank of the value space
87
virtual std::size_t value_rank() const
92
/// Return the dimension of the value space for axis i
93
virtual std::size_t value_dimension(std::size_t i) const
98
/// Evaluate basis function i at given point x in cell
99
virtual void evaluate_basis(std::size_t i,
102
const double* vertex_coordinates,
103
int cell_orientation) const
107
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
109
// Compute Jacobian inverse and determinant
112
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
117
// Compute subdeterminants
119
// Get coordinates and map to the reference (FIAT) element
125
// Array of basisvalues
126
double basisvalues[1] = {0.0};
128
// Declare helper variables
130
// Compute basisvalues
131
basisvalues[0] = 1.0;
133
// Table(s) of coefficients
134
static const double coefficients0[1] = \
138
for (unsigned int r = 0; r < 1; r++)
140
*values += coefficients0[r]*basisvalues[r];
141
}// end loop over 'r'
144
/// Evaluate all basis functions at given point x in cell
145
virtual void evaluate_basis_all(double* values,
147
const double* vertex_coordinates,
148
int cell_orientation) const
150
// Element is constant, calling evaluate_basis.
151
evaluate_basis(0, values, x, vertex_coordinates, cell_orientation);
154
/// Evaluate order n derivatives of basis function i at given point x in cell
155
virtual void evaluate_basis_derivatives(std::size_t i,
159
const double* vertex_coordinates,
160
int cell_orientation) const
164
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
166
// Compute Jacobian inverse and determinant
169
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
174
// Compute subdeterminants
176
// Get coordinates and map to the reference (FIAT) element
179
// Compute number of derivatives.
180
unsigned int num_derivatives = 1;
181
for (unsigned int r = 0; r < n; r++)
183
num_derivatives *= 3;
184
}// end loop over 'r'
186
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
187
unsigned int **combinations = new unsigned int *[num_derivatives];
188
for (unsigned int row = 0; row < num_derivatives; row++)
190
combinations[row] = new unsigned int [n];
191
for (unsigned int col = 0; col < n; col++)
192
combinations[row][col] = 0;
195
// Generate combinations of derivatives
196
for (unsigned int row = 1; row < num_derivatives; row++)
198
for (unsigned int num = 0; num < row; num++)
200
for (unsigned int col = n-1; col+1 > 0; col--)
202
if (combinations[row][col] + 1 > 2)
203
combinations[row][col] = 0;
206
combinations[row][col] += 1;
213
// Compute inverse of Jacobian
214
const double Jinv[3][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}, {K[6], K[7], K[8]}};
216
// Declare transformation matrix
217
// Declare pointer to two dimensional array and initialise
218
double **transform = new double *[num_derivatives];
220
for (unsigned int j = 0; j < num_derivatives; j++)
222
transform[j] = new double [num_derivatives];
223
for (unsigned int k = 0; k < num_derivatives; k++)
227
// Construct transformation matrix
228
for (unsigned int row = 0; row < num_derivatives; row++)
230
for (unsigned int col = 0; col < num_derivatives; col++)
232
for (unsigned int k = 0; k < n; k++)
233
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
237
// Reset values. Assuming that values is always an array.
238
for (unsigned int r = 0; r < num_derivatives; r++)
241
}// end loop over 'r'
244
// Array of basisvalues
245
double basisvalues[1] = {0.0};
247
// Declare helper variables
249
// Compute basisvalues
250
basisvalues[0] = 1.0;
252
// Table(s) of coefficients
253
static const double coefficients0[1] = \
256
// Tables of derivatives of the polynomial base (transpose).
257
static const double dmats0[1][1] = \
260
static const double dmats1[1][1] = \
263
static const double dmats2[1][1] = \
266
// Compute reference derivatives.
267
// Declare pointer to array of derivatives on FIAT element.
268
double *derivatives = new double[num_derivatives];
269
for (unsigned int r = 0; r < num_derivatives; r++)
271
derivatives[r] = 0.0;
272
}// end loop over 'r'
274
// Declare derivative matrix (of polynomial basis).
275
double dmats[1][1] = \
278
// Declare (auxiliary) derivative matrix (of polynomial basis).
279
double dmats_old[1][1] = \
282
// Loop possible derivatives.
283
for (unsigned int r = 0; r < num_derivatives; r++)
285
// Resetting dmats values to compute next derivative.
286
for (unsigned int t = 0; t < 1; t++)
288
for (unsigned int u = 0; u < 1; u++)
296
}// end loop over 'u'
297
}// end loop over 't'
299
// Looping derivative order to generate dmats.
300
for (unsigned int s = 0; s < n; s++)
302
// Updating dmats_old with new values and resetting dmats.
303
for (unsigned int t = 0; t < 1; t++)
305
for (unsigned int u = 0; u < 1; u++)
307
dmats_old[t][u] = dmats[t][u];
309
}// end loop over 'u'
310
}// end loop over 't'
312
// Update dmats using an inner product.
313
if (combinations[r][s] == 0)
315
for (unsigned int t = 0; t < 1; t++)
317
for (unsigned int u = 0; u < 1; u++)
319
for (unsigned int tu = 0; tu < 1; tu++)
321
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
322
}// end loop over 'tu'
323
}// end loop over 'u'
324
}// end loop over 't'
327
if (combinations[r][s] == 1)
329
for (unsigned int t = 0; t < 1; t++)
331
for (unsigned int u = 0; u < 1; u++)
333
for (unsigned int tu = 0; tu < 1; tu++)
335
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
336
}// end loop over 'tu'
337
}// end loop over 'u'
338
}// end loop over 't'
341
if (combinations[r][s] == 2)
343
for (unsigned int t = 0; t < 1; t++)
345
for (unsigned int u = 0; u < 1; u++)
347
for (unsigned int tu = 0; tu < 1; tu++)
349
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
350
}// end loop over 'tu'
351
}// end loop over 'u'
352
}// end loop over 't'
355
}// end loop over 's'
356
for (unsigned int s = 0; s < 1; s++)
358
for (unsigned int t = 0; t < 1; t++)
360
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
361
}// end loop over 't'
362
}// end loop over 's'
363
}// end loop over 'r'
365
// Transform derivatives back to physical element
366
for (unsigned int r = 0; r < num_derivatives; r++)
368
for (unsigned int s = 0; s < num_derivatives; s++)
370
values[r] += transform[r][s]*derivatives[s];
371
}// end loop over 's'
372
}// end loop over 'r'
374
// Delete pointer to array of derivatives on FIAT element
375
delete [] derivatives;
377
// Delete pointer to array of combinations of derivatives and transform
378
for (unsigned int r = 0; r < num_derivatives; r++)
380
delete [] combinations[r];
381
}// end loop over 'r'
382
delete [] combinations;
383
for (unsigned int r = 0; r < num_derivatives; r++)
385
delete [] transform[r];
386
}// end loop over 'r'
390
/// Evaluate order n derivatives of all basis functions at given point x in cell
391
virtual void evaluate_basis_derivatives_all(std::size_t n,
394
const double* vertex_coordinates,
395
int cell_orientation) const
397
// Element is constant, calling evaluate_basis_derivatives.
398
evaluate_basis_derivatives(0, n, values, x, vertex_coordinates, cell_orientation);
401
/// Evaluate linear functional for dof i on the function f
402
virtual double evaluate_dof(std::size_t i,
403
const ufc::function& f,
404
const double* vertex_coordinates,
405
int cell_orientation,
406
const ufc::cell& c) const
408
// Declare variables for result of evaluation
411
// Declare variable for physical coordinates
417
y[0] = 0.25*vertex_coordinates[0] + 0.25*vertex_coordinates[3] + 0.25*vertex_coordinates[6] + 0.25*vertex_coordinates[9];
418
y[1] = 0.25*vertex_coordinates[1] + 0.25*vertex_coordinates[4] + 0.25*vertex_coordinates[7] + 0.25*vertex_coordinates[10];
419
y[2] = 0.25*vertex_coordinates[2] + 0.25*vertex_coordinates[5] + 0.25*vertex_coordinates[8] + 0.25*vertex_coordinates[11];
420
f.evaluate(vals, y, c);
429
/// Evaluate linear functionals for all dofs on the function f
430
virtual void evaluate_dofs(double* values,
431
const ufc::function& f,
432
const double* vertex_coordinates,
433
int cell_orientation,
434
const ufc::cell& c) const
436
// Declare variables for result of evaluation
439
// Declare variable for physical coordinates
441
y[0] = 0.25*vertex_coordinates[0] + 0.25*vertex_coordinates[3] + 0.25*vertex_coordinates[6] + 0.25*vertex_coordinates[9];
442
y[1] = 0.25*vertex_coordinates[1] + 0.25*vertex_coordinates[4] + 0.25*vertex_coordinates[7] + 0.25*vertex_coordinates[10];
443
y[2] = 0.25*vertex_coordinates[2] + 0.25*vertex_coordinates[5] + 0.25*vertex_coordinates[8] + 0.25*vertex_coordinates[11];
444
f.evaluate(vals, y, c);
448
/// Interpolate vertex values from dof values
449
virtual void interpolate_vertex_values(double* vertex_values,
450
const double* dof_values,
451
const double* vertex_coordinates,
452
int cell_orientation,
453
const ufc::cell& c) const
455
// Evaluate function and change variables
456
vertex_values[0] = dof_values[0];
457
vertex_values[1] = dof_values[0];
458
vertex_values[2] = dof_values[0];
459
vertex_values[3] = dof_values[0];
462
/// Map coordinate xhat from reference cell to coordinate x in cell
463
virtual void map_from_reference_cell(double* x,
465
const ufc::cell& c) const
467
throw std::runtime_error("map_from_reference_cell not yet implemented.");
470
/// Map from coordinate x in cell to coordinate xhat in reference cell
471
virtual void map_to_reference_cell(double* xhat,
473
const ufc::cell& c) const
475
throw std::runtime_error("map_to_reference_cell not yet implemented.");
478
/// Return the number of sub elements (for a mixed element)
479
virtual std::size_t num_sub_elements() const
484
/// Create a new finite element for sub element i (for a mixed element)
485
virtual ufc::finite_element* create_sub_element(std::size_t i) const
490
/// Create a new class instance
491
virtual ufc::finite_element* create() const
493
return new navierstokes_finite_element_0();
498
/// This class defines the interface for a finite element.
500
class navierstokes_finite_element_1: public ufc::finite_element
505
navierstokes_finite_element_1() : ufc::finite_element()
511
virtual ~navierstokes_finite_element_1()
516
/// Return a string identifying the finite element
517
virtual const char* signature() const
519
return "FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, None)";
522
/// Return the cell shape
523
virtual ufc::shape cell_shape() const
525
return ufc::tetrahedron;
528
/// Return the topological dimension of the cell shape
529
virtual std::size_t topological_dimension() const
534
/// Return the geometric dimension of the cell shape
535
virtual std::size_t geometric_dimension() const
540
/// Return the dimension of the finite element function space
541
virtual std::size_t space_dimension() const
546
/// Return the rank of the value space
547
virtual std::size_t value_rank() const
552
/// Return the dimension of the value space for axis i
553
virtual std::size_t value_dimension(std::size_t i) const
558
/// Evaluate basis function i at given point x in cell
559
virtual void evaluate_basis(std::size_t i,
562
const double* vertex_coordinates,
563
int cell_orientation) const
567
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
569
// Compute Jacobian inverse and determinant
572
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
576
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
577
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
578
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
580
// Compute subdeterminants
581
const double d_00 = J[4]*J[8] - J[5]*J[7];
582
const double d_01 = J[5]*J[6] - J[3]*J[8];
583
const double d_02 = J[3]*J[7] - J[4]*J[6];
584
const double d_10 = J[2]*J[7] - J[1]*J[8];
585
const double d_11 = J[0]*J[8] - J[2]*J[6];
586
const double d_12 = J[1]*J[6] - J[0]*J[7];
587
const double d_20 = J[1]*J[5] - J[2]*J[4];
588
const double d_21 = J[2]*J[3] - J[0]*J[5];
589
const double d_22 = J[0]*J[4] - J[1]*J[3];
591
// Get coordinates and map to the reference (FIAT) element
592
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
593
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
594
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
604
// Array of basisvalues
605
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
607
// Declare helper variables
608
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
610
// Compute basisvalues
611
basisvalues[0] = 1.0;
612
basisvalues[1] = tmp0;
613
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
614
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
615
basisvalues[0] *= std::sqrt(0.75);
616
basisvalues[3] *= std::sqrt(1.25);
617
basisvalues[2] *= std::sqrt(2.5);
618
basisvalues[1] *= std::sqrt(7.5);
620
// Table(s) of coefficients
621
static const double coefficients0[4] = \
622
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
625
for (unsigned int r = 0; r < 4; r++)
627
*values += coefficients0[r]*basisvalues[r];
628
}// end loop over 'r'
634
// Array of basisvalues
635
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
637
// Declare helper variables
638
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
640
// Compute basisvalues
641
basisvalues[0] = 1.0;
642
basisvalues[1] = tmp0;
643
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
644
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
645
basisvalues[0] *= std::sqrt(0.75);
646
basisvalues[3] *= std::sqrt(1.25);
647
basisvalues[2] *= std::sqrt(2.5);
648
basisvalues[1] *= std::sqrt(7.5);
650
// Table(s) of coefficients
651
static const double coefficients0[4] = \
652
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
655
for (unsigned int r = 0; r < 4; r++)
657
*values += coefficients0[r]*basisvalues[r];
658
}// end loop over 'r'
664
// Array of basisvalues
665
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
667
// Declare helper variables
668
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
670
// Compute basisvalues
671
basisvalues[0] = 1.0;
672
basisvalues[1] = tmp0;
673
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
674
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
675
basisvalues[0] *= std::sqrt(0.75);
676
basisvalues[3] *= std::sqrt(1.25);
677
basisvalues[2] *= std::sqrt(2.5);
678
basisvalues[1] *= std::sqrt(7.5);
680
// Table(s) of coefficients
681
static const double coefficients0[4] = \
682
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
685
for (unsigned int r = 0; r < 4; r++)
687
*values += coefficients0[r]*basisvalues[r];
688
}// end loop over 'r'
694
// Array of basisvalues
695
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
697
// Declare helper variables
698
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
700
// Compute basisvalues
701
basisvalues[0] = 1.0;
702
basisvalues[1] = tmp0;
703
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
704
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
705
basisvalues[0] *= std::sqrt(0.75);
706
basisvalues[3] *= std::sqrt(1.25);
707
basisvalues[2] *= std::sqrt(2.5);
708
basisvalues[1] *= std::sqrt(7.5);
710
// Table(s) of coefficients
711
static const double coefficients0[4] = \
712
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
715
for (unsigned int r = 0; r < 4; r++)
717
*values += coefficients0[r]*basisvalues[r];
718
}// end loop over 'r'
725
/// Evaluate all basis functions at given point x in cell
726
virtual void evaluate_basis_all(double* values,
728
const double* vertex_coordinates,
729
int cell_orientation) const
731
// Helper variable to hold values of a single dof.
732
double dof_values = 0.0;
734
// Loop dofs and call evaluate_basis
735
for (unsigned int r = 0; r < 4; r++)
737
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
738
values[r] = dof_values;
739
}// end loop over 'r'
742
/// Evaluate order n derivatives of basis function i at given point x in cell
743
virtual void evaluate_basis_derivatives(std::size_t i,
747
const double* vertex_coordinates,
748
int cell_orientation) const
752
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
754
// Compute Jacobian inverse and determinant
757
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
761
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
762
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
763
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
765
// Compute subdeterminants
766
const double d_00 = J[4]*J[8] - J[5]*J[7];
767
const double d_01 = J[5]*J[6] - J[3]*J[8];
768
const double d_02 = J[3]*J[7] - J[4]*J[6];
769
const double d_10 = J[2]*J[7] - J[1]*J[8];
770
const double d_11 = J[0]*J[8] - J[2]*J[6];
771
const double d_12 = J[1]*J[6] - J[0]*J[7];
772
const double d_20 = J[1]*J[5] - J[2]*J[4];
773
const double d_21 = J[2]*J[3] - J[0]*J[5];
774
const double d_22 = J[0]*J[4] - J[1]*J[3];
776
// Get coordinates and map to the reference (FIAT) element
777
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
778
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
779
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
782
// Compute number of derivatives.
783
unsigned int num_derivatives = 1;
784
for (unsigned int r = 0; r < n; r++)
786
num_derivatives *= 3;
787
}// end loop over 'r'
789
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
790
unsigned int **combinations = new unsigned int *[num_derivatives];
791
for (unsigned int row = 0; row < num_derivatives; row++)
793
combinations[row] = new unsigned int [n];
794
for (unsigned int col = 0; col < n; col++)
795
combinations[row][col] = 0;
798
// Generate combinations of derivatives
799
for (unsigned int row = 1; row < num_derivatives; row++)
801
for (unsigned int num = 0; num < row; num++)
803
for (unsigned int col = n-1; col+1 > 0; col--)
805
if (combinations[row][col] + 1 > 2)
806
combinations[row][col] = 0;
809
combinations[row][col] += 1;
816
// Compute inverse of Jacobian
817
const double Jinv[3][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}, {K[6], K[7], K[8]}};
819
// Declare transformation matrix
820
// Declare pointer to two dimensional array and initialise
821
double **transform = new double *[num_derivatives];
823
for (unsigned int j = 0; j < num_derivatives; j++)
825
transform[j] = new double [num_derivatives];
826
for (unsigned int k = 0; k < num_derivatives; k++)
830
// Construct transformation matrix
831
for (unsigned int row = 0; row < num_derivatives; row++)
833
for (unsigned int col = 0; col < num_derivatives; col++)
835
for (unsigned int k = 0; k < n; k++)
836
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
840
// Reset values. Assuming that values is always an array.
841
for (unsigned int r = 0; r < num_derivatives; r++)
844
}// end loop over 'r'
851
// Array of basisvalues
852
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
854
// Declare helper variables
855
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
857
// Compute basisvalues
858
basisvalues[0] = 1.0;
859
basisvalues[1] = tmp0;
860
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
861
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
862
basisvalues[0] *= std::sqrt(0.75);
863
basisvalues[3] *= std::sqrt(1.25);
864
basisvalues[2] *= std::sqrt(2.5);
865
basisvalues[1] *= std::sqrt(7.5);
867
// Table(s) of coefficients
868
static const double coefficients0[4] = \
869
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
871
// Tables of derivatives of the polynomial base (transpose).
872
static const double dmats0[4][4] = \
873
{{0.0, 0.0, 0.0, 0.0},
874
{6.32455532033676, 0.0, 0.0, 0.0},
875
{0.0, 0.0, 0.0, 0.0},
876
{0.0, 0.0, 0.0, 0.0}};
878
static const double dmats1[4][4] = \
879
{{0.0, 0.0, 0.0, 0.0},
880
{3.16227766016838, 0.0, 0.0, 0.0},
881
{5.47722557505166, 0.0, 0.0, 0.0},
882
{0.0, 0.0, 0.0, 0.0}};
884
static const double dmats2[4][4] = \
885
{{0.0, 0.0, 0.0, 0.0},
886
{3.16227766016838, 0.0, 0.0, 0.0},
887
{1.82574185835055, 0.0, 0.0, 0.0},
888
{5.16397779494322, 0.0, 0.0, 0.0}};
890
// Compute reference derivatives.
891
// Declare pointer to array of derivatives on FIAT element.
892
double *derivatives = new double[num_derivatives];
893
for (unsigned int r = 0; r < num_derivatives; r++)
895
derivatives[r] = 0.0;
896
}// end loop over 'r'
898
// Declare derivative matrix (of polynomial basis).
899
double dmats[4][4] = \
900
{{1.0, 0.0, 0.0, 0.0},
901
{0.0, 1.0, 0.0, 0.0},
902
{0.0, 0.0, 1.0, 0.0},
903
{0.0, 0.0, 0.0, 1.0}};
905
// Declare (auxiliary) derivative matrix (of polynomial basis).
906
double dmats_old[4][4] = \
907
{{1.0, 0.0, 0.0, 0.0},
908
{0.0, 1.0, 0.0, 0.0},
909
{0.0, 0.0, 1.0, 0.0},
910
{0.0, 0.0, 0.0, 1.0}};
912
// Loop possible derivatives.
913
for (unsigned int r = 0; r < num_derivatives; r++)
915
// Resetting dmats values to compute next derivative.
916
for (unsigned int t = 0; t < 4; t++)
918
for (unsigned int u = 0; u < 4; u++)
926
}// end loop over 'u'
927
}// end loop over 't'
929
// Looping derivative order to generate dmats.
930
for (unsigned int s = 0; s < n; s++)
932
// Updating dmats_old with new values and resetting dmats.
933
for (unsigned int t = 0; t < 4; t++)
935
for (unsigned int u = 0; u < 4; u++)
937
dmats_old[t][u] = dmats[t][u];
939
}// end loop over 'u'
940
}// end loop over 't'
942
// Update dmats using an inner product.
943
if (combinations[r][s] == 0)
945
for (unsigned int t = 0; t < 4; t++)
947
for (unsigned int u = 0; u < 4; u++)
949
for (unsigned int tu = 0; tu < 4; tu++)
951
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
952
}// end loop over 'tu'
953
}// end loop over 'u'
954
}// end loop over 't'
957
if (combinations[r][s] == 1)
959
for (unsigned int t = 0; t < 4; t++)
961
for (unsigned int u = 0; u < 4; u++)
963
for (unsigned int tu = 0; tu < 4; tu++)
965
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
966
}// end loop over 'tu'
967
}// end loop over 'u'
968
}// end loop over 't'
971
if (combinations[r][s] == 2)
973
for (unsigned int t = 0; t < 4; t++)
975
for (unsigned int u = 0; u < 4; u++)
977
for (unsigned int tu = 0; tu < 4; tu++)
979
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
980
}// end loop over 'tu'
981
}// end loop over 'u'
982
}// end loop over 't'
985
}// end loop over 's'
986
for (unsigned int s = 0; s < 4; s++)
988
for (unsigned int t = 0; t < 4; t++)
990
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
991
}// end loop over 't'
992
}// end loop over 's'
993
}// end loop over 'r'
995
// Transform derivatives back to physical element
996
for (unsigned int r = 0; r < num_derivatives; r++)
998
for (unsigned int s = 0; s < num_derivatives; s++)
1000
values[r] += transform[r][s]*derivatives[s];
1001
}// end loop over 's'
1002
}// end loop over 'r'
1004
// Delete pointer to array of derivatives on FIAT element
1005
delete [] derivatives;
1007
// Delete pointer to array of combinations of derivatives and transform
1008
for (unsigned int r = 0; r < num_derivatives; r++)
1010
delete [] combinations[r];
1011
}// end loop over 'r'
1012
delete [] combinations;
1013
for (unsigned int r = 0; r < num_derivatives; r++)
1015
delete [] transform[r];
1016
}// end loop over 'r'
1017
delete [] transform;
1023
// Array of basisvalues
1024
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1026
// Declare helper variables
1027
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1029
// Compute basisvalues
1030
basisvalues[0] = 1.0;
1031
basisvalues[1] = tmp0;
1032
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1033
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1034
basisvalues[0] *= std::sqrt(0.75);
1035
basisvalues[3] *= std::sqrt(1.25);
1036
basisvalues[2] *= std::sqrt(2.5);
1037
basisvalues[1] *= std::sqrt(7.5);
1039
// Table(s) of coefficients
1040
static const double coefficients0[4] = \
1041
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
1043
// Tables of derivatives of the polynomial base (transpose).
1044
static const double dmats0[4][4] = \
1045
{{0.0, 0.0, 0.0, 0.0},
1046
{6.32455532033676, 0.0, 0.0, 0.0},
1047
{0.0, 0.0, 0.0, 0.0},
1048
{0.0, 0.0, 0.0, 0.0}};
1050
static const double dmats1[4][4] = \
1051
{{0.0, 0.0, 0.0, 0.0},
1052
{3.16227766016838, 0.0, 0.0, 0.0},
1053
{5.47722557505166, 0.0, 0.0, 0.0},
1054
{0.0, 0.0, 0.0, 0.0}};
1056
static const double dmats2[4][4] = \
1057
{{0.0, 0.0, 0.0, 0.0},
1058
{3.16227766016838, 0.0, 0.0, 0.0},
1059
{1.82574185835055, 0.0, 0.0, 0.0},
1060
{5.16397779494322, 0.0, 0.0, 0.0}};
1062
// Compute reference derivatives.
1063
// Declare pointer to array of derivatives on FIAT element.
1064
double *derivatives = new double[num_derivatives];
1065
for (unsigned int r = 0; r < num_derivatives; r++)
1067
derivatives[r] = 0.0;
1068
}// end loop over 'r'
1070
// Declare derivative matrix (of polynomial basis).
1071
double dmats[4][4] = \
1072
{{1.0, 0.0, 0.0, 0.0},
1073
{0.0, 1.0, 0.0, 0.0},
1074
{0.0, 0.0, 1.0, 0.0},
1075
{0.0, 0.0, 0.0, 1.0}};
1077
// Declare (auxiliary) derivative matrix (of polynomial basis).
1078
double dmats_old[4][4] = \
1079
{{1.0, 0.0, 0.0, 0.0},
1080
{0.0, 1.0, 0.0, 0.0},
1081
{0.0, 0.0, 1.0, 0.0},
1082
{0.0, 0.0, 0.0, 1.0}};
1084
// Loop possible derivatives.
1085
for (unsigned int r = 0; r < num_derivatives; r++)
1087
// Resetting dmats values to compute next derivative.
1088
for (unsigned int t = 0; t < 4; t++)
1090
for (unsigned int u = 0; u < 4; u++)
1098
}// end loop over 'u'
1099
}// end loop over 't'
1101
// Looping derivative order to generate dmats.
1102
for (unsigned int s = 0; s < n; s++)
1104
// Updating dmats_old with new values and resetting dmats.
1105
for (unsigned int t = 0; t < 4; t++)
1107
for (unsigned int u = 0; u < 4; u++)
1109
dmats_old[t][u] = dmats[t][u];
1111
}// end loop over 'u'
1112
}// end loop over 't'
1114
// Update dmats using an inner product.
1115
if (combinations[r][s] == 0)
1117
for (unsigned int t = 0; t < 4; t++)
1119
for (unsigned int u = 0; u < 4; u++)
1121
for (unsigned int tu = 0; tu < 4; tu++)
1123
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1124
}// end loop over 'tu'
1125
}// end loop over 'u'
1126
}// end loop over 't'
1129
if (combinations[r][s] == 1)
1131
for (unsigned int t = 0; t < 4; t++)
1133
for (unsigned int u = 0; u < 4; u++)
1135
for (unsigned int tu = 0; tu < 4; tu++)
1137
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1138
}// end loop over 'tu'
1139
}// end loop over 'u'
1140
}// end loop over 't'
1143
if (combinations[r][s] == 2)
1145
for (unsigned int t = 0; t < 4; t++)
1147
for (unsigned int u = 0; u < 4; u++)
1149
for (unsigned int tu = 0; tu < 4; tu++)
1151
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1152
}// end loop over 'tu'
1153
}// end loop over 'u'
1154
}// end loop over 't'
1157
}// end loop over 's'
1158
for (unsigned int s = 0; s < 4; s++)
1160
for (unsigned int t = 0; t < 4; t++)
1162
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1163
}// end loop over 't'
1164
}// end loop over 's'
1165
}// end loop over 'r'
1167
// Transform derivatives back to physical element
1168
for (unsigned int r = 0; r < num_derivatives; r++)
1170
for (unsigned int s = 0; s < num_derivatives; s++)
1172
values[r] += transform[r][s]*derivatives[s];
1173
}// end loop over 's'
1174
}// end loop over 'r'
1176
// Delete pointer to array of derivatives on FIAT element
1177
delete [] derivatives;
1179
// Delete pointer to array of combinations of derivatives and transform
1180
for (unsigned int r = 0; r < num_derivatives; r++)
1182
delete [] combinations[r];
1183
}// end loop over 'r'
1184
delete [] combinations;
1185
for (unsigned int r = 0; r < num_derivatives; r++)
1187
delete [] transform[r];
1188
}// end loop over 'r'
1189
delete [] transform;
1195
// Array of basisvalues
1196
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1198
// Declare helper variables
1199
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1201
// Compute basisvalues
1202
basisvalues[0] = 1.0;
1203
basisvalues[1] = tmp0;
1204
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1205
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1206
basisvalues[0] *= std::sqrt(0.75);
1207
basisvalues[3] *= std::sqrt(1.25);
1208
basisvalues[2] *= std::sqrt(2.5);
1209
basisvalues[1] *= std::sqrt(7.5);
1211
// Table(s) of coefficients
1212
static const double coefficients0[4] = \
1213
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
1215
// Tables of derivatives of the polynomial base (transpose).
1216
static const double dmats0[4][4] = \
1217
{{0.0, 0.0, 0.0, 0.0},
1218
{6.32455532033676, 0.0, 0.0, 0.0},
1219
{0.0, 0.0, 0.0, 0.0},
1220
{0.0, 0.0, 0.0, 0.0}};
1222
static const double dmats1[4][4] = \
1223
{{0.0, 0.0, 0.0, 0.0},
1224
{3.16227766016838, 0.0, 0.0, 0.0},
1225
{5.47722557505166, 0.0, 0.0, 0.0},
1226
{0.0, 0.0, 0.0, 0.0}};
1228
static const double dmats2[4][4] = \
1229
{{0.0, 0.0, 0.0, 0.0},
1230
{3.16227766016838, 0.0, 0.0, 0.0},
1231
{1.82574185835055, 0.0, 0.0, 0.0},
1232
{5.16397779494322, 0.0, 0.0, 0.0}};
1234
// Compute reference derivatives.
1235
// Declare pointer to array of derivatives on FIAT element.
1236
double *derivatives = new double[num_derivatives];
1237
for (unsigned int r = 0; r < num_derivatives; r++)
1239
derivatives[r] = 0.0;
1240
}// end loop over 'r'
1242
// Declare derivative matrix (of polynomial basis).
1243
double dmats[4][4] = \
1244
{{1.0, 0.0, 0.0, 0.0},
1245
{0.0, 1.0, 0.0, 0.0},
1246
{0.0, 0.0, 1.0, 0.0},
1247
{0.0, 0.0, 0.0, 1.0}};
1249
// Declare (auxiliary) derivative matrix (of polynomial basis).
1250
double dmats_old[4][4] = \
1251
{{1.0, 0.0, 0.0, 0.0},
1252
{0.0, 1.0, 0.0, 0.0},
1253
{0.0, 0.0, 1.0, 0.0},
1254
{0.0, 0.0, 0.0, 1.0}};
1256
// Loop possible derivatives.
1257
for (unsigned int r = 0; r < num_derivatives; r++)
1259
// Resetting dmats values to compute next derivative.
1260
for (unsigned int t = 0; t < 4; t++)
1262
for (unsigned int u = 0; u < 4; u++)
1270
}// end loop over 'u'
1271
}// end loop over 't'
1273
// Looping derivative order to generate dmats.
1274
for (unsigned int s = 0; s < n; s++)
1276
// Updating dmats_old with new values and resetting dmats.
1277
for (unsigned int t = 0; t < 4; t++)
1279
for (unsigned int u = 0; u < 4; u++)
1281
dmats_old[t][u] = dmats[t][u];
1283
}// end loop over 'u'
1284
}// end loop over 't'
1286
// Update dmats using an inner product.
1287
if (combinations[r][s] == 0)
1289
for (unsigned int t = 0; t < 4; t++)
1291
for (unsigned int u = 0; u < 4; u++)
1293
for (unsigned int tu = 0; tu < 4; tu++)
1295
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1296
}// end loop over 'tu'
1297
}// end loop over 'u'
1298
}// end loop over 't'
1301
if (combinations[r][s] == 1)
1303
for (unsigned int t = 0; t < 4; t++)
1305
for (unsigned int u = 0; u < 4; u++)
1307
for (unsigned int tu = 0; tu < 4; tu++)
1309
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1310
}// end loop over 'tu'
1311
}// end loop over 'u'
1312
}// end loop over 't'
1315
if (combinations[r][s] == 2)
1317
for (unsigned int t = 0; t < 4; t++)
1319
for (unsigned int u = 0; u < 4; u++)
1321
for (unsigned int tu = 0; tu < 4; tu++)
1323
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1324
}// end loop over 'tu'
1325
}// end loop over 'u'
1326
}// end loop over 't'
1329
}// end loop over 's'
1330
for (unsigned int s = 0; s < 4; s++)
1332
for (unsigned int t = 0; t < 4; t++)
1334
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1335
}// end loop over 't'
1336
}// end loop over 's'
1337
}// end loop over 'r'
1339
// Transform derivatives back to physical element
1340
for (unsigned int r = 0; r < num_derivatives; r++)
1342
for (unsigned int s = 0; s < num_derivatives; s++)
1344
values[r] += transform[r][s]*derivatives[s];
1345
}// end loop over 's'
1346
}// end loop over 'r'
1348
// Delete pointer to array of derivatives on FIAT element
1349
delete [] derivatives;
1351
// Delete pointer to array of combinations of derivatives and transform
1352
for (unsigned int r = 0; r < num_derivatives; r++)
1354
delete [] combinations[r];
1355
}// end loop over 'r'
1356
delete [] combinations;
1357
for (unsigned int r = 0; r < num_derivatives; r++)
1359
delete [] transform[r];
1360
}// end loop over 'r'
1361
delete [] transform;
1367
// Array of basisvalues
1368
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1370
// Declare helper variables
1371
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1373
// Compute basisvalues
1374
basisvalues[0] = 1.0;
1375
basisvalues[1] = tmp0;
1376
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1377
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1378
basisvalues[0] *= std::sqrt(0.75);
1379
basisvalues[3] *= std::sqrt(1.25);
1380
basisvalues[2] *= std::sqrt(2.5);
1381
basisvalues[1] *= std::sqrt(7.5);
1383
// Table(s) of coefficients
1384
static const double coefficients0[4] = \
1385
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
1387
// Tables of derivatives of the polynomial base (transpose).
1388
static const double dmats0[4][4] = \
1389
{{0.0, 0.0, 0.0, 0.0},
1390
{6.32455532033676, 0.0, 0.0, 0.0},
1391
{0.0, 0.0, 0.0, 0.0},
1392
{0.0, 0.0, 0.0, 0.0}};
1394
static const double dmats1[4][4] = \
1395
{{0.0, 0.0, 0.0, 0.0},
1396
{3.16227766016838, 0.0, 0.0, 0.0},
1397
{5.47722557505166, 0.0, 0.0, 0.0},
1398
{0.0, 0.0, 0.0, 0.0}};
1400
static const double dmats2[4][4] = \
1401
{{0.0, 0.0, 0.0, 0.0},
1402
{3.16227766016838, 0.0, 0.0, 0.0},
1403
{1.82574185835055, 0.0, 0.0, 0.0},
1404
{5.16397779494322, 0.0, 0.0, 0.0}};
1406
// Compute reference derivatives.
1407
// Declare pointer to array of derivatives on FIAT element.
1408
double *derivatives = new double[num_derivatives];
1409
for (unsigned int r = 0; r < num_derivatives; r++)
1411
derivatives[r] = 0.0;
1412
}// end loop over 'r'
1414
// Declare derivative matrix (of polynomial basis).
1415
double dmats[4][4] = \
1416
{{1.0, 0.0, 0.0, 0.0},
1417
{0.0, 1.0, 0.0, 0.0},
1418
{0.0, 0.0, 1.0, 0.0},
1419
{0.0, 0.0, 0.0, 1.0}};
1421
// Declare (auxiliary) derivative matrix (of polynomial basis).
1422
double dmats_old[4][4] = \
1423
{{1.0, 0.0, 0.0, 0.0},
1424
{0.0, 1.0, 0.0, 0.0},
1425
{0.0, 0.0, 1.0, 0.0},
1426
{0.0, 0.0, 0.0, 1.0}};
1428
// Loop possible derivatives.
1429
for (unsigned int r = 0; r < num_derivatives; r++)
1431
// Resetting dmats values to compute next derivative.
1432
for (unsigned int t = 0; t < 4; t++)
1434
for (unsigned int u = 0; u < 4; u++)
1442
}// end loop over 'u'
1443
}// end loop over 't'
1445
// Looping derivative order to generate dmats.
1446
for (unsigned int s = 0; s < n; s++)
1448
// Updating dmats_old with new values and resetting dmats.
1449
for (unsigned int t = 0; t < 4; t++)
1451
for (unsigned int u = 0; u < 4; u++)
1453
dmats_old[t][u] = dmats[t][u];
1455
}// end loop over 'u'
1456
}// end loop over 't'
1458
// Update dmats using an inner product.
1459
if (combinations[r][s] == 0)
1461
for (unsigned int t = 0; t < 4; t++)
1463
for (unsigned int u = 0; u < 4; u++)
1465
for (unsigned int tu = 0; tu < 4; tu++)
1467
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1468
}// end loop over 'tu'
1469
}// end loop over 'u'
1470
}// end loop over 't'
1473
if (combinations[r][s] == 1)
1475
for (unsigned int t = 0; t < 4; t++)
1477
for (unsigned int u = 0; u < 4; u++)
1479
for (unsigned int tu = 0; tu < 4; tu++)
1481
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1482
}// end loop over 'tu'
1483
}// end loop over 'u'
1484
}// end loop over 't'
1487
if (combinations[r][s] == 2)
1489
for (unsigned int t = 0; t < 4; t++)
1491
for (unsigned int u = 0; u < 4; u++)
1493
for (unsigned int tu = 0; tu < 4; tu++)
1495
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1496
}// end loop over 'tu'
1497
}// end loop over 'u'
1498
}// end loop over 't'
1501
}// end loop over 's'
1502
for (unsigned int s = 0; s < 4; s++)
1504
for (unsigned int t = 0; t < 4; t++)
1506
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1507
}// end loop over 't'
1508
}// end loop over 's'
1509
}// end loop over 'r'
1511
// Transform derivatives back to physical element
1512
for (unsigned int r = 0; r < num_derivatives; r++)
1514
for (unsigned int s = 0; s < num_derivatives; s++)
1516
values[r] += transform[r][s]*derivatives[s];
1517
}// end loop over 's'
1518
}// end loop over 'r'
1520
// Delete pointer to array of derivatives on FIAT element
1521
delete [] derivatives;
1523
// Delete pointer to array of combinations of derivatives and transform
1524
for (unsigned int r = 0; r < num_derivatives; r++)
1526
delete [] combinations[r];
1527
}// end loop over 'r'
1528
delete [] combinations;
1529
for (unsigned int r = 0; r < num_derivatives; r++)
1531
delete [] transform[r];
1532
}// end loop over 'r'
1533
delete [] transform;
1540
/// Evaluate order n derivatives of all basis functions at given point x in cell
1541
virtual void evaluate_basis_derivatives_all(std::size_t n,
1544
const double* vertex_coordinates,
1545
int cell_orientation) const
1547
// Compute number of derivatives.
1548
unsigned int num_derivatives = 1;
1549
for (unsigned int r = 0; r < n; r++)
1551
num_derivatives *= 3;
1552
}// end loop over 'r'
1554
// Helper variable to hold values of a single dof.
1555
double *dof_values = new double[num_derivatives];
1556
for (unsigned int r = 0; r < num_derivatives; r++)
1558
dof_values[r] = 0.0;
1559
}// end loop over 'r'
1561
// Loop dofs and call evaluate_basis_derivatives.
1562
for (unsigned int r = 0; r < 4; r++)
1564
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1565
for (unsigned int s = 0; s < num_derivatives; s++)
1567
values[r*num_derivatives + s] = dof_values[s];
1568
}// end loop over 's'
1569
}// end loop over 'r'
1572
delete [] dof_values;
1575
/// Evaluate linear functional for dof i on the function f
1576
virtual double evaluate_dof(std::size_t i,
1577
const ufc::function& f,
1578
const double* vertex_coordinates,
1579
int cell_orientation,
1580
const ufc::cell& c) const
1582
// Declare variables for result of evaluation
1585
// Declare variable for physical coordinates
1591
y[0] = vertex_coordinates[0];
1592
y[1] = vertex_coordinates[1];
1593
y[2] = vertex_coordinates[2];
1594
f.evaluate(vals, y, c);
1600
y[0] = vertex_coordinates[3];
1601
y[1] = vertex_coordinates[4];
1602
y[2] = vertex_coordinates[5];
1603
f.evaluate(vals, y, c);
1609
y[0] = vertex_coordinates[6];
1610
y[1] = vertex_coordinates[7];
1611
y[2] = vertex_coordinates[8];
1612
f.evaluate(vals, y, c);
1618
y[0] = vertex_coordinates[9];
1619
y[1] = vertex_coordinates[10];
1620
y[2] = vertex_coordinates[11];
1621
f.evaluate(vals, y, c);
1630
/// Evaluate linear functionals for all dofs on the function f
1631
virtual void evaluate_dofs(double* values,
1632
const ufc::function& f,
1633
const double* vertex_coordinates,
1634
int cell_orientation,
1635
const ufc::cell& c) const
1637
// Declare variables for result of evaluation
1640
// Declare variable for physical coordinates
1642
y[0] = vertex_coordinates[0];
1643
y[1] = vertex_coordinates[1];
1644
y[2] = vertex_coordinates[2];
1645
f.evaluate(vals, y, c);
1646
values[0] = vals[0];
1647
y[0] = vertex_coordinates[3];
1648
y[1] = vertex_coordinates[4];
1649
y[2] = vertex_coordinates[5];
1650
f.evaluate(vals, y, c);
1651
values[1] = vals[0];
1652
y[0] = vertex_coordinates[6];
1653
y[1] = vertex_coordinates[7];
1654
y[2] = vertex_coordinates[8];
1655
f.evaluate(vals, y, c);
1656
values[2] = vals[0];
1657
y[0] = vertex_coordinates[9];
1658
y[1] = vertex_coordinates[10];
1659
y[2] = vertex_coordinates[11];
1660
f.evaluate(vals, y, c);
1661
values[3] = vals[0];
1664
/// Interpolate vertex values from dof values
1665
virtual void interpolate_vertex_values(double* vertex_values,
1666
const double* dof_values,
1667
const double* vertex_coordinates,
1668
int cell_orientation,
1669
const ufc::cell& c) const
1671
// Evaluate function and change variables
1672
vertex_values[0] = dof_values[0];
1673
vertex_values[1] = dof_values[1];
1674
vertex_values[2] = dof_values[2];
1675
vertex_values[3] = dof_values[3];
1678
/// Map coordinate xhat from reference cell to coordinate x in cell
1679
virtual void map_from_reference_cell(double* x,
1681
const ufc::cell& c) const
1683
throw std::runtime_error("map_from_reference_cell not yet implemented.");
1686
/// Map from coordinate x in cell to coordinate xhat in reference cell
1687
virtual void map_to_reference_cell(double* xhat,
1689
const ufc::cell& c) const
1691
throw std::runtime_error("map_to_reference_cell not yet implemented.");
1694
/// Return the number of sub elements (for a mixed element)
1695
virtual std::size_t num_sub_elements() const
1700
/// Create a new finite element for sub element i (for a mixed element)
1701
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1706
/// Create a new class instance
1707
virtual ufc::finite_element* create() const
1709
return new navierstokes_finite_element_1();
1714
/// This class defines the interface for a finite element.
1716
class navierstokes_finite_element_2: public ufc::finite_element
1721
navierstokes_finite_element_2() : ufc::finite_element()
1727
virtual ~navierstokes_finite_element_2()
1732
/// Return a string identifying the finite element
1733
virtual const char* signature() const
1735
return "VectorElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, 3, None)";
1738
/// Return the cell shape
1739
virtual ufc::shape cell_shape() const
1741
return ufc::tetrahedron;
1744
/// Return the topological dimension of the cell shape
1745
virtual std::size_t topological_dimension() const
1750
/// Return the geometric dimension of the cell shape
1751
virtual std::size_t geometric_dimension() const
1756
/// Return the dimension of the finite element function space
1757
virtual std::size_t space_dimension() const
1762
/// Return the rank of the value space
1763
virtual std::size_t value_rank() const
1768
/// Return the dimension of the value space for axis i
1769
virtual std::size_t value_dimension(std::size_t i) const
1783
/// Evaluate basis function i at given point x in cell
1784
virtual void evaluate_basis(std::size_t i,
1787
const double* vertex_coordinates,
1788
int cell_orientation) const
1792
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
1794
// Compute Jacobian inverse and determinant
1797
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
1800
// Compute constants
1801
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
1802
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
1803
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
1805
// Compute subdeterminants
1806
const double d_00 = J[4]*J[8] - J[5]*J[7];
1807
const double d_01 = J[5]*J[6] - J[3]*J[8];
1808
const double d_02 = J[3]*J[7] - J[4]*J[6];
1809
const double d_10 = J[2]*J[7] - J[1]*J[8];
1810
const double d_11 = J[0]*J[8] - J[2]*J[6];
1811
const double d_12 = J[1]*J[6] - J[0]*J[7];
1812
const double d_20 = J[1]*J[5] - J[2]*J[4];
1813
const double d_21 = J[2]*J[3] - J[0]*J[5];
1814
const double d_22 = J[0]*J[4] - J[1]*J[3];
1816
// Get coordinates and map to the reference (FIAT) element
1817
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
1818
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
1819
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
1831
// Array of basisvalues
1832
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1834
// Declare helper variables
1835
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1837
// Compute basisvalues
1838
basisvalues[0] = 1.0;
1839
basisvalues[1] = tmp0;
1840
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1841
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1842
basisvalues[0] *= std::sqrt(0.75);
1843
basisvalues[3] *= std::sqrt(1.25);
1844
basisvalues[2] *= std::sqrt(2.5);
1845
basisvalues[1] *= std::sqrt(7.5);
1847
// Table(s) of coefficients
1848
static const double coefficients0[4] = \
1849
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
1852
for (unsigned int r = 0; r < 4; r++)
1854
values[0] += coefficients0[r]*basisvalues[r];
1855
}// end loop over 'r'
1861
// Array of basisvalues
1862
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1864
// Declare helper variables
1865
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1867
// Compute basisvalues
1868
basisvalues[0] = 1.0;
1869
basisvalues[1] = tmp0;
1870
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1871
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1872
basisvalues[0] *= std::sqrt(0.75);
1873
basisvalues[3] *= std::sqrt(1.25);
1874
basisvalues[2] *= std::sqrt(2.5);
1875
basisvalues[1] *= std::sqrt(7.5);
1877
// Table(s) of coefficients
1878
static const double coefficients0[4] = \
1879
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
1882
for (unsigned int r = 0; r < 4; r++)
1884
values[0] += coefficients0[r]*basisvalues[r];
1885
}// end loop over 'r'
1891
// Array of basisvalues
1892
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1894
// Declare helper variables
1895
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1897
// Compute basisvalues
1898
basisvalues[0] = 1.0;
1899
basisvalues[1] = tmp0;
1900
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1901
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1902
basisvalues[0] *= std::sqrt(0.75);
1903
basisvalues[3] *= std::sqrt(1.25);
1904
basisvalues[2] *= std::sqrt(2.5);
1905
basisvalues[1] *= std::sqrt(7.5);
1907
// Table(s) of coefficients
1908
static const double coefficients0[4] = \
1909
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
1912
for (unsigned int r = 0; r < 4; r++)
1914
values[0] += coefficients0[r]*basisvalues[r];
1915
}// end loop over 'r'
1921
// Array of basisvalues
1922
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1924
// Declare helper variables
1925
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1927
// Compute basisvalues
1928
basisvalues[0] = 1.0;
1929
basisvalues[1] = tmp0;
1930
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1931
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1932
basisvalues[0] *= std::sqrt(0.75);
1933
basisvalues[3] *= std::sqrt(1.25);
1934
basisvalues[2] *= std::sqrt(2.5);
1935
basisvalues[1] *= std::sqrt(7.5);
1937
// Table(s) of coefficients
1938
static const double coefficients0[4] = \
1939
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
1942
for (unsigned int r = 0; r < 4; r++)
1944
values[0] += coefficients0[r]*basisvalues[r];
1945
}// end loop over 'r'
1951
// Array of basisvalues
1952
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1954
// Declare helper variables
1955
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1957
// Compute basisvalues
1958
basisvalues[0] = 1.0;
1959
basisvalues[1] = tmp0;
1960
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1961
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1962
basisvalues[0] *= std::sqrt(0.75);
1963
basisvalues[3] *= std::sqrt(1.25);
1964
basisvalues[2] *= std::sqrt(2.5);
1965
basisvalues[1] *= std::sqrt(7.5);
1967
// Table(s) of coefficients
1968
static const double coefficients0[4] = \
1969
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
1972
for (unsigned int r = 0; r < 4; r++)
1974
values[1] += coefficients0[r]*basisvalues[r];
1975
}// end loop over 'r'
1981
// Array of basisvalues
1982
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
1984
// Declare helper variables
1985
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1987
// Compute basisvalues
1988
basisvalues[0] = 1.0;
1989
basisvalues[1] = tmp0;
1990
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1991
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1992
basisvalues[0] *= std::sqrt(0.75);
1993
basisvalues[3] *= std::sqrt(1.25);
1994
basisvalues[2] *= std::sqrt(2.5);
1995
basisvalues[1] *= std::sqrt(7.5);
1997
// Table(s) of coefficients
1998
static const double coefficients0[4] = \
1999
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
2002
for (unsigned int r = 0; r < 4; r++)
2004
values[1] += coefficients0[r]*basisvalues[r];
2005
}// end loop over 'r'
2011
// Array of basisvalues
2012
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2014
// Declare helper variables
2015
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2017
// Compute basisvalues
2018
basisvalues[0] = 1.0;
2019
basisvalues[1] = tmp0;
2020
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2021
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2022
basisvalues[0] *= std::sqrt(0.75);
2023
basisvalues[3] *= std::sqrt(1.25);
2024
basisvalues[2] *= std::sqrt(2.5);
2025
basisvalues[1] *= std::sqrt(7.5);
2027
// Table(s) of coefficients
2028
static const double coefficients0[4] = \
2029
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
2032
for (unsigned int r = 0; r < 4; r++)
2034
values[1] += coefficients0[r]*basisvalues[r];
2035
}// end loop over 'r'
2041
// Array of basisvalues
2042
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2044
// Declare helper variables
2045
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2047
// Compute basisvalues
2048
basisvalues[0] = 1.0;
2049
basisvalues[1] = tmp0;
2050
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2051
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2052
basisvalues[0] *= std::sqrt(0.75);
2053
basisvalues[3] *= std::sqrt(1.25);
2054
basisvalues[2] *= std::sqrt(2.5);
2055
basisvalues[1] *= std::sqrt(7.5);
2057
// Table(s) of coefficients
2058
static const double coefficients0[4] = \
2059
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
2062
for (unsigned int r = 0; r < 4; r++)
2064
values[1] += coefficients0[r]*basisvalues[r];
2065
}// end loop over 'r'
2071
// Array of basisvalues
2072
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2074
// Declare helper variables
2075
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2077
// Compute basisvalues
2078
basisvalues[0] = 1.0;
2079
basisvalues[1] = tmp0;
2080
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2081
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2082
basisvalues[0] *= std::sqrt(0.75);
2083
basisvalues[3] *= std::sqrt(1.25);
2084
basisvalues[2] *= std::sqrt(2.5);
2085
basisvalues[1] *= std::sqrt(7.5);
2087
// Table(s) of coefficients
2088
static const double coefficients0[4] = \
2089
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
2092
for (unsigned int r = 0; r < 4; r++)
2094
values[2] += coefficients0[r]*basisvalues[r];
2095
}// end loop over 'r'
2101
// Array of basisvalues
2102
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2104
// Declare helper variables
2105
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2107
// Compute basisvalues
2108
basisvalues[0] = 1.0;
2109
basisvalues[1] = tmp0;
2110
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2111
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2112
basisvalues[0] *= std::sqrt(0.75);
2113
basisvalues[3] *= std::sqrt(1.25);
2114
basisvalues[2] *= std::sqrt(2.5);
2115
basisvalues[1] *= std::sqrt(7.5);
2117
// Table(s) of coefficients
2118
static const double coefficients0[4] = \
2119
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
2122
for (unsigned int r = 0; r < 4; r++)
2124
values[2] += coefficients0[r]*basisvalues[r];
2125
}// end loop over 'r'
2131
// Array of basisvalues
2132
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2134
// Declare helper variables
2135
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2137
// Compute basisvalues
2138
basisvalues[0] = 1.0;
2139
basisvalues[1] = tmp0;
2140
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2141
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2142
basisvalues[0] *= std::sqrt(0.75);
2143
basisvalues[3] *= std::sqrt(1.25);
2144
basisvalues[2] *= std::sqrt(2.5);
2145
basisvalues[1] *= std::sqrt(7.5);
2147
// Table(s) of coefficients
2148
static const double coefficients0[4] = \
2149
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
2152
for (unsigned int r = 0; r < 4; r++)
2154
values[2] += coefficients0[r]*basisvalues[r];
2155
}// end loop over 'r'
2161
// Array of basisvalues
2162
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2164
// Declare helper variables
2165
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2167
// Compute basisvalues
2168
basisvalues[0] = 1.0;
2169
basisvalues[1] = tmp0;
2170
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2171
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2172
basisvalues[0] *= std::sqrt(0.75);
2173
basisvalues[3] *= std::sqrt(1.25);
2174
basisvalues[2] *= std::sqrt(2.5);
2175
basisvalues[1] *= std::sqrt(7.5);
2177
// Table(s) of coefficients
2178
static const double coefficients0[4] = \
2179
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
2182
for (unsigned int r = 0; r < 4; r++)
2184
values[2] += coefficients0[r]*basisvalues[r];
2185
}// end loop over 'r'
2192
/// Evaluate all basis functions at given point x in cell
2193
virtual void evaluate_basis_all(double* values,
2195
const double* vertex_coordinates,
2196
int cell_orientation) const
2198
// Helper variable to hold values of a single dof.
2199
double dof_values[3] = {0.0, 0.0, 0.0};
2201
// Loop dofs and call evaluate_basis
2202
for (unsigned int r = 0; r < 12; r++)
2204
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
2205
for (unsigned int s = 0; s < 3; s++)
2207
values[r*3 + s] = dof_values[s];
2208
}// end loop over 's'
2209
}// end loop over 'r'
2212
/// Evaluate order n derivatives of basis function i at given point x in cell
2213
virtual void evaluate_basis_derivatives(std::size_t i,
2217
const double* vertex_coordinates,
2218
int cell_orientation) const
2222
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
2224
// Compute Jacobian inverse and determinant
2227
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
2230
// Compute constants
2231
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
2232
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
2233
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
2235
// Compute subdeterminants
2236
const double d_00 = J[4]*J[8] - J[5]*J[7];
2237
const double d_01 = J[5]*J[6] - J[3]*J[8];
2238
const double d_02 = J[3]*J[7] - J[4]*J[6];
2239
const double d_10 = J[2]*J[7] - J[1]*J[8];
2240
const double d_11 = J[0]*J[8] - J[2]*J[6];
2241
const double d_12 = J[1]*J[6] - J[0]*J[7];
2242
const double d_20 = J[1]*J[5] - J[2]*J[4];
2243
const double d_21 = J[2]*J[3] - J[0]*J[5];
2244
const double d_22 = J[0]*J[4] - J[1]*J[3];
2246
// Get coordinates and map to the reference (FIAT) element
2247
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
2248
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
2249
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
2252
// Compute number of derivatives.
2253
unsigned int num_derivatives = 1;
2254
for (unsigned int r = 0; r < n; r++)
2256
num_derivatives *= 3;
2257
}// end loop over 'r'
2259
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
2260
unsigned int **combinations = new unsigned int *[num_derivatives];
2261
for (unsigned int row = 0; row < num_derivatives; row++)
2263
combinations[row] = new unsigned int [n];
2264
for (unsigned int col = 0; col < n; col++)
2265
combinations[row][col] = 0;
2268
// Generate combinations of derivatives
2269
for (unsigned int row = 1; row < num_derivatives; row++)
2271
for (unsigned int num = 0; num < row; num++)
2273
for (unsigned int col = n-1; col+1 > 0; col--)
2275
if (combinations[row][col] + 1 > 2)
2276
combinations[row][col] = 0;
2279
combinations[row][col] += 1;
2286
// Compute inverse of Jacobian
2287
const double Jinv[3][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}, {K[6], K[7], K[8]}};
2289
// Declare transformation matrix
2290
// Declare pointer to two dimensional array and initialise
2291
double **transform = new double *[num_derivatives];
2293
for (unsigned int j = 0; j < num_derivatives; j++)
2295
transform[j] = new double [num_derivatives];
2296
for (unsigned int k = 0; k < num_derivatives; k++)
2297
transform[j][k] = 1;
2300
// Construct transformation matrix
2301
for (unsigned int row = 0; row < num_derivatives; row++)
2303
for (unsigned int col = 0; col < num_derivatives; col++)
2305
for (unsigned int k = 0; k < n; k++)
2306
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2310
// Reset values. Assuming that values is always an array.
2311
for (unsigned int r = 0; r < 3*num_derivatives; r++)
2314
}// end loop over 'r'
2321
// Array of basisvalues
2322
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2324
// Declare helper variables
2325
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2327
// Compute basisvalues
2328
basisvalues[0] = 1.0;
2329
basisvalues[1] = tmp0;
2330
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2331
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2332
basisvalues[0] *= std::sqrt(0.75);
2333
basisvalues[3] *= std::sqrt(1.25);
2334
basisvalues[2] *= std::sqrt(2.5);
2335
basisvalues[1] *= std::sqrt(7.5);
2337
// Table(s) of coefficients
2338
static const double coefficients0[4] = \
2339
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
2341
// Tables of derivatives of the polynomial base (transpose).
2342
static const double dmats0[4][4] = \
2343
{{0.0, 0.0, 0.0, 0.0},
2344
{6.32455532033676, 0.0, 0.0, 0.0},
2345
{0.0, 0.0, 0.0, 0.0},
2346
{0.0, 0.0, 0.0, 0.0}};
2348
static const double dmats1[4][4] = \
2349
{{0.0, 0.0, 0.0, 0.0},
2350
{3.16227766016838, 0.0, 0.0, 0.0},
2351
{5.47722557505166, 0.0, 0.0, 0.0},
2352
{0.0, 0.0, 0.0, 0.0}};
2354
static const double dmats2[4][4] = \
2355
{{0.0, 0.0, 0.0, 0.0},
2356
{3.16227766016838, 0.0, 0.0, 0.0},
2357
{1.82574185835055, 0.0, 0.0, 0.0},
2358
{5.16397779494322, 0.0, 0.0, 0.0}};
2360
// Compute reference derivatives.
2361
// Declare pointer to array of derivatives on FIAT element.
2362
double *derivatives = new double[num_derivatives];
2363
for (unsigned int r = 0; r < num_derivatives; r++)
2365
derivatives[r] = 0.0;
2366
}// end loop over 'r'
2368
// Declare derivative matrix (of polynomial basis).
2369
double dmats[4][4] = \
2370
{{1.0, 0.0, 0.0, 0.0},
2371
{0.0, 1.0, 0.0, 0.0},
2372
{0.0, 0.0, 1.0, 0.0},
2373
{0.0, 0.0, 0.0, 1.0}};
2375
// Declare (auxiliary) derivative matrix (of polynomial basis).
2376
double dmats_old[4][4] = \
2377
{{1.0, 0.0, 0.0, 0.0},
2378
{0.0, 1.0, 0.0, 0.0},
2379
{0.0, 0.0, 1.0, 0.0},
2380
{0.0, 0.0, 0.0, 1.0}};
2382
// Loop possible derivatives.
2383
for (unsigned int r = 0; r < num_derivatives; r++)
2385
// Resetting dmats values to compute next derivative.
2386
for (unsigned int t = 0; t < 4; t++)
2388
for (unsigned int u = 0; u < 4; u++)
2396
}// end loop over 'u'
2397
}// end loop over 't'
2399
// Looping derivative order to generate dmats.
2400
for (unsigned int s = 0; s < n; s++)
2402
// Updating dmats_old with new values and resetting dmats.
2403
for (unsigned int t = 0; t < 4; t++)
2405
for (unsigned int u = 0; u < 4; u++)
2407
dmats_old[t][u] = dmats[t][u];
2409
}// end loop over 'u'
2410
}// end loop over 't'
2412
// Update dmats using an inner product.
2413
if (combinations[r][s] == 0)
2415
for (unsigned int t = 0; t < 4; t++)
2417
for (unsigned int u = 0; u < 4; u++)
2419
for (unsigned int tu = 0; tu < 4; tu++)
2421
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2422
}// end loop over 'tu'
2423
}// end loop over 'u'
2424
}// end loop over 't'
2427
if (combinations[r][s] == 1)
2429
for (unsigned int t = 0; t < 4; t++)
2431
for (unsigned int u = 0; u < 4; u++)
2433
for (unsigned int tu = 0; tu < 4; tu++)
2435
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2436
}// end loop over 'tu'
2437
}// end loop over 'u'
2438
}// end loop over 't'
2441
if (combinations[r][s] == 2)
2443
for (unsigned int t = 0; t < 4; t++)
2445
for (unsigned int u = 0; u < 4; u++)
2447
for (unsigned int tu = 0; tu < 4; tu++)
2449
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2450
}// end loop over 'tu'
2451
}// end loop over 'u'
2452
}// end loop over 't'
2455
}// end loop over 's'
2456
for (unsigned int s = 0; s < 4; s++)
2458
for (unsigned int t = 0; t < 4; t++)
2460
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2461
}// end loop over 't'
2462
}// end loop over 's'
2463
}// end loop over 'r'
2465
// Transform derivatives back to physical element
2466
for (unsigned int r = 0; r < num_derivatives; r++)
2468
for (unsigned int s = 0; s < num_derivatives; s++)
2470
values[r] += transform[r][s]*derivatives[s];
2471
}// end loop over 's'
2472
}// end loop over 'r'
2474
// Delete pointer to array of derivatives on FIAT element
2475
delete [] derivatives;
2477
// Delete pointer to array of combinations of derivatives and transform
2478
for (unsigned int r = 0; r < num_derivatives; r++)
2480
delete [] combinations[r];
2481
}// end loop over 'r'
2482
delete [] combinations;
2483
for (unsigned int r = 0; r < num_derivatives; r++)
2485
delete [] transform[r];
2486
}// end loop over 'r'
2487
delete [] transform;
2493
// Array of basisvalues
2494
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2496
// Declare helper variables
2497
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2499
// Compute basisvalues
2500
basisvalues[0] = 1.0;
2501
basisvalues[1] = tmp0;
2502
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2503
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2504
basisvalues[0] *= std::sqrt(0.75);
2505
basisvalues[3] *= std::sqrt(1.25);
2506
basisvalues[2] *= std::sqrt(2.5);
2507
basisvalues[1] *= std::sqrt(7.5);
2509
// Table(s) of coefficients
2510
static const double coefficients0[4] = \
2511
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
2513
// Tables of derivatives of the polynomial base (transpose).
2514
static const double dmats0[4][4] = \
2515
{{0.0, 0.0, 0.0, 0.0},
2516
{6.32455532033676, 0.0, 0.0, 0.0},
2517
{0.0, 0.0, 0.0, 0.0},
2518
{0.0, 0.0, 0.0, 0.0}};
2520
static const double dmats1[4][4] = \
2521
{{0.0, 0.0, 0.0, 0.0},
2522
{3.16227766016838, 0.0, 0.0, 0.0},
2523
{5.47722557505166, 0.0, 0.0, 0.0},
2524
{0.0, 0.0, 0.0, 0.0}};
2526
static const double dmats2[4][4] = \
2527
{{0.0, 0.0, 0.0, 0.0},
2528
{3.16227766016838, 0.0, 0.0, 0.0},
2529
{1.82574185835055, 0.0, 0.0, 0.0},
2530
{5.16397779494322, 0.0, 0.0, 0.0}};
2532
// Compute reference derivatives.
2533
// Declare pointer to array of derivatives on FIAT element.
2534
double *derivatives = new double[num_derivatives];
2535
for (unsigned int r = 0; r < num_derivatives; r++)
2537
derivatives[r] = 0.0;
2538
}// end loop over 'r'
2540
// Declare derivative matrix (of polynomial basis).
2541
double dmats[4][4] = \
2542
{{1.0, 0.0, 0.0, 0.0},
2543
{0.0, 1.0, 0.0, 0.0},
2544
{0.0, 0.0, 1.0, 0.0},
2545
{0.0, 0.0, 0.0, 1.0}};
2547
// Declare (auxiliary) derivative matrix (of polynomial basis).
2548
double dmats_old[4][4] = \
2549
{{1.0, 0.0, 0.0, 0.0},
2550
{0.0, 1.0, 0.0, 0.0},
2551
{0.0, 0.0, 1.0, 0.0},
2552
{0.0, 0.0, 0.0, 1.0}};
2554
// Loop possible derivatives.
2555
for (unsigned int r = 0; r < num_derivatives; r++)
2557
// Resetting dmats values to compute next derivative.
2558
for (unsigned int t = 0; t < 4; t++)
2560
for (unsigned int u = 0; u < 4; u++)
2568
}// end loop over 'u'
2569
}// end loop over 't'
2571
// Looping derivative order to generate dmats.
2572
for (unsigned int s = 0; s < n; s++)
2574
// Updating dmats_old with new values and resetting dmats.
2575
for (unsigned int t = 0; t < 4; t++)
2577
for (unsigned int u = 0; u < 4; u++)
2579
dmats_old[t][u] = dmats[t][u];
2581
}// end loop over 'u'
2582
}// end loop over 't'
2584
// Update dmats using an inner product.
2585
if (combinations[r][s] == 0)
2587
for (unsigned int t = 0; t < 4; t++)
2589
for (unsigned int u = 0; u < 4; u++)
2591
for (unsigned int tu = 0; tu < 4; tu++)
2593
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2594
}// end loop over 'tu'
2595
}// end loop over 'u'
2596
}// end loop over 't'
2599
if (combinations[r][s] == 1)
2601
for (unsigned int t = 0; t < 4; t++)
2603
for (unsigned int u = 0; u < 4; u++)
2605
for (unsigned int tu = 0; tu < 4; tu++)
2607
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2608
}// end loop over 'tu'
2609
}// end loop over 'u'
2610
}// end loop over 't'
2613
if (combinations[r][s] == 2)
2615
for (unsigned int t = 0; t < 4; t++)
2617
for (unsigned int u = 0; u < 4; u++)
2619
for (unsigned int tu = 0; tu < 4; tu++)
2621
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2622
}// end loop over 'tu'
2623
}// end loop over 'u'
2624
}// end loop over 't'
2627
}// end loop over 's'
2628
for (unsigned int s = 0; s < 4; s++)
2630
for (unsigned int t = 0; t < 4; t++)
2632
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2633
}// end loop over 't'
2634
}// end loop over 's'
2635
}// end loop over 'r'
2637
// Transform derivatives back to physical element
2638
for (unsigned int r = 0; r < num_derivatives; r++)
2640
for (unsigned int s = 0; s < num_derivatives; s++)
2642
values[r] += transform[r][s]*derivatives[s];
2643
}// end loop over 's'
2644
}// end loop over 'r'
2646
// Delete pointer to array of derivatives on FIAT element
2647
delete [] derivatives;
2649
// Delete pointer to array of combinations of derivatives and transform
2650
for (unsigned int r = 0; r < num_derivatives; r++)
2652
delete [] combinations[r];
2653
}// end loop over 'r'
2654
delete [] combinations;
2655
for (unsigned int r = 0; r < num_derivatives; r++)
2657
delete [] transform[r];
2658
}// end loop over 'r'
2659
delete [] transform;
2665
// Array of basisvalues
2666
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2668
// Declare helper variables
2669
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2671
// Compute basisvalues
2672
basisvalues[0] = 1.0;
2673
basisvalues[1] = tmp0;
2674
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2675
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2676
basisvalues[0] *= std::sqrt(0.75);
2677
basisvalues[3] *= std::sqrt(1.25);
2678
basisvalues[2] *= std::sqrt(2.5);
2679
basisvalues[1] *= std::sqrt(7.5);
2681
// Table(s) of coefficients
2682
static const double coefficients0[4] = \
2683
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
2685
// Tables of derivatives of the polynomial base (transpose).
2686
static const double dmats0[4][4] = \
2687
{{0.0, 0.0, 0.0, 0.0},
2688
{6.32455532033676, 0.0, 0.0, 0.0},
2689
{0.0, 0.0, 0.0, 0.0},
2690
{0.0, 0.0, 0.0, 0.0}};
2692
static const double dmats1[4][4] = \
2693
{{0.0, 0.0, 0.0, 0.0},
2694
{3.16227766016838, 0.0, 0.0, 0.0},
2695
{5.47722557505166, 0.0, 0.0, 0.0},
2696
{0.0, 0.0, 0.0, 0.0}};
2698
static const double dmats2[4][4] = \
2699
{{0.0, 0.0, 0.0, 0.0},
2700
{3.16227766016838, 0.0, 0.0, 0.0},
2701
{1.82574185835055, 0.0, 0.0, 0.0},
2702
{5.16397779494322, 0.0, 0.0, 0.0}};
2704
// Compute reference derivatives.
2705
// Declare pointer to array of derivatives on FIAT element.
2706
double *derivatives = new double[num_derivatives];
2707
for (unsigned int r = 0; r < num_derivatives; r++)
2709
derivatives[r] = 0.0;
2710
}// end loop over 'r'
2712
// Declare derivative matrix (of polynomial basis).
2713
double dmats[4][4] = \
2714
{{1.0, 0.0, 0.0, 0.0},
2715
{0.0, 1.0, 0.0, 0.0},
2716
{0.0, 0.0, 1.0, 0.0},
2717
{0.0, 0.0, 0.0, 1.0}};
2719
// Declare (auxiliary) derivative matrix (of polynomial basis).
2720
double dmats_old[4][4] = \
2721
{{1.0, 0.0, 0.0, 0.0},
2722
{0.0, 1.0, 0.0, 0.0},
2723
{0.0, 0.0, 1.0, 0.0},
2724
{0.0, 0.0, 0.0, 1.0}};
2726
// Loop possible derivatives.
2727
for (unsigned int r = 0; r < num_derivatives; r++)
2729
// Resetting dmats values to compute next derivative.
2730
for (unsigned int t = 0; t < 4; t++)
2732
for (unsigned int u = 0; u < 4; u++)
2740
}// end loop over 'u'
2741
}// end loop over 't'
2743
// Looping derivative order to generate dmats.
2744
for (unsigned int s = 0; s < n; s++)
2746
// Updating dmats_old with new values and resetting dmats.
2747
for (unsigned int t = 0; t < 4; t++)
2749
for (unsigned int u = 0; u < 4; u++)
2751
dmats_old[t][u] = dmats[t][u];
2753
}// end loop over 'u'
2754
}// end loop over 't'
2756
// Update dmats using an inner product.
2757
if (combinations[r][s] == 0)
2759
for (unsigned int t = 0; t < 4; t++)
2761
for (unsigned int u = 0; u < 4; u++)
2763
for (unsigned int tu = 0; tu < 4; tu++)
2765
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2766
}// end loop over 'tu'
2767
}// end loop over 'u'
2768
}// end loop over 't'
2771
if (combinations[r][s] == 1)
2773
for (unsigned int t = 0; t < 4; t++)
2775
for (unsigned int u = 0; u < 4; u++)
2777
for (unsigned int tu = 0; tu < 4; tu++)
2779
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2780
}// end loop over 'tu'
2781
}// end loop over 'u'
2782
}// end loop over 't'
2785
if (combinations[r][s] == 2)
2787
for (unsigned int t = 0; t < 4; t++)
2789
for (unsigned int u = 0; u < 4; u++)
2791
for (unsigned int tu = 0; tu < 4; tu++)
2793
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2794
}// end loop over 'tu'
2795
}// end loop over 'u'
2796
}// end loop over 't'
2799
}// end loop over 's'
2800
for (unsigned int s = 0; s < 4; s++)
2802
for (unsigned int t = 0; t < 4; t++)
2804
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2805
}// end loop over 't'
2806
}// end loop over 's'
2807
}// end loop over 'r'
2809
// Transform derivatives back to physical element
2810
for (unsigned int r = 0; r < num_derivatives; r++)
2812
for (unsigned int s = 0; s < num_derivatives; s++)
2814
values[r] += transform[r][s]*derivatives[s];
2815
}// end loop over 's'
2816
}// end loop over 'r'
2818
// Delete pointer to array of derivatives on FIAT element
2819
delete [] derivatives;
2821
// Delete pointer to array of combinations of derivatives and transform
2822
for (unsigned int r = 0; r < num_derivatives; r++)
2824
delete [] combinations[r];
2825
}// end loop over 'r'
2826
delete [] combinations;
2827
for (unsigned int r = 0; r < num_derivatives; r++)
2829
delete [] transform[r];
2830
}// end loop over 'r'
2831
delete [] transform;
2837
// Array of basisvalues
2838
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
2840
// Declare helper variables
2841
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2843
// Compute basisvalues
2844
basisvalues[0] = 1.0;
2845
basisvalues[1] = tmp0;
2846
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2847
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2848
basisvalues[0] *= std::sqrt(0.75);
2849
basisvalues[3] *= std::sqrt(1.25);
2850
basisvalues[2] *= std::sqrt(2.5);
2851
basisvalues[1] *= std::sqrt(7.5);
2853
// Table(s) of coefficients
2854
static const double coefficients0[4] = \
2855
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
2857
// Tables of derivatives of the polynomial base (transpose).
2858
static const double dmats0[4][4] = \
2859
{{0.0, 0.0, 0.0, 0.0},
2860
{6.32455532033676, 0.0, 0.0, 0.0},
2861
{0.0, 0.0, 0.0, 0.0},
2862
{0.0, 0.0, 0.0, 0.0}};
2864
static const double dmats1[4][4] = \
2865
{{0.0, 0.0, 0.0, 0.0},
2866
{3.16227766016838, 0.0, 0.0, 0.0},
2867
{5.47722557505166, 0.0, 0.0, 0.0},
2868
{0.0, 0.0, 0.0, 0.0}};
2870
static const double dmats2[4][4] = \
2871
{{0.0, 0.0, 0.0, 0.0},
2872
{3.16227766016838, 0.0, 0.0, 0.0},
2873
{1.82574185835055, 0.0, 0.0, 0.0},
2874
{5.16397779494322, 0.0, 0.0, 0.0}};
2876
// Compute reference derivatives.
2877
// Declare pointer to array of derivatives on FIAT element.
2878
double *derivatives = new double[num_derivatives];
2879
for (unsigned int r = 0; r < num_derivatives; r++)
2881
derivatives[r] = 0.0;
2882
}// end loop over 'r'
2884
// Declare derivative matrix (of polynomial basis).
2885
double dmats[4][4] = \
2886
{{1.0, 0.0, 0.0, 0.0},
2887
{0.0, 1.0, 0.0, 0.0},
2888
{0.0, 0.0, 1.0, 0.0},
2889
{0.0, 0.0, 0.0, 1.0}};
2891
// Declare (auxiliary) derivative matrix (of polynomial basis).
2892
double dmats_old[4][4] = \
2893
{{1.0, 0.0, 0.0, 0.0},
2894
{0.0, 1.0, 0.0, 0.0},
2895
{0.0, 0.0, 1.0, 0.0},
2896
{0.0, 0.0, 0.0, 1.0}};
2898
// Loop possible derivatives.
2899
for (unsigned int r = 0; r < num_derivatives; r++)
2901
// Resetting dmats values to compute next derivative.
2902
for (unsigned int t = 0; t < 4; t++)
2904
for (unsigned int u = 0; u < 4; u++)
2912
}// end loop over 'u'
2913
}// end loop over 't'
2915
// Looping derivative order to generate dmats.
2916
for (unsigned int s = 0; s < n; s++)
2918
// Updating dmats_old with new values and resetting dmats.
2919
for (unsigned int t = 0; t < 4; t++)
2921
for (unsigned int u = 0; u < 4; u++)
2923
dmats_old[t][u] = dmats[t][u];
2925
}// end loop over 'u'
2926
}// end loop over 't'
2928
// Update dmats using an inner product.
2929
if (combinations[r][s] == 0)
2931
for (unsigned int t = 0; t < 4; t++)
2933
for (unsigned int u = 0; u < 4; u++)
2935
for (unsigned int tu = 0; tu < 4; tu++)
2937
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2938
}// end loop over 'tu'
2939
}// end loop over 'u'
2940
}// end loop over 't'
2943
if (combinations[r][s] == 1)
2945
for (unsigned int t = 0; t < 4; t++)
2947
for (unsigned int u = 0; u < 4; u++)
2949
for (unsigned int tu = 0; tu < 4; tu++)
2951
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2952
}// end loop over 'tu'
2953
}// end loop over 'u'
2954
}// end loop over 't'
2957
if (combinations[r][s] == 2)
2959
for (unsigned int t = 0; t < 4; t++)
2961
for (unsigned int u = 0; u < 4; u++)
2963
for (unsigned int tu = 0; tu < 4; tu++)
2965
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2966
}// end loop over 'tu'
2967
}// end loop over 'u'
2968
}// end loop over 't'
2971
}// end loop over 's'
2972
for (unsigned int s = 0; s < 4; s++)
2974
for (unsigned int t = 0; t < 4; t++)
2976
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2977
}// end loop over 't'
2978
}// end loop over 's'
2979
}// end loop over 'r'
2981
// Transform derivatives back to physical element
2982
for (unsigned int r = 0; r < num_derivatives; r++)
2984
for (unsigned int s = 0; s < num_derivatives; s++)
2986
values[r] += transform[r][s]*derivatives[s];
2987
}// end loop over 's'
2988
}// end loop over 'r'
2990
// Delete pointer to array of derivatives on FIAT element
2991
delete [] derivatives;
2993
// Delete pointer to array of combinations of derivatives and transform
2994
for (unsigned int r = 0; r < num_derivatives; r++)
2996
delete [] combinations[r];
2997
}// end loop over 'r'
2998
delete [] combinations;
2999
for (unsigned int r = 0; r < num_derivatives; r++)
3001
delete [] transform[r];
3002
}// end loop over 'r'
3003
delete [] transform;
3009
// Array of basisvalues
3010
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3012
// Declare helper variables
3013
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3015
// Compute basisvalues
3016
basisvalues[0] = 1.0;
3017
basisvalues[1] = tmp0;
3018
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3019
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3020
basisvalues[0] *= std::sqrt(0.75);
3021
basisvalues[3] *= std::sqrt(1.25);
3022
basisvalues[2] *= std::sqrt(2.5);
3023
basisvalues[1] *= std::sqrt(7.5);
3025
// Table(s) of coefficients
3026
static const double coefficients0[4] = \
3027
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
3029
// Tables of derivatives of the polynomial base (transpose).
3030
static const double dmats0[4][4] = \
3031
{{0.0, 0.0, 0.0, 0.0},
3032
{6.32455532033676, 0.0, 0.0, 0.0},
3033
{0.0, 0.0, 0.0, 0.0},
3034
{0.0, 0.0, 0.0, 0.0}};
3036
static const double dmats1[4][4] = \
3037
{{0.0, 0.0, 0.0, 0.0},
3038
{3.16227766016838, 0.0, 0.0, 0.0},
3039
{5.47722557505166, 0.0, 0.0, 0.0},
3040
{0.0, 0.0, 0.0, 0.0}};
3042
static const double dmats2[4][4] = \
3043
{{0.0, 0.0, 0.0, 0.0},
3044
{3.16227766016838, 0.0, 0.0, 0.0},
3045
{1.82574185835055, 0.0, 0.0, 0.0},
3046
{5.16397779494322, 0.0, 0.0, 0.0}};
3048
// Compute reference derivatives.
3049
// Declare pointer to array of derivatives on FIAT element.
3050
double *derivatives = new double[num_derivatives];
3051
for (unsigned int r = 0; r < num_derivatives; r++)
3053
derivatives[r] = 0.0;
3054
}// end loop over 'r'
3056
// Declare derivative matrix (of polynomial basis).
3057
double dmats[4][4] = \
3058
{{1.0, 0.0, 0.0, 0.0},
3059
{0.0, 1.0, 0.0, 0.0},
3060
{0.0, 0.0, 1.0, 0.0},
3061
{0.0, 0.0, 0.0, 1.0}};
3063
// Declare (auxiliary) derivative matrix (of polynomial basis).
3064
double dmats_old[4][4] = \
3065
{{1.0, 0.0, 0.0, 0.0},
3066
{0.0, 1.0, 0.0, 0.0},
3067
{0.0, 0.0, 1.0, 0.0},
3068
{0.0, 0.0, 0.0, 1.0}};
3070
// Loop possible derivatives.
3071
for (unsigned int r = 0; r < num_derivatives; r++)
3073
// Resetting dmats values to compute next derivative.
3074
for (unsigned int t = 0; t < 4; t++)
3076
for (unsigned int u = 0; u < 4; u++)
3084
}// end loop over 'u'
3085
}// end loop over 't'
3087
// Looping derivative order to generate dmats.
3088
for (unsigned int s = 0; s < n; s++)
3090
// Updating dmats_old with new values and resetting dmats.
3091
for (unsigned int t = 0; t < 4; t++)
3093
for (unsigned int u = 0; u < 4; u++)
3095
dmats_old[t][u] = dmats[t][u];
3097
}// end loop over 'u'
3098
}// end loop over 't'
3100
// Update dmats using an inner product.
3101
if (combinations[r][s] == 0)
3103
for (unsigned int t = 0; t < 4; t++)
3105
for (unsigned int u = 0; u < 4; u++)
3107
for (unsigned int tu = 0; tu < 4; tu++)
3109
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3110
}// end loop over 'tu'
3111
}// end loop over 'u'
3112
}// end loop over 't'
3115
if (combinations[r][s] == 1)
3117
for (unsigned int t = 0; t < 4; t++)
3119
for (unsigned int u = 0; u < 4; u++)
3121
for (unsigned int tu = 0; tu < 4; tu++)
3123
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3124
}// end loop over 'tu'
3125
}// end loop over 'u'
3126
}// end loop over 't'
3129
if (combinations[r][s] == 2)
3131
for (unsigned int t = 0; t < 4; t++)
3133
for (unsigned int u = 0; u < 4; u++)
3135
for (unsigned int tu = 0; tu < 4; tu++)
3137
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3138
}// end loop over 'tu'
3139
}// end loop over 'u'
3140
}// end loop over 't'
3143
}// end loop over 's'
3144
for (unsigned int s = 0; s < 4; s++)
3146
for (unsigned int t = 0; t < 4; t++)
3148
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3149
}// end loop over 't'
3150
}// end loop over 's'
3151
}// end loop over 'r'
3153
// Transform derivatives back to physical element
3154
for (unsigned int r = 0; r < num_derivatives; r++)
3156
for (unsigned int s = 0; s < num_derivatives; s++)
3158
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3159
}// end loop over 's'
3160
}// end loop over 'r'
3162
// Delete pointer to array of derivatives on FIAT element
3163
delete [] derivatives;
3165
// Delete pointer to array of combinations of derivatives and transform
3166
for (unsigned int r = 0; r < num_derivatives; r++)
3168
delete [] combinations[r];
3169
}// end loop over 'r'
3170
delete [] combinations;
3171
for (unsigned int r = 0; r < num_derivatives; r++)
3173
delete [] transform[r];
3174
}// end loop over 'r'
3175
delete [] transform;
3181
// Array of basisvalues
3182
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3184
// Declare helper variables
3185
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3187
// Compute basisvalues
3188
basisvalues[0] = 1.0;
3189
basisvalues[1] = tmp0;
3190
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3191
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3192
basisvalues[0] *= std::sqrt(0.75);
3193
basisvalues[3] *= std::sqrt(1.25);
3194
basisvalues[2] *= std::sqrt(2.5);
3195
basisvalues[1] *= std::sqrt(7.5);
3197
// Table(s) of coefficients
3198
static const double coefficients0[4] = \
3199
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
3201
// Tables of derivatives of the polynomial base (transpose).
3202
static const double dmats0[4][4] = \
3203
{{0.0, 0.0, 0.0, 0.0},
3204
{6.32455532033676, 0.0, 0.0, 0.0},
3205
{0.0, 0.0, 0.0, 0.0},
3206
{0.0, 0.0, 0.0, 0.0}};
3208
static const double dmats1[4][4] = \
3209
{{0.0, 0.0, 0.0, 0.0},
3210
{3.16227766016838, 0.0, 0.0, 0.0},
3211
{5.47722557505166, 0.0, 0.0, 0.0},
3212
{0.0, 0.0, 0.0, 0.0}};
3214
static const double dmats2[4][4] = \
3215
{{0.0, 0.0, 0.0, 0.0},
3216
{3.16227766016838, 0.0, 0.0, 0.0},
3217
{1.82574185835055, 0.0, 0.0, 0.0},
3218
{5.16397779494322, 0.0, 0.0, 0.0}};
3220
// Compute reference derivatives.
3221
// Declare pointer to array of derivatives on FIAT element.
3222
double *derivatives = new double[num_derivatives];
3223
for (unsigned int r = 0; r < num_derivatives; r++)
3225
derivatives[r] = 0.0;
3226
}// end loop over 'r'
3228
// Declare derivative matrix (of polynomial basis).
3229
double dmats[4][4] = \
3230
{{1.0, 0.0, 0.0, 0.0},
3231
{0.0, 1.0, 0.0, 0.0},
3232
{0.0, 0.0, 1.0, 0.0},
3233
{0.0, 0.0, 0.0, 1.0}};
3235
// Declare (auxiliary) derivative matrix (of polynomial basis).
3236
double dmats_old[4][4] = \
3237
{{1.0, 0.0, 0.0, 0.0},
3238
{0.0, 1.0, 0.0, 0.0},
3239
{0.0, 0.0, 1.0, 0.0},
3240
{0.0, 0.0, 0.0, 1.0}};
3242
// Loop possible derivatives.
3243
for (unsigned int r = 0; r < num_derivatives; r++)
3245
// Resetting dmats values to compute next derivative.
3246
for (unsigned int t = 0; t < 4; t++)
3248
for (unsigned int u = 0; u < 4; u++)
3256
}// end loop over 'u'
3257
}// end loop over 't'
3259
// Looping derivative order to generate dmats.
3260
for (unsigned int s = 0; s < n; s++)
3262
// Updating dmats_old with new values and resetting dmats.
3263
for (unsigned int t = 0; t < 4; t++)
3265
for (unsigned int u = 0; u < 4; u++)
3267
dmats_old[t][u] = dmats[t][u];
3269
}// end loop over 'u'
3270
}// end loop over 't'
3272
// Update dmats using an inner product.
3273
if (combinations[r][s] == 0)
3275
for (unsigned int t = 0; t < 4; t++)
3277
for (unsigned int u = 0; u < 4; u++)
3279
for (unsigned int tu = 0; tu < 4; tu++)
3281
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3282
}// end loop over 'tu'
3283
}// end loop over 'u'
3284
}// end loop over 't'
3287
if (combinations[r][s] == 1)
3289
for (unsigned int t = 0; t < 4; t++)
3291
for (unsigned int u = 0; u < 4; u++)
3293
for (unsigned int tu = 0; tu < 4; tu++)
3295
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3296
}// end loop over 'tu'
3297
}// end loop over 'u'
3298
}// end loop over 't'
3301
if (combinations[r][s] == 2)
3303
for (unsigned int t = 0; t < 4; t++)
3305
for (unsigned int u = 0; u < 4; u++)
3307
for (unsigned int tu = 0; tu < 4; tu++)
3309
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3310
}// end loop over 'tu'
3311
}// end loop over 'u'
3312
}// end loop over 't'
3315
}// end loop over 's'
3316
for (unsigned int s = 0; s < 4; s++)
3318
for (unsigned int t = 0; t < 4; t++)
3320
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3321
}// end loop over 't'
3322
}// end loop over 's'
3323
}// end loop over 'r'
3325
// Transform derivatives back to physical element
3326
for (unsigned int r = 0; r < num_derivatives; r++)
3328
for (unsigned int s = 0; s < num_derivatives; s++)
3330
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3331
}// end loop over 's'
3332
}// end loop over 'r'
3334
// Delete pointer to array of derivatives on FIAT element
3335
delete [] derivatives;
3337
// Delete pointer to array of combinations of derivatives and transform
3338
for (unsigned int r = 0; r < num_derivatives; r++)
3340
delete [] combinations[r];
3341
}// end loop over 'r'
3342
delete [] combinations;
3343
for (unsigned int r = 0; r < num_derivatives; r++)
3345
delete [] transform[r];
3346
}// end loop over 'r'
3347
delete [] transform;
3353
// Array of basisvalues
3354
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3356
// Declare helper variables
3357
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3359
// Compute basisvalues
3360
basisvalues[0] = 1.0;
3361
basisvalues[1] = tmp0;
3362
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3363
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3364
basisvalues[0] *= std::sqrt(0.75);
3365
basisvalues[3] *= std::sqrt(1.25);
3366
basisvalues[2] *= std::sqrt(2.5);
3367
basisvalues[1] *= std::sqrt(7.5);
3369
// Table(s) of coefficients
3370
static const double coefficients0[4] = \
3371
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
3373
// Tables of derivatives of the polynomial base (transpose).
3374
static const double dmats0[4][4] = \
3375
{{0.0, 0.0, 0.0, 0.0},
3376
{6.32455532033676, 0.0, 0.0, 0.0},
3377
{0.0, 0.0, 0.0, 0.0},
3378
{0.0, 0.0, 0.0, 0.0}};
3380
static const double dmats1[4][4] = \
3381
{{0.0, 0.0, 0.0, 0.0},
3382
{3.16227766016838, 0.0, 0.0, 0.0},
3383
{5.47722557505166, 0.0, 0.0, 0.0},
3384
{0.0, 0.0, 0.0, 0.0}};
3386
static const double dmats2[4][4] = \
3387
{{0.0, 0.0, 0.0, 0.0},
3388
{3.16227766016838, 0.0, 0.0, 0.0},
3389
{1.82574185835055, 0.0, 0.0, 0.0},
3390
{5.16397779494322, 0.0, 0.0, 0.0}};
3392
// Compute reference derivatives.
3393
// Declare pointer to array of derivatives on FIAT element.
3394
double *derivatives = new double[num_derivatives];
3395
for (unsigned int r = 0; r < num_derivatives; r++)
3397
derivatives[r] = 0.0;
3398
}// end loop over 'r'
3400
// Declare derivative matrix (of polynomial basis).
3401
double dmats[4][4] = \
3402
{{1.0, 0.0, 0.0, 0.0},
3403
{0.0, 1.0, 0.0, 0.0},
3404
{0.0, 0.0, 1.0, 0.0},
3405
{0.0, 0.0, 0.0, 1.0}};
3407
// Declare (auxiliary) derivative matrix (of polynomial basis).
3408
double dmats_old[4][4] = \
3409
{{1.0, 0.0, 0.0, 0.0},
3410
{0.0, 1.0, 0.0, 0.0},
3411
{0.0, 0.0, 1.0, 0.0},
3412
{0.0, 0.0, 0.0, 1.0}};
3414
// Loop possible derivatives.
3415
for (unsigned int r = 0; r < num_derivatives; r++)
3417
// Resetting dmats values to compute next derivative.
3418
for (unsigned int t = 0; t < 4; t++)
3420
for (unsigned int u = 0; u < 4; u++)
3428
}// end loop over 'u'
3429
}// end loop over 't'
3431
// Looping derivative order to generate dmats.
3432
for (unsigned int s = 0; s < n; s++)
3434
// Updating dmats_old with new values and resetting dmats.
3435
for (unsigned int t = 0; t < 4; t++)
3437
for (unsigned int u = 0; u < 4; u++)
3439
dmats_old[t][u] = dmats[t][u];
3441
}// end loop over 'u'
3442
}// end loop over 't'
3444
// Update dmats using an inner product.
3445
if (combinations[r][s] == 0)
3447
for (unsigned int t = 0; t < 4; t++)
3449
for (unsigned int u = 0; u < 4; u++)
3451
for (unsigned int tu = 0; tu < 4; tu++)
3453
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3454
}// end loop over 'tu'
3455
}// end loop over 'u'
3456
}// end loop over 't'
3459
if (combinations[r][s] == 1)
3461
for (unsigned int t = 0; t < 4; t++)
3463
for (unsigned int u = 0; u < 4; u++)
3465
for (unsigned int tu = 0; tu < 4; tu++)
3467
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3468
}// end loop over 'tu'
3469
}// end loop over 'u'
3470
}// end loop over 't'
3473
if (combinations[r][s] == 2)
3475
for (unsigned int t = 0; t < 4; t++)
3477
for (unsigned int u = 0; u < 4; u++)
3479
for (unsigned int tu = 0; tu < 4; tu++)
3481
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3482
}// end loop over 'tu'
3483
}// end loop over 'u'
3484
}// end loop over 't'
3487
}// end loop over 's'
3488
for (unsigned int s = 0; s < 4; s++)
3490
for (unsigned int t = 0; t < 4; t++)
3492
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3493
}// end loop over 't'
3494
}// end loop over 's'
3495
}// end loop over 'r'
3497
// Transform derivatives back to physical element
3498
for (unsigned int r = 0; r < num_derivatives; r++)
3500
for (unsigned int s = 0; s < num_derivatives; s++)
3502
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3503
}// end loop over 's'
3504
}// end loop over 'r'
3506
// Delete pointer to array of derivatives on FIAT element
3507
delete [] derivatives;
3509
// Delete pointer to array of combinations of derivatives and transform
3510
for (unsigned int r = 0; r < num_derivatives; r++)
3512
delete [] combinations[r];
3513
}// end loop over 'r'
3514
delete [] combinations;
3515
for (unsigned int r = 0; r < num_derivatives; r++)
3517
delete [] transform[r];
3518
}// end loop over 'r'
3519
delete [] transform;
3525
// Array of basisvalues
3526
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3528
// Declare helper variables
3529
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3531
// Compute basisvalues
3532
basisvalues[0] = 1.0;
3533
basisvalues[1] = tmp0;
3534
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3535
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3536
basisvalues[0] *= std::sqrt(0.75);
3537
basisvalues[3] *= std::sqrt(1.25);
3538
basisvalues[2] *= std::sqrt(2.5);
3539
basisvalues[1] *= std::sqrt(7.5);
3541
// Table(s) of coefficients
3542
static const double coefficients0[4] = \
3543
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
3545
// Tables of derivatives of the polynomial base (transpose).
3546
static const double dmats0[4][4] = \
3547
{{0.0, 0.0, 0.0, 0.0},
3548
{6.32455532033676, 0.0, 0.0, 0.0},
3549
{0.0, 0.0, 0.0, 0.0},
3550
{0.0, 0.0, 0.0, 0.0}};
3552
static const double dmats1[4][4] = \
3553
{{0.0, 0.0, 0.0, 0.0},
3554
{3.16227766016838, 0.0, 0.0, 0.0},
3555
{5.47722557505166, 0.0, 0.0, 0.0},
3556
{0.0, 0.0, 0.0, 0.0}};
3558
static const double dmats2[4][4] = \
3559
{{0.0, 0.0, 0.0, 0.0},
3560
{3.16227766016838, 0.0, 0.0, 0.0},
3561
{1.82574185835055, 0.0, 0.0, 0.0},
3562
{5.16397779494322, 0.0, 0.0, 0.0}};
3564
// Compute reference derivatives.
3565
// Declare pointer to array of derivatives on FIAT element.
3566
double *derivatives = new double[num_derivatives];
3567
for (unsigned int r = 0; r < num_derivatives; r++)
3569
derivatives[r] = 0.0;
3570
}// end loop over 'r'
3572
// Declare derivative matrix (of polynomial basis).
3573
double dmats[4][4] = \
3574
{{1.0, 0.0, 0.0, 0.0},
3575
{0.0, 1.0, 0.0, 0.0},
3576
{0.0, 0.0, 1.0, 0.0},
3577
{0.0, 0.0, 0.0, 1.0}};
3579
// Declare (auxiliary) derivative matrix (of polynomial basis).
3580
double dmats_old[4][4] = \
3581
{{1.0, 0.0, 0.0, 0.0},
3582
{0.0, 1.0, 0.0, 0.0},
3583
{0.0, 0.0, 1.0, 0.0},
3584
{0.0, 0.0, 0.0, 1.0}};
3586
// Loop possible derivatives.
3587
for (unsigned int r = 0; r < num_derivatives; r++)
3589
// Resetting dmats values to compute next derivative.
3590
for (unsigned int t = 0; t < 4; t++)
3592
for (unsigned int u = 0; u < 4; u++)
3600
}// end loop over 'u'
3601
}// end loop over 't'
3603
// Looping derivative order to generate dmats.
3604
for (unsigned int s = 0; s < n; s++)
3606
// Updating dmats_old with new values and resetting dmats.
3607
for (unsigned int t = 0; t < 4; t++)
3609
for (unsigned int u = 0; u < 4; u++)
3611
dmats_old[t][u] = dmats[t][u];
3613
}// end loop over 'u'
3614
}// end loop over 't'
3616
// Update dmats using an inner product.
3617
if (combinations[r][s] == 0)
3619
for (unsigned int t = 0; t < 4; t++)
3621
for (unsigned int u = 0; u < 4; u++)
3623
for (unsigned int tu = 0; tu < 4; tu++)
3625
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3626
}// end loop over 'tu'
3627
}// end loop over 'u'
3628
}// end loop over 't'
3631
if (combinations[r][s] == 1)
3633
for (unsigned int t = 0; t < 4; t++)
3635
for (unsigned int u = 0; u < 4; u++)
3637
for (unsigned int tu = 0; tu < 4; tu++)
3639
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3640
}// end loop over 'tu'
3641
}// end loop over 'u'
3642
}// end loop over 't'
3645
if (combinations[r][s] == 2)
3647
for (unsigned int t = 0; t < 4; t++)
3649
for (unsigned int u = 0; u < 4; u++)
3651
for (unsigned int tu = 0; tu < 4; tu++)
3653
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3654
}// end loop over 'tu'
3655
}// end loop over 'u'
3656
}// end loop over 't'
3659
}// end loop over 's'
3660
for (unsigned int s = 0; s < 4; s++)
3662
for (unsigned int t = 0; t < 4; t++)
3664
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3665
}// end loop over 't'
3666
}// end loop over 's'
3667
}// end loop over 'r'
3669
// Transform derivatives back to physical element
3670
for (unsigned int r = 0; r < num_derivatives; r++)
3672
for (unsigned int s = 0; s < num_derivatives; s++)
3674
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3675
}// end loop over 's'
3676
}// end loop over 'r'
3678
// Delete pointer to array of derivatives on FIAT element
3679
delete [] derivatives;
3681
// Delete pointer to array of combinations of derivatives and transform
3682
for (unsigned int r = 0; r < num_derivatives; r++)
3684
delete [] combinations[r];
3685
}// end loop over 'r'
3686
delete [] combinations;
3687
for (unsigned int r = 0; r < num_derivatives; r++)
3689
delete [] transform[r];
3690
}// end loop over 'r'
3691
delete [] transform;
3697
// Array of basisvalues
3698
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3700
// Declare helper variables
3701
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3703
// Compute basisvalues
3704
basisvalues[0] = 1.0;
3705
basisvalues[1] = tmp0;
3706
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3707
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3708
basisvalues[0] *= std::sqrt(0.75);
3709
basisvalues[3] *= std::sqrt(1.25);
3710
basisvalues[2] *= std::sqrt(2.5);
3711
basisvalues[1] *= std::sqrt(7.5);
3713
// Table(s) of coefficients
3714
static const double coefficients0[4] = \
3715
{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
3717
// Tables of derivatives of the polynomial base (transpose).
3718
static const double dmats0[4][4] = \
3719
{{0.0, 0.0, 0.0, 0.0},
3720
{6.32455532033676, 0.0, 0.0, 0.0},
3721
{0.0, 0.0, 0.0, 0.0},
3722
{0.0, 0.0, 0.0, 0.0}};
3724
static const double dmats1[4][4] = \
3725
{{0.0, 0.0, 0.0, 0.0},
3726
{3.16227766016838, 0.0, 0.0, 0.0},
3727
{5.47722557505166, 0.0, 0.0, 0.0},
3728
{0.0, 0.0, 0.0, 0.0}};
3730
static const double dmats2[4][4] = \
3731
{{0.0, 0.0, 0.0, 0.0},
3732
{3.16227766016838, 0.0, 0.0, 0.0},
3733
{1.82574185835055, 0.0, 0.0, 0.0},
3734
{5.16397779494322, 0.0, 0.0, 0.0}};
3736
// Compute reference derivatives.
3737
// Declare pointer to array of derivatives on FIAT element.
3738
double *derivatives = new double[num_derivatives];
3739
for (unsigned int r = 0; r < num_derivatives; r++)
3741
derivatives[r] = 0.0;
3742
}// end loop over 'r'
3744
// Declare derivative matrix (of polynomial basis).
3745
double dmats[4][4] = \
3746
{{1.0, 0.0, 0.0, 0.0},
3747
{0.0, 1.0, 0.0, 0.0},
3748
{0.0, 0.0, 1.0, 0.0},
3749
{0.0, 0.0, 0.0, 1.0}};
3751
// Declare (auxiliary) derivative matrix (of polynomial basis).
3752
double dmats_old[4][4] = \
3753
{{1.0, 0.0, 0.0, 0.0},
3754
{0.0, 1.0, 0.0, 0.0},
3755
{0.0, 0.0, 1.0, 0.0},
3756
{0.0, 0.0, 0.0, 1.0}};
3758
// Loop possible derivatives.
3759
for (unsigned int r = 0; r < num_derivatives; r++)
3761
// Resetting dmats values to compute next derivative.
3762
for (unsigned int t = 0; t < 4; t++)
3764
for (unsigned int u = 0; u < 4; u++)
3772
}// end loop over 'u'
3773
}// end loop over 't'
3775
// Looping derivative order to generate dmats.
3776
for (unsigned int s = 0; s < n; s++)
3778
// Updating dmats_old with new values and resetting dmats.
3779
for (unsigned int t = 0; t < 4; t++)
3781
for (unsigned int u = 0; u < 4; u++)
3783
dmats_old[t][u] = dmats[t][u];
3785
}// end loop over 'u'
3786
}// end loop over 't'
3788
// Update dmats using an inner product.
3789
if (combinations[r][s] == 0)
3791
for (unsigned int t = 0; t < 4; t++)
3793
for (unsigned int u = 0; u < 4; u++)
3795
for (unsigned int tu = 0; tu < 4; tu++)
3797
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3798
}// end loop over 'tu'
3799
}// end loop over 'u'
3800
}// end loop over 't'
3803
if (combinations[r][s] == 1)
3805
for (unsigned int t = 0; t < 4; t++)
3807
for (unsigned int u = 0; u < 4; u++)
3809
for (unsigned int tu = 0; tu < 4; tu++)
3811
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3812
}// end loop over 'tu'
3813
}// end loop over 'u'
3814
}// end loop over 't'
3817
if (combinations[r][s] == 2)
3819
for (unsigned int t = 0; t < 4; t++)
3821
for (unsigned int u = 0; u < 4; u++)
3823
for (unsigned int tu = 0; tu < 4; tu++)
3825
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3826
}// end loop over 'tu'
3827
}// end loop over 'u'
3828
}// end loop over 't'
3831
}// end loop over 's'
3832
for (unsigned int s = 0; s < 4; s++)
3834
for (unsigned int t = 0; t < 4; t++)
3836
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3837
}// end loop over 't'
3838
}// end loop over 's'
3839
}// end loop over 'r'
3841
// Transform derivatives back to physical element
3842
for (unsigned int r = 0; r < num_derivatives; r++)
3844
for (unsigned int s = 0; s < num_derivatives; s++)
3846
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
3847
}// end loop over 's'
3848
}// end loop over 'r'
3850
// Delete pointer to array of derivatives on FIAT element
3851
delete [] derivatives;
3853
// Delete pointer to array of combinations of derivatives and transform
3854
for (unsigned int r = 0; r < num_derivatives; r++)
3856
delete [] combinations[r];
3857
}// end loop over 'r'
3858
delete [] combinations;
3859
for (unsigned int r = 0; r < num_derivatives; r++)
3861
delete [] transform[r];
3862
}// end loop over 'r'
3863
delete [] transform;
3869
// Array of basisvalues
3870
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
3872
// Declare helper variables
3873
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
3875
// Compute basisvalues
3876
basisvalues[0] = 1.0;
3877
basisvalues[1] = tmp0;
3878
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
3879
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
3880
basisvalues[0] *= std::sqrt(0.75);
3881
basisvalues[3] *= std::sqrt(1.25);
3882
basisvalues[2] *= std::sqrt(2.5);
3883
basisvalues[1] *= std::sqrt(7.5);
3885
// Table(s) of coefficients
3886
static const double coefficients0[4] = \
3887
{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
3889
// Tables of derivatives of the polynomial base (transpose).
3890
static const double dmats0[4][4] = \
3891
{{0.0, 0.0, 0.0, 0.0},
3892
{6.32455532033676, 0.0, 0.0, 0.0},
3893
{0.0, 0.0, 0.0, 0.0},
3894
{0.0, 0.0, 0.0, 0.0}};
3896
static const double dmats1[4][4] = \
3897
{{0.0, 0.0, 0.0, 0.0},
3898
{3.16227766016838, 0.0, 0.0, 0.0},
3899
{5.47722557505166, 0.0, 0.0, 0.0},
3900
{0.0, 0.0, 0.0, 0.0}};
3902
static const double dmats2[4][4] = \
3903
{{0.0, 0.0, 0.0, 0.0},
3904
{3.16227766016838, 0.0, 0.0, 0.0},
3905
{1.82574185835055, 0.0, 0.0, 0.0},
3906
{5.16397779494322, 0.0, 0.0, 0.0}};
3908
// Compute reference derivatives.
3909
// Declare pointer to array of derivatives on FIAT element.
3910
double *derivatives = new double[num_derivatives];
3911
for (unsigned int r = 0; r < num_derivatives; r++)
3913
derivatives[r] = 0.0;
3914
}// end loop over 'r'
3916
// Declare derivative matrix (of polynomial basis).
3917
double dmats[4][4] = \
3918
{{1.0, 0.0, 0.0, 0.0},
3919
{0.0, 1.0, 0.0, 0.0},
3920
{0.0, 0.0, 1.0, 0.0},
3921
{0.0, 0.0, 0.0, 1.0}};
3923
// Declare (auxiliary) derivative matrix (of polynomial basis).
3924
double dmats_old[4][4] = \
3925
{{1.0, 0.0, 0.0, 0.0},
3926
{0.0, 1.0, 0.0, 0.0},
3927
{0.0, 0.0, 1.0, 0.0},
3928
{0.0, 0.0, 0.0, 1.0}};
3930
// Loop possible derivatives.
3931
for (unsigned int r = 0; r < num_derivatives; r++)
3933
// Resetting dmats values to compute next derivative.
3934
for (unsigned int t = 0; t < 4; t++)
3936
for (unsigned int u = 0; u < 4; u++)
3944
}// end loop over 'u'
3945
}// end loop over 't'
3947
// Looping derivative order to generate dmats.
3948
for (unsigned int s = 0; s < n; s++)
3950
// Updating dmats_old with new values and resetting dmats.
3951
for (unsigned int t = 0; t < 4; t++)
3953
for (unsigned int u = 0; u < 4; u++)
3955
dmats_old[t][u] = dmats[t][u];
3957
}// end loop over 'u'
3958
}// end loop over 't'
3960
// Update dmats using an inner product.
3961
if (combinations[r][s] == 0)
3963
for (unsigned int t = 0; t < 4; t++)
3965
for (unsigned int u = 0; u < 4; u++)
3967
for (unsigned int tu = 0; tu < 4; tu++)
3969
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3970
}// end loop over 'tu'
3971
}// end loop over 'u'
3972
}// end loop over 't'
3975
if (combinations[r][s] == 1)
3977
for (unsigned int t = 0; t < 4; t++)
3979
for (unsigned int u = 0; u < 4; u++)
3981
for (unsigned int tu = 0; tu < 4; tu++)
3983
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3984
}// end loop over 'tu'
3985
}// end loop over 'u'
3986
}// end loop over 't'
3989
if (combinations[r][s] == 2)
3991
for (unsigned int t = 0; t < 4; t++)
3993
for (unsigned int u = 0; u < 4; u++)
3995
for (unsigned int tu = 0; tu < 4; tu++)
3997
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3998
}// end loop over 'tu'
3999
}// end loop over 'u'
4000
}// end loop over 't'
4003
}// end loop over 's'
4004
for (unsigned int s = 0; s < 4; s++)
4006
for (unsigned int t = 0; t < 4; t++)
4008
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4009
}// end loop over 't'
4010
}// end loop over 's'
4011
}// end loop over 'r'
4013
// Transform derivatives back to physical element
4014
for (unsigned int r = 0; r < num_derivatives; r++)
4016
for (unsigned int s = 0; s < num_derivatives; s++)
4018
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
4019
}// end loop over 's'
4020
}// end loop over 'r'
4022
// Delete pointer to array of derivatives on FIAT element
4023
delete [] derivatives;
4025
// Delete pointer to array of combinations of derivatives and transform
4026
for (unsigned int r = 0; r < num_derivatives; r++)
4028
delete [] combinations[r];
4029
}// end loop over 'r'
4030
delete [] combinations;
4031
for (unsigned int r = 0; r < num_derivatives; r++)
4033
delete [] transform[r];
4034
}// end loop over 'r'
4035
delete [] transform;
4041
// Array of basisvalues
4042
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
4044
// Declare helper variables
4045
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
4047
// Compute basisvalues
4048
basisvalues[0] = 1.0;
4049
basisvalues[1] = tmp0;
4050
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
4051
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
4052
basisvalues[0] *= std::sqrt(0.75);
4053
basisvalues[3] *= std::sqrt(1.25);
4054
basisvalues[2] *= std::sqrt(2.5);
4055
basisvalues[1] *= std::sqrt(7.5);
4057
// Table(s) of coefficients
4058
static const double coefficients0[4] = \
4059
{0.288675134594813, 0.0, 0.210818510677892, -0.074535599249993};
4061
// Tables of derivatives of the polynomial base (transpose).
4062
static const double dmats0[4][4] = \
4063
{{0.0, 0.0, 0.0, 0.0},
4064
{6.32455532033676, 0.0, 0.0, 0.0},
4065
{0.0, 0.0, 0.0, 0.0},
4066
{0.0, 0.0, 0.0, 0.0}};
4068
static const double dmats1[4][4] = \
4069
{{0.0, 0.0, 0.0, 0.0},
4070
{3.16227766016838, 0.0, 0.0, 0.0},
4071
{5.47722557505166, 0.0, 0.0, 0.0},
4072
{0.0, 0.0, 0.0, 0.0}};
4074
static const double dmats2[4][4] = \
4075
{{0.0, 0.0, 0.0, 0.0},
4076
{3.16227766016838, 0.0, 0.0, 0.0},
4077
{1.82574185835055, 0.0, 0.0, 0.0},
4078
{5.16397779494322, 0.0, 0.0, 0.0}};
4080
// Compute reference derivatives.
4081
// Declare pointer to array of derivatives on FIAT element.
4082
double *derivatives = new double[num_derivatives];
4083
for (unsigned int r = 0; r < num_derivatives; r++)
4085
derivatives[r] = 0.0;
4086
}// end loop over 'r'
4088
// Declare derivative matrix (of polynomial basis).
4089
double dmats[4][4] = \
4090
{{1.0, 0.0, 0.0, 0.0},
4091
{0.0, 1.0, 0.0, 0.0},
4092
{0.0, 0.0, 1.0, 0.0},
4093
{0.0, 0.0, 0.0, 1.0}};
4095
// Declare (auxiliary) derivative matrix (of polynomial basis).
4096
double dmats_old[4][4] = \
4097
{{1.0, 0.0, 0.0, 0.0},
4098
{0.0, 1.0, 0.0, 0.0},
4099
{0.0, 0.0, 1.0, 0.0},
4100
{0.0, 0.0, 0.0, 1.0}};
4102
// Loop possible derivatives.
4103
for (unsigned int r = 0; r < num_derivatives; r++)
4105
// Resetting dmats values to compute next derivative.
4106
for (unsigned int t = 0; t < 4; t++)
4108
for (unsigned int u = 0; u < 4; u++)
4116
}// end loop over 'u'
4117
}// end loop over 't'
4119
// Looping derivative order to generate dmats.
4120
for (unsigned int s = 0; s < n; s++)
4122
// Updating dmats_old with new values and resetting dmats.
4123
for (unsigned int t = 0; t < 4; t++)
4125
for (unsigned int u = 0; u < 4; u++)
4127
dmats_old[t][u] = dmats[t][u];
4129
}// end loop over 'u'
4130
}// end loop over 't'
4132
// Update dmats using an inner product.
4133
if (combinations[r][s] == 0)
4135
for (unsigned int t = 0; t < 4; t++)
4137
for (unsigned int u = 0; u < 4; u++)
4139
for (unsigned int tu = 0; tu < 4; tu++)
4141
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4142
}// end loop over 'tu'
4143
}// end loop over 'u'
4144
}// end loop over 't'
4147
if (combinations[r][s] == 1)
4149
for (unsigned int t = 0; t < 4; t++)
4151
for (unsigned int u = 0; u < 4; u++)
4153
for (unsigned int tu = 0; tu < 4; tu++)
4155
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4156
}// end loop over 'tu'
4157
}// end loop over 'u'
4158
}// end loop over 't'
4161
if (combinations[r][s] == 2)
4163
for (unsigned int t = 0; t < 4; t++)
4165
for (unsigned int u = 0; u < 4; u++)
4167
for (unsigned int tu = 0; tu < 4; tu++)
4169
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
4170
}// end loop over 'tu'
4171
}// end loop over 'u'
4172
}// end loop over 't'
4175
}// end loop over 's'
4176
for (unsigned int s = 0; s < 4; s++)
4178
for (unsigned int t = 0; t < 4; t++)
4180
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4181
}// end loop over 't'
4182
}// end loop over 's'
4183
}// end loop over 'r'
4185
// Transform derivatives back to physical element
4186
for (unsigned int r = 0; r < num_derivatives; r++)
4188
for (unsigned int s = 0; s < num_derivatives; s++)
4190
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
4191
}// end loop over 's'
4192
}// end loop over 'r'
4194
// Delete pointer to array of derivatives on FIAT element
4195
delete [] derivatives;
4197
// Delete pointer to array of combinations of derivatives and transform
4198
for (unsigned int r = 0; r < num_derivatives; r++)
4200
delete [] combinations[r];
4201
}// end loop over 'r'
4202
delete [] combinations;
4203
for (unsigned int r = 0; r < num_derivatives; r++)
4205
delete [] transform[r];
4206
}// end loop over 'r'
4207
delete [] transform;
4213
// Array of basisvalues
4214
double basisvalues[4] = {0.0, 0.0, 0.0, 0.0};
4216
// Declare helper variables
4217
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
4219
// Compute basisvalues
4220
basisvalues[0] = 1.0;
4221
basisvalues[1] = tmp0;
4222
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
4223
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
4224
basisvalues[0] *= std::sqrt(0.75);
4225
basisvalues[3] *= std::sqrt(1.25);
4226
basisvalues[2] *= std::sqrt(2.5);
4227
basisvalues[1] *= std::sqrt(7.5);
4229
// Table(s) of coefficients
4230
static const double coefficients0[4] = \
4231
{0.288675134594813, 0.0, 0.0, 0.223606797749979};
4233
// Tables of derivatives of the polynomial base (transpose).
4234
static const double dmats0[4][4] = \
4235
{{0.0, 0.0, 0.0, 0.0},
4236
{6.32455532033676, 0.0, 0.0, 0.0},
4237
{0.0, 0.0, 0.0, 0.0},
4238
{0.0, 0.0, 0.0, 0.0}};
4240
static const double dmats1[4][4] = \
4241
{{0.0, 0.0, 0.0, 0.0},
4242
{3.16227766016838, 0.0, 0.0, 0.0},
4243
{5.47722557505166, 0.0, 0.0, 0.0},
4244
{0.0, 0.0, 0.0, 0.0}};
4246
static const double dmats2[4][4] = \
4247
{{0.0, 0.0, 0.0, 0.0},
4248
{3.16227766016838, 0.0, 0.0, 0.0},
4249
{1.82574185835055, 0.0, 0.0, 0.0},
4250
{5.16397779494322, 0.0, 0.0, 0.0}};
4252
// Compute reference derivatives.
4253
// Declare pointer to array of derivatives on FIAT element.
4254
double *derivatives = new double[num_derivatives];
4255
for (unsigned int r = 0; r < num_derivatives; r++)
4257
derivatives[r] = 0.0;
4258
}// end loop over 'r'
4260
// Declare derivative matrix (of polynomial basis).
4261
double dmats[4][4] = \
4262
{{1.0, 0.0, 0.0, 0.0},
4263
{0.0, 1.0, 0.0, 0.0},
4264
{0.0, 0.0, 1.0, 0.0},
4265
{0.0, 0.0, 0.0, 1.0}};
4267
// Declare (auxiliary) derivative matrix (of polynomial basis).
4268
double dmats_old[4][4] = \
4269
{{1.0, 0.0, 0.0, 0.0},
4270
{0.0, 1.0, 0.0, 0.0},
4271
{0.0, 0.0, 1.0, 0.0},
4272
{0.0, 0.0, 0.0, 1.0}};
4274
// Loop possible derivatives.
4275
for (unsigned int r = 0; r < num_derivatives; r++)
4277
// Resetting dmats values to compute next derivative.
4278
for (unsigned int t = 0; t < 4; t++)
4280
for (unsigned int u = 0; u < 4; u++)
4288
}// end loop over 'u'
4289
}// end loop over 't'
4291
// Looping derivative order to generate dmats.
4292
for (unsigned int s = 0; s < n; s++)
4294
// Updating dmats_old with new values and resetting dmats.
4295
for (unsigned int t = 0; t < 4; t++)
4297
for (unsigned int u = 0; u < 4; u++)
4299
dmats_old[t][u] = dmats[t][u];
4301
}// end loop over 'u'
4302
}// end loop over 't'
4304
// Update dmats using an inner product.
4305
if (combinations[r][s] == 0)
4307
for (unsigned int t = 0; t < 4; t++)
4309
for (unsigned int u = 0; u < 4; u++)
4311
for (unsigned int tu = 0; tu < 4; tu++)
4313
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4314
}// end loop over 'tu'
4315
}// end loop over 'u'
4316
}// end loop over 't'
4319
if (combinations[r][s] == 1)
4321
for (unsigned int t = 0; t < 4; t++)
4323
for (unsigned int u = 0; u < 4; u++)
4325
for (unsigned int tu = 0; tu < 4; tu++)
4327
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4328
}// end loop over 'tu'
4329
}// end loop over 'u'
4330
}// end loop over 't'
4333
if (combinations[r][s] == 2)
4335
for (unsigned int t = 0; t < 4; t++)
4337
for (unsigned int u = 0; u < 4; u++)
4339
for (unsigned int tu = 0; tu < 4; tu++)
4341
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
4342
}// end loop over 'tu'
4343
}// end loop over 'u'
4344
}// end loop over 't'
4347
}// end loop over 's'
4348
for (unsigned int s = 0; s < 4; s++)
4350
for (unsigned int t = 0; t < 4; t++)
4352
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4353
}// end loop over 't'
4354
}// end loop over 's'
4355
}// end loop over 'r'
4357
// Transform derivatives back to physical element
4358
for (unsigned int r = 0; r < num_derivatives; r++)
4360
for (unsigned int s = 0; s < num_derivatives; s++)
4362
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
4363
}// end loop over 's'
4364
}// end loop over 'r'
4366
// Delete pointer to array of derivatives on FIAT element
4367
delete [] derivatives;
4369
// Delete pointer to array of combinations of derivatives and transform
4370
for (unsigned int r = 0; r < num_derivatives; r++)
4372
delete [] combinations[r];
4373
}// end loop over 'r'
4374
delete [] combinations;
4375
for (unsigned int r = 0; r < num_derivatives; r++)
4377
delete [] transform[r];
4378
}// end loop over 'r'
4379
delete [] transform;
4386
/// Evaluate order n derivatives of all basis functions at given point x in cell
4387
virtual void evaluate_basis_derivatives_all(std::size_t n,
4390
const double* vertex_coordinates,
4391
int cell_orientation) const
4393
// Compute number of derivatives.
4394
unsigned int num_derivatives = 1;
4395
for (unsigned int r = 0; r < n; r++)
4397
num_derivatives *= 3;
4398
}// end loop over 'r'
4400
// Helper variable to hold values of a single dof.
4401
double *dof_values = new double[3*num_derivatives];
4402
for (unsigned int r = 0; r < 3*num_derivatives; r++)
4404
dof_values[r] = 0.0;
4405
}// end loop over 'r'
4407
// Loop dofs and call evaluate_basis_derivatives.
4408
for (unsigned int r = 0; r < 12; r++)
4410
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
4411
for (unsigned int s = 0; s < 3*num_derivatives; s++)
4413
values[r*3*num_derivatives + s] = dof_values[s];
4414
}// end loop over 's'
4415
}// end loop over 'r'
4418
delete [] dof_values;
4421
/// Evaluate linear functional for dof i on the function f
4422
virtual double evaluate_dof(std::size_t i,
4423
const ufc::function& f,
4424
const double* vertex_coordinates,
4425
int cell_orientation,
4426
const ufc::cell& c) const
4428
// Declare variables for result of evaluation
4431
// Declare variable for physical coordinates
4437
y[0] = vertex_coordinates[0];
4438
y[1] = vertex_coordinates[1];
4439
y[2] = vertex_coordinates[2];
4440
f.evaluate(vals, y, c);
4446
y[0] = vertex_coordinates[3];
4447
y[1] = vertex_coordinates[4];
4448
y[2] = vertex_coordinates[5];
4449
f.evaluate(vals, y, c);
4455
y[0] = vertex_coordinates[6];
4456
y[1] = vertex_coordinates[7];
4457
y[2] = vertex_coordinates[8];
4458
f.evaluate(vals, y, c);
4464
y[0] = vertex_coordinates[9];
4465
y[1] = vertex_coordinates[10];
4466
y[2] = vertex_coordinates[11];
4467
f.evaluate(vals, y, c);
4473
y[0] = vertex_coordinates[0];
4474
y[1] = vertex_coordinates[1];
4475
y[2] = vertex_coordinates[2];
4476
f.evaluate(vals, y, c);
4482
y[0] = vertex_coordinates[3];
4483
y[1] = vertex_coordinates[4];
4484
y[2] = vertex_coordinates[5];
4485
f.evaluate(vals, y, c);
4491
y[0] = vertex_coordinates[6];
4492
y[1] = vertex_coordinates[7];
4493
y[2] = vertex_coordinates[8];
4494
f.evaluate(vals, y, c);
4500
y[0] = vertex_coordinates[9];
4501
y[1] = vertex_coordinates[10];
4502
y[2] = vertex_coordinates[11];
4503
f.evaluate(vals, y, c);
4509
y[0] = vertex_coordinates[0];
4510
y[1] = vertex_coordinates[1];
4511
y[2] = vertex_coordinates[2];
4512
f.evaluate(vals, y, c);
4518
y[0] = vertex_coordinates[3];
4519
y[1] = vertex_coordinates[4];
4520
y[2] = vertex_coordinates[5];
4521
f.evaluate(vals, y, c);
4527
y[0] = vertex_coordinates[6];
4528
y[1] = vertex_coordinates[7];
4529
y[2] = vertex_coordinates[8];
4530
f.evaluate(vals, y, c);
4536
y[0] = vertex_coordinates[9];
4537
y[1] = vertex_coordinates[10];
4538
y[2] = vertex_coordinates[11];
4539
f.evaluate(vals, y, c);
4548
/// Evaluate linear functionals for all dofs on the function f
4549
virtual void evaluate_dofs(double* values,
4550
const ufc::function& f,
4551
const double* vertex_coordinates,
4552
int cell_orientation,
4553
const ufc::cell& c) const
4555
// Declare variables for result of evaluation
4558
// Declare variable for physical coordinates
4560
y[0] = vertex_coordinates[0];
4561
y[1] = vertex_coordinates[1];
4562
y[2] = vertex_coordinates[2];
4563
f.evaluate(vals, y, c);
4564
values[0] = vals[0];
4565
y[0] = vertex_coordinates[3];
4566
y[1] = vertex_coordinates[4];
4567
y[2] = vertex_coordinates[5];
4568
f.evaluate(vals, y, c);
4569
values[1] = vals[0];
4570
y[0] = vertex_coordinates[6];
4571
y[1] = vertex_coordinates[7];
4572
y[2] = vertex_coordinates[8];
4573
f.evaluate(vals, y, c);
4574
values[2] = vals[0];
4575
y[0] = vertex_coordinates[9];
4576
y[1] = vertex_coordinates[10];
4577
y[2] = vertex_coordinates[11];
4578
f.evaluate(vals, y, c);
4579
values[3] = vals[0];
4580
y[0] = vertex_coordinates[0];
4581
y[1] = vertex_coordinates[1];
4582
y[2] = vertex_coordinates[2];
4583
f.evaluate(vals, y, c);
4584
values[4] = vals[1];
4585
y[0] = vertex_coordinates[3];
4586
y[1] = vertex_coordinates[4];
4587
y[2] = vertex_coordinates[5];
4588
f.evaluate(vals, y, c);
4589
values[5] = vals[1];
4590
y[0] = vertex_coordinates[6];
4591
y[1] = vertex_coordinates[7];
4592
y[2] = vertex_coordinates[8];
4593
f.evaluate(vals, y, c);
4594
values[6] = vals[1];
4595
y[0] = vertex_coordinates[9];
4596
y[1] = vertex_coordinates[10];
4597
y[2] = vertex_coordinates[11];
4598
f.evaluate(vals, y, c);
4599
values[7] = vals[1];
4600
y[0] = vertex_coordinates[0];
4601
y[1] = vertex_coordinates[1];
4602
y[2] = vertex_coordinates[2];
4603
f.evaluate(vals, y, c);
4604
values[8] = vals[2];
4605
y[0] = vertex_coordinates[3];
4606
y[1] = vertex_coordinates[4];
4607
y[2] = vertex_coordinates[5];
4608
f.evaluate(vals, y, c);
4609
values[9] = vals[2];
4610
y[0] = vertex_coordinates[6];
4611
y[1] = vertex_coordinates[7];
4612
y[2] = vertex_coordinates[8];
4613
f.evaluate(vals, y, c);
4614
values[10] = vals[2];
4615
y[0] = vertex_coordinates[9];
4616
y[1] = vertex_coordinates[10];
4617
y[2] = vertex_coordinates[11];
4618
f.evaluate(vals, y, c);
4619
values[11] = vals[2];
4622
/// Interpolate vertex values from dof values
4623
virtual void interpolate_vertex_values(double* vertex_values,
4624
const double* dof_values,
4625
const double* vertex_coordinates,
4626
int cell_orientation,
4627
const ufc::cell& c) const
4629
// Evaluate function and change variables
4630
vertex_values[0] = dof_values[0];
4631
vertex_values[3] = dof_values[1];
4632
vertex_values[6] = dof_values[2];
4633
vertex_values[9] = dof_values[3];
4634
// Evaluate function and change variables
4635
vertex_values[1] = dof_values[4];
4636
vertex_values[4] = dof_values[5];
4637
vertex_values[7] = dof_values[6];
4638
vertex_values[10] = dof_values[7];
4639
// Evaluate function and change variables
4640
vertex_values[2] = dof_values[8];
4641
vertex_values[5] = dof_values[9];
4642
vertex_values[8] = dof_values[10];
4643
vertex_values[11] = dof_values[11];
4646
/// Map coordinate xhat from reference cell to coordinate x in cell
4647
virtual void map_from_reference_cell(double* x,
4649
const ufc::cell& c) const
4651
throw std::runtime_error("map_from_reference_cell not yet implemented.");
4654
/// Map from coordinate x in cell to coordinate xhat in reference cell
4655
virtual void map_to_reference_cell(double* xhat,
4657
const ufc::cell& c) const
4659
throw std::runtime_error("map_to_reference_cell not yet implemented.");
4662
/// Return the number of sub elements (for a mixed element)
4663
virtual std::size_t num_sub_elements() const
4668
/// Create a new finite element for sub element i (for a mixed element)
4669
virtual ufc::finite_element* create_sub_element(std::size_t i) const
4675
return new navierstokes_finite_element_1();
4680
return new navierstokes_finite_element_1();
4685
return new navierstokes_finite_element_1();
4693
/// Create a new class instance
4694
virtual ufc::finite_element* create() const
4696
return new navierstokes_finite_element_2();
4701
/// This class defines the interface for a local-to-global mapping of
4702
/// degrees of freedom (dofs).
4704
class navierstokes_dofmap_0: public ufc::dofmap
4709
navierstokes_dofmap_0() : ufc::dofmap()
4715
virtual ~navierstokes_dofmap_0()
4720
/// Return a string identifying the dofmap
4721
virtual const char* signature() const
4723
return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 0, None)";
4726
/// Return true iff mesh entities of topological dimension d are needed
4727
virtual bool needs_mesh_entities(std::size_t d) const
4756
/// Return the topological dimension of the associated cell shape
4757
virtual std::size_t topological_dimension() const
4762
/// Return the geometric dimension of the associated cell shape
4763
virtual std::size_t geometric_dimension() const
4768
/// Return the dimension of the global finite element function space
4769
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4770
num_global_entities) const
4772
return num_global_entities[3];
4775
/// Return the dimension of the local finite element function space for a cell
4776
virtual std::size_t local_dimension() const
4781
/// Return the number of dofs on each cell facet
4782
virtual std::size_t num_facet_dofs() const
4787
/// Return the number of dofs associated with each cell entity of dimension d
4788
virtual std::size_t num_entity_dofs(std::size_t d) const
4817
/// Tabulate the local-to-global mapping of dofs on a cell
4818
virtual void tabulate_dofs(std::size_t* dofs,
4819
const std::vector<std::size_t>& num_global_entities,
4820
const ufc::cell& c) const
4822
dofs[0] = c.entity_indices[3][0];
4825
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4826
virtual void tabulate_facet_dofs(std::size_t* dofs,
4827
std::size_t facet) const
4855
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4856
virtual void tabulate_entity_dofs(std::size_t* dofs,
4857
std::size_t d, std::size_t i) const
4861
throw std::runtime_error("d is larger than dimension (3)");
4885
throw std::runtime_error("i is larger than number of entities (0)");
4895
/// Tabulate the coordinates of all dofs on a cell
4896
virtual void tabulate_coordinates(double** dof_coordinates,
4897
const double* vertex_coordinates) const
4899
dof_coordinates[0][0] = 0.25*vertex_coordinates[0] + 0.25*vertex_coordinates[3] + 0.25*vertex_coordinates[6] + 0.25*vertex_coordinates[9];
4900
dof_coordinates[0][1] = 0.25*vertex_coordinates[1] + 0.25*vertex_coordinates[4] + 0.25*vertex_coordinates[7] + 0.25*vertex_coordinates[10];
4901
dof_coordinates[0][2] = 0.25*vertex_coordinates[2] + 0.25*vertex_coordinates[5] + 0.25*vertex_coordinates[8] + 0.25*vertex_coordinates[11];
4904
/// Return the number of sub dofmaps (for a mixed element)
4905
virtual std::size_t num_sub_dofmaps() const
4910
/// Create a new dofmap for sub dofmap i (for a mixed element)
4911
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4916
/// Create a new class instance
4917
virtual ufc::dofmap* create() const
4919
return new navierstokes_dofmap_0();
4924
/// This class defines the interface for a local-to-global mapping of
4925
/// degrees of freedom (dofs).
4927
class navierstokes_dofmap_1: public ufc::dofmap
4932
navierstokes_dofmap_1() : ufc::dofmap()
4938
virtual ~navierstokes_dofmap_1()
4943
/// Return a string identifying the dofmap
4944
virtual const char* signature() const
4946
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, None)";
4949
/// Return true iff mesh entities of topological dimension d are needed
4950
virtual bool needs_mesh_entities(std::size_t d) const
4979
/// Return the topological dimension of the associated cell shape
4980
virtual std::size_t topological_dimension() const
4985
/// Return the geometric dimension of the associated cell shape
4986
virtual std::size_t geometric_dimension() const
4991
/// Return the dimension of the global finite element function space
4992
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4993
num_global_entities) const
4995
return num_global_entities[0];
4998
/// Return the dimension of the local finite element function space for a cell
4999
virtual std::size_t local_dimension() const
5004
/// Return the number of dofs on each cell facet
5005
virtual std::size_t num_facet_dofs() const
5010
/// Return the number of dofs associated with each cell entity of dimension d
5011
virtual std::size_t num_entity_dofs(std::size_t d) const
5040
/// Tabulate the local-to-global mapping of dofs on a cell
5041
virtual void tabulate_dofs(std::size_t* dofs,
5042
const std::vector<std::size_t>& num_global_entities,
5043
const ufc::cell& c) const
5045
dofs[0] = c.entity_indices[0][0];
5046
dofs[1] = c.entity_indices[0][1];
5047
dofs[2] = c.entity_indices[0][2];
5048
dofs[3] = c.entity_indices[0][3];
5051
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
5052
virtual void tabulate_facet_dofs(std::size_t* dofs,
5053
std::size_t facet) const
5089
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
5090
virtual void tabulate_entity_dofs(std::size_t* dofs,
5091
std::size_t d, std::size_t i) const
5095
throw std::runtime_error("d is larger than dimension (3)");
5104
throw std::runtime_error("i is larger than number of entities (3)");
5152
/// Tabulate the coordinates of all dofs on a cell
5153
virtual void tabulate_coordinates(double** dof_coordinates,
5154
const double* vertex_coordinates) const
5156
dof_coordinates[0][0] = vertex_coordinates[0];
5157
dof_coordinates[0][1] = vertex_coordinates[1];
5158
dof_coordinates[0][2] = vertex_coordinates[2];
5159
dof_coordinates[1][0] = vertex_coordinates[3];
5160
dof_coordinates[1][1] = vertex_coordinates[4];
5161
dof_coordinates[1][2] = vertex_coordinates[5];
5162
dof_coordinates[2][0] = vertex_coordinates[6];
5163
dof_coordinates[2][1] = vertex_coordinates[7];
5164
dof_coordinates[2][2] = vertex_coordinates[8];
5165
dof_coordinates[3][0] = vertex_coordinates[9];
5166
dof_coordinates[3][1] = vertex_coordinates[10];
5167
dof_coordinates[3][2] = vertex_coordinates[11];
5170
/// Return the number of sub dofmaps (for a mixed element)
5171
virtual std::size_t num_sub_dofmaps() const
5176
/// Create a new dofmap for sub dofmap i (for a mixed element)
5177
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
5182
/// Create a new class instance
5183
virtual ufc::dofmap* create() const
5185
return new navierstokes_dofmap_1();
5190
/// This class defines the interface for a local-to-global mapping of
5191
/// degrees of freedom (dofs).
5193
class navierstokes_dofmap_2: public ufc::dofmap
5198
navierstokes_dofmap_2() : ufc::dofmap()
5204
virtual ~navierstokes_dofmap_2()
5209
/// Return a string identifying the dofmap
5210
virtual const char* signature() const
5212
return "FFC dofmap for VectorElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 1, 3, None)";
5215
/// Return true iff mesh entities of topological dimension d are needed
5216
virtual bool needs_mesh_entities(std::size_t d) const
5245
/// Return the topological dimension of the associated cell shape
5246
virtual std::size_t topological_dimension() const
5251
/// Return the geometric dimension of the associated cell shape
5252
virtual std::size_t geometric_dimension() const
5257
/// Return the dimension of the global finite element function space
5258
virtual std::size_t global_dimension(const std::vector<std::size_t>&
5259
num_global_entities) const
5261
return 3*num_global_entities[0];
5264
/// Return the dimension of the local finite element function space for a cell
5265
virtual std::size_t local_dimension() const
5270
/// Return the number of dofs on each cell facet
5271
virtual std::size_t num_facet_dofs() const
5276
/// Return the number of dofs associated with each cell entity of dimension d
5277
virtual std::size_t num_entity_dofs(std::size_t d) const
5306
/// Tabulate the local-to-global mapping of dofs on a cell
5307
virtual void tabulate_dofs(std::size_t* dofs,
5308
const std::vector<std::size_t>& num_global_entities,
5309
const ufc::cell& c) const
5311
unsigned int offset = 0;
5312
dofs[0] = offset + c.entity_indices[0][0];
5313
dofs[1] = offset + c.entity_indices[0][1];
5314
dofs[2] = offset + c.entity_indices[0][2];
5315
dofs[3] = offset + c.entity_indices[0][3];
5316
offset += num_global_entities[0];
5317
dofs[4] = offset + c.entity_indices[0][0];
5318
dofs[5] = offset + c.entity_indices[0][1];
5319
dofs[6] = offset + c.entity_indices[0][2];
5320
dofs[7] = offset + c.entity_indices[0][3];
5321
offset += num_global_entities[0];
5322
dofs[8] = offset + c.entity_indices[0][0];
5323
dofs[9] = offset + c.entity_indices[0][1];
5324
dofs[10] = offset + c.entity_indices[0][2];
5325
dofs[11] = offset + c.entity_indices[0][3];
5326
offset += num_global_entities[0];
5329
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
5330
virtual void tabulate_facet_dofs(std::size_t* dofs,
5331
std::size_t facet) const
5391
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
5392
virtual void tabulate_entity_dofs(std::size_t* dofs,
5393
std::size_t d, std::size_t i) const
5397
throw std::runtime_error("d is larger than dimension (3)");
5406
throw std::runtime_error("i is larger than number of entities (3)");
5462
/// Tabulate the coordinates of all dofs on a cell
5463
virtual void tabulate_coordinates(double** dof_coordinates,
5464
const double* vertex_coordinates) const
5466
dof_coordinates[0][0] = vertex_coordinates[0];
5467
dof_coordinates[0][1] = vertex_coordinates[1];
5468
dof_coordinates[0][2] = vertex_coordinates[2];
5469
dof_coordinates[1][0] = vertex_coordinates[3];
5470
dof_coordinates[1][1] = vertex_coordinates[4];
5471
dof_coordinates[1][2] = vertex_coordinates[5];
5472
dof_coordinates[2][0] = vertex_coordinates[6];
5473
dof_coordinates[2][1] = vertex_coordinates[7];
5474
dof_coordinates[2][2] = vertex_coordinates[8];
5475
dof_coordinates[3][0] = vertex_coordinates[9];
5476
dof_coordinates[3][1] = vertex_coordinates[10];
5477
dof_coordinates[3][2] = vertex_coordinates[11];
5478
dof_coordinates[4][0] = vertex_coordinates[0];
5479
dof_coordinates[4][1] = vertex_coordinates[1];
5480
dof_coordinates[4][2] = vertex_coordinates[2];
5481
dof_coordinates[5][0] = vertex_coordinates[3];
5482
dof_coordinates[5][1] = vertex_coordinates[4];
5483
dof_coordinates[5][2] = vertex_coordinates[5];
5484
dof_coordinates[6][0] = vertex_coordinates[6];
5485
dof_coordinates[6][1] = vertex_coordinates[7];
5486
dof_coordinates[6][2] = vertex_coordinates[8];
5487
dof_coordinates[7][0] = vertex_coordinates[9];
5488
dof_coordinates[7][1] = vertex_coordinates[10];
5489
dof_coordinates[7][2] = vertex_coordinates[11];
5490
dof_coordinates[8][0] = vertex_coordinates[0];
5491
dof_coordinates[8][1] = vertex_coordinates[1];
5492
dof_coordinates[8][2] = vertex_coordinates[2];
5493
dof_coordinates[9][0] = vertex_coordinates[3];
5494
dof_coordinates[9][1] = vertex_coordinates[4];
5495
dof_coordinates[9][2] = vertex_coordinates[5];
5496
dof_coordinates[10][0] = vertex_coordinates[6];
5497
dof_coordinates[10][1] = vertex_coordinates[7];
5498
dof_coordinates[10][2] = vertex_coordinates[8];
5499
dof_coordinates[11][0] = vertex_coordinates[9];
5500
dof_coordinates[11][1] = vertex_coordinates[10];
5501
dof_coordinates[11][2] = vertex_coordinates[11];
5504
/// Return the number of sub dofmaps (for a mixed element)
5505
virtual std::size_t num_sub_dofmaps() const
5510
/// Create a new dofmap for sub dofmap i (for a mixed element)
5511
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
5517
return new navierstokes_dofmap_1();
5522
return new navierstokes_dofmap_1();
5527
return new navierstokes_dofmap_1();
5535
/// Create a new class instance
5536
virtual ufc::dofmap* create() const
5538
return new navierstokes_dofmap_2();
5543
/// This class defines the interface for the tabulation of the cell
5544
/// tensor corresponding to the local contribution to a form from
5545
/// the integral over a cell.
5547
class navierstokes_cell_integral_0_otherwise: public ufc::cell_integral
5552
navierstokes_cell_integral_0_otherwise() : ufc::cell_integral()
5558
virtual ~navierstokes_cell_integral_0_otherwise()
5563
/// Tabulate the tensor for the contribution from a local cell
5564
virtual void tabulate_tensor(double* A,
5565
const double * const * w,
5566
const double* vertex_coordinates,
5567
int cell_orientation) const
5571
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
5573
// Compute Jacobian inverse and determinant
5576
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
5579
const double det = std::abs(detJ);
5583
// Compute circumradius
5586
// Facet area (divide by two because 'det' is scaled by area of reference triangle)
5588
// Array of quadrature weights.
5589
static const double W4[4] = {0.0416666666666667, 0.0416666666666667, 0.0416666666666667, 0.0416666666666667};
5590
// Quadrature points on the UFC reference element: (0.585410196624969, 0.138196601125011, 0.138196601125011), (0.138196601125011, 0.585410196624969, 0.138196601125011), (0.138196601125011, 0.138196601125011, 0.585410196624969), (0.138196601125011, 0.138196601125011, 0.138196601125011)
5592
// Value of basis functions at quadrature points.
5593
static const double FE1_C0[4][4] = \
5594
{{0.138196601125009, 0.585410196624969, 0.138196601125011, 0.138196601125011},
5595
{0.138196601125009, 0.138196601125011, 0.585410196624969, 0.138196601125011},
5596
{0.138196601125009, 0.138196601125011, 0.138196601125011, 0.585410196624969},
5597
{0.585410196624967, 0.138196601125011, 0.138196601125011, 0.138196601125011}};
5599
// Array of non-zero columns
5600
static const unsigned int nzc0[4] = {0, 1, 2, 3};
5602
// Array of non-zero columns
5603
static const unsigned int nzc4[4] = {4, 5, 6, 7};
5605
// Array of non-zero columns
5606
static const unsigned int nzc8[4] = {8, 9, 10, 11};
5608
static const double FE1_C0_D001[4][2] = \
5614
// Array of non-zero columns
5615
static const unsigned int nzc1[2] = {0, 3};
5617
// Array of non-zero columns
5618
static const unsigned int nzc3[2] = {0, 1};
5620
// Array of non-zero columns
5621
static const unsigned int nzc9[2] = {8, 11};
5623
// Array of non-zero columns
5624
static const unsigned int nzc2[2] = {0, 2};
5626
// Array of non-zero columns
5627
static const unsigned int nzc6[2] = {4, 6};
5629
// Array of non-zero columns
5630
static const unsigned int nzc7[2] = {4, 5};
5632
// Array of non-zero columns
5633
static const unsigned int nzc11[2] = {8, 9};
5635
// Array of non-zero columns
5636
static const unsigned int nzc10[2] = {8, 10};
5638
// Array of non-zero columns
5639
static const unsigned int nzc5[2] = {4, 7};
5641
// Reset values in the element tensor.
5642
for (unsigned int r = 0; r < 144; r++)
5645
}// end loop over 'r'
5646
// Number of operations to compute geometry constants: 567.
5648
G[0] = 0.5*det*w[3][0]*(K[5]*K[5]*w[2][0] + w[4][0]*(K[3]*K[3] + K[4]*K[4] + K[5]*K[5]));
5649
G[1] = 0.5*K[3]*K[3]*det*w[1][0]*w[3][0];
5650
G[2] = K[3]*K[4]*det*w[1][0]*w[3][0];
5651
G[3] = K[3]*K[5]*det*w[1][0]*w[3][0];
5652
G[4] = 0.5*K[4]*K[4]*det*w[1][0]*w[3][0];
5653
G[5] = K[4]*K[5]*det*w[1][0]*w[3][0];
5654
G[6] = 0.5*K[5]*K[5]*det*w[1][0]*w[3][0];
5655
G[7] = 0.5*det*w[3][0]*(K[2]*K[5]*w[2][0] + w[4][0]*(K[0]*K[3] + K[1]*K[4] + K[2]*K[5]));
5656
G[8] = 0.5*K[0]*K[3]*det*w[1][0]*w[3][0];
5657
G[9] = 0.5*det*w[1][0]*w[3][0]*(K[0]*K[4] + K[1]*K[3]);
5658
G[10] = 0.5*det*w[1][0]*w[3][0]*(K[0]*K[5] + K[2]*K[3]);
5659
G[11] = 0.5*K[1]*K[4]*det*w[1][0]*w[3][0];
5660
G[12] = 0.5*det*w[1][0]*w[3][0]*(K[1]*K[5] + K[2]*K[4]);
5661
G[13] = 0.5*K[2]*K[5]*det*w[1][0]*w[3][0];
5662
G[14] = 0.5*K[5]*K[6]*det*w[2][0]*w[3][0];
5663
G[15] = 0.5*K[3]*K[5]*det*w[2][0]*w[3][0];
5664
G[16] = 0.5*K[0]*K[5]*det*w[2][0]*w[3][0];
5665
G[17] = 0.5*K[5]*K[7]*det*w[2][0]*w[3][0];
5666
G[18] = 0.5*K[4]*K[5]*det*w[2][0]*w[3][0];
5667
G[19] = 0.5*K[1]*K[5]*det*w[2][0]*w[3][0];
5668
G[20] = 0.5*det*w[3][0]*(K[5]*K[8]*w[2][0] + w[4][0]*(K[3]*K[6] + K[4]*K[7] + K[5]*K[8]));
5669
G[21] = 0.5*K[3]*K[6]*det*w[1][0]*w[3][0];
5670
G[22] = 0.5*det*w[1][0]*w[3][0]*(K[3]*K[7] + K[4]*K[6]);
5671
G[23] = 0.5*det*w[1][0]*w[3][0]*(K[3]*K[8] + K[5]*K[6]);
5672
G[24] = 0.5*K[4]*K[7]*det*w[1][0]*w[3][0];
5673
G[25] = 0.5*det*w[1][0]*w[3][0]*(K[4]*K[8] + K[5]*K[7]);
5674
G[26] = 0.5*K[5]*K[8]*det*w[1][0]*w[3][0];
5675
G[27] = 0.5*det*w[3][0]*(K[2]*K[2]*w[2][0] + w[4][0]*(K[0]*K[0] + K[1]*K[1] + K[2]*K[2]));
5676
G[28] = 0.5*K[0]*K[0]*det*w[1][0]*w[3][0];
5677
G[29] = K[0]*K[1]*det*w[1][0]*w[3][0];
5678
G[30] = K[0]*K[2]*det*w[1][0]*w[3][0];
5679
G[31] = 0.5*K[1]*K[1]*det*w[1][0]*w[3][0];
5680
G[32] = K[1]*K[2]*det*w[1][0]*w[3][0];
5681
G[33] = 0.5*K[2]*K[2]*det*w[1][0]*w[3][0];
5682
G[34] = 0.5*K[2]*K[6]*det*w[2][0]*w[3][0];
5683
G[35] = 0.5*K[2]*K[3]*det*w[2][0]*w[3][0];
5684
G[36] = 0.5*K[0]*K[2]*det*w[2][0]*w[3][0];
5685
G[37] = 0.5*K[2]*K[7]*det*w[2][0]*w[3][0];
5686
G[38] = 0.5*K[2]*K[4]*det*w[2][0]*w[3][0];
5687
G[39] = 0.5*K[1]*K[2]*det*w[2][0]*w[3][0];
5688
G[40] = 0.5*det*w[3][0]*(K[2]*K[8]*w[2][0] + w[4][0]*(K[0]*K[6] + K[1]*K[7] + K[2]*K[8]));
5689
G[41] = 0.5*K[0]*K[6]*det*w[1][0]*w[3][0];
5690
G[42] = 0.5*det*w[1][0]*w[3][0]*(K[0]*K[7] + K[1]*K[6]);
5691
G[43] = 0.5*det*w[1][0]*w[3][0]*(K[0]*K[8] + K[2]*K[6]);
5692
G[44] = 0.5*K[1]*K[7]*det*w[1][0]*w[3][0];
5693
G[45] = 0.5*det*w[1][0]*w[3][0]*(K[1]*K[8] + K[2]*K[7]);
5694
G[46] = 0.5*K[2]*K[8]*det*w[1][0]*w[3][0];
5695
G[47] = 0.5*det*w[3][0]*(K[6]*K[6]*w[2][0] + w[4][0]*(K[6]*K[6] + K[7]*K[7] + K[8]*K[8]));
5696
G[48] = 0.5*K[6]*K[6]*det*w[1][0]*w[3][0];
5697
G[49] = K[6]*K[7]*det*w[1][0]*w[3][0];
5698
G[50] = K[6]*K[8]*det*w[1][0]*w[3][0];
5699
G[51] = 0.5*K[7]*K[7]*det*w[1][0]*w[3][0];
5700
G[52] = K[7]*K[8]*det*w[1][0]*w[3][0];
5701
G[53] = 0.5*K[8]*K[8]*det*w[1][0]*w[3][0];
5702
G[54] = 0.5*det*w[3][0]*(K[3]*K[6]*w[2][0] + w[4][0]*(K[3]*K[6] + K[4]*K[7] + K[5]*K[8]));
5703
G[55] = 0.5*det*w[3][0]*(K[0]*K[6]*w[2][0] + w[4][0]*(K[0]*K[6] + K[1]*K[7] + K[2]*K[8]));
5704
G[56] = 0.5*K[6]*K[7]*det*w[2][0]*w[3][0];
5705
G[57] = 0.5*K[4]*K[6]*det*w[2][0]*w[3][0];
5706
G[58] = 0.5*K[1]*K[6]*det*w[2][0]*w[3][0];
5707
G[59] = 0.5*K[6]*K[8]*det*w[2][0]*w[3][0];
5708
G[60] = 0.5*det*w[3][0]*(K[3]*K[3]*w[2][0] + w[4][0]*(K[3]*K[3] + K[4]*K[4] + K[5]*K[5]));
5709
G[61] = 0.5*det*w[3][0]*(K[0]*K[3]*w[2][0] + w[4][0]*(K[0]*K[3] + K[1]*K[4] + K[2]*K[5]));
5710
G[62] = 0.5*K[3]*K[7]*det*w[2][0]*w[3][0];
5711
G[63] = 0.5*K[3]*K[4]*det*w[2][0]*w[3][0];
5712
G[64] = 0.5*K[1]*K[3]*det*w[2][0]*w[3][0];
5713
G[65] = 0.5*K[3]*K[8]*det*w[2][0]*w[3][0];
5714
G[66] = 0.5*det*w[3][0]*(K[0]*K[0]*w[2][0] + w[4][0]*(K[0]*K[0] + K[1]*K[1] + K[2]*K[2]));
5715
G[67] = 0.5*K[0]*K[7]*det*w[2][0]*w[3][0];
5716
G[68] = 0.5*K[0]*K[4]*det*w[2][0]*w[3][0];
5717
G[69] = 0.5*K[0]*K[1]*det*w[2][0]*w[3][0];
5718
G[70] = 0.5*K[0]*K[8]*det*w[2][0]*w[3][0];
5719
G[71] = 0.5*det*w[3][0]*(K[7]*K[7]*w[2][0] + w[4][0]*(K[6]*K[6] + K[7]*K[7] + K[8]*K[8]));
5720
G[72] = 0.5*det*w[3][0]*(K[4]*K[7]*w[2][0] + w[4][0]*(K[3]*K[6] + K[4]*K[7] + K[5]*K[8]));
5721
G[73] = 0.5*det*w[3][0]*(K[1]*K[7]*w[2][0] + w[4][0]*(K[0]*K[6] + K[1]*K[7] + K[2]*K[8]));
5722
G[74] = 0.5*K[7]*K[8]*det*w[2][0]*w[3][0];
5723
G[75] = 0.5*det*w[3][0]*(K[4]*K[4]*w[2][0] + w[4][0]*(K[3]*K[3] + K[4]*K[4] + K[5]*K[5]));
5724
G[76] = 0.5*det*w[3][0]*(K[1]*K[4]*w[2][0] + w[4][0]*(K[0]*K[3] + K[1]*K[4] + K[2]*K[5]));
5725
G[77] = 0.5*K[4]*K[8]*det*w[2][0]*w[3][0];
5726
G[78] = 0.5*det*w[3][0]*(K[1]*K[1]*w[2][0] + w[4][0]*(K[0]*K[0] + K[1]*K[1] + K[2]*K[2]));
5727
G[79] = 0.5*K[1]*K[8]*det*w[2][0]*w[3][0];
5728
G[80] = 0.5*det*w[3][0]*(K[8]*K[8]*w[2][0] + w[4][0]*(K[6]*K[6] + K[7]*K[7] + K[8]*K[8]));
5729
G[81] = 0.5*K[6]*det*w[3][0];
5730
G[82] = 0.5*K[7]*det*w[3][0];
5731
G[83] = 0.5*K[8]*det*w[3][0];
5732
G[84] = 0.5*K[3]*det*w[3][0];
5733
G[85] = 0.5*K[4]*det*w[3][0];
5734
G[86] = 0.5*K[5]*det*w[3][0];
5735
G[87] = 0.5*K[0]*det*w[3][0];
5736
G[88] = 0.5*K[1]*det*w[3][0];
5737
G[89] = 0.5*K[2]*det*w[3][0];
5739
// Compute element tensor using UFL quadrature representation
5740
// Optimisations: ('eliminate zeros', True), ('ignore ones', True), ('ignore zero tables', True), ('optimisation', 'simplify_expressions'), ('remove zero terms', True)
5742
// Loop quadrature points for integral.
5743
// Number of operations to compute element tensor for following IP loop = 6760
5744
for (unsigned int ip = 0; ip < 4; ip++)
5747
// Coefficient declarations.
5752
// Total number of operations to compute function values = 24
5753
for (unsigned int r = 0; r < 4; r++)
5755
F0 += FE1_C0[ip][r]*w[0][nzc0[r]];
5756
F1 += FE1_C0[ip][r]*w[0][nzc4[r]];
5757
F2 += FE1_C0[ip][r]*w[0][nzc8[r]];
5758
}// end loop over 'r'
5760
// Number of operations to compute ip constants: 334
5762
// Number of operations: 16
5763
I[0] = W4[ip]*(G[0] + F0*F0*G[1] + F1*(F0*G[2] + F1*G[4]) + F2*(F0*G[3] + F1*G[5] + F2*G[6]));
5765
// Number of operations: 16
5766
I[1] = W4[ip]*(G[7] + F0*F0*G[8] + F1*(F0*G[9] + F1*G[11]) + F2*(F0*G[10] + F1*G[12] + F2*G[13]));
5768
// Number of operations: 1
5769
I[2] = G[14]*W4[ip];
5771
// Number of operations: 1
5772
I[3] = G[15]*W4[ip];
5774
// Number of operations: 1
5775
I[4] = G[16]*W4[ip];
5777
// Number of operations: 1
5778
I[5] = G[17]*W4[ip];
5780
// Number of operations: 1
5781
I[6] = G[18]*W4[ip];
5783
// Number of operations: 1
5784
I[7] = G[19]*W4[ip];
5786
// Number of operations: 16
5787
I[8] = W4[ip]*(G[20] + F0*F0*G[21] + F1*(F0*G[22] + F1*G[24]) + F2*(F0*G[23] + F1*G[25] + F2*G[26]));
5789
// Number of operations: 16
5790
I[9] = W4[ip]*(G[27] + F0*F0*G[28] + F1*(F0*G[29] + F1*G[31]) + F2*(F0*G[30] + F1*G[32] + F2*G[33]));
5792
// Number of operations: 1
5793
I[10] = G[34]*W4[ip];
5795
// Number of operations: 1
5796
I[11] = G[35]*W4[ip];
5798
// Number of operations: 1
5799
I[12] = G[36]*W4[ip];
5801
// Number of operations: 1
5802
I[13] = G[37]*W4[ip];
5804
// Number of operations: 1
5805
I[14] = G[38]*W4[ip];
5807
// Number of operations: 1
5808
I[15] = G[39]*W4[ip];
5810
// Number of operations: 16
5811
I[16] = W4[ip]*(G[40] + F0*F0*G[41] + F1*(F0*G[42] + F1*G[44]) + F2*(F0*G[43] + F1*G[45] + F2*G[46]));
5813
// Number of operations: 16
5814
I[17] = W4[ip]*(G[47] + F0*F0*G[48] + F1*(F0*G[49] + F1*G[51]) + F2*(F0*G[50] + F1*G[52] + F2*G[53]));
5816
// Number of operations: 16
5817
I[18] = W4[ip]*(G[54] + F0*F0*G[21] + F1*(F0*G[22] + F1*G[24]) + F2*(F0*G[23] + F1*G[25] + F2*G[26]));
5819
// Number of operations: 16
5820
I[19] = W4[ip]*(G[55] + F0*F0*G[41] + F1*(F0*G[42] + F1*G[44]) + F2*(F0*G[43] + F1*G[45] + F2*G[46]));
5822
// Number of operations: 1
5823
I[20] = G[56]*W4[ip];
5825
// Number of operations: 1
5826
I[21] = G[57]*W4[ip];
5828
// Number of operations: 1
5829
I[22] = G[58]*W4[ip];
5831
// Number of operations: 1
5832
I[23] = G[59]*W4[ip];
5834
// Number of operations: 16
5835
I[24] = W4[ip]*(G[60] + F0*F0*G[1] + F1*(F0*G[2] + F1*G[4]) + F2*(F0*G[3] + F1*G[5] + F2*G[6]));
5837
// Number of operations: 16
5838
I[25] = W4[ip]*(G[61] + F0*F0*G[8] + F1*(F0*G[9] + F1*G[11]) + F2*(F0*G[10] + F1*G[12] + F2*G[13]));
5840
// Number of operations: 1
5841
I[26] = G[62]*W4[ip];
5843
// Number of operations: 1
5844
I[27] = G[63]*W4[ip];
5846
// Number of operations: 1
5847
I[28] = G[64]*W4[ip];
5849
// Number of operations: 1
5850
I[29] = G[65]*W4[ip];
5852
// Number of operations: 16
5853
I[30] = W4[ip]*(G[66] + F0*F0*G[28] + F1*(F0*G[29] + F1*G[31]) + F2*(F0*G[30] + F1*G[32] + F2*G[33]));
5855
// Number of operations: 1
5856
I[31] = G[67]*W4[ip];
5858
// Number of operations: 1
5859
I[32] = G[68]*W4[ip];
5861
// Number of operations: 1
5862
I[33] = G[69]*W4[ip];
5864
// Number of operations: 1
5865
I[34] = G[70]*W4[ip];
5867
// Number of operations: 16
5868
I[35] = W4[ip]*(G[71] + F0*F0*G[48] + F1*(F0*G[49] + F1*G[51]) + F2*(F0*G[50] + F1*G[52] + F2*G[53]));
5870
// Number of operations: 16
5871
I[36] = W4[ip]*(G[72] + F0*F0*G[21] + F1*(F0*G[22] + F1*G[24]) + F2*(F0*G[23] + F1*G[25] + F2*G[26]));
5873
// Number of operations: 16
5874
I[37] = W4[ip]*(G[73] + F0*F0*G[41] + F1*(F0*G[42] + F1*G[44]) + F2*(F0*G[43] + F1*G[45] + F2*G[46]));
5876
// Number of operations: 1
5877
I[38] = G[74]*W4[ip];
5879
// Number of operations: 16
5880
I[39] = W4[ip]*(G[75] + F0*F0*G[1] + F1*(F0*G[2] + F1*G[4]) + F2*(F0*G[3] + F1*G[5] + F2*G[6]));
5882
// Number of operations: 16
5883
I[40] = W4[ip]*(G[76] + F0*F0*G[8] + F1*(F0*G[9] + F1*G[11]) + F2*(F0*G[10] + F1*G[12] + F2*G[13]));
5885
// Number of operations: 1
5886
I[41] = G[77]*W4[ip];
5888
// Number of operations: 16
5889
I[42] = W4[ip]*(G[78] + F0*F0*G[28] + F1*(F0*G[29] + F1*G[31]) + F2*(F0*G[30] + F1*G[32] + F2*G[33]));
5891
// Number of operations: 1
5892
I[43] = G[79]*W4[ip];
5894
// Number of operations: 16
5895
I[44] = W4[ip]*(G[80] + F0*F0*G[48] + F1*(F0*G[49] + F1*G[51]) + F2*(F0*G[50] + F1*G[52] + F2*G[53]));
5897
// Number of operations: 6
5898
I[45] = W4[ip]*(F0*G[81] + F1*G[82] + F2*G[83]);
5900
// Number of operations: 6
5901
I[46] = W4[ip]*(F0*G[84] + F1*G[85] + F2*G[86]);
5903
// Number of operations: 6
5904
I[47] = W4[ip]*(F0*G[87] + F1*G[88] + F2*G[89]);
5906
// Number of operations: 1
5910
// Number of operations for primary indices: 216
5911
for (unsigned int j = 0; j < 4; j++)
5913
for (unsigned int k = 0; k < 2; k++)
5915
// Number of operations to compute entry: 3
5916
A[nzc0[j]*12 + nzc1[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[45];
5917
// Number of operations to compute entry: 3
5918
A[nzc0[j]*12 + nzc2[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[46];
5919
// Number of operations to compute entry: 3
5920
A[nzc0[j]*12 + nzc3[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[47];
5921
// Number of operations to compute entry: 3
5922
A[nzc4[j]*12 + nzc5[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[45];
5923
// Number of operations to compute entry: 3
5924
A[nzc4[j]*12 + nzc6[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[46];
5925
// Number of operations to compute entry: 3
5926
A[nzc4[j]*12 + nzc7[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[47];
5927
// Number of operations to compute entry: 3
5928
A[nzc8[j]*12 + nzc10[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[46];
5929
// Number of operations to compute entry: 3
5930
A[nzc8[j]*12 + nzc11[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[47];
5931
// Number of operations to compute entry: 3
5932
A[nzc8[j]*12 + nzc9[k]] += FE1_C0[ip][j]*FE1_C0_D001[ip][k]*I[45];
5933
}// end loop over 'k'
5934
}// end loop over 'j'
5936
// Number of operations for primary indices: 972
5937
for (unsigned int j = 0; j < 2; j++)
5939
for (unsigned int k = 0; k < 2; k++)
5941
// Number of operations to compute entry: 3
5942
A[nzc10[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[0];
5943
// Number of operations to compute entry: 3
5944
A[nzc10[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[1];
5945
// Number of operations to compute entry: 3
5946
A[nzc10[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[2];
5947
// Number of operations to compute entry: 3
5948
A[nzc10[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[3];
5949
// Number of operations to compute entry: 3
5950
A[nzc10[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[4];
5951
// Number of operations to compute entry: 3
5952
A[nzc10[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[5];
5953
// Number of operations to compute entry: 3
5954
A[nzc10[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[6];
5955
// Number of operations to compute entry: 3
5956
A[nzc10[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[7];
5957
// Number of operations to compute entry: 3
5958
A[nzc10[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[8];
5959
// Number of operations to compute entry: 3
5960
A[nzc11[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[1];
5961
// Number of operations to compute entry: 3
5962
A[nzc11[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[9];
5963
// Number of operations to compute entry: 3
5964
A[nzc11[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[10];
5965
// Number of operations to compute entry: 3
5966
A[nzc11[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[11];
5967
// Number of operations to compute entry: 3
5968
A[nzc11[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[12];
5969
// Number of operations to compute entry: 3
5970
A[nzc11[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[13];
5971
// Number of operations to compute entry: 3
5972
A[nzc11[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[14];
5973
// Number of operations to compute entry: 3
5974
A[nzc11[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[15];
5975
// Number of operations to compute entry: 3
5976
A[nzc11[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[16];
5977
// Number of operations to compute entry: 3
5978
A[nzc1[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[2];
5979
// Number of operations to compute entry: 3
5980
A[nzc1[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[10];
5981
// Number of operations to compute entry: 3
5982
A[nzc1[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[17];
5983
// Number of operations to compute entry: 3
5984
A[nzc1[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[18];
5985
// Number of operations to compute entry: 3
5986
A[nzc1[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[19];
5987
// Number of operations to compute entry: 3
5988
A[nzc1[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[20];
5989
// Number of operations to compute entry: 3
5990
A[nzc1[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[21];
5991
// Number of operations to compute entry: 3
5992
A[nzc1[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[22];
5993
// Number of operations to compute entry: 3
5994
A[nzc1[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[23];
5995
// Number of operations to compute entry: 3
5996
A[nzc2[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[3];
5997
// Number of operations to compute entry: 3
5998
A[nzc2[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[11];
5999
// Number of operations to compute entry: 3
6000
A[nzc2[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[18];
6001
// Number of operations to compute entry: 3
6002
A[nzc2[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[24];
6003
// Number of operations to compute entry: 3
6004
A[nzc2[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[25];
6005
// Number of operations to compute entry: 3
6006
A[nzc2[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[26];
6007
// Number of operations to compute entry: 3
6008
A[nzc2[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[27];
6009
// Number of operations to compute entry: 3
6010
A[nzc2[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[28];
6011
// Number of operations to compute entry: 3
6012
A[nzc2[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[29];
6013
// Number of operations to compute entry: 3
6014
A[nzc3[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[4];
6015
// Number of operations to compute entry: 3
6016
A[nzc3[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[12];
6017
// Number of operations to compute entry: 3
6018
A[nzc3[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[19];
6019
// Number of operations to compute entry: 3
6020
A[nzc3[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[25];
6021
// Number of operations to compute entry: 3
6022
A[nzc3[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[30];
6023
// Number of operations to compute entry: 3
6024
A[nzc3[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[31];
6025
// Number of operations to compute entry: 3
6026
A[nzc3[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[32];
6027
// Number of operations to compute entry: 3
6028
A[nzc3[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[33];
6029
// Number of operations to compute entry: 3
6030
A[nzc3[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[34];
6031
// Number of operations to compute entry: 3
6032
A[nzc5[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[5];
6033
// Number of operations to compute entry: 3
6034
A[nzc5[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[13];
6035
// Number of operations to compute entry: 3
6036
A[nzc5[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[20];
6037
// Number of operations to compute entry: 3
6038
A[nzc5[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[26];
6039
// Number of operations to compute entry: 3
6040
A[nzc5[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[31];
6041
// Number of operations to compute entry: 3
6042
A[nzc5[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[35];
6043
// Number of operations to compute entry: 3
6044
A[nzc5[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[36];
6045
// Number of operations to compute entry: 3
6046
A[nzc5[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[37];
6047
// Number of operations to compute entry: 3
6048
A[nzc5[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[38];
6049
// Number of operations to compute entry: 3
6050
A[nzc6[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[6];
6051
// Number of operations to compute entry: 3
6052
A[nzc6[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[14];
6053
// Number of operations to compute entry: 3
6054
A[nzc6[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[21];
6055
// Number of operations to compute entry: 3
6056
A[nzc6[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[27];
6057
// Number of operations to compute entry: 3
6058
A[nzc6[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[32];
6059
// Number of operations to compute entry: 3
6060
A[nzc6[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[36];
6061
// Number of operations to compute entry: 3
6062
A[nzc6[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[39];
6063
// Number of operations to compute entry: 3
6064
A[nzc6[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[40];
6065
// Number of operations to compute entry: 3
6066
A[nzc6[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[41];
6067
// Number of operations to compute entry: 3
6068
A[nzc7[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[7];
6069
// Number of operations to compute entry: 3
6070
A[nzc7[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[15];
6071
// Number of operations to compute entry: 3
6072
A[nzc7[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[22];
6073
// Number of operations to compute entry: 3
6074
A[nzc7[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[28];
6075
// Number of operations to compute entry: 3
6076
A[nzc7[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[33];
6077
// Number of operations to compute entry: 3
6078
A[nzc7[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[37];
6079
// Number of operations to compute entry: 3
6080
A[nzc7[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[40];
6081
// Number of operations to compute entry: 3
6082
A[nzc7[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[42];
6083
// Number of operations to compute entry: 3
6084
A[nzc7[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[43];
6085
// Number of operations to compute entry: 3
6086
A[nzc9[j]*12 + nzc10[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[8];
6087
// Number of operations to compute entry: 3
6088
A[nzc9[j]*12 + nzc11[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[16];
6089
// Number of operations to compute entry: 3
6090
A[nzc9[j]*12 + nzc1[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[23];
6091
// Number of operations to compute entry: 3
6092
A[nzc9[j]*12 + nzc2[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[29];
6093
// Number of operations to compute entry: 3
6094
A[nzc9[j]*12 + nzc3[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[34];
6095
// Number of operations to compute entry: 3
6096
A[nzc9[j]*12 + nzc5[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[38];
6097
// Number of operations to compute entry: 3
6098
A[nzc9[j]*12 + nzc6[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[41];
6099
// Number of operations to compute entry: 3
6100
A[nzc9[j]*12 + nzc7[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[43];
6101
// Number of operations to compute entry: 3
6102
A[nzc9[j]*12 + nzc9[k]] += FE1_C0_D001[ip][j]*FE1_C0_D001[ip][k]*I[44];
6103
}// end loop over 'k'
6104
}// end loop over 'j'
6106
// Number of operations for primary indices: 144
6107
for (unsigned int j = 0; j < 4; j++)
6109
for (unsigned int k = 0; k < 4; k++)
6111
// Number of operations to compute entry: 3
6112
A[nzc0[j]*12 + nzc0[k]] += FE1_C0[ip][j]*FE1_C0[ip][k]*I[48];
6113
// Number of operations to compute entry: 3
6114
A[nzc4[j]*12 + nzc4[k]] += FE1_C0[ip][j]*FE1_C0[ip][k]*I[48];
6115
// Number of operations to compute entry: 3
6116
A[nzc8[j]*12 + nzc8[k]] += FE1_C0[ip][j]*FE1_C0[ip][k]*I[48];
6117
}// end loop over 'k'
6118
}// end loop over 'j'
6119
}// end loop over 'ip'
6124
/// This class defines the interface for the assembly of the global
6125
/// tensor corresponding to a form with r + n arguments, that is, a
6128
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
6130
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
6131
/// global tensor A is defined by
6133
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
6135
/// where each argument Vj represents the application to the
6136
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
6137
/// fixed functions (coefficients).
6139
class navierstokes_form_0: public ufc::form
6144
navierstokes_form_0() : ufc::form()
6150
virtual ~navierstokes_form_0()
6155
/// Return a string identifying the form
6156
virtual const char* signature() const
6158
return "119156fa8b43c2d4219ae3a35407b2cc76c4df05cd875a026e131301f5136d4bdf5b6cdc72fb64902c0d7de2096e2c7d6a4ace396a7d0bdbe95636d07b8bf6b7";
6161
/// Return the rank of the global tensor (r)
6162
virtual std::size_t rank() const
6167
/// Return the number of coefficients (n)
6168
virtual std::size_t num_coefficients() const
6173
/// Return the number of cell domains
6174
virtual std::size_t num_cell_domains() const
6179
/// Return the number of exterior facet domains
6180
virtual std::size_t num_exterior_facet_domains() const
6185
/// Return the number of interior facet domains
6186
virtual std::size_t num_interior_facet_domains() const
6191
/// Return the number of point domains
6192
virtual std::size_t num_point_domains() const
6197
/// Return whether the form has any cell integrals
6198
virtual bool has_cell_integrals() const
6203
/// Return whether the form has any exterior facet integrals
6204
virtual bool has_exterior_facet_integrals() const
6209
/// Return whether the form has any interior facet integrals
6210
virtual bool has_interior_facet_integrals() const
6215
/// Return whether the form has any point integrals
6216
virtual bool has_point_integrals() const
6221
/// Create a new finite element for argument function i
6222
virtual ufc::finite_element* create_finite_element(std::size_t i) const
6228
return new navierstokes_finite_element_2();
6233
return new navierstokes_finite_element_2();
6238
return new navierstokes_finite_element_2();
6243
return new navierstokes_finite_element_0();
6248
return new navierstokes_finite_element_0();
6253
return new navierstokes_finite_element_0();
6258
return new navierstokes_finite_element_0();
6266
/// Create a new dofmap for argument function i
6267
virtual ufc::dofmap* create_dofmap(std::size_t i) const
6273
return new navierstokes_dofmap_2();
6278
return new navierstokes_dofmap_2();
6283
return new navierstokes_dofmap_2();
6288
return new navierstokes_dofmap_0();
6293
return new navierstokes_dofmap_0();
6298
return new navierstokes_dofmap_0();
6303
return new navierstokes_dofmap_0();
6311
/// Create a new cell integral on sub domain i
6312
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
6317
/// Create a new exterior facet integral on sub domain i
6318
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
6323
/// Create a new interior facet integral on sub domain i
6324
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
6329
/// Create a new point integral on sub domain i
6330
virtual ufc::point_integral* create_point_integral(std::size_t i) const
6335
/// Create a new cell integral on everywhere else
6336
virtual ufc::cell_integral* create_default_cell_integral() const
6338
return new navierstokes_cell_integral_0_otherwise();
6341
/// Create a new exterior facet integral on everywhere else
6342
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
6347
/// Create a new interior facet integral on everywhere else
6348
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
6353
/// Create a new point integral on everywhere else
6354
virtual ufc::point_integral* create_default_point_integral() const
6363
// Standard library includes
6367
#include <dolfin/common/NoDeleter.h>
6368
#include <dolfin/mesh/Restriction.h>
6369
#include <dolfin/fem/FiniteElement.h>
6370
#include <dolfin/fem/DofMap.h>
6371
#include <dolfin/fem/Form.h>
6372
#include <dolfin/function/FunctionSpace.h>
6373
#include <dolfin/function/GenericFunction.h>
6374
#include <dolfin/function/CoefficientAssigner.h>
6375
#include <dolfin/adaptivity/ErrorControl.h>
6376
#include <dolfin/adaptivity/GoalFunctional.h>
6378
namespace NavierStokes
6381
class CoefficientSpace_d1: public dolfin::FunctionSpace
6385
//--- Constructors for standard function space, 2 different versions ---
6387
// Create standard function space (reference version)
6388
CoefficientSpace_d1(const dolfin::Mesh& mesh):
6389
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6390
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6391
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh)))
6396
// Create standard function space (shared pointer version)
6397
CoefficientSpace_d1(boost::shared_ptr<const dolfin::Mesh> mesh):
6398
dolfin::FunctionSpace(mesh,
6399
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6400
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh)))
6405
//--- Constructors for constrained function space, 2 different versions ---
6407
// Create standard function space (reference version)
6408
CoefficientSpace_d1(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
6409
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6410
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6411
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh,
6412
dolfin::reference_to_no_delete_pointer(constrained_domain))))
6417
// Create standard function space (shared pointer version)
6418
CoefficientSpace_d1(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
6419
dolfin::FunctionSpace(mesh,
6420
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6421
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh, constrained_domain)))
6426
//--- Constructors for restricted function space, 2 different versions ---
6428
// Create restricted function space (reference version)
6429
CoefficientSpace_d1(const dolfin::Restriction& restriction):
6430
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
6431
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6432
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6433
reference_to_no_delete_pointer(restriction))))
6438
// Create restricted function space (shared pointer version)
6439
CoefficientSpace_d1(boost::shared_ptr<const dolfin::Restriction> restriction):
6440
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
6441
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6442
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6449
~CoefficientSpace_d1()
6455
class CoefficientSpace_d2: public dolfin::FunctionSpace
6459
//--- Constructors for standard function space, 2 different versions ---
6461
// Create standard function space (reference version)
6462
CoefficientSpace_d2(const dolfin::Mesh& mesh):
6463
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6464
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6465
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh)))
6470
// Create standard function space (shared pointer version)
6471
CoefficientSpace_d2(boost::shared_ptr<const dolfin::Mesh> mesh):
6472
dolfin::FunctionSpace(mesh,
6473
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6474
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh)))
6479
//--- Constructors for constrained function space, 2 different versions ---
6481
// Create standard function space (reference version)
6482
CoefficientSpace_d2(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
6483
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6484
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6485
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh,
6486
dolfin::reference_to_no_delete_pointer(constrained_domain))))
6491
// Create standard function space (shared pointer version)
6492
CoefficientSpace_d2(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
6493
dolfin::FunctionSpace(mesh,
6494
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6495
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh, constrained_domain)))
6500
//--- Constructors for restricted function space, 2 different versions ---
6502
// Create restricted function space (reference version)
6503
CoefficientSpace_d2(const dolfin::Restriction& restriction):
6504
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
6505
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6506
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6507
reference_to_no_delete_pointer(restriction))))
6512
// Create restricted function space (shared pointer version)
6513
CoefficientSpace_d2(boost::shared_ptr<const dolfin::Restriction> restriction):
6514
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
6515
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6516
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6523
~CoefficientSpace_d2()
6529
class CoefficientSpace_k: public dolfin::FunctionSpace
6533
//--- Constructors for standard function space, 2 different versions ---
6535
// Create standard function space (reference version)
6536
CoefficientSpace_k(const dolfin::Mesh& mesh):
6537
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6538
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6539
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh)))
6544
// Create standard function space (shared pointer version)
6545
CoefficientSpace_k(boost::shared_ptr<const dolfin::Mesh> mesh):
6546
dolfin::FunctionSpace(mesh,
6547
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6548
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh)))
6553
//--- Constructors for constrained function space, 2 different versions ---
6555
// Create standard function space (reference version)
6556
CoefficientSpace_k(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
6557
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6558
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6559
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh,
6560
dolfin::reference_to_no_delete_pointer(constrained_domain))))
6565
// Create standard function space (shared pointer version)
6566
CoefficientSpace_k(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
6567
dolfin::FunctionSpace(mesh,
6568
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6569
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh, constrained_domain)))
6574
//--- Constructors for restricted function space, 2 different versions ---
6576
// Create restricted function space (reference version)
6577
CoefficientSpace_k(const dolfin::Restriction& restriction):
6578
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
6579
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6580
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6581
reference_to_no_delete_pointer(restriction))))
6586
// Create restricted function space (shared pointer version)
6587
CoefficientSpace_k(boost::shared_ptr<const dolfin::Restriction> restriction):
6588
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
6589
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6590
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6597
~CoefficientSpace_k()
6603
class CoefficientSpace_nu: public dolfin::FunctionSpace
6607
//--- Constructors for standard function space, 2 different versions ---
6609
// Create standard function space (reference version)
6610
CoefficientSpace_nu(const dolfin::Mesh& mesh):
6611
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6612
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6613
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh)))
6618
// Create standard function space (shared pointer version)
6619
CoefficientSpace_nu(boost::shared_ptr<const dolfin::Mesh> mesh):
6620
dolfin::FunctionSpace(mesh,
6621
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6622
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh)))
6627
//--- Constructors for constrained function space, 2 different versions ---
6629
// Create standard function space (reference version)
6630
CoefficientSpace_nu(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
6631
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6632
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6633
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), mesh,
6634
dolfin::reference_to_no_delete_pointer(constrained_domain))))
6639
// Create standard function space (shared pointer version)
6640
CoefficientSpace_nu(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
6641
dolfin::FunctionSpace(mesh,
6642
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6643
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()), *mesh, constrained_domain)))
6648
//--- Constructors for restricted function space, 2 different versions ---
6650
// Create restricted function space (reference version)
6651
CoefficientSpace_nu(const dolfin::Restriction& restriction):
6652
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
6653
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6654
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6655
reference_to_no_delete_pointer(restriction))))
6660
// Create restricted function space (shared pointer version)
6661
CoefficientSpace_nu(boost::shared_ptr<const dolfin::Restriction> restriction):
6662
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
6663
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_0()))),
6664
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_0()),
6671
~CoefficientSpace_nu()
6677
class CoefficientSpace_w: public dolfin::FunctionSpace
6681
//--- Constructors for standard function space, 2 different versions ---
6683
// Create standard function space (reference version)
6684
CoefficientSpace_w(const dolfin::Mesh& mesh):
6685
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6686
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6687
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), mesh)))
6692
// Create standard function space (shared pointer version)
6693
CoefficientSpace_w(boost::shared_ptr<const dolfin::Mesh> mesh):
6694
dolfin::FunctionSpace(mesh,
6695
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6696
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), *mesh)))
6701
//--- Constructors for constrained function space, 2 different versions ---
6703
// Create standard function space (reference version)
6704
CoefficientSpace_w(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
6705
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6706
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6707
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), mesh,
6708
dolfin::reference_to_no_delete_pointer(constrained_domain))))
6713
// Create standard function space (shared pointer version)
6714
CoefficientSpace_w(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
6715
dolfin::FunctionSpace(mesh,
6716
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6717
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), *mesh, constrained_domain)))
6722
//--- Constructors for restricted function space, 2 different versions ---
6724
// Create restricted function space (reference version)
6725
CoefficientSpace_w(const dolfin::Restriction& restriction):
6726
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
6727
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6728
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()),
6729
reference_to_no_delete_pointer(restriction))))
6734
// Create restricted function space (shared pointer version)
6735
CoefficientSpace_w(boost::shared_ptr<const dolfin::Restriction> restriction):
6736
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
6737
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6738
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()),
6745
~CoefficientSpace_w()
6751
class Form_a_FunctionSpace_0: public dolfin::FunctionSpace
6755
//--- Constructors for standard function space, 2 different versions ---
6757
// Create standard function space (reference version)
6758
Form_a_FunctionSpace_0(const dolfin::Mesh& mesh):
6759
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6760
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6761
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), mesh)))
6766
// Create standard function space (shared pointer version)
6767
Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
6768
dolfin::FunctionSpace(mesh,
6769
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6770
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), *mesh)))
6775
//--- Constructors for constrained function space, 2 different versions ---
6777
// Create standard function space (reference version)
6778
Form_a_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
6779
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6780
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6781
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), mesh,
6782
dolfin::reference_to_no_delete_pointer(constrained_domain))))
6787
// Create standard function space (shared pointer version)
6788
Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
6789
dolfin::FunctionSpace(mesh,
6790
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6791
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), *mesh, constrained_domain)))
6796
//--- Constructors for restricted function space, 2 different versions ---
6798
// Create restricted function space (reference version)
6799
Form_a_FunctionSpace_0(const dolfin::Restriction& restriction):
6800
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
6801
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6802
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()),
6803
reference_to_no_delete_pointer(restriction))))
6808
// Create restricted function space (shared pointer version)
6809
Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Restriction> restriction):
6810
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
6811
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6812
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()),
6819
~Form_a_FunctionSpace_0()
6825
class Form_a_FunctionSpace_1: public dolfin::FunctionSpace
6829
//--- Constructors for standard function space, 2 different versions ---
6831
// Create standard function space (reference version)
6832
Form_a_FunctionSpace_1(const dolfin::Mesh& mesh):
6833
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6834
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6835
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), mesh)))
6840
// Create standard function space (shared pointer version)
6841
Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh):
6842
dolfin::FunctionSpace(mesh,
6843
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6844
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), *mesh)))
6849
//--- Constructors for constrained function space, 2 different versions ---
6851
// Create standard function space (reference version)
6852
Form_a_FunctionSpace_1(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
6853
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
6854
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6855
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), mesh,
6856
dolfin::reference_to_no_delete_pointer(constrained_domain))))
6861
// Create standard function space (shared pointer version)
6862
Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
6863
dolfin::FunctionSpace(mesh,
6864
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6865
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()), *mesh, constrained_domain)))
6870
//--- Constructors for restricted function space, 2 different versions ---
6872
// Create restricted function space (reference version)
6873
Form_a_FunctionSpace_1(const dolfin::Restriction& restriction):
6874
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
6875
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6876
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()),
6877
reference_to_no_delete_pointer(restriction))))
6882
// Create restricted function space (shared pointer version)
6883
Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Restriction> restriction):
6884
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
6885
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new navierstokes_finite_element_2()))),
6886
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new navierstokes_dofmap_2()),
6893
~Form_a_FunctionSpace_1()
6899
typedef CoefficientSpace_w Form_a_FunctionSpace_2;
6901
typedef CoefficientSpace_d1 Form_a_FunctionSpace_3;
6903
typedef CoefficientSpace_d2 Form_a_FunctionSpace_4;
6905
typedef CoefficientSpace_k Form_a_FunctionSpace_5;
6907
typedef CoefficientSpace_nu Form_a_FunctionSpace_6;
6909
class Form_a: public dolfin::Form
6914
Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0):
6915
dolfin::Form(2, 5), w(*this, 0), d1(*this, 1), d2(*this, 2), k(*this, 3), nu(*this, 4)
6917
_function_spaces[0] = reference_to_no_delete_pointer(V0);
6918
_function_spaces[1] = reference_to_no_delete_pointer(V1);
6920
_ufc_form = boost::shared_ptr<const ufc::form>(new navierstokes_form_0());
6924
Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0, const dolfin::GenericFunction& w, const dolfin::GenericFunction& d1, const dolfin::GenericFunction& d2, const dolfin::GenericFunction& k, const dolfin::GenericFunction& nu):
6925
dolfin::Form(2, 5), w(*this, 0), d1(*this, 1), d2(*this, 2), k(*this, 3), nu(*this, 4)
6927
_function_spaces[0] = reference_to_no_delete_pointer(V0);
6928
_function_spaces[1] = reference_to_no_delete_pointer(V1);
6936
_ufc_form = boost::shared_ptr<const ufc::form>(new navierstokes_form_0());
6940
Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0, boost::shared_ptr<const dolfin::GenericFunction> w, boost::shared_ptr<const dolfin::GenericFunction> d1, boost::shared_ptr<const dolfin::GenericFunction> d2, boost::shared_ptr<const dolfin::GenericFunction> k, boost::shared_ptr<const dolfin::GenericFunction> nu):
6941
dolfin::Form(2, 5), w(*this, 0), d1(*this, 1), d2(*this, 2), k(*this, 3), nu(*this, 4)
6943
_function_spaces[0] = reference_to_no_delete_pointer(V0);
6944
_function_spaces[1] = reference_to_no_delete_pointer(V1);
6952
_ufc_form = boost::shared_ptr<const ufc::form>(new navierstokes_form_0());
6956
Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0):
6957
dolfin::Form(2, 5), w(*this, 0), d1(*this, 1), d2(*this, 2), k(*this, 3), nu(*this, 4)
6959
_function_spaces[0] = V0;
6960
_function_spaces[1] = V1;
6962
_ufc_form = boost::shared_ptr<const ufc::form>(new navierstokes_form_0());
6966
Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0, const dolfin::GenericFunction& w, const dolfin::GenericFunction& d1, const dolfin::GenericFunction& d2, const dolfin::GenericFunction& k, const dolfin::GenericFunction& nu):
6967
dolfin::Form(2, 5), w(*this, 0), d1(*this, 1), d2(*this, 2), k(*this, 3), nu(*this, 4)
6969
_function_spaces[0] = V0;
6970
_function_spaces[1] = V1;
6978
_ufc_form = boost::shared_ptr<const ufc::form>(new navierstokes_form_0());
6982
Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::GenericFunction> w, boost::shared_ptr<const dolfin::GenericFunction> d1, boost::shared_ptr<const dolfin::GenericFunction> d2, boost::shared_ptr<const dolfin::GenericFunction> k, boost::shared_ptr<const dolfin::GenericFunction> nu):
6983
dolfin::Form(2, 5), w(*this, 0), d1(*this, 1), d2(*this, 2), k(*this, 3), nu(*this, 4)
6985
_function_spaces[0] = V0;
6986
_function_spaces[1] = V1;
6994
_ufc_form = boost::shared_ptr<const ufc::form>(new navierstokes_form_0());
7001
/// Return the number of the coefficient with this name
7002
virtual std::size_t coefficient_number(const std::string& name) const
7006
else if (name == "d1")
7008
else if (name == "d2")
7010
else if (name == "k")
7012
else if (name == "nu")
7015
dolfin::dolfin_error("generated code for class Form",
7016
"access coefficient data",
7017
"Invalid coefficient");
7021
/// Return the name of the coefficient with this number
7022
virtual std::string coefficient_name(std::size_t i) const
7038
dolfin::dolfin_error("generated code for class Form",
7039
"access coefficient data",
7040
"Invalid coefficient");
7045
typedef Form_a_FunctionSpace_0 TestSpace;
7046
typedef Form_a_FunctionSpace_1 TrialSpace;
7047
typedef Form_a_FunctionSpace_2 CoefficientSpace_w;
7048
typedef Form_a_FunctionSpace_3 CoefficientSpace_d1;
7049
typedef Form_a_FunctionSpace_4 CoefficientSpace_d2;
7050
typedef Form_a_FunctionSpace_5 CoefficientSpace_k;
7051
typedef Form_a_FunctionSpace_6 CoefficientSpace_nu;
7054
dolfin::CoefficientAssigner w;
7055
dolfin::CoefficientAssigner d1;
7056
dolfin::CoefficientAssigner d2;
7057
dolfin::CoefficientAssigner k;
7058
dolfin::CoefficientAssigner nu;
7062
typedef Form_a BilinearForm;
7063
typedef Form_a JacobianForm;
7064
typedef Form_a::TestSpace FunctionSpace;