Lagrange.h 24.6 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut für Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
// ==                                                                        ==
// ============================================================================

/** \file Lagrange.h */

#ifndef AMDIS_LAGRANGE_H
#define AMDIS_LAGRANGE_H

#include "BasisFunction.h"
#include "FixVec.h"
#include "AbstractFunction.h"
#include <list>

namespace AMDiS {

#define MAX_DIM 3
#define MAX_DEGREE 4

  template<typename ReturnType, typename ArgumentType> class AbstractFunction;

  /** \ingroup FEMSpace
   * \brief
   * Lagrange basis functions. Sub class of BasisFunction
   */
  class Lagrange : public BasisFunction
  {
  protected:
    /** \brief
     * Constructs lagrange basis functions with the given dim and degree.
     * Constructor is protected to avoid multiple instantiation of identical
     * basis functions. Use \ref getLagrange instead.
     */
    Lagrange(int dim_, int degree_);

    /** \brief
     * destructor
     */
    virtual ~Lagrange();

  public:
    /** \brief
     * Returns a pointer to lagrange basis functions with the given dim and
     * degree. Multiple instantiation of identical basis functions is avoided
     * by rembering once created basis functions in \ref allBasFcts.
     */
62
    static Lagrange* getLagrange(int dim, int degree);
63

64
    /// Implements BasisFunction::interpol
65
66
67
    const double *interpol(const ElInfo *, int, const int *, 
			   AbstractFunction<double, WorldVector<double> >*, 
			   double *);
68

69
    /// Implements BasisFunction::interpol
70
71
72
73
    const WorldVector<double> *interpol(const ElInfo *, int, 
					const int *b_no,
					AbstractFunction<WorldVector<double>, 
					WorldVector<double> >*, 
74
75
					WorldVector<double> *);

76
    /// Returns the barycentric coordinates of the i-th basis function.
77
78
    DimVec<double> *getCoords(int i) const;

79
    /// Implements BasisFunction::getDOFIndices
80
81
    const DegreeOfFreedom* getDOFIndices(const Element*,
					 const DOFAdmin&, 
82
83
					 DegreeOfFreedom*) const;

84
    /// Implements BasisFunction::getBound
Thomas Witkowski's avatar
Thomas Witkowski committed
85
    void getBound(const ElInfo*, BoundaryType *) const;
86
87
88
89
90
91
92
93
94
95
96
97
98
99

    /** \brief
     * Calculates the local vertex indices which are involved in evaluating
     * the nodeIndex-th DOF at the positionIndex-th part of type position 
     * (VERTEX/EDGE/FACE/CENTER). nodeIndex determines the permutation
     * of the involved vertices. So in 1d for lagrange4 there are two DOFs at
     * the CENTER (which is an edge in this case). Then vertices[0] = {0, 1} and
     * vertices[1] = {1, 0}. This allows to use the same local basis function 
     * for all DOFs at the same position.
     */
    static void setVertices(int dim, int degree, 
			    GeoIndex position, int positionIndex, int nodeIndex, 
			    int** vertices);

100
101
    /// Implements BasisFunction::refineInter
    inline void refineInter(DOFIndexed<double> *drv, RCNeighbourList* list, int n) 
102
    {
103
      if (refineInter_fct)
104
	(*refineInter_fct)(drv, list, n, this);
105
    }
106

107
108
    /// Implements BasisFunction::coarseRestrict
    inline void coarseRestr(DOFIndexed<double> *drv, RCNeighbourList* list, int n)
109
    {
110
      if (coarseRestr_fct)
111
	(*coarseRestr_fct)(drv, list, n, this);
112
    }
113
  
114
115
    /// Implements BasisFunction::coarseInter
    inline void coarseInter(DOFIndexed<double> *drv, RCNeighbourList* list, int n) 
116
    {
117
      if (coarseInter_fct)
118
	(*coarseInter_fct)(drv, list, n, this);
Thomas Witkowski's avatar
Thomas Witkowski committed
119
    }
120
  
121
    /// Implements BasisFunction::getLocalIndices().
Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
124
    const DegreeOfFreedom *getLocalIndices(const Element *el,
					   const DOFAdmin *admin,
					   DegreeOfFreedom *dofs) const;
125

Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
128
    void getLocalDofPtrVec(const Element *el, 
			   const DOFAdmin *admin,
			   std::vector<const DegreeOfFreedom*>& vec) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
129

130
    /// Implements BasisFunction::l2ScpFctBas
131
132
133
134
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<double, WorldVector<double> >* f,
		     DOFVector<double>* fh);

135
    /// Implements BasisFunction::l2ScpFctBas
136
137
138
139
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
		     DOFVector<WorldVector<double> >* fh);

