Lagrange.h 25.2 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
// ============================================================================
// ==                                                                        ==
// == 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 <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
65
66

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


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

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

    /** \brief
     * Implements BasisFunction::getDOFIndices
     */
89
90
    const DegreeOfFreedom* getDOFIndices(const Element*,
					 const DOFAdmin&, 
91
92
93
94
95
					 DegreeOfFreedom*) const;

    /** \brief
     * Implements BasisFunction::getBound
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
96
    void getBound(const ElInfo*, BoundaryType *) const;
97
98
99
100
101
102
103
104
105
106
107
108
109
110

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

111
112
    /// Implements BasisFunction::refineInter
    inline void refineInter(DOFIndexed<double> *drv, RCNeighbourList* list, int n) 
113
    {
114
      if (refineInter_fct)
115
	(*refineInter_fct)(drv, list, n, this);
116
    }
117

118
119
    /// Implements BasisFunction::coarseRestrict
    inline void coarseRestr(DOFIndexed<double> *drv, RCNeighbourList* list, int n)
120
    {
121
      if (coarseRestr_fct)
122
	(*coarseRestr_fct)(drv, list, n, this);
123
    }
124
  
125
126
    /// Implements BasisFunction::coarseInter
    inline void coarseInter(DOFIndexed<double> *drv, RCNeighbourList* list, int n) 
127
    {
128
      if (coarseInter_fct)
129
	(*coarseInter_fct)(drv, list, n, this);
Thomas Witkowski's avatar
Thomas Witkowski committed
130
    }
131
  
132
    /// Implements BasisFunction::getLocalIndices().
133
134
135
136
    const DegreeOfFreedom *getLocalIndices(const Element*,
					   const DOFAdmin*,
					   DegreeOfFreedom*) const;

137
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
140
141
    void getLocalIndicesVec(const Element*,
			    const DOFAdmin*,
			    Vector<DegreeOfFreedom>*) const;

142
    /// Implements BasisFunction::l2ScpFctBas
143
144
145
146
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<double, WorldVector<double> >* f,
		     DOFVector<double>* fh);

147
    /// Implements BasisFunction::l2ScpFctBas
148
149
150
151
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
		     DOFVector<WorldVector<double> >* fh);

152
153
    static void clear();

154
155
156
157
158
159
160
161
162
163
164
  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, 
165
		      DimVec<double>* vec = NULL);
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184

    /** \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
     */
185
186
    int* orderOfPositionIndices(const Element* el, 
				GeoIndex position, 
187
188
189
190
191
192
193
194
195
196
197
198
				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
     */
199
    std::vector<DimVec<double>* > *bary;
200
201
202
203

    /** \name static dim-degree-arrays
     * \{
     */
204
205
206
207
208
209
    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];
210
211
212
213
214
215
    /** \} */

    /** \brief
     * List of all used BasisFunctions in the whole program. Avoids duplicate
     * instantiation of identical BasisFunctions. 
     */
216
    static std::list<Lagrange*> allBasFcts;
217
218
219
220
221
222
223
224
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
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


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

319
320


321
322
323
324
325
326
327
328
329
330
331
332
333
    // ===== Phi =========================================================

    /** \brief
     * AbstractFunction which implements lagrange basis functions
     */
    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.
       */
334
      Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349

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

    private:
      /** \brief
       * vertices needed for evaluation of this function
       */
      int* vertices;
    
      /** \brief
       * Pointer to the evaluating function
       */
350
      double (*func)(const DimVec<double>& lambda, int* vert);
351
352
353
354

      /** \brief
       * Returns \ref func(lambda, vertices)
       */
355
      inline double operator()(const DimVec<double>& lambda) const {
356
	return func(lambda, vertices);
357
      }
358
359
360
361
362

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

363
      // ====== Lagrange, degree = 0 =====================================
364
      // center
365
366
      inline static double phi0c(const DimVec<double>&, int*) {
	return 1.0;
367
      }
368

369
      // ====== Lagrange, degree = 1 =====================================
370
      // vertex
371
372
      inline static double phi1v(const DimVec<double>& lambda, int* vertices) {
	return lambda[vertices[0]]; 
373
      }
374

375
      // ====== Lagrange, degree = 2 =====================================
376
      // vertex
377
378
      inline static double phi2v(const DimVec<double>& lambda, int* vertices) {
	return lambda[vertices[0]] * (2.0 * lambda[vertices[0]] - 1.0);
379
      }    
380
381

      // edge
382
383
      inline static double phi2e(const DimVec<double>& lambda, int* vertices) {
	return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]);
384
      }
385

386
      // ====== Lagrange, degree = 3 =====================================
387
      // vertex
388
389
      inline static double phi3v(const DimVec<double>& lambda, int* vertices) {
	return (4.5 * (lambda[vertices[0]] - 1.0) * lambda[vertices[0]] + 1.0) * 
390
	  lambda[vertices[0]];
391
      }
