Lagrange.h 23.9 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
#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:
Thomas Witkowski's avatar
Thomas Witkowski committed
46
47
48
    /// 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.
49
50
51
52
53
54
55
56
    Lagrange(int dim_, int degree_);

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

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
57
58
59
    /// 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.
60
    static Lagrange* getLagrange(int dim, int degree);
61

62
    /// Implements BasisFunction::interpol
63
64
65
    void interpol(const ElInfo *, int, const int *, 
		  AbstractFunction<double, WorldVector<double> >*, 
		  mtl::dense_vector<double>&) const;
66

67
    /// Implements BasisFunction::interpol
68
69
70
71
    void interpol(const ElInfo *, int, 
		  const int *b_no,
		  AbstractFunction<WorldVector<double>, WorldVector<double> >*, 
		  mtl::dense_vector<WorldVector<double> >&) const;
72

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

76
    /// Implements BasisFunction::getBound
Thomas Witkowski's avatar
Thomas Witkowski committed
77
    void getBound(const ElInfo*, BoundaryType *) const;
78
79
80
81
82
83
84
85
86
87
88
89
90
91

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

92
93
    /// Implements BasisFunction::refineInter
    inline void refineInter(DOFIndexed<double> *drv, RCNeighbourList* list, int n) 
94
    {
95
      if (refineInter_fct)
96
	(*refineInter_fct)(drv, list, n, this);
97
    }
98

99
100
    /// Implements BasisFunction::coarseRestrict
    inline void coarseRestr(DOFIndexed<double> *drv, RCNeighbourList* list, int n)
101
    {
102
      if (coarseRestr_fct)
103
	(*coarseRestr_fct)(drv, list, n, this);
104
    }
105
  
106
107
    /// Implements BasisFunction::coarseInter
    inline void coarseInter(DOFIndexed<double> *drv, RCNeighbourList* list, int n) 
108
    {
109
      if (coarseInter_fct)
110
	(*coarseInter_fct)(drv, list, n, this);
Thomas Witkowski's avatar
Thomas Witkowski committed
111
    }
112
  
113
    /// Implements BasisFunction::getLocalIndices().
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
114
115
116
    void getLocalIndices(const Element *el,
			 const DOFAdmin *admin,
			 std::vector<DegreeOfFreedom> &dofs) const;
117

Thomas Witkowski's avatar
Thomas Witkowski committed
118
119
    void getLocalDofPtrVec(const Element *el, 
			   const DOFAdmin *admin,
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
120
			   vector<const DegreeOfFreedom*>& vec) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
121

122
    /// Implements BasisFunction::l2ScpFctBas
123
124
125
126
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<double, WorldVector<double> >* f,
		     DOFVector<double>* fh);

127
    /// Implements BasisFunction::l2ScpFctBas
128
129
130
131
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
		     DOFVector<WorldVector<double> >* fh);

132
133
    static void clear();

134
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
    /// sets the barycentric coordinates (stored in \ref bary) of the local 
    /// basis functions.
137
138
    void setBary();

139
    /// Recursive calculation of coordinates. Used by \ref setBary
140
    void createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
141
		      DimVec<double>* vec = NULL);
142

143
    /// Used by \ref setBary
144
145
    int** getIndexPermutations(int numIndices) const;

146
    /// Implements BasisFunction::setNDOF
147
148
    void setNDOF();

149
    /// Sets used function pointers
150
151
    void setFunctionPointer();

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
152
    /// Used by \ref getVec
153
154
    int* orderOfPositionIndices(const Element* el, 
				GeoIndex position, 
155
156
				int positionIndex) const;

157
158
    /// Calculates the number of DOFs needed for Lagrange of the given dim 
    /// and degree.
159
    static int getNumberOfDofs(int dim, int degree);
160
161

  private:
162
    /// barycentric coordinates of the locations of all basis functions
163
    std::vector<DimVec<double>* > *bary;
164
165
166
167

    /** \name static dim-degree-arrays
     * \{
     */
168
169
170
171
172
173
    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];
174
175
    /** \} */

Thomas Witkowski's avatar
Thomas Witkowski committed
176
177
    /// List of all used BasisFunctions in the whole program. Avoids duplicate
    /// instantiation of identical BasisFunctions. 
178
    static std::list<Lagrange*> allBasFcts;
179
180
181


  protected:
182
183
    /// Pointer to the used refineInter function
    void (*refineInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211

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

212
213
    /// Pointer to the used coarseRestr function
    void (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);
214
215
216
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

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

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

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

272

273
    /// AbstractFunction which implements lagrange basis functions
274
275
276
    class Phi : public BasFctType
    {
    public:
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
279
      /// 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.
280
      Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
281

282
      /// Destructor
283
284
285
      virtual ~Phi();

    private:
286
      /// vertices needed for evaluation of this function
287
288
      int* vertices;
    
289
      /// Pointer to the evaluating function
290
      double (*func)(const DimVec<double>& lambda, int* vert);
291

292
293
294
      /// Returns \ref func(lambda, vertices)
      inline double operator()(const DimVec<double>& lambda) const 
      {
295
	return func(lambda, vertices);
296
      }
297
298
299
300
301

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

302
      // ====== Lagrange, degree = 0 =====================================
303
      // center
304
305
      inline static double phi0c(const DimVec<double>&, int*) 
      {
306
	return 1.0;
307
      }
308

309
      // ====== Lagrange, degree = 1 =====================================
310
      // vertex
311
312
      inline static double phi1v(const DimVec<double>& lambda, int* vertices) 
      {
313
	return lambda[vertices[0]]; 
314
      }
315

316
      // ====== Lagrange, degree = 2 =====================================
317
      // vertex
318
319
      inline static double phi2v(const DimVec<double>& lambda, int* vertices) 
      {
320
	return lambda[vertices[0]] * (2.0 * lambda[vertices[0]] - 1.0);
321
      }    
322
323

      // edge
324
325
      inline static double phi2e(const DimVec<double>& lambda, int* vertices) 
      {
326
	return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]);
327
      }
328

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

      // edge
338
339
      inline static double phi3e(const DimVec<double>& lambda, int* vertices) 
      {
340
341
	return (13.5 * lambda[vertices[0]] - 4.5) * 
	  lambda[vertices[0]] * lambda[vertices[1]];
342
      }
343
344

      // face
345
346
      inline static double phi3f(const DimVec<double>& lambda, int* vertices) 
      {
347
348
	return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]];
349
      }
350

351
      // ====== Lagrange, degree = 4 ======================================
352
      // vertex
353
354
      inline static double phi4v(const DimVec<double>& lambda, int* vertices) 
      {
355
356
	return (((32.0 * lambda[vertices[0]] - 48.0) * lambda[vertices[0]] + 22.0)
	       * lambda[vertices[0]] - 3.0) * lambda[vertices[0]] / 3.0;
357
      }
358
359

      // edge
360
361
      inline static double phi4e0(const DimVec<double>& lambda, int* vertices) 
      {
362
363
	return ((128.0 * lambda[vertices[0]] - 96.0) * lambda[vertices[0]] + 16.0)
	  * lambda[vertices[0]] * lambda[vertices[1]] / 3.0;
364
      }
365

366
367
      inline static double phi4e1(const DimVec<double>& lambda, int* vertices) 
      {
368
369
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  (4.0 * lambda[vertices[1]] - 1.0) * lambda[vertices[1]] * 4.0;
370
      }
371
372

      // face
373
374
      inline static double phi4f(const DimVec<double>& lambda,  int* vertices) 
      {
375
376
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  lambda[vertices[1]] * lambda[vertices[2]] * 32.0;
377
      }
378
379

      // center
380
381
      inline static double phi4c(const DimVec<double>& lambda, int* vertices) 
      {
382
383
	return 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]] * lambda[vertices[3]];
384
      }
385
386
387
388
389

    };
  
    /** \} */

390
391


Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
    /// AbstractFunction which implements gradients of lagrange basis functions.
    /// See \ref Phi
394
395
396
    class GrdPhi : public GrdBasFctType
    {
    public:
397
398
      GrdPhi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);

399
400
401
      virtual ~GrdPhi();
    private:
      int* vertices;
402

403
404
405
      void (*func)(const DimVec<double>& lambda, 
		   int* vertices_, 
		   mtl::dense_vector<double>& result);
406

407
408
      inline void operator()(const DimVec<double>& lambda, 
			     mtl::dense_vector<double>& result) const 
409
      {
410
	func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
411
      }
412
413
414

      // ====== Lagrange0 ================================================
      // center
415
416
      inline static void grdPhi0c(const DimVec<double>&, 
				  int*, 
417
				  mtl::dense_vector<double>& result) 
418
      {
419
	result = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
420
      }
421
422
423

      // ====== Lagrange1 ================================================
      // vertex
424
425
      inline static void grdPhi1v(const DimVec<double>&, 
				  int* vertices, 
426
				  mtl::dense_vector<double>& result) 
427
      {
428
	result = 0.0;
429
	result[vertices[0]] = 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
430
      }
431
432
433

      // ====== Lagrange2 ================================================
      // vertex
