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