← Back to branch summary

~ubuntu-branches/ubuntu/raring/ffc/raring

« back to all changes in this revision

Viewing changes to test/regression/references/Poisson.h

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2010-07-01 19:54:32 UTC
  • mfrom: (1.1.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20100701195432-xz3gw5nprdj79jcb
Tags: 0.9.3-1
* New upstream release.
* debian/control:
  - Minor fix in Vcs fields.
  - Bump Standards-Version to 3.9.0 (no changes needed).
  - Update version for python-ufc, python-fiat, and python-ufl in
    Depends field.
* Switch to dpkg-source 3.0 (quilt) format.
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/regression/references/Poisson.h
1
 
// This code conforms with the UFC specification version 1.4
2
 
// and was automatically generated by FFC version 0.9.2.
 
1
// This code conforms with the UFC specification version 1.4.1
 
2
// and was automatically generated by FFC version 0.9.3.
3
3
// 
4
4
// This code was generated with the following parameters:
5
5
// 
105
105
    
106
106
    // Reset values.
107
107
    *values = 0.00000000;
108
 
    
109
 
    // Map degree of freedom to element degree of freedom
110
 
    const unsigned int dof = i;
111
 
    
112
 
    // Array of basisvalues.
113
 
    double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
114
 
    
115
 
    // Declare helper variables.
116
 
    unsigned int rr = 0;
117
 
    unsigned int ss = 0;
118
 
    double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
119
 
    
120
 
    // Compute basisvalues.
121
 
    basisvalues[0] = 1.00000000;
122
 
    basisvalues[1] = tmp0;
123
 
    for (unsigned int r = 0; r < 1; r++)
124
 
    {
125
 
      rr = (r + 1)*(r + 1 + 1)/2 + 1;
126
 
      ss = r*(r + 1)/2;
127
 
      basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
128
 
    }// end loop over 'r'
129
 
    for (unsigned int r = 0; r < 2; r++)
130
 
    {
131
 
      for (unsigned int s = 0; s < 2 - r; s++)
132
 
      {
133
 
        rr = (r + s)*(r + s + 1)/2 + s;
134
 
        basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
135
 
      }// end loop over 's'
136
 
    }// end loop over 'r'
137
 
    
138
 
    // Table(s) of coefficients.
139
 
    static const double coefficients0[3][3] = \
140
 
    {{0.47140452, -0.28867513, -0.16666667},
141
 
    {0.47140452, 0.28867513, -0.16666667},
142
 
    {0.47140452, 0.00000000, 0.33333333}};
143
 
    
144
 
    // Compute value(s).
145
 
    for (unsigned int r = 0; r < 3; r++)
146
 
    {
147
 
      *values += coefficients0[dof][r]*basisvalues[r];
148
 
    }// end loop over 'r'
 
108
    switch (i)
 
109
    {
 
110
    case 0:
 
111
      {
 
112
        
 
113
      // Array of basisvalues.
 
114
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
115
      
 
116
      // Declare helper variables.
 
117
      unsigned int rr = 0;
 
118
      unsigned int ss = 0;
 
119
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
120
      
 
121
      // Compute basisvalues.
 
122
      basisvalues[0] = 1.00000000;
 
123
      basisvalues[1] = tmp0;
 
124
      for (unsigned int r = 0; r < 1; r++)
 
125
      {
 
126
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
127
        ss = r*(r + 1)/2;
 
128
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
129
      }// end loop over 'r'
 
130
      for (unsigned int r = 0; r < 2; r++)
 
131
      {
 
132
        for (unsigned int s = 0; s < 2 - r; s++)
 
133
        {
 
134
          rr = (r + s)*(r + s + 1)/2 + s;
 
135
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
136
        }// end loop over 's'
 
137
      }// end loop over 'r'
 
138
      
 
139
      // Table(s) of coefficients.
 
140
      static const double coefficients0[3] = \
 
141
      {0.47140452, -0.28867513, -0.16666667};
 
142
      
 
143
      // Compute value(s).
 
144
      for (unsigned int r = 0; r < 3; r++)
 
145
      {
 
146
        *values += coefficients0[r]*basisvalues[r];
 
147
      }// end loop over 'r'
 
148
        break;
 
149
      }
 
150
    case 1:
 
151
      {
 
152
        
 
153
      // Array of basisvalues.
 
154
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
155
      
 
156
      // Declare helper variables.
 
157
      unsigned int rr = 0;
 
158
      unsigned int ss = 0;
 
159
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
160
      
 
161
      // Compute basisvalues.
 
162
      basisvalues[0] = 1.00000000;
 
163
      basisvalues[1] = tmp0;
 
164
      for (unsigned int r = 0; r < 1; r++)
 
165
      {
 
166
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
167
        ss = r*(r + 1)/2;
 
168
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
169
      }// end loop over 'r'
 
170
      for (unsigned int r = 0; r < 2; r++)
 
171
      {
 
172
        for (unsigned int s = 0; s < 2 - r; s++)
 
173
        {
 
174
          rr = (r + s)*(r + s + 1)/2 + s;
 
175
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
176
        }// end loop over 's'
 
177
      }// end loop over 'r'
 
178
      
 
179
      // Table(s) of coefficients.
 
180
      static const double coefficients0[3] = \
 
181
      {0.47140452, 0.28867513, -0.16666667};
 
182
      
 
183
      // Compute value(s).
 
184
      for (unsigned int r = 0; r < 3; r++)
 
185
      {
 
186
        *values += coefficients0[r]*basisvalues[r];
 
187
      }// end loop over 'r'
 
188
        break;
 
189
      }
 
190
    case 2:
 
191
      {
 
192
        
 
193
      // Array of basisvalues.
 
194
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
195
      
 
196
      // Declare helper variables.
 
197
      unsigned int rr = 0;
 
198
      unsigned int ss = 0;
 
199
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
200
      
 
201
      // Compute basisvalues.
 
202
      basisvalues[0] = 1.00000000;
 
203
      basisvalues[1] = tmp0;
 
204
      for (unsigned int r = 0; r < 1; r++)
 
205
      {
 
206
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
207
        ss = r*(r + 1)/2;
 
208
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
209
      }// end loop over 'r'
 
210
      for (unsigned int r = 0; r < 2; r++)
 
211
      {
 
212
        for (unsigned int s = 0; s < 2 - r; s++)
 
213
        {
 
214
          rr = (r + s)*(r + s + 1)/2 + s;
 
215
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
216
        }// end loop over 's'
 
217
      }// end loop over 'r'
 
218
      
 
219
      // Table(s) of coefficients.
 
220
      static const double coefficients0[3] = \
 
221
      {0.47140452, 0.00000000, 0.33333333};
 
222
      
 
223
      // Compute value(s).
 
224
      for (unsigned int r = 0; r < 3; r++)
 
225
      {
 
226
        *values += coefficients0[r]*basisvalues[r];
 
227
      }// end loop over 'r'
 
228
        break;
 
229
      }
 
230
    }
 
231
    
149
232
  }
