Lagrange.h 25.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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
62
63
64
65
66
67
68
69
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file Lagrange.h */

#ifndef AMDIS_LAGRANGE_H
#define AMDIS_LAGRANGE_H

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

namespace AMDiS {

#define MAX_DIM 3
#define MAX_DEGREE 4

  template<typename ReturnType, typename ArgumentType> class AbstractFunction;

  // ============================================================================
  // ===== class Lagrange =======================================================
  // ============================================================================

  /** \ingroup FEMSpace
   * \brief
   * Lagrange basis functions. Sub class of BasisFunction
   */
  class Lagrange : public BasisFunction
  {
  public:
    MEMORY_MANAGED(Lagrange);

  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.
     */
70
    static Lagrange* getLagrange(int dim, int degree);
71
72
73
74

    /** \brief
     * Implements BasisFunction::interpol
     */
75
76
77
    const double *interpol(const ElInfo *, int, const int *, 
			   AbstractFunction<double, WorldVector<double> >*, 
			   double *);
78
79
80
81
82


    /** \brief
     * Implements BasisFunction::interpol
     */
83
84
85
86
    const WorldVector<double> *interpol(const ElInfo *, int, 
					const int *b_no,
					AbstractFunction<WorldVector<double>, 
					WorldVector<double> >*, 
87
88
89
90
91
92
93
94
95
96
					WorldVector<double> *);

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

    /** \brief
     * Implements BasisFunction::getDOFIndices
     */
97
98
    const DegreeOfFreedom* getDOFIndices(const Element*,
					 const DOFAdmin&, 
99
100
101
102
103
					 DegreeOfFreedom*) const;

    /** \brief
     * Implements BasisFunction::getBound
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
104
    void getBound(const ElInfo*, BoundaryType *) const;
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

    /** \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);

    /** \brief
     * Implements BasisFunction::refineInter
     */
    inline void  refineInter(DOFIndexed<double> *drv,
			     RCNeighbourList* list, 
			     int n) 
    {
126
      if (refineInter_fct)
127
	(*refineInter_fct)(drv, list, n, this);
128
    }
129
130
131
132
133
134
135
136

    /** \brief
     * Implements BasisFunction::coarseRestrict
     */
    inline void  coarseRestr(DOFIndexed<double> *drv,
			     RCNeighbourList* list,
			     int n)
    {
137
      if (coarseRestr_fct)
138
	(*coarseRestr_fct)(drv, list, n, this);
139
    }
140
141
142
143
144
145
146
147
  
    /** \brief
     * Implements BasisFunction::coarseInter
     */
    inline void  coarseInter(DOFIndexed<double> *drv, 
			     RCNeighbourList* list, 
			     int n) 
    {
148
      if (coarseInter_fct)
149
150
151
152
153
154
155
156
157
158
	(*coarseInter_fct)(drv, list, n, this);
    };
  
    /** \brief
     * Implements BasisFunction::getLocalIndices().
     */
    const DegreeOfFreedom *getLocalIndices(const Element*,
					   const DOFAdmin*,
					   DegreeOfFreedom*) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
161
162
163

    void getLocalIndicesVec(const Element*,
			    const DOFAdmin*,
			    Vector<DegreeOfFreedom>*) const;

164
165
166
167
168
169
170
171
172
173
174
175
176
177
    /** \brief
     * Implements BasisFunction::l2ScpFctBas
     */
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<double, WorldVector<double> >* f,
		     DOFVector<double>* fh);

    /** \brief
     * Implements BasisFunction::l2ScpFctBas
     */
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
		     DOFVector<WorldVector<double> >* fh);

178
179
    static void clear();

180
181
182
183
184
185
186
187
188
189
190
  protected:
    /** \brief
     * sets the barycentric coordinates (stored in \ref bary) of the local 
     * basis functions.
     */
    void setBary();

    /** \brief
     * Recursive calculation of coordinates. Used by \ref setBary
     */
    void createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
191
		      DimVec<double>* vec = NULL);
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210

    /** \brief
     * Used by \ref setBary
     */
    int** getIndexPermutations(int numIndices) const;

    /** \brief
     * Implements BasisFunction::setNDOF
     */
    void setNDOF();

    /** \brief
     * Sets used function pointers
     */
    void setFunctionPointer();

    /** \brief
     * Used by \ref getDOFIndices and \ref getVec
     */