140
141
    static void clear();

142
143
144
145
146
147
148
  protected:
    /** \brief
     * sets the barycentric coordinates (stored in \ref bary) of the local 
     * basis functions.
     */
    void setBary();

149
    /// Recursive calculation of coordinates. Used by \ref setBary
150
    void createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
151
		      DimVec<double>* vec = NULL);
152

153
    /// Used by \ref setBary
154
155
    int** getIndexPermutations(int numIndices) const;

156
    /// Implements BasisFunction::setNDOF
157
158
    void setNDOF();

159
    /// Sets used function pointers
160
161
    void setFunctionPointer();

162
    /// Used by \ref getDOFIndices and \ref getVec
163
164
    int* orderOfPositionIndices(const Element* el, 
				GeoIndex position, 
165
166
				int positionIndex) const;

167
    /// Calculates the number of DOFs needed for Lagrange of the given dim and degree.
168
    static int getNumberOfDofs(int dim, int degree);
169
170

  private:
171
    /// barycentric coordinates of the locations of all basis functions
172
    std::vector<DimVec<double>* > *bary;
173
174
175
176

    /** \name static dim-degree-arrays
     * \{
     */
177
178
179
180
181
182
    static std::vector<DimVec<double>* > baryDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
    static DimVec<int>* ndofDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
    static int nBasFctsDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
    static std::vector<BasFctType*> phiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
    static std::vector<GrdBasFctType*> grdPhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
    static std::vector<D2BasFctType*> D2PhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
183
184
185
186
187
188
    /** \} */

    /** \brief
     * List of all used BasisFunctions in the whole program. Avoids duplicate
     * instantiation of identical BasisFunctions. 
     */
189
    static std::list<Lagrange*> allBasFcts;
190
191
192


  protected:
193
194
    /// Pointer to the used refineInter function
    void (*refineInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222

    /** \name refineInter functions
     * \{
     */  
    static void  refineInter0(DOFIndexed<double> *, RCNeighbourList*, int, 
			      BasisFunction*);
    static void  refineInter1(DOFIndexed<double> *, RCNeighbourList*, int, 
			      BasisFunction*);
    static void  refineInter2_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter2_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter2_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter3_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter3_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter3_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter4_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter4_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  refineInter4_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    /** \} */

223
224
    /// Pointer to the used coarseRestr function
    void (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252

    /** \name coarseRestr functions
     * \{
     */  
    static void  coarseRestr0(DOFIndexed<double> *, RCNeighbourList*, int, 
			      BasisFunction*);
    static void  coarseRestr1(DOFIndexed<double> *, RCNeighbourList*, int, 
			      BasisFunction*);
    static void  coarseRestr2_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr2_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr2_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr3_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr3_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr3_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr4_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr4_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseRestr4_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    /** \} */

253
254
    /// Pointer to the used coarseInter function
    void (*coarseInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);  
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282

    /** \name coarseInter functions
     * \{
     */
    static void  coarseInter0(DOFIndexed<double> *, RCNeighbourList*, int, 
			      BasisFunction*);
    static void  coarseInter1(DOFIndexed<double> *, RCNeighbourList*, int, 
			      BasisFunction*);
    static void  coarseInter2_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter2_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter2_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter3_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter3_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter3_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter4_1d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter4_2d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    static void  coarseInter4_3d(DOFIndexed<double> *, RCNeighbourList*, int, 
				 BasisFunction*);
    /** \} */

283

284
    /// AbstractFunction which implements lagrange basis functions
285
286
287
288
289
290
291
292
    class Phi : public BasFctType
    {
    public:
      /** \brief
       * Constructs the local lagrange basis function for the given position,
       * positionIndex and nodeIndex. owner_ is a pointer to the Lagrange object
       * this basis function belongs to.
       */
293
      Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
294

295
      /// Destructor
296
297
298
      virtual ~Phi();

    private:
299
      /// vertices needed for evaluation of this function
300
301
      int* vertices;
    
302
      /// Pointer to the evaluating function
303
      double (*func)(const DimVec<double>& lambda, int* vert);
304

305
306
307
      /// Returns \ref func(lambda, vertices)
      inline double operator()(const DimVec<double>& lambda) const 
      {
308
	return func(lambda, vertices);
309
      }
310
311
312
313
314

      /** \name basis functions for different degrees
       * \{
       */

315
      // ====== Lagrange, degree = 0 =====================================
316
      // center
317
318
      inline static double phi0c(const DimVec<double>&, int*) 
      {
319
	return 1.0;
320
      }
321

322
      // ====== Lagrange, degree = 1 =====================================
323
      // vertex
324
325
      inline static double phi1v(const DimVec<double>& lambda, int* vertices) 
      {
326
	return lambda[vertices[0]]; 
327
      }
328

329
      // ====== Lagrange, degree = 2 =====================================
330
      // vertex
331
332
      inline static double phi2v(const DimVec<double>& lambda, int* vertices) 
      {
333
	return lambda[vertices[0]] * (2.0 * lambda[vertices[0]] - 1.0);
334
      }    
335
336

      // edge
337
338
      inline static double phi2e(const DimVec<double>& lambda, int* vertices) 
      {
339
	return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]);
340
      }
341

342
      // ====== Lagrange, degree = 3 =====================================
343
      // vertex
344
345
      inline static double phi3v(const DimVec<double>& lambda, int* vertices) 
      {
346
	return (4.5 * (lambda[vertices[0]] - 1.0) * lambda[vertices[0]] + 1.0) * 
347
	  lambda[vertices[0]];
348
      }
349
350

      // edge
351
352
      inline static double phi3e(const DimVec<double>& lambda, int* vertices) 
      {
353
354
	return (13.5 * lambda[vertices[0]] - 4.5) * 
	  lambda[vertices[0]] * lambda[vertices[1]];
355
      }
356
357

      // face
358
359
      inline static double phi3f(const DimVec<double>& lambda, int* vertices) 
      {
360
361
	return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]];
362
      }
363

364
      // ====== Lagrange, degree = 4 ======================================
365
      // vertex
366
367
      inline static double phi4v(const DimVec<double>& lambda, int* vertices) 
      {
368
369
	return (((32.0 * lambda[vertices[0]] - 48.0) * lambda[vertices[0]] + 22.0)
	       * lambda[vertices[0]] - 3.0) * lambda[vertices[0]] / 3.0;
370
      }
371
372

      // edge
373
374
      inline static double phi4e0(const DimVec<double>& lambda, int* vertices) 
      {
375
376
	return ((128.0 * lambda[vertices[0]] - 96.0) * lambda[vertices[0]] + 16.0)
	  * lambda[vertices[0]] * lambda[vertices[1]] / 3.0;
377
      }
378

379
380
      inline static double phi4e1(const DimVec<double>& lambda, int* vertices) 
      {
381
382
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  (4.0 * lambda[vertices[1]] - 1.0) * lambda[vertices[1]] * 4.0;
383
      }
384
385

      // face
386
387
      inline static double phi4f(const DimVec<double>& lambda,  int* vertices) 
      {
388
389
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  lambda[vertices[1]] * lambda[vertices[2]] * 32.0;
390
      }