150
233
 
151
234
  /// Evaluate all basis functions at given point in cell
261
344
      values[r] = 0.00000000;
262
345
    }// end loop over 'r'
263
346
    
264
 
    // Map degree of freedom to element degree of freedom
265
 
    const unsigned int dof = i;
266
 
    
267
 
    // Array of basisvalues.
268
 
    double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
269
 
    
270
 
    // Declare helper variables.
271
 
    unsigned int rr = 0;
272
 
    unsigned int ss = 0;
273
 
    double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
274
 
    
275
 
    // Compute basisvalues.
276
 
    basisvalues[0] = 1.00000000;
277
 
    basisvalues[1] = tmp0;
278
 
    for (unsigned int r = 0; r < 1; r++)
279
 
    {
280
 
      rr = (r + 1)*(r + 1 + 1)/2 + 1;
281
 
      ss = r*(r + 1)/2;
282
 
      basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
283
 
    }// end loop over 'r'
284
 
    for (unsigned int r = 0; r < 2; r++)
285
 
    {
286
 
      for (unsigned int s = 0; s < 2 - r; s++)
287
 
      {
288
 
        rr = (r + s)*(r + s + 1)/2 + s;
289
 
        basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
290
 
      }// end loop over 's'
291
 
    }// end loop over 'r'