392
393

      // edge
394
395
396
      inline static double phi3e(const DimVec<double>& lambda, int* vertices) {
	return (13.5 * lambda[vertices[0]] - 4.5) * 
	  lambda[vertices[0]] * lambda[vertices[1]];
397
      }
398
399

      // face
400
401
402
      inline static double phi3f(const DimVec<double>& lambda, int* vertices) {
	return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]];
403
      }
404

405
      // ====== Lagrange, degree = 4 ======================================
406
      // vertex
407
408
409
      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;
410
      }
411
412

      // edge
413
414
415
      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;
416
      }
417

418
419
420
      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;
421
      }
422
423

      // face
424
425
426
      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;
427
      }
428
429

      // center
430
431
432
      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]];
433
      }
434
435
436
437
438

    };
  
    /** \} */

439
440


441
442
443
444
445
446
447
448
449
    // ===== GrdPhi ======================================================

    /** \brief
     * AbstractFunction which implements gradients of lagrange basis functions.
     * See \ref Phi
     */
    class GrdPhi : public GrdBasFctType
    {
    public:
450
451
      GrdPhi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);

452
453
454
      virtual ~GrdPhi();
    private:
      int* vertices;
455
456
457
      void (*func)(const DimVec<double>& lambda,
		   int* vertices_,
		   DimVec<double>& result);
458

459
460
      inline void operator()(const DimVec<double>& lambda, DimVec<double>& result) const {
	func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
461
      }
462
463
464

      // ====== Lagrange0 ================================================
      // center
465
466
467
468
469
      inline static void grdPhi0c(const DimVec<double>&, 
				  int*, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
470
      }
471
472
473

      // ====== Lagrange1 ================================================
      // vertex
474
475
476
477
478
479
      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
480
      }
481
482
483

      // ====== Lagrange2 ================================================
      // vertex
484
485
486
487
488
489
      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
490
      }
491
492

      // edge
493
494
495
496
497
498
499
      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
500
      }
501
502
503

      // ===== Lagrange3 ================================================
      // vertex
504
505
506
507
508
509
510
      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
511
      }
512
513

      // edge
514
515
516
517
518
519
      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) * 
520
	  lambda[vertices[1]];
521
	result[vertices[1]] = (13.5 * lambda[vertices[0]] - 4.5) * 
522
	  lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
523
      }
524
525

      // face
526
527
528
529
530
531
532
533
      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
534
      }
535
536
537
538

    
      // ===== Lagrange4 ================================================
      // vertex
539
540
541
542
543
544
545
546
      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
547
      }
548
549
    
      // edge
550
551
552
553
554
555
556
557
558
      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
559
      }
560
561
562
563
564
565
566
567
568
569

      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
570
      }
571
572

      // face
573
574
575
576
577
578
579
580
581
582
583
      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
584
      }
585
586
    
      // center
587
588
589
590
591
592
593
594
595
      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
596
      }
597
598
    };

599
600


601
602
603
604
605
606
607
608
609
610
    
    // ===== D2Phi =======================================================

    /** \brief
     * AbstractFunction which implements second derivatives of Lagrange basis
     * functions. See \ref Phi
     */
    class D2Phi : public D2BasFctType
    {
    public:
611
      D2Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
612
613
614
615
616

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

619
620
      inline void operator()(const DimVec<double>& lambda, DimMat<double>& result) const {
	return func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
621
      }
622
623
624

      // ===== Lagrange0 ================================================
      // center
625
626
627
      inline static void D2Phi0c(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);       
Thomas Witkowski's avatar
Thomas Witkowski committed
628
      }
629
630
631

      // ===== Lagrange1 ================================================
      // vertex
632
633
634
      inline static void D2Phi1v(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
635
      }
636
637
638

      // ===== Lagrange2 ================================================
      // vertex
639
640
641
642
643
      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
644
      }
645
646

      // edge
647
648
649
650
651
652
      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
653
      }
654
655
656
657


      // ===== Lagrange3 ================================================
      // vertex
658
659
660
661
662
      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
663
      }
664
665

      // edge
666
667
668
669
670
671
672
      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
673
      }
674
675

      // face
676
677
678
679
680
681
682
683
684
685
      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
686
      }
687
688
689
690


      // ===== Lagrange4 ================================================
      // vertex
691
692
693
694
695
696
      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
697
      }
698
699

      // edge
700
701
702
703
704
705
706
707
708
      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
709
      }
710
711
712
713
714
715
716
717
718
719
720
721

      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
722
      }
723
724

      // face
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
      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
740
      }
741
742

      // center
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
      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
765
      }
766
767
768
769
770
771
    };
  };

}

#endif // AMDIS_LAGRANGE_H