211
212
    int* orderOfPositionIndices(const Element* el, 
				GeoIndex position, 
213
214
215
216
217
218
219
220
221
222
223
224
				int positionIndex) const;

    /** \brief
     * Calculates the number of DOFs needed for Lagrange of the given dim and
     * degree.
     */
    static int getNumberOfDOFs(int dim, int degree);

  private:
    /** \brief
     * barycentric coordinates of the locations of all basis functions
     */
225
    std::vector<DimVec<double>* > *bary;
226
227
228
229

    /** \name static dim-degree-arrays
     * \{
     */
230
231
232
233
234
235
    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];
236
237
238
239
240
241
    /** \} */

    /** \brief
     * List of all used BasisFunctions in the whole program. Avoids duplicate
     * instantiation of identical BasisFunctions. 
     */
242
    static std::list<Lagrange*> allBasFcts;
243
244
245
246
247
248
249
250
251
252
253
254
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
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344


  protected:
    /** \brief
     * Pointer to the used refineInter function
     */
    void  (*refineInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, 
			     BasisFunction*);

    /** \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*);
    /** \} */

    /** \brief
     * Pointer to the used coarseRestr function
     */
    void  (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, 
			     BasisFunction*);

    /** \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*);
    /** \} */

    /** \brief
     * Pointer to the used coarseInter function
     */
    void  (*coarseInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, 
			     BasisFunction*);  

    /** \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*);
    /** \} */

345
346


347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
    // ===== Phi =========================================================

    /** \brief
     * AbstractFunction which implements lagrange basis functions
     */
    class Phi : public BasFctType
    {
    public:
      MEMORY_MANAGED(Phi);

      /** \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.
       */
362
      Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377

      /** \brief
       * Destructor
       */
      virtual ~Phi();

    private:
      /** \brief
       * vertices needed for evaluation of this function
       */
      int* vertices;
    
      /** \brief
       * Pointer to the evaluating function
       */
378
      double (*func)(const DimVec<double>& lambda, int* vert);
379
380
381
382

      /** \brief
       * Returns \ref func(lambda, vertices)
       */
383
      inline double operator()(const DimVec<double>& lambda) const {
384
	return func(lambda, vertices);
385
      }
386
387
388
389
390

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

391
      // ====== Lagrange, degree = 0 =====================================
392
      // center
393
394
      inline static double phi0c(const DimVec<double>&, int*) {
	return 1.0;
395
      }
396

397
      // ====== Lagrange, degree = 1 =====================================
398
      // vertex
399
400
      inline static double phi1v(const DimVec<double>& lambda, int* vertices) {
	return lambda[vertices[0]]; 
401
      }
402

403
      // ====== Lagrange, degree = 2 =====================================
404
      // vertex
405
406
      inline static double phi2v(const DimVec<double>& lambda, int* vertices) {
	return lambda[vertices[0]] * (2.0 * lambda[vertices[0]] - 1.0);
407
      }    
408
409

      // edge
410
411
      inline static double phi2e(const DimVec<double>& lambda, int* vertices) {
	return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]);
412
      }
413

414
      // ====== Lagrange, degree = 3 =====================================
415
      // vertex
416
417
      inline static double phi3v(const DimVec<double>& lambda, int* vertices) {
	return (4.5 * (lambda[vertices[0]] - 1.0) * lambda[vertices[0]] + 1.0) * 
418
	  lambda[vertices[0]];
419
      }
420
421

      // edge
422
423
424
      inline static double phi3e(const DimVec<double>& lambda, int* vertices) {
	return (13.5 * lambda[vertices[0]] - 4.5) * 
	  lambda[vertices[0]] * lambda[vertices[1]];
425
      }
426
427

      // face
428
429
430
      inline static double phi3f(const DimVec<double>& lambda, int* vertices) {
	return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]];
431
      }
432

433
      // ====== Lagrange, degree = 4 ======================================
434
      // vertex
435
436
437
      inline static double phi4v(const DimVec<double>& lambda, int* vertices) {
	return (((32.0 * lambda[vertices[0]] - 48.0) * lambda[vertices[0]] + 22.0)
	       * lambda[vertices[0]] - 3.0) * lambda[vertices[0]] / 3.0;
438
      }
439
440

      // edge
441
442
443
      inline static double phi4e0(const DimVec<double>& lambda, int* vertices) {
	return ((128.0 * lambda[vertices[0]] - 96.0) * lambda[vertices[0]] + 16.0)
	  * lambda[vertices[0]] * lambda[vertices[1]] / 3.0;
444
      }
445

446
447
448
      inline static double phi4e1(const DimVec<double>& lambda, int* vertices) {
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  (4.0 * lambda[vertices[1]] - 1.0) * lambda[vertices[1]] * 4.0;
449
      }