391
392

      // center
393
394
      inline static double phi4c(const DimVec<double>& lambda, int* vertices) 
      {
395
396
	return 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]] * lambda[vertices[3]];
397
      }
398
399
400
401
402

    };
  
    /** \} */

403
404


405
406
407
408
409
410
411
    /** \brief
     * AbstractFunction which implements gradients of lagrange basis functions.
     * See \ref Phi
     */
    class GrdPhi : public GrdBasFctType
    {
    public:
412
413
      GrdPhi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);

414
415
416
      virtual ~GrdPhi();
    private:
      int* vertices;
417
418

      void (*func)(const DimVec<double>& lambda, int* vertices_, 
419
		   DimVec<double>& result);
420

421
422
      inline void operator()(const DimVec<double>& lambda, DimVec<double>& result) const 
      {
423
	func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
424
      }
425
426
427

      // ====== Lagrange0 ================================================
      // center
428
429
430
431
432
      inline static void grdPhi0c(const DimVec<double>&, 
				  int*, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
433
      }
434
435
436

      // ====== Lagrange1 ================================================
      // vertex
437
438
439
440
441
442
      inline static void grdPhi1v(const DimVec<double>&, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
443
      }
444
445
446

      // ====== Lagrange2 ================================================
      // vertex
447
448
449
450
451
452
      inline static void grdPhi2v(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 4.0 * lambda[vertices[0]] - 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
453
      }
454
455

      // edge
456
457
458
459
460
461
462
      inline static void grdPhi2e(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 4.0 * lambda[vertices[1]];
	result[vertices[1]] = 4.0 * lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
463
      }
464
465
466

      // ===== Lagrange3 ================================================
      // vertex
467
468
469
470
471
472
473
      inline static void grdPhi3v(const DimVec<double>& lambda,
				  int* vertices,
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = (13.5 * lambda[vertices[0]] - 9.0) * 
	  lambda[vertices[0]] + 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
474
      }
475
476

      // edge
477
478
479
480
481
482
      inline static void grdPhi3e(const DimVec<double>& lambda,
				  int* vertices, 
				  DimVec<double>& result)
      {
	result.set(0.0);
	result[vertices[0]] = (27.0 * lambda[vertices[0]] - 4.5) * 
483
	  lambda[vertices[1]];
484
	result[vertices[1]] = (13.5 * lambda[vertices[0]] - 4.5) * 
485
	  lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
486
      }
487
488

      // face
489
490
491
492
493
494
495
496
      inline static void grdPhi3f(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 27.0 * lambda[vertices[1]] * lambda[vertices[2]];
	result[vertices[1]] = 27.0 * lambda[vertices[0]] * lambda[vertices[2]];
	result[vertices[2]] = 27.0 * lambda[vertices[0]] * lambda[vertices[1]];
Thomas Witkowski's avatar
Thomas Witkowski committed
497
      }
498
499
500
501

    
      // ===== Lagrange4 ================================================
      // vertex
502
503
504
505
506
507
508
509
      inline static void grdPhi4v(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 
	  ((128.0 * lambda[vertices[0]] - 144.0) * lambda[vertices[0]] + 44.0) * 
	  lambda[vertices[0]] / 3.0 - 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
510
      }
511
512
    
      // edge
513
514
515
516
517
518
519
520
521
      inline static void grdPhi4e0(const DimVec<double>& lambda,
				   int* vertices, 
				   DimVec<double>& result)
      {
	result.set(0.0);
	result[vertices[0]] = ((128.0 * lambda[vertices[0]] - 64.0) * 
			       lambda[vertices[0]] + 16.0 / 3.0) * lambda[vertices[1]];
	result[vertices[1]] = ((128.0 * lambda[vertices[0]] - 96.0) * 
			       lambda[vertices[0]] + 16.0)*lambda[vertices[0]] / 3.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
522
      }
