1
// This code conforms with the UFC specification version 2.2.0
2
// and was automatically generated by FFC version 1.2.0.
4
// This code was generated with the following parameters:
7
// convert_exceptions_to_warnings: True
9
// cpp_optimize_flags: '-O2'
11
// error_control: False
19
// quadrature_degree: 'auto'
20
// quadrature_rule: 'auto'
21
// representation: 'quadrature'
23
// swig_binary: 'swig'
26
#ifndef __X_ELEMENT5_H
27
#define __X_ELEMENT5_H
34
/// This class defines the interface for a finite element.
36
class x_element5_finite_element_0: public ufc::finite_element
41
x_element5_finite_element_0() : ufc::finite_element()
47
virtual ~x_element5_finite_element_0()
52
/// Return a string identifying the finite element
53
virtual const char* signature() const
55
return "FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None)";
58
/// Return the cell shape
59
virtual ufc::shape cell_shape() const
64
/// Return the topological dimension of the cell shape
65
virtual std::size_t topological_dimension() const
70
/// Return the geometric dimension of the cell shape
71
virtual std::size_t geometric_dimension() const
76
/// Return the dimension of the finite element function space
77
virtual std::size_t space_dimension() const
82
/// Return the rank of the value space
83
virtual std::size_t value_rank() const
88
/// Return the dimension of the value space for axis i
89
virtual std::size_t value_dimension(std::size_t i) const
94
/// Evaluate basis function i at given point x in cell
95
virtual void evaluate_basis(std::size_t i,
98
const double* vertex_coordinates,
99
int cell_orientation) const
103
compute_jacobian_triangle_3d(J, vertex_coordinates);
105
// Compute Jacobian inverse and determinant
108
compute_jacobian_inverse_triangle_3d(K, detJ, J);
111
const double b0 = vertex_coordinates[0];
112
const double b1 = vertex_coordinates[1];
113
const double b2 = vertex_coordinates[2];
115
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
116
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
117
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
127
// Array of basisvalues
128
double basisvalues[3] = {0.0, 0.0, 0.0};
130
// Declare helper variables
131
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
133
// Compute basisvalues
134
basisvalues[0] = 1.0;
135
basisvalues[1] = tmp0;
136
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
137
basisvalues[0] *= std::sqrt(0.5);
138
basisvalues[2] *= std::sqrt(1.0);
139
basisvalues[1] *= std::sqrt(3.0);
141
// Table(s) of coefficients
142
static const double coefficients0[3] = \
143
{0.47140452, -0.28867513, -0.16666667};
146
for (unsigned int r = 0; r < 3; r++)
148
*values += coefficients0[r]*basisvalues[r];
149
}// end loop over 'r'
155
// Array of basisvalues
156
double basisvalues[3] = {0.0, 0.0, 0.0};
158
// Declare helper variables
159
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
161
// Compute basisvalues
162
basisvalues[0] = 1.0;
163
basisvalues[1] = tmp0;
164
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
165
basisvalues[0] *= std::sqrt(0.5);
166
basisvalues[2] *= std::sqrt(1.0);
167
basisvalues[1] *= std::sqrt(3.0);
169
// Table(s) of coefficients
170
static const double coefficients0[3] = \
171
{0.47140452, 0.28867513, -0.16666667};
174
for (unsigned int r = 0; r < 3; r++)
176
*values += coefficients0[r]*basisvalues[r];
177
}// end loop over 'r'
183
// Array of basisvalues
184
double basisvalues[3] = {0.0, 0.0, 0.0};
186
// Declare helper variables
187
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
189
// Compute basisvalues
190
basisvalues[0] = 1.0;
191
basisvalues[1] = tmp0;
192
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
193
basisvalues[0] *= std::sqrt(0.5);
194
basisvalues[2] *= std::sqrt(1.0);
195
basisvalues[1] *= std::sqrt(3.0);
197
// Table(s) of coefficients
198
static const double coefficients0[3] = \
199
{0.47140452, 0.0, 0.33333333};
202
for (unsigned int r = 0; r < 3; r++)
204
*values += coefficients0[r]*basisvalues[r];
205
}// end loop over 'r'
212
/// Evaluate all basis functions at given point x in cell
213
virtual void evaluate_basis_all(double* values,
215
const double* vertex_coordinates,
216
int cell_orientation) const
218
// Helper variable to hold values of a single dof.
219
double dof_values = 0.0;
221
// Loop dofs and call evaluate_basis
222
for (unsigned int r = 0; r < 3; r++)
224
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
225
values[r] = dof_values;
226
}// end loop over 'r'
229
/// Evaluate order n derivatives of basis function i at given point x in cell
230
virtual void evaluate_basis_derivatives(std::size_t i,
234
const double* vertex_coordinates,
235
int cell_orientation) const
239
compute_jacobian_triangle_3d(J, vertex_coordinates);
241
// Compute Jacobian inverse and determinant
244
compute_jacobian_inverse_triangle_3d(K, detJ, J);
247
const double b0 = vertex_coordinates[0];
248
const double b1 = vertex_coordinates[1];
249
const double b2 = vertex_coordinates[2];
251
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
252
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
253
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
256
// Compute number of derivatives.
257
unsigned int num_derivatives_t = 1;
258
for (unsigned int r = 0; r < n; r++)
260
num_derivatives_t *= 2;
261
}// end loop over 'r'
263
// Compute number of derivatives.
264
unsigned int num_derivatives_g = 1;
265
for (unsigned int r = 0; r < n; r++)
267
num_derivatives_g *= 3;
268
}// end loop over 'r'
270
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
271
unsigned int **combinations_t = new unsigned int *[num_derivatives_t];
272
for (unsigned int row = 0; row < num_derivatives_t; row++)
274
combinations_t[row] = new unsigned int [n];
275
for (unsigned int col = 0; col < n; col++)
276
combinations_t[row][col] = 0;
279
// Generate combinations of derivatives
280
for (unsigned int row = 1; row < num_derivatives_t; row++)
282
for (unsigned int num = 0; num < row; num++)
284
for (unsigned int col = n-1; col+1 > 0; col--)
286
if (combinations_t[row][col] + 1 > 1)
287
combinations_t[row][col] = 0;
290
combinations_t[row][col] += 1;
297
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
298
unsigned int **combinations_g = new unsigned int *[num_derivatives_g];
299
for (unsigned int row = 0; row < num_derivatives_g; row++)
301
combinations_g[row] = new unsigned int [n];
302
for (unsigned int col = 0; col < n; col++)
303
combinations_g[row][col] = 0;
306
// Generate combinations of derivatives
307
for (unsigned int row = 1; row < num_derivatives_g; row++)
309
for (unsigned int num = 0; num < row; num++)
311
for (unsigned int col = n-1; col+1 > 0; col--)
313
if (combinations_g[row][col] + 1 > 2)
314
combinations_g[row][col] = 0;
317
combinations_g[row][col] += 1;
324
// Compute inverse of Jacobian
325
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
327
// Declare transformation matrix
328
// Declare pointer to two dimensional array and initialise
329
double **transform = new double *[num_derivatives_g];
331
for (unsigned int j = 0; j < num_derivatives_g; j++)
333
transform[j] = new double [num_derivatives_t];
334
for (unsigned int k = 0; k < num_derivatives_t; k++)
338
// Construct transformation matrix
339
for (unsigned int row = 0; row < num_derivatives_g; row++)
341
for (unsigned int col = 0; col < num_derivatives_t; col++)
343
for (unsigned int k = 0; k < n; k++)
344
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
348
// Reset values. Assuming that values is always an array.
349
for (unsigned int r = 0; r < num_derivatives_g; r++)
352
}// end loop over 'r'
359
// Array of basisvalues
360
double basisvalues[3] = {0.0, 0.0, 0.0};
362
// Declare helper variables
363
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
365
// Compute basisvalues
366
basisvalues[0] = 1.0;
367
basisvalues[1] = tmp0;
368
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
369
basisvalues[0] *= std::sqrt(0.5);
370
basisvalues[2] *= std::sqrt(1.0);
371
basisvalues[1] *= std::sqrt(3.0);
373
// Table(s) of coefficients
374
static const double coefficients0[3] = \
375
{0.47140452, -0.28867513, -0.16666667};
377
// Tables of derivatives of the polynomial base (transpose).
378
static const double dmats0[3][3] = \
380
{4.8989795, 0.0, 0.0},
383
static const double dmats1[3][3] = \
385
{2.4494897, 0.0, 0.0},
386
{4.2426407, 0.0, 0.0}};
388
// Compute reference derivatives.
389
// Declare pointer to array of derivatives on FIAT element.
390
double *derivatives = new double[num_derivatives_t];
391
for (unsigned int r = 0; r < num_derivatives_t; r++)
393
derivatives[r] = 0.0;
394
}// end loop over 'r'
396
// Declare derivative matrix (of polynomial basis).
397
double dmats[3][3] = \
402
// Declare (auxiliary) derivative matrix (of polynomial basis).
403
double dmats_old[3][3] = \
408
// Loop possible derivatives.
409
for (unsigned int r = 0; r < num_derivatives_t; r++)
411
// Resetting dmats values to compute next derivative.
412
for (unsigned int t = 0; t < 3; t++)
414
for (unsigned int u = 0; u < 3; u++)
422
}// end loop over 'u'
423
}// end loop over 't'
425
// Looping derivative order to generate dmats.
426
for (unsigned int s = 0; s < n; s++)
428
// Updating dmats_old with new values and resetting dmats.
429
for (unsigned int t = 0; t < 3; t++)
431
for (unsigned int u = 0; u < 3; u++)
433
dmats_old[t][u] = dmats[t][u];
435
}// end loop over 'u'
436
}// end loop over 't'
438
// Update dmats using an inner product.
439
if (combinations_t[r][s] == 0)
441
for (unsigned int t = 0; t < 3; t++)
443
for (unsigned int u = 0; u < 3; u++)
445
for (unsigned int tu = 0; tu < 3; tu++)
447
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
448
}// end loop over 'tu'
449
}// end loop over 'u'
450
}// end loop over 't'
453
if (combinations_t[r][s] == 1)
455
for (unsigned int t = 0; t < 3; t++)
457
for (unsigned int u = 0; u < 3; u++)
459
for (unsigned int tu = 0; tu < 3; tu++)
461
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
462
}// end loop over 'tu'
463
}// end loop over 'u'
464
}// end loop over 't'
467
}// end loop over 's'
468
for (unsigned int s = 0; s < 3; s++)
470
for (unsigned int t = 0; t < 3; t++)
472
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
473
}// end loop over 't'
474
}// end loop over 's'
475
}// end loop over 'r'
477
// Transform derivatives back to physical element
478
for (unsigned int r = 0; r < num_derivatives_g; r++)
480
for (unsigned int s = 0; s < num_derivatives_t; s++)
482
values[r] += transform[r][s]*derivatives[s];
483
}// end loop over 's'
484
}// end loop over 'r'
486
// Delete pointer to array of derivatives on FIAT element
487
delete [] derivatives;
489
// Delete pointer to array of combinations of derivatives and transform
490
for (unsigned int r = 0; r < num_derivatives_t; r++)
492
delete [] combinations_t[r];
493
}// end loop over 'r'
494
delete [] combinations_t;
495
for (unsigned int r = 0; r < num_derivatives_g; r++)
497
delete [] combinations_g[r];
498
}// end loop over 'r'
499
delete [] combinations_g;
500
for (unsigned int r = 0; r < num_derivatives_g; r++)
502
delete [] transform[r];
503
}// end loop over 'r'
510
// Array of basisvalues
511
double basisvalues[3] = {0.0, 0.0, 0.0};
513
// Declare helper variables
514
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
516
// Compute basisvalues
517
basisvalues[0] = 1.0;
518
basisvalues[1] = tmp0;
519
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
520
basisvalues[0] *= std::sqrt(0.5);
521
basisvalues[2] *= std::sqrt(1.0);
522
basisvalues[1] *= std::sqrt(3.0);
524
// Table(s) of coefficients
525
static const double coefficients0[3] = \
526
{0.47140452, 0.28867513, -0.16666667};
528
// Tables of derivatives of the polynomial base (transpose).
529
static const double dmats0[3][3] = \
531
{4.8989795, 0.0, 0.0},
534
static const double dmats1[3][3] = \
536
{2.4494897, 0.0, 0.0},
537
{4.2426407, 0.0, 0.0}};
539
// Compute reference derivatives.
540
// Declare pointer to array of derivatives on FIAT element.
541
double *derivatives = new double[num_derivatives_t];
542
for (unsigned int r = 0; r < num_derivatives_t; r++)
544
derivatives[r] = 0.0;
545
}// end loop over 'r'
547
// Declare derivative matrix (of polynomial basis).
548
double dmats[3][3] = \
553
// Declare (auxiliary) derivative matrix (of polynomial basis).
554
double dmats_old[3][3] = \
559
// Loop possible derivatives.
560
for (unsigned int r = 0; r < num_derivatives_t; r++)
562
// Resetting dmats values to compute next derivative.
563
for (unsigned int t = 0; t < 3; t++)
565
for (unsigned int u = 0; u < 3; u++)
573
}// end loop over 'u'
574
}// end loop over 't'
576
// Looping derivative order to generate dmats.
577
for (unsigned int s = 0; s < n; s++)
579
// Updating dmats_old with new values and resetting dmats.
580
for (unsigned int t = 0; t < 3; t++)
582
for (unsigned int u = 0; u < 3; u++)
584
dmats_old[t][u] = dmats[t][u];
586
}// end loop over 'u'
587
}// end loop over 't'
589
// Update dmats using an inner product.
590
if (combinations_t[r][s] == 0)
592
for (unsigned int t = 0; t < 3; t++)
594
for (unsigned int u = 0; u < 3; u++)
596
for (unsigned int tu = 0; tu < 3; tu++)
598
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
599
}// end loop over 'tu'
600
}// end loop over 'u'
601
}// end loop over 't'
604
if (combinations_t[r][s] == 1)
606
for (unsigned int t = 0; t < 3; t++)
608
for (unsigned int u = 0; u < 3; u++)
610
for (unsigned int tu = 0; tu < 3; tu++)
612
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
613
}// end loop over 'tu'
614
}// end loop over 'u'
615
}// end loop over 't'
618
}// end loop over 's'
619
for (unsigned int s = 0; s < 3; s++)
621
for (unsigned int t = 0; t < 3; t++)
623
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
624
}// end loop over 't'
625
}// end loop over 's'
626
}// end loop over 'r'
628
// Transform derivatives back to physical element
629
for (unsigned int r = 0; r < num_derivatives_g; r++)
631
for (unsigned int s = 0; s < num_derivatives_t; s++)
633
values[r] += transform[r][s]*derivatives[s];
634
}// end loop over 's'
635
}// end loop over 'r'
637
// Delete pointer to array of derivatives on FIAT element
638
delete [] derivatives;
640
// Delete pointer to array of combinations of derivatives and transform
641
for (unsigned int r = 0; r < num_derivatives_t; r++)
643
delete [] combinations_t[r];
644
}// end loop over 'r'
645
delete [] combinations_t;
646
for (unsigned int r = 0; r < num_derivatives_g; r++)
648
delete [] combinations_g[r];
649
}// end loop over 'r'
650
delete [] combinations_g;
651
for (unsigned int r = 0; r < num_derivatives_g; r++)
653
delete [] transform[r];
654
}// end loop over 'r'
661
// Array of basisvalues
662
double basisvalues[3] = {0.0, 0.0, 0.0};
664
// Declare helper variables
665
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
667
// Compute basisvalues
668
basisvalues[0] = 1.0;
669
basisvalues[1] = tmp0;
670
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
671
basisvalues[0] *= std::sqrt(0.5);
672
basisvalues[2] *= std::sqrt(1.0);
673
basisvalues[1] *= std::sqrt(3.0);
675
// Table(s) of coefficients
676
static const double coefficients0[3] = \
677
{0.47140452, 0.0, 0.33333333};
679
// Tables of derivatives of the polynomial base (transpose).
680
static const double dmats0[3][3] = \
682
{4.8989795, 0.0, 0.0},
685
static const double dmats1[3][3] = \
687
{2.4494897, 0.0, 0.0},
688
{4.2426407, 0.0, 0.0}};
690
// Compute reference derivatives.
691
// Declare pointer to array of derivatives on FIAT element.
692
double *derivatives = new double[num_derivatives_t];
693
for (unsigned int r = 0; r < num_derivatives_t; r++)
695
derivatives[r] = 0.0;
696
}// end loop over 'r'
698
// Declare derivative matrix (of polynomial basis).
699
double dmats[3][3] = \
704
// Declare (auxiliary) derivative matrix (of polynomial basis).
705
double dmats_old[3][3] = \
710
// Loop possible derivatives.
711
for (unsigned int r = 0; r < num_derivatives_t; r++)
713
// Resetting dmats values to compute next derivative.
714
for (unsigned int t = 0; t < 3; t++)
716
for (unsigned int u = 0; u < 3; u++)
724
}// end loop over 'u'
725
}// end loop over 't'
727
// Looping derivative order to generate dmats.
728
for (unsigned int s = 0; s < n; s++)
730
// Updating dmats_old with new values and resetting dmats.
731
for (unsigned int t = 0; t < 3; t++)
733
for (unsigned int u = 0; u < 3; u++)
735
dmats_old[t][u] = dmats[t][u];
737
}// end loop over 'u'
738
}// end loop over 't'
740
// Update dmats using an inner product.
741
if (combinations_t[r][s] == 0)
743
for (unsigned int t = 0; t < 3; t++)
745
for (unsigned int u = 0; u < 3; u++)
747
for (unsigned int tu = 0; tu < 3; tu++)
749
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
750
}// end loop over 'tu'
751
}// end loop over 'u'
752
}// end loop over 't'
755
if (combinations_t[r][s] == 1)
757
for (unsigned int t = 0; t < 3; t++)
759
for (unsigned int u = 0; u < 3; u++)
761
for (unsigned int tu = 0; tu < 3; tu++)
763
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
764
}// end loop over 'tu'
765
}// end loop over 'u'
766
}// end loop over 't'
769
}// end loop over 's'
770
for (unsigned int s = 0; s < 3; s++)
772
for (unsigned int t = 0; t < 3; t++)
774
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
775
}// end loop over 't'
776
}// end loop over 's'
777
}// end loop over 'r'
779
// Transform derivatives back to physical element
780
for (unsigned int r = 0; r < num_derivatives_g; r++)
782
for (unsigned int s = 0; s < num_derivatives_t; s++)
784
values[r] += transform[r][s]*derivatives[s];
785
}// end loop over 's'
786
}// end loop over 'r'
788
// Delete pointer to array of derivatives on FIAT element
789
delete [] derivatives;
791
// Delete pointer to array of combinations of derivatives and transform
792
for (unsigned int r = 0; r < num_derivatives_t; r++)
794
delete [] combinations_t[r];
795
}// end loop over 'r'
796
delete [] combinations_t;
797
for (unsigned int r = 0; r < num_derivatives_g; r++)
799
delete [] combinations_g[r];
800
}// end loop over 'r'
801
delete [] combinations_g;
802
for (unsigned int r = 0; r < num_derivatives_g; r++)
804
delete [] transform[r];
805
}// end loop over 'r'
813
/// Evaluate order n derivatives of all basis functions at given point x in cell
814
virtual void evaluate_basis_derivatives_all(std::size_t n,
817
const double* vertex_coordinates,
818
int cell_orientation) const
820
// Compute number of derivatives.
821
unsigned int num_derivatives_g = 1;
822
for (unsigned int r = 0; r < n; r++)
824
num_derivatives_g *= 3;
825
}// end loop over 'r'
827
// Helper variable to hold values of a single dof.
828
double *dof_values = new double[num_derivatives_g];
829
for (unsigned int r = 0; r < num_derivatives_g; r++)
832
}// end loop over 'r'
834
// Loop dofs and call evaluate_basis_derivatives.
835
for (unsigned int r = 0; r < 3; r++)
837
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
838
for (unsigned int s = 0; s < num_derivatives_g; s++)
840
values[r*num_derivatives_g + s] = dof_values[s];
841
}// end loop over 's'
842
}// end loop over 'r'
845
delete [] dof_values;
848
/// Evaluate linear functional for dof i on the function f
849
virtual double evaluate_dof(std::size_t i,
850
const ufc::function& f,
851
const double* vertex_coordinates,
852
int cell_orientation,
853
const ufc::cell& c) const
855
// Declare variables for result of evaluation
858
// Declare variable for physical coordinates
864
y[0] = vertex_coordinates[0];
865
y[1] = vertex_coordinates[1];
866
y[2] = vertex_coordinates[2];
867
f.evaluate(vals, y, c);
873
y[0] = vertex_coordinates[3];
874
y[1] = vertex_coordinates[4];
875
y[2] = vertex_coordinates[5];
876
f.evaluate(vals, y, c);
882
y[0] = vertex_coordinates[6];
883
y[1] = vertex_coordinates[7];
884
y[2] = vertex_coordinates[8];
885
f.evaluate(vals, y, c);
894
/// Evaluate linear functionals for all dofs on the function f
895
virtual void evaluate_dofs(double* values,
896
const ufc::function& f,
897
const double* vertex_coordinates,
898
int cell_orientation,
899
const ufc::cell& c) const
901
// Declare variables for result of evaluation
904
// Declare variable for physical coordinates
906
y[0] = vertex_coordinates[0];
907
y[1] = vertex_coordinates[1];
908
y[2] = vertex_coordinates[2];
909
f.evaluate(vals, y, c);
911
y[0] = vertex_coordinates[3];
912
y[1] = vertex_coordinates[4];
913
y[2] = vertex_coordinates[5];
914
f.evaluate(vals, y, c);
916
y[0] = vertex_coordinates[6];
917
y[1] = vertex_coordinates[7];
918
y[2] = vertex_coordinates[8];
919
f.evaluate(vals, y, c);
923
/// Interpolate vertex values from dof values
924
virtual void interpolate_vertex_values(double* vertex_values,
925
const double* dof_values,
926
const double* vertex_coordinates,
927
int cell_orientation,
928
const ufc::cell& c) const
930
// Evaluate function and change variables
931
vertex_values[0] = dof_values[0];
932
vertex_values[1] = dof_values[1];
933
vertex_values[2] = dof_values[2];
936
/// Map coordinate xhat from reference cell to coordinate x in cell
937
virtual void map_from_reference_cell(double* x,
939
const ufc::cell& c) const
941
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
944
/// Map from coordinate x in cell to coordinate xhat in reference cell
945
virtual void map_to_reference_cell(double* xhat,
947
const ufc::cell& c) const
949
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
952
/// Return the number of sub elements (for a mixed element)
953
virtual std::size_t num_sub_elements() const
958
/// Create a new finite element for sub element i (for a mixed element)
959
virtual ufc::finite_element* create_sub_element(std::size_t i) const
964
/// Create a new class instance
965
virtual ufc::finite_element* create() const
967
return new x_element5_finite_element_0();
972
/// This class defines the interface for a finite element.
974
class x_element5_finite_element_1: public ufc::finite_element
979
x_element5_finite_element_1() : ufc::finite_element()
985
virtual ~x_element5_finite_element_1()
990
/// Return a string identifying the finite element
991
virtual const char* signature() const
993
return "VectorElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, 3, None)";
996
/// Return the cell shape
997
virtual ufc::shape cell_shape() const
999
return ufc::triangle;
1002
/// Return the topological dimension of the cell shape
1003
virtual std::size_t topological_dimension() const
1008
/// Return the geometric dimension of the cell shape
1009
virtual std::size_t geometric_dimension() const
1014
/// Return the dimension of the finite element function space
1015
virtual std::size_t space_dimension() const
1020
/// Return the rank of the value space
1021
virtual std::size_t value_rank() const
1026
/// Return the dimension of the value space for axis i
1027
virtual std::size_t value_dimension(std::size_t i) const
1041
/// Evaluate basis function i at given point x in cell
1042
virtual void evaluate_basis(std::size_t i,
1045
const double* vertex_coordinates,
1046
int cell_orientation) const
1050
compute_jacobian_triangle_3d(J, vertex_coordinates);
1052
// Compute Jacobian inverse and determinant
1055
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1058
const double b0 = vertex_coordinates[0];
1059
const double b1 = vertex_coordinates[1];
1060
const double b2 = vertex_coordinates[2];
1062
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1063
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
1064
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
1076
// Array of basisvalues
1077
double basisvalues[3] = {0.0, 0.0, 0.0};
1079
// Declare helper variables
1080
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1082
// Compute basisvalues
1083
basisvalues[0] = 1.0;
1084
basisvalues[1] = tmp0;
1085
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1086
basisvalues[0] *= std::sqrt(0.5);
1087
basisvalues[2] *= std::sqrt(1.0);
1088
basisvalues[1] *= std::sqrt(3.0);
1090
// Table(s) of coefficients
1091
static const double coefficients0[3] = \
1092
{0.47140452, -0.28867513, -0.16666667};
1095
for (unsigned int r = 0; r < 3; r++)
1097
values[0] += coefficients0[r]*basisvalues[r];
1098
}// end loop over 'r'
1104
// Array of basisvalues
1105
double basisvalues[3] = {0.0, 0.0, 0.0};
1107
// Declare helper variables
1108
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1110
// Compute basisvalues
1111
basisvalues[0] = 1.0;
1112
basisvalues[1] = tmp0;
1113
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1114
basisvalues[0] *= std::sqrt(0.5);
1115
basisvalues[2] *= std::sqrt(1.0);
1116
basisvalues[1] *= std::sqrt(3.0);
1118
// Table(s) of coefficients
1119
static const double coefficients0[3] = \
1120
{0.47140452, 0.28867513, -0.16666667};
1123
for (unsigned int r = 0; r < 3; r++)
1125
values[0] += coefficients0[r]*basisvalues[r];
1126
}// end loop over 'r'
1132
// Array of basisvalues
1133
double basisvalues[3] = {0.0, 0.0, 0.0};
1135
// Declare helper variables
1136
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1138
// Compute basisvalues
1139
basisvalues[0] = 1.0;
1140
basisvalues[1] = tmp0;
1141
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1142
basisvalues[0] *= std::sqrt(0.5);
1143
basisvalues[2] *= std::sqrt(1.0);
1144
basisvalues[1] *= std::sqrt(3.0);
1146
// Table(s) of coefficients
1147
static const double coefficients0[3] = \
1148
{0.47140452, 0.0, 0.33333333};
1151
for (unsigned int r = 0; r < 3; r++)
1153
values[0] += coefficients0[r]*basisvalues[r];
1154
}// end loop over 'r'
1160
// Array of basisvalues
1161
double basisvalues[3] = {0.0, 0.0, 0.0};
1163
// Declare helper variables
1164
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1166
// Compute basisvalues
1167
basisvalues[0] = 1.0;
1168
basisvalues[1] = tmp0;
1169
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1170
basisvalues[0] *= std::sqrt(0.5);
1171
basisvalues[2] *= std::sqrt(1.0);
1172
basisvalues[1] *= std::sqrt(3.0);
1174
// Table(s) of coefficients
1175
static const double coefficients0[3] = \
1176
{0.47140452, -0.28867513, -0.16666667};
1179
for (unsigned int r = 0; r < 3; r++)
1181
values[1] += coefficients0[r]*basisvalues[r];
1182
}// end loop over 'r'
1188
// Array of basisvalues
1189
double basisvalues[3] = {0.0, 0.0, 0.0};
1191
// Declare helper variables
1192
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1194
// Compute basisvalues
1195
basisvalues[0] = 1.0;
1196
basisvalues[1] = tmp0;
1197
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1198
basisvalues[0] *= std::sqrt(0.5);
1199
basisvalues[2] *= std::sqrt(1.0);
1200
basisvalues[1] *= std::sqrt(3.0);
1202
// Table(s) of coefficients
1203
static const double coefficients0[3] = \
1204
{0.47140452, 0.28867513, -0.16666667};
1207
for (unsigned int r = 0; r < 3; r++)
1209
values[1] += coefficients0[r]*basisvalues[r];
1210
}// end loop over 'r'
1216
// Array of basisvalues
1217
double basisvalues[3] = {0.0, 0.0, 0.0};
1219
// Declare helper variables
1220
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1222
// Compute basisvalues
1223
basisvalues[0] = 1.0;
1224
basisvalues[1] = tmp0;
1225
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1226
basisvalues[0] *= std::sqrt(0.5);
1227
basisvalues[2] *= std::sqrt(1.0);
1228
basisvalues[1] *= std::sqrt(3.0);
1230
// Table(s) of coefficients
1231
static const double coefficients0[3] = \
1232
{0.47140452, 0.0, 0.33333333};
1235
for (unsigned int r = 0; r < 3; r++)
1237
values[1] += coefficients0[r]*basisvalues[r];
1238
}// end loop over 'r'
1244
// Array of basisvalues
1245
double basisvalues[3] = {0.0, 0.0, 0.0};
1247
// Declare helper variables
1248
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1250
// Compute basisvalues
1251
basisvalues[0] = 1.0;
1252
basisvalues[1] = tmp0;
1253
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1254
basisvalues[0] *= std::sqrt(0.5);
1255
basisvalues[2] *= std::sqrt(1.0);
1256
basisvalues[1] *= std::sqrt(3.0);
1258
// Table(s) of coefficients
1259
static const double coefficients0[3] = \
1260
{0.47140452, -0.28867513, -0.16666667};
1263
for (unsigned int r = 0; r < 3; r++)
1265
values[2] += coefficients0[r]*basisvalues[r];
1266
}// end loop over 'r'
1272
// Array of basisvalues
1273
double basisvalues[3] = {0.0, 0.0, 0.0};
1275
// Declare helper variables
1276
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1278
// Compute basisvalues
1279
basisvalues[0] = 1.0;
1280
basisvalues[1] = tmp0;
1281
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1282
basisvalues[0] *= std::sqrt(0.5);
1283
basisvalues[2] *= std::sqrt(1.0);
1284
basisvalues[1] *= std::sqrt(3.0);
1286
// Table(s) of coefficients
1287
static const double coefficients0[3] = \
1288
{0.47140452, 0.28867513, -0.16666667};
1291
for (unsigned int r = 0; r < 3; r++)
1293
values[2] += coefficients0[r]*basisvalues[r];
1294
}// end loop over 'r'
1300
// Array of basisvalues
1301
double basisvalues[3] = {0.0, 0.0, 0.0};
1303
// Declare helper variables
1304
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1306
// Compute basisvalues
1307
basisvalues[0] = 1.0;
1308
basisvalues[1] = tmp0;
1309
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1310
basisvalues[0] *= std::sqrt(0.5);
1311
basisvalues[2] *= std::sqrt(1.0);
1312
basisvalues[1] *= std::sqrt(3.0);
1314
// Table(s) of coefficients
1315
static const double coefficients0[3] = \
1316
{0.47140452, 0.0, 0.33333333};
1319
for (unsigned int r = 0; r < 3; r++)
1321
values[2] += coefficients0[r]*basisvalues[r];
1322
}// end loop over 'r'
1329
/// Evaluate all basis functions at given point x in cell
1330
virtual void evaluate_basis_all(double* values,
1332
const double* vertex_coordinates,
1333
int cell_orientation) const
1335
// Helper variable to hold values of a single dof.
1336
double dof_values[3] = {0.0, 0.0, 0.0};
1338
// Loop dofs and call evaluate_basis
1339
for (unsigned int r = 0; r < 9; r++)
1341
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
1342
for (unsigned int s = 0; s < 3; s++)
1344
values[r*3 + s] = dof_values[s];
1345
}// end loop over 's'
1346
}// end loop over 'r'
1349
/// Evaluate order n derivatives of basis function i at given point x in cell
1350
virtual void evaluate_basis_derivatives(std::size_t i,
1354
const double* vertex_coordinates,
1355
int cell_orientation) const
1359
compute_jacobian_triangle_3d(J, vertex_coordinates);
1361
// Compute Jacobian inverse and determinant
1364
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1367
const double b0 = vertex_coordinates[0];
1368
const double b1 = vertex_coordinates[1];
1369
const double b2 = vertex_coordinates[2];
1371
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1372
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
1373
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
1376
// Compute number of derivatives.
1377
unsigned int num_derivatives_t = 1;
1378
for (unsigned int r = 0; r < n; r++)
1380
num_derivatives_t *= 2;
1381
}// end loop over 'r'
1383
// Compute number of derivatives.
1384
unsigned int num_derivatives_g = 1;
1385
for (unsigned int r = 0; r < n; r++)
1387
num_derivatives_g *= 3;
1388
}// end loop over 'r'
1390
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1391
unsigned int **combinations_t = new unsigned int *[num_derivatives_t];
1392
for (unsigned int row = 0; row < num_derivatives_t; row++)
1394
combinations_t[row] = new unsigned int [n];
1395
for (unsigned int col = 0; col < n; col++)
1396
combinations_t[row][col] = 0;
1399
// Generate combinations of derivatives
1400
for (unsigned int row = 1; row < num_derivatives_t; row++)
1402
for (unsigned int num = 0; num < row; num++)
1404
for (unsigned int col = n-1; col+1 > 0; col--)
1406
if (combinations_t[row][col] + 1 > 1)
1407
combinations_t[row][col] = 0;
1410
combinations_t[row][col] += 1;
1417
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1418
unsigned int **combinations_g = new unsigned int *[num_derivatives_g];
1419
for (unsigned int row = 0; row < num_derivatives_g; row++)
1421
combinations_g[row] = new unsigned int [n];
1422
for (unsigned int col = 0; col < n; col++)
1423
combinations_g[row][col] = 0;
1426
// Generate combinations of derivatives
1427
for (unsigned int row = 1; row < num_derivatives_g; row++)
1429
for (unsigned int num = 0; num < row; num++)
1431
for (unsigned int col = n-1; col+1 > 0; col--)
1433
if (combinations_g[row][col] + 1 > 2)
1434
combinations_g[row][col] = 0;
1437
combinations_g[row][col] += 1;
1444
// Compute inverse of Jacobian
1445
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
1447
// Declare transformation matrix
1448
// Declare pointer to two dimensional array and initialise
1449
double **transform = new double *[num_derivatives_g];
1451
for (unsigned int j = 0; j < num_derivatives_g; j++)
1453
transform[j] = new double [num_derivatives_t];
1454
for (unsigned int k = 0; k < num_derivatives_t; k++)
1455
transform[j][k] = 1;
1458
// Construct transformation matrix
1459
for (unsigned int row = 0; row < num_derivatives_g; row++)
1461
for (unsigned int col = 0; col < num_derivatives_t; col++)
1463
for (unsigned int k = 0; k < n; k++)
1464
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
1468
// Reset values. Assuming that values is always an array.
1469
for (unsigned int r = 0; r < 3*num_derivatives_g; r++)
1472
}// end loop over 'r'
1479
// Array of basisvalues
1480
double basisvalues[3] = {0.0, 0.0, 0.0};
1482
// Declare helper variables
1483
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1485
// Compute basisvalues
1486
basisvalues[0] = 1.0;
1487
basisvalues[1] = tmp0;
1488
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1489
basisvalues[0] *= std::sqrt(0.5);
1490
basisvalues[2] *= std::sqrt(1.0);
1491
basisvalues[1] *= std::sqrt(3.0);
1493
// Table(s) of coefficients
1494
static const double coefficients0[3] = \
1495
{0.47140452, -0.28867513, -0.16666667};
1497
// Tables of derivatives of the polynomial base (transpose).
1498
static const double dmats0[3][3] = \
1500
{4.8989795, 0.0, 0.0},
1503
static const double dmats1[3][3] = \
1505
{2.4494897, 0.0, 0.0},
1506
{4.2426407, 0.0, 0.0}};
1508
// Compute reference derivatives.
1509
// Declare pointer to array of derivatives on FIAT element.
1510
double *derivatives = new double[num_derivatives_t];
1511
for (unsigned int r = 0; r < num_derivatives_t; r++)
1513
derivatives[r] = 0.0;
1514
}// end loop over 'r'
1516
// Declare derivative matrix (of polynomial basis).
1517
double dmats[3][3] = \
1522
// Declare (auxiliary) derivative matrix (of polynomial basis).
1523
double dmats_old[3][3] = \
1528
// Loop possible derivatives.
1529
for (unsigned int r = 0; r < num_derivatives_t; r++)
1531
// Resetting dmats values to compute next derivative.
1532
for (unsigned int t = 0; t < 3; t++)
1534
for (unsigned int u = 0; u < 3; u++)
1542
}// end loop over 'u'
1543
}// end loop over 't'
1545
// Looping derivative order to generate dmats.
1546
for (unsigned int s = 0; s < n; s++)
1548
// Updating dmats_old with new values and resetting dmats.
1549
for (unsigned int t = 0; t < 3; t++)
1551
for (unsigned int u = 0; u < 3; u++)
1553
dmats_old[t][u] = dmats[t][u];
1555
}// end loop over 'u'
1556
}// end loop over 't'
1558
// Update dmats using an inner product.
1559
if (combinations_t[r][s] == 0)
1561
for (unsigned int t = 0; t < 3; t++)
1563
for (unsigned int u = 0; u < 3; u++)
1565
for (unsigned int tu = 0; tu < 3; tu++)
1567
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1568
}// end loop over 'tu'
1569
}// end loop over 'u'
1570
}// end loop over 't'
1573
if (combinations_t[r][s] == 1)
1575
for (unsigned int t = 0; t < 3; t++)
1577
for (unsigned int u = 0; u < 3; u++)
1579
for (unsigned int tu = 0; tu < 3; tu++)
1581
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1582
}// end loop over 'tu'
1583
}// end loop over 'u'
1584
}// end loop over 't'
1587
}// end loop over 's'
1588
for (unsigned int s = 0; s < 3; s++)
1590
for (unsigned int t = 0; t < 3; t++)
1592
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1593
}// end loop over 't'
1594
}// end loop over 's'
1595
}// end loop over 'r'
1597
// Transform derivatives back to physical element
1598
for (unsigned int r = 0; r < num_derivatives_g; r++)
1600
for (unsigned int s = 0; s < num_derivatives_t; s++)
1602
values[r] += transform[r][s]*derivatives[s];
1603
}// end loop over 's'
1604
}// end loop over 'r'
1606
// Delete pointer to array of derivatives on FIAT element
1607
delete [] derivatives;
1609
// Delete pointer to array of combinations of derivatives and transform
1610
for (unsigned int r = 0; r < num_derivatives_t; r++)
1612
delete [] combinations_t[r];
1613
}// end loop over 'r'
1614
delete [] combinations_t;
1615
for (unsigned int r = 0; r < num_derivatives_g; r++)
1617
delete [] combinations_g[r];
1618
}// end loop over 'r'
1619
delete [] combinations_g;
1620
for (unsigned int r = 0; r < num_derivatives_g; r++)
1622
delete [] transform[r];
1623
}// end loop over 'r'
1624
delete [] transform;
1630
// Array of basisvalues
1631
double basisvalues[3] = {0.0, 0.0, 0.0};
1633
// Declare helper variables
1634
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1636
// Compute basisvalues
1637
basisvalues[0] = 1.0;
1638
basisvalues[1] = tmp0;
1639
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1640
basisvalues[0] *= std::sqrt(0.5);
1641
basisvalues[2] *= std::sqrt(1.0);
1642
basisvalues[1] *= std::sqrt(3.0);
1644
// Table(s) of coefficients
1645
static const double coefficients0[3] = \
1646
{0.47140452, 0.28867513, -0.16666667};
1648
// Tables of derivatives of the polynomial base (transpose).
1649
static const double dmats0[3][3] = \
1651
{4.8989795, 0.0, 0.0},
1654
static const double dmats1[3][3] = \
1656
{2.4494897, 0.0, 0.0},
1657
{4.2426407, 0.0, 0.0}};
1659
// Compute reference derivatives.
1660
// Declare pointer to array of derivatives on FIAT element.
1661
double *derivatives = new double[num_derivatives_t];
1662
for (unsigned int r = 0; r < num_derivatives_t; r++)
1664
derivatives[r] = 0.0;
1665
}// end loop over 'r'
1667
// Declare derivative matrix (of polynomial basis).
1668
double dmats[3][3] = \
1673
// Declare (auxiliary) derivative matrix (of polynomial basis).
1674
double dmats_old[3][3] = \
1679
// Loop possible derivatives.
1680
for (unsigned int r = 0; r < num_derivatives_t; r++)
1682
// Resetting dmats values to compute next derivative.
1683
for (unsigned int t = 0; t < 3; t++)
1685
for (unsigned int u = 0; u < 3; u++)
1693
}// end loop over 'u'
1694
}// end loop over 't'
1696
// Looping derivative order to generate dmats.
1697
for (unsigned int s = 0; s < n; s++)
1699
// Updating dmats_old with new values and resetting dmats.
1700
for (unsigned int t = 0; t < 3; t++)
1702
for (unsigned int u = 0; u < 3; u++)
1704
dmats_old[t][u] = dmats[t][u];
1706
}// end loop over 'u'
1707
}// end loop over 't'
1709
// Update dmats using an inner product.
1710
if (combinations_t[r][s] == 0)
1712
for (unsigned int t = 0; t < 3; t++)
1714
for (unsigned int u = 0; u < 3; u++)
1716
for (unsigned int tu = 0; tu < 3; tu++)
1718
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1719
}// end loop over 'tu'
1720
}// end loop over 'u'
1721
}// end loop over 't'
1724
if (combinations_t[r][s] == 1)
1726
for (unsigned int t = 0; t < 3; t++)
1728
for (unsigned int u = 0; u < 3; u++)
1730
for (unsigned int tu = 0; tu < 3; tu++)
1732
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1733
}// end loop over 'tu'
1734
}// end loop over 'u'
1735
}// end loop over 't'
1738
}// end loop over 's'
1739
for (unsigned int s = 0; s < 3; s++)
1741
for (unsigned int t = 0; t < 3; t++)
1743
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1744
}// end loop over 't'
1745
}// end loop over 's'
1746
}// end loop over 'r'
1748
// Transform derivatives back to physical element
1749
for (unsigned int r = 0; r < num_derivatives_g; r++)
1751
for (unsigned int s = 0; s < num_derivatives_t; s++)
1753
values[r] += transform[r][s]*derivatives[s];
1754
}// end loop over 's'
1755
}// end loop over 'r'
1757
// Delete pointer to array of derivatives on FIAT element
1758
delete [] derivatives;
1760
// Delete pointer to array of combinations of derivatives and transform
1761
for (unsigned int r = 0; r < num_derivatives_t; r++)
1763
delete [] combinations_t[r];
1764
}// end loop over 'r'
1765
delete [] combinations_t;
1766
for (unsigned int r = 0; r < num_derivatives_g; r++)
1768
delete [] combinations_g[r];
1769
}// end loop over 'r'
1770
delete [] combinations_g;
1771
for (unsigned int r = 0; r < num_derivatives_g; r++)
1773
delete [] transform[r];
1774
}// end loop over 'r'
1775
delete [] transform;
1781
// Array of basisvalues
1782
double basisvalues[3] = {0.0, 0.0, 0.0};
1784
// Declare helper variables
1785
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1787
// Compute basisvalues
1788
basisvalues[0] = 1.0;
1789
basisvalues[1] = tmp0;
1790
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1791
basisvalues[0] *= std::sqrt(0.5);
1792
basisvalues[2] *= std::sqrt(1.0);
1793
basisvalues[1] *= std::sqrt(3.0);
1795
// Table(s) of coefficients
1796
static const double coefficients0[3] = \
1797
{0.47140452, 0.0, 0.33333333};
1799
// Tables of derivatives of the polynomial base (transpose).
1800
static const double dmats0[3][3] = \
1802
{4.8989795, 0.0, 0.0},
1805
static const double dmats1[3][3] = \
1807
{2.4494897, 0.0, 0.0},
1808
{4.2426407, 0.0, 0.0}};
1810
// Compute reference derivatives.
1811
// Declare pointer to array of derivatives on FIAT element.
1812
double *derivatives = new double[num_derivatives_t];
1813
for (unsigned int r = 0; r < num_derivatives_t; r++)
1815
derivatives[r] = 0.0;
1816
}// end loop over 'r'
1818
// Declare derivative matrix (of polynomial basis).
1819
double dmats[3][3] = \
1824
// Declare (auxiliary) derivative matrix (of polynomial basis).
1825
double dmats_old[3][3] = \
1830
// Loop possible derivatives.
1831
for (unsigned int r = 0; r < num_derivatives_t; r++)
1833
// Resetting dmats values to compute next derivative.
1834
for (unsigned int t = 0; t < 3; t++)
1836
for (unsigned int u = 0; u < 3; u++)
1844
}// end loop over 'u'
1845
}// end loop over 't'
1847
// Looping derivative order to generate dmats.
1848
for (unsigned int s = 0; s < n; s++)
1850
// Updating dmats_old with new values and resetting dmats.
1851
for (unsigned int t = 0; t < 3; t++)
1853
for (unsigned int u = 0; u < 3; u++)
1855
dmats_old[t][u] = dmats[t][u];
1857
}// end loop over 'u'
1858
}// end loop over 't'
1860
// Update dmats using an inner product.
1861
if (combinations_t[r][s] == 0)
1863
for (unsigned int t = 0; t < 3; t++)
1865
for (unsigned int u = 0; u < 3; u++)
1867
for (unsigned int tu = 0; tu < 3; tu++)
1869
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1870
}// end loop over 'tu'
1871
}// end loop over 'u'
1872
}// end loop over 't'
1875
if (combinations_t[r][s] == 1)
1877
for (unsigned int t = 0; t < 3; t++)
1879
for (unsigned int u = 0; u < 3; u++)
1881
for (unsigned int tu = 0; tu < 3; tu++)
1883
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1884
}// end loop over 'tu'
1885
}// end loop over 'u'
1886
}// end loop over 't'
1889
}// end loop over 's'
1890
for (unsigned int s = 0; s < 3; s++)
1892
for (unsigned int t = 0; t < 3; t++)
1894
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1895
}// end loop over 't'
1896
}// end loop over 's'
1897
}// end loop over 'r'
1899
// Transform derivatives back to physical element
1900
for (unsigned int r = 0; r < num_derivatives_g; r++)
1902
for (unsigned int s = 0; s < num_derivatives_t; s++)
1904
values[r] += transform[r][s]*derivatives[s];
1905
}// end loop over 's'
1906
}// end loop over 'r'
1908
// Delete pointer to array of derivatives on FIAT element
1909
delete [] derivatives;
1911
// Delete pointer to array of combinations of derivatives and transform
1912
for (unsigned int r = 0; r < num_derivatives_t; r++)
1914
delete [] combinations_t[r];
1915
}// end loop over 'r'
1916
delete [] combinations_t;
1917
for (unsigned int r = 0; r < num_derivatives_g; r++)
1919
delete [] combinations_g[r];
1920
}// end loop over 'r'
1921
delete [] combinations_g;
1922
for (unsigned int r = 0; r < num_derivatives_g; r++)
1924
delete [] transform[r];
1925
}// end loop over 'r'
1926
delete [] transform;
1932
// Array of basisvalues
1933
double basisvalues[3] = {0.0, 0.0, 0.0};
1935
// Declare helper variables
1936
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1938
// Compute basisvalues
1939
basisvalues[0] = 1.0;
1940
basisvalues[1] = tmp0;
1941
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1942
basisvalues[0] *= std::sqrt(0.5);
1943
basisvalues[2] *= std::sqrt(1.0);
1944
basisvalues[1] *= std::sqrt(3.0);
1946
// Table(s) of coefficients
1947
static const double coefficients0[3] = \
1948
{0.47140452, -0.28867513, -0.16666667};
1950
// Tables of derivatives of the polynomial base (transpose).
1951
static const double dmats0[3][3] = \
1953
{4.8989795, 0.0, 0.0},
1956
static const double dmats1[3][3] = \
1958
{2.4494897, 0.0, 0.0},
1959
{4.2426407, 0.0, 0.0}};
1961
// Compute reference derivatives.
1962
// Declare pointer to array of derivatives on FIAT element.
1963
double *derivatives = new double[num_derivatives_t];
1964
for (unsigned int r = 0; r < num_derivatives_t; r++)
1966
derivatives[r] = 0.0;
1967
}// end loop over 'r'
1969
// Declare derivative matrix (of polynomial basis).
1970
double dmats[3][3] = \
1975
// Declare (auxiliary) derivative matrix (of polynomial basis).
1976
double dmats_old[3][3] = \
1981
// Loop possible derivatives.
1982
for (unsigned int r = 0; r < num_derivatives_t; r++)
1984
// Resetting dmats values to compute next derivative.
1985
for (unsigned int t = 0; t < 3; t++)
1987
for (unsigned int u = 0; u < 3; u++)
1995
}// end loop over 'u'
1996
}// end loop over 't'
1998
// Looping derivative order to generate dmats.
1999
for (unsigned int s = 0; s < n; s++)
2001
// Updating dmats_old with new values and resetting dmats.
2002
for (unsigned int t = 0; t < 3; t++)
2004
for (unsigned int u = 0; u < 3; u++)
2006
dmats_old[t][u] = dmats[t][u];
2008
}// end loop over 'u'
2009
}// end loop over 't'
2011
// Update dmats using an inner product.
2012
if (combinations_t[r][s] == 0)
2014
for (unsigned int t = 0; t < 3; t++)
2016
for (unsigned int u = 0; u < 3; u++)
2018
for (unsigned int tu = 0; tu < 3; tu++)
2020
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2021
}// end loop over 'tu'
2022
}// end loop over 'u'
2023
}// end loop over 't'
2026
if (combinations_t[r][s] == 1)
2028
for (unsigned int t = 0; t < 3; t++)
2030
for (unsigned int u = 0; u < 3; u++)
2032
for (unsigned int tu = 0; tu < 3; tu++)
2034
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2035
}// end loop over 'tu'
2036
}// end loop over 'u'
2037
}// end loop over 't'
2040
}// end loop over 's'
2041
for (unsigned int s = 0; s < 3; s++)
2043
for (unsigned int t = 0; t < 3; t++)
2045
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2046
}// end loop over 't'
2047
}// end loop over 's'
2048
}// end loop over 'r'
2050
// Transform derivatives back to physical element
2051
for (unsigned int r = 0; r < num_derivatives_g; r++)
2053
for (unsigned int s = 0; s < num_derivatives_t; s++)
2055
values[num_derivatives_g + r] += transform[r][s]*derivatives[s];
2056
}// end loop over 's'
2057
}// end loop over 'r'
2059
// Delete pointer to array of derivatives on FIAT element
2060
delete [] derivatives;
2062
// Delete pointer to array of combinations of derivatives and transform
2063
for (unsigned int r = 0; r < num_derivatives_t; r++)
2065
delete [] combinations_t[r];
2066
}// end loop over 'r'
2067
delete [] combinations_t;
2068
for (unsigned int r = 0; r < num_derivatives_g; r++)
2070
delete [] combinations_g[r];
2071
}// end loop over 'r'
2072
delete [] combinations_g;
2073
for (unsigned int r = 0; r < num_derivatives_g; r++)
2075
delete [] transform[r];
2076
}// end loop over 'r'
2077
delete [] transform;
2083
// Array of basisvalues
2084
double basisvalues[3] = {0.0, 0.0, 0.0};
2086
// Declare helper variables
2087
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2089
// Compute basisvalues
2090
basisvalues[0] = 1.0;
2091
basisvalues[1] = tmp0;
2092
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2093
basisvalues[0] *= std::sqrt(0.5);
2094
basisvalues[2] *= std::sqrt(1.0);
2095
basisvalues[1] *= std::sqrt(3.0);
2097
// Table(s) of coefficients
2098
static const double coefficients0[3] = \
2099
{0.47140452, 0.28867513, -0.16666667};
2101
// Tables of derivatives of the polynomial base (transpose).
2102
static const double dmats0[3][3] = \
2104
{4.8989795, 0.0, 0.0},
2107
static const double dmats1[3][3] = \
2109
{2.4494897, 0.0, 0.0},
2110
{4.2426407, 0.0, 0.0}};
2112
// Compute reference derivatives.
2113
// Declare pointer to array of derivatives on FIAT element.
2114
double *derivatives = new double[num_derivatives_t];
2115
for (unsigned int r = 0; r < num_derivatives_t; r++)
2117
derivatives[r] = 0.0;
2118
}// end loop over 'r'
2120
// Declare derivative matrix (of polynomial basis).
2121
double dmats[3][3] = \
2126
// Declare (auxiliary) derivative matrix (of polynomial basis).
2127
double dmats_old[3][3] = \
2132
// Loop possible derivatives.
2133
for (unsigned int r = 0; r < num_derivatives_t; r++)
2135
// Resetting dmats values to compute next derivative.
2136
for (unsigned int t = 0; t < 3; t++)
2138
for (unsigned int u = 0; u < 3; u++)
2146
}// end loop over 'u'
2147
}// end loop over 't'
2149
// Looping derivative order to generate dmats.
2150
for (unsigned int s = 0; s < n; s++)
2152
// Updating dmats_old with new values and resetting dmats.
2153
for (unsigned int t = 0; t < 3; t++)
2155
for (unsigned int u = 0; u < 3; u++)
2157
dmats_old[t][u] = dmats[t][u];
2159
}// end loop over 'u'
2160
}// end loop over 't'
2162
// Update dmats using an inner product.
2163
if (combinations_t[r][s] == 0)
2165
for (unsigned int t = 0; t < 3; t++)
2167
for (unsigned int u = 0; u < 3; u++)
2169
for (unsigned int tu = 0; tu < 3; tu++)
2171
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2172
}// end loop over 'tu'
2173
}// end loop over 'u'
2174
}// end loop over 't'
2177
if (combinations_t[r][s] == 1)
2179
for (unsigned int t = 0; t < 3; t++)
2181
for (unsigned int u = 0; u < 3; u++)
2183
for (unsigned int tu = 0; tu < 3; tu++)
2185
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2186
}// end loop over 'tu'
2187
}// end loop over 'u'
2188
}// end loop over 't'
2191
}// end loop over 's'
2192
for (unsigned int s = 0; s < 3; s++)
2194
for (unsigned int t = 0; t < 3; t++)
2196
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2197
}// end loop over 't'
2198
}// end loop over 's'
2199
}// end loop over 'r'
2201
// Transform derivatives back to physical element
2202
for (unsigned int r = 0; r < num_derivatives_g; r++)
2204
for (unsigned int s = 0; s < num_derivatives_t; s++)
2206
values[num_derivatives_g + r] += transform[r][s]*derivatives[s];
2207
}// end loop over 's'
2208
}// end loop over 'r'
2210
// Delete pointer to array of derivatives on FIAT element
2211
delete [] derivatives;
2213
// Delete pointer to array of combinations of derivatives and transform
2214
for (unsigned int r = 0; r < num_derivatives_t; r++)
2216
delete [] combinations_t[r];
2217
}// end loop over 'r'
2218
delete [] combinations_t;
2219
for (unsigned int r = 0; r < num_derivatives_g; r++)
2221
delete [] combinations_g[r];
2222
}// end loop over 'r'
2223
delete [] combinations_g;
2224
for (unsigned int r = 0; r < num_derivatives_g; r++)
2226
delete [] transform[r];
2227
}// end loop over 'r'
2228
delete [] transform;
2234
// Array of basisvalues
2235
double basisvalues[3] = {0.0, 0.0, 0.0};
2237
// Declare helper variables
2238
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2240
// Compute basisvalues
2241
basisvalues[0] = 1.0;
2242
basisvalues[1] = tmp0;
2243
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2244
basisvalues[0] *= std::sqrt(0.5);
2245
basisvalues[2] *= std::sqrt(1.0);
2246
basisvalues[1] *= std::sqrt(3.0);
2248
// Table(s) of coefficients
2249
static const double coefficients0[3] = \
2250
{0.47140452, 0.0, 0.33333333};
2252
// Tables of derivatives of the polynomial base (transpose).
2253
static const double dmats0[3][3] = \
2255
{4.8989795, 0.0, 0.0},
2258
static const double dmats1[3][3] = \
2260
{2.4494897, 0.0, 0.0},
2261
{4.2426407, 0.0, 0.0}};
2263
// Compute reference derivatives.
2264
// Declare pointer to array of derivatives on FIAT element.
2265
double *derivatives = new double[num_derivatives_t];
2266
for (unsigned int r = 0; r < num_derivatives_t; r++)
2268
derivatives[r] = 0.0;
2269
}// end loop over 'r'
2271
// Declare derivative matrix (of polynomial basis).
2272
double dmats[3][3] = \
2277
// Declare (auxiliary) derivative matrix (of polynomial basis).
2278
double dmats_old[3][3] = \
2283
// Loop possible derivatives.
2284
for (unsigned int r = 0; r < num_derivatives_t; r++)
2286
// Resetting dmats values to compute next derivative.
2287
for (unsigned int t = 0; t < 3; t++)
2289
for (unsigned int u = 0; u < 3; u++)
2297
}// end loop over 'u'
2298
}// end loop over 't'
2300
// Looping derivative order to generate dmats.
2301
for (unsigned int s = 0; s < n; s++)
2303
// Updating dmats_old with new values and resetting dmats.
2304
for (unsigned int t = 0; t < 3; t++)
2306
for (unsigned int u = 0; u < 3; u++)
2308
dmats_old[t][u] = dmats[t][u];
2310
}// end loop over 'u'
2311
}// end loop over 't'
2313
// Update dmats using an inner product.
2314
if (combinations_t[r][s] == 0)
2316
for (unsigned int t = 0; t < 3; t++)
2318
for (unsigned int u = 0; u < 3; u++)
2320
for (unsigned int tu = 0; tu < 3; tu++)
2322
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2323
}// end loop over 'tu'
2324
}// end loop over 'u'
2325
}// end loop over 't'
2328
if (combinations_t[r][s] == 1)
2330
for (unsigned int t = 0; t < 3; t++)
2332
for (unsigned int u = 0; u < 3; u++)
2334
for (unsigned int tu = 0; tu < 3; tu++)
2336
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2337
}// end loop over 'tu'
2338
}// end loop over 'u'
2339
}// end loop over 't'
2342
}// end loop over 's'
2343
for (unsigned int s = 0; s < 3; s++)
2345
for (unsigned int t = 0; t < 3; t++)
2347
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2348
}// end loop over 't'
2349
}// end loop over 's'
2350
}// end loop over 'r'
2352
// Transform derivatives back to physical element
2353
for (unsigned int r = 0; r < num_derivatives_g; r++)
2355
for (unsigned int s = 0; s < num_derivatives_t; s++)
2357
values[num_derivatives_g + r] += transform[r][s]*derivatives[s];
2358
}// end loop over 's'
2359
}// end loop over 'r'
2361
// Delete pointer to array of derivatives on FIAT element
2362
delete [] derivatives;
2364
// Delete pointer to array of combinations of derivatives and transform
2365
for (unsigned int r = 0; r < num_derivatives_t; r++)
2367
delete [] combinations_t[r];
2368
}// end loop over 'r'
2369
delete [] combinations_t;
2370
for (unsigned int r = 0; r < num_derivatives_g; r++)
2372
delete [] combinations_g[r];
2373
}// end loop over 'r'
2374
delete [] combinations_g;
2375
for (unsigned int r = 0; r < num_derivatives_g; r++)
2377
delete [] transform[r];
2378
}// end loop over 'r'
2379
delete [] transform;
2385
// Array of basisvalues
2386
double basisvalues[3] = {0.0, 0.0, 0.0};
2388
// Declare helper variables
2389
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2391
// Compute basisvalues
2392
basisvalues[0] = 1.0;
2393
basisvalues[1] = tmp0;
2394
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2395
basisvalues[0] *= std::sqrt(0.5);
2396
basisvalues[2] *= std::sqrt(1.0);
2397
basisvalues[1] *= std::sqrt(3.0);
2399
// Table(s) of coefficients
2400
static const double coefficients0[3] = \
2401
{0.47140452, -0.28867513, -0.16666667};
2403
// Tables of derivatives of the polynomial base (transpose).
2404
static const double dmats0[3][3] = \
2406
{4.8989795, 0.0, 0.0},
2409
static const double dmats1[3][3] = \
2411
{2.4494897, 0.0, 0.0},
2412
{4.2426407, 0.0, 0.0}};
2414
// Compute reference derivatives.
2415
// Declare pointer to array of derivatives on FIAT element.
2416
double *derivatives = new double[num_derivatives_t];
2417
for (unsigned int r = 0; r < num_derivatives_t; r++)
2419
derivatives[r] = 0.0;
2420
}// end loop over 'r'
2422
// Declare derivative matrix (of polynomial basis).
2423
double dmats[3][3] = \
2428
// Declare (auxiliary) derivative matrix (of polynomial basis).
2429
double dmats_old[3][3] = \
2434
// Loop possible derivatives.
2435
for (unsigned int r = 0; r < num_derivatives_t; r++)
2437
// Resetting dmats values to compute next derivative.
2438
for (unsigned int t = 0; t < 3; t++)
2440
for (unsigned int u = 0; u < 3; u++)
2448
}// end loop over 'u'
2449
}// end loop over 't'
2451
// Looping derivative order to generate dmats.
2452
for (unsigned int s = 0; s < n; s++)
2454
// Updating dmats_old with new values and resetting dmats.
2455
for (unsigned int t = 0; t < 3; t++)
2457
for (unsigned int u = 0; u < 3; u++)
2459
dmats_old[t][u] = dmats[t][u];
2461
}// end loop over 'u'
2462
}// end loop over 't'
2464
// Update dmats using an inner product.
2465
if (combinations_t[r][s] == 0)
2467
for (unsigned int t = 0; t < 3; t++)
2469
for (unsigned int u = 0; u < 3; u++)
2471
for (unsigned int tu = 0; tu < 3; tu++)
2473
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2474
}// end loop over 'tu'
2475
}// end loop over 'u'
2476
}// end loop over 't'
2479
if (combinations_t[r][s] == 1)
2481
for (unsigned int t = 0; t < 3; t++)
2483
for (unsigned int u = 0; u < 3; u++)
2485
for (unsigned int tu = 0; tu < 3; tu++)
2487
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2488
}// end loop over 'tu'
2489
}// end loop over 'u'
2490
}// end loop over 't'
2493
}// end loop over 's'
2494
for (unsigned int s = 0; s < 3; s++)
2496
for (unsigned int t = 0; t < 3; t++)
2498
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2499
}// end loop over 't'
2500
}// end loop over 's'
2501
}// end loop over 'r'
2503
// Transform derivatives back to physical element
2504
for (unsigned int r = 0; r < num_derivatives_g; r++)
2506
for (unsigned int s = 0; s < num_derivatives_t; s++)
2508
values[2*num_derivatives_g + r] += transform[r][s]*derivatives[s];
2509
}// end loop over 's'
2510
}// end loop over 'r'
2512
// Delete pointer to array of derivatives on FIAT element
2513
delete [] derivatives;
2515
// Delete pointer to array of combinations of derivatives and transform
2516
for (unsigned int r = 0; r < num_derivatives_t; r++)
2518
delete [] combinations_t[r];
2519
}// end loop over 'r'
2520
delete [] combinations_t;
2521
for (unsigned int r = 0; r < num_derivatives_g; r++)
2523
delete [] combinations_g[r];
2524
}// end loop over 'r'
2525
delete [] combinations_g;
2526
for (unsigned int r = 0; r < num_derivatives_g; r++)
2528
delete [] transform[r];
2529
}// end loop over 'r'
2530
delete [] transform;
2536
// Array of basisvalues
2537
double basisvalues[3] = {0.0, 0.0, 0.0};
2539
// Declare helper variables
2540
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2542
// Compute basisvalues
2543
basisvalues[0] = 1.0;
2544
basisvalues[1] = tmp0;
2545
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2546
basisvalues[0] *= std::sqrt(0.5);
2547
basisvalues[2] *= std::sqrt(1.0);
2548
basisvalues[1] *= std::sqrt(3.0);
2550
// Table(s) of coefficients
2551
static const double coefficients0[3] = \
2552
{0.47140452, 0.28867513, -0.16666667};
2554
// Tables of derivatives of the polynomial base (transpose).
2555
static const double dmats0[3][3] = \
2557
{4.8989795, 0.0, 0.0},
2560
static const double dmats1[3][3] = \
2562
{2.4494897, 0.0, 0.0},
2563
{4.2426407, 0.0, 0.0}};
2565
// Compute reference derivatives.
2566
// Declare pointer to array of derivatives on FIAT element.
2567
double *derivatives = new double[num_derivatives_t];
2568
for (unsigned int r = 0; r < num_derivatives_t; r++)
2570
derivatives[r] = 0.0;
2571
}// end loop over 'r'
2573
// Declare derivative matrix (of polynomial basis).
2574
double dmats[3][3] = \
2579
// Declare (auxiliary) derivative matrix (of polynomial basis).
2580
double dmats_old[3][3] = \
2585
// Loop possible derivatives.
2586
for (unsigned int r = 0; r < num_derivatives_t; r++)
2588
// Resetting dmats values to compute next derivative.
2589
for (unsigned int t = 0; t < 3; t++)
2591
for (unsigned int u = 0; u < 3; u++)
2599
}// end loop over 'u'
2600
}// end loop over 't'
2602
// Looping derivative order to generate dmats.
2603
for (unsigned int s = 0; s < n; s++)
2605
// Updating dmats_old with new values and resetting dmats.
2606
for (unsigned int t = 0; t < 3; t++)
2608
for (unsigned int u = 0; u < 3; u++)
2610
dmats_old[t][u] = dmats[t][u];
2612
}// end loop over 'u'
2613
}// end loop over 't'
2615
// Update dmats using an inner product.
2616
if (combinations_t[r][s] == 0)
2618
for (unsigned int t = 0; t < 3; t++)
2620
for (unsigned int u = 0; u < 3; u++)
2622
for (unsigned int tu = 0; tu < 3; tu++)
2624
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2625
}// end loop over 'tu'
2626
}// end loop over 'u'
2627
}// end loop over 't'
2630
if (combinations_t[r][s] == 1)
2632
for (unsigned int t = 0; t < 3; t++)
2634
for (unsigned int u = 0; u < 3; u++)
2636
for (unsigned int tu = 0; tu < 3; tu++)
2638
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2639
}// end loop over 'tu'
2640
}// end loop over 'u'
2641
}// end loop over 't'
2644
}// end loop over 's'
2645
for (unsigned int s = 0; s < 3; s++)
2647
for (unsigned int t = 0; t < 3; t++)
2649
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2650
}// end loop over 't'
2651
}// end loop over 's'
2652
}// end loop over 'r'
2654
// Transform derivatives back to physical element
2655
for (unsigned int r = 0; r < num_derivatives_g; r++)
2657
for (unsigned int s = 0; s < num_derivatives_t; s++)
2659
values[2*num_derivatives_g + r] += transform[r][s]*derivatives[s];
2660
}// end loop over 's'
2661
}// end loop over 'r'
2663
// Delete pointer to array of derivatives on FIAT element
2664
delete [] derivatives;
2666
// Delete pointer to array of combinations of derivatives and transform
2667
for (unsigned int r = 0; r < num_derivatives_t; r++)
2669
delete [] combinations_t[r];
2670
}// end loop over 'r'
2671
delete [] combinations_t;
2672
for (unsigned int r = 0; r < num_derivatives_g; r++)
2674
delete [] combinations_g[r];
2675
}// end loop over 'r'
2676
delete [] combinations_g;
2677
for (unsigned int r = 0; r < num_derivatives_g; r++)
2679
delete [] transform[r];
2680
}// end loop over 'r'
2681
delete [] transform;
2687
// Array of basisvalues
2688
double basisvalues[3] = {0.0, 0.0, 0.0};
2690
// Declare helper variables
2691
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2693
// Compute basisvalues
2694
basisvalues[0] = 1.0;
2695
basisvalues[1] = tmp0;
2696
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2697
basisvalues[0] *= std::sqrt(0.5);
2698
basisvalues[2] *= std::sqrt(1.0);
2699
basisvalues[1] *= std::sqrt(3.0);
2701
// Table(s) of coefficients
2702
static const double coefficients0[3] = \
2703
{0.47140452, 0.0, 0.33333333};
2705
// Tables of derivatives of the polynomial base (transpose).
2706
static const double dmats0[3][3] = \
2708
{4.8989795, 0.0, 0.0},
2711
static const double dmats1[3][3] = \
2713
{2.4494897, 0.0, 0.0},
2714
{4.2426407, 0.0, 0.0}};
2716
// Compute reference derivatives.
2717
// Declare pointer to array of derivatives on FIAT element.
2718
double *derivatives = new double[num_derivatives_t];
2719
for (unsigned int r = 0; r < num_derivatives_t; r++)
2721
derivatives[r] = 0.0;
2722
}// end loop over 'r'
2724
// Declare derivative matrix (of polynomial basis).
2725
double dmats[3][3] = \
2730
// Declare (auxiliary) derivative matrix (of polynomial basis).
2731
double dmats_old[3][3] = \
2736
// Loop possible derivatives.
2737
for (unsigned int r = 0; r < num_derivatives_t; r++)
2739
// Resetting dmats values to compute next derivative.
2740
for (unsigned int t = 0; t < 3; t++)
2742
for (unsigned int u = 0; u < 3; u++)
2750
}// end loop over 'u'
2751
}// end loop over 't'
2753
// Looping derivative order to generate dmats.
2754
for (unsigned int s = 0; s < n; s++)
2756
// Updating dmats_old with new values and resetting dmats.
2757
for (unsigned int t = 0; t < 3; t++)
2759
for (unsigned int u = 0; u < 3; u++)
2761
dmats_old[t][u] = dmats[t][u];
2763
}// end loop over 'u'
2764
}// end loop over 't'
2766
// Update dmats using an inner product.
2767
if (combinations_t[r][s] == 0)
2769
for (unsigned int t = 0; t < 3; t++)
2771
for (unsigned int u = 0; u < 3; u++)
2773
for (unsigned int tu = 0; tu < 3; tu++)
2775
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2776
}// end loop over 'tu'
2777
}// end loop over 'u'
2778
}// end loop over 't'
2781
if (combinations_t[r][s] == 1)
2783
for (unsigned int t = 0; t < 3; t++)
2785
for (unsigned int u = 0; u < 3; u++)
2787
for (unsigned int tu = 0; tu < 3; tu++)
2789
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2790
}// end loop over 'tu'
2791
}// end loop over 'u'
2792
}// end loop over 't'
2795
}// end loop over 's'
2796
for (unsigned int s = 0; s < 3; s++)
2798
for (unsigned int t = 0; t < 3; t++)
2800
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2801
}// end loop over 't'
2802
}// end loop over 's'
2803
}// end loop over 'r'
2805
// Transform derivatives back to physical element
2806
for (unsigned int r = 0; r < num_derivatives_g; r++)
2808
for (unsigned int s = 0; s < num_derivatives_t; s++)
2810
values[2*num_derivatives_g + r] += transform[r][s]*derivatives[s];
2811
}// end loop over 's'
2812
}// end loop over 'r'
2814
// Delete pointer to array of derivatives on FIAT element
2815
delete [] derivatives;
2817
// Delete pointer to array of combinations of derivatives and transform
2818
for (unsigned int r = 0; r < num_derivatives_t; r++)
2820
delete [] combinations_t[r];
2821
}// end loop over 'r'
2822
delete [] combinations_t;
2823
for (unsigned int r = 0; r < num_derivatives_g; r++)
2825
delete [] combinations_g[r];
2826
}// end loop over 'r'
2827
delete [] combinations_g;
2828
for (unsigned int r = 0; r < num_derivatives_g; r++)
2830
delete [] transform[r];
2831
}// end loop over 'r'
2832
delete [] transform;
2839
/// Evaluate order n derivatives of all basis functions at given point x in cell
2840
virtual void evaluate_basis_derivatives_all(std::size_t n,
2843
const double* vertex_coordinates,
2844
int cell_orientation) const
2846
// Compute number of derivatives.
2847
unsigned int num_derivatives_g = 1;
2848
for (unsigned int r = 0; r < n; r++)
2850
num_derivatives_g *= 3;
2851
}// end loop over 'r'
2853
// Helper variable to hold values of a single dof.
2854
double *dof_values = new double[3*num_derivatives_g];
2855
for (unsigned int r = 0; r < 3*num_derivatives_g; r++)
2857
dof_values[r] = 0.0;
2858
}// end loop over 'r'
2860
// Loop dofs and call evaluate_basis_derivatives.
2861
for (unsigned int r = 0; r < 9; r++)
2863
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2864
for (unsigned int s = 0; s < 3*num_derivatives_g; s++)
2866
values[r*3*num_derivatives_g + s] = dof_values[s];
2867
}// end loop over 's'
2868
}// end loop over 'r'
2871
delete [] dof_values;
2874
/// Evaluate linear functional for dof i on the function f
2875
virtual double evaluate_dof(std::size_t i,
2876
const ufc::function& f,
2877
const double* vertex_coordinates,
2878
int cell_orientation,
2879
const ufc::cell& c) const
2881
// Declare variables for result of evaluation
2884
// Declare variable for physical coordinates
2890
y[0] = vertex_coordinates[0];
2891
y[1] = vertex_coordinates[1];
2892
y[2] = vertex_coordinates[2];
2893
f.evaluate(vals, y, c);
2899
y[0] = vertex_coordinates[3];
2900
y[1] = vertex_coordinates[4];
2901
y[2] = vertex_coordinates[5];
2902
f.evaluate(vals, y, c);
2908
y[0] = vertex_coordinates[6];
2909
y[1] = vertex_coordinates[7];
2910
y[2] = vertex_coordinates[8];
2911
f.evaluate(vals, y, c);
2917
y[0] = vertex_coordinates[0];
2918
y[1] = vertex_coordinates[1];
2919
y[2] = vertex_coordinates[2];
2920
f.evaluate(vals, y, c);
2926
y[0] = vertex_coordinates[3];
2927
y[1] = vertex_coordinates[4];
2928
y[2] = vertex_coordinates[5];
2929
f.evaluate(vals, y, c);
2935
y[0] = vertex_coordinates[6];
2936
y[1] = vertex_coordinates[7];
2937
y[2] = vertex_coordinates[8];
2938
f.evaluate(vals, y, c);
2944
y[0] = vertex_coordinates[0];
2945
y[1] = vertex_coordinates[1];
2946
y[2] = vertex_coordinates[2];
2947
f.evaluate(vals, y, c);
2953
y[0] = vertex_coordinates[3];
2954
y[1] = vertex_coordinates[4];
2955
y[2] = vertex_coordinates[5];
2956
f.evaluate(vals, y, c);
2962
y[0] = vertex_coordinates[6];
2963
y[1] = vertex_coordinates[7];
2964
y[2] = vertex_coordinates[8];
2965
f.evaluate(vals, y, c);
2974
/// Evaluate linear functionals for all dofs on the function f
2975
virtual void evaluate_dofs(double* values,
2976
const ufc::function& f,
2977
const double* vertex_coordinates,
2978
int cell_orientation,
2979
const ufc::cell& c) const
2981
// Declare variables for result of evaluation
2984
// Declare variable for physical coordinates
2986
y[0] = vertex_coordinates[0];
2987
y[1] = vertex_coordinates[1];
2988
y[2] = vertex_coordinates[2];
2989
f.evaluate(vals, y, c);
2990
values[0] = vals[0];
2991
y[0] = vertex_coordinates[3];
2992
y[1] = vertex_coordinates[4];
2993
y[2] = vertex_coordinates[5];
2994
f.evaluate(vals, y, c);
2995
values[1] = vals[0];
2996
y[0] = vertex_coordinates[6];
2997
y[1] = vertex_coordinates[7];
2998
y[2] = vertex_coordinates[8];
2999
f.evaluate(vals, y, c);
3000
values[2] = vals[0];
3001
y[0] = vertex_coordinates[0];
3002
y[1] = vertex_coordinates[1];
3003
y[2] = vertex_coordinates[2];
3004
f.evaluate(vals, y, c);
3005
values[3] = vals[1];
3006
y[0] = vertex_coordinates[3];
3007
y[1] = vertex_coordinates[4];
3008
y[2] = vertex_coordinates[5];
3009
f.evaluate(vals, y, c);
3010
values[4] = vals[1];
3011
y[0] = vertex_coordinates[6];
3012
y[1] = vertex_coordinates[7];
3013
y[2] = vertex_coordinates[8];
3014
f.evaluate(vals, y, c);
3015
values[5] = vals[1];
3016
y[0] = vertex_coordinates[0];
3017
y[1] = vertex_coordinates[1];
3018
y[2] = vertex_coordinates[2];
3019
f.evaluate(vals, y, c);
3020
values[6] = vals[2];
3021
y[0] = vertex_coordinates[3];
3022
y[1] = vertex_coordinates[4];
3023
y[2] = vertex_coordinates[5];
3024
f.evaluate(vals, y, c);
3025
values[7] = vals[2];
3026
y[0] = vertex_coordinates[6];
3027
y[1] = vertex_coordinates[7];
3028
y[2] = vertex_coordinates[8];
3029
f.evaluate(vals, y, c);
3030
values[8] = vals[2];
3033
/// Interpolate vertex values from dof values
3034
virtual void interpolate_vertex_values(double* vertex_values,
3035
const double* dof_values,
3036
const double* vertex_coordinates,
3037
int cell_orientation,
3038
const ufc::cell& c) const
3040
// Evaluate function and change variables
3041
vertex_values[0] = dof_values[0];
3042
vertex_values[3] = dof_values[1];
3043
vertex_values[6] = dof_values[2];
3044
// Evaluate function and change variables
3045
vertex_values[1] = dof_values[3];
3046
vertex_values[4] = dof_values[4];
3047
vertex_values[7] = dof_values[5];
3048
// Evaluate function and change variables
3049
vertex_values[2] = dof_values[6];
3050
vertex_values[5] = dof_values[7];
3051
vertex_values[8] = dof_values[8];
3054
/// Map coordinate xhat from reference cell to coordinate x in cell
3055
virtual void map_from_reference_cell(double* x,
3057
const ufc::cell& c) const
3059
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
3062
/// Map from coordinate x in cell to coordinate xhat in reference cell
3063
virtual void map_to_reference_cell(double* xhat,
3065
const ufc::cell& c) const
3067
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
3070
/// Return the number of sub elements (for a mixed element)
3071
virtual std::size_t num_sub_elements() const
3076
/// Create a new finite element for sub element i (for a mixed element)
3077
virtual ufc::finite_element* create_sub_element(std::size_t i) const
3083
return new x_element5_finite_element_0();
3088
return new x_element5_finite_element_0();
3093
return new x_element5_finite_element_0();
3101
/// Create a new class instance
3102
virtual ufc::finite_element* create() const
3104
return new x_element5_finite_element_1();
3109
/// This class defines the interface for a local-to-global mapping of
3110
/// degrees of freedom (dofs).
3112
class x_element5_dofmap_0: public ufc::dofmap
3117
x_element5_dofmap_0() : ufc::dofmap()
3123
virtual ~x_element5_dofmap_0()
3128
/// Return a string identifying the dofmap
3129
virtual const char* signature() const
3131
return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None)";
3134
/// Return true iff mesh entities of topological dimension d are needed
3135
virtual bool needs_mesh_entities(std::size_t d) const
3159
/// Return the topological dimension of the associated cell shape
3160
virtual std::size_t topological_dimension() const
3165
/// Return the geometric dimension of the associated cell shape
3166
virtual std::size_t geometric_dimension() const
3171
/// Return the dimension of the global finite element function space
3172
virtual std::size_t global_dimension(const std::vector<std::size_t>&
3173
num_global_entities) const
3175
return 3*num_global_entities[2];
3178
/// Return the dimension of the local finite element function space for a cell
3179
virtual std::size_t local_dimension(const ufc::cell& c) const
3184
/// Return the maximum dimension of the local finite element function space
3185
virtual std::size_t max_local_dimension() const
3190
/// Return the number of dofs on each cell facet
3191
virtual std::size_t num_facet_dofs() const
3196
/// Return the number of dofs associated with each cell entity of dimension d
3197
virtual std::size_t num_entity_dofs(std::size_t d) const
3221
/// Tabulate the local-to-global mapping of dofs on a cell
3222
virtual void tabulate_dofs(std::size_t* dofs,
3223
const std::vector<std::size_t>& num_global_entities,
3224
const ufc::cell& c) const
3226
dofs[0] = 3*c.entity_indices[2][0];
3227
dofs[1] = 3*c.entity_indices[2][0] + 1;
3228
dofs[2] = 3*c.entity_indices[2][0] + 2;
3231
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3232
virtual void tabulate_facet_dofs(std::size_t* dofs,
3233
std::size_t facet) const
3256
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3257
virtual void tabulate_entity_dofs(std::size_t* dofs,
3258
std::size_t d, std::size_t i) const
3262
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3281
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
3293
/// Tabulate the coordinates of all dofs on a cell
3294
virtual void tabulate_coordinates(double** dof_coordinates,
3295
const double* vertex_coordinates) const
3297
dof_coordinates[0][0] = vertex_coordinates[0];
3298
dof_coordinates[0][1] = vertex_coordinates[1];
3299
dof_coordinates[0][2] = vertex_coordinates[2];
3300
dof_coordinates[1][0] = vertex_coordinates[3];
3301
dof_coordinates[1][1] = vertex_coordinates[4];
3302
dof_coordinates[1][2] = vertex_coordinates[5];
3303
dof_coordinates[2][0] = vertex_coordinates[6];
3304
dof_coordinates[2][1] = vertex_coordinates[7];
3305
dof_coordinates[2][2] = vertex_coordinates[8];
3308
/// Return the number of sub dofmaps (for a mixed element)
3309
virtual std::size_t num_sub_dofmaps() const
3314
/// Create a new dofmap for sub dofmap i (for a mixed element)
3315
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3320
/// Create a new class instance
3321
virtual ufc::dofmap* create() const
3323
return new x_element5_dofmap_0();
3328
/// This class defines the interface for a local-to-global mapping of
3329
/// degrees of freedom (dofs).
3331
class x_element5_dofmap_1: public ufc::dofmap
3336
x_element5_dofmap_1() : ufc::dofmap()
3342
virtual ~x_element5_dofmap_1()
3347
/// Return a string identifying the dofmap
3348
virtual const char* signature() const
3350
return "FFC dofmap for VectorElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, 3, None)";
3353
/// Return true iff mesh entities of topological dimension d are needed
3354
virtual bool needs_mesh_entities(std::size_t d) const
3378
/// Return the topological dimension of the associated cell shape
3379
virtual std::size_t topological_dimension() const
3384
/// Return the geometric dimension of the associated cell shape
3385
virtual std::size_t geometric_dimension() const
3390
/// Return the dimension of the global finite element function space
3391
virtual std::size_t global_dimension(const std::vector<std::size_t>&
3392
num_global_entities) const
3394
return 9*num_global_entities[2];
3397
/// Return the dimension of the local finite element function space for a cell
3398
virtual std::size_t local_dimension(const ufc::cell& c) const
3403
/// Return the maximum dimension of the local finite element function space
3404
virtual std::size_t max_local_dimension() const
3409
/// Return the number of dofs on each cell facet
3410
virtual std::size_t num_facet_dofs() const
3415
/// Return the number of dofs associated with each cell entity of dimension d
3416
virtual std::size_t num_entity_dofs(std::size_t d) const
3440
/// Tabulate the local-to-global mapping of dofs on a cell
3441
virtual void tabulate_dofs(std::size_t* dofs,
3442
const std::vector<std::size_t>& num_global_entities,
3443
const ufc::cell& c) const
3445
unsigned int offset = 0;
3446
dofs[0] = offset + 3*c.entity_indices[2][0];
3447
dofs[1] = offset + 3*c.entity_indices[2][0] + 1;
3448
dofs[2] = offset + 3*c.entity_indices[2][0] + 2;
3449
offset += 3*num_global_entities[2];
3450
dofs[3] = offset + 3*c.entity_indices[2][0];
3451
dofs[4] = offset + 3*c.entity_indices[2][0] + 1;
3452
dofs[5] = offset + 3*c.entity_indices[2][0] + 2;
3453
offset += 3*num_global_entities[2];
3454
dofs[6] = offset + 3*c.entity_indices[2][0];
3455
dofs[7] = offset + 3*c.entity_indices[2][0] + 1;
3456
dofs[8] = offset + 3*c.entity_indices[2][0] + 2;
3457
offset += 3*num_global_entities[2];
3460
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3461
virtual void tabulate_facet_dofs(std::size_t* dofs,
3462
std::size_t facet) const
3485
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3486
virtual void tabulate_entity_dofs(std::size_t* dofs,
3487
std::size_t d, std::size_t i) const
3491
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3510
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
3528
/// Tabulate the coordinates of all dofs on a cell
3529
virtual void tabulate_coordinates(double** dof_coordinates,
3530
const double* vertex_coordinates) const
3532
dof_coordinates[0][0] = vertex_coordinates[0];
3533
dof_coordinates[0][1] = vertex_coordinates[1];
3534
dof_coordinates[0][2] = vertex_coordinates[2];
3535
dof_coordinates[1][0] = vertex_coordinates[3];
3536
dof_coordinates[1][1] = vertex_coordinates[4];
3537
dof_coordinates[1][2] = vertex_coordinates[5];
3538
dof_coordinates[2][0] = vertex_coordinates[6];
3539
dof_coordinates[2][1] = vertex_coordinates[7];
3540
dof_coordinates[2][2] = vertex_coordinates[8];
3541
dof_coordinates[3][0] = vertex_coordinates[0];
3542
dof_coordinates[3][1] = vertex_coordinates[1];
3543
dof_coordinates[3][2] = vertex_coordinates[2];
3544
dof_coordinates[4][0] = vertex_coordinates[3];
3545
dof_coordinates[4][1] = vertex_coordinates[4];
3546
dof_coordinates[4][2] = vertex_coordinates[5];
3547
dof_coordinates[5][0] = vertex_coordinates[6];
3548
dof_coordinates[5][1] = vertex_coordinates[7];
3549
dof_coordinates[5][2] = vertex_coordinates[8];
3550
dof_coordinates[6][0] = vertex_coordinates[0];
3551
dof_coordinates[6][1] = vertex_coordinates[1];
3552
dof_coordinates[6][2] = vertex_coordinates[2];
3553
dof_coordinates[7][0] = vertex_coordinates[3];
3554
dof_coordinates[7][1] = vertex_coordinates[4];
3555
dof_coordinates[7][2] = vertex_coordinates[5];
3556
dof_coordinates[8][0] = vertex_coordinates[6];
3557
dof_coordinates[8][1] = vertex_coordinates[7];
3558
dof_coordinates[8][2] = vertex_coordinates[8];
3561
/// Return the number of sub dofmaps (for a mixed element)
3562
virtual std::size_t num_sub_dofmaps() const
3567
/// Create a new dofmap for sub dofmap i (for a mixed element)
3568
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3574
return new x_element5_dofmap_0();
3579
return new x_element5_dofmap_0();
3584
return new x_element5_dofmap_0();
3592
/// Create a new class instance
3593
virtual ufc::dofmap* create() const
3595
return new x_element5_dofmap_1();