Lagrange.h 24.2 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25

/** \file Lagrange.h */

#ifndef AMDIS_LAGRANGE_H
#define AMDIS_LAGRANGE_H

26
27
28
#include <list>
#include <boost/numeric/mtl/mtl.hpp>
#include "AbstractFunction.h"
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
#include "BasisFunction.h"
#include "FixVec.h"

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.
     */
64
    static Lagrange* getLagrange(int dim, int degree);
65

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

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

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

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

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

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

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

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

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

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

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

142
143
    static void clear();

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

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

155
    /// Used by \ref setBary
156
157
    int** getIndexPermutations(int numIndices) const;

158
    /// Implements BasisFunction::setNDOF
159
160
    void setNDOF();

161
    /// Sets used function pointers
162
163
    void setFunctionPointer();

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

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

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

    /** \name static dim-degree-arrays
     * \{
     */
179
180
181
182
183
184
    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];
185
186
187
188
189
190
    /** \} */

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


  protected:
195
196
    /// Pointer to the used refineInter function
    void (*refineInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);
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
223
224

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

225
226
    /// Pointer to the used coarseRestr function
    void (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);
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

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

255
256
    /// Pointer to the used coarseInter function
    void (*coarseInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);  
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

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

285

286
    /// AbstractFunction which implements lagrange basis functions
287
288
289
290
291
292
293
294
    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.
       */
295
      Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
296

297
      /// Destructor
298
299
300
      virtual ~Phi();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    };
  
    /** \} */

405
406


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

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

420
421
422
      void (*func)(const DimVec<double>& lambda, 
		   int* vertices_, 
		   mtl::dense_vector<double>& result);
423

424
425
      inline void operator()(const DimVec<double>& lambda, 
			     mtl::dense_vector<double>& result) const 
426
      {
427
	func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
428
      }
429
430
431

      // ====== Lagrange0 ================================================
      // center
432
433
      inline static void grdPhi0c(const DimVec<double>&, 
				  int*, 
434
				  mtl::dense_vector<double>& result) 
435
      {
436
	result = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
437
      }
438
439
440

      // ====== Lagrange1 ================================================
      // vertex
441
442
      inline static void grdPhi1v(const DimVec<double>&, 
				  int* vertices, 
443
				  mtl::dense_vector<double>& result) 
444
      {
445
	result = 0.0;
446
	result[vertices[0]] = 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
447
      }
448
449
450

      // ====== Lagrange2 ================================================
      // vertex
451
452
      inline static void grdPhi2v(const DimVec<double>& lambda, 
				  int* vertices, 
453
				  mtl::dense_vector<double>& result) 
454
      {
455
	result = 0.0;
456
	result[vertices[0]] = 4.0 * lambda[vertices[0]] - 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
457
      }
458
459

      // edge
460
461
      inline static void grdPhi2e(const DimVec<double>& lambda, 
				  int* vertices, 
462
				  mtl::dense_vector<double>& result) 
463
      {
464
	result = 0.0;
465
466
	result[vertices[0]] = 4.0 * lambda[vertices[1]];
	result[vertices[1]] = 4.0 * lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
467
      }
468
469
470

      // ===== Lagrange3 ================================================
      // vertex
471
472
      inline static void grdPhi3v(const DimVec<double>& lambda,
				  int* vertices,
473
				  mtl::dense_vector<double>& result) 
474
      {
475
	result = 0.0;
476
477
	result[vertices[0]] = (13.5 * lambda[vertices[0]] - 9.0) * 
	  lambda[vertices[0]] + 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
478
      }
479
480

      // edge
481
482
      inline static void grdPhi3e(const DimVec<double>& lambda,
				  int* vertices, 
483
				  mtl::dense_vector<double>& result)
484
      {
485
	result = 0.0;
486
	result[vertices[0]] = (27.0 * lambda[vertices[0]] - 4.5) * 
487
	  lambda[vertices[1]];
488
	result[vertices[1]] = (13.5 * lambda[vertices[0]] - 4.5) * 
489
	  lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
490
      }
491
492

      // face
493
494
      inline static void grdPhi3f(const DimVec<double>& lambda, 
				  int* vertices, 
495
				  mtl::dense_vector<double>& result) 
496
      {
497
	result = 0.0;
498
499
500
	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
501
      }
502
503
504
505

    
      // ===== Lagrange4 ================================================
      // vertex
506
507
      inline static void grdPhi4v(const DimVec<double>& lambda, 
				  int* vertices, 
508
				  mtl::dense_vector<double>& result) 
509
      {
510
	result = 0.0;
511
512
513
	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
514
      }
515
516
    
      // edge
517
518
      inline static void grdPhi4e0(const DimVec<double>& lambda,
				   int* vertices, 
519
				   mtl::dense_vector<double>& result)
520
      {
521
	result = 0.0;
522
523
524
525
	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
526
      }
527
528
529

      inline static void grdPhi4e1(const DimVec<double>& lambda,
				   int* vertices,
530
				   mtl::dense_vector<double>& result)
531
      {
532
	result = 0.0;
533
534
535
536
	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
537
      }
538
539

      // face
540
541
      inline static void grdPhi4f(const DimVec<double>& lambda, 
				  int* vertices, 
542
				  mtl::dense_vector<double>& result)
543
      {
544
	result = 0.0;
545
546
547
548
549
550
	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
551
      }
552
553
    
      // center
554
555
      inline static void grdPhi4c(const DimVec<double>& lambda, 
				  int* vertices,
556
				  mtl::dense_vector<double>& result) 
557
      {
558
559
560
561
562
563
564
565
566
	result = 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
567
      }
568
569
    };

570
571


572
573
574
575
576
577
578
    /** \brief
     * AbstractFunction which implements second derivatives of Lagrange basis
     * functions. See \ref Phi
     */
    class D2Phi : public D2BasFctType
    {
    public:
579
      D2Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
580
581
582
583
584

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

587
588
      inline void operator()(const DimVec<double>& lambda, DimMat<double>& result) const {
	return func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
589
      }
590
591
592

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

      // ===== Lagrange1 ================================================
      // vertex
600
601
602
      inline static void D2Phi1v(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
603
      }
604
605
606

      // ===== Lagrange2 ================================================
      // vertex
607
608
609
610
611
      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
612
      }
613
614

      // edge
615
616
617
618
619
620
      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
621
      }
622
623
624
625


      // ===== Lagrange3 ================================================
      // vertex
626
627
628
629
630
      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
631
      }
632
633

      // edge
634
635
636
637
638
639
640
      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
641
      }
642
643

      // face
644
645
646
647
648
649
650
651
652
653
      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
654
      }
655
656
657
658


      // ===== Lagrange4 ================================================
      // vertex
659
660
661
662
663
664
      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
665
      }
666
667

      // edge
668
669
670
671
672
673
674
675
676
      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
677
      }
678
679
680
681
682
683
684
685
686
687
688
689

      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
690
      }
691
692

      // face
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
      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
708
      }
709
710

      // center
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
      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
733
      }
734
735
736
737
738
739
    };
  };

}

#endif // AMDIS_LAGRANGE_H