~ubuntu-branches/ubuntu/utopic/ffc/utopic
» Revision
17
: test/regression/references/r_auto/StabilisedStokes.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_auto/StabilisedStokes.h
- browse files at revision 17
- compare with another revision
- download diff
- download tarball
- view history from revision 17
- Committer: Package Import Robot
- Author(s): Johannes Ring
- Date: 2013-06-26 14:48:32 UTC
- mfrom: (1.1.13)
- Revision ID: package-import@ubuntu.com-20130626144832-1xd8htax4s3utybz
Tags: 1.2.0-1
* New upstream release.
* debian/control:
- Bump required version for python-ufc, python-fiat, python-instant
and python-ufl in Depends field.
- Bump Standards-Version to 3.9.4.
- Remove DM-Upload-Allowed field.
- Bump required debhelper version in Build-Depends.
- Remove cdbs from Build-Depends.
- Use canonical URIs for Vcs-* fields.
* debian/compat: Bump to compatibility level 9.
* debian/rules: Rewrite for debhelper (drop cdbs).
* debian/control:
- Bump required version for python-ufc, python-fiat, python-instant
and python-ufl in Depends field.
- Bump Standards-Version to 3.9.4.
- Remove DM-Upload-Allowed field.
- Bump required debhelper version in Build-Depends.
- Remove cdbs from Build-Depends.
- Use canonical URIs for Vcs-* fields.
* debian/compat: Bump to compatibility level 9.
* debian/rules: Rewrite for debhelper (drop cdbs).
- files added:
- .pc
- ffc/uflacsrepr
- sandbox/features
- sandbox/multidomain
- files modified:
- ChangeLog
added
removed
test/regression/references/r_auto/StabilisedStokes.h
1
// This code conforms with the UFC specification version 2.0.5
2
// and was automatically generated by FFC version 1.0.0.
1
// This code conforms with the UFC specification version 2.2.0
2
// and was automatically generated by FFC version 1.2.0.
3
3
//
4
4
// This code was generated with the following parameters:
5
5
//
52
52
/// Return a string identifying the finite element
53
53
virtual const char* signature() const
54
54
{
55
return "FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)";
55
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
56
56
}
57
57
58
58
/// Return the cell shape
62
62
}
63
63
64
64
/// Return the topological dimension of the cell shape
65
virtual unsigned int topological_dimension() const
65
virtual std::size_t topological_dimension() const
66
66
{
67
67
return 2;
68
68
}
69
69
70
70
/// Return the geometric dimension of the cell shape
71
virtual unsigned int geometric_dimension() const
71
virtual std::size_t geometric_dimension() const
72
72
{
73
73
return 2;
74
74
}
75
75
76
76
/// Return the dimension of the finite element function space
77
virtual unsigned int space_dimension() const
77
virtual std::size_t space_dimension() const
78
78
{
79
79
return 3;
80
80
}
81
81
82
82
/// Return the rank of the value space
83
virtual unsigned int value_rank() const
83
virtual std::size_t value_rank() const
84
84
{
85
85
return 0;
86
86
}
87
87
88
88
/// Return the dimension of the value space for axis i
89
virtual unsigned int value_dimension(unsigned int i) const
89
virtual std::size_t value_dimension(std::size_t i) const
90
90
{
91
91
return 1;
92
92
}
93
93
94
/// Evaluate basis function i at given point in cell
95
virtual void evaluate_basis(unsigned int i,
94
/// Evaluate basis function i at given point x in cell
95
virtual void evaluate_basis(std::size_t i,
96
96
double* values,
97
const double* coordinates,
98
const ufc::cell& c) const
97
const double* x,
98
const double* vertex_coordinates,
99
int cell_orientation) const
99
100
{
100
// Extract vertex coordinates
101
const double * const * x = c.coordinates;
102
103
// Compute Jacobian of affine map from reference cell
104
const double J_00 = x[1][0] - x[0][0];
105
const double J_01 = x[2][0] - x[0][0];
106
const double J_10 = x[1][1] - x[0][1];
107
const double J_11 = x[2][1] - x[0][1];
108
109
// Compute determinant of Jacobian
110
double detJ = J_00*J_11 - J_01*J_10;
111
112
// Compute inverse of Jacobian
101
// Compute Jacobian
102
double J[4];
103
compute_jacobian_triangle_2d(J, vertex_coordinates);
104
105
// Compute Jacobian inverse and determinant
106
double K[4];
107
double detJ;
108
compute_jacobian_inverse_triangle_2d(K, detJ, J);
109
113
110
114
111
// Compute constants
115
const double C0 = x[1][0] + x[2][0];
116
const double C1 = x[1][1] + x[2][1];
112
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
113
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
117
114
118
115
// Get coordinates and map to the reference (FIAT) element
119
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
120
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
116
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
117
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
121
118
122
// Reset values.
119
// Reset values
123
120
*values = 0.0;
124
121
switch (i)
125
122
{
126
123
case 0:
127
124
{
128
125
129
// Array of basisvalues.
126
// Array of basisvalues
130
127
double basisvalues[3] = {0.0, 0.0, 0.0};
131
128
132
// Declare helper variables.
129
// Declare helper variables
133
130
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
134
131
135
// Compute basisvalues.
132
// Compute basisvalues
136
133
basisvalues[0] = 1.0;
137
134
basisvalues[1] = tmp0;
138
135
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
140
137
basisvalues[2] *= std::sqrt(1.0);
141
138
basisvalues[1] *= std::sqrt(3.0);
142
139
143
// Table(s) of coefficients.
140
// Table(s) of coefficients
144
141
static const double coefficients0[3] = \
145
142
{0.47140452, -0.28867513, -0.16666667};
146
143
147
// Compute value(s).
144
// Compute value(s)
148
145
for (unsigned int r = 0; r < 3; r++)
149
146
{
150
147
*values += coefficients0[r]*basisvalues[r];
154
151
case 1:
155
152
{
156
153
157
// Array of basisvalues.
154
// Array of basisvalues
158
155
double basisvalues[3] = {0.0, 0.0, 0.0};
159
156
160
// Declare helper variables.
157
// Declare helper variables
161
158
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
162
159
163
// Compute basisvalues.
160
// Compute basisvalues
164
161
basisvalues[0] = 1.0;
165
162
basisvalues[1] = tmp0;
166
163
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
168
165
basisvalues[2] *= std::sqrt(1.0);
169
166
basisvalues[1] *= std::sqrt(3.0);
170
167
171
// Table(s) of coefficients.
168
// Table(s) of coefficients
172
169
static const double coefficients0[3] = \
173
170
{0.47140452, 0.28867513, -0.16666667};
174
171
175
// Compute value(s).
172
// Compute value(s)
176
173
for (unsigned int r = 0; r < 3; r++)
177
174
{
178
175
*values += coefficients0[r]*basisvalues[r];
182
179
case 2:
183
180
{
184
181
185
// Array of basisvalues.
182
// Array of basisvalues
186
183
double basisvalues[3] = {0.0, 0.0, 0.0};
187
184
188
// Declare helper variables.
185
// Declare helper variables
189
186
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
190
187
191
// Compute basisvalues.
188
// Compute basisvalues
192
189
basisvalues[0] = 1.0;
193
190
basisvalues[1] = tmp0;
194
191
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
196
193
basisvalues[2] *= std::sqrt(1.0);
197
194
basisvalues[1] *= std::sqrt(3.0);
198
195
199
// Table(s) of coefficients.
196
// Table(s) of coefficients
200
197
static const double coefficients0[3] = \
201
198
{0.47140452, 0.0, 0.33333333};
202
199
203
// Compute value(s).
200
// Compute value(s)
204
201
for (unsigned int r = 0; r < 3; r++)
205
202
{
206
203
*values += coefficients0[r]*basisvalues[r];
211
208
212
209
}
213
210
214
/// Evaluate all basis functions at given point in cell
211
/// Evaluate all basis functions at given point x in cell
215
212
virtual void evaluate_basis_all(double* values,
216
const double* coordinates,
217
const ufc::cell& c) const
213
const double* x,
214
const double* vertex_coordinates,
215
int cell_orientation) const
218
216
{
219
217
// Helper variable to hold values of a single dof.
220
218
double dof_values = 0.0;
221
219
222
// Loop dofs and call evaluate_basis.
220
// Loop dofs and call evaluate_basis
223
221
for (unsigned int r = 0; r < 3; r++)
224
222
{
225
evaluate_basis(r, &dof_values, coordinates, c);
223
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
226
224
values[r] = dof_values;
227
225
}// end loop over 'r'
228
226
}
229
227
230
/// Evaluate order n derivatives of basis function i at given point in cell
231
virtual void evaluate_basis_derivatives(unsigned int i,
232
unsigned int n,
228
/// Evaluate order n derivatives of basis function i at given point x in cell
229
virtual void evaluate_basis_derivatives(std::size_t i,
230
std::size_t n,
233
231
double* values,
234
const double* coordinates,
235
const ufc::cell& c) const
232
const double* x,
233
const double* vertex_coordinates,
234
int cell_orientation) const
236
235
{
237
// Extract vertex coordinates
238
const double * const * x = c.coordinates;
239
240
// Compute Jacobian of affine map from reference cell
241
const double J_00 = x[1][0] - x[0][0];
242
const double J_01 = x[2][0] - x[0][0];
243
const double J_10 = x[1][1] - x[0][1];
244
const double J_11 = x[2][1] - x[0][1];
245
246
// Compute determinant of Jacobian
247
double detJ = J_00*J_11 - J_01*J_10;
248
249
// Compute inverse of Jacobian
250
const double K_00 = J_11 / detJ;
251
const double K_01 = -J_01 / detJ;
252
const double K_10 = -J_10 / detJ;
253
const double K_11 = J_00 / detJ;
236
// Compute Jacobian
237
double J[4];
238
compute_jacobian_triangle_2d(J, vertex_coordinates);
239
240
// Compute Jacobian inverse and determinant
241
double K[4];
242
double detJ;
243
compute_jacobian_inverse_triangle_2d(K, detJ, J);
244
254
245
255
246
// Compute constants
256
const double C0 = x[1][0] + x[2][0];
257
const double C1 = x[1][1] + x[2][1];
247
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
248
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
258
249
259
250
// Get coordinates and map to the reference (FIAT) element
260
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
261
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
251
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
252
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
262
253
263
254
// Compute number of derivatives.
264
255
unsigned int num_derivatives = 1;
295
286
}
296
287
297
288
// Compute inverse of Jacobian
298
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
289
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
299
290
300
291
// Declare transformation matrix
301
292
// Declare pointer to two dimensional array and initialise
329
320
case 0:
330
321
{
331
322
332
// Array of basisvalues.
323
// Array of basisvalues
333
324
double basisvalues[3] = {0.0, 0.0, 0.0};
334
325
335
// Declare helper variables.
326
// Declare helper variables
336
327
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
337
328
338
// Compute basisvalues.
329
// Compute basisvalues
339
330
basisvalues[0] = 1.0;
340
331
basisvalues[1] = tmp0;
341
332
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
343
334
basisvalues[2] *= std::sqrt(1.0);
344
335
basisvalues[1] *= std::sqrt(3.0);
345
336
346
// Table(s) of coefficients.
337
// Table(s) of coefficients
347
338
static const double coefficients0[3] = \
348
339
{0.47140452, -0.28867513, -0.16666667};
349
340
475
466
case 1:
476
467
{
477
468
478
// Array of basisvalues.
469
// Array of basisvalues
479
470
double basisvalues[3] = {0.0, 0.0, 0.0};
480
471
481
// Declare helper variables.
472
// Declare helper variables
482
473
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
483
474
484
// Compute basisvalues.
475
// Compute basisvalues
485
476
basisvalues[0] = 1.0;
486
477
basisvalues[1] = tmp0;
487
478
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
489
480
basisvalues[2] *= std::sqrt(1.0);
490
481
basisvalues[1] *= std::sqrt(3.0);
491
482
492
// Table(s) of coefficients.
483
// Table(s) of coefficients
493
484
static const double coefficients0[3] = \
494
485
{0.47140452, 0.28867513, -0.16666667};
495
486
621
612
case 2:
622
613
{
623
614
624
// Array of basisvalues.
615
// Array of basisvalues
625
616
double basisvalues[3] = {0.0, 0.0, 0.0};
626
617
627
// Declare helper variables.
618
// Declare helper variables
628
619
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
629
620
630
// Compute basisvalues.
621
// Compute basisvalues
631
622
basisvalues[0] = 1.0;
632
623
basisvalues[1] = tmp0;
633
624
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
635
626
basisvalues[2] *= std::sqrt(1.0);
636
627
basisvalues[1] *= std::sqrt(3.0);
637
628
638
// Table(s) of coefficients.
629
// Table(s) of coefficients
639
630
static const double coefficients0[3] = \
640
631
{0.47140452, 0.0, 0.33333333};
641
632
768
759
769
760
}
770
761
771
/// Evaluate order n derivatives of all basis functions at given point in cell
772
virtual void evaluate_basis_derivatives_all(unsigned int n,
762
/// Evaluate order n derivatives of all basis functions at given point x in cell
763
virtual void evaluate_basis_derivatives_all(std::size_t n,
773
764
double* values,
774
const double* coordinates,
775
const ufc::cell& c) const
765
const double* x,
766
const double* vertex_coordinates,
767
int cell_orientation) const
776
768
{
777
769
// Compute number of derivatives.
778
770
unsigned int num_derivatives = 1;
791
783
// Loop dofs and call evaluate_basis_derivatives.
792
784
for (unsigned int r = 0; r < 3; r++)
793
785
{
794
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
786
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
795
787
for (unsigned int s = 0; s < num_derivatives; s++)
796
788
{
797
789
values[r*num_derivatives + s] = dof_values[s];
803
795
}
804
796
805
797
/// Evaluate linear functional for dof i on the function f
806
virtual double evaluate_dof(unsigned int i,
798
virtual double evaluate_dof(std::size_t i,
807
799
const ufc::function& f,
800
const double* vertex_coordinates,
801
int cell_orientation,
808
802
const ufc::cell& c) const
809
803
{
810
// Declare variables for result of evaluation.
804
// Declare variables for result of evaluation
811
805
double vals[1];
812
806
813
// Declare variable for physical coordinates.
807
// Declare variable for physical coordinates
814
808
double y[2];
815
const double * const * x = c.coordinates;
816
809
switch (i)
817
810
{
818
811
case 0:
819
812
{
820
y[0] = x[0][0];
821
y[1] = x[0][1];
813
y[0] = vertex_coordinates[0];
814
y[1] = vertex_coordinates[1];
822
815
f.evaluate(vals, y, c);
823
816
return vals[0];
824
817
break;
825
818
}
826
819
case 1:
827
820
{
828
y[0] = x[1][0];
829
y[1] = x[1][1];
821
y[0] = vertex_coordinates[2];
822
y[1] = vertex_coordinates[3];
830
823
f.evaluate(vals, y, c);
831
824
return vals[0];
832
825
break;
833
826
}
834
827
case 2:
835
828
{
836
y[0] = x[2][0];
837
y[1] = x[2][1];
829
y[0] = vertex_coordinates[4];
830
y[1] = vertex_coordinates[5];
838
831
f.evaluate(vals, y, c);
839
832
return vals[0];
840
833
break;
847
840
/// Evaluate linear functionals for all dofs on the function f
848
841
virtual void evaluate_dofs(double* values,
849
842
const ufc::function& f,
843
const double* vertex_coordinates,
844
int cell_orientation,
850
845
const ufc::cell& c) const
851
846
{
852
// Declare variables for result of evaluation.
847
// Declare variables for result of evaluation
853
848
double vals[1];
854
849
855
// Declare variable for physical coordinates.
850
// Declare variable for physical coordinates
856
851
double y[2];
857
const double * const * x = c.coordinates;
858
y[0] = x[0][0];
859
y[1] = x[0][1];
852
y[0] = vertex_coordinates[0];
853
y[1] = vertex_coordinates[1];
860
854
f.evaluate(vals, y, c);
861
855
values[0] = vals[0];
862
y[0] = x[1][0];
863
y[1] = x[1][1];
856
y[0] = vertex_coordinates[2];
857
y[1] = vertex_coordinates[3];
864
858
f.evaluate(vals, y, c);
865
859
values[1] = vals[0];
866
y[0] = x[2][0];
867
y[1] = x[2][1];
860
y[0] = vertex_coordinates[4];
861
y[1] = vertex_coordinates[5];
868
862
f.evaluate(vals, y, c);
869
863
values[2] = vals[0];
870
864
}
872
866
/// Interpolate vertex values from dof values
873
867
virtual void interpolate_vertex_values(double* vertex_values,
874
868
const double* dof_values,
869
const double* vertex_coordinates,
870
int cell_orientation,
875
871
const ufc::cell& c) const
876
872
{
877
873
// Evaluate function and change variables
885
881
const double* xhat,
886
882
const ufc::cell& c) const
887
883
{
888
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
884
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
889
885
}
890
886
891
887
/// Map from coordinate x in cell to coordinate xhat in reference cell
893
889
const double* x,
894
890
const ufc::cell& c) const
895
891
{
896
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
892
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
897
893
}
898
894
899
895
/// Return the number of sub elements (for a mixed element)
900
virtual unsigned int num_sub_elements() const
896
virtual std::size_t num_sub_elements() const
901
897
{
902
898
return 0;
903
899
}
904
900
905
901
/// Create a new finite element for sub element i (for a mixed element)
906
virtual ufc::finite_element* create_sub_element(unsigned int i) const
902
virtual ufc::finite_element* create_sub_element(std::size_t i) const
907
903
{
908
904
return 0;
909
905
}
937
933
/// Return a string identifying the finite element
938
934
virtual const char* signature() const
939
935
{
940
return "VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None)";
936
return "VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, 2, None)";
941
937
}
942
938
943
939
/// Return the cell shape
947
943
}
948
944
949
945
/// Return the topological dimension of the cell shape
950
virtual unsigned int topological_dimension() const
946
virtual std::size_t topological_dimension() const
951
947
{
952
948
return 2;
953
949
}
954
950
955
951
/// Return the geometric dimension of the cell shape
956
virtual unsigned int geometric_dimension() const
952
virtual std::size_t geometric_dimension() const
957
953
{
958
954
return 2;
959
955
}
960
956
961
957
/// Return the dimension of the finite element function space
962
virtual unsigned int space_dimension() const
958
virtual std::size_t space_dimension() const
963
959
{
964
960
return 6;
965
961
}
966
962
967
963
/// Return the rank of the value space
968
virtual unsigned int value_rank() const
964
virtual std::size_t value_rank() const
969
965
{
970
966
return 1;
971
967
}
972
968
973
969
/// Return the dimension of the value space for axis i
974
virtual unsigned int value_dimension(unsigned int i) const
970
virtual std::size_t value_dimension(std::size_t i) const
975
971
{
976
972
switch (i)
977
973
{
985
981
return 0;
986
982
}
987
983
988
/// Evaluate basis function i at given point in cell
989
virtual void evaluate_basis(unsigned int i,
984
/// Evaluate basis function i at given point x in cell
985
virtual void evaluate_basis(std::size_t i,
990
986
double* values,
991
const double* coordinates,
992
const ufc::cell& c) const
987
const double* x,
988
const double* vertex_coordinates,
989
int cell_orientation) const
993
990
{
994
// Extract vertex coordinates
995
const double * const * x = c.coordinates;
996
997
// Compute Jacobian of affine map from reference cell
998
const double J_00 = x[1][0] - x[0][0];
999
const double J_01 = x[2][0] - x[0][0];
1000
const double J_10 = x[1][1] - x[0][1];
1001
const double J_11 = x[2][1] - x[0][1];
1002
1003
// Compute determinant of Jacobian
1004
double detJ = J_00*J_11 - J_01*J_10;
1005
1006
// Compute inverse of Jacobian
991
// Compute Jacobian
992
double J[4];
993
compute_jacobian_triangle_2d(J, vertex_coordinates);
994
995
// Compute Jacobian inverse and determinant
996
double K[4];
997
double detJ;
998
compute_jacobian_inverse_triangle_2d(K, detJ, J);
999
1007
1000
1008
1001
// Compute constants
1009
const double C0 = x[1][0] + x[2][0];
1010
const double C1 = x[1][1] + x[2][1];
1002
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1003
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1011
1004
1012
1005
// Get coordinates and map to the reference (FIAT) element
1013
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1014
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
1006
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1007
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1015
1008
1016
// Reset values.
1009
// Reset values
1017
1010
values[0] = 0.0;
1018
1011
values[1] = 0.0;
1019
1012
switch (i)
1021
1014
case 0:
1022
1015
{
1023
1016
1024
// Array of basisvalues.
1017
// Array of basisvalues
1025
1018
double basisvalues[3] = {0.0, 0.0, 0.0};
1026
1019
1027
// Declare helper variables.
1020
// Declare helper variables
1028
1021
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1029
1022
1030
// Compute basisvalues.
1023
// Compute basisvalues
1031
1024
basisvalues[0] = 1.0;
1032
1025
basisvalues[1] = tmp0;
1033
1026
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1035
1028
basisvalues[2] *= std::sqrt(1.0);
1036
1029
basisvalues[1] *= std::sqrt(3.0);
1037
1030
1038
// Table(s) of coefficients.
1031
// Table(s) of coefficients
1039
1032
static const double coefficients0[3] = \
1040
1033
{0.47140452, -0.28867513, -0.16666667};
1041
1034
1042
// Compute value(s).
1035
// Compute value(s)
1043
1036
for (unsigned int r = 0; r < 3; r++)
1044
1037
{
1045
1038
values[0] += coefficients0[r]*basisvalues[r];
1049
1042
case 1:
1050
1043
{
1051
1044
1052
// Array of basisvalues.
1045
// Array of basisvalues
1053
1046
double basisvalues[3] = {0.0, 0.0, 0.0};
1054
1047
1055
// Declare helper variables.
1048
// Declare helper variables
1056
1049
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1057
1050
1058
// Compute basisvalues.
1051
// Compute basisvalues
1059
1052
basisvalues[0] = 1.0;
1060
1053
basisvalues[1] = tmp0;
1061
1054
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1063
1056
basisvalues[2] *= std::sqrt(1.0);
1064
1057
basisvalues[1] *= std::sqrt(3.0);
1065
1058
1066
// Table(s) of coefficients.
1059
// Table(s) of coefficients
1067
1060
static const double coefficients0[3] = \
1068
1061
{0.47140452, 0.28867513, -0.16666667};
1069
1062
1070
// Compute value(s).
1063
// Compute value(s)
1071
1064
for (unsigned int r = 0; r < 3; r++)
1072
1065
{
1073
1066
values[0] += coefficients0[r]*basisvalues[r];
1077
1070
case 2:
1078
1071
{
1079
1072
1080
// Array of basisvalues.
1073
// Array of basisvalues
1081
1074
double basisvalues[3] = {0.0, 0.0, 0.0};
1082
1075
1083
// Declare helper variables.
1076
// Declare helper variables
1084
1077
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1085
1078
1086
// Compute basisvalues.
1079
// Compute basisvalues
1087
1080
basisvalues[0] = 1.0;
1088
1081
basisvalues[1] = tmp0;
1089
1082
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1091
1084
basisvalues[2] *= std::sqrt(1.0);
1092
1085
basisvalues[1] *= std::sqrt(3.0);
1093
1086
1094
// Table(s) of coefficients.
1087
// Table(s) of coefficients
1095
1088
static const double coefficients0[3] = \
1096
1089
{0.47140452, 0.0, 0.33333333};
1097
1090
1098
// Compute value(s).
1091
// Compute value(s)
1099
1092
for (unsigned int r = 0; r < 3; r++)
1100
1093
{
1101
1094
values[0] += coefficients0[r]*basisvalues[r];
1105
1098
case 3:
1106
1099
{
1107
1100
1108
// Array of basisvalues.
1101
// Array of basisvalues
1109
1102
double basisvalues[3] = {0.0, 0.0, 0.0};
1110
1103
1111
// Declare helper variables.
1104
// Declare helper variables
1112
1105
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1113
1106
1114
// Compute basisvalues.
1107
// Compute basisvalues
1115
1108
basisvalues[0] = 1.0;
1116
1109
basisvalues[1] = tmp0;
1117
1110
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1119
1112
basisvalues[2] *= std::sqrt(1.0);
1120
1113
basisvalues[1] *= std::sqrt(3.0);
1121
1114
1122
// Table(s) of coefficients.
1115
// Table(s) of coefficients
1123
1116
static const double coefficients0[3] = \
1124
1117
{0.47140452, -0.28867513, -0.16666667};
1125
1118
1126
// Compute value(s).
1119
// Compute value(s)
1127
1120
for (unsigned int r = 0; r < 3; r++)
1128
1121
{
1129
1122
values[1] += coefficients0[r]*basisvalues[r];
1133
1126
case 4:
1134
1127
{
1135
1128
1136
// Array of basisvalues.
1129
// Array of basisvalues
1137
1130
double basisvalues[3] = {0.0, 0.0, 0.0};
1138
1131
1139
// Declare helper variables.
1132
// Declare helper variables
1140
1133
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1141
1134
1142
// Compute basisvalues.
1135
// Compute basisvalues
1143
1136
basisvalues[0] = 1.0;
1144
1137
basisvalues[1] = tmp0;
1145
1138
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1147
1140
basisvalues[2] *= std::sqrt(1.0);
1148
1141
basisvalues[1] *= std::sqrt(3.0);
1149
1142
1150
// Table(s) of coefficients.
1143
// Table(s) of coefficients
1151
1144
static const double coefficients0[3] = \
1152
1145
{0.47140452, 0.28867513, -0.16666667};
1153
1146
1154
// Compute value(s).
1147
// Compute value(s)
1155
1148
for (unsigned int r = 0; r < 3; r++)
1156
1149
{
1157
1150
values[1] += coefficients0[r]*basisvalues[r];
1161
1154
case 5:
1162
1155
{
1163
1156
1164
// Array of basisvalues.
1157
// Array of basisvalues
1165
1158
double basisvalues[3] = {0.0, 0.0, 0.0};
1166
1159
1167
// Declare helper variables.
1160
// Declare helper variables
1168
1161
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1169
1162
1170
// Compute basisvalues.
1163
// Compute basisvalues
1171
1164
basisvalues[0] = 1.0;
1172
1165
basisvalues[1] = tmp0;
1173
1166
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1175
1168
basisvalues[2] *= std::sqrt(1.0);
1176
1169
basisvalues[1] *= std::sqrt(3.0);
1177
1170
1178
// Table(s) of coefficients.
1171
// Table(s) of coefficients
1179
1172
static const double coefficients0[3] = \
1180
1173
{0.47140452, 0.0, 0.33333333};
1181
1174
1182
// Compute value(s).
1175
// Compute value(s)
1183
1176
for (unsigned int r = 0; r < 3; r++)
1184
1177
{
1185
1178
values[1] += coefficients0[r]*basisvalues[r];
1190
1183
1191
1184
}
1192
1185
1193
/// Evaluate all basis functions at given point in cell
1186
/// Evaluate all basis functions at given point x in cell
1194
1187
virtual void evaluate_basis_all(double* values,
1195
const double* coordinates,
1196
const ufc::cell& c) const
1188
const double* x,
1189
const double* vertex_coordinates,
1190
int cell_orientation) const
1197
1191
{
1198
1192
// Helper variable to hold values of a single dof.
1199
1193
double dof_values[2] = {0.0, 0.0};
1200
1194
1201
// Loop dofs and call evaluate_basis.
1195
// Loop dofs and call evaluate_basis
1202
1196
for (unsigned int r = 0; r < 6; r++)
1203
1197
{
1204
evaluate_basis(r, dof_values, coordinates, c);
1198
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
1205
1199
for (unsigned int s = 0; s < 2; s++)
1206
1200
{
1207
1201
values[r*2 + s] = dof_values[s];
1209
1203
}// end loop over 'r'
1210
1204
}
1211
1205
1212
/// Evaluate order n derivatives of basis function i at given point in cell
1213
virtual void evaluate_basis_derivatives(unsigned int i,
1214
unsigned int n,
1206
/// Evaluate order n derivatives of basis function i at given point x in cell
1207
virtual void evaluate_basis_derivatives(std::size_t i,
1208
std::size_t n,
1215
1209
double* values,
1216
const double* coordinates,
1217
const ufc::cell& c) const
1210
const double* x,
1211
const double* vertex_coordinates,
1212
int cell_orientation) const
1218
1213
{
1219
// Extract vertex coordinates
1220
const double * const * x = c.coordinates;
1221
1222
// Compute Jacobian of affine map from reference cell
1223
const double J_00 = x[1][0] - x[0][0];
1224
const double J_01 = x[2][0] - x[0][0];
1225
const double J_10 = x[1][1] - x[0][1];
1226
const double J_11 = x[2][1] - x[0][1];
1227
1228
// Compute determinant of Jacobian
1229
double detJ = J_00*J_11 - J_01*J_10;
1230
1231
// Compute inverse of Jacobian
1232
const double K_00 = J_11 / detJ;
1233
const double K_01 = -J_01 / detJ;
1234
const double K_10 = -J_10 / detJ;
1235
const double K_11 = J_00 / detJ;
1214
// Compute Jacobian
1215
double J[4];
1216
compute_jacobian_triangle_2d(J, vertex_coordinates);
1217
1218
// Compute Jacobian inverse and determinant
1219
double K[4];
1220
double detJ;
1221
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1222
1236
1223
1237
1224
// Compute constants
1238
const double C0 = x[1][0] + x[2][0];
1239
const double C1 = x[1][1] + x[2][1];
1225
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1226
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1240
1227
1241
1228
// Get coordinates and map to the reference (FIAT) element
1242
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1243
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
1229
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1230
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1244
1231
1245
1232
// Compute number of derivatives.
1246
1233
unsigned int num_derivatives = 1;
1277
1264
}
1278
1265
1279
1266
// Compute inverse of Jacobian
1280
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
1267
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
1281
1268
1282
1269
// Declare transformation matrix
1283
1270
// Declare pointer to two dimensional array and initialise
1311
1298
case 0:
1312
1299
{
1313
1300
1314
// Array of basisvalues.
1301
// Array of basisvalues
1315
1302
double basisvalues[3] = {0.0, 0.0, 0.0};
1316
1303
1317
// Declare helper variables.
1304
// Declare helper variables
1318
1305
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1319
1306
1320
// Compute basisvalues.
1307
// Compute basisvalues
1321
1308
basisvalues[0] = 1.0;
1322
1309
basisvalues[1] = tmp0;
1323
1310
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1325
1312
basisvalues[2] *= std::sqrt(1.0);
1326
1313
basisvalues[1] *= std::sqrt(3.0);
1327
1314
1328
// Table(s) of coefficients.
1315
// Table(s) of coefficients
1329
1316
static const double coefficients0[3] = \
1330
1317
{0.47140452, -0.28867513, -0.16666667};
1331
1318
1457
1444
case 1:
1458
1445
{
1459
1446
1460
// Array of basisvalues.
1447
// Array of basisvalues
1461
1448
double basisvalues[3] = {0.0, 0.0, 0.0};
1462
1449
1463
// Declare helper variables.
1450
// Declare helper variables
1464
1451
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1465
1452
1466
// Compute basisvalues.
1453
// Compute basisvalues
1467
1454
basisvalues[0] = 1.0;
1468
1455
basisvalues[1] = tmp0;
1469
1456
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1471
1458
basisvalues[2] *= std::sqrt(1.0);
1472
1459
basisvalues[1] *= std::sqrt(3.0);
1473
1460
1474
// Table(s) of coefficients.
1461
// Table(s) of coefficients
1475
1462
static const double coefficients0[3] = \
1476
1463
{0.47140452, 0.28867513, -0.16666667};
1477
1464
1603
1590
case 2:
1604
1591
{
1605
1592
1606
// Array of basisvalues.
1593
// Array of basisvalues
1607
1594
double basisvalues[3] = {0.0, 0.0, 0.0};
1608
1595
1609
// Declare helper variables.
1596
// Declare helper variables
1610
1597
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1611
1598
1612
// Compute basisvalues.
1599
// Compute basisvalues
1613
1600
basisvalues[0] = 1.0;
1614
1601
basisvalues[1] = tmp0;
1615
1602
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1617
1604
basisvalues[2] *= std::sqrt(1.0);
1618
1605
basisvalues[1] *= std::sqrt(3.0);
1619
1606
1620
// Table(s) of coefficients.
1607
// Table(s) of coefficients
1621
1608
static const double coefficients0[3] = \
1622
1609
{0.47140452, 0.0, 0.33333333};
1623
1610
1749
1736
case 3:
1750
1737
{
1751
1738
1752
// Array of basisvalues.
1739
// Array of basisvalues
1753
1740
double basisvalues[3] = {0.0, 0.0, 0.0};
1754
1741
1755
// Declare helper variables.
1742
// Declare helper variables
1756
1743
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1757
1744
1758
// Compute basisvalues.
1745
// Compute basisvalues
1759
1746
basisvalues[0] = 1.0;
1760
1747
basisvalues[1] = tmp0;
1761
1748
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1763
1750
basisvalues[2] *= std::sqrt(1.0);
1764
1751
basisvalues[1] *= std::sqrt(3.0);
1765
1752
1766
// Table(s) of coefficients.
1753
// Table(s) of coefficients
1767
1754
static const double coefficients0[3] = \
1768
1755
{0.47140452, -0.28867513, -0.16666667};
1769
1756
1895
1882
case 4:
1896
1883
{
1897
1884
1898
// Array of basisvalues.
1885
// Array of basisvalues
1899
1886
double basisvalues[3] = {0.0, 0.0, 0.0};
1900
1887
1901
// Declare helper variables.
1888
// Declare helper variables
1902
1889
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1903
1890
1904
// Compute basisvalues.
1891
// Compute basisvalues
1905
1892
basisvalues[0] = 1.0;
1906
1893
basisvalues[1] = tmp0;
1907
1894
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1909
1896
basisvalues[2] *= std::sqrt(1.0);
1910
1897
basisvalues[1] *= std::sqrt(3.0);
1911
1898
1912
// Table(s) of coefficients.
1899
// Table(s) of coefficients
1913
1900
static const double coefficients0[3] = \
1914
1901
{0.47140452, 0.28867513, -0.16666667};
1915
1902
2041
2028
case 5:
2042
2029
{
2043
2030
2044
// Array of basisvalues.
2031
// Array of basisvalues
2045
2032
double basisvalues[3] = {0.0, 0.0, 0.0};
2046
2033
2047
// Declare helper variables.
2034
// Declare helper variables
2048
2035
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2049
2036
2050
// Compute basisvalues.
2037
// Compute basisvalues
2051
2038
basisvalues[0] = 1.0;
2052
2039
basisvalues[1] = tmp0;
2053
2040
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2055
2042
basisvalues[2] *= std::sqrt(1.0);
2056
2043
basisvalues[1] *= std::sqrt(3.0);
2057
2044
2058
// Table(s) of coefficients.
2045
// Table(s) of coefficients
2059
2046
static const double coefficients0[3] = \
2060
2047
{0.47140452, 0.0, 0.33333333};
2061
2048
2188
2175
2189
2176
}
2190
2177
2191
/// Evaluate order n derivatives of all basis functions at given point in cell
2192
virtual void evaluate_basis_derivatives_all(unsigned int n,
2178
/// Evaluate order n derivatives of all basis functions at given point x in cell
2179
virtual void evaluate_basis_derivatives_all(std::size_t n,
2193
2180
double* values,
2194
const double* coordinates,
2195
const ufc::cell& c) const
2181
const double* x,
2182
const double* vertex_coordinates,
2183
int cell_orientation) const
2196
2184
{
2197
2185
// Compute number of derivatives.
2198
2186
unsigned int num_derivatives = 1;
2211
2199
// Loop dofs and call evaluate_basis_derivatives.
2212
2200
for (unsigned int r = 0; r < 6; r++)
2213
2201
{
2214
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
2202
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2215
2203
for (unsigned int s = 0; s < 2*num_derivatives; s++)
2216
2204
{
2217
2205
values[r*2*num_derivatives + s] = dof_values[s];
2223
2211
}
2224
2212
2225
2213
/// Evaluate linear functional for dof i on the function f
2226
virtual double evaluate_dof(unsigned int i,
2214
virtual double evaluate_dof(std::size_t i,
2227
2215
const ufc::function& f,
2216
const double* vertex_coordinates,
2217
int cell_orientation,
2228
2218
const ufc::cell& c) const
2229
2219
{
2230
// Declare variables for result of evaluation.
2220
// Declare variables for result of evaluation
2231
2221
double vals[2];
2232
2222
2233
// Declare variable for physical coordinates.
2223
// Declare variable for physical coordinates
2234
2224
double y[2];
2235
const double * const * x = c.coordinates;
2236
2225
switch (i)
2237
2226
{
2238
2227
case 0:
2239
2228
{
2240
y[0] = x[0][0];
2241
y[1] = x[0][1];
2229
y[0] = vertex_coordinates[0];
2230
y[1] = vertex_coordinates[1];
2242
2231
f.evaluate(vals, y, c);
2243
2232
return vals[0];
2244
2233
break;
2245
2234
}
2246
2235
case 1:
2247
2236
{
2248
y[0] = x[1][0];
2249
y[1] = x[1][1];
2237
y[0] = vertex_coordinates[2];
2238
y[1] = vertex_coordinates[3];
2250
2239
f.evaluate(vals, y, c);
2251
2240
return vals[0];
2252
2241
break;
2253
2242
}
2254
2243
case 2:
2255
2244
{
2256
y[0] = x[2][0];
2257
y[1] = x[2][1];
2245
y[0] = vertex_coordinates[4];
2246
y[1] = vertex_coordinates[5];
2258
2247
f.evaluate(vals, y, c);
2259
2248
return vals[0];
2260
2249
break;
2261
2250
}
2262
2251
case 3:
2263
2252
{
2264
y[0] = x[0][0];
2265
y[1] = x[0][1];
2253
y[0] = vertex_coordinates[0];
2254
y[1] = vertex_coordinates[1];
2266
2255
f.evaluate(vals, y, c);
2267
2256
return vals[1];
2268
2257
break;
2269
2258
}
2270
2259
case 4:
2271
2260
{
2272
y[0] = x[1][0];
2273
y[1] = x[1][1];
2261
y[0] = vertex_coordinates[2];
2262
y[1] = vertex_coordinates[3];
2274
2263
f.evaluate(vals, y, c);
2275
2264
return vals[1];
2276
2265
break;
2277
2266
}
2278
2267
case 5:
2279
2268
{
2280
y[0] = x[2][0];
2281
y[1] = x[2][1];
2269
y[0] = vertex_coordinates[4];
2270
y[1] = vertex_coordinates[5];
2282
2271
f.evaluate(vals, y, c);
2283
2272
return vals[1];
2284
2273
break;
2291
2280
/// Evaluate linear functionals for all dofs on the function f
2292
2281
virtual void evaluate_dofs(double* values,
2293
2282
const ufc::function& f,
2283
const double* vertex_coordinates,
2284
int cell_orientation,
2294
2285
const ufc::cell& c) const
2295
2286
{
2296
// Declare variables for result of evaluation.
2287
// Declare variables for result of evaluation
2297
2288
double vals[2];
2298
2289
2299
// Declare variable for physical coordinates.
2290
// Declare variable for physical coordinates
2300
2291
double y[2];
2301
const double * const * x = c.coordinates;
2302
y[0] = x[0][0];
2303
y[1] = x[0][1];
2292
y[0] = vertex_coordinates[0];
2293
y[1] = vertex_coordinates[1];
2304
2294
f.evaluate(vals, y, c);
2305
2295
values[0] = vals[0];
2306
y[0] = x[1][0];
2307
y[1] = x[1][1];
2296
y[0] = vertex_coordinates[2];
2297
y[1] = vertex_coordinates[3];
2308
2298
f.evaluate(vals, y, c);
2309
2299
values[1] = vals[0];
2310
y[0] = x[2][0];
2311
y[1] = x[2][1];
2300
y[0] = vertex_coordinates[4];
2301
y[1] = vertex_coordinates[5];
2312
2302
f.evaluate(vals, y, c);
2313
2303
values[2] = vals[0];
2314
y[0] = x[0][0];
2315
y[1] = x[0][1];
2304
y[0] = vertex_coordinates[0];
2305
y[1] = vertex_coordinates[1];
2316
2306
f.evaluate(vals, y, c);
2317
2307
values[3] = vals[1];
2318
y[0] = x[1][0];
2319
y[1] = x[1][1];
2308
y[0] = vertex_coordinates[2];
2309
y[1] = vertex_coordinates[3];
2320
2310
f.evaluate(vals, y, c);
2321
2311
values[4] = vals[1];
2322
y[0] = x[2][0];
2323
y[1] = x[2][1];
2312
y[0] = vertex_coordinates[4];
2313
y[1] = vertex_coordinates[5];
2324
2314
f.evaluate(vals, y, c);
2325
2315
values[5] = vals[1];
2326
2316
}
2328
2318
/// Interpolate vertex values from dof values
2329
2319
virtual void interpolate_vertex_values(double* vertex_values,
2330
2320
const double* dof_values,
2321
const double* vertex_coordinates,
2322
int cell_orientation,
2331
2323
const ufc::cell& c) const
2332
2324
{
2333
2325
// Evaluate function and change variables
2345
2337
const double* xhat,
2346
2338
const ufc::cell& c) const
2347
2339
{
2348
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
2340
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
2349
2341
}
2350
2342
2351
2343
/// Map from coordinate x in cell to coordinate xhat in reference cell
2353
2345
const double* x,
2354
2346
const ufc::cell& c) const
2355
2347
{
2356
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
2348
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
2357
2349
}
2358
2350
2359
2351
/// Return the number of sub elements (for a mixed element)
2360
virtual unsigned int num_sub_elements() const
2352
virtual std::size_t num_sub_elements() const
2361
2353
{
2362
2354
return 2;
2363
2355
}
2364
2356
2365
2357
/// Create a new finite element for sub element i (for a mixed element)
2366
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2358
virtual ufc::finite_element* create_sub_element(std::size_t i) const
2367
2359
{
2368
2360
switch (i)
2369
2361
{
2411
2403
/// Return a string identifying the finite element
2412
2404
virtual const char* signature() const
2413
2405
{
2414
return "MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) })";
2406
return "MixedElement(*[VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (3,) })";
2415
2407
}
2416
2408
2417
2409
/// Return the cell shape
2421
2413
}
2422
2414
2423
2415
/// Return the topological dimension of the cell shape
2424
virtual unsigned int topological_dimension() const
2416
virtual std::size_t topological_dimension() const
2425
2417
{
2426
2418
return 2;
2427
2419
}
2428
2420
2429
2421
/// Return the geometric dimension of the cell shape
2430
virtual unsigned int geometric_dimension() const
2422
virtual std::size_t geometric_dimension() const
2431
2423
{
2432
2424
return 2;
2433
2425
}
2434
2426
2435
2427
/// Return the dimension of the finite element function space
2436
virtual unsigned int space_dimension() const
2428
virtual std::size_t space_dimension() const
2437
2429
{
2438
2430
return 9;
2439
2431
}
2440
2432
2441
2433
/// Return the rank of the value space
2442
virtual unsigned int value_rank() const
2434
virtual std::size_t value_rank() const
2443
2435
{
2444
2436
return 1;
2445
2437
}
2446
2438
2447
2439
/// Return the dimension of the value space for axis i
2448
virtual unsigned int value_dimension(unsigned int i) const
2440
virtual std::size_t value_dimension(std::size_t i) const
2449
2441
{
2450
2442
switch (i)
2451
2443
{
2459
2451
return 0;
2460
2452
}
2461
2453
2462
/// Evaluate basis function i at given point in cell
2463
virtual void evaluate_basis(unsigned int i,
2454
/// Evaluate basis function i at given point x in cell
2455
virtual void evaluate_basis(std::size_t i,
2464
2456
double* values,
2465
const double* coordinates,
2466
const ufc::cell& c) const
2457
const double* x,
2458
const double* vertex_coordinates,
2459
int cell_orientation) const
2467
2460
{
2468
// Extract vertex coordinates
2469
const double * const * x = c.coordinates;
2470
2471
// Compute Jacobian of affine map from reference cell
2472
const double J_00 = x[1][0] - x[0][0];
2473
const double J_01 = x[2][0] - x[0][0];
2474
const double J_10 = x[1][1] - x[0][1];
2475
const double J_11 = x[2][1] - x[0][1];
2476
2477
// Compute determinant of Jacobian
2478
double detJ = J_00*J_11 - J_01*J_10;
2479
2480
// Compute inverse of Jacobian
2461
// Compute Jacobian
2462
double J[4];
2463
compute_jacobian_triangle_2d(J, vertex_coordinates);
2464
2465
// Compute Jacobian inverse and determinant
2466
double K[4];
2467
double detJ;
2468
compute_jacobian_inverse_triangle_2d(K, detJ, J);
2469
2481
2470
2482
2471
// Compute constants
2483
const double C0 = x[1][0] + x[2][0];
2484
const double C1 = x[1][1] + x[2][1];
2472
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
2473
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
2485
2474
2486
2475
// Get coordinates and map to the reference (FIAT) element
2487
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
2488
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
2476
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
2477
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
2489
2478
2490
// Reset values.
2479
// Reset values
2491
2480
values[0] = 0.0;
2492
2481
values[1] = 0.0;
2493
2482
values[2] = 0.0;
2496
2485
case 0:
2497
2486
{
2498
2487
2499
// Array of basisvalues.
2488
// Array of basisvalues
2500
2489
double basisvalues[3] = {0.0, 0.0, 0.0};
2501
2490
2502
// Declare helper variables.
2491
// Declare helper variables
2503
2492
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2504
2493
2505
// Compute basisvalues.
2494
// Compute basisvalues
2506
2495
basisvalues[0] = 1.0;
2507
2496
basisvalues[1] = tmp0;
2508
2497
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2510
2499
basisvalues[2] *= std::sqrt(1.0);
2511
2500
basisvalues[1] *= std::sqrt(3.0);
2512
2501
2513
// Table(s) of coefficients.
2502
// Table(s) of coefficients
2514
2503
static const double coefficients0[3] = \
2515
2504
{0.47140452, -0.28867513, -0.16666667};
2516
2505
2517
// Compute value(s).
2506
// Compute value(s)
2518
2507
for (unsigned int r = 0; r < 3; r++)
2519
2508
{
2520
2509
values[0] += coefficients0[r]*basisvalues[r];
2524
2513
case 1:
2525
2514
{
2526
2515
2527
// Array of basisvalues.
2516
// Array of basisvalues
2528
2517
double basisvalues[3] = {0.0, 0.0, 0.0};
2529
2518
2530
// Declare helper variables.
2519
// Declare helper variables
2531
2520
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2532
2521
2533
// Compute basisvalues.
2522
// Compute basisvalues
2534
2523
basisvalues[0] = 1.0;
2535
2524
basisvalues[1] = tmp0;
2536
2525
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2538
2527
basisvalues[2] *= std::sqrt(1.0);
2539
2528
basisvalues[1] *= std::sqrt(3.0);
2540
2529
2541
// Table(s) of coefficients.
2530
// Table(s) of coefficients
2542
2531
static const double coefficients0[3] = \
2543
2532
{0.47140452, 0.28867513, -0.16666667};
2544
2533
2545
// Compute value(s).
2534
// Compute value(s)
2546
2535
for (unsigned int r = 0; r < 3; r++)
2547
2536
{
2548
2537
values[0] += coefficients0[r]*basisvalues[r];
2552
2541
case 2:
2553
2542
{
2554
2543
2555
// Array of basisvalues.
2544
// Array of basisvalues
2556
2545
double basisvalues[3] = {0.0, 0.0, 0.0};
2557
2546
2558
// Declare helper variables.
2547
// Declare helper variables
2559
2548
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2560
2549
2561
// Compute basisvalues.
2550
// Compute basisvalues
2562
2551
basisvalues[0] = 1.0;
2563
2552
basisvalues[1] = tmp0;
2564
2553
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2566
2555
basisvalues[2] *= std::sqrt(1.0);
2567
2556
basisvalues[1] *= std::sqrt(3.0);
2568
2557
2569
// Table(s) of coefficients.
2558
// Table(s) of coefficients
2570
2559
static const double coefficients0[3] = \
2571
2560
{0.47140452, 0.0, 0.33333333};
2572
2561
2573
// Compute value(s).
2562
// Compute value(s)
2574
2563
for (unsigned int r = 0; r < 3; r++)
2575
2564
{
2576
2565
values[0] += coefficients0[r]*basisvalues[r];
2580
2569
case 3:
2581
2570
{
2582
2571
2583
// Array of basisvalues.
2572
// Array of basisvalues
2584
2573
double basisvalues[3] = {0.0, 0.0, 0.0};
2585
2574
2586
// Declare helper variables.
2575
// Declare helper variables
2587
2576
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2588
2577
2589
// Compute basisvalues.
2578
// Compute basisvalues
2590
2579
basisvalues[0] = 1.0;
2591
2580
basisvalues[1] = tmp0;
2592
2581
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2594
2583
basisvalues[2] *= std::sqrt(1.0);
2595
2584
basisvalues[1] *= std::sqrt(3.0);
2596
2585
2597
// Table(s) of coefficients.
2586
// Table(s) of coefficients
2598
2587
static const double coefficients0[3] = \
2599
2588
{0.47140452, -0.28867513, -0.16666667};
2600
2589
2601
// Compute value(s).
2590
// Compute value(s)
2602
2591
for (unsigned int r = 0; r < 3; r++)
2603
2592
{
2604
2593
values[1] += coefficients0[r]*basisvalues[r];
2608
2597
case 4:
2609
2598
{
2610
2599
2611
// Array of basisvalues.
2600
// Array of basisvalues
2612
2601
double basisvalues[3] = {0.0, 0.0, 0.0};
2613
2602
2614
// Declare helper variables.
2603
// Declare helper variables
2615
2604
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2616
2605
2617
// Compute basisvalues.
2606
// Compute basisvalues
2618
2607
basisvalues[0] = 1.0;
2619
2608
basisvalues[1] = tmp0;
2620
2609
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2622
2611
basisvalues[2] *= std::sqrt(1.0);
2623
2612
basisvalues[1] *= std::sqrt(3.0);
2624
2613
2625
// Table(s) of coefficients.
2614
// Table(s) of coefficients
2626
2615
static const double coefficients0[3] = \
2627
2616
{0.47140452, 0.28867513, -0.16666667};
2628
2617
2629
// Compute value(s).
2618
// Compute value(s)
2630
2619
for (unsigned int r = 0; r < 3; r++)
2631
2620
{
2632
2621
values[1] += coefficients0[r]*basisvalues[r];
2636
2625
case 5:
2637
2626
{
2638
2627
2639
// Array of basisvalues.
2628
// Array of basisvalues
2640
2629
double basisvalues[3] = {0.0, 0.0, 0.0};
2641
2630
2642
// Declare helper variables.
2631
// Declare helper variables
2643
2632
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2644
2633
2645
// Compute basisvalues.
2634
// Compute basisvalues
2646
2635
basisvalues[0] = 1.0;
2647
2636
basisvalues[1] = tmp0;
2648
2637
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2650
2639
basisvalues[2] *= std::sqrt(1.0);
2651
2640
basisvalues[1] *= std::sqrt(3.0);
2652
2641
2653
// Table(s) of coefficients.
2642
// Table(s) of coefficients
2654
2643
static const double coefficients0[3] = \
2655
2644
{0.47140452, 0.0, 0.33333333};
2656
2645
2657
// Compute value(s).
2646
// Compute value(s)
2658
2647
for (unsigned int r = 0; r < 3; r++)
2659
2648
{
2660
2649
values[1] += coefficients0[r]*basisvalues[r];
2664
2653
case 6:
2665
2654
{
2666
2655
2667
// Array of basisvalues.
2656
// Array of basisvalues
2668
2657
double basisvalues[3] = {0.0, 0.0, 0.0};
2669
2658
2670
// Declare helper variables.
2659
// Declare helper variables
2671
2660
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2672
2661
2673
// Compute basisvalues.
2662
// Compute basisvalues
2674
2663
basisvalues[0] = 1.0;
2675
2664
basisvalues[1] = tmp0;
2676
2665
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2678
2667
basisvalues[2] *= std::sqrt(1.0);
2679
2668
basisvalues[1] *= std::sqrt(3.0);
2680
2669
2681
// Table(s) of coefficients.
2670
// Table(s) of coefficients
2682
2671
static const double coefficients0[3] = \
2683
2672
{0.47140452, -0.28867513, -0.16666667};
2684
2673
2685
// Compute value(s).
2674
// Compute value(s)
2686
2675
for (unsigned int r = 0; r < 3; r++)
2687
2676
{
2688
2677
values[2] += coefficients0[r]*basisvalues[r];
2692
2681
case 7:
2693
2682
{
2694
2683
2695
// Array of basisvalues.
2684
// Array of basisvalues
2696
2685
double basisvalues[3] = {0.0, 0.0, 0.0};
2697
2686
2698
// Declare helper variables.
2687
// Declare helper variables
2699
2688
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2700
2689
2701
// Compute basisvalues.
2690
// Compute basisvalues
2702
2691
basisvalues[0] = 1.0;
2703
2692
basisvalues[1] = tmp0;
2704
2693
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2706
2695
basisvalues[2] *= std::sqrt(1.0);
2707
2696
basisvalues[1] *= std::sqrt(3.0);
2708
2697
2709
// Table(s) of coefficients.
2698
// Table(s) of coefficients
2710
2699
static const double coefficients0[3] = \
2711
2700
{0.47140452, 0.28867513, -0.16666667};
2712
2701
2713
// Compute value(s).
2702
// Compute value(s)
2714
2703
for (unsigned int r = 0; r < 3; r++)
2715
2704
{
2716
2705
values[2] += coefficients0[r]*basisvalues[r];
2720
2709
case 8:
2721
2710
{
2722
2711
2723
// Array of basisvalues.
2712
// Array of basisvalues
2724
2713
double basisvalues[3] = {0.0, 0.0, 0.0};
2725
2714
2726
// Declare helper variables.
2715
// Declare helper variables
2727
2716
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2728
2717
2729
// Compute basisvalues.
2718
// Compute basisvalues
2730
2719
basisvalues[0] = 1.0;
2731
2720
basisvalues[1] = tmp0;
2732
2721
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2734
2723
basisvalues[2] *= std::sqrt(1.0);
2735
2724
basisvalues[1] *= std::sqrt(3.0);
2736
2725
2737
// Table(s) of coefficients.
2726
// Table(s) of coefficients
2738
2727
static const double coefficients0[3] = \
2739
2728
{0.47140452, 0.0, 0.33333333};
2740
2729
2741
// Compute value(s).
2730
// Compute value(s)
2742
2731
for (unsigned int r = 0; r < 3; r++)
2743
2732
{
2744
2733
values[2] += coefficients0[r]*basisvalues[r];
2749
2738
2750
2739
}
2751
2740
2752
/// Evaluate all basis functions at given point in cell
2741
/// Evaluate all basis functions at given point x in cell
2753
2742
virtual void evaluate_basis_all(double* values,
2754
const double* coordinates,
2755
const ufc::cell& c) const
2743
const double* x,
2744
const double* vertex_coordinates,
2745
int cell_orientation) const
2756
2746
{
2757
2747
// Helper variable to hold values of a single dof.
2758
2748
double dof_values[3] = {0.0, 0.0, 0.0};
2759
2749
2760
// Loop dofs and call evaluate_basis.
2750
// Loop dofs and call evaluate_basis
2761
2751
for (unsigned int r = 0; r < 9; r++)
2762
2752
{
2763
evaluate_basis(r, dof_values, coordinates, c);
2753
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
2764
2754
for (unsigned int s = 0; s < 3; s++)
2765
2755
{
2766
2756
values[r*3 + s] = dof_values[s];
2768
2758
}// end loop over 'r'
2769
2759
}
2770
2760
2771
/// Evaluate order n derivatives of basis function i at given point in cell
2772
virtual void evaluate_basis_derivatives(unsigned int i,
2773
unsigned int n,
2761
/// Evaluate order n derivatives of basis function i at given point x in cell
2762
virtual void evaluate_basis_derivatives(std::size_t i,
2763
std::size_t n,
2774
2764
double* values,
2775
const double* coordinates,
2776
const ufc::cell& c) const
2765
const double* x,
2766
const double* vertex_coordinates,
2767
int cell_orientation) const
2777
2768
{
2778
// Extract vertex coordinates
2779
const double * const * x = c.coordinates;
2780
2781
// Compute Jacobian of affine map from reference cell
2782
const double J_00 = x[1][0] - x[0][0];
2783
const double J_01 = x[2][0] - x[0][0];
2784
const double J_10 = x[1][1] - x[0][1];
2785
const double J_11 = x[2][1] - x[0][1];
2786
2787
// Compute determinant of Jacobian
2788
double detJ = J_00*J_11 - J_01*J_10;
2789
2790
// Compute inverse of Jacobian
2791
const double K_00 = J_11 / detJ;
2792
const double K_01 = -J_01 / detJ;
2793
const double K_10 = -J_10 / detJ;
2794
const double K_11 = J_00 / detJ;
2769
// Compute Jacobian
2770
double J[4];
2771
compute_jacobian_triangle_2d(J, vertex_coordinates);
2772
2773
// Compute Jacobian inverse and determinant
2774
double K[4];
2775
double detJ;
2776
compute_jacobian_inverse_triangle_2d(K, detJ, J);
2777
2795
2778
2796
2779
// Compute constants
2797
const double C0 = x[1][0] + x[2][0];
2798
const double C1 = x[1][1] + x[2][1];
2780
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
2781
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
2799
2782
2800
2783
// Get coordinates and map to the reference (FIAT) element
2801
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
2802
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
2784
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
2785
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
2803
2786
2804
2787
// Compute number of derivatives.
2805
2788
unsigned int num_derivatives = 1;
2836
2819
}
2837
2820
2838
2821
// Compute inverse of Jacobian
2839
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
2822
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
2840
2823
2841
2824
// Declare transformation matrix
2842
2825
// Declare pointer to two dimensional array and initialise
2870
2853
case 0:
2871
2854
{
2872
2855
2873
// Array of basisvalues.
2856
// Array of basisvalues
2874
2857
double basisvalues[3] = {0.0, 0.0, 0.0};
2875
2858
2876
// Declare helper variables.
2859
// Declare helper variables
2877
2860
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2878
2861
2879
// Compute basisvalues.
2862
// Compute basisvalues
2880
2863
basisvalues[0] = 1.0;
2881
2864
basisvalues[1] = tmp0;
2882
2865
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2884
2867
basisvalues[2] *= std::sqrt(1.0);
2885
2868
basisvalues[1] *= std::sqrt(3.0);
2886
2869
2887
// Table(s) of coefficients.
2870
// Table(s) of coefficients
2888
2871
static const double coefficients0[3] = \
2889
2872
{0.47140452, -0.28867513, -0.16666667};
2890
2873
3016
2999
case 1:
3017
3000
{
3018
3001
3019
// Array of basisvalues.
3002
// Array of basisvalues
3020
3003
double basisvalues[3] = {0.0, 0.0, 0.0};
3021
3004
3022
// Declare helper variables.
3005
// Declare helper variables
3023
3006
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3024
3007
3025
// Compute basisvalues.
3008
// Compute basisvalues
3026
3009
basisvalues[0] = 1.0;
3027
3010
basisvalues[1] = tmp0;
3028
3011
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3030
3013
basisvalues[2] *= std::sqrt(1.0);
3031
3014
basisvalues[1] *= std::sqrt(3.0);
3032
3015
3033
// Table(s) of coefficients.
3016
// Table(s) of coefficients
3034
3017
static const double coefficients0[3] = \
3035
3018
{0.47140452, 0.28867513, -0.16666667};
3036
3019
3162
3145
case 2:
3163
3146
{
3164
3147
3165
// Array of basisvalues.
3148
// Array of basisvalues
3166
3149
double basisvalues[3] = {0.0, 0.0, 0.0};
3167
3150
3168
// Declare helper variables.
3151
// Declare helper variables
3169
3152
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3170
3153
3171
// Compute basisvalues.
3154
// Compute basisvalues
3172
3155
basisvalues[0] = 1.0;
3173
3156
basisvalues[1] = tmp0;
3174
3157
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3176
3159
basisvalues[2] *= std::sqrt(1.0);
3177
3160
basisvalues[1] *= std::sqrt(3.0);
3178
3161
3179
// Table(s) of coefficients.
3162
// Table(s) of coefficients
3180
3163
static const double coefficients0[3] = \
3181
3164
{0.47140452, 0.0, 0.33333333};
3182
3165
3308
3291
case 3:
3309
3292
{
3310
3293
3311
// Array of basisvalues.
3294
// Array of basisvalues
3312
3295
double basisvalues[3] = {0.0, 0.0, 0.0};
3313
3296
3314
// Declare helper variables.
3297
// Declare helper variables
3315
3298
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3316
3299
3317
// Compute basisvalues.
3300
// Compute basisvalues
3318
3301
basisvalues[0] = 1.0;
3319
3302
basisvalues[1] = tmp0;
3320
3303
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3322
3305
basisvalues[2] *= std::sqrt(1.0);
3323
3306
basisvalues[1] *= std::sqrt(3.0);
3324
3307
3325
// Table(s) of coefficients.
3308
// Table(s) of coefficients
3326
3309
static const double coefficients0[3] = \
3327
3310
{0.47140452, -0.28867513, -0.16666667};
3328
3311
3454
3437
case 4:
3455
3438
{
3456
3439
3457
// Array of basisvalues.
3440
// Array of basisvalues
3458
3441
double basisvalues[3] = {0.0, 0.0, 0.0};
3459
3442
3460
// Declare helper variables.
3443
// Declare helper variables
3461
3444
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3462
3445
3463
// Compute basisvalues.
3446
// Compute basisvalues
3464
3447
basisvalues[0] = 1.0;
3465
3448
basisvalues[1] = tmp0;
3466
3449
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3468
3451
basisvalues[2] *= std::sqrt(1.0);
3469
3452
basisvalues[1] *= std::sqrt(3.0);
3470
3453
3471
// Table(s) of coefficients.
3454
// Table(s) of coefficients
3472
3455
static const double coefficients0[3] = \
3473
3456
{0.47140452, 0.28867513, -0.16666667};
3474
3457
3600
3583
case 5:
3601
3584
{
3602
3585
3603
// Array of basisvalues.
3586
// Array of basisvalues
3604
3587
double basisvalues[3] = {0.0, 0.0, 0.0};
3605
3588
3606
// Declare helper variables.
3589
// Declare helper variables
3607
3590
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3608
3591
3609
// Compute basisvalues.
3592
// Compute basisvalues
3610
3593
basisvalues[0] = 1.0;
3611
3594
basisvalues[1] = tmp0;
3612
3595
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3614
3597
basisvalues[2] *= std::sqrt(1.0);
3615
3598
basisvalues[1] *= std::sqrt(3.0);
3616
3599
3617
// Table(s) of coefficients.
3600
// Table(s) of coefficients
3618
3601
static const double coefficients0[3] = \
3619
3602
{0.47140452, 0.0, 0.33333333};
3620
3603
3746
3729
case 6:
3747
3730
{
3748
3731
3749
// Array of basisvalues.
3732
// Array of basisvalues
3750
3733
double basisvalues[3] = {0.0, 0.0, 0.0};
3751
3734
3752
// Declare helper variables.
3735
// Declare helper variables
3753
3736
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3754
3737
3755
// Compute basisvalues.
3738
// Compute basisvalues
3756
3739
basisvalues[0] = 1.0;
3757
3740
basisvalues[1] = tmp0;
3758
3741
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3760
3743
basisvalues[2] *= std::sqrt(1.0);
3761
3744
basisvalues[1] *= std::sqrt(3.0);
3762
3745
3763
// Table(s) of coefficients.
3746
// Table(s) of coefficients
3764
3747
static const double coefficients0[3] = \
3765
3748
{0.47140452, -0.28867513, -0.16666667};
3766
3749
3892
3875
case 7:
3893
3876
{
3894
3877
3895
// Array of basisvalues.
3878
// Array of basisvalues
3896
3879
double basisvalues[3] = {0.0, 0.0, 0.0};
3897
3880
3898
// Declare helper variables.
3881
// Declare helper variables
3899
3882
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3900
3883
3901
// Compute basisvalues.
3884
// Compute basisvalues
3902
3885
basisvalues[0] = 1.0;
3903
3886
basisvalues[1] = tmp0;
3904
3887
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3906
3889
basisvalues[2] *= std::sqrt(1.0);
3907
3890
basisvalues[1] *= std::sqrt(3.0);
3908
3891
3909
// Table(s) of coefficients.
3892
// Table(s) of coefficients
3910
3893
static const double coefficients0[3] = \
3911
3894
{0.47140452, 0.28867513, -0.16666667};
3912
3895
4038
4021
case 8:
4039
4022
{
4040
4023
4041
// Array of basisvalues.
4024
// Array of basisvalues
4042
4025
double basisvalues[3] = {0.0, 0.0, 0.0};
4043
4026
4044
// Declare helper variables.
4027
// Declare helper variables
4045
4028
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4046
4029
4047
// Compute basisvalues.
4030
// Compute basisvalues
4048
4031
basisvalues[0] = 1.0;
4049
4032
basisvalues[1] = tmp0;
4050
4033
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
4052
4035
basisvalues[2] *= std::sqrt(1.0);
4053
4036
basisvalues[1] *= std::sqrt(3.0);
4054
4037
4055
// Table(s) of coefficients.
4038
// Table(s) of coefficients
4056
4039
static const double coefficients0[3] = \
4057
4040
{0.47140452, 0.0, 0.33333333};
4058
4041
4185
4168
4186
4169
}
4187
4170
4188
/// Evaluate order n derivatives of all basis functions at given point in cell
4189
virtual void evaluate_basis_derivatives_all(unsigned int n,
4171
/// Evaluate order n derivatives of all basis functions at given point x in cell
4172
virtual void evaluate_basis_derivatives_all(std::size_t n,
4190
4173
double* values,
4191
const double* coordinates,
4192
const ufc::cell& c) const
4174
const double* x,
4175
const double* vertex_coordinates,
4176
int cell_orientation) const
4193
4177
{
4194
4178
// Compute number of derivatives.
4195
4179
unsigned int num_derivatives = 1;
4208
4192
// Loop dofs and call evaluate_basis_derivatives.
4209
4193
for (unsigned int r = 0; r < 9; r++)
4210
4194
{
4211
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
4195
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
4212
4196
for (unsigned int s = 0; s < 3*num_derivatives; s++)
4213
4197
{
4214
4198
values[r*3*num_derivatives + s] = dof_values[s];
4220
4204
}
4221
4205
4222
4206
/// Evaluate linear functional for dof i on the function f
4223
virtual double evaluate_dof(unsigned int i,
4207
virtual double evaluate_dof(std::size_t i,
4224
4208
const ufc::function& f,
4209
const double* vertex_coordinates,
4210
int cell_orientation,
4225
4211
const ufc::cell& c) const
4226
4212
{
4227
// Declare variables for result of evaluation.
4213
// Declare variables for result of evaluation
4228
4214
double vals[3];
4229
4215
4230
// Declare variable for physical coordinates.
4216
// Declare variable for physical coordinates
4231
4217
double y[2];
4232
const double * const * x = c.coordinates;
4233
4218
switch (i)
4234
4219
{
4235
4220
case 0:
4236
4221
{
4237
y[0] = x[0][0];
4238
y[1] = x[0][1];
4222
y[0] = vertex_coordinates[0];
4223
y[1] = vertex_coordinates[1];
4239
4224
f.evaluate(vals, y, c);
4240
4225
return vals[0];
4241
4226
break;
4242
4227
}
4243
4228
case 1:
4244
4229
{
4245
y[0] = x[1][0];
4246
y[1] = x[1][1];
4230
y[0] = vertex_coordinates[2];
4231
y[1] = vertex_coordinates[3];
4247
4232
f.evaluate(vals, y, c);
4248
4233
return vals[0];
4249
4234
break;
4250
4235
}
4251
4236
case 2:
4252
4237
{
4253
y[0] = x[2][0];
4254
y[1] = x[2][1];
4238
y[0] = vertex_coordinates[4];
4239
y[1] = vertex_coordinates[5];
4255
4240
f.evaluate(vals, y, c);
4256
4241
return vals[0];
4257
4242
break;
4258
4243
}
4259
4244
case 3:
4260
4245
{
4261
y[0] = x[0][0];
4262
y[1] = x[0][1];
4246
y[0] = vertex_coordinates[0];
4247
y[1] = vertex_coordinates[1];
4263
4248
f.evaluate(vals, y, c);
4264
4249
return vals[1];
4265
4250
break;
4266
4251
}
4267
4252
case 4:
4268
4253
{
4269
y[0] = x[1][0];
4270
y[1] = x[1][1];
4254
y[0] = vertex_coordinates[2];
4255
y[1] = vertex_coordinates[3];
4271
4256
f.evaluate(vals, y, c);
4272
4257
return vals[1];
4273
4258
break;
4274
4259
}
4275
4260
case 5:
4276
4261
{
4277
y[0] = x[2][0];
4278
y[1] = x[2][1];
4262
y[0] = vertex_coordinates[4];
4263
y[1] = vertex_coordinates[5];
4279
4264
f.evaluate(vals, y, c);
4280
4265
return vals[1];
4281
4266
break;
4282
4267
}
4283
4268
case 6:
4284
4269
{
4285
y[0] = x[0][0];
4286
y[1] = x[0][1];
4270
y[0] = vertex_coordinates[0];
4271
y[1] = vertex_coordinates[1];
4287
4272
f.evaluate(vals, y, c);
4288
4273
return vals[2];
4289
4274
break;
4290
4275
}
4291
4276
case 7:
4292
4277
{
4293
y[0] = x[1][0];
4294
y[1] = x[1][1];
4278
y[0] = vertex_coordinates[2];
4279
y[1] = vertex_coordinates[3];
4295
4280
f.evaluate(vals, y, c);
4296
4281
return vals[2];
4297
4282
break;
4298
4283
}
4299
4284
case 8:
4300
4285
{
4301
y[0] = x[2][0];
4302
y[1] = x[2][1];
4286
y[0] = vertex_coordinates[4];
4287
y[1] = vertex_coordinates[5];
4303
4288
f.evaluate(vals, y, c);
4304
4289
return vals[2];
4305
4290
break;
4312
4297
/// Evaluate linear functionals for all dofs on the function f
4313
4298
virtual void evaluate_dofs(double* values,
4314
4299
const ufc::function& f,
4300
const double* vertex_coordinates,
4301
int cell_orientation,
4315
4302
const ufc::cell& c) const
4316
4303
{
4317
// Declare variables for result of evaluation.
4304
// Declare variables for result of evaluation
4318
4305
double vals[3];
4319
4306
4320
// Declare variable for physical coordinates.
4307
// Declare variable for physical coordinates
4321
4308
double y[2];
4322
const double * const * x = c.coordinates;
4323
y[0] = x[0][0];
4324
y[1] = x[0][1];
4309
y[0] = vertex_coordinates[0];
4310
y[1] = vertex_coordinates[1];
4325
4311
f.evaluate(vals, y, c);
4326
4312
values[0] = vals[0];
4327
y[0] = x[1][0];
4328
y[1] = x[1][1];
4313
y[0] = vertex_coordinates[2];
4314
y[1] = vertex_coordinates[3];
4329
4315
f.evaluate(vals, y, c);
4330
4316
values[1] = vals[0];
4331
y[0] = x[2][0];
4332
y[1] = x[2][1];
4317
y[0] = vertex_coordinates[4];
4318
y[1] = vertex_coordinates[5];
4333
4319
f.evaluate(vals, y, c);
4334
4320
values[2] = vals[0];
4335
y[0] = x[0][0];
4336
y[1] = x[0][1];
4321
y[0] = vertex_coordinates[0];
4322
y[1] = vertex_coordinates[1];
4337
4323
f.evaluate(vals, y, c);
4338
4324
values[3] = vals[1];
4339
y[0] = x[1][0];
4340
y[1] = x[1][1];
4325
y[0] = vertex_coordinates[2];
4326
y[1] = vertex_coordinates[3];
4341
4327
f.evaluate(vals, y, c);
4342
4328
values[4] = vals[1];
4343
y[0] = x[2][0];
4344
y[1] = x[2][1];
4329
y[0] = vertex_coordinates[4];
4330
y[1] = vertex_coordinates[5];
4345
4331
f.evaluate(vals, y, c);
4346
4332
values[5] = vals[1];
4347
y[0] = x[0][0];
4348
y[1] = x[0][1];
4333
y[0] = vertex_coordinates[0];
4334
y[1] = vertex_coordinates[1];
4349
4335
f.evaluate(vals, y, c);
4350
4336
values[6] = vals[2];
4351
y[0] = x[1][0];
4352
y[1] = x[1][1];
4337
y[0] = vertex_coordinates[2];
4338
y[1] = vertex_coordinates[3];
4353
4339
f.evaluate(vals, y, c);
4354
4340
values[7] = vals[2];
4355
y[0] = x[2][0];
4356
y[1] = x[2][1];
4341
y[0] = vertex_coordinates[4];
4342
y[1] = vertex_coordinates[5];
4357
4343
f.evaluate(vals, y, c);
4358
4344
values[8] = vals[2];
4359
4345
}
4361
4347
/// Interpolate vertex values from dof values
4362
4348
virtual void interpolate_vertex_values(double* vertex_values,
4363
4349
const double* dof_values,
4350
const double* vertex_coordinates,
4351
int cell_orientation,
4364
4352
const ufc::cell& c) const
4365
4353
{
4366
4354
// Evaluate function and change variables
4382
4370
const double* xhat,
4383
4371
const ufc::cell& c) const
4384
4372
{
4385
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4373
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
4386
4374
}
4387
4375
4388
4376
/// Map from coordinate x in cell to coordinate xhat in reference cell
4390
4378
const double* x,
4391
4379
const ufc::cell& c) const
4392
4380
{
4393
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
4381
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
4394
4382
}
4395
4383
4396
4384
/// Return the number of sub elements (for a mixed element)
4397
virtual unsigned int num_sub_elements() const
4385
virtual std::size_t num_sub_elements() const
4398
4386
{
4399
4387
return 2;
4400
4388
}
4401
4389
4402
4390
/// Create a new finite element for sub element i (for a mixed element)
4403
virtual ufc::finite_element* create_sub_element(unsigned int i) const
4391
virtual ufc::finite_element* create_sub_element(std::size_t i) const
4404
4392
{
4405
4393
switch (i)
4406
4394
{
4432
4420
4433
4421
class stabilisedstokes_dofmap_0: public ufc::dofmap
4434
4422
{
4435
private:
4436
4437
unsigned int _global_dimension;
4438
4423
public:
4439
4424
4440
4425
/// Constructor
4441
4426
stabilisedstokes_dofmap_0() : ufc::dofmap()
4442
4427
{
4443
_global_dimension = 0;
4428
// Do nothing
4444
4429
}
4445
4430
4446
4431
/// Destructor
4452
4437
/// Return a string identifying the dofmap
4453
4438
virtual const char* signature() const
4454
4439
{
4455
return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)";
4440
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
4456
4441
}
4457
4442
4458
4443
/// Return true iff mesh entities of topological dimension d are needed
4459
virtual bool needs_mesh_entities(unsigned int d) const
4444
virtual bool needs_mesh_entities(std::size_t d) const
4460
4445
{
4461
4446
switch (d)
4462
4447
{
4480
4465
return false;
4481
4466
}
4482
4467
4483
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
4484
virtual bool init_mesh(const ufc::mesh& m)
4485
{
4486
_global_dimension = m.num_entities[0];
4487
return false;
4488
}
4489
4490
/// Initialize dofmap for given cell
4491
virtual void init_cell(const ufc::mesh& m,
4492
const ufc::cell& c)
4493
{
4494
// Do nothing
4495
}
4496
4497
/// Finish initialization of dofmap for cells
4498
virtual void init_cell_finalize()
4499
{
4500
// Do nothing
4501
}
4502
4503
4468
/// Return the topological dimension of the associated cell shape
4504
virtual unsigned int topological_dimension() const
4469
virtual std::size_t topological_dimension() const
4505
4470
{
4506
4471
return 2;
4507
4472
}
4508
4473
4509
4474
/// Return the geometric dimension of the associated cell shape
4510
virtual unsigned int geometric_dimension() const
4475
virtual std::size_t geometric_dimension() const
4511
4476
{
4512
4477
return 2;
4513
4478
}
4514
4479
4515
4480
/// Return the dimension of the global finite element function space
4516
virtual unsigned int global_dimension() const
4481
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4482
num_global_entities) const
4517
4483
{
4518
return _global_dimension;
4484
return num_global_entities[0];
4519
4485
}
4520
4486
4521
4487
/// Return the dimension of the local finite element function space for a cell
4522
virtual unsigned int local_dimension(const ufc::cell& c) const
4488
virtual std::size_t local_dimension(const ufc::cell& c) const
4523
4489
{
4524
4490
return 3;
4525
4491
}
4526
4492
4527
4493
/// Return the maximum dimension of the local finite element function space
4528
virtual unsigned int max_local_dimension() const
4494
virtual std::size_t max_local_dimension() const
4529
4495
{
4530
4496
return 3;
4531
4497
}
4532
4498
4533
4499
/// Return the number of dofs on each cell facet
4534
virtual unsigned int num_facet_dofs() const
4500
virtual std::size_t num_facet_dofs() const
4535
4501
{
4536
4502
return 2;
4537
4503
}
4538
4504
4539
4505
/// Return the number of dofs associated with each cell entity of dimension d
4540
virtual unsigned int num_entity_dofs(unsigned int d) const
4506
virtual std::size_t num_entity_dofs(std::size_t d) const
4541
4507
{
4542
4508
switch (d)
4543
4509
{
4562
4528
}
4563
4529
4564
4530
/// Tabulate the local-to-global mapping of dofs on a cell
4565
virtual void tabulate_dofs(unsigned int* dofs,
4566
const ufc::mesh& m,
4531
virtual void tabulate_dofs(std::size_t* dofs,
4532
const std::vector<std::size_t>& num_global_entities,
4567
4533
const ufc::cell& c) const
4568
4534
{
4569
4535
dofs[0] = c.entity_indices[0][0];
4572
4538
}
4573
4539
4574
4540
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4575
virtual void tabulate_facet_dofs(unsigned int* dofs,
4576
unsigned int facet) const
4541
virtual void tabulate_facet_dofs(std::size_t* dofs,
4542
std::size_t facet) const
4577
4543
{
4578
4544
switch (facet)
4579
4545
{
4600
4566
}
4601
4567
4602
4568
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4603
virtual void tabulate_entity_dofs(unsigned int* dofs,
4604
unsigned int d, unsigned int i) const
4569
virtual void tabulate_entity_dofs(std::size_t* dofs,
4570
std::size_t d, std::size_t i) const
4605
4571
{
4606
4572
if (d > 2)
4607
4573
{
4653
4619
}
4654
4620
4655
4621
/// Tabulate the coordinates of all dofs on a cell
4656
virtual void tabulate_coordinates(double** coordinates,
4657
const ufc::cell& c) const
4622
virtual void tabulate_coordinates(double** dof_coordinates,
4623
const double* vertex_coordinates) const
4658
4624
{
4659
const double * const * x = c.coordinates;
4660
4661
coordinates[0][0] = x[0][0];
4662
coordinates[0][1] = x[0][1];
4663
coordinates[1][0] = x[1][0];
4664
coordinates[1][1] = x[1][1];
4665
coordinates[2][0] = x[2][0];
4666
coordinates[2][1] = x[2][1];
4625
dof_coordinates[0][0] = vertex_coordinates[0];
4626
dof_coordinates[0][1] = vertex_coordinates[1];
4627
dof_coordinates[1][0] = vertex_coordinates[2];
4628
dof_coordinates[1][1] = vertex_coordinates[3];
4629
dof_coordinates[2][0] = vertex_coordinates[4];
4630
dof_coordinates[2][1] = vertex_coordinates[5];
4667
4631
}
4668
4632
4669
4633
/// Return the number of sub dofmaps (for a mixed element)
4670
virtual unsigned int num_sub_dofmaps() const
4634
virtual std::size_t num_sub_dofmaps() const
4671
4635
{
4672
4636
return 0;
4673
4637
}
4674
4638
4675
4639
/// Create a new dofmap for sub dofmap i (for a mixed element)
4676
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
4640
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4677
4641
{
4678
4642
return 0;
4679
4643
}
4691
4655
4692
4656
class stabilisedstokes_dofmap_1: public ufc::dofmap
4693
4657
{
4694
private:
4695
4696
unsigned int _global_dimension;
4697
4658
public:
4698
4659
4699
4660
/// Constructor
4700
4661
stabilisedstokes_dofmap_1() : ufc::dofmap()
4701
4662
{
4702
_global_dimension = 0;
4663
// Do nothing
4703
4664
}
4704
4665
4705
4666
/// Destructor
4711
4672
/// Return a string identifying the dofmap
4712
4673
virtual const char* signature() const
4713
4674
{
4714
return "FFC dofmap for VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None)";
4675
return "FFC dofmap for VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, 2, None)";
4715
4676
}
4716
4677
4717
4678
/// Return true iff mesh entities of topological dimension d are needed
4718
virtual bool needs_mesh_entities(unsigned int d) const
4679
virtual bool needs_mesh_entities(std::size_t d) const
4719
4680
{
4720
4681
switch (d)
4721
4682
{
4739
4700
return false;
4740
4701
}
4741
4702
4742
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
4743
virtual bool init_mesh(const ufc::mesh& m)
4744
{
4745
_global_dimension = 2*m.num_entities[0];
4746
return false;
4747
}
4748
4749
/// Initialize dofmap for given cell
4750
virtual void init_cell(const ufc::mesh& m,
4751
const ufc::cell& c)
4752
{
4753
// Do nothing
4754
}
4755
4756
/// Finish initialization of dofmap for cells
4757
virtual void init_cell_finalize()
4758
{
4759
// Do nothing
4760
}
4761
4762
4703
/// Return the topological dimension of the associated cell shape
4763
virtual unsigned int topological_dimension() const
4704
virtual std::size_t topological_dimension() const
4764
4705
{
4765
4706
return 2;
4766
4707
}
4767
4708
4768
4709
/// Return the geometric dimension of the associated cell shape
4769
virtual unsigned int geometric_dimension() const
4710
virtual std::size_t geometric_dimension() const
4770
4711
{
4771
4712
return 2;
4772
4713
}
4773
4714
4774
4715
/// Return the dimension of the global finite element function space
4775
virtual unsigned int global_dimension() const
4716
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4717
num_global_entities) const
4776
4718
{
4777
return _global_dimension;
4719
return 2*num_global_entities[0];
4778
4720
}
4779
4721
4780
4722
/// Return the dimension of the local finite element function space for a cell
4781
virtual unsigned int local_dimension(const ufc::cell& c) const
4723
virtual std::size_t local_dimension(const ufc::cell& c) const
4782
4724
{
4783
4725
return 6;
4784
4726
}
4785
4727
4786
4728
/// Return the maximum dimension of the local finite element function space
4787
virtual unsigned int max_local_dimension() const
4729
virtual std::size_t max_local_dimension() const
4788
4730
{
4789
4731
return 6;
4790
4732
}
4791
4733
4792
4734
/// Return the number of dofs on each cell facet
4793
virtual unsigned int num_facet_dofs() const
4735
virtual std::size_t num_facet_dofs() const
4794
4736
{
4795
4737
return 4;
4796
4738
}
4797
4739
4798
4740
/// Return the number of dofs associated with each cell entity of dimension d
4799
virtual unsigned int num_entity_dofs(unsigned int d) const
4741
virtual std::size_t num_entity_dofs(std::size_t d) const
4800
4742
{
4801
4743
switch (d)
4802
4744
{
4821
4763
}
4822
4764
4823
4765
/// Tabulate the local-to-global mapping of dofs on a cell
4824
virtual void tabulate_dofs(unsigned int* dofs,
4825
const ufc::mesh& m,
4766
virtual void tabulate_dofs(std::size_t* dofs,
4767
const std::vector<std::size_t>& num_global_entities,
4826
4768
const ufc::cell& c) const
4827
4769
{
4828
4770
unsigned int offset = 0;
4829
4771
dofs[0] = offset + c.entity_indices[0][0];
4830
4772
dofs[1] = offset + c.entity_indices[0][1];
4831
4773
dofs[2] = offset + c.entity_indices[0][2];
4832
offset += m.num_entities[0];
4774
offset += num_global_entities[0];
4833
4775
dofs[3] = offset + c.entity_indices[0][0];
4834
4776
dofs[4] = offset + c.entity_indices[0][1];
4835
4777
dofs[5] = offset + c.entity_indices[0][2];
4836
offset += m.num_entities[0];
4778
offset += num_global_entities[0];
4837
4779
}
4838
4780
4839
4781
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4840
virtual void tabulate_facet_dofs(unsigned int* dofs,
4841
unsigned int facet) const
4782
virtual void tabulate_facet_dofs(std::size_t* dofs,
4783
std::size_t facet) const
4842
4784
{
4843
4785
switch (facet)
4844
4786
{
4871
4813
}
4872
4814
4873
4815
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4874
virtual void tabulate_entity_dofs(unsigned int* dofs,
4875
unsigned int d, unsigned int i) const
4816
virtual void tabulate_entity_dofs(std::size_t* dofs,
4817
std::size_t d, std::size_t i) const
4876
4818
{
4877
4819
if (d > 2)
4878
4820
{
4927
4869
}
4928
4870
4929
4871
/// Tabulate the coordinates of all dofs on a cell
4930
virtual void tabulate_coordinates(double** coordinates,
4931
const ufc::cell& c) const
4872
virtual void tabulate_coordinates(double** dof_coordinates,
4873
const double* vertex_coordinates) const
4932
4874
{
4933
const double * const * x = c.coordinates;
4934
4935
coordinates[0][0] = x[0][0];
4936
coordinates[0][1] = x[0][1];
4937
coordinates[1][0] = x[1][0];
4938
coordinates[1][1] = x[1][1];
4939
coordinates[2][0] = x[2][0];
4940
coordinates[2][1] = x[2][1];
4941
coordinates[3][0] = x[0][0];
4942
coordinates[3][1] = x[0][1];
4943
coordinates[4][0] = x[1][0];
4944
coordinates[4][1] = x[1][1];
4945
coordinates[5][0] = x[2][0];
4946
coordinates[5][1] = x[2][1];
4875
dof_coordinates[0][0] = vertex_coordinates[0];
4876
dof_coordinates[0][1] = vertex_coordinates[1];
4877
dof_coordinates[1][0] = vertex_coordinates[2];
4878
dof_coordinates[1][1] = vertex_coordinates[3];
4879
dof_coordinates[2][0] = vertex_coordinates[4];
4880
dof_coordinates[2][1] = vertex_coordinates[5];
4881
dof_coordinates[3][0] = vertex_coordinates[0];
4882
dof_coordinates[3][1] = vertex_coordinates[1];
4883
dof_coordinates[4][0] = vertex_coordinates[2];
4884
dof_coordinates[4][1] = vertex_coordinates[3];
4885
dof_coordinates[5][0] = vertex_coordinates[4];
4886
dof_coordinates[5][1] = vertex_coordinates[5];
4947
4887
}
4948
4888
4949
4889
/// Return the number of sub dofmaps (for a mixed element)
4950
virtual unsigned int num_sub_dofmaps() const
4890
virtual std::size_t num_sub_dofmaps() const
4951
4891
{
4952
4892
return 2;
4953
4893
}
4954
4894
4955
4895
/// Create a new dofmap for sub dofmap i (for a mixed element)
4956
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
4896
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4957
4897
{
4958
4898
switch (i)
4959
4899
{
4985
4925
4986
4926
class stabilisedstokes_dofmap_2: public ufc::dofmap
4987
4927
{
4988
private:
4989
4990
unsigned int _global_dimension;
4991
4928
public:
4992
4929
4993
4930
/// Constructor
4994
4931
stabilisedstokes_dofmap_2() : ufc::dofmap()
4995
4932
{
4996
_global_dimension = 0;
4933
// Do nothing
4997
4934
}
4998
4935
4999
4936
/// Destructor
5005
4942
/// Return a string identifying the dofmap
5006
4943
virtual const char* signature() const
5007
4944
{
5008
return "FFC dofmap for MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) })";
4945
return "FFC dofmap for MixedElement(*[VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (3,) })";
5009
4946
}
5010
4947
5011
4948
/// Return true iff mesh entities of topological dimension d are needed
5012
virtual bool needs_mesh_entities(unsigned int d) const
4949
virtual bool needs_mesh_entities(std::size_t d) const
5013
4950
{
5014
4951
switch (d)
5015
4952
{
5033
4970
return false;
5034
4971
}
5035
4972
5036
/// Initialize dofmap for mesh (return true iff init_cell() is needed)
5037
virtual bool init_mesh(const ufc::mesh& m)
5038
{
5039
_global_dimension = 3*m.num_entities[0];
5040
return false;
5041
}
5042
5043
/// Initialize dofmap for given cell
5044
virtual void init_cell(const ufc::mesh& m,
5045
const ufc::cell& c)
5046
{
5047
// Do nothing
5048
}
5049
5050
/// Finish initialization of dofmap for cells
5051
virtual void init_cell_finalize()
5052
{
5053
// Do nothing
5054
}
5055
5056
4973
/// Return the topological dimension of the associated cell shape
5057
virtual unsigned int topological_dimension() const
4974
virtual std::size_t topological_dimension() const
5058
4975
{
5059
4976
return 2;
5060
4977
}
5061
4978
5062
4979
/// Return the geometric dimension of the associated cell shape
5063
virtual unsigned int geometric_dimension() const
4980
virtual std::size_t geometric_dimension() const
5064
4981
{
5065
4982
return 2;
5066
4983
}
5067
4984
5068
4985
/// Return the dimension of the global finite element function space
5069
virtual unsigned int global_dimension() const
4986
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4987
num_global_entities) const
5070
4988
{
5071
return _global_dimension;
4989
return 3*num_global_entities[0];
5072
4990
}
5073
4991
5074
4992
/// Return the dimension of the local finite element function space for a cell
5075
virtual unsigned int local_dimension(const ufc::cell& c) const
4993
virtual std::size_t local_dimension(const ufc::cell& c) const
5076
4994
{
5077
4995
return 9;
5078
4996
}
5079
4997
5080
4998
/// Return the maximum dimension of the local finite element function space
5081
virtual unsigned int max_local_dimension() const
4999
virtual std::size_t max_local_dimension() const
5082
5000
{
5083
5001
return 9;
5084
5002
}
5085
5003
5086
5004
/// Return the number of dofs on each cell facet
5087
virtual unsigned int num_facet_dofs() const
5005
virtual std::size_t num_facet_dofs() const
5088
5006
{
5089
5007
return 6;
5090
5008
}
5091
5009
5092
5010
/// Return the number of dofs associated with each cell entity of dimension d
5093
virtual unsigned int num_entity_dofs(unsigned int d) const
5011
virtual std::size_t num_entity_dofs(std::size_t d) const
5094
5012
{
5095
5013
switch (d)
5096
5014
{
5115
5033
}
5116
5034
5117
5035
/// Tabulate the local-to-global mapping of dofs on a cell
5118
virtual void tabulate_dofs(unsigned int* dofs,
5119
const ufc::mesh& m,
5036
virtual void tabulate_dofs(std::size_t* dofs,
5037
const std::vector<std::size_t>& num_global_entities,
5120
5038
const ufc::cell& c) const
5121
5039
{
5122
5040
unsigned int offset = 0;
5123
5041
dofs[0] = offset + c.entity_indices[0][0];
5124
5042
dofs[1] = offset + c.entity_indices[0][1];
5125
5043
dofs[2] = offset + c.entity_indices[0][2];
5126
offset += m.num_entities[0];
5044
offset += num_global_entities[0];
5127
5045
dofs[3] = offset + c.entity_indices[0][0];
5128
5046
dofs[4] = offset + c.entity_indices[0][1];
5129
5047
dofs[5] = offset + c.entity_indices[0][2];
5130
offset += m.num_entities[0];
5048
offset += num_global_entities[0];
5131
5049
dofs[6] = offset + c.entity_indices[0][0];
5132
5050
dofs[7] = offset + c.entity_indices[0][1];
5133
5051
dofs[8] = offset + c.entity_indices[0][2];
5134
offset += m.num_entities[0];
5052
offset += num_global_entities[0];
5135
5053
}
5136
5054
5137
5055
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
5138
virtual void tabulate_facet_dofs(unsigned int* dofs,
5139
unsigned int facet) const
5056
virtual void tabulate_facet_dofs(std::size_t* dofs,
5057
std::size_t facet) const
5140
5058
{
5141
5059
switch (facet)
5142
5060
{
5175
5093
}
5176
5094
5177
5095
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
5178
virtual void tabulate_entity_dofs(unsigned int* dofs,
5179
unsigned int d, unsigned int i) const
5096
virtual void tabulate_entity_dofs(std::size_t* dofs,
5097
std::size_t d, std::size_t i) const
5180
5098
{
5181
5099
if (d > 2)
5182
5100
{
5234
5152
}
5235
5153
5236
5154
/// Tabulate the coordinates of all dofs on a cell
5237
virtual void tabulate_coordinates(double** coordinates,
5238
const ufc::cell& c) const
5155
virtual void tabulate_coordinates(double** dof_coordinates,
5156
const double* vertex_coordinates) const
5239
5157
{
5240
const double * const * x = c.coordinates;
5241
5242
coordinates[0][0] = x[0][0];
5243
coordinates[0][1] = x[0][1];
5244
coordinates[1][0] = x[1][0];
5245
coordinates[1][1] = x[1][1];
5246
coordinates[2][0] = x[2][0];
5247
coordinates[2][1] = x[2][1];
5248
coordinates[3][0] = x[0][0];
5249
coordinates[3][1] = x[0][1];
5250
coordinates[4][0] = x[1][0];
5251
coordinates[4][1] = x[1][1];
5252
coordinates[5][0] = x[2][0];
5253
coordinates[5][1] = x[2][1];
5254
coordinates[6][0] = x[0][0];
5255
coordinates[6][1] = x[0][1];
5256
coordinates[7][0] = x[1][0];
5257
coordinates[7][1] = x[1][1];
5258
coordinates[8][0] = x[2][0];
5259
coordinates[8][1] = x[2][1];
5158
dof_coordinates[0][0] = vertex_coordinates[0];
5159
dof_coordinates[0][1] = vertex_coordinates[1];
5160
dof_coordinates[1][0] = vertex_coordinates[2];
5161
dof_coordinates[1][1] = vertex_coordinates[3];
5162
dof_coordinates[2][0] = vertex_coordinates[4];
5163
dof_coordinates[2][1] = vertex_coordinates[5];
5164
dof_coordinates[3][0] = vertex_coordinates[0];
5165
dof_coordinates[3][1] = vertex_coordinates[1];
5166
dof_coordinates[4][0] = vertex_coordinates[2];
5167
dof_coordinates[4][1] = vertex_coordinates[3];
5168
dof_coordinates[5][0] = vertex_coordinates[4];
5169
dof_coordinates[5][1] = vertex_coordinates[5];
5170
dof_coordinates[6][0] = vertex_coordinates[0];
5171
dof_coordinates[6][1] = vertex_coordinates[1];
5172
dof_coordinates[7][0] = vertex_coordinates[2];
5173
dof_coordinates[7][1] = vertex_coordinates[3];
5174
dof_coordinates[8][0] = vertex_coordinates[4];
5175
dof_coordinates[8][1] = vertex_coordinates[5];
5260
5176
}
5261
5177
5262
5178
/// Return the number of sub dofmaps (for a mixed element)
5263
virtual unsigned int num_sub_dofmaps() const
5179
virtual std::size_t num_sub_dofmaps() const
5264
5180
{
5265
5181
return 2;
5266
5182
}
5267
5183
5268
5184
/// Create a new dofmap for sub dofmap i (for a mixed element)
5269
virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
5185
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
5270
5186
{
5271
5187
switch (i)
5272
5188
{
5297
5213
/// tensor corresponding to the local contribution to a form from
5298
5214
/// the integral over a cell.
5299
5215
5300
class stabilisedstokes_cell_integral_0_0: public ufc::cell_integral
5216
class stabilisedstokes_cell_integral_0_otherwise: public ufc::cell_integral
5301
5217
{
5302
5218
public:
5303
5219
5304
5220
/// Constructor
5305
stabilisedstokes_cell_integral_0_0() : ufc::cell_integral()
5221
stabilisedstokes_cell_integral_0_otherwise() : ufc::cell_integral()
5306
5222
{
5307
5223
// Do nothing
5308
5224
}
5309
5225
5310
5226
/// Destructor
5311
virtual ~stabilisedstokes_cell_integral_0_0()
5227
virtual ~stabilisedstokes_cell_integral_0_otherwise()
5312
5228
{
5313
5229
// Do nothing
5314
5230
}
5316
5232
/// Tabulate the tensor for the contribution from a local cell
5317
5233
virtual void tabulate_tensor(double* A,
5318
5234
const double * const * w,
5319
const ufc::cell& c) const
5235
const double* vertex_coordinates,
5236
int cell_orientation) const
5320
5237
{
5321
// Extract vertex coordinates
5322
const double * const * x = c.coordinates;
5323
5324
// Compute Jacobian of affine map from reference cell
5325
const double J_00 = x[1][0] - x[0][0];
5326
const double J_01 = x[2][0] - x[0][0];
5327
const double J_10 = x[1][1] - x[0][1];
5328
const double J_11 = x[2][1] - x[0][1];
5329
5330
// Compute determinant of Jacobian
5331
double detJ = J_00*J_11 - J_01*J_10;
5332
5333
// Compute inverse of Jacobian
5334
const double K_00 = J_11 / detJ;
5335
const double K_01 = -J_01 / detJ;
5336
const double K_10 = -J_10 / detJ;
5337
const double K_11 = J_00 / detJ;
5238
// Compute Jacobian
5239
double J[4];
5240
compute_jacobian_triangle_2d(J, vertex_coordinates);
5241
5242
// Compute Jacobian inverse and determinant
5243
double K[4];
5244
double detJ;
5245
compute_jacobian_inverse_triangle_2d(K, detJ, J);
5338
5246
5339
5247
// Set scale factor
5340
5248
const double det = std::abs(detJ);
5341
5249
5342
// Cell Volume.
5343
5344
// Compute circumradius, assuming triangle is embedded in 2D.
5345
5346
5347
// Facet Area.
5250
// Cell volume
5251
5252
// Compute circumradius of triangle in 2D
5253
5254
5255
// Facet area
5348
5256
5349
5257
// Array of quadrature weights.
5350
5258
static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
5420
5328
for (unsigned int k = 0; k < 9; k++)
5421
5329
{
5422
5330
// Number of operations to compute entry: 72
5423
A[j*9 + k] += (((((K_01*FE1_C2_D10[ip][j] + K_11*FE1_C2_D01[ip][j]))*((K_01*FE1_C2_D10[ip][k] + K_11*FE1_C2_D01[ip][k])) + ((K_00*FE1_C2_D10[ip][j] + K_10*FE1_C2_D01[ip][j]))*((K_00*FE1_C2_D10[ip][k] + K_10*FE1_C2_D01[ip][k]))))*F0*0.2*F0 + ((((K_00*FE1_C0_D10[ip][k] + K_10*FE1_C0_D01[ip][k]) + (K_01*FE1_C1_D10[ip][k] + K_11*FE1_C1_D01[ip][k])))*FE1_C2[ip][j] + (((((K_00*FE1_C0_D10[ip][j] + K_10*FE1_C0_D01[ip][j]) + (K_01*FE1_C1_D10[ip][j] + K_11*FE1_C1_D01[ip][j])))*FE1_C2[ip][k])*(-1.0) + ((((K_01*FE1_C1_D10[ip][j] + K_11*FE1_C1_D01[ip][j]))*((K_01*FE1_C1_D10[ip][k] + K_11*FE1_C1_D01[ip][k])) + ((K_01*FE1_C0_D10[ip][j] + K_11*FE1_C0_D01[ip][j]))*((K_01*FE1_C0_D10[ip][k] + K_11*FE1_C0_D01[ip][k]))) + (((K_00*FE1_C1_D10[ip][j] + K_10*FE1_C1_D01[ip][j]))*((K_00*FE1_C1_D10[ip][k] + K_10*FE1_C1_D01[ip][k])) + ((K_00*FE1_C0_D10[ip][j] + K_10*FE1_C0_D01[ip][j]))*((K_00*FE1_C0_D10[ip][k] + K_10*FE1_C0_D01[ip][k])))))))*W3[ip]*det;
5331
A[j*9 + k] += (((((((K[1]*FE1_C1_D10[ip][j] + K[3]*FE1_C1_D01[ip][j]) + (K[0]*FE1_C0_D10[ip][j] + K[2]*FE1_C0_D01[ip][j])))*FE1_C2[ip][k])*(-1.0) + ((((K[1]*FE1_C1_D10[ip][j] + K[3]*FE1_C1_D01[ip][j]))*((K[1]*FE1_C1_D10[ip][k] + K[3]*FE1_C1_D01[ip][k])) + ((K[1]*FE1_C0_D10[ip][j] + K[3]*FE1_C0_D01[ip][j]))*((K[1]*FE1_C0_D10[ip][k] + K[3]*FE1_C0_D01[ip][k]))) + (((K[0]*FE1_C1_D10[ip][j] + K[2]*FE1_C1_D01[ip][j]))*((K[0]*FE1_C1_D10[ip][k] + K[2]*FE1_C1_D01[ip][k])) + ((K[0]*FE1_C0_D10[ip][j] + K[2]*FE1_C0_D01[ip][j]))*((K[0]*FE1_C0_D10[ip][k] + K[2]*FE1_C0_D01[ip][k]))))) + (((K[0]*FE1_C0_D10[ip][k] + K[2]*FE1_C0_D01[ip][k]) + (K[1]*FE1_C1_D10[ip][k] + K[3]*FE1_C1_D01[ip][k])))*FE1_C2[ip][j]) + ((((K[1]*FE1_C2_D10[ip][j] + K[3]*FE1_C2_D01[ip][j]))*((K[1]*FE1_C2_D10[ip][k] + K[3]*FE1_C2_D01[ip][k])) + ((K[0]*FE1_C2_D10[ip][j] + K[2]*FE1_C2_D01[ip][j]))*((K[0]*FE1_C2_D10[ip][k] + K[2]*FE1_C2_D01[ip][k]))))*F0*0.2*F0)*W3[ip]*det;
5424
5332
}// end loop over 'k'
5425
5333
}// end loop over 'j'
5426
5334
}// end loop over 'ip'
5427
5335
}
5428
5336
5429
/// Tabulate the tensor for the contribution from a local cell
5430
/// using the specified reference cell quadrature points/weights
5431
virtual void tabulate_tensor(double* A,
5432
const double * const * w,
5433
const ufc::cell& c,
5434
unsigned int num_quadrature_points,
5435
const double * const * quadrature_points,
5436
const double* quadrature_weights) const
5437
{
5438
std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
5439
}
5440
5441
5337
};
5442
5338
5443
5339
/// This class defines the interface for the tabulation of the cell
5444
5340
/// tensor corresponding to the local contribution to a form from
5445
5341
/// the integral over a cell.
5446
5342
5447
class stabilisedstokes_cell_integral_1_0: public ufc::cell_integral
5343
class stabilisedstokes_cell_integral_1_otherwise: public ufc::cell_integral
5448
5344
{
5449
5345
public:
5450
5346
5451
5347
/// Constructor
5452
stabilisedstokes_cell_integral_1_0() : ufc::cell_integral()
5348
stabilisedstokes_cell_integral_1_otherwise() : ufc::cell_integral()
5453
5349
{
5454
5350
// Do nothing
5455
5351
}
5456
5352
5457
5353
/// Destructor
5458
virtual ~stabilisedstokes_cell_integral_1_0()
5354
virtual ~stabilisedstokes_cell_integral_1_otherwise()
5459
5355
{
5460
5356
// Do nothing
5461
5357
}
5463
5359
/// Tabulate the tensor for the contribution from a local cell
5464
5360
virtual void tabulate_tensor(double* A,
5465
5361
const double * const * w,
5466
const ufc::cell& c) const
5362
const double* vertex_coordinates,
5363
int cell_orientation) const
5467
5364
{
5468
// Extract vertex coordinates
5469
const double * const * x = c.coordinates;
5470
5471
// Compute Jacobian of affine map from reference cell
5472
const double J_00 = x[1][0] - x[0][0];
5473
const double J_01 = x[2][0] - x[0][0];
5474
const double J_10 = x[1][1] - x[0][1];
5475
const double J_11 = x[2][1] - x[0][1];
5476
5477
// Compute determinant of Jacobian
5478
double detJ = J_00*J_11 - J_01*J_10;
5479
5480
// Compute inverse of Jacobian
5481
const double K_00 = J_11 / detJ;
5482
const double K_01 = -J_01 / detJ;
5483
const double K_10 = -J_10 / detJ;
5484
const double K_11 = J_00 / detJ;
5365
// Compute Jacobian
5366
double J[4];
5367
compute_jacobian_triangle_2d(J, vertex_coordinates);
5368
5369
// Compute Jacobian inverse and determinant
5370
double K[4];
5371
double detJ;
5372
compute_jacobian_inverse_triangle_2d(K, detJ, J);
5485
5373
5486
5374
// Set scale factor
5487
5375
const double det = std::abs(detJ);
5488
5376
5489
// Cell Volume.
5490
5491
// Compute circumradius, assuming triangle is embedded in 2D.
5492
5493
5494
// Facet Area.
5377
// Cell volume
5378
5379
// Compute circumradius of triangle in 2D
5380
5381
5382
// Facet area
5495
5383
5496
5384
// Array of quadrature weights.
5497
5385
static const double W6[6] = {0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333, 0.083333333};
5590
5478
for (unsigned int j = 0; j < 9; j++)
5591
5479
{
5592
5480
// Number of operations to compute entry: 20
5593
A[j] += (((FE2_C1[ip][j] + ((K_01*FE2_C2_D10[ip][j] + K_11*FE2_C2_D01[ip][j]))*F1*0.2*F1))*F2 + ((((K_00*FE2_C2_D10[ip][j] + K_10*FE2_C2_D01[ip][j]))*F1*0.2*F1 + FE2_C0[ip][j]))*F0)*W6[ip]*det;
5481
A[j] += (((((K[0]*FE2_C2_D10[ip][j] + K[2]*FE2_C2_D01[ip][j]))*F1*0.2*F1 + FE2_C0[ip][j]))*F0 + ((FE2_C1[ip][j] + ((K[1]*FE2_C2_D10[ip][j] + K[3]*FE2_C2_D01[ip][j]))*F1*0.2*F1))*F2)*W6[ip]*det;
5594
5482
}// end loop over 'j'
5595
5483
}// end loop over 'ip'
5596
5484
}
5597
5485
5598
/// Tabulate the tensor for the contribution from a local cell
5599
/// using the specified reference cell quadrature points/weights
5600
virtual void tabulate_tensor(double* A,
5601
const double * const * w,
5602
const ufc::cell& c,
5603
unsigned int num_quadrature_points,
5604
const double * const * quadrature_points,
5605
const double* quadrature_weights) const
5606
{
5607
std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
5608
}
5609
5610
5486
};
5611
5487
5612
5488
/// This class defines the interface for the assembly of the global
5643
5519
/// Return a string identifying the form
5644
5520
virtual const char* signature() const
5645
5521
{
5646
return "Form([Integral(Sum(Product(IndexSum(Product(Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((FixedIndex(2),), {})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((FixedIndex(2),), {})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0), Product(FloatValue(0.2, (), (), {}), Coefficient(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0)))), Sum(Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((Index(3),), {Index(3): 2})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(2),), {}))), Sum(IndexSum(IndexSum(Product(Indexed(ComponentTensor(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(4),), {Index(4): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(4),), {Index(4): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(5),), {Index(5): 2})), MultiIndex((Index(5), Index(4)), {Index(4): 2, Index(5): 2})), MultiIndex((Index(6), Index(7)), {Index(7): 2, Index(6): 2})), Indexed(ComponentTensor(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(8),), {Index(8): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((Index(8),), {Index(8): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(9),), {Index(9): 2})), MultiIndex((Index(9), Index(8)), {Index(8): 2, Index(9): 2})), MultiIndex((Index(6), Index(7)), {Index(7): 2, Index(6): 2}))), MultiIndex((Index(6),), {Index(6): 2})), MultiIndex((Index(7),), {Index(7): 2})), Product(IntValue(-1, (), (), {}), Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((Index(10),), {Index(10): 2})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 1), MultiIndex((FixedIndex(2),), {}))))))), Measure('cell', 0, None))])";
5522
return "f60265baac0db00f7c1f0f27eebe47270c7ba139cd94da79afe0ad3e65f48e700510e32ef2f7b053ba9c6069c0035dad185c278323f066e40944f4e009cdcf21";
5647
5523
}
5648
5524
5649
5525
/// Return the rank of the global tensor (r)
5650
virtual unsigned int rank() const
5526
virtual std::size_t rank() const
5651
5527
{
5652
5528
return 2;
5653
5529
}
5654
5530
5655
5531
/// Return the number of coefficients (n)
5656
virtual unsigned int num_coefficients() const
5532
virtual std::size_t num_coefficients() const
5657
5533
{
5658
5534
return 1;
5659
5535
}
5660
5536
5661
5537
/// Return the number of cell domains
5662
virtual unsigned int num_cell_domains() const
5538
virtual std::size_t num_cell_domains() const
5663
5539
{
5664
return 1;
5540
return 0;
5665
5541
}
5666
5542
5667
5543
/// Return the number of exterior facet domains
5668
virtual unsigned int num_exterior_facet_domains() const
5544
virtual std::size_t num_exterior_facet_domains() const
5669
5545
{
5670
5546
return 0;
5671
5547
}
5672
5548
5673
5549
/// Return the number of interior facet domains
5674
virtual unsigned int num_interior_facet_domains() const
5675
{
5676
return 0;
5550
virtual std::size_t num_interior_facet_domains() const
5551
{
5552
return 0;
5553
}
5554
5555
/// Return the number of point domains
5556
virtual std::size_t num_point_domains() const
5557
{
5558
return 0;
5559
}
5560
5561
/// Return whether the form has any cell integrals
5562
virtual bool has_cell_integrals() const
5563
{
5564
return true;
5565
}
5566
5567
/// Return whether the form has any exterior facet integrals
5568
virtual bool has_exterior_facet_integrals() const
5569
{
5570
return false;
5571
}
5572
5573
/// Return whether the form has any interior facet integrals
5574
virtual bool has_interior_facet_integrals() const
5575
{
5576
return false;
5577
}
5578
5579
/// Return whether the form has any point integrals
5580
virtual bool has_point_integrals() const
5581
{
5582
return false;
5677
5583
}
5678
5584
5679
5585
/// Create a new finite element for argument function i
5680
virtual ufc::finite_element* create_finite_element(unsigned int i) const
5586
virtual ufc::finite_element* create_finite_element(std::size_t i) const
5681
5587
{
5682
5588
switch (i)
5683
5589
{
5702
5608
}
5703
5609
5704
5610
/// Create a new dofmap for argument function i
5705
virtual ufc::dofmap* create_dofmap(unsigned int i) const
5611
virtual ufc::dofmap* create_dofmap(std::size_t i) const
5706
5612
{
5707
5613
switch (i)
5708
5614
{
5727
5633
}
5728
5634
5729
5635
/// Create a new cell integral on sub domain i
5730
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
5636
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
5731
5637
{
5732
switch (i)
5733
{
5734
case 0:
5735
{
5736
return new stabilisedstokes_cell_integral_0_0();
5737
break;
5738
}
5739
}
5740
5741
5638
return 0;
5742
5639
}
5743
5640
5744
5641
/// Create a new exterior facet integral on sub domain i
5745
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
5642
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
5746
5643
{
5747
5644
return 0;
5748
5645
}
5749
5646
5750
5647
/// Create a new interior facet integral on sub domain i
5751
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
5648
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
5649
{
5650
return 0;
5651
}
5652
5653
/// Create a new point integral on sub domain i
5654
virtual ufc::point_integral* create_point_integral(std::size_t i) const
5655
{
5656
return 0;
5657
}
5658
5659
/// Create a new cell integral on everywhere else
5660
virtual ufc::cell_integral* create_default_cell_integral() const
5661
{
5662
return new stabilisedstokes_cell_integral_0_otherwise();
5663
}
5664
5665
/// Create a new exterior facet integral on everywhere else
5666
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
5667
{
5668
return 0;
5669
}
5670
5671
/// Create a new interior facet integral on everywhere else
5672
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
5673
{
5674
return 0;
5675
}
5676
5677
/// Create a new point integral on everywhere else
5678
virtual ufc::point_integral* create_default_point_integral() const
5752
5679
{
5753
5680
return 0;
5754
5681
}
5789
5716
/// Return a string identifying the form
5790
5717
virtual const char* signature() const
5791
5718
{
5792
return "Form([Integral(IndexSum(Product(Indexed(Coefficient(VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), 0), MultiIndex((Index(0),), {Index(0): 2})), Indexed(Sum(ComponentTensor(Product(Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((FixedIndex(2),), {})), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((Index(2),), {Index(2): 2})), Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 1), Product(FloatValue(0.2, (), (), {}), Coefficient(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 1)))), MultiIndex((Index(2),), {Index(2): 2})), ListTensor(Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(0),), {})), Indexed(Argument(MixedElement(*[VectorElement('Lagrange', Cell('triangle', Space(2)), 1, 2, None), FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(1),), {})))), MultiIndex((Index(0),), {Index(0): 2}))), MultiIndex((Index(0),), {Index(0): 2})), Measure('cell', 0, None))])";
5719
return "f191f10a1ec4b9dfd5fd0fb1f0bf10c58b799d565995c8a739d80cf6853d9218abaf909aa5d4c1e1b0ca79350344217265d273fa9e12df4e877b433840b206d2";
5793
5720
}
5794
5721
5795
5722
/// Return the rank of the global tensor (r)
5796
virtual unsigned int rank() const
5723
virtual std::size_t rank() const
5797
5724
{
5798
5725
return 1;
5799
5726
}
5800
5727
5801
5728
/// Return the number of coefficients (n)
5802
virtual unsigned int num_coefficients() const
5729
virtual std::size_t num_coefficients() const
5803
5730
{
5804
5731
return 2;
5805
5732
}
5806
5733
5807
5734
/// Return the number of cell domains
5808
virtual unsigned int num_cell_domains() const
5735
virtual std::size_t num_cell_domains() const
5809
5736
{
5810
return 1;
5737
return 0;
5811
5738
}
5812
5739
5813
5740
/// Return the number of exterior facet domains
5814
virtual unsigned int num_exterior_facet_domains() const
5741
virtual std::size_t num_exterior_facet_domains() const
5815
5742
{
5816
5743
return 0;
5817
5744
}
5818
5745
5819
5746
/// Return the number of interior facet domains
5820
virtual unsigned int num_interior_facet_domains() const
5821
{
5822
return 0;
5747
virtual std::size_t num_interior_facet_domains() const
5748
{
5749
return 0;
5750
}
5751
5752
/// Return the number of point domains
5753
virtual std::size_t num_point_domains() const
5754
{
5755
return 0;
5756
}
5757
5758
/// Return whether the form has any cell integrals
5759
virtual bool has_cell_integrals() const
5760
{
5761
return true;
5762
}
5763
5764
/// Return whether the form has any exterior facet integrals
5765
virtual bool has_exterior_facet_integrals() const
5766
{
5767
return false;
5768
}
5769
5770
/// Return whether the form has any interior facet integrals
5771
virtual bool has_interior_facet_integrals() const
5772
{
5773
return false;
5774
}
5775
5776
/// Return whether the form has any point integrals
5777
virtual bool has_point_integrals() const
5778
{
5779
return false;
5823
5780
}
5824
5781
5825
5782
/// Create a new finite element for argument function i
5826
virtual ufc::finite_element* create_finite_element(unsigned int i) const
5783
virtual ufc::finite_element* create_finite_element(std::size_t i) const
5827
5784
{
5828
5785
switch (i)
5829
5786
{
5848
5805
}
5849
5806
5850
5807
/// Create a new dofmap for argument function i
5851
virtual ufc::dofmap* create_dofmap(unsigned int i) const
5808
virtual ufc::dofmap* create_dofmap(std::size_t i) const
5852
5809
{
5853
5810
switch (i)
5854
5811
{
5873
5830
}
5874
5831
5875
5832
/// Create a new cell integral on sub domain i
5876
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
5833
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
5877
5834
{
5878
switch (i)
5879
{
5880
case 0:
5881
{
5882
return new stabilisedstokes_cell_integral_1_0();
5883
break;
5884
}
5885
}
5886
5887
5835
return 0;
5888
5836
}
5889
5837
5890
5838
/// Create a new exterior facet integral on sub domain i
5891
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
5839
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
5892
5840
{
5893
5841
return 0;
5894
5842
}
5895
5843
5896
5844
/// Create a new interior facet integral on sub domain i
5897
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
5845
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
5846
{
5847
return 0;
5848
}
5849
5850
/// Create a new point integral on sub domain i
5851
virtual ufc::point_integral* create_point_integral(std::size_t i) const
5852
{
5853
return 0;
5854
}
5855
5856
/// Create a new cell integral on everywhere else
5857
virtual ufc::cell_integral* create_default_cell_integral() const
5858
{
5859
return new stabilisedstokes_cell_integral_1_otherwise();
5860
}
5861
5862
/// Create a new exterior facet integral on everywhere else
5863
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
5864
{
5865
return 0;
5866
}
5867
5868
/// Create a new interior facet integral on everywhere else
5869
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
5870
{
5871
return 0;
5872
}
5873
5874
/// Create a new point integral on everywhere else
5875
virtual ufc::point_integral* create_default_point_integral() const
5898
5876
{
5899
5877
return 0;
5900
5878
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1