Lagrange.h 25.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file Lagrange.h */

#ifndef AMDIS_LAGRANGE_H
#define AMDIS_LAGRANGE_H

#include "BasisFunction.h"
#include "FixVec.h"
#include "AbstractFunction.h"
#include <list>

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

    /** \brief
     * Implements BasisFunction::interpol
     */
67
68
69
    const double *interpol(const ElInfo *, int, const int *, 
			   AbstractFunction<double, WorldVector<double> >*, 
			   double *);
70
71
72
73
74


    /** \brief
     * Implements BasisFunction::interpol
     */
75
76
77
78
    const WorldVector<double> *interpol(const ElInfo *, int, 
					const int *b_no,
					AbstractFunction<WorldVector<double>, 
					WorldVector<double> >*, 
79
80
81
82
83
84
85
86
87
88
					WorldVector<double> *);

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

    /** \brief
     * Implements BasisFunction::getDOFIndices
     */
89
90
    const DegreeOfFreedom* getDOFIndices(const Element*,
					 const DOFAdmin&, 
91
92
93
94
95
					 DegreeOfFreedom*) const;

    /** \brief
     * Implements BasisFunction::getBound
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
96
    void getBound(const ElInfo*, BoundaryType *) const;
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

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

    /** \brief
     * Implements BasisFunction::refineInter
     */
    inline void  refineInter(DOFIndexed<double> *drv,
			     RCNeighbourList* list, 
			     int n) 
    {
118
      if (refineInter_fct)
119
	(*refineInter_fct)(drv, list, n, this);
120
    }
121
122
123
124
125
126
127
128

    /** \brief
     * Implements BasisFunction::coarseRestrict
     */
    inline void  coarseRestr(DOFIndexed<double> *drv,
			     RCNeighbourList* list,
			     int n)
    {
129
      if (coarseRestr_fct)
130
	(*coarseRestr_fct)(drv, list, n, this);
131
    }
132
133
134
135
136
137
138
139
  
    /** \brief
     * Implements BasisFunction::coarseInter
     */
    inline void  coarseInter(DOFIndexed<double> *drv, 
			     RCNeighbourList* list, 
			     int n) 
    {
140
      if (coarseInter_fct)
141
	(*coarseInter_fct)(drv, list, n, this);
Thomas Witkowski's avatar
Thomas Witkowski committed
142
    }
143
144
145
146
147
148
149
150
  
    /** \brief
     * Implements BasisFunction::getLocalIndices().
     */
    const DegreeOfFreedom *getLocalIndices(const Element*,
					   const DOFAdmin*,
					   DegreeOfFreedom*) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
154
155

    void getLocalIndicesVec(const Element*,
			    const DOFAdmin*,
			    Vector<DegreeOfFreedom>*) const;

156
157
158
159
160
161
162
163
164
165
166
167
168
169
    /** \brief
     * Implements BasisFunction::l2ScpFctBas
     */
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<double, WorldVector<double> >* f,
		     DOFVector<double>* fh);

    /** \brief
     * Implements BasisFunction::l2ScpFctBas
     */
    void l2ScpFctBas(Quadrature* q,
		     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
		     DOFVector<WorldVector<double> >* fh);

170
171
    static void clear();

172
173
174
175
176
177
178
179
180
181
182
  protected:
    /** \brief
     * sets the barycentric coordinates (stored in \ref bary) of the local 
     * basis functions.
     */
    void setBary();

    /** \brief
     * Recursive calculation of coordinates. Used by \ref setBary
     */
    void createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
183
		      DimVec<double>* vec = NULL);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202

    /** \brief
     * Used by \ref setBary
     */
    int** getIndexPermutations(int numIndices) const;

    /** \brief
     * Implements BasisFunction::setNDOF
     */
    void setNDOF();

    /** \brief
     * Sets used function pointers
     */
    void setFunctionPointer();

    /** \brief
     * Used by \ref getDOFIndices and \ref getVec
     */
203
204
    int* orderOfPositionIndices(const Element* el, 
				GeoIndex position, 
205
206
207
208
209
210
211
212
213
214
215
216
				int positionIndex) const;

    /** \brief
     * Calculates the number of DOFs needed for Lagrange of the given dim and
     * degree.
     */
    static int getNumberOfDOFs(int dim, int degree);

  private:
    /** \brief
     * barycentric coordinates of the locations of all basis functions
     */
217
    std::vector<DimVec<double>* > *bary;
218
219
220
221

    /** \name static dim-degree-arrays
     * \{
     */
222
223
224
225
226
227
    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];
228
229
230
231
232
233
    /** \} */

    /** \brief
     * List of all used BasisFunctions in the whole program. Avoids duplicate
     * instantiation of identical BasisFunctions. 
     */
234
    static std::list<Lagrange*> allBasFcts;