450
451

      // face
452
453
454
      inline static double phi4f(const DimVec<double>& lambda,  int* vertices) {
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  lambda[vertices[1]] * lambda[vertices[2]] * 32.0;
455
      }
456
457

      // center
458
459
460
      inline static double phi4c(const DimVec<double>& lambda, int* vertices) {
	return 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]] * lambda[vertices[3]];
461
      }
462
463
464
465
466

    };
  
    /** \} */

467
468


469
470
471
472
473
474
475
476
477
478
479
    // ===== GrdPhi ======================================================

    /** \brief
     * AbstractFunction which implements gradients of lagrange basis functions.
     * See \ref Phi
     */
    class GrdPhi : public GrdBasFctType
    {
    public:
      MEMORY_MANAGED(GrdPhi);

480
481
      GrdPhi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);

482
483
484
      virtual ~GrdPhi();
    private:
      int* vertices;
485
486
487
      void (*func)(const DimVec<double>& lambda,
		   int* vertices_,
		   DimVec<double>& result);
488

489
490
      inline void operator()(const DimVec<double>& lambda, DimVec<double>& result) const {
	func(lambda, vertices, result);
491
492
493
494
      };

      // ====== Lagrange0 ================================================
      // center
495
496
497
498
499
      inline static void grdPhi0c(const DimVec<double>&, 
				  int*, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
500
501
502
503
      };

      // ====== Lagrange1 ================================================
      // vertex
504
505
506
507
508
509
      inline static void grdPhi1v(const DimVec<double>&, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 1.0;
510
511
512
513
      };

      // ====== Lagrange2 ================================================
      // vertex
514
515
516
517
518
519
      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;
520
521
522
      };

      // edge
523
524
525
526
527
528
529
      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]];
530
531
532
533
      };

      // ===== Lagrange3 ================================================
      // vertex
534
535
536
537
538
539
540
      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;
541
542
543
      };

      // edge
544
545
546
547
548
549
      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) * 
550
	  lambda[vertices[1]];
551
	result[vertices[1]] = (13.5 * lambda[vertices[0]] - 4.5) * 
552
553
554
555
	  lambda[vertices[0]];
      };

      // face
556
557
558
559
560
561
562
563
      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]];
564
565
566
567
568
      };

    
      // ===== Lagrange4 ================================================
      // vertex
569
570
571
572
573
574
575
576
      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;
577
578
579
      };
    
      // edge
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
      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;
      };

      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);
600
601
602
      };

      // face
603
604
605
606
607
608
609
610
611
612
613
      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]];	
614
615
616
      };
    
      // center
617
618
619
620
621
622
623
624
625
      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]];
626
627
628
      };
    };

629
630


631
632
633
634
635
636
637
638
639
640
641
642
    
    // ===== D2Phi =======================================================

    /** \brief
     * AbstractFunction which implements second derivatives of Lagrange basis
     * functions. See \ref Phi
     */
    class D2Phi : public D2BasFctType
    {
    public:
      MEMORY_MANAGED(D2Phi);

643
      D2Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
644
645
646
647
648

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

651
652
      inline void operator()(const DimVec<double>& lambda, DimMat<double>& result) const {
	return func(lambda, vertices, result);
653
654
655
656
      };

      // ===== Lagrange0 ================================================
      // center
657
658
659
      inline static void D2Phi0c(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);       
660
661
662
663
      };

      // ===== Lagrange1 ================================================
      // vertex
664
665
666
      inline static void D2Phi1v(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);
667
668
669
670
      };

      // ===== Lagrange2 ================================================
      // vertex
671
672
673
674
675
      inline static void D2Phi2v(const DimVec<double>&, int* vertices, 
				 DimMat<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]][vertices[0]] = 4.0;
676
677
678
      };

      // edge
679
680
681
682
683
684
      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;
685
686
687
688
689
      };


      // ===== Lagrange3 ================================================
      // vertex
690
691
692
693
694
      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;
695
696
697
      };

      // edge
698
699
700
701
702
703
704
      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;
705
706
707
      };

      // face
708
709
710
711
712
713
714
715
716
717
      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]];
718
719
720
721
722
      };


      // ===== Lagrange4 ================================================
      // vertex
723
724
725
726
727
728
      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;
729
730
731
      };

      // edge
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
      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;
      };

      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);
754
755
756
      };

      // face
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
      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]];
772
773
774
      };

      // center
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
      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]];
797
798
799
800
801
802
803
      };
    };
  };

}

#endif // AMDIS_LAGRANGE_H