292
 
    
293
 
    // Table(s) of coefficients.
294
 
    static const double coefficients0[3][3] = \
295
 
    {{0.47140452, -0.28867513, -0.16666667},
296
 
    {0.47140452, 0.28867513, -0.16666667},
297
 
    {0.47140452, 0.00000000, 0.33333333}};
298
 
    
299
 
    // Tables of derivatives of the polynomial base (transpose).
300
 
    static const double dmats0[3][3] = \
301
 
    {{0.00000000, 0.00000000, 0.00000000},
302
 
    {4.89897949, 0.00000000, 0.00000000},
303
 
    {0.00000000, 0.00000000, 0.00000000}};
304
 
    
305
 
    static const double dmats1[3][3] = \
306
 
    {{0.00000000, 0.00000000, 0.00000000},
307
 
    {2.44948974, 0.00000000, 0.00000000},
308
 
    {4.24264069, 0.00000000, 0.00000000}};
309
 
    
310
 
    // Compute reference derivatives.
311
 
    // Declare pointer to array of derivatives on FIAT element.
312
 
    double *derivatives = new double[num_derivatives];
313
 
    for (unsigned int r = 0; r < num_derivatives; r++)
314
 
    {
315
 
      derivatives[r] = 0.00000000;
316
 
    }// end loop over 'r'
317
 
    
318
 
    // Declare derivative matrix (of polynomial basis).
319
 
    double dmats[3][3] = \
320
 
    {{1.00000000, 0.00000000, 0.00000000},
321
 
    {0.00000000, 1.00000000, 0.00000000},
322
 
    {0.00000000, 0.00000000, 1.00000000}};
323
 
    
324
 
    // Declare (auxiliary) derivative matrix (of polynomial basis).
325
 
    double dmats_old[3][3] = \
326
 
    {{1.00000000, 0.00000000, 0.00000000},
327
 
    {0.00000000, 1.00000000, 0.00000000},
328
 
    {0.00000000, 0.00000000, 1.00000000}};
329
 
    
330
 
    // Loop possible derivatives.
331
 
    for (unsigned int r = 0; r < num_derivatives; r++)
332
 
    {
333
 
      // Resetting dmats values to compute next derivative.
334
 
      for (unsigned int t = 0; t < 3; t++)
335
 
      {
336
 
        for (unsigned int u = 0; u < 3; u++)
337
 
        {
338
 
          dmats[t][u] = 0.00000000;
339
 
          if (t == u)
340
 
          {
341
 
          dmats[t][u] = 1.00000000;
342
 
          }
343
 
          
344
 
        }// end loop over 'u'
345
 
      }// end loop over 't'
346
 
      
347
 
      // Looping derivative order to generate dmats.
348
 
      for (unsigned int s = 0; s < n; s++)
349
 
      {
350
 
        // Updating dmats_old with new values and resetting dmats.
351
 
        for (unsigned int t = 0; t < 3; t++)
352
 
        {
353
 
          for (unsigned int u = 0; u < 3; u++)
354
 
          {
355
 
            dmats_old[t][u] = dmats[t][u];
356
 
            dmats[t][u] = 0.00000000;
357
 
          }// end loop over 'u'
358
 
        }// end loop over 't'
359
 
        
360
 
        // Update dmats using an inner product.
361
 
        if (combinations[r][s] == 0)
362
 
        {
363
 
        for (unsigned int t = 0; t < 3; t++)
364
 
        {
365
 
          for (unsigned int u = 0; u < 3; u++)
366
 
          {
367
 
            for (unsigned int tu = 0; tu < 3; tu++)
368
 
            {
369
 
              dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
370
 
            }// end loop over 'tu'
371
 
          }// end loop over 'u'
372
 
        }// end loop over 't'
373
 
        }
374
 
        
375
 
        if (combinations[r][s] == 1)
376
 
        {
377
 
        for (unsigned int t = 0; t < 3; t++)
378
 
        {
379
 
          for (unsigned int u = 0; u < 3; u++)
380
 
          {
381
 
            for (unsigned int tu = 0; tu < 3; tu++)
382
 
            {
383
 
              dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
384
 
            }// end loop over 'tu'
385
 
          }// end loop over 'u'
386
 
        }// end loop over 't'
387
 
        }
388
 
        
389
 
      }// end loop over 's'
390
 
      for (unsigned int s = 0; s < 3; s++)
391
 
      {
392
 
        for (unsigned int t = 0; t < 3; t++)
393
 
        {
394
 
          derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
395
 
        }// end loop over 't'
396
 
      }// end loop over 's'
397
 
    }// end loop over 'r'
