Lagrange.h 24.1 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::getDOFIndices
77
78
    const DegreeOfFreedom* getDOFIndices(const Element*,
					 const DOFAdmin&, 
79
80
					 DegreeOfFreedom*) const;

81
    /// Implements BasisFunction::getBound
Thomas Witkowski's avatar
Thomas Witkowski committed
82
    void getBound(const ElInfo*, BoundaryType *) const;
83
84
85
86
87
88
89
90
91
92
93
94
95
96

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

97
98
    /// Implements BasisFunction::refineInter
    inline void refineInter(DOFIndexed<double> *drv, RCNeighbourList* list, int n) 
99
    {
100
      if (refineInter_fct)
101
	(*refineInter_fct)(drv, list, n, this);
102
    }
103

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

Thomas Witkowski's avatar
Thomas Witkowski committed
123
124
125
    void getLocalDofPtrVec(const Element *el, 
			   const DOFAdmin *admin,
			   std::vector<const DegreeOfFreedom*>& vec) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
126

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

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

137
138
    static void clear();

139
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
    /// sets the barycentric coordinates (stored in \ref bary) of the local 
    /// basis functions.
142
143
    void setBary();

144
    /// Recursive calculation of coordinates. Used by \ref setBary
145
    void createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
146
		      DimVec<double>* vec = NULL);
147

148
    /// Used by \ref setBary
149
150
    int** getIndexPermutations(int numIndices) const;

151
    /// Implements BasisFunction::setNDOF
152
153
    void setNDOF();

154
    /// Sets used function pointers
155
156
    void setFunctionPointer();

157
    /// Used by \ref getDOFIndices and \ref getVec
158
159
    int* orderOfPositionIndices(const Element* el, 
				GeoIndex position, 
160
161
				int positionIndex) const;

162
    /// Calculates the number of DOFs needed for Lagrange of the given dim and degree.
163
    static int getNumberOfDofs(int dim, int degree);
164
165

  private:
166
    /// barycentric coordinates of the locations of all basis functions
167
    std::vector<DimVec<double>* > *bary;
168
169
170
171

    /** \name static dim-degree-arrays
     * \{
     */
172
173
174
175
176
177
    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];
178
179
    /** \} */

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


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

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

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

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

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

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

276

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

286
      /// Destructor
287
288
289
      virtual ~Phi();

    private:
290
      /// vertices needed for evaluation of this function
291
292
      int* vertices;
    
293
      /// Pointer to the evaluating function
294
      double (*func)(const DimVec<double>& lambda, int* vert);
295

296
297
298
      /// Returns \ref func(lambda, vertices)
      inline double operator()(const DimVec<double>& lambda) const 
      {
299
	return func(lambda, vertices);
300
      }
301
302
303
304
305

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

306
      // ====== Lagrange, degree = 0 =====================================
307
      // center
308
309
      inline static double phi0c(const DimVec<double>&, int*) 
      {
310
	return 1.0;
311
      }
312

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

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

      // edge
328
329
      inline static double phi2e(const DimVec<double>& lambda, int* vertices) 
      {
330
	return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]);
331
      }
332

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

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

      // face
349
350
      inline static double phi3f(const DimVec<double>& lambda, int* vertices) 
      {
351
352
	return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]];
353
      }
354

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

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

370
371
      inline static double phi4e1(const DimVec<double>& lambda, int* vertices) 
      {
372
373
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  (4.0 * lambda[vertices[1]] - 1.0) * lambda[vertices[1]] * 4.0;
374
      }
375
376

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

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

    };
  
    /** \} */

394
395


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

403
404
405
      virtual ~GrdPhi();
    private:
      int* vertices;
406

407
408
409
      void (*func)(const DimVec<double>& lambda, 
		   int* vertices_, 
		   mtl::dense_vector<double>& result);
410

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

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

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

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

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

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

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

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

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

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

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

557
558


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

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

572
573
      inline void operator()(const DimVec<double>& lambda, DimMat<double>& result) const {
	return func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
574
      }
575
576
577

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

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

      // ===== Lagrange2 ================================================
      // vertex
592
593
594
595
596
      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
597
      }
598
599

      // edge
600
601
602
603
604
605
      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
606
      }
607
608
609
610


      // ===== Lagrange3 ================================================
      // vertex
611
612
613
614
615
      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
616
      }
617
618

      // edge
619
620
621
622
623
624
625
      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
626
      }
627
628

      // face
629
630
631
632
633
634
635
636
637
638
      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
639
      }
640
641
642
643


      // ===== Lagrange4 ================================================
      // vertex
644
645
646
647
648
649
      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
650
      }
651
652

      // edge
653
654
655
656
657
658
659
660
661
      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
662
      }
663
664
665
666
667
668
669
670
671
672
673
674

      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
675
      }
676
677

      // face
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
      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
693
      }
694
695

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

}

#endif // AMDIS_LAGRANGE_H