523
524
525
526
527
528
529
530
531
532

      inline static void grdPhi4e1(const DimVec<double>& lambda,
				   int* vertices,
				   DimVec<double>& result)
      {
	result.set(0.0);
	result[vertices[0]] = 4.0 * (8.0 * lambda[vertices[0]] - 1.0) * 
	  lambda[vertices[1]] * (4.0 * lambda[vertices[1]] - 1.0);
	result[vertices[1]] = 4.0 * lambda[vertices[0]] * 
	  (4.0 * lambda[vertices[0]] - 1.0) * (8.0 * lambda[vertices[1]] - 1.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
533
      }
534
535

      // face
536
537
538
539
540
541
542
543
544
545
546
      inline static void grdPhi4f(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result)
      {
	result.set(0.0);
	result[vertices[0]] = 32.0 * (8.0 * lambda[vertices[0]] - 1.0) * 
	  lambda[vertices[1]] * lambda[vertices[2]];
	result[vertices[1]] = 32.0 * (4.0 * lambda[vertices[0]] - 1.0) * 
	  lambda[vertices[0]] * lambda[vertices[2]];
	result[vertices[2]] = 32.0 * (4.0 * lambda[vertices[0]] - 1.0) * 
	  lambda[vertices[0]] * lambda[vertices[1]];	
Thomas Witkowski's avatar
Thomas Witkowski committed
547
      }
548
549
    
      // center
550
551
552
553
554
555
556
557
558
      inline static void grdPhi4c(const DimVec<double>& lambda, 
				  int* vertices,
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[0] = 256.0 * lambda[vertices[1]] * lambda[vertices[2]] * lambda[vertices[3]];
	result[1] = 256.0 * lambda[vertices[0]] * lambda[vertices[2]] * lambda[vertices[3]];
	result[2] = 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[3]];
	result[3] = 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[2]];
Thomas Witkowski's avatar
Thomas Witkowski committed
559
      }
560
561
    };

562
563


564
565
566
567
568
569
570
    /** \brief
     * AbstractFunction which implements second derivatives of Lagrange basis
     * functions. See \ref Phi
     */
    class D2Phi : public D2BasFctType
    {
    public:
571
      D2Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
572
573
574
575
576

      virtual ~D2Phi();
    private:
      int* vertices;
    
577
      void (*func)(const DimVec<double>& lambda, int* vertices_, DimMat<double>& result);
578

579
580
      inline void operator()(const DimVec<double>& lambda, DimMat<double>& result) const {
	return func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
581
      }
582
583
584

      // ===== Lagrange0 ================================================
      // center
585
586
587
      inline static void D2Phi0c(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);       
Thomas Witkowski's avatar
Thomas Witkowski committed
588
      }
589
590
591

      // ===== Lagrange1 ================================================
      // vertex
592
593
594
      inline static void D2Phi1v(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
595
      }
596
597
598

      // ===== Lagrange2 ================================================
      // vertex
599
600
601
602
603
      inline static void D2Phi2v(const DimVec<double>&, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 4.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
604
      }
605
606

      // edge
607
608
609
610
611
612
      inline static void D2Phi2e(const DimVec<double>&, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[1]] = 4.0;
	result[vertices[1]][vertices[0]] = 4.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
613
      }
614
615
616
617


      // ===== Lagrange3 ================================================
      // vertex
618
619
620
621
622
      inline static void D2Phi3v(const DimVec<double>& lambda, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 27.0 * lambda[vertices[0]] - 9.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
623
      }
624
625

      // edge
626
627
628
629
630
631
632
      inline static void D2Phi3e(const DimVec<double>& lambda, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 27.0 * lambda[vertices[1]];
	result[vertices[0]][vertices[1]] = 
	  result[vertices[1]][vertices[0]] = 27.0 * lambda[vertices[0]] - 4.5;
Thomas Witkowski's avatar
Thomas Witkowski committed
633
      }
634
635

      // face
