~ubuntu-branches/ubuntu/wily/ffc/wily
» Revision
17
: test/regression/references/r_auto/PointMeasure.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_auto/PointMeasure.h
- browse files at revision 17
- compare with another revision
- download diff
- download tarball
- view history from revision 17
- Committer: Package Import Robot
- Author(s): Johannes Ring
- Date: 2013-06-26 14:48:32 UTC
- mfrom: (1.1.13)
- Revision ID: package-import@ubuntu.com-20130626144832-1xd8htax4s3utybz
Tags: 1.2.0-1
* New upstream release.
* debian/control:
- Bump required version for python-ufc, python-fiat, python-instant
and python-ufl in Depends field.
- Bump Standards-Version to 3.9.4.
- Remove DM-Upload-Allowed field.
- Bump required debhelper version in Build-Depends.
- Remove cdbs from Build-Depends.
- Use canonical URIs for Vcs-* fields.
* debian/compat: Bump to compatibility level 9.
* debian/rules: Rewrite for debhelper (drop cdbs).
* debian/control:
- Bump required version for python-ufc, python-fiat, python-instant
and python-ufl in Depends field.
- Bump Standards-Version to 3.9.4.
- Remove DM-Upload-Allowed field.
- Bump required debhelper version in Build-Depends.
- Remove cdbs from Build-Depends.
- Use canonical URIs for Vcs-* fields.
* debian/compat: Bump to compatibility level 9.
* debian/rules: Rewrite for debhelper (drop cdbs).
- files added:
- .pc
- ffc/uflacsrepr
- sandbox/features
- sandbox/multidomain
- files modified:
- ChangeLog
added
removed
test/regression/references/r_auto/PointMeasure.h
1
// This code conforms with the UFC specification version 2.2.0
2
// and was automatically generated by FFC version 1.2.0.
3
//
4
// This code was generated with the following parameters:
5
//
6
// cache_dir: ''
7
// convert_exceptions_to_warnings: True
8
// cpp_optimize: False
9
// cpp_optimize_flags: '-O2'
10
// epsilon: 1e-14
11
// error_control: False
12
// form_postfix: True
13
// format: 'ufc'
14
// log_level: 20
15
// log_prefix: ''
16
// optimize: False
17
// output_dir: '.'
18
// precision: '8'
19
// quadrature_degree: 'auto'
20
// quadrature_rule: 'auto'
21
// representation: 'auto'
22
// split: False
23
// swig_binary: 'swig'
24
// swig_path: ''
25
26
#ifndef __POINTMEASURE_H
27
#define __POINTMEASURE_H
28
29
#include <cmath>
30
#include <stdexcept>
31
#include <fstream>
32
#include <ufc.h>
33
34
/// This class defines the interface for a finite element.
35
36
class pointmeasure_finite_element_0: public ufc::finite_element
37
{
38
public:
39
40
/// Constructor
41
pointmeasure_finite_element_0() : ufc::finite_element()
42
{
43
// Do nothing
44
}
45
46
/// Destructor
47
virtual ~pointmeasure_finite_element_0()
48
{
49
// Do nothing
50
}
51
52
/// Return a string identifying the finite element
53
virtual const char* signature() const
54
{
55
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
56
}
57
58
/// Return the cell shape
59
virtual ufc::shape cell_shape() const
60
{
61
return ufc::triangle;
62
}
63
64
/// Return the topological dimension of the cell shape
65
virtual std::size_t topological_dimension() const
66
{
67
return 2;
68
}
69
70
/// Return the geometric dimension of the cell shape
71
virtual std::size_t geometric_dimension() const
72
{
73
return 2;
74
}
75
76
/// Return the dimension of the finite element function space
77
virtual std::size_t space_dimension() const
78
{
79
return 6;
80
}
81
82
/// Return the rank of the value space
83
virtual std::size_t value_rank() const
84
{
85
return 0;
86
}
87
88
/// Return the dimension of the value space for axis i
89
virtual std::size_t value_dimension(std::size_t i) const
90
{
91
return 1;
92
}
93
94
/// Evaluate basis function i at given point x in cell
95
virtual void evaluate_basis(std::size_t i,
96
double* values,
97
const double* x,
98
const double* vertex_coordinates,
99
int cell_orientation) const
100
{
101
// Compute Jacobian
102
double J[4];
103
compute_jacobian_triangle_2d(J, vertex_coordinates);
104
105
// Compute Jacobian inverse and determinant
106
double K[4];
107
double detJ;
108
compute_jacobian_inverse_triangle_2d(K, detJ, J);
109
110
111
// Compute constants
112
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
113
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
114
115
// Get coordinates and map to the reference (FIAT) element
116
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
117
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
118
119
// Reset values
120
*values = 0.0;
121
switch (i)
122
{
123
case 0:
124
{
125
126
// Array of basisvalues
127
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
128
129
// Declare helper variables
130
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
131
double tmp1 = (1.0 - Y)/2.0;
132
double tmp2 = tmp1*tmp1;
133
134
// Compute basisvalues
135
basisvalues[0] = 1.0;
136
basisvalues[1] = tmp0;
137
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
138
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
139
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
140
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
141
basisvalues[0] *= std::sqrt(0.5);
142
basisvalues[2] *= std::sqrt(1.0);
143
basisvalues[5] *= std::sqrt(1.5);
144
basisvalues[1] *= std::sqrt(3.0);
145
basisvalues[4] *= std::sqrt(4.5);
146
basisvalues[3] *= std::sqrt(7.5);
147
148
// Table(s) of coefficients
149
static const double coefficients0[6] = \
150
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
151
152
// Compute value(s)
153
for (unsigned int r = 0; r < 6; r++)
154
{
155
*values += coefficients0[r]*basisvalues[r];
156
}// end loop over 'r'
157
break;
158
}
159
case 1:
160
{
161
162
// Array of basisvalues
163
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
164
165
// Declare helper variables
166
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
167
double tmp1 = (1.0 - Y)/2.0;
168
double tmp2 = tmp1*tmp1;
169
170
// Compute basisvalues
171
basisvalues[0] = 1.0;
172
basisvalues[1] = tmp0;
173
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
174
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
175
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
176
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
177
basisvalues[0] *= std::sqrt(0.5);
178
basisvalues[2] *= std::sqrt(1.0);
179
basisvalues[5] *= std::sqrt(1.5);
180
basisvalues[1] *= std::sqrt(3.0);
181
basisvalues[4] *= std::sqrt(4.5);
182
basisvalues[3] *= std::sqrt(7.5);
183
184
// Table(s) of coefficients
185
static const double coefficients0[6] = \
186
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
187
188
// Compute value(s)
189
for (unsigned int r = 0; r < 6; r++)
190
{
191
*values += coefficients0[r]*basisvalues[r];
192
}// end loop over 'r'
193
break;
194
}
195
case 2:
196
{
197
198
// Array of basisvalues
199
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
200
201
// Declare helper variables
202
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
203
double tmp1 = (1.0 - Y)/2.0;
204
double tmp2 = tmp1*tmp1;
205
206
// Compute basisvalues
207
basisvalues[0] = 1.0;
208
basisvalues[1] = tmp0;
209
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
210
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
211
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
212
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
213
basisvalues[0] *= std::sqrt(0.5);
214
basisvalues[2] *= std::sqrt(1.0);
215
basisvalues[5] *= std::sqrt(1.5);
216
basisvalues[1] *= std::sqrt(3.0);
217
basisvalues[4] *= std::sqrt(4.5);
218
basisvalues[3] *= std::sqrt(7.5);
219
220
// Table(s) of coefficients
221
static const double coefficients0[6] = \
222
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
223
224
// Compute value(s)
225
for (unsigned int r = 0; r < 6; r++)
226
{
227
*values += coefficients0[r]*basisvalues[r];
228
}// end loop over 'r'
229
break;
230
}
231
case 3:
232
{
233
234
// Array of basisvalues
235
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
236
237
// Declare helper variables
238
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
239
double tmp1 = (1.0 - Y)/2.0;
240
double tmp2 = tmp1*tmp1;
241
242
// Compute basisvalues
243
basisvalues[0] = 1.0;
244
basisvalues[1] = tmp0;
245
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
246
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
247
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
248
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
249
basisvalues[0] *= std::sqrt(0.5);
250
basisvalues[2] *= std::sqrt(1.0);
251
basisvalues[5] *= std::sqrt(1.5);
252
basisvalues[1] *= std::sqrt(3.0);
253
basisvalues[4] *= std::sqrt(4.5);
254
basisvalues[3] *= std::sqrt(7.5);
255
256
// Table(s) of coefficients
257
static const double coefficients0[6] = \
258
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
259
260
// Compute value(s)
261
for (unsigned int r = 0; r < 6; r++)
262
{
263
*values += coefficients0[r]*basisvalues[r];
264
}// end loop over 'r'
265
break;
266
}
267
case 4:
268
{
269
270
// Array of basisvalues
271
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
272
273
// Declare helper variables
274
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
275
double tmp1 = (1.0 - Y)/2.0;
276
double tmp2 = tmp1*tmp1;
277
278
// Compute basisvalues
279
basisvalues[0] = 1.0;
280
basisvalues[1] = tmp0;
281
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
282
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
283
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
284
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
285
basisvalues[0] *= std::sqrt(0.5);
286
basisvalues[2] *= std::sqrt(1.0);
287
basisvalues[5] *= std::sqrt(1.5);
288
basisvalues[1] *= std::sqrt(3.0);
289
basisvalues[4] *= std::sqrt(4.5);
290
basisvalues[3] *= std::sqrt(7.5);
291
292
// Table(s) of coefficients
293
static const double coefficients0[6] = \
294
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
295
296
// Compute value(s)
297
for (unsigned int r = 0; r < 6; r++)
298
{
299
*values += coefficients0[r]*basisvalues[r];
300
}// end loop over 'r'
301
break;
302
}
303
case 5:
304
{
305
306
// Array of basisvalues
307
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
308
309
// Declare helper variables
310
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
311
double tmp1 = (1.0 - Y)/2.0;
312
double tmp2 = tmp1*tmp1;
313
314
// Compute basisvalues
315
basisvalues[0] = 1.0;
316
basisvalues[1] = tmp0;
317
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
318
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
319
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
320
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
321
basisvalues[0] *= std::sqrt(0.5);
322
basisvalues[2] *= std::sqrt(1.0);
323
basisvalues[5] *= std::sqrt(1.5);
324
basisvalues[1] *= std::sqrt(3.0);
325
basisvalues[4] *= std::sqrt(4.5);
326
basisvalues[3] *= std::sqrt(7.5);
327
328
// Table(s) of coefficients
329
static const double coefficients0[6] = \
330
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
331
332
// Compute value(s)
333
for (unsigned int r = 0; r < 6; r++)
334
{
335
*values += coefficients0[r]*basisvalues[r];
336
}// end loop over 'r'
337
break;
338
}
339
}
340
341
}
342
343
/// Evaluate all basis functions at given point x in cell
344
virtual void evaluate_basis_all(double* values,
345
const double* x,
346
const double* vertex_coordinates,
347
int cell_orientation) const
348
{
349
// Helper variable to hold values of a single dof.
350
double dof_values = 0.0;
351
352
// Loop dofs and call evaluate_basis
353
for (unsigned int r = 0; r < 6; r++)
354
{
355
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
356
values[r] = dof_values;
357
}// end loop over 'r'
358
}
359
360
/// Evaluate order n derivatives of basis function i at given point x in cell
361
virtual void evaluate_basis_derivatives(std::size_t i,
362
std::size_t n,
363
double* values,
364
const double* x,
365
const double* vertex_coordinates,
366
int cell_orientation) const
367
{
368
// Compute Jacobian
369
double J[4];
370
compute_jacobian_triangle_2d(J, vertex_coordinates);
371
372
// Compute Jacobian inverse and determinant
373
double K[4];
374
double detJ;
375
compute_jacobian_inverse_triangle_2d(K, detJ, J);
376
377
378
// Compute constants
379
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
380
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
381
382
// Get coordinates and map to the reference (FIAT) element
383
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
384
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
385
386
// Compute number of derivatives.
387
unsigned int num_derivatives = 1;
388
for (unsigned int r = 0; r < n; r++)
389
{
390
num_derivatives *= 2;
391
}// end loop over 'r'
392
393
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
394
unsigned int **combinations = new unsigned int *[num_derivatives];
395
for (unsigned int row = 0; row < num_derivatives; row++)
396
{
397
combinations[row] = new unsigned int [n];
398
for (unsigned int col = 0; col < n; col++)
399
combinations[row][col] = 0;
400
}
401
402
// Generate combinations of derivatives
403
for (unsigned int row = 1; row < num_derivatives; row++)
404
{
405
for (unsigned int num = 0; num < row; num++)
406
{
407
for (unsigned int col = n-1; col+1 > 0; col--)
408
{
409
if (combinations[row][col] + 1 > 1)
410
combinations[row][col] = 0;
411
else
412
{
413
combinations[row][col] += 1;
414
break;
415
}
416
}
417
}
418
}
419
420
// Compute inverse of Jacobian
421
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
422
423
// Declare transformation matrix
424
// Declare pointer to two dimensional array and initialise
425
double **transform = new double *[num_derivatives];
426
427
for (unsigned int j = 0; j < num_derivatives; j++)
428
{
429
transform[j] = new double [num_derivatives];
430
for (unsigned int k = 0; k < num_derivatives; k++)
431
transform[j][k] = 1;
432
}
433
434
// Construct transformation matrix
435
for (unsigned int row = 0; row < num_derivatives; row++)
436
{
437
for (unsigned int col = 0; col < num_derivatives; col++)
438
{
439
for (unsigned int k = 0; k < n; k++)
440
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
441
}
442
}
443
444
// Reset values. Assuming that values is always an array.
445
for (unsigned int r = 0; r < num_derivatives; r++)
446
{
447
values[r] = 0.0;
448
}// end loop over 'r'
449
450
switch (i)
451
{
452
case 0:
453
{
454
455
// Array of basisvalues
456
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
457
458
// Declare helper variables
459
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
460
double tmp1 = (1.0 - Y)/2.0;
461
double tmp2 = tmp1*tmp1;
462
463
// Compute basisvalues
464
basisvalues[0] = 1.0;
465
basisvalues[1] = tmp0;
466
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
467
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
468
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
469
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
470
basisvalues[0] *= std::sqrt(0.5);
471
basisvalues[2] *= std::sqrt(1.0);
472
basisvalues[5] *= std::sqrt(1.5);
473
basisvalues[1] *= std::sqrt(3.0);
474
basisvalues[4] *= std::sqrt(4.5);
475
basisvalues[3] *= std::sqrt(7.5);
476
477
// Table(s) of coefficients
478
static const double coefficients0[6] = \
479
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
480
481
// Tables of derivatives of the polynomial base (transpose).
482
static const double dmats0[6][6] = \
483
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
484
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
485
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
486
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
487
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
488
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
489
490
static const double dmats1[6][6] = \
491
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
492
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
493
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
494
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
495
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
496
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
497
498
// Compute reference derivatives.
499
// Declare pointer to array of derivatives on FIAT element.
500
double *derivatives = new double[num_derivatives];
501
for (unsigned int r = 0; r < num_derivatives; r++)
502
{
503
derivatives[r] = 0.0;
504
}// end loop over 'r'
505
506
// Declare derivative matrix (of polynomial basis).
507
double dmats[6][6] = \
508
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
509
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
510
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
511
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
512
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
513
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
514
515
// Declare (auxiliary) derivative matrix (of polynomial basis).
516
double dmats_old[6][6] = \
517
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
518
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
519
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
520
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
521
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
522
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
523
524
// Loop possible derivatives.
525
for (unsigned int r = 0; r < num_derivatives; r++)
526
{
527
// Resetting dmats values to compute next derivative.
528
for (unsigned int t = 0; t < 6; t++)
529
{
530
for (unsigned int u = 0; u < 6; u++)
531
{
532
dmats[t][u] = 0.0;
533
if (t == u)
534
{
535
dmats[t][u] = 1.0;
536
}
537
538
}// end loop over 'u'
539
}// end loop over 't'
540
541
// Looping derivative order to generate dmats.
542
for (unsigned int s = 0; s < n; s++)
543
{
544
// Updating dmats_old with new values and resetting dmats.
545
for (unsigned int t = 0; t < 6; t++)
546
{
547
for (unsigned int u = 0; u < 6; u++)
548
{
549
dmats_old[t][u] = dmats[t][u];
550
dmats[t][u] = 0.0;
551
}// end loop over 'u'
552
}// end loop over 't'
553
554
// Update dmats using an inner product.
555
if (combinations[r][s] == 0)
556
{
557
for (unsigned int t = 0; t < 6; t++)
558
{
559
for (unsigned int u = 0; u < 6; u++)
560
{
561
for (unsigned int tu = 0; tu < 6; tu++)
562
{
563
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
564
}// end loop over 'tu'
565
}// end loop over 'u'
566
}// end loop over 't'
567
}
568
569
if (combinations[r][s] == 1)
570
{
571
for (unsigned int t = 0; t < 6; t++)
572
{
573
for (unsigned int u = 0; u < 6; u++)
574
{
575
for (unsigned int tu = 0; tu < 6; tu++)
576
{
577
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
578
}// end loop over 'tu'
579
}// end loop over 'u'
580
}// end loop over 't'
581
}
582
583
}// end loop over 's'
584
for (unsigned int s = 0; s < 6; s++)
585
{
586
for (unsigned int t = 0; t < 6; t++)
587
{
588
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
589
}// end loop over 't'
590
}// end loop over 's'
591
}// end loop over 'r'
592
593
// Transform derivatives back to physical element
594
for (unsigned int r = 0; r < num_derivatives; r++)
595
{
596
for (unsigned int s = 0; s < num_derivatives; s++)
597
{
598
values[r] += transform[r][s]*derivatives[s];
599
}// end loop over 's'
600
}// end loop over 'r'
601
602
// Delete pointer to array of derivatives on FIAT element
603
delete [] derivatives;
604
605
// Delete pointer to array of combinations of derivatives and transform
606
for (unsigned int r = 0; r < num_derivatives; r++)
607
{
608
delete [] combinations[r];
609
}// end loop over 'r'
610
delete [] combinations;
611
for (unsigned int r = 0; r < num_derivatives; r++)
612
{
613
delete [] transform[r];
614
}// end loop over 'r'
615
delete [] transform;
616
break;
617
}
618
case 1:
619
{
620
621
// Array of basisvalues
622
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
623
624
// Declare helper variables
625
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
626
double tmp1 = (1.0 - Y)/2.0;
627
double tmp2 = tmp1*tmp1;
628
629
// Compute basisvalues
630
basisvalues[0] = 1.0;
631
basisvalues[1] = tmp0;
632
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
633
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
634
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
635
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
636
basisvalues[0] *= std::sqrt(0.5);
637
basisvalues[2] *= std::sqrt(1.0);
638
basisvalues[5] *= std::sqrt(1.5);
639
basisvalues[1] *= std::sqrt(3.0);
640
basisvalues[4] *= std::sqrt(4.5);
641
basisvalues[3] *= std::sqrt(7.5);
642
643
// Table(s) of coefficients
644
static const double coefficients0[6] = \
645
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
646
647
// Tables of derivatives of the polynomial base (transpose).
648
static const double dmats0[6][6] = \
649
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
650
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
651
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
652
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
653
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
654
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
655
656
static const double dmats1[6][6] = \
657
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
658
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
659
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
660
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
661
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
662
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
663
664
// Compute reference derivatives.
665
// Declare pointer to array of derivatives on FIAT element.
666
double *derivatives = new double[num_derivatives];
667
for (unsigned int r = 0; r < num_derivatives; r++)
668
{
669
derivatives[r] = 0.0;
670
}// end loop over 'r'
671
672
// Declare derivative matrix (of polynomial basis).
673
double dmats[6][6] = \
674
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
675
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
676
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
677
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
678
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
679
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
680
681
// Declare (auxiliary) derivative matrix (of polynomial basis).
682
double dmats_old[6][6] = \
683
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
684
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
685
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
686
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
687
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
688
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
689
690
// Loop possible derivatives.
691
for (unsigned int r = 0; r < num_derivatives; r++)
692
{
693
// Resetting dmats values to compute next derivative.
694
for (unsigned int t = 0; t < 6; t++)
695
{
696
for (unsigned int u = 0; u < 6; u++)
697
{
698
dmats[t][u] = 0.0;
699
if (t == u)
700
{
701
dmats[t][u] = 1.0;
702
}
703
704
}// end loop over 'u'
705
}// end loop over 't'
706
707
// Looping derivative order to generate dmats.
708
for (unsigned int s = 0; s < n; s++)
709
{
710
// Updating dmats_old with new values and resetting dmats.
711
for (unsigned int t = 0; t < 6; t++)
712
{
713
for (unsigned int u = 0; u < 6; u++)
714
{
715
dmats_old[t][u] = dmats[t][u];
716
dmats[t][u] = 0.0;
717
}// end loop over 'u'
718
}// end loop over 't'
719
720
// Update dmats using an inner product.
721
if (combinations[r][s] == 0)
722
{
723
for (unsigned int t = 0; t < 6; t++)
724
{
725
for (unsigned int u = 0; u < 6; u++)
726
{
727
for (unsigned int tu = 0; tu < 6; tu++)
728
{
729
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
730
}// end loop over 'tu'
731
}// end loop over 'u'
732
}// end loop over 't'
733
}
734
735
if (combinations[r][s] == 1)
736
{
737
for (unsigned int t = 0; t < 6; t++)
738
{
739
for (unsigned int u = 0; u < 6; u++)
740
{
741
for (unsigned int tu = 0; tu < 6; tu++)
742
{
743
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
744
}// end loop over 'tu'
745
}// end loop over 'u'
746
}// end loop over 't'
747
}
748
749
}// end loop over 's'
750
for (unsigned int s = 0; s < 6; s++)
751
{
752
for (unsigned int t = 0; t < 6; t++)
753
{
754
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
755
}// end loop over 't'
756
}// end loop over 's'
757
}// end loop over 'r'
758
759
// Transform derivatives back to physical element
760
for (unsigned int r = 0; r < num_derivatives; r++)
761
{
762
for (unsigned int s = 0; s < num_derivatives; s++)
763
{
764
values[r] += transform[r][s]*derivatives[s];
765
}// end loop over 's'
766
}// end loop over 'r'
767
768
// Delete pointer to array of derivatives on FIAT element
769
delete [] derivatives;
770
771
// Delete pointer to array of combinations of derivatives and transform
772
for (unsigned int r = 0; r < num_derivatives; r++)
773
{
774
delete [] combinations[r];
775
}// end loop over 'r'
776
delete [] combinations;
777
for (unsigned int r = 0; r < num_derivatives; r++)
778
{
779
delete [] transform[r];
780
}// end loop over 'r'
781
delete [] transform;
782
break;
783
}
784
case 2:
785
{
786
787
// Array of basisvalues
788
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
789
790
// Declare helper variables
791
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
792
double tmp1 = (1.0 - Y)/2.0;
793
double tmp2 = tmp1*tmp1;
794
795
// Compute basisvalues
796
basisvalues[0] = 1.0;
797
basisvalues[1] = tmp0;
798
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
799
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
800
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
801
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
802
basisvalues[0] *= std::sqrt(0.5);
803
basisvalues[2] *= std::sqrt(1.0);
804
basisvalues[5] *= std::sqrt(1.5);
805
basisvalues[1] *= std::sqrt(3.0);
806
basisvalues[4] *= std::sqrt(4.5);
807
basisvalues[3] *= std::sqrt(7.5);
808
809
// Table(s) of coefficients
810
static const double coefficients0[6] = \
811
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
812
813
// Tables of derivatives of the polynomial base (transpose).
814
static const double dmats0[6][6] = \
815
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
816
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
817
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
818
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
819
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
820
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
821
822
static const double dmats1[6][6] = \
823
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
824
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
825
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
826
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
827
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
828
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
829
830
// Compute reference derivatives.
831
// Declare pointer to array of derivatives on FIAT element.
832
double *derivatives = new double[num_derivatives];
833
for (unsigned int r = 0; r < num_derivatives; r++)
834
{
835
derivatives[r] = 0.0;
836
}// end loop over 'r'
837
838
// Declare derivative matrix (of polynomial basis).
839
double dmats[6][6] = \
840
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
841
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
842
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
843
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
844
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
845
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
846
847
// Declare (auxiliary) derivative matrix (of polynomial basis).
848
double dmats_old[6][6] = \
849
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
850
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
851
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
852
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
853
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
854
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
855
856
// Loop possible derivatives.
857
for (unsigned int r = 0; r < num_derivatives; r++)
858
{
859
// Resetting dmats values to compute next derivative.
860
for (unsigned int t = 0; t < 6; t++)
861
{
862
for (unsigned int u = 0; u < 6; u++)
863
{
864
dmats[t][u] = 0.0;
865
if (t == u)
866
{
867
dmats[t][u] = 1.0;
868
}
869
870
}// end loop over 'u'
871
}// end loop over 't'
872
873
// Looping derivative order to generate dmats.
874
for (unsigned int s = 0; s < n; s++)
875
{
876
// Updating dmats_old with new values and resetting dmats.
877
for (unsigned int t = 0; t < 6; t++)
878
{
879
for (unsigned int u = 0; u < 6; u++)
880
{
881
dmats_old[t][u] = dmats[t][u];
882
dmats[t][u] = 0.0;
883
}// end loop over 'u'
884
}// end loop over 't'
885
886
// Update dmats using an inner product.
887
if (combinations[r][s] == 0)
888
{
889
for (unsigned int t = 0; t < 6; t++)
890
{
891
for (unsigned int u = 0; u < 6; u++)
892
{
893
for (unsigned int tu = 0; tu < 6; tu++)
894
{
895
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
896
}// end loop over 'tu'
897
}// end loop over 'u'
898
}// end loop over 't'
899
}
900
901
if (combinations[r][s] == 1)
902
{
903
for (unsigned int t = 0; t < 6; t++)
904
{
905
for (unsigned int u = 0; u < 6; u++)
906
{
907
for (unsigned int tu = 0; tu < 6; tu++)
908
{
909
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
910
}// end loop over 'tu'
911
}// end loop over 'u'
912
}// end loop over 't'
913
}
914
915
}// end loop over 's'
916
for (unsigned int s = 0; s < 6; s++)
917
{
918
for (unsigned int t = 0; t < 6; t++)
919
{
920
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
921
}// end loop over 't'
922
}// end loop over 's'
923
}// end loop over 'r'
924
925
// Transform derivatives back to physical element
926
for (unsigned int r = 0; r < num_derivatives; r++)
927
{
928
for (unsigned int s = 0; s < num_derivatives; s++)
929
{
930
values[r] += transform[r][s]*derivatives[s];
931
}// end loop over 's'
932
}// end loop over 'r'
933
934
// Delete pointer to array of derivatives on FIAT element
935
delete [] derivatives;
936
937
// Delete pointer to array of combinations of derivatives and transform
938
for (unsigned int r = 0; r < num_derivatives; r++)
939
{
940
delete [] combinations[r];
941
}// end loop over 'r'
942
delete [] combinations;
943
for (unsigned int r = 0; r < num_derivatives; r++)
944
{
945
delete [] transform[r];
946
}// end loop over 'r'
947
delete [] transform;
948
break;
949
}
950
case 3:
951
{
952
953
// Array of basisvalues
954
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
955
956
// Declare helper variables
957
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
958
double tmp1 = (1.0 - Y)/2.0;
959
double tmp2 = tmp1*tmp1;
960
961
// Compute basisvalues
962
basisvalues[0] = 1.0;
963
basisvalues[1] = tmp0;
964
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
965
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
966
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
967
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
968
basisvalues[0] *= std::sqrt(0.5);
969
basisvalues[2] *= std::sqrt(1.0);
970
basisvalues[5] *= std::sqrt(1.5);
971
basisvalues[1] *= std::sqrt(3.0);
972
basisvalues[4] *= std::sqrt(4.5);
973
basisvalues[3] *= std::sqrt(7.5);
974
975
// Table(s) of coefficients
976
static const double coefficients0[6] = \
977
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
978
979
// Tables of derivatives of the polynomial base (transpose).
980
static const double dmats0[6][6] = \
981
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
982
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
983
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
984
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
985
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
986
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
987
988
static const double dmats1[6][6] = \
989
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
990
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
991
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
992
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
993
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
994
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
995
996
// Compute reference derivatives.
997
// Declare pointer to array of derivatives on FIAT element.
998
double *derivatives = new double[num_derivatives];
999
for (unsigned int r = 0; r < num_derivatives; r++)
1000
{
1001
derivatives[r] = 0.0;
1002
}// end loop over 'r'
1003
1004
// Declare derivative matrix (of polynomial basis).
1005
double dmats[6][6] = \
1006
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1007
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1008
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1009
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1010
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1011
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1012
1013
// Declare (auxiliary) derivative matrix (of polynomial basis).
1014
double dmats_old[6][6] = \
1015
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1016
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1017
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1018
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1019
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1020
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1021
1022
// Loop possible derivatives.
1023
for (unsigned int r = 0; r < num_derivatives; r++)
1024
{
1025
// Resetting dmats values to compute next derivative.
1026
for (unsigned int t = 0; t < 6; t++)
1027
{
1028
for (unsigned int u = 0; u < 6; u++)
1029
{
1030
dmats[t][u] = 0.0;
1031
if (t == u)
1032
{
1033
dmats[t][u] = 1.0;
1034
}
1035
1036
}// end loop over 'u'
1037
}// end loop over 't'
1038
1039
// Looping derivative order to generate dmats.
1040
for (unsigned int s = 0; s < n; s++)
1041
{
1042
// Updating dmats_old with new values and resetting dmats.
1043
for (unsigned int t = 0; t < 6; t++)
1044
{
1045
for (unsigned int u = 0; u < 6; u++)
1046
{
1047
dmats_old[t][u] = dmats[t][u];
1048
dmats[t][u] = 0.0;
1049
}// end loop over 'u'
1050
}// end loop over 't'
1051
1052
// Update dmats using an inner product.
1053
if (combinations[r][s] == 0)
1054
{
1055
for (unsigned int t = 0; t < 6; t++)
1056
{
1057
for (unsigned int u = 0; u < 6; u++)
1058
{
1059
for (unsigned int tu = 0; tu < 6; tu++)
1060
{
1061
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1062
}// end loop over 'tu'
1063
}// end loop over 'u'
1064
}// end loop over 't'
1065
}
1066
1067
if (combinations[r][s] == 1)
1068
{
1069
for (unsigned int t = 0; t < 6; t++)
1070
{
1071
for (unsigned int u = 0; u < 6; u++)
1072
{
1073
for (unsigned int tu = 0; tu < 6; tu++)
1074
{
1075
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1076
}// end loop over 'tu'
1077
}// end loop over 'u'
1078
}// end loop over 't'
1079
}
1080
1081
}// end loop over 's'
1082
for (unsigned int s = 0; s < 6; s++)
1083
{
1084
for (unsigned int t = 0; t < 6; t++)
1085
{
1086
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1087
}// end loop over 't'
1088
}// end loop over 's'
1089
}// end loop over 'r'
1090
1091
// Transform derivatives back to physical element
1092
for (unsigned int r = 0; r < num_derivatives; r++)
1093
{
1094
for (unsigned int s = 0; s < num_derivatives; s++)
1095
{
1096
values[r] += transform[r][s]*derivatives[s];
1097
}// end loop over 's'
1098
}// end loop over 'r'
1099
1100
// Delete pointer to array of derivatives on FIAT element
1101
delete [] derivatives;
1102
1103
// Delete pointer to array of combinations of derivatives and transform
1104
for (unsigned int r = 0; r < num_derivatives; r++)
1105
{
1106
delete [] combinations[r];
1107
}// end loop over 'r'
1108
delete [] combinations;
1109
for (unsigned int r = 0; r < num_derivatives; r++)
1110
{
1111
delete [] transform[r];
1112
}// end loop over 'r'
1113
delete [] transform;
1114
break;
1115
}
1116
case 4:
1117
{
1118
1119
// Array of basisvalues
1120
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1121
1122
// Declare helper variables
1123
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1124
double tmp1 = (1.0 - Y)/2.0;
1125
double tmp2 = tmp1*tmp1;
1126
1127
// Compute basisvalues
1128
basisvalues[0] = 1.0;
1129
basisvalues[1] = tmp0;
1130
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1131
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1132
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1133
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
1134
basisvalues[0] *= std::sqrt(0.5);
1135
basisvalues[2] *= std::sqrt(1.0);
1136
basisvalues[5] *= std::sqrt(1.5);
1137
basisvalues[1] *= std::sqrt(3.0);
1138
basisvalues[4] *= std::sqrt(4.5);
1139
basisvalues[3] *= std::sqrt(7.5);
1140
1141
// Table(s) of coefficients
1142
static const double coefficients0[6] = \
1143
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
1144
1145
// Tables of derivatives of the polynomial base (transpose).
1146
static const double dmats0[6][6] = \
1147
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1148
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
1149
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1150
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
1151
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
1152
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1153
1154
static const double dmats1[6][6] = \
1155
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1156
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
1157
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
1158
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
1159
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
1160
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
1161
1162
// Compute reference derivatives.
1163
// Declare pointer to array of derivatives on FIAT element.
1164
double *derivatives = new double[num_derivatives];
1165
for (unsigned int r = 0; r < num_derivatives; r++)
1166
{
1167
derivatives[r] = 0.0;
1168
}// end loop over 'r'
1169
1170
// Declare derivative matrix (of polynomial basis).
1171
double dmats[6][6] = \
1172
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1173
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1174
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1175
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1176
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1177
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1178
1179
// Declare (auxiliary) derivative matrix (of polynomial basis).
1180
double dmats_old[6][6] = \
1181
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1182
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1183
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1184
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1185
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1186
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1187
1188
// Loop possible derivatives.
1189
for (unsigned int r = 0; r < num_derivatives; r++)
1190
{
1191
// Resetting dmats values to compute next derivative.
1192
for (unsigned int t = 0; t < 6; t++)
1193
{
1194
for (unsigned int u = 0; u < 6; u++)
1195
{
1196
dmats[t][u] = 0.0;
1197
if (t == u)
1198
{
1199
dmats[t][u] = 1.0;
1200
}
1201
1202
}// end loop over 'u'
1203
}// end loop over 't'
1204
1205
// Looping derivative order to generate dmats.
1206
for (unsigned int s = 0; s < n; s++)
1207
{
1208
// Updating dmats_old with new values and resetting dmats.
1209
for (unsigned int t = 0; t < 6; t++)
1210
{
1211
for (unsigned int u = 0; u < 6; u++)
1212
{
1213
dmats_old[t][u] = dmats[t][u];
1214
dmats[t][u] = 0.0;
1215
}// end loop over 'u'
1216
}// end loop over 't'
1217
1218
// Update dmats using an inner product.
1219
if (combinations[r][s] == 0)
1220
{
1221
for (unsigned int t = 0; t < 6; t++)
1222
{
1223
for (unsigned int u = 0; u < 6; u++)
1224
{
1225
for (unsigned int tu = 0; tu < 6; tu++)
1226
{
1227
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1228
}// end loop over 'tu'
1229
}// end loop over 'u'
1230
}// end loop over 't'
1231
}
1232
1233
if (combinations[r][s] == 1)
1234
{
1235
for (unsigned int t = 0; t < 6; t++)
1236
{
1237
for (unsigned int u = 0; u < 6; u++)
1238
{
1239
for (unsigned int tu = 0; tu < 6; tu++)
1240
{
1241
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1242
}// end loop over 'tu'
1243
}// end loop over 'u'
1244
}// end loop over 't'
1245
}
1246
1247
}// end loop over 's'
1248
for (unsigned int s = 0; s < 6; s++)
1249
{
1250
for (unsigned int t = 0; t < 6; t++)
1251
{
1252
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1253
}// end loop over 't'
1254
}// end loop over 's'
1255
}// end loop over 'r'
1256
1257
// Transform derivatives back to physical element
1258
for (unsigned int r = 0; r < num_derivatives; r++)
1259
{
1260
for (unsigned int s = 0; s < num_derivatives; s++)
1261
{
1262
values[r] += transform[r][s]*derivatives[s];
1263
}// end loop over 's'
1264
}// end loop over 'r'
1265
1266
// Delete pointer to array of derivatives on FIAT element
1267
delete [] derivatives;
1268
1269
// Delete pointer to array of combinations of derivatives and transform
1270
for (unsigned int r = 0; r < num_derivatives; r++)
1271
{
1272
delete [] combinations[r];
1273
}// end loop over 'r'
1274
delete [] combinations;
1275
for (unsigned int r = 0; r < num_derivatives; r++)
1276
{
1277
delete [] transform[r];
1278
}// end loop over 'r'
1279
delete [] transform;
1280
break;
1281
}
1282
case 5:
1283
{
1284
1285
// Array of basisvalues
1286
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1287
1288
// Declare helper variables
1289
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1290
double tmp1 = (1.0 - Y)/2.0;
1291
double tmp2 = tmp1*tmp1;
1292
1293
// Compute basisvalues
1294
basisvalues[0] = 1.0;
1295
basisvalues[1] = tmp0;
1296
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1297
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1298
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1299
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
1300
basisvalues[0] *= std::sqrt(0.5);
1301
basisvalues[2] *= std::sqrt(1.0);
1302
basisvalues[5] *= std::sqrt(1.5);
1303
basisvalues[1] *= std::sqrt(3.0);
1304
basisvalues[4] *= std::sqrt(4.5);
1305
basisvalues[3] *= std::sqrt(7.5);
1306
1307
// Table(s) of coefficients
1308
static const double coefficients0[6] = \
1309
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
1310
1311
// Tables of derivatives of the polynomial base (transpose).
1312
static const double dmats0[6][6] = \
1313
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1314
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
1315
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1316
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
1317
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
1318
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1319
1320
static const double dmats1[6][6] = \
1321
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1322
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
1323
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
1324
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
1325
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
1326
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
1327
1328
// Compute reference derivatives.
1329
// Declare pointer to array of derivatives on FIAT element.
1330
double *derivatives = new double[num_derivatives];
1331
for (unsigned int r = 0; r < num_derivatives; r++)
1332
{
1333
derivatives[r] = 0.0;
1334
}// end loop over 'r'
1335
1336
// Declare derivative matrix (of polynomial basis).
1337
double dmats[6][6] = \
1338
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1339
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1340
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1341
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1342
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1343
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1344
1345
// Declare (auxiliary) derivative matrix (of polynomial basis).
1346
double dmats_old[6][6] = \
1347
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1348
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1349
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1350
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1351
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1352
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1353
1354
// Loop possible derivatives.
1355
for (unsigned int r = 0; r < num_derivatives; r++)
1356
{
1357
// Resetting dmats values to compute next derivative.
1358
for (unsigned int t = 0; t < 6; t++)
1359
{
1360
for (unsigned int u = 0; u < 6; u++)
1361
{
1362
dmats[t][u] = 0.0;
1363
if (t == u)
1364
{
1365
dmats[t][u] = 1.0;
1366
}
1367
1368
}// end loop over 'u'
1369
}// end loop over 't'
1370
1371
// Looping derivative order to generate dmats.
1372
for (unsigned int s = 0; s < n; s++)
1373
{
1374
// Updating dmats_old with new values and resetting dmats.
1375
for (unsigned int t = 0; t < 6; t++)
1376
{
1377
for (unsigned int u = 0; u < 6; u++)
1378
{
1379
dmats_old[t][u] = dmats[t][u];
1380
dmats[t][u] = 0.0;
1381
}// end loop over 'u'
1382
}// end loop over 't'
1383
1384
// Update dmats using an inner product.
1385
if (combinations[r][s] == 0)
1386
{
1387
for (unsigned int t = 0; t < 6; t++)
1388
{
1389
for (unsigned int u = 0; u < 6; u++)
1390
{
1391
for (unsigned int tu = 0; tu < 6; tu++)
1392
{
1393
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1394
}// end loop over 'tu'
1395
}// end loop over 'u'
1396
}// end loop over 't'
1397
}
1398
1399
if (combinations[r][s] == 1)
1400
{
1401
for (unsigned int t = 0; t < 6; t++)
1402
{
1403
for (unsigned int u = 0; u < 6; u++)
1404
{
1405
for (unsigned int tu = 0; tu < 6; tu++)
1406
{
1407
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1408
}// end loop over 'tu'
1409
}// end loop over 'u'
1410
}// end loop over 't'
1411
}
1412
1413
}// end loop over 's'
1414
for (unsigned int s = 0; s < 6; s++)
1415
{
1416
for (unsigned int t = 0; t < 6; t++)
1417
{
1418
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1419
}// end loop over 't'
1420
}// end loop over 's'
1421
}// end loop over 'r'
1422
1423
// Transform derivatives back to physical element
1424
for (unsigned int r = 0; r < num_derivatives; r++)
1425
{
1426
for (unsigned int s = 0; s < num_derivatives; s++)
1427
{
1428
values[r] += transform[r][s]*derivatives[s];
1429
}// end loop over 's'
1430
}// end loop over 'r'
1431
1432
// Delete pointer to array of derivatives on FIAT element
1433
delete [] derivatives;
1434
1435
// Delete pointer to array of combinations of derivatives and transform
1436
for (unsigned int r = 0; r < num_derivatives; r++)
1437
{
1438
delete [] combinations[r];
1439
}// end loop over 'r'
1440
delete [] combinations;
1441
for (unsigned int r = 0; r < num_derivatives; r++)
1442
{
1443
delete [] transform[r];
1444
}// end loop over 'r'
1445
delete [] transform;
1446
break;
1447
}
1448
}
1449
1450
}
1451
1452
/// Evaluate order n derivatives of all basis functions at given point x in cell
1453
virtual void evaluate_basis_derivatives_all(std::size_t n,
1454
double* values,
1455
const double* x,
1456
const double* vertex_coordinates,
1457
int cell_orientation) const
1458
{
1459
// Compute number of derivatives.
1460
unsigned int num_derivatives = 1;
1461
for (unsigned int r = 0; r < n; r++)
1462
{
1463
num_derivatives *= 2;
1464
}// end loop over 'r'
1465
1466
// Helper variable to hold values of a single dof.
1467
double *dof_values = new double[num_derivatives];
1468
for (unsigned int r = 0; r < num_derivatives; r++)
1469
{
1470
dof_values[r] = 0.0;
1471
}// end loop over 'r'
1472
1473
// Loop dofs and call evaluate_basis_derivatives.
1474
for (unsigned int r = 0; r < 6; r++)
1475
{
1476
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1477
for (unsigned int s = 0; s < num_derivatives; s++)
1478
{
1479
values[r*num_derivatives + s] = dof_values[s];
1480
}// end loop over 's'
1481
}// end loop over 'r'
1482
1483
// Delete pointer.
1484
delete [] dof_values;
1485
}
1486
1487
/// Evaluate linear functional for dof i on the function f
1488
virtual double evaluate_dof(std::size_t i,
1489
const ufc::function& f,
1490
const double* vertex_coordinates,
1491
int cell_orientation,
1492
const ufc::cell& c) const
1493
{
1494
// Declare variables for result of evaluation
1495
double vals[1];
1496
1497
// Declare variable for physical coordinates
1498
double y[2];
1499
switch (i)
1500
{
1501
case 0:
1502
{
1503
y[0] = vertex_coordinates[0];
1504
y[1] = vertex_coordinates[1];
1505
f.evaluate(vals, y, c);
1506
return vals[0];
1507
break;
1508
}
1509
case 1:
1510
{
1511
y[0] = vertex_coordinates[2];
1512
y[1] = vertex_coordinates[3];
1513
f.evaluate(vals, y, c);
1514
return vals[0];
1515
break;
1516
}
1517
case 2:
1518
{
1519
y[0] = vertex_coordinates[4];
1520
y[1] = vertex_coordinates[5];
1521
f.evaluate(vals, y, c);
1522
return vals[0];
1523
break;
1524
}
1525
case 3:
1526
{
1527
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
1528
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1529
f.evaluate(vals, y, c);
1530
return vals[0];
1531
break;
1532
}
1533
case 4:
1534
{
1535
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
1536
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1537
f.evaluate(vals, y, c);
1538
return vals[0];
1539
break;
1540
}
1541
case 5:
1542
{
1543
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
1544
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1545
f.evaluate(vals, y, c);
1546
return vals[0];
1547
break;
1548
}
1549
}
1550
1551
return 0.0;
1552
}
1553
1554
/// Evaluate linear functionals for all dofs on the function f
1555
virtual void evaluate_dofs(double* values,
1556
const ufc::function& f,
1557
const double* vertex_coordinates,
1558
int cell_orientation,
1559
const ufc::cell& c) const
1560
{
1561
// Declare variables for result of evaluation
1562
double vals[1];
1563
1564
// Declare variable for physical coordinates
1565
double y[2];
1566
y[0] = vertex_coordinates[0];
1567
y[1] = vertex_coordinates[1];
1568
f.evaluate(vals, y, c);
1569
values[0] = vals[0];
1570
y[0] = vertex_coordinates[2];
1571
y[1] = vertex_coordinates[3];
1572
f.evaluate(vals, y, c);
1573
values[1] = vals[0];
1574
y[0] = vertex_coordinates[4];
1575
y[1] = vertex_coordinates[5];
1576
f.evaluate(vals, y, c);
1577
values[2] = vals[0];
1578
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
1579
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1580
f.evaluate(vals, y, c);
1581
values[3] = vals[0];
1582
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
1583
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1584
f.evaluate(vals, y, c);
1585
values[4] = vals[0];
1586
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
1587
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1588
f.evaluate(vals, y, c);
1589
values[5] = vals[0];
1590
}
1591
1592
/// Interpolate vertex values from dof values
1593
virtual void interpolate_vertex_values(double* vertex_values,
1594
const double* dof_values,
1595
const double* vertex_coordinates,
1596
int cell_orientation,
1597
const ufc::cell& c) const
1598
{
1599
// Evaluate function and change variables
1600
vertex_values[0] = dof_values[0];
1601
vertex_values[1] = dof_values[1];
1602
vertex_values[2] = dof_values[2];
1603
}
1604
1605
/// Map coordinate xhat from reference cell to coordinate x in cell
1606
virtual void map_from_reference_cell(double* x,
1607
const double* xhat,
1608
const ufc::cell& c) const
1609
{
1610
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1611
}
1612
1613
/// Map from coordinate x in cell to coordinate xhat in reference cell
1614
virtual void map_to_reference_cell(double* xhat,
1615
const double* x,
1616
const ufc::cell& c) const
1617
{
1618
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1619
}
1620
1621
/// Return the number of sub elements (for a mixed element)
1622
virtual std::size_t num_sub_elements() const
1623
{
1624
return 0;
1625
}
1626
1627
/// Create a new finite element for sub element i (for a mixed element)
1628
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1629
{
1630
return 0;
1631
}
1632
1633
/// Create a new class instance
1634
virtual ufc::finite_element* create() const
1635
{
1636
return new pointmeasure_finite_element_0();
1637
}
1638
1639
};
1640
1641
/// This class defines the interface for a finite element.
1642
1643
class pointmeasure_finite_element_1: public ufc::finite_element
1644
{
1645
public:
1646
1647
/// Constructor
1648
pointmeasure_finite_element_1() : ufc::finite_element()
1649
{
1650
// Do nothing
1651
}
1652
1653
/// Destructor
1654
virtual ~pointmeasure_finite_element_1()
1655
{
1656
// Do nothing
1657
}
1658
1659
/// Return a string identifying the finite element
1660
virtual const char* signature() const
1661
{
1662
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
1663
}
1664
1665
/// Return the cell shape
1666
virtual ufc::shape cell_shape() const
1667
{
1668
return ufc::triangle;
1669
}
1670
1671
/// Return the topological dimension of the cell shape
1672
virtual std::size_t topological_dimension() const
1673
{
1674
return 2;
1675
}
1676
1677
/// Return the geometric dimension of the cell shape
1678
virtual std::size_t geometric_dimension() const
1679
{
1680
return 2;
1681
}
1682
1683
/// Return the dimension of the finite element function space
1684
virtual std::size_t space_dimension() const
1685
{
1686
return 3;
1687
}
1688
1689
/// Return the rank of the value space
1690
virtual std::size_t value_rank() const
1691
{
1692
return 0;
1693
}
1694
1695
/// Return the dimension of the value space for axis i
1696
virtual std::size_t value_dimension(std::size_t i) const
1697
{
1698
return 1;
1699
}
1700
1701
/// Evaluate basis function i at given point x in cell
1702
virtual void evaluate_basis(std::size_t i,
1703
double* values,
1704
const double* x,
1705
const double* vertex_coordinates,
1706
int cell_orientation) const
1707
{
1708
// Compute Jacobian
1709
double J[4];
1710
compute_jacobian_triangle_2d(J, vertex_coordinates);
1711
1712
// Compute Jacobian inverse and determinant
1713
double K[4];
1714
double detJ;
1715
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1716
1717
1718
// Compute constants
1719
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1720
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1721
1722
// Get coordinates and map to the reference (FIAT) element
1723
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1724
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1725
1726
// Reset values
1727
*values = 0.0;
1728
switch (i)
1729
{
1730
case 0:
1731
{
1732
1733
// Array of basisvalues
1734
double basisvalues[3] = {0.0, 0.0, 0.0};
1735
1736
// Declare helper variables
1737
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1738
1739
// Compute basisvalues
1740
basisvalues[0] = 1.0;
1741
basisvalues[1] = tmp0;
1742
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1743
basisvalues[0] *= std::sqrt(0.5);
1744
basisvalues[2] *= std::sqrt(1.0);
1745
basisvalues[1] *= std::sqrt(3.0);
1746
1747
// Table(s) of coefficients
1748
static const double coefficients0[3] = \
1749
{0.47140452, -0.28867513, -0.16666667};
1750
1751
// Compute value(s)
1752
for (unsigned int r = 0; r < 3; r++)
1753
{
1754
*values += coefficients0[r]*basisvalues[r];
1755
}// end loop over 'r'
1756
break;
1757
}
1758
case 1:
1759
{
1760
1761
// Array of basisvalues
1762
double basisvalues[3] = {0.0, 0.0, 0.0};
1763
1764
// Declare helper variables
1765
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1766
1767
// Compute basisvalues
1768
basisvalues[0] = 1.0;
1769
basisvalues[1] = tmp0;
1770
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1771
basisvalues[0] *= std::sqrt(0.5);
1772
basisvalues[2] *= std::sqrt(1.0);
1773
basisvalues[1] *= std::sqrt(3.0);
1774
1775
// Table(s) of coefficients
1776
static const double coefficients0[3] = \
1777
{0.47140452, 0.28867513, -0.16666667};
1778
1779
// Compute value(s)
1780
for (unsigned int r = 0; r < 3; r++)
1781
{
1782
*values += coefficients0[r]*basisvalues[r];
1783
}// end loop over 'r'
1784
break;
1785
}
1786
case 2:
1787
{
1788
1789
// Array of basisvalues
1790
double basisvalues[3] = {0.0, 0.0, 0.0};
1791
1792
// Declare helper variables
1793
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1794
1795
// Compute basisvalues
1796
basisvalues[0] = 1.0;
1797
basisvalues[1] = tmp0;
1798
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1799
basisvalues[0] *= std::sqrt(0.5);
1800
basisvalues[2] *= std::sqrt(1.0);
1801
basisvalues[1] *= std::sqrt(3.0);
1802
1803
// Table(s) of coefficients
1804
static const double coefficients0[3] = \
1805
{0.47140452, 0.0, 0.33333333};
1806
1807
// Compute value(s)
1808
for (unsigned int r = 0; r < 3; r++)
1809
{
1810
*values += coefficients0[r]*basisvalues[r];
1811
}// end loop over 'r'
1812
break;
1813
}
1814
}
1815
1816
}
1817
1818
/// Evaluate all basis functions at given point x in cell
1819
virtual void evaluate_basis_all(double* values,
1820
const double* x,
1821
const double* vertex_coordinates,
1822
int cell_orientation) const
1823
{
1824
// Helper variable to hold values of a single dof.
1825
double dof_values = 0.0;
1826
1827
// Loop dofs and call evaluate_basis
1828
for (unsigned int r = 0; r < 3; r++)
1829
{
1830
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
1831
values[r] = dof_values;
1832
}// end loop over 'r'
1833
}
1834
1835
/// Evaluate order n derivatives of basis function i at given point x in cell
1836
virtual void evaluate_basis_derivatives(std::size_t i,
1837
std::size_t n,
1838
double* values,
1839
const double* x,
1840
const double* vertex_coordinates,
1841
int cell_orientation) const
1842
{
1843
// Compute Jacobian
1844
double J[4];
1845
compute_jacobian_triangle_2d(J, vertex_coordinates);
1846
1847
// Compute Jacobian inverse and determinant
1848
double K[4];
1849
double detJ;
1850
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1851
1852
1853
// Compute constants
1854
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1855
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1856
1857
// Get coordinates and map to the reference (FIAT) element
1858
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1859
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1860
1861
// Compute number of derivatives.
1862
unsigned int num_derivatives = 1;
1863
for (unsigned int r = 0; r < n; r++)
1864
{
1865
num_derivatives *= 2;
1866
}// end loop over 'r'
1867
1868
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1869
unsigned int **combinations = new unsigned int *[num_derivatives];
1870
for (unsigned int row = 0; row < num_derivatives; row++)
1871
{
1872
combinations[row] = new unsigned int [n];
1873
for (unsigned int col = 0; col < n; col++)
1874
combinations[row][col] = 0;
1875
}
1876
1877
// Generate combinations of derivatives
1878
for (unsigned int row = 1; row < num_derivatives; row++)
1879
{
1880
for (unsigned int num = 0; num < row; num++)
1881
{
1882
for (unsigned int col = n-1; col+1 > 0; col--)
1883
{
1884
if (combinations[row][col] + 1 > 1)
1885
combinations[row][col] = 0;
1886
else
1887
{
1888
combinations[row][col] += 1;
1889
break;
1890
}
1891
}
1892
}
1893
}
1894
1895
// Compute inverse of Jacobian
1896
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
1897
1898
// Declare transformation matrix
1899
// Declare pointer to two dimensional array and initialise
1900
double **transform = new double *[num_derivatives];
1901
1902
for (unsigned int j = 0; j < num_derivatives; j++)
1903
{
1904
transform[j] = new double [num_derivatives];
1905
for (unsigned int k = 0; k < num_derivatives; k++)
1906
transform[j][k] = 1;
1907
}
1908
1909
// Construct transformation matrix
1910
for (unsigned int row = 0; row < num_derivatives; row++)
1911
{
1912
for (unsigned int col = 0; col < num_derivatives; col++)
1913
{
1914
for (unsigned int k = 0; k < n; k++)
1915
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1916
}
1917
}
1918
1919
// Reset values. Assuming that values is always an array.
1920
for (unsigned int r = 0; r < num_derivatives; r++)
1921
{
1922
values[r] = 0.0;
1923
}// end loop over 'r'
1924
1925
switch (i)
1926
{
1927
case 0:
1928
{
1929
1930
// Array of basisvalues
1931
double basisvalues[3] = {0.0, 0.0, 0.0};
1932
1933
// Declare helper variables
1934
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1935
1936
// Compute basisvalues
1937
basisvalues[0] = 1.0;
1938
basisvalues[1] = tmp0;
1939
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1940
basisvalues[0] *= std::sqrt(0.5);
1941
basisvalues[2] *= std::sqrt(1.0);
1942
basisvalues[1] *= std::sqrt(3.0);
1943
1944
// Table(s) of coefficients
1945
static const double coefficients0[3] = \
1946
{0.47140452, -0.28867513, -0.16666667};
1947
1948
// Tables of derivatives of the polynomial base (transpose).
1949
static const double dmats0[3][3] = \
1950
{{0.0, 0.0, 0.0},
1951
{4.8989795, 0.0, 0.0},
1952
{0.0, 0.0, 0.0}};
1953
1954
static const double dmats1[3][3] = \
1955
{{0.0, 0.0, 0.0},
1956
{2.4494897, 0.0, 0.0},
1957
{4.2426407, 0.0, 0.0}};
1958
1959
// Compute reference derivatives.
1960
// Declare pointer to array of derivatives on FIAT element.
1961
double *derivatives = new double[num_derivatives];
1962
for (unsigned int r = 0; r < num_derivatives; r++)
1963
{
1964
derivatives[r] = 0.0;
1965
}// end loop over 'r'
1966
1967
// Declare derivative matrix (of polynomial basis).
1968
double dmats[3][3] = \
1969
{{1.0, 0.0, 0.0},
1970
{0.0, 1.0, 0.0},
1971
{0.0, 0.0, 1.0}};
1972
1973
// Declare (auxiliary) derivative matrix (of polynomial basis).
1974
double dmats_old[3][3] = \
1975
{{1.0, 0.0, 0.0},
1976
{0.0, 1.0, 0.0},
1977
{0.0, 0.0, 1.0}};
1978
1979
// Loop possible derivatives.
1980
for (unsigned int r = 0; r < num_derivatives; r++)
1981
{
1982
// Resetting dmats values to compute next derivative.
1983
for (unsigned int t = 0; t < 3; t++)
1984
{
1985
for (unsigned int u = 0; u < 3; u++)
1986
{
1987
dmats[t][u] = 0.0;
1988
if (t == u)
1989
{
1990
dmats[t][u] = 1.0;
1991
}
1992
1993
}// end loop over 'u'
1994
}// end loop over 't'
1995
1996
// Looping derivative order to generate dmats.
1997
for (unsigned int s = 0; s < n; s++)
1998
{
1999
// Updating dmats_old with new values and resetting dmats.
2000
for (unsigned int t = 0; t < 3; t++)
2001
{
2002
for (unsigned int u = 0; u < 3; u++)
2003
{
2004
dmats_old[t][u] = dmats[t][u];
2005
dmats[t][u] = 0.0;
2006
}// end loop over 'u'
2007
}// end loop over 't'
2008
2009
// Update dmats using an inner product.
2010
if (combinations[r][s] == 0)
2011
{
2012
for (unsigned int t = 0; t < 3; t++)
2013
{
2014
for (unsigned int u = 0; u < 3; u++)
2015
{
2016
for (unsigned int tu = 0; tu < 3; tu++)
2017
{
2018
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2019
}// end loop over 'tu'
2020
}// end loop over 'u'
2021
}// end loop over 't'
2022
}
2023
2024
if (combinations[r][s] == 1)
2025
{
2026
for (unsigned int t = 0; t < 3; t++)
2027
{
2028
for (unsigned int u = 0; u < 3; u++)
2029
{
2030
for (unsigned int tu = 0; tu < 3; tu++)
2031
{
2032
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2033
}// end loop over 'tu'
2034
}// end loop over 'u'
2035
}// end loop over 't'
2036
}
2037
2038
}// end loop over 's'
2039
for (unsigned int s = 0; s < 3; s++)
2040
{
2041
for (unsigned int t = 0; t < 3; t++)
2042
{
2043
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2044
}// end loop over 't'
2045
}// end loop over 's'
2046
}// end loop over 'r'
2047
2048
// Transform derivatives back to physical element
2049
for (unsigned int r = 0; r < num_derivatives; r++)
2050
{
2051
for (unsigned int s = 0; s < num_derivatives; s++)
2052
{
2053
values[r] += transform[r][s]*derivatives[s];
2054
}// end loop over 's'
2055
}// end loop over 'r'
2056
2057
// Delete pointer to array of derivatives on FIAT element
2058
delete [] derivatives;
2059
2060
// Delete pointer to array of combinations of derivatives and transform
2061
for (unsigned int r = 0; r < num_derivatives; r++)
2062
{
2063
delete [] combinations[r];
2064
}// end loop over 'r'
2065
delete [] combinations;
2066
for (unsigned int r = 0; r < num_derivatives; r++)
2067
{
2068
delete [] transform[r];
2069
}// end loop over 'r'
2070
delete [] transform;
2071
break;
2072
}
2073
case 1:
2074
{
2075
2076
// Array of basisvalues
2077
double basisvalues[3] = {0.0, 0.0, 0.0};
2078
2079
// Declare helper variables
2080
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2081
2082
// Compute basisvalues
2083
basisvalues[0] = 1.0;
2084
basisvalues[1] = tmp0;
2085
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2086
basisvalues[0] *= std::sqrt(0.5);
2087
basisvalues[2] *= std::sqrt(1.0);
2088
basisvalues[1] *= std::sqrt(3.0);
2089
2090
// Table(s) of coefficients
2091
static const double coefficients0[3] = \
2092
{0.47140452, 0.28867513, -0.16666667};
2093
2094
// Tables of derivatives of the polynomial base (transpose).
2095
static const double dmats0[3][3] = \
2096
{{0.0, 0.0, 0.0},
2097
{4.8989795, 0.0, 0.0},
2098
{0.0, 0.0, 0.0}};
2099
2100
static const double dmats1[3][3] = \
2101
{{0.0, 0.0, 0.0},
2102
{2.4494897, 0.0, 0.0},
2103
{4.2426407, 0.0, 0.0}};
2104
2105
// Compute reference derivatives.
2106
// Declare pointer to array of derivatives on FIAT element.
2107
double *derivatives = new double[num_derivatives];
2108
for (unsigned int r = 0; r < num_derivatives; r++)
2109
{
2110
derivatives[r] = 0.0;
2111
}// end loop over 'r'
2112
2113
// Declare derivative matrix (of polynomial basis).
2114
double dmats[3][3] = \
2115
{{1.0, 0.0, 0.0},
2116
{0.0, 1.0, 0.0},
2117
{0.0, 0.0, 1.0}};
2118
2119
// Declare (auxiliary) derivative matrix (of polynomial basis).
2120
double dmats_old[3][3] = \
2121
{{1.0, 0.0, 0.0},
2122
{0.0, 1.0, 0.0},
2123
{0.0, 0.0, 1.0}};
2124
2125
// Loop possible derivatives.
2126
for (unsigned int r = 0; r < num_derivatives; r++)
2127
{
2128
// Resetting dmats values to compute next derivative.
2129
for (unsigned int t = 0; t < 3; t++)
2130
{
2131
for (unsigned int u = 0; u < 3; u++)
2132
{
2133
dmats[t][u] = 0.0;
2134
if (t == u)
2135
{
2136
dmats[t][u] = 1.0;
2137
}
2138
2139
}// end loop over 'u'
2140
}// end loop over 't'
2141
2142
// Looping derivative order to generate dmats.
2143
for (unsigned int s = 0; s < n; s++)
2144
{
2145
// Updating dmats_old with new values and resetting dmats.
2146
for (unsigned int t = 0; t < 3; t++)
2147
{
2148
for (unsigned int u = 0; u < 3; u++)
2149
{
2150
dmats_old[t][u] = dmats[t][u];
2151
dmats[t][u] = 0.0;
2152
}// end loop over 'u'
2153
}// end loop over 't'
2154
2155
// Update dmats using an inner product.
2156
if (combinations[r][s] == 0)
2157
{
2158
for (unsigned int t = 0; t < 3; t++)
2159
{
2160
for (unsigned int u = 0; u < 3; u++)
2161
{
2162
for (unsigned int tu = 0; tu < 3; tu++)
2163
{
2164
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2165
}// end loop over 'tu'
2166
}// end loop over 'u'
2167
}// end loop over 't'
2168
}
2169
2170
if (combinations[r][s] == 1)
2171
{
2172
for (unsigned int t = 0; t < 3; t++)
2173
{
2174
for (unsigned int u = 0; u < 3; u++)
2175
{
2176
for (unsigned int tu = 0; tu < 3; tu++)
2177
{
2178
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2179
}// end loop over 'tu'
2180
}// end loop over 'u'
2181
}// end loop over 't'
2182
}
2183
2184
}// end loop over 's'
2185
for (unsigned int s = 0; s < 3; s++)
2186
{
2187
for (unsigned int t = 0; t < 3; t++)
2188
{
2189
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2190
}// end loop over 't'
2191
}// end loop over 's'
2192
}// end loop over 'r'
2193
2194
// Transform derivatives back to physical element
2195
for (unsigned int r = 0; r < num_derivatives; r++)
2196
{
2197
for (unsigned int s = 0; s < num_derivatives; s++)
2198
{
2199
values[r] += transform[r][s]*derivatives[s];
2200
}// end loop over 's'
2201
}// end loop over 'r'
2202
2203
// Delete pointer to array of derivatives on FIAT element
2204
delete [] derivatives;
2205
2206
// Delete pointer to array of combinations of derivatives and transform
2207
for (unsigned int r = 0; r < num_derivatives; r++)
2208
{
2209
delete [] combinations[r];
2210
}// end loop over 'r'
2211
delete [] combinations;
2212
for (unsigned int r = 0; r < num_derivatives; r++)
2213
{
2214
delete [] transform[r];
2215
}// end loop over 'r'
2216
delete [] transform;
2217
break;
2218
}
2219
case 2:
2220
{
2221
2222
// Array of basisvalues
2223
double basisvalues[3] = {0.0, 0.0, 0.0};
2224
2225
// Declare helper variables
2226
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2227
2228
// Compute basisvalues
2229
basisvalues[0] = 1.0;
2230
basisvalues[1] = tmp0;
2231
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2232
basisvalues[0] *= std::sqrt(0.5);
2233
basisvalues[2] *= std::sqrt(1.0);
2234
basisvalues[1] *= std::sqrt(3.0);
2235
2236
// Table(s) of coefficients
2237
static const double coefficients0[3] = \
2238
{0.47140452, 0.0, 0.33333333};
2239
2240
// Tables of derivatives of the polynomial base (transpose).
2241
static const double dmats0[3][3] = \
2242
{{0.0, 0.0, 0.0},
2243
{4.8989795, 0.0, 0.0},
2244
{0.0, 0.0, 0.0}};
2245
2246
static const double dmats1[3][3] = \
2247
{{0.0, 0.0, 0.0},
2248
{2.4494897, 0.0, 0.0},
2249
{4.2426407, 0.0, 0.0}};
2250
2251
// Compute reference derivatives.
2252
// Declare pointer to array of derivatives on FIAT element.
2253
double *derivatives = new double[num_derivatives];
2254
for (unsigned int r = 0; r < num_derivatives; r++)
2255
{
2256
derivatives[r] = 0.0;
2257
}// end loop over 'r'
2258
2259
// Declare derivative matrix (of polynomial basis).
2260
double dmats[3][3] = \
2261
{{1.0, 0.0, 0.0},
2262
{0.0, 1.0, 0.0},
2263
{0.0, 0.0, 1.0}};
2264
2265
// Declare (auxiliary) derivative matrix (of polynomial basis).
2266
double dmats_old[3][3] = \
2267
{{1.0, 0.0, 0.0},
2268
{0.0, 1.0, 0.0},
2269
{0.0, 0.0, 1.0}};
2270
2271
// Loop possible derivatives.
2272
for (unsigned int r = 0; r < num_derivatives; r++)
2273
{
2274
// Resetting dmats values to compute next derivative.
2275
for (unsigned int t = 0; t < 3; t++)
2276
{
2277
for (unsigned int u = 0; u < 3; u++)
2278
{
2279
dmats[t][u] = 0.0;
2280
if (t == u)
2281
{
2282
dmats[t][u] = 1.0;
2283
}
2284
2285
}// end loop over 'u'
2286
}// end loop over 't'
2287
2288
// Looping derivative order to generate dmats.
2289
for (unsigned int s = 0; s < n; s++)
2290
{
2291
// Updating dmats_old with new values and resetting dmats.
2292
for (unsigned int t = 0; t < 3; t++)
2293
{
2294
for (unsigned int u = 0; u < 3; u++)
2295
{
2296
dmats_old[t][u] = dmats[t][u];
2297
dmats[t][u] = 0.0;
2298
}// end loop over 'u'
2299
}// end loop over 't'
2300
2301
// Update dmats using an inner product.
2302
if (combinations[r][s] == 0)
2303
{
2304
for (unsigned int t = 0; t < 3; t++)
2305
{
2306
for (unsigned int u = 0; u < 3; u++)
2307
{
2308
for (unsigned int tu = 0; tu < 3; tu++)
2309
{
2310
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2311
}// end loop over 'tu'
2312
}// end loop over 'u'
2313
}// end loop over 't'
2314
}
2315
2316
if (combinations[r][s] == 1)
2317
{
2318
for (unsigned int t = 0; t < 3; t++)
2319
{
2320
for (unsigned int u = 0; u < 3; u++)
2321
{
2322
for (unsigned int tu = 0; tu < 3; tu++)
2323
{
2324
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2325
}// end loop over 'tu'
2326
}// end loop over 'u'
2327
}// end loop over 't'
2328
}
2329
2330
}// end loop over 's'
2331
for (unsigned int s = 0; s < 3; s++)
2332
{
2333
for (unsigned int t = 0; t < 3; t++)
2334
{
2335
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2336
}// end loop over 't'
2337
}// end loop over 's'
2338
}// end loop over 'r'
2339
2340
// Transform derivatives back to physical element
2341
for (unsigned int r = 0; r < num_derivatives; r++)
2342
{
2343
for (unsigned int s = 0; s < num_derivatives; s++)
2344
{
2345
values[r] += transform[r][s]*derivatives[s];
2346
}// end loop over 's'
2347
}// end loop over 'r'
2348
2349
// Delete pointer to array of derivatives on FIAT element
2350
delete [] derivatives;
2351
2352
// Delete pointer to array of combinations of derivatives and transform
2353
for (unsigned int r = 0; r < num_derivatives; r++)
2354
{
2355
delete [] combinations[r];
2356
}// end loop over 'r'
2357
delete [] combinations;
2358
for (unsigned int r = 0; r < num_derivatives; r++)
2359
{
2360
delete [] transform[r];
2361
}// end loop over 'r'
2362
delete [] transform;
2363
break;
2364
}
2365
}
2366
2367
}
2368
2369
/// Evaluate order n derivatives of all basis functions at given point x in cell
2370
virtual void evaluate_basis_derivatives_all(std::size_t n,
2371
double* values,
2372
const double* x,
2373
const double* vertex_coordinates,
2374
int cell_orientation) const
2375
{
2376
// Compute number of derivatives.
2377
unsigned int num_derivatives = 1;
2378
for (unsigned int r = 0; r < n; r++)
2379
{
2380
num_derivatives *= 2;
2381
}// end loop over 'r'
2382
2383
// Helper variable to hold values of a single dof.
2384
double *dof_values = new double[num_derivatives];
2385
for (unsigned int r = 0; r < num_derivatives; r++)
2386
{
2387
dof_values[r] = 0.0;
2388
}// end loop over 'r'
2389
2390
// Loop dofs and call evaluate_basis_derivatives.
2391
for (unsigned int r = 0; r < 3; r++)
2392
{
2393
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2394
for (unsigned int s = 0; s < num_derivatives; s++)
2395
{
2396
values[r*num_derivatives + s] = dof_values[s];
2397
}// end loop over 's'
2398
}// end loop over 'r'
2399
2400
// Delete pointer.
2401
delete [] dof_values;
2402
}
2403
2404
/// Evaluate linear functional for dof i on the function f
2405
virtual double evaluate_dof(std::size_t i,
2406
const ufc::function& f,
2407
const double* vertex_coordinates,
2408
int cell_orientation,
2409
const ufc::cell& c) const
2410
{
2411
// Declare variables for result of evaluation
2412
double vals[1];
2413
2414
// Declare variable for physical coordinates
2415
double y[2];
2416
switch (i)
2417
{
2418
case 0:
2419
{
2420
y[0] = vertex_coordinates[0];
2421
y[1] = vertex_coordinates[1];
2422
f.evaluate(vals, y, c);
2423
return vals[0];
2424
break;
2425
}
2426
case 1:
2427
{
2428
y[0] = vertex_coordinates[2];
2429
y[1] = vertex_coordinates[3];
2430
f.evaluate(vals, y, c);
2431
return vals[0];
2432
break;
2433
}
2434
case 2:
2435
{
2436
y[0] = vertex_coordinates[4];
2437
y[1] = vertex_coordinates[5];
2438
f.evaluate(vals, y, c);
2439
return vals[0];
2440
break;
2441
}
2442
}
2443
2444
return 0.0;
2445
}
2446
2447
/// Evaluate linear functionals for all dofs on the function f
2448
virtual void evaluate_dofs(double* values,
2449
const ufc::function& f,
2450
const double* vertex_coordinates,
2451
int cell_orientation,
2452
const ufc::cell& c) const
2453
{
2454
// Declare variables for result of evaluation
2455
double vals[1];
2456
2457
// Declare variable for physical coordinates
2458
double y[2];
2459
y[0] = vertex_coordinates[0];
2460
y[1] = vertex_coordinates[1];
2461
f.evaluate(vals, y, c);
2462
values[0] = vals[0];
2463
y[0] = vertex_coordinates[2];
2464
y[1] = vertex_coordinates[3];
2465
f.evaluate(vals, y, c);
2466
values[1] = vals[0];
2467
y[0] = vertex_coordinates[4];
2468
y[1] = vertex_coordinates[5];
2469
f.evaluate(vals, y, c);
2470
values[2] = vals[0];
2471
}
2472
2473
/// Interpolate vertex values from dof values
2474
virtual void interpolate_vertex_values(double* vertex_values,
2475
const double* dof_values,
2476
const double* vertex_coordinates,
2477
int cell_orientation,
2478
const ufc::cell& c) const
2479
{
2480
// Evaluate function and change variables
2481
vertex_values[0] = dof_values[0];
2482
vertex_values[1] = dof_values[1];
2483
vertex_values[2] = dof_values[2];
2484
}
2485
2486
/// Map coordinate xhat from reference cell to coordinate x in cell
2487
virtual void map_from_reference_cell(double* x,
2488
const double* xhat,
2489
const ufc::cell& c) const
2490
{
2491
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
2492
}
2493
2494
/// Map from coordinate x in cell to coordinate xhat in reference cell
2495
virtual void map_to_reference_cell(double* xhat,
2496
const double* x,
2497
const ufc::cell& c) const
2498
{
2499
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
2500
}
2501
2502
/// Return the number of sub elements (for a mixed element)
2503
virtual std::size_t num_sub_elements() const
2504
{
2505
return 0;
2506
}
2507
2508
/// Create a new finite element for sub element i (for a mixed element)
2509
virtual ufc::finite_element* create_sub_element(std::size_t i) const
2510
{
2511
return 0;
2512
}
2513
2514
/// Create a new class instance
2515
virtual ufc::finite_element* create() const
2516
{
2517
return new pointmeasure_finite_element_1();
2518
}
2519
2520
};
2521
2522
/// This class defines the interface for a local-to-global mapping of
2523
/// degrees of freedom (dofs).
2524
2525
class pointmeasure_dofmap_0: public ufc::dofmap
2526
{
2527
public:
2528
2529
/// Constructor
2530
pointmeasure_dofmap_0() : ufc::dofmap()
2531
{
2532
// Do nothing
2533
}
2534
2535
/// Destructor
2536
virtual ~pointmeasure_dofmap_0()
2537
{
2538
// Do nothing
2539
}
2540
2541
/// Return a string identifying the dofmap
2542
virtual const char* signature() const
2543
{
2544
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
2545
}
2546
2547
/// Return true iff mesh entities of topological dimension d are needed
2548
virtual bool needs_mesh_entities(std::size_t d) const
2549
{
2550
switch (d)
2551
{
2552
case 0:
2553
{
2554
return true;
2555
break;
2556
}
2557
case 1:
2558
{
2559
return true;
2560
break;
2561
}
2562
case 2:
2563
{
2564
return false;
2565
break;
2566
}
2567
}
2568
2569
return false;
2570
}
2571
2572
/// Return the topological dimension of the associated cell shape
2573
virtual std::size_t topological_dimension() const
2574
{
2575
return 2;
2576
}
2577
2578
/// Return the geometric dimension of the associated cell shape
2579
virtual std::size_t geometric_dimension() const
2580
{
2581
return 2;
2582
}
2583
2584
/// Return the dimension of the global finite element function space
2585
virtual std::size_t global_dimension(const std::vector<std::size_t>&
2586
num_global_entities) const
2587
{
2588
return num_global_entities[0] + num_global_entities[1];
2589
}
2590
2591
/// Return the dimension of the local finite element function space for a cell
2592
virtual std::size_t local_dimension(const ufc::cell& c) const
2593
{
2594
return 6;
2595
}
2596
2597
/// Return the maximum dimension of the local finite element function space
2598
virtual std::size_t max_local_dimension() const
2599
{
2600
return 6;
2601
}
2602
2603
/// Return the number of dofs on each cell facet
2604
virtual std::size_t num_facet_dofs() const
2605
{
2606
return 3;
2607
}
2608
2609
/// Return the number of dofs associated with each cell entity of dimension d
2610
virtual std::size_t num_entity_dofs(std::size_t d) const
2611
{
2612
switch (d)
2613
{
2614
case 0:
2615
{
2616
return 1;
2617
break;
2618
}
2619
case 1:
2620
{
2621
return 1;
2622
break;
2623
}
2624
case 2:
2625
{
2626
return 0;
2627
break;
2628
}
2629
}
2630
2631
return 0;
2632
}
2633
2634
/// Tabulate the local-to-global mapping of dofs on a cell
2635
virtual void tabulate_dofs(std::size_t* dofs,
2636
const std::vector<std::size_t>& num_global_entities,
2637
const ufc::cell& c) const
2638
{
2639
unsigned int offset = 0;
2640
dofs[0] = offset + c.entity_indices[0][0];
2641
dofs[1] = offset + c.entity_indices[0][1];
2642
dofs[2] = offset + c.entity_indices[0][2];
2643
offset += num_global_entities[0];
2644
dofs[3] = offset + c.entity_indices[1][0];
2645
dofs[4] = offset + c.entity_indices[1][1];
2646
dofs[5] = offset + c.entity_indices[1][2];
2647
offset += num_global_entities[1];
2648
}
2649
2650
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2651
virtual void tabulate_facet_dofs(std::size_t* dofs,
2652
std::size_t facet) const
2653
{
2654
switch (facet)
2655
{
2656
case 0:
2657
{
2658
dofs[0] = 1;
2659
dofs[1] = 2;
2660
dofs[2] = 3;
2661
break;
2662
}
2663
case 1:
2664
{
2665
dofs[0] = 0;
2666
dofs[1] = 2;
2667
dofs[2] = 4;
2668
break;
2669
}
2670
case 2:
2671
{
2672
dofs[0] = 0;
2673
dofs[1] = 1;
2674
dofs[2] = 5;
2675
break;
2676
}
2677
}
2678
2679
}
2680
2681
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2682
virtual void tabulate_entity_dofs(std::size_t* dofs,
2683
std::size_t d, std::size_t i) const
2684
{
2685
if (d > 2)
2686
{
2687
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2688
}
2689
2690
switch (d)
2691
{
2692
case 0:
2693
{
2694
if (i > 2)
2695
{
2696
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2697
}
2698
2699
switch (i)
2700
{
2701
case 0:
2702
{
2703
dofs[0] = 0;
2704
break;
2705
}
2706
case 1:
2707
{
2708
dofs[0] = 1;
2709
break;
2710
}
2711
case 2:
2712
{
2713
dofs[0] = 2;
2714
break;
2715
}
2716
}
2717
2718
break;
2719
}
2720
case 1:
2721
{
2722
if (i > 2)
2723
{
2724
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2725
}
2726
2727
switch (i)
2728
{
2729
case 0:
2730
{
2731
dofs[0] = 3;
2732
break;
2733
}
2734
case 1:
2735
{
2736
dofs[0] = 4;
2737
break;
2738
}
2739
case 2:
2740
{
2741
dofs[0] = 5;
2742
break;
2743
}
2744
}
2745
2746
break;
2747
}
2748
case 2:
2749
{
2750
2751
break;
2752
}
2753
}
2754
2755
}
2756
2757
/// Tabulate the coordinates of all dofs on a cell
2758
virtual void tabulate_coordinates(double** dof_coordinates,
2759
const double* vertex_coordinates) const
2760
{
2761
dof_coordinates[0][0] = vertex_coordinates[0];
2762
dof_coordinates[0][1] = vertex_coordinates[1];
2763
dof_coordinates[1][0] = vertex_coordinates[2];
2764
dof_coordinates[1][1] = vertex_coordinates[3];
2765
dof_coordinates[2][0] = vertex_coordinates[4];
2766
dof_coordinates[2][1] = vertex_coordinates[5];
2767
dof_coordinates[3][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
2768
dof_coordinates[3][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
2769
dof_coordinates[4][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
2770
dof_coordinates[4][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
2771
dof_coordinates[5][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
2772
dof_coordinates[5][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
2773
}
2774
2775
/// Return the number of sub dofmaps (for a mixed element)
2776
virtual std::size_t num_sub_dofmaps() const
2777
{
2778
return 0;
2779
}
2780
2781
/// Create a new dofmap for sub dofmap i (for a mixed element)
2782
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
2783
{
2784
return 0;
2785
}
2786
2787
/// Create a new class instance
2788
virtual ufc::dofmap* create() const
2789
{
2790
return new pointmeasure_dofmap_0();
2791
}
2792
2793
};
2794
2795
/// This class defines the interface for a local-to-global mapping of
2796
/// degrees of freedom (dofs).
2797
2798
class pointmeasure_dofmap_1: public ufc::dofmap
2799
{
2800
public:
2801
2802
/// Constructor
2803
pointmeasure_dofmap_1() : ufc::dofmap()
2804
{
2805
// Do nothing
2806
}
2807
2808
/// Destructor
2809
virtual ~pointmeasure_dofmap_1()
2810
{
2811
// Do nothing
2812
}
2813
2814
/// Return a string identifying the dofmap
2815
virtual const char* signature() const
2816
{
2817
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
2818
}
2819
2820
/// Return true iff mesh entities of topological dimension d are needed
2821
virtual bool needs_mesh_entities(std::size_t d) const
2822
{
2823
switch (d)
2824
{
2825
case 0:
2826
{
2827
return true;
2828
break;
2829
}
2830
case 1:
2831
{
2832
return false;
2833
break;
2834
}
2835
case 2:
2836
{
2837
return false;
2838
break;
2839
}
2840
}
2841
2842
return false;
2843
}
2844
2845
/// Return the topological dimension of the associated cell shape
2846
virtual std::size_t topological_dimension() const
2847
{
2848
return 2;
2849
}
2850
2851
/// Return the geometric dimension of the associated cell shape
2852
virtual std::size_t geometric_dimension() const
2853
{
2854
return 2;
2855
}
2856
2857
/// Return the dimension of the global finite element function space
2858
virtual std::size_t global_dimension(const std::vector<std::size_t>&
2859
num_global_entities) const
2860
{
2861
return num_global_entities[0];
2862
}
2863
2864
/// Return the dimension of the local finite element function space for a cell
2865
virtual std::size_t local_dimension(const ufc::cell& c) const
2866
{
2867
return 3;
2868
}
2869
2870
/// Return the maximum dimension of the local finite element function space
2871
virtual std::size_t max_local_dimension() const
2872
{
2873
return 3;
2874
}
2875
2876
/// Return the number of dofs on each cell facet
2877
virtual std::size_t num_facet_dofs() const
2878
{
2879
return 2;
2880
}
2881
2882
/// Return the number of dofs associated with each cell entity of dimension d
2883
virtual std::size_t num_entity_dofs(std::size_t d) const
2884
{
2885
switch (d)
2886
{
2887
case 0:
2888
{
2889
return 1;
2890
break;
2891
}
2892
case 1:
2893
{
2894
return 0;
2895
break;
2896
}
2897
case 2:
2898
{
2899
return 0;
2900
break;
2901
}
2902
}
2903
2904
return 0;
2905
}
2906
2907
/// Tabulate the local-to-global mapping of dofs on a cell
2908
virtual void tabulate_dofs(std::size_t* dofs,
2909
const std::vector<std::size_t>& num_global_entities,
2910
const ufc::cell& c) const
2911
{
2912
dofs[0] = c.entity_indices[0][0];
2913
dofs[1] = c.entity_indices[0][1];
2914
dofs[2] = c.entity_indices[0][2];
2915
}
2916
2917
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2918
virtual void tabulate_facet_dofs(std::size_t* dofs,
2919
std::size_t facet) const
2920
{
2921
switch (facet)
2922
{
2923
case 0:
2924
{
2925
dofs[0] = 1;
2926
dofs[1] = 2;
2927
break;
2928
}
2929
case 1:
2930
{
2931
dofs[0] = 0;
2932
dofs[1] = 2;
2933
break;
2934
}
2935
case 2:
2936
{
2937
dofs[0] = 0;
2938
dofs[1] = 1;
2939
break;
2940
}
2941
}
2942
2943
}
2944
2945
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2946
virtual void tabulate_entity_dofs(std::size_t* dofs,
2947
std::size_t d, std::size_t i) const
2948
{
2949
if (d > 2)
2950
{
2951
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2952
}
2953
2954
switch (d)
2955
{
2956
case 0:
2957
{
2958
if (i > 2)
2959
{
2960
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2961
}
2962
2963
switch (i)
2964
{
2965
case 0:
2966
{
2967
dofs[0] = 0;
2968
break;
2969
}
2970
case 1:
2971
{
2972
dofs[0] = 1;
2973
break;
2974
}
2975
case 2:
2976
{
2977
dofs[0] = 2;
2978
break;
2979
}
2980
}
2981
2982
break;
2983
}
2984
case 1:
2985
{
2986
2987
break;
2988
}
2989
case 2:
2990
{
2991
2992
break;
2993
}
2994
}
2995
2996
}
2997
2998
/// Tabulate the coordinates of all dofs on a cell
2999
virtual void tabulate_coordinates(double** dof_coordinates,
3000
const double* vertex_coordinates) const
3001
{
3002
dof_coordinates[0][0] = vertex_coordinates[0];
3003
dof_coordinates[0][1] = vertex_coordinates[1];
3004
dof_coordinates[1][0] = vertex_coordinates[2];
3005
dof_coordinates[1][1] = vertex_coordinates[3];
3006
dof_coordinates[2][0] = vertex_coordinates[4];
3007
dof_coordinates[2][1] = vertex_coordinates[5];
3008
}
3009
3010
/// Return the number of sub dofmaps (for a mixed element)
3011
virtual std::size_t num_sub_dofmaps() const
3012
{
3013
return 0;
3014
}
3015
3016
/// Create a new dofmap for sub dofmap i (for a mixed element)
3017
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3018
{
3019
return 0;
3020
}
3021
3022
/// Create a new class instance
3023
virtual ufc::dofmap* create() const
3024
{
3025
return new pointmeasure_dofmap_1();
3026
}
3027
3028
};
3029
3030
/// This class defines the interface for the tabulation of the cell
3031
/// tensor corresponding to the local contribution to a form from
3032
/// the integral over a cell.
3033
3034
class pointmeasure_cell_integral_0_otherwise: public ufc::cell_integral
3035
{
3036
public:
3037
3038
/// Constructor
3039
pointmeasure_cell_integral_0_otherwise() : ufc::cell_integral()
3040
{
3041
// Do nothing
3042
}
3043
3044
/// Destructor
3045
virtual ~pointmeasure_cell_integral_0_otherwise()
3046
{
3047
// Do nothing
3048
}
3049
3050
/// Tabulate the tensor for the contribution from a local cell
3051
virtual void tabulate_tensor(double* A,
3052
const double * const * w,
3053
const double* vertex_coordinates,
3054
int cell_orientation) const
3055
{
3056
// Number of operations (multiply-add pairs) for Jacobian data: 3
3057
// Number of operations (multiply-add pairs) for geometry tensor: 0
3058
// Number of operations (multiply-add pairs) for tensor contraction: 4
3059
// Total number of operations (multiply-add pairs): 7
3060
3061
// Compute Jacobian
3062
double J[4];
3063
compute_jacobian_triangle_2d(J, vertex_coordinates);
3064
3065
// Compute Jacobian inverse and determinant
3066
double K[4];
3067
double detJ;
3068
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3069
3070
// Set scale factor
3071
const double det = std::abs(detJ);
3072
3073
// Compute geometry tensor
3074
const double G0_ = det;
3075
3076
// Compute element tensor
3077
A[0] = 0.083333333*G0_;
3078
A[1] = 0.041666667*G0_;
3079
A[2] = 0.041666667*G0_;
3080
A[3] = 0.041666667*G0_;
3081
A[4] = 0.083333333*G0_;
3082
A[5] = 0.041666667*G0_;
3083
A[6] = 0.041666667*G0_;
3084
A[7] = 0.041666667*G0_;
3085
A[8] = 0.083333333*G0_;
3086
}
3087
3088
};
3089
3090
/// This class defines the interface for the tabulation of
3091
/// an expression evaluated at exactly one point.
3092
3093
class pointmeasure_point_integral_0_1: public ufc::point_integral
3094
{
3095
public:
3096
3097
/// Constructor
3098
pointmeasure_point_integral_0_1() : ufc::point_integral()
3099
{
3100
// Do nothing
3101
}
3102
3103
/// Destructor
3104
virtual ~pointmeasure_point_integral_0_1()
3105
{
3106
// Do nothing
3107
}
3108
3109
/// Tabulate the tensor for the contribution from the local vertex
3110
virtual void tabulate_tensor(double* A,
3111
const double * const * w,
3112
const double* vertex_coordinates,
3113
std::size_t vertex) const
3114
{
3115
// Compute Jacobian
3116
double J[4];
3117
compute_jacobian_triangle_2d(J, vertex_coordinates);
3118
3119
// Compute Jacobian inverse and determinant
3120
double K[4];
3121
double detJ;
3122
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3123
3124
3125
3126
// Get vertices on edge
3127
3128
// Compute scale factor (length of edge scaled by length of reference interval)
3129
3130
3131
// Array of quadrature weights.
3132
static const double W1 = 1.0;
3133
// Quadrature points on the UFC reference element: ()
3134
3135
// Value of basis functions at quadrature points.
3136
static const double FE0_v0[1][3] = \
3137
{{1.0, 0.0, 0.0}};
3138
3139
static const double FE0_v1[1][3] = \
3140
{{0.0, 1.0, 0.0}};
3141
3142
static const double FE0_v2[1][3] = \
3143
{{0.0, 0.0, 1.0}};
3144
3145
// Reset values in the element tensor.
3146
for (unsigned int r = 0; r < 9; r++)
3147
{
3148
A[r] = 0.0;
3149
}// end loop over 'r'
3150
3151
// Compute element tensor using UFL quadrature representation
3152
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3153
switch (vertex)
3154
{
3155
case 0:
3156
{
3157
// Total number of operations to compute element tensor (from this point): 69
3158
3159
// Loop quadrature points for integral.
3160
// Number of operations to compute element tensor for following IP loop = 69
3161
// Only 1 integration point, omitting IP loop.
3162
3163
// Coefficient declarations.
3164
double F0 = 0.0;
3165
3166
// Total number of operations to compute function values = 6
3167
for (unsigned int r = 0; r < 3; r++)
3168
{
3169
F0 += FE0_v0[0][r]*w[0][r];
3170
}// end loop over 'r'
3171
3172
// Number of operations for primary indices: 63
3173
for (unsigned int j = 0; j < 3; j++)
3174
{
3175
for (unsigned int k = 0; k < 3; k++)
3176
{
3177
// Number of operations to compute entry: 7
3178
A[j*3 + k] += (FE0_v0[0][j]*FE0_v0[0][k] + FE0_v0[0][j]*FE0_v0[0][k]*F0*F0)*W1;
3179
}// end loop over 'k'
3180
}// end loop over 'j'
3181
break;
3182
}
3183
case 1:
3184
{
3185
// Total number of operations to compute element tensor (from this point): 69
3186
3187
// Loop quadrature points for integral.
3188
// Number of operations to compute element tensor for following IP loop = 69
3189
// Only 1 integration point, omitting IP loop.
3190
3191
// Coefficient declarations.
3192
double F0 = 0.0;
3193
3194
// Total number of operations to compute function values = 6
3195
for (unsigned int r = 0; r < 3; r++)
3196
{
3197
F0 += FE0_v1[0][r]*w[0][r];
3198
}// end loop over 'r'
3199
3200
// Number of operations for primary indices: 63
3201
for (unsigned int j = 0; j < 3; j++)
3202
{
3203
for (unsigned int k = 0; k < 3; k++)
3204
{
3205
// Number of operations to compute entry: 7
3206
A[j*3 + k] += (FE0_v1[0][j]*FE0_v1[0][k]*F0*F0 + FE0_v1[0][j]*FE0_v1[0][k])*W1;
3207
}// end loop over 'k'
3208
}// end loop over 'j'
3209
break;
3210
}
3211
case 2:
3212
{
3213
// Total number of operations to compute element tensor (from this point): 69
3214
3215
// Loop quadrature points for integral.
3216
// Number of operations to compute element tensor for following IP loop = 69
3217
// Only 1 integration point, omitting IP loop.
3218
3219
// Coefficient declarations.
3220
double F0 = 0.0;
3221
3222
// Total number of operations to compute function values = 6
3223
for (unsigned int r = 0; r < 3; r++)
3224
{
3225
F0 += FE0_v2[0][r]*w[0][r];
3226
}// end loop over 'r'
3227
3228
// Number of operations for primary indices: 63
3229
for (unsigned int j = 0; j < 3; j++)
3230
{
3231
for (unsigned int k = 0; k < 3; k++)
3232
{
3233
// Number of operations to compute entry: 7
3234
A[j*3 + k] += (FE0_v2[0][j]*FE0_v2[0][k]*F0*F0 + FE0_v2[0][j]*FE0_v2[0][k])*W1;
3235
}// end loop over 'k'
3236
}// end loop over 'j'
3237
break;
3238
}
3239
}
3240
3241
}
3242
3243
};
3244
3245
/// This class defines the interface for the tabulation of
3246
/// an expression evaluated at exactly one point.
3247
3248
class pointmeasure_point_integral_0_otherwise: public ufc::point_integral
3249
{
3250
public:
3251
3252
/// Constructor
3253
pointmeasure_point_integral_0_otherwise() : ufc::point_integral()
3254
{
3255
// Do nothing
3256
}
3257
3258
/// Destructor
3259
virtual ~pointmeasure_point_integral_0_otherwise()
3260
{
3261
// Do nothing
3262
}
3263
3264
/// Tabulate the tensor for the contribution from the local vertex
3265
virtual void tabulate_tensor(double* A,
3266
const double * const * w,
3267
const double* vertex_coordinates,
3268
std::size_t vertex) const
3269
{
3270
// Compute Jacobian
3271
double J[4];
3272
compute_jacobian_triangle_2d(J, vertex_coordinates);
3273
3274
// Compute Jacobian inverse and determinant
3275
double K[4];
3276
double detJ;
3277
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3278
3279
3280
3281
// Get vertices on edge
3282
3283
// Compute scale factor (length of edge scaled by length of reference interval)
3284
3285
3286
// Array of quadrature weights.
3287
static const double W1 = 1.0;
3288
// Quadrature points on the UFC reference element: ()
3289
3290
// Value of basis functions at quadrature points.
3291
static const double FE0_v0[1][3] = \
3292
{{1.0, 0.0, 0.0}};
3293
3294
static const double FE0_v1[1][3] = \
3295
{{0.0, 1.0, 0.0}};
3296
3297
static const double FE0_v2[1][3] = \
3298
{{0.0, 0.0, 1.0}};
3299
3300
// Reset values in the element tensor.
3301
for (unsigned int r = 0; r < 9; r++)
3302
{
3303
A[r] = 0.0;
3304
}// end loop over 'r'
3305
3306
// Compute element tensor using UFL quadrature representation
3307
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3308
switch (vertex)
3309
{
3310
case 0:
3311
{
3312
// Total number of operations to compute element tensor (from this point): 27
3313
3314
// Loop quadrature points for integral.
3315
// Number of operations to compute element tensor for following IP loop = 27
3316
// Only 1 integration point, omitting IP loop.
3317
3318
// Number of operations for primary indices: 27
3319
for (unsigned int j = 0; j < 3; j++)
3320
{
3321
for (unsigned int k = 0; k < 3; k++)
3322
{
3323
// Number of operations to compute entry: 3
3324
A[j*3 + k] += FE0_v0[0][j]*FE0_v0[0][k]*W1;
3325
}// end loop over 'k'
3326
}// end loop over 'j'
3327
break;
3328
}
3329
case 1:
3330
{
3331
// Total number of operations to compute element tensor (from this point): 27
3332
3333
// Loop quadrature points for integral.
3334
// Number of operations to compute element tensor for following IP loop = 27
3335
// Only 1 integration point, omitting IP loop.
3336
3337
// Number of operations for primary indices: 27
3338
for (unsigned int j = 0; j < 3; j++)
3339
{
3340
for (unsigned int k = 0; k < 3; k++)
3341
{
3342
// Number of operations to compute entry: 3
3343
A[j*3 + k] += FE0_v1[0][j]*FE0_v1[0][k]*W1;
3344
}// end loop over 'k'
3345
}// end loop over 'j'
3346
break;
3347
}
3348
case 2:
3349
{
3350
// Total number of operations to compute element tensor (from this point): 27
3351
3352
// Loop quadrature points for integral.
3353
// Number of operations to compute element tensor for following IP loop = 27
3354
// Only 1 integration point, omitting IP loop.
3355
3356
// Number of operations for primary indices: 27
3357
for (unsigned int j = 0; j < 3; j++)
3358
{
3359
for (unsigned int k = 0; k < 3; k++)
3360
{
3361
// Number of operations to compute entry: 3
3362
A[j*3 + k] += FE0_v2[0][j]*FE0_v2[0][k]*W1;
3363
}// end loop over 'k'
3364
}// end loop over 'j'
3365
break;
3366
}
3367
}
3368
3369
}
3370
3371
};
3372
3373
/// This class defines the interface for the tabulation of
3374
/// an expression evaluated at exactly one point.
3375
3376
class pointmeasure_point_integral_1_otherwise: public ufc::point_integral
3377
{
3378
public:
3379
3380
/// Constructor
3381
pointmeasure_point_integral_1_otherwise() : ufc::point_integral()
3382
{
3383
// Do nothing
3384
}
3385
3386
/// Destructor
3387
virtual ~pointmeasure_point_integral_1_otherwise()
3388
{
3389
// Do nothing
3390
}
3391
3392
/// Tabulate the tensor for the contribution from the local vertex
3393
virtual void tabulate_tensor(double* A,
3394
const double * const * w,
3395
const double* vertex_coordinates,
3396
std::size_t vertex) const
3397
{
3398
// Compute Jacobian
3399
double J[4];
3400
compute_jacobian_triangle_2d(J, vertex_coordinates);
3401
3402
// Compute Jacobian inverse and determinant
3403
double K[4];
3404
double detJ;
3405
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3406
3407
3408
3409
// Get vertices on edge
3410
3411
// Compute scale factor (length of edge scaled by length of reference interval)
3412
3413
3414
// Array of quadrature weights.
3415
static const double W1 = 1.0;
3416
// Quadrature points on the UFC reference element: ()
3417
3418
// Value of basis functions at quadrature points.
3419
static const double FE0_v0[1][3] = \
3420
{{1.0, 0.0, 0.0}};
3421
3422
static const double FE0_v1[1][3] = \
3423
{{0.0, 1.0, 0.0}};
3424
3425
static const double FE0_v2[1][3] = \
3426
{{0.0, 0.0, 1.0}};
3427
3428
static const double FE1_v0[1][6] = \
3429
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
3430
3431
static const double FE1_v1[1][6] = \
3432
{{0.0, 1.0, 0.0, 0.0, 0.0, 0.0}};
3433
3434
static const double FE1_v2[1][6] = \
3435
{{0.0, 0.0, 1.0, 0.0, 0.0, 0.0}};
3436
3437
// Reset values in the element tensor.
3438
for (unsigned int r = 0; r < 3; r++)
3439
{
3440
A[r] = 0.0;
3441
}// end loop over 'r'
3442
3443
// Compute element tensor using UFL quadrature representation
3444
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3445
switch (vertex)
3446
{
3447
case 0:
3448
{
3449
// Total number of operations to compute element tensor (from this point): 21
3450
3451
// Loop quadrature points for integral.
3452
// Number of operations to compute element tensor for following IP loop = 21
3453
// Only 1 integration point, omitting IP loop.
3454
3455
// Coefficient declarations.
3456
double F0 = 0.0;
3457
3458
// Total number of operations to compute function values = 12
3459
for (unsigned int r = 0; r < 6; r++)
3460
{
3461
F0 += FE1_v0[0][r]*w[0][r];
3462
}// end loop over 'r'
3463
3464
// Number of operations for primary indices: 9
3465
for (unsigned int j = 0; j < 3; j++)
3466
{
3467
// Number of operations to compute entry: 3
3468
A[j] += FE0_v0[0][j]*F0*W1;
3469
}// end loop over 'j'
3470
break;
3471
}
3472
case 1:
3473
{
3474
// Total number of operations to compute element tensor (from this point): 21
3475
3476
// Loop quadrature points for integral.
3477
// Number of operations to compute element tensor for following IP loop = 21
3478
// Only 1 integration point, omitting IP loop.
3479
3480
// Coefficient declarations.
3481
double F0 = 0.0;
3482
3483
// Total number of operations to compute function values = 12
3484
for (unsigned int r = 0; r < 6; r++)
3485
{
3486
F0 += FE1_v1[0][r]*w[0][r];
3487
}// end loop over 'r'
3488
3489
// Number of operations for primary indices: 9
3490
for (unsigned int j = 0; j < 3; j++)
3491
{
3492
// Number of operations to compute entry: 3
3493
A[j] += FE0_v1[0][j]*F0*W1;
3494
}// end loop over 'j'
3495
break;
3496
}
3497
case 2:
3498
{
3499
// Total number of operations to compute element tensor (from this point): 21
3500
3501
// Loop quadrature points for integral.
3502
// Number of operations to compute element tensor for following IP loop = 21
3503
// Only 1 integration point, omitting IP loop.
3504
3505
// Coefficient declarations.
3506
double F0 = 0.0;
3507
3508
// Total number of operations to compute function values = 12
3509
for (unsigned int r = 0; r < 6; r++)
3510
{
3511
F0 += FE1_v2[0][r]*w[0][r];
3512
}// end loop over 'r'
3513
3514
// Number of operations for primary indices: 9
3515
for (unsigned int j = 0; j < 3; j++)
3516
{
3517
// Number of operations to compute entry: 3
3518
A[j] += FE0_v2[0][j]*F0*W1;
3519
}// end loop over 'j'
3520
break;
3521
}
3522
}
3523
3524
}
3525
3526
};
3527
3528
/// This class defines the interface for the assembly of the global
3529
/// tensor corresponding to a form with r + n arguments, that is, a
3530
/// mapping
3531
///
3532
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3533
///
3534
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3535
/// global tensor A is defined by
3536
///
3537
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3538
///
3539
/// where each argument Vj represents the application to the
3540
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3541
/// fixed functions (coefficients).
3542
3543
class pointmeasure_form_0: public ufc::form
3544
{
3545
public:
3546
3547
/// Constructor
3548
pointmeasure_form_0() : ufc::form()
3549
{
3550
// Do nothing
3551
}
3552
3553
/// Destructor
3554
virtual ~pointmeasure_form_0()
3555
{
3556
// Do nothing
3557
}
3558
3559
/// Return a string identifying the form
3560
virtual const char* signature() const
3561
{
3562
return "856f5aa1a8f7253511ba78602090dc3cf582e71ee4c859de2e80d2fd37c92108a283e09a2ae620c72f13f280a38368ed9c5e78e015492db3cbebd6d8f66016fe";
3563
}
3564
3565
/// Return the rank of the global tensor (r)
3566
virtual std::size_t rank() const
3567
{
3568
return 2;
3569
}
3570
3571
/// Return the number of coefficients (n)
3572
virtual std::size_t num_coefficients() const
3573
{
3574
return 1;
3575
}
3576
3577
/// Return the number of cell domains
3578
virtual std::size_t num_cell_domains() const
3579
{
3580
return 0;
3581
}
3582
3583
/// Return the number of exterior facet domains
3584
virtual std::size_t num_exterior_facet_domains() const
3585
{
3586
return 0;
3587
}
3588
3589
/// Return the number of interior facet domains
3590
virtual std::size_t num_interior_facet_domains() const
3591
{
3592
return 0;
3593
}
3594
3595
/// Return the number of point domains
3596
virtual std::size_t num_point_domains() const
3597
{
3598
return 2;
3599
}
3600
3601
/// Return whether the form has any cell integrals
3602
virtual bool has_cell_integrals() const
3603
{
3604
return true;
3605
}
3606
3607
/// Return whether the form has any exterior facet integrals
3608
virtual bool has_exterior_facet_integrals() const
3609
{
3610
return false;
3611
}
3612
3613
/// Return whether the form has any interior facet integrals
3614
virtual bool has_interior_facet_integrals() const
3615
{
3616
return false;
3617
}
3618
3619
/// Return whether the form has any point integrals
3620
virtual bool has_point_integrals() const
3621
{
3622
return true;
3623
}
3624
3625
/// Create a new finite element for argument function i
3626
virtual ufc::finite_element* create_finite_element(std::size_t i) const
3627
{
3628
switch (i)
3629
{
3630
case 0:
3631
{
3632
return new pointmeasure_finite_element_1();
3633
break;
3634
}
3635
case 1:
3636
{
3637
return new pointmeasure_finite_element_1();
3638
break;
3639
}
3640
case 2:
3641
{
3642
return new pointmeasure_finite_element_1();
3643
break;
3644
}
3645
}
3646
3647
return 0;
3648
}
3649
3650
/// Create a new dofmap for argument function i
3651
virtual ufc::dofmap* create_dofmap(std::size_t i) const
3652
{
3653
switch (i)
3654
{
3655
case 0:
3656
{
3657
return new pointmeasure_dofmap_1();
3658
break;
3659
}
3660
case 1:
3661
{
3662
return new pointmeasure_dofmap_1();
3663
break;
3664
}
3665
case 2:
3666
{
3667
return new pointmeasure_dofmap_1();
3668
break;
3669
}
3670
}
3671
3672
return 0;
3673
}
3674
3675
/// Create a new cell integral on sub domain i
3676
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
3677
{
3678
return 0;
3679
}
3680
3681
/// Create a new exterior facet integral on sub domain i
3682
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
3683
{
3684
return 0;
3685
}
3686
3687
/// Create a new interior facet integral on sub domain i
3688
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
3689
{
3690
return 0;
3691
}
3692
3693
/// Create a new point integral on sub domain i
3694
virtual ufc::point_integral* create_point_integral(std::size_t i) const
3695
{
3696
switch (i)
3697
{
3698
case 1:
3699
{
3700
return new pointmeasure_point_integral_0_1();
3701
break;
3702
}
3703
}
3704
3705
return 0;
3706
}
3707
3708
/// Create a new cell integral on everywhere else
3709
virtual ufc::cell_integral* create_default_cell_integral() const
3710
{
3711
return new pointmeasure_cell_integral_0_otherwise();
3712
}
3713
3714
/// Create a new exterior facet integral on everywhere else
3715
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
3716
{
3717
return 0;
3718
}
3719
3720
/// Create a new interior facet integral on everywhere else
3721
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
3722
{
3723
return 0;
3724
}
3725
3726
/// Create a new point integral on everywhere else
3727
virtual ufc::point_integral* create_default_point_integral() const
3728
{
3729
return new pointmeasure_point_integral_0_otherwise();
3730
}
3731
3732
};
3733
3734
/// This class defines the interface for the assembly of the global
3735
/// tensor corresponding to a form with r + n arguments, that is, a
3736
/// mapping
3737
///
3738
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3739
///
3740
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3741
/// global tensor A is defined by
3742
///
3743
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3744
///
3745
/// where each argument Vj represents the application to the
3746
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3747
/// fixed functions (coefficients).
3748
3749
class pointmeasure_form_1: public ufc::form
3750
{
3751
public:
3752
3753
/// Constructor
3754
pointmeasure_form_1() : ufc::form()
3755
{
3756
// Do nothing
3757
}
3758
3759
/// Destructor
3760
virtual ~pointmeasure_form_1()
3761
{
3762
// Do nothing
3763
}
3764
3765
/// Return a string identifying the form
3766
virtual const char* signature() const
3767
{
3768
return "cbc400c8cb006d3574d6a58223a2c80d90b230d781feaabb65968d87ac0550f09264f2d1f32621512f1d42a7bc048c03d38e4482a1abf0fba8fad9c1cab8cc38";
3769
}
3770
3771
/// Return the rank of the global tensor (r)
3772
virtual std::size_t rank() const
3773
{
3774
return 1;
3775
}
3776
3777
/// Return the number of coefficients (n)
3778
virtual std::size_t num_coefficients() const
3779
{
3780
return 1;
3781
}
3782
3783
/// Return the number of cell domains
3784
virtual std::size_t num_cell_domains() const
3785
{
3786
return 0;
3787
}
3788
3789
/// Return the number of exterior facet domains
3790
virtual std::size_t num_exterior_facet_domains() const
3791
{
3792
return 0;
3793
}
3794
3795
/// Return the number of interior facet domains
3796
virtual std::size_t num_interior_facet_domains() const
3797
{
3798
return 0;
3799
}
3800
3801
/// Return the number of point domains
3802
virtual std::size_t num_point_domains() const
3803
{
3804
return 0;
3805
}
3806
3807
/// Return whether the form has any cell integrals
3808
virtual bool has_cell_integrals() const
3809
{
3810
return false;
3811
}
3812
3813
/// Return whether the form has any exterior facet integrals
3814
virtual bool has_exterior_facet_integrals() const
3815
{
3816
return false;
3817
}
3818
3819
/// Return whether the form has any interior facet integrals
3820
virtual bool has_interior_facet_integrals() const
3821
{
3822
return false;
3823
}
3824
3825
/// Return whether the form has any point integrals
3826
virtual bool has_point_integrals() const
3827
{
3828
return true;
3829
}
3830
3831
/// Create a new finite element for argument function i
3832
virtual ufc::finite_element* create_finite_element(std::size_t i) const
3833
{
3834
switch (i)
3835
{
3836
case 0:
3837
{
3838
return new pointmeasure_finite_element_1();
3839
break;
3840
}
3841
case 1:
3842
{
3843
return new pointmeasure_finite_element_0();
3844
break;
3845
}
3846
}
3847
3848
return 0;
3849
}
3850
3851
/// Create a new dofmap for argument function i
3852
virtual ufc::dofmap* create_dofmap(std::size_t i) const
3853
{
3854
switch (i)
3855
{
3856
case 0:
3857
{
3858
return new pointmeasure_dofmap_1();
3859
break;
3860
}
3861
case 1:
3862
{
3863
return new pointmeasure_dofmap_0();
3864
break;
3865
}
3866
}
3867
3868
return 0;
3869
}
3870
3871
/// Create a new cell integral on sub domain i
3872
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
3873
{
3874
return 0;
3875
}
3876
3877
/// Create a new exterior facet integral on sub domain i
3878
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
3879
{
3880
return 0;
3881
}
3882
3883
/// Create a new interior facet integral on sub domain i
3884
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
3885
{
3886
return 0;
3887
}
3888
3889
/// Create a new point integral on sub domain i
3890
virtual ufc::point_integral* create_point_integral(std::size_t i) const
3891
{
3892
return 0;
3893
}
3894
3895
/// Create a new cell integral on everywhere else
3896
virtual ufc::cell_integral* create_default_cell_integral() const
3897
{
3898
return 0;
3899
}
3900
3901
/// Create a new exterior facet integral on everywhere else
3902
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
3903
{
3904
return 0;
3905
}
3906
3907
/// Create a new interior facet integral on everywhere else
3908
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
3909
{
3910
return 0;
3911
}
3912
3913
/// Create a new point integral on everywhere else
3914
virtual ufc::point_integral* create_default_point_integral() const
3915
{
3916
return new pointmeasure_point_integral_1_otherwise();
3917
}
3918
3919
};
3920
3921
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1