398
 
    
399
 
    // Transform derivatives back to physical element
400
 
    for (unsigned int r = 0; r < num_derivatives; r++)
401
 
    {
402
 
      for (unsigned int s = 0; s < num_derivatives; s++)
403
 
      {
404
 
        values[r] += transform[r][s]*derivatives[s];
405
 
      }// end loop over 's'
406
 
    }// end loop over 'r'
407
 
    
408
 
    // Delete pointer to array of derivatives on FIAT element
409
 
    delete [] derivatives;
410
 
    
411
 
    // Delete pointer to array of combinations of derivatives and transform
412
 
    for (unsigned int r = 0; r < num_derivatives; r++)
413
 
    {
414
 
      delete [] combinations[r];
415
 
    }// end loop over 'r'
416
 
    delete [] combinations;
417
 
    for (unsigned int r = 0; r < num_derivatives; r++)
418
 
    {
419
 
      delete [] transform[r];
420
 
    }// end loop over 'r'
421
 
    delete [] transform;
 
347
    switch (i)
 
348
    {
 
349
    case 0:
 
350
      {
 
351
        
 
352
      // Array of basisvalues.
 
353
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
354
      
 
355
      // Declare helper variables.
 
356
      unsigned int rr = 0;
 
357
      unsigned int ss = 0;
 
358
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
359
      
 
360
      // Compute basisvalues.
 
361
      basisvalues[0] = 1.00000000;
 
362
      basisvalues[1] = tmp0;
 
363
      for (unsigned int r = 0; r < 1; r++)
 
364
      {
 
365
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
366
        ss = r*(r + 1)/2;
 
367
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
368
      }// end loop over 'r'
 
369
      for (unsigned int r = 0; r < 2; r++)
 
370
      {
 
371
        for (unsigned int s = 0; s < 2 - r; s++)
 
372
        {
 
373
          rr = (r + s)*(r + s + 1)/2 + s;
 
374
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
375
        }// end loop over 's'
 
376
      }// end loop over 'r'
 
377
      
 
378
      // Table(s) of coefficients.
 
379
      static const double coefficients0[3] = \
 
380
      {0.47140452, -0.28867513, -0.16666667};
 
381
      
 
382
      // Tables of derivatives of the polynomial base (transpose).
 
383
      static const double dmats0[3][3] = \
 
384
      {{0.00000000, 0.00000000, 0.00000000},
 
385
      {4.89897949, 0.00000000, 0.00000000},
 
386
      {0.00000000, 0.00000000, 0.00000000}};
 
387
      
 
388
      static const double dmats1[3][3] = \
 
389
      {{0.00000000, 0.00000000, 0.00000000},
 
390
      {2.44948974, 0.00000000, 0.00000000},
 
391
      {4.24264069, 0.00000000, 0.00000000}};
 
392
      
 
393
      // Compute reference derivatives.
 
394
      // Declare pointer to array of derivatives on FIAT element.
 
395
      double *derivatives = new double[num_derivatives];
 
396
      for (unsigned int r = 0; r < num_derivatives; r++)
 
397
      {
 
398
        derivatives[r] = 0.00000000;
 
399
      }// end loop over 'r'
 
400
      
 
401
      // Declare derivative matrix (of polynomial basis).
 
402
      double dmats[3][3] = \
 
403
      {{1.00000000, 0.00000000, 0.00000000},
 
404
      {0.00000000, 1.00000000, 0.00000000},
 
405
      {0.00000000, 0.00000000, 1.00000000}};
 
406
      
 
407
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
408
      double dmats_old[3][3] = \
 
409
      {{1.00000000, 0.00000000, 0.00000000},
 
410
      {0.00000000, 1.00000000, 0.00000000},
 
411
      {0.00000000, 0.00000000, 1.00000000}};
 
412
      
 
413
      // Loop possible derivatives.
 
414
      for (unsigned int r = 0; r < num_derivatives; r++)
 
415
      {
 
416
        // Resetting dmats values to compute next derivative.
 
417
        for (unsigned int t = 0; t < 3; t++)
 
418
        {
 
419
          for (unsigned int u = 0; u < 3; u++)
 
420
          {
 
421
            dmats[t][u] = 0.00000000;
 
422
            if (t == u)
 
423
            {
 
424
            dmats[t][u] = 1.00000000;
 
425
            }
 
426
            
 
427
          }// end loop over 'u'
 
428
        }// end loop over 't'
 
429
        
 
430
        // Looping derivative order to generate dmats.
 
431
        for (unsigned int s = 0; s < n; s++)
 
432
        {
 
433
          // Updating dmats_old with new values and resetting dmats.
 
434
          for (unsigned int t = 0; t < 3; t++)
 
435
          {
 
436
            for (unsigned int u = 0; u < 3; u++)
 
437
            {
 
438
              dmats_old[t][u] = dmats[t][u];
 
439
              dmats[t][u] = 0.00000000;
 
440
            }// end loop over 'u'
 
441
          }// end loop over 't'
 
442
          
 
443
          // Update dmats using an inner product.
 
444
          if (combinations[r][s] == 0)
 
445
          {
 
446
          for (unsigned int t = 0; t < 3; t++)
 
447
          {
 
448
            for (unsigned int u = 0; u < 3; u++)
 
449
            {
 
450
              for (unsigned int tu = 0; tu < 3; tu++)
 
451
              {
 
452
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
453
              }// end loop over 'tu'
 
454
            }// end loop over 'u'
 
455
          }// end loop over 't'
 
456
          }
 
457
          
 
458
          if (combinations[r][s] == 1)
 
459
          {
 
460
          for (unsigned int t = 0; t < 3; t++)
 
461
          {
 
462
            for (unsigned int u = 0; u < 3; u++)
 
463
            {
 
464
              for (unsigned int tu = 0; tu < 3; tu++)
 
465
              {
 
466
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
467
              }// end loop over 'tu'
 
468
            }// end loop over 'u'
 
469
          }// end loop over 't'
 
470
          }
 
471
          
 
472
        }// end loop over 's'
 
473
        for (unsigned int s = 0; s < 3; s++)
 
474
        {
 
475
          for (unsigned int t = 0; t < 3; t++)
 
476
          {
 
477
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
478
          }// end loop over 't'
 
479
        }// end loop over 's'
 
480
      }// end loop over 'r'
 
481
      
 
482
      // Transform derivatives back to physical element
 
483
      for (unsigned int r = 0; r < num_derivatives; r++)
 
484
      {
 
485
        for (unsigned int s = 0; s < num_derivatives; s++)
 
486
        {
 
487
          values[r] += transform[r][s]*derivatives[s];
 
488
        }// end loop over 's'
 
489
      }// end loop over 'r'
 
490
      
 
491
      // Delete pointer to array of derivatives on FIAT element
 
492
      delete [] derivatives;
 
493
      
 
494
      // Delete pointer to array of combinations of derivatives and transform
 
495
      for (unsigned int r = 0; r < num_derivatives; r++)
 
496
      {
 
497
        delete [] combinations[r];
 
498
      }// end loop over 'r'
 
499
      delete [] combinations;
 
500
      for (unsigned int r = 0; r < num_derivatives; r++)
 
501
      {
 
502
        delete [] transform[r];
 
503
      }// end loop over 'r'
 
504
      delete [] transform;
 
505
        break;
 
506
      }
 
507
    case 1:
 
508
      {
 
509
        
 
510
      // Array of basisvalues.
 
511
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
512
      
 
513
      // Declare helper variables.
 
514
      unsigned int rr = 0;
 
515
      unsigned int ss = 0;
 
516
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
517
      
 
518
      // Compute basisvalues.
 
519
      basisvalues[0] = 1.00000000;
 
520
      basisvalues[1] = tmp0;
 
521
      for (unsigned int r = 0; r < 1; r++)
 
522
      {
 
523
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
524
        ss = r*(r + 1)/2;
 
525
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
526
      }// end loop over 'r'
 
527
      for (unsigned int r = 0; r < 2; r++)
 
528
      {
 
529
        for (unsigned int s = 0; s < 2 - r; s++)
 
530
        {
 
531
          rr = (r + s)*(r + s + 1)/2 + s;
 
532
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
533
        }// end loop over 's'
 
534
      }// end loop over 'r'
 
535
      
 
536
      // Table(s) of coefficients.
 
537
      static const double coefficients0[3] = \
 
538
      {0.47140452, 0.28867513, -0.16666667};
 
539
      
 
540
      // Tables of derivatives of the polynomial base (transpose).
 
541
      static const double dmats0[3][3] = \
 
542
      {{0.00000000, 0.00000000, 0.00000000},
 
543
      {4.89897949, 0.00000000, 0.00000000},
 
544
      {0.00000000, 0.00000000, 0.00000000}};
 
545
      
 
546
      static const double dmats1[3][3] = \
 
547
      {{0.00000000, 0.00000000, 0.00000000},
 
548
      {2.44948974, 0.00000000, 0.00000000},
 
549
      {4.24264069, 0.00000000, 0.00000000}};
 
550
      
 
551
      // Compute reference derivatives.
 
552
      // Declare pointer to array of derivatives on FIAT element.
 
553
      double *derivatives = new double[num_derivatives];
 
554
      for (unsigned int r = 0; r < num_derivatives; r++)
 
555
      {
 
556
        derivatives[r] = 0.00000000;
 
557
      }// end loop over 'r'
 
558
      
 
559
      // Declare derivative matrix (of polynomial basis).
 
560
      double dmats[3][3] = \
 
561
      {{1.00000000, 0.00000000, 0.00000000},
 
562
      {0.00000000, 1.00000000, 0.00000000},
 
563
      {0.00000000, 0.00000000, 1.00000000}};
 
564
      
 
565
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
566
      double dmats_old[3][3] = \
 
567
      {{1.00000000, 0.00000000, 0.00000000},
 
568
      {0.00000000, 1.00000000, 0.00000000},
 
569
      {0.00000000, 0.00000000, 1.00000000}};
 
570
      
 
571
      // Loop possible derivatives.
 
572
      for (unsigned int r = 0; r < num_derivatives; r++)
 
573
      {
 
574
        // Resetting dmats values to compute next derivative.
 
575
        for (unsigned int t = 0; t < 3; t++)
 
576
        {
 
577
          for (unsigned int u = 0; u < 3; u++)
 
578
          {
 
579
            dmats[t][u] = 0.00000000;
 
580
            if (t == u)
 
581
            {
 
582
            dmats[t][u] = 1.00000000;
 
583
            }
 
584
            
 
585
          }// end loop over 'u'
 
586
        }// end loop over 't'
 
587
        
 
588
        // Looping derivative order to generate dmats.
 
589
        for (unsigned int s = 0; s < n; s++)
 
590
        {
 
591
          // Updating dmats_old with new values and resetting dmats.
 
592
          for (unsigned int t = 0; t < 3; t++)
 
593
          {
 
594
            for (unsigned int u = 0; u < 3; u++)
 
595
            {
 
596
              dmats_old[t][u] = dmats[t][u];
 
597
              dmats[t][u] = 0.00000000;
 
598
            }// end loop over 'u'
 
599
          }// end loop over 't'
 
600
          
 
601
          // Update dmats using an inner product.
 
602
          if (combinations[r][s] == 0)
 
603
          {
 
604
          for (unsigned int t = 0; t < 3; t++)
 
605
          {
 
606
            for (unsigned int u = 0; u < 3; u++)
 
607
            {
 
608
              for (unsigned int tu = 0; tu < 3; tu++)
 
609
              {
 
610
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
611
              }// end loop over 'tu'
 
612
            }// end loop over 'u'
 
613
          }// end loop over 't'
 
614
          }
 
615
          
 
616
          if (combinations[r][s] == 1)
 
617
          {
 
618
          for (unsigned int t = 0; t < 3; t++)
 
619
          {
 
620
            for (unsigned int u = 0; u < 3; u++)
 
621
            {
 
622
              for (unsigned int tu = 0; tu < 3; tu++)
 
623
              {
 
624
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
625
              }// end loop over 'tu'
 
626
            }// end loop over 'u'
 
627
          }// end loop over 't'
 
628
          }
 
629
          
 
630
        }// end loop over 's'
 
631
        for (unsigned int s = 0; s < 3; s++)
 
632
        {
 
633
          for (unsigned int t = 0; t < 3; t++)
 
634
          {
 
635
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
636
          }// end loop over 't'
 
637
        }// end loop over 's'
 
638
      }// end loop over 'r'
 
639
      
 
640
      // Transform derivatives back to physical element
 
641
      for (unsigned int r = 0; r < num_derivatives; r++)
 
642
      {
 
643
        for (unsigned int s = 0; s < num_derivatives; s++)
 
644
        {
 
645
          values[r] += transform[r][s]*derivatives[s];
 
646
        }// end loop over 's'
 
647
      }// end loop over 'r'
 
648
      
 
649
      // Delete pointer to array of derivatives on FIAT element
 
650
      delete [] derivatives;
 
651
      
 
652
      // Delete pointer to array of combinations of derivatives and transform
 
653
      for (unsigned int r = 0; r < num_derivatives; r++)
 
654
      {
 
655
        delete [] combinations[r];
 
656
      }// end loop over 'r'
 
657
      delete [] combinations;
 
658
      for (unsigned int r = 0; r < num_derivatives; r++)
 
659
      {
 
660
        delete [] transform[r];
 
661
      }// end loop over 'r'
 
662
      delete [] transform;
 
663
        break;
 
664
      }
 
665
    case 2:
 
666
      {
 
667
        
 
668
      // Array of basisvalues.
 
669
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
670
      
 
671
      // Declare helper variables.
 
672
      unsigned int rr = 0;
 
673
      unsigned int ss = 0;
 
674
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
675
      
 
676
      // Compute basisvalues.
 
677
      basisvalues[0] = 1.00000000;
 
678
      basisvalues[1] = tmp0;
 
679
      for (unsigned int r = 0; r < 1; r++)
 
680
      {
 
681
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
682
        ss = r*(r + 1)/2;
 
683
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
684
      }// end loop over 'r'
 
685
      for (unsigned int r = 0; r < 2; r++)
 
686
      {
 
687
        for (unsigned int s = 0; s < 2 - r; s++)
 
688
        {
 
689
          rr = (r + s)*(r + s + 1)/2 + s;
 
690
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
691
        }// end loop over 's'
 
692
      }// end loop over 'r'
 
693
      
 
694
      // Table(s) of coefficients.
 
695
      static const double coefficients0[3] = \
 
696
      {0.47140452, 0.00000000, 0.33333333};
 
697
      
 
698
      // Tables of derivatives of the polynomial base (transpose).
 
699
      static const double dmats0[3][3] = \
 
700
      {{0.00000000, 0.00000000, 0.00000000},
 
701
      {4.89897949, 0.00000000, 0.00000000},
 
702
      {0.00000000, 0.00000000, 0.00000000}};
 
703
      
 
704
      static const double dmats1[3][3] = \
 
705
      {{0.00000000, 0.00000000, 0.00000000},
 
706
      {2.44948974, 0.00000000, 0.00000000},
 
707
      {4.24264069, 0.00000000, 0.00000000}};
 
708
      
 
709
      // Compute reference derivatives.
 
710
      // Declare pointer to array of derivatives on FIAT element.
 
711
      double *derivatives = new double[num_derivatives];
 
712
      for (unsigned int r = 0; r < num_derivatives; r++)
 
713
      {
 
714
        derivatives[r] = 0.00000000;
 
715
      }// end loop over 'r'
 
716
      
 
717
      // Declare derivative matrix (of polynomial basis).
 
718
      double dmats[3][3] = \
 
719
      {{1.00000000, 0.00000000, 0.00000000},
 
720
      {0.00000000, 1.00000000, 0.00000000},
 
721
      {0.00000000, 0.00000000, 1.00000000}};
 
722
      
 
723
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
724
      double dmats_old[3][3] = \
 
725
      {{1.00000000, 0.00000000, 0.00000000},
 
726
      {0.00000000, 1.00000000, 0.00000000},
 
727
      {0.00000000, 0.00000000, 1.00000000}};
 
728
      
 
729
      // Loop possible derivatives.
 
730
      for (unsigned int r = 0; r < num_derivatives; r++)
 
731
      {
 
732
        // Resetting dmats values to compute next derivative.
 
733
        for (unsigned int t = 0; t < 3; t++)
 
734
        {
 
735
          for (unsigned int u = 0; u < 3; u++)
 
736
          {
 
737
            dmats[t][u] = 0.00000000;
 
738
            if (t == u)
 
739
            {
 
740
            dmats[t][u] = 1.00000000;
 
741
            }
 
742
            
 
743
          }// end loop over 'u'
 
744
        }// end loop over 't'
 
745
        
 
746
        // Looping derivative order to generate dmats.
 
747
        for (unsigned int s = 0; s < n; s++)
 
748
        {
 
749
          // Updating dmats_old with new values and resetting dmats.
 
750
          for (unsigned int t = 0; t < 3; t++)
 
751
          {
 
752
            for (unsigned int u = 0; u < 3; u++)
 
753
            {
 
754
              dmats_old[t][u] = dmats[t][u];
 
755
              dmats[t][u] = 0.00000000;
 
756
            }// end loop over 'u'
 
757
          }// end loop over 't'
 
758
          
 
759
          // Update dmats using an inner product.
 
760
          if (combinations[r][s] == 0)
 
761
          {
 
762
          for (unsigned int t = 0; t < 3; t++)
 
763
          {
 
764
            for (unsigned int u = 0; u < 3; u++)
 
765
            {
 
766
              for (unsigned int tu = 0; tu < 3; tu++)
 
767
              {
 
768
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
769
              }// end loop over 'tu'
 
770
            }// end loop over 'u'
 
771
          }// end loop over 't'
 
772
          }
 
773
          
 
774
          if (combinations[r][s] == 1)
 
775
          {
 
776
          for (unsigned int t = 0; t < 3; t++)
 
777
          {
 
778
            for (unsigned int u = 0; u < 3; u++)
 
779
            {
 
780
              for (unsigned int tu = 0; tu < 3; tu++)
 
781
              {
 
782
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
783
              }// end loop over 'tu'
 
784
            }// end loop over 'u'
 
785
          }// end loop over 't'
 
786
          }
 
787
          
 
788
        }// end loop over 's'
 
789
        for (unsigned int s = 0; s < 3; s++)
 
790
        {
 
791
          for (unsigned int t = 0; t < 3; t++)
 
792
          {
 
793
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
794
          }// end loop over 't'
 
795
        }// end loop over 's'
 
796
      }// end loop over 'r'
 
797
      
 
798
      // Transform derivatives back to physical element
 
799
      for (unsigned int r = 0; r < num_derivatives; r++)
 
800
      {
 
801
        for (unsigned int s = 0; s < num_derivatives; s++)
 
802
        {
 
803
          values[r] += transform[r][s]*derivatives[s];
 
804
        }// end loop over 's'
 
805
      }// end loop over 'r'
 
806
      
 
807
      // Delete pointer to array of derivatives on FIAT element
 
808
      delete [] derivatives;
 
809
      
 
810
      // Delete pointer to array of combinations of derivatives and transform
 
811
      for (unsigned int r = 0; r < num_derivatives; r++)
 
812
      {
 
813
        delete [] combinations[r];
 
814
      }// end loop over 'r'
 
815
      delete [] combinations;
 
816
      for (unsigned int r = 0; r < num_derivatives; r++)
 
817
      {
 
818
        delete [] transform[r];
 
819
      }// end loop over 'r'
 
820
      delete [] transform;
 
821
        break;
 
822
      }
 
823
    }
 
824
    
422
825
  }
423
826
 
424
827
  /// Evaluate order n derivatives of all basis functions at given point in cell