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
    /// Calculates the number of DOFs needed for Lagrange of the given dim and degree.
158
    static int getNumberOfDofs(int dim, int degree);
159
160

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

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

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


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

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

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

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

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

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

271

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    };
  
    /** \} */

389
390


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

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

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

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

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

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

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

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

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

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

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

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

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

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

552
553


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

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

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

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

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

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

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


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

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

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


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

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

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

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

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

}

#endif // AMDIS_LAGRANGE_H