235
236
237
238
239
240
241
242
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
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336


  protected:
    /** \brief
     * Pointer to the used refineInter function
     */
    void  (*refineInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, 
			     BasisFunction*);

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

    /** \brief
     * Pointer to the used coarseRestr function
     */
    void  (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, 
			     BasisFunction*);

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

    /** \brief
     * Pointer to the used coarseInter function
     */
    void  (*coarseInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, 
			     BasisFunction*);  

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

337
338


339
340
341
342
343
344
345
346
347
348
349
350
351
    // ===== Phi =========================================================

    /** \brief
     * AbstractFunction which implements lagrange basis functions
     */
    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.
       */
352
      Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367

      /** \brief
       * Destructor
       */
      virtual ~Phi();

    private:
      /** \brief
       * vertices needed for evaluation of this function
       */
      int* vertices;
    
      /** \brief
       * Pointer to the evaluating function
       */
368
      double (*func)(const DimVec<double>& lambda, int* vert);
369
370
371
372

      /** \brief
       * Returns \ref func(lambda, vertices)
       */
373
      inline double operator()(const DimVec<double>& lambda) const {
374
	return func(lambda, vertices);
375
      }
376
377
378
379
380

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

381
      // ====== Lagrange, degree = 0 =====================================
382
      // center
383
384
      inline static double phi0c(const DimVec<double>&, int*) {
	return 1.0;
385
      }
386

387
      // ====== Lagrange, degree = 1 =====================================
388
      // vertex
389
390
      inline static double phi1v(const DimVec<double>& lambda, int* vertices) {
	return lambda[vertices[0]]; 
391
      }
392

393
      // ====== Lagrange, degree = 2 =====================================
394
      // vertex
395
396
      inline static double phi2v(const DimVec<double>& lambda, int* vertices) {
	return lambda[vertices[0]] * (2.0 * lambda[vertices[0]] - 1.0);
397
      }    
398
399

      // edge
400
401
      inline static double phi2e(const DimVec<double>& lambda, int* vertices) {
	return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]);
402
      }
403

404
      // ====== Lagrange, degree = 3 =====================================
405
      // vertex
406
407
      inline static double phi3v(const DimVec<double>& lambda, int* vertices) {
	return (4.5 * (lambda[vertices[0]] - 1.0) * lambda[vertices[0]] + 1.0) * 
408
	  lambda[vertices[0]];
409
      }
410
411

      // edge
412
413
414
      inline static double phi3e(const DimVec<double>& lambda, int* vertices) {
	return (13.5 * lambda[vertices[0]] - 4.5) * 
	  lambda[vertices[0]] * lambda[vertices[1]];
415
      }
416
417

      // face
418
419
420
      inline static double phi3f(const DimVec<double>& lambda, int* vertices) {
	return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]];
421
      }
422

423
      // ====== Lagrange, degree = 4 ======================================
424
      // vertex
425
426
427
      inline static double phi4v(const DimVec<double>& lambda, int* vertices) {
	return (((32.0 * lambda[vertices[0]] - 48.0) * lambda[vertices[0]] + 22.0)
	       * lambda[vertices[0]] - 3.0) * lambda[vertices[0]] / 3.0;
428
      }
429
430

      // edge
431
432
433
      inline static double phi4e0(const DimVec<double>& lambda, int* vertices) {
	return ((128.0 * lambda[vertices[0]] - 96.0) * lambda[vertices[0]] + 16.0)
	  * lambda[vertices[0]] * lambda[vertices[1]] / 3.0;
434
      }
435

436
437
438
      inline static double phi4e1(const DimVec<double>& lambda, int* vertices) {
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  (4.0 * lambda[vertices[1]] - 1.0) * lambda[vertices[1]] * 4.0;
439
      }
440
441

      // face
442
443
444
      inline static double phi4f(const DimVec<double>& lambda,  int* vertices) {
	return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * 
	  lambda[vertices[1]] * lambda[vertices[2]] * 32.0;
445
      }
446
447

      // center
448
449
450
      inline static double phi4c(const DimVec<double>& lambda, int* vertices) {
	return 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * 
	  lambda[vertices[2]] * lambda[vertices[3]];
451
      }
452
453
454
455
456

    };
  
    /** \} */

457
458


459
460
461
462
463
464
465
466
467
    // ===== GrdPhi ======================================================

    /** \brief
     * AbstractFunction which implements gradients of lagrange basis functions.
     * See \ref Phi
     */
    class GrdPhi : public GrdBasFctType
    {
    public:
468
469
      GrdPhi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);

470
471
472
      virtual ~GrdPhi();
    private:
      int* vertices;
473
474
475
      void (*func)(const DimVec<double>& lambda,
		   int* vertices_,
		   DimVec<double>& result);
476

477
478
      inline void operator()(const DimVec<double>& lambda, DimVec<double>& result) const {
	func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
479
      }
