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
    void interpol(const ElInfo *, int, const int *, 
		  AbstractFunction<double, WorldVector<double> >*, 
		  mtl::dense_vector<double>&) const;
70

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

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

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

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

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

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

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

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

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

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

141
142
    static void clear();

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

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

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

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

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

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

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

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

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

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


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

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

224
225
    /// Pointer to the used coarseRestr function
    void (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

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

254
255
    /// Pointer to the used coarseInter function
    void (*coarseInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*);  
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283

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

284

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    };
  
    /** \} */

404
405


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

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

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

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

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

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

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

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

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

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

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

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

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

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

569
570


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

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

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

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

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

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

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


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

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

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


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

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

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

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

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

}

#endif // AMDIS_LAGRANGE_H