636
637
638
639
640
641
642
643
644
645
      inline static void D2Phi3f(const DimVec<double>& lambda, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[1]] = 
	  result[vertices[1]][vertices[0]] = 27.0 * lambda[vertices[2]];
	result[vertices[0]][vertices[2]] = 
	  result[vertices[2]][vertices[0]] = 27.0 * lambda[vertices[1]];
	result[vertices[1]][vertices[2]] = 
	  result[vertices[2]][vertices[1]] = 27.0 * lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
646
      }
647
648
649
650


      // ===== Lagrange4 ================================================
      // vertex
651
652
653
654
655
656
      inline static void D2Phi4v(const DimVec<double>& lambda, 
				 int* vertices, 
				 DimMat<double>& result) {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 
	  (128.0 * lambda[vertices[0]] - 96.0) * lambda[vertices[0]] + 44.0 / 3.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
657
      }
658
659

      // edge
660
661
662
663
664
665
666
667
668
      inline static void D2Phi4e0(const DimVec<double>& lambda, int* vertices, 
				  DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 
	  (256.0 * lambda[vertices[0]] - 64.0) * lambda[vertices[1]];
	result[vertices[0]][vertices[1]] = 
	  result[vertices[1]][vertices[0]] =
	  (128.0 * lambda[vertices[0]] - 64.0) * lambda[vertices[0]] + 16.0 / 3.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
669
      }
670
671
672
673
674
675
676
677
678
679
680
681

      inline static void D2Phi4e1(const DimVec<double>& lambda, int* vertices, 
				  DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 
	  32.0 * lambda[vertices[1]] * (4.0 * lambda[vertices[1]] - 1.0);
	result[vertices[0]][vertices[1]] = 
	  result[vertices[1]][vertices[0]] = 
	  4.0 * (8.0 * lambda[vertices[0]] - 1.0) * (8.0 * lambda[vertices[1]] - 1.0);
	result[vertices[1]][vertices[1]] = 
	  32.0 * lambda[vertices[0]] * (4.0 * lambda[vertices[0]] - 1.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
682
      }
683
684

      // face
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
      inline static void D2Phi4f(const DimVec<double>& lambda, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 
	  256.0 * lambda[vertices[1]] * lambda[vertices[2]];
	result[vertices[0]][vertices[1]] = 
	  result[vertices[1]][vertices[0]] = 
	  32.0 * (8.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[2]];
	result[vertices[0]][vertices[2]] = 
	  result[vertices[2]][vertices[0]] = 
	  32.0 * (8.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[1]];
	result[vertices[1]][vertices[2]] = 
	  result[vertices[2]][vertices[1]] = 
	  32.0 * (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
700
      }
701
702

      // center
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
      inline static void D2Phi4c(const DimVec<double>& lambda, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[1]] = 
	  result[vertices[1]][vertices[0]] = 
	  256.0 * lambda[vertices[2]] * lambda[vertices[3]];
	result[vertices[0]][vertices[2]] = 
	  result[vertices[2]][vertices[0]] = 
	  256.0 * lambda[vertices[1]] * lambda[vertices[3]];
	result[vertices[0]][vertices[3]] = 
	  result[vertices[3]][vertices[0]] = 
	  256.0 * lambda[vertices[1]] * lambda[vertices[2]];
	result[vertices[1]][vertices[2]] = 
	  result[vertices[2]][vertices[1]] = 
	  256.0 * lambda[vertices[0]] * lambda[vertices[3]];
	result[vertices[1]][vertices[3]] = 
	  result[vertices[3]][vertices[1]] = 
	  256.0 * lambda[vertices[0]] * lambda[vertices[2]];
	result[vertices[2]][vertices[3]] = 
	  result[vertices[3]][vertices[2]] = 
	  256.0 * lambda[vertices[0]] * lambda[vertices[1]];
Thomas Witkowski's avatar
Thomas Witkowski committed
725
      }
726
727
728
729
730
731
    };
  };

}

#endif // AMDIS_LAGRANGE_H