480
481
482

      // ====== Lagrange0 ================================================
      // center
483
484
485
486
487
      inline static void grdPhi0c(const DimVec<double>&, 
				  int*, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
488
      }
489
490
491

      // ====== Lagrange1 ================================================
      // vertex
492
493
494
495
496
497
      inline static void grdPhi1v(const DimVec<double>&, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
498
      }
499
500
501

      // ====== Lagrange2 ================================================
      // vertex
502
503
504
505
506
507
      inline static void grdPhi2v(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 4.0 * lambda[vertices[0]] - 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
508
      }
509
510

      // edge
511
512
513
514
515
516
517
      inline static void grdPhi2e(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = 4.0 * lambda[vertices[1]];
	result[vertices[1]] = 4.0 * lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
518
      }
519
520
521

      // ===== Lagrange3 ================================================
      // vertex
522
523
524
525
526
527
528
      inline static void grdPhi3v(const DimVec<double>& lambda,
				  int* vertices,
				  DimVec<double>& result) 
      {
	result.set(0.0);
	result[vertices[0]] = (13.5 * lambda[vertices[0]] - 9.0) * 
	  lambda[vertices[0]] + 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
529
      }
530
531

      // edge
532
533
534
535
536
537
      inline static void grdPhi3e(const DimVec<double>& lambda,
				  int* vertices, 
				  DimVec<double>& result)
      {
	result.set(0.0);
	result[vertices[0]] = (27.0 * lambda[vertices[0]] - 4.5) * 
538
	  lambda[vertices[1]];
539
	result[vertices[1]] = (13.5 * lambda[vertices[0]] - 4.5) * 
540
	  lambda[vertices[0]];
Thomas Witkowski's avatar
Thomas Witkowski committed
541
      }
542
543

      // face
544
545
546
547
548
549
550
551
      inline static void grdPhi3f(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	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
552
      }
553
554
555
556

    
      // ===== Lagrange4 ================================================
      // vertex
557
558
559
560
561
562
563
564
      inline static void grdPhi4v(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result) 
      {
	result.set(0.0);
	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
565
      }
566
567
    
      // edge
568
569
570
571
572
573
574
575
576
      inline static void grdPhi4e0(const DimVec<double>& lambda,
				   int* vertices, 
				   DimVec<double>& result)
      {
	result.set(0.0);
	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
577
      }
578
579
580
581
582
583
584
585
586
587

      inline static void grdPhi4e1(const DimVec<double>& lambda,
				   int* vertices,
				   DimVec<double>& result)
      {
	result.set(0.0);
	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
588
      }
589
590

      // face
591
592
593
594
595
596
597
598
599
600
601
      inline static void grdPhi4f(const DimVec<double>& lambda, 
				  int* vertices, 
				  DimVec<double>& result)
      {
	result.set(0.0);
	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
602
      }
603
604
    
      // center
605
606
607
608
609
610
611
612
613
      inline static void grdPhi4c(const DimVec<double>& lambda, 
				  int* vertices,
				  DimVec<double>& result) 
      {
	result.set(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
614
      }
615
616
    };

617
618


619
620
621
622
623
624
625
626
627
628
    
    // ===== D2Phi =======================================================

    /** \brief
     * AbstractFunction which implements second derivatives of Lagrange basis
     * functions. See \ref Phi
     */
    class D2Phi : public D2BasFctType
    {
    public:
629
      D2Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex);
630
631
632
633
634

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

637
638
      inline void operator()(const DimVec<double>& lambda, DimMat<double>& result) const {
	return func(lambda, vertices, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
639
      }
640
641
642

      // ===== Lagrange0 ================================================
      // center
643
644
645
      inline static void D2Phi0c(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);       
Thomas Witkowski's avatar
Thomas Witkowski committed
646
      }
647
648
649

      // ===== Lagrange1 ================================================
      // vertex
650
651
652
      inline static void D2Phi1v(const DimVec<double>&, int*, DimMat<double>& result) 
      {
	result.set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
653
      }
654
655
656

      // ===== Lagrange2 ================================================
      // vertex
657
658
659
660
661
      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
662
      }
663
664

      // edge
665
666
667
668
669
670
      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
671
      }
672
673
674
675


      // ===== Lagrange3 ================================================
      // vertex
676
677
678
679
680
      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
681
      }
682
683

      // edge
684
685
686
687
688
689
690
      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
691
      }
692
693

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


      // ===== Lagrange4 ================================================
      // vertex
709
710
711
712
713
714
      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
715
      }
716
717

      // edge
718
719
720
721
722
723
724
725
726
      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
727
      }
728
729
730
731
732
733
734
735
736
737
738
739

      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
740
      }
741
742

      // face
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
      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
758
      }
759
760

      // center
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
      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
783
      }
784
785
786
787
788
789
    };
  };

}

#endif // AMDIS_LAGRANGE_H