434
435
      inline static void grdPhi2v(const DimVec<double>& lambda, 
				  int* vertices, 
436
				  mtl::dense_vector<double>& result) 
437
      {
438
	result = 0.0;
439
	result[vertices[0]] = 4.0 * lambda[vertices[0]] - 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
440
      }
441
442

      // edge
443
444
      inline static void grdPhi2e(const DimVec<double>& lambda, 
				  int* vertices, 
445
				  mtl::dense_vector<double>& result) 
446
      {
447
	result = 0.0;
448
449
	result[vertices[0]] = 4.0 * lambda[vertices[1]];
	result[vertices[1]] = 4.0 * lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
450
      }
451
452
453

      // ===== Lagrange3 ================================================
      // vertex
454
455
      inline static void grdPhi3v(const DimVec<double>& lambda,
				  int* vertices,
456
				  mtl::dense_vector<double>& result) 
457
      {
458
	result = 0.0;
459
460
	result[vertices[0]] = (13.5 * lambda[vertices[0]] - 9.0) * 
	  lambda[vertices[0]] + 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
461
      }
462
463

      // edge
464
465
      inline static void grdPhi3e(const DimVec<double>& lambda,
				  int* vertices, 
466
				  mtl::dense_vector<double>& result)
467
      {
468
	result = 0.0;
469
	result[vertices[0]] = (27.0 * lambda[vertices[0]] - 4.5) * 
470
	  lambda[vertices[1]];
471
	result[vertices[1]] = (13.5 * lambda[vertices[0]] - 4.5) * 
472
	  lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
473
      }
474
475

      // face
476
477
      inline static void grdPhi3f(const DimVec<double>& lambda, 
				  int* vertices, 
478
				  mtl::dense_vector<double>& result) 
479
      {
480
	result = 0.0;
481
482
483
	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
484
      }
485
486
487
488

    
      // ===== Lagrange4 ================================================
      // vertex
489
490
      inline static void grdPhi4v(const DimVec<double>& lambda, 
				  int* vertices, 
491
				  mtl::dense_vector<double>& result) 
492
      {
493
	result = 0.0;
494
495
496
	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
497
      }
498
499
    
      // edge
500
501
      inline static void grdPhi4e0(const DimVec<double>& lambda,
				   int* vertices, 
502
				   mtl::dense_vector<double>& result)
503
      {
504
	result = 0.0;
505
506
507
508
	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
509
      }
510
511
512

      inline static void grdPhi4e1(const DimVec<double>& lambda,
				   int* vertices,
513
				   mtl::dense_vector<double>& result)
514
      {
515
	result = 0.0;
516
517
518
519
	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
520
      }
521
522

      // face
523
524
      inline static void grdPhi4f(const DimVec<double>& lambda, 
				  int* vertices, 
525
				  mtl::dense_vector<double>& result)
526
      {
527
	result = 0.0;
528
529
530
531
532
533
	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
534
      }
535
536
    
      // center
537
538
      inline static void grdPhi4c(const DimVec<double>& lambda, 
				  int* vertices,
539
				  mtl::dense_vector<double>& result) 
540
      {
541
542
543
544
545
546
547
548
549
	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
550
      }
551
552
    };

553
554


Thomas Witkowski's avatar
Thomas Witkowski committed
555
556
    /// AbstractFunction which implements second derivatives of Lagrange basis
    /// functions. See \ref Phi
557
558
559
    class D2Phi : public D2BasFctType
    {
    public:
560
      D2Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
561
562
563
564
565

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

568
569
      inline void operator()(const DimVec<double>& lambda, DimMat<double>& result) const {
	return func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
570
      }
571
572
573

      // ===== Lagrange0 ================================================
      // center
574
575
576
      inline static void D2Phi0c(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);       
Thomas Witkowski's avatar
Thomas Witkowski committed
577
      }
578
579
580

      // ===== Lagrange1 ================================================
      // vertex
581
582
583
      inline static void D2Phi1v(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
584
      }
585
586
587

      // ===== Lagrange2 ================================================
      // vertex
588
589
590
591
592
      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
593
      }
594
595

      // edge
596
597
598
599
600
601
      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
602
      }
603
604
605
606


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

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

      // face
625
626
627
628
629
630
631
632
633
634
      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
635
      }
636
637
638
639


      // ===== Lagrange4 ================================================
      // vertex
640
641
642
643
644
645
      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
646
      }
647
648

      // edge
649
650
651
652
653
654
655
656
657
      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
658
      }
659
660
661
662
663
664
665
666
667
668
669
670

      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
671
      }
672
673

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

      // center
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
      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
714
      }
715
716
717
718
719
720
    };
  };

}

#endif // AMDIS_LAGRANGE_H