Lagrange.cc 179 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include "Mesh.h"
#include "Element.h"
#include "Lagrange.h"
#include "DOFAdmin.h"
#include "RCNeighbourList.h"
#include "DOFVector.h"
#include "Traverse.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "Parametric.h"

#include <stdio.h>
#include <algorithm>
#include <list>

namespace AMDiS {

19
  std::vector<DimVec<double>* > Lagrange::baryDimDegree[MAX_DIM+1][MAX_DEGREE+1];
20
21
  DimVec<int>* Lagrange::ndofDimDegree[MAX_DIM+1][MAX_DEGREE+1];
  int Lagrange::nBasFctsDimDegree[MAX_DIM+1][MAX_DEGREE+1];
22
23
24
25
  std::vector<BasFctType*> Lagrange::phiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
  std::vector<GrdBasFctType*> Lagrange::grdPhiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
  std::vector<D2BasFctType*> Lagrange::D2PhiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
  std::list<Lagrange*> Lagrange::allBasFcts;
26
27
28


  Lagrange::Lagrange(int dim, int degree)
29
    : BasisFunction(std::string("Lagrange"), dim, degree)
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
  {
    // set name
    char dummy[3];
    sprintf(dummy, "%d", dim);
    name.append(dummy);
    name.append(" ");
    sprintf(dummy, "%d", degree);  
    name.append(dummy);

    // set nDOF
    setNDOF();

    // set barycentric coordinates
    setBary();

    // set function pointer
    setFunctionPointer();

  }

  Lagrange::~Lagrange()
  {
  }

  Lagrange* Lagrange::getLagrange(int dim, int degree)
  {
56
    std::list<Lagrange*>::iterator it;
57
58
59
60
61
62
63
64
65
66
67
68
69
70

    for(it = allBasFcts.begin(); it != allBasFcts.end(); it++) {
      if(((*it)->dim == dim) && ((*it)->degree == degree)) {
	return (*it);
      } 
    }

    Lagrange* newLagrange = NEW Lagrange(dim, degree);
    allBasFcts.push_back(newLagrange);
    return newLagrange;
  }

  void Lagrange::setFunctionPointer()
  {
71
    if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) {
72
      // for all positions
73
      for (int i = 0; i < dim+1; i++) {
74
	// no vertex dofs for degree 0 ?
75
76
	if (degree == 0 && i != dim) 
	  continue;
77
	// for all vertices/edges/...
78
	for (int j = 0; j < Global::getGeo(INDEX_OF_DIM(i, dim),dim); j++) {
79
	  // for all dofs at this position
80
	  for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
81
82
	    // basis function
	    phiDimDegree[dim][degree].push_back(NEW Phi(this, 
83
							INDEX_OF_DIM(i ,dim),j, k));
84
85
	    // gradients
	    grdPhiDimDegree[dim][degree].push_back(NEW GrdPhi(this, 
86
87
							      INDEX_OF_DIM(i, dim),
							      j, k));
88
89
	    // D2-Matrices
	    D2PhiDimDegree[dim][degree].push_back(NEW D2Phi(this, 
90
91
							    INDEX_OF_DIM(i, dim),
							    j, k));
92
93
94
95
	  }
	}
      }
    }
96
    phi = &phiDimDegree[dim][degree];
97
    grdPhi = &grdPhiDimDegree[dim][degree];
98
    d2Phi = &D2PhiDimDegree[dim][degree];
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177

    switch(degree) {
    case 0:
      refineInter_fct = refineInter0;
      coarseRestr_fct = coarseRestr0;
      coarseInter_fct = coarseInter0;  // not yet implemented
      break;
    case 1:
      refineInter_fct = refineInter1;
      coarseRestr_fct = coarseRestr1;
      coarseInter_fct = NULL;  // not yet implemented
      break;
    case 2:
      switch(dim) {
      case 1:
	refineInter_fct = refineInter2_1d;
	coarseRestr_fct = NULL; // not yet implemented
	coarseInter_fct = coarseInter2_1d;
	break;
      case 2: 
	refineInter_fct = refineInter2_2d;
	coarseRestr_fct = coarseRestr2_2d;
	coarseInter_fct = coarseInter2_2d;
	break;
      case 3:
	refineInter_fct = refineInter2_3d;
	coarseRestr_fct = coarseRestr2_3d;
	coarseInter_fct = coarseInter2_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 3:
      switch(dim) {
      case 1:
	refineInter_fct = refineInter3_1d;
	coarseRestr_fct = coarseRestr3_1d;
	coarseInter_fct = coarseInter3_1d;
	break;
      case 2: 
	refineInter_fct = refineInter3_2d;
	coarseRestr_fct = coarseRestr3_2d;
	coarseInter_fct = coarseInter3_2d;
	break;
      case 3:
	refineInter_fct = refineInter3_3d;
	coarseRestr_fct = coarseRestr3_3d;
	coarseInter_fct = coarseInter3_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 4:
      switch(dim) {
      case 1: 
	refineInter_fct = refineInter4_1d;
	coarseRestr_fct = coarseRestr4_1d;
	coarseInter_fct = coarseInter4_1d;
	break;
      case 2:
	refineInter_fct = refineInter4_2d;
	coarseRestr_fct = coarseRestr4_2d;
	coarseInter_fct = coarseInter4_2d;
	break;
      case 3:
	refineInter_fct = refineInter4_3d;
	coarseRestr_fct = coarseRestr4_3d;
	coarseInter_fct = coarseInter4_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

  void Lagrange::setNDOF()
  {
178
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
179
180
      ndofDimDegree[dim][degree] = NEW DimVec<int>(dim, DEFAULT_VALUE, 0);

181
182
183
184
185
      if (degree != 0) {
	(*ndofDimDegree[dim][degree])[VERTEX] = 1;
      } else {
	(*ndofDimDegree[dim][degree])[VERTEX] = 0;
      }
186

187
      for (int i = 1; i < dim + 1; i++) {
188
189
	nBasFcts = getNumberOfDOFs(i, degree);
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
190
	for (int j = 0; j < i; j++) {
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
	  (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] -= 
	    Global::getGeo(INDEX_OF_DIM(j, dim), i) * 
	    (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(j, dim)];
	}
      }
      nBasFctsDimDegree[dim][degree] = nBasFcts;
    }
    nBasFcts = nBasFctsDimDegree[dim][degree];
    nDOF = ndofDimDegree[dim][degree];
  } 

  DimVec<double> *Lagrange::getCoords(int i) const
  {
    return (*bary)[i];
  }

  void Lagrange::setVertices(int dim, int degree, 
			     GeoIndex position, int positionIndex, int nodeIndex, 
			     int** vertices)
  {
211
    TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");
212
213
214
215
    int dimOfPosition = DIM_OF_INDEX(position, dim);

    *vertices = GET_MEMORY(int, dimOfPosition + 1);

216
    if ((degree == 4) && (dimOfPosition==1)) {
217
218
219
220
221
222
223
224
225
      (*vertices)[(nodeIndex != 2) ? 0 : 1] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex,
							      0);
      (*vertices)[(nodeIndex != 2) ? 1 : 0] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex, 
							      1);
    } else if ((degree==4) && (dimOfPosition==2)) {
226
227
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+dimOfPosition*nodeIndex) % (dimOfPosition + 1)] = 
228
229
230
231
232
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }    
    } else {
233
234
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+nodeIndex) % (dimOfPosition + 1)] = 
235
236
237
238
239
240
241
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }
    }
  }

242
  Lagrange::Phi::Phi(Lagrange* owner, 
243
244
245
		     GeoIndex position, 
		     int positionIndex, 
		     int nodeIndex)
246
    : vertices(NULL)
247
248
249
250
251
252
253
254
255
256
257
  {
    FUNCNAME("Lagrange::Phi::Phi")
      // get relevant vertices
      Lagrange::setVertices(owner->getDim(), 
			    owner->getDegree(), 
			    position, 
			    positionIndex, 
			    nodeIndex, 
			    &vertices);
  
    // set function pointer
258
    switch (owner->getDegree()) {
259
    case 0:
260
      switch (position) {
261
262
263
264
265
266
267
268
      case CENTER:
	func = phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
269
      switch (position) {
270
271
272
273
274
275
276
277
      case VERTEX:
	func = phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
278
      switch (position) {
279
280
281
282
      case VERTEX:
	func = phi2v;
	break;
      case EDGE:
283
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
284
285
286
	func = phi2e;
	break;
      case CENTER:
287
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
288
289
290
291
292
293
294
	func = phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
295
      switch (position) {
296
297
298
299
300
301
302
      case VERTEX:
	func = phi3v;
	break;
      case EDGE:
	func = phi3e;
	break;
      case FACE:
303
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
304
305
306
	func = phi3f;
	break;
      case CENTER:
307
	switch (owner->getDim()) {
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
	case 1:
	  func = phi3e;
	  break;
	case 2:
	  func = phi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
324
      switch (position) {
325
326
327
328
      case VERTEX:
	func = phi4v;
	break;
      case EDGE:
329
	if (nodeIndex == 1)
330
331
332
333
334
	  func = phi4e1;
	else 
	  func = phi4e0;
	break;
      case FACE:
335
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
336
337
338
	func = phi4f;      
	break;
      case CENTER:
339
	switch (owner->getDim()) {
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
	case 1:
	  if(nodeIndex == 1)
	    func = phi4e1;
	  else 
	    func = phi4e0;
	  break;
	case 2:
	  func = phi4f;
	  break;
	case 3:
	  func = phi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

366
367
368
  Lagrange::Phi::~Phi() { 
    DELETE [] vertices; 
  }
369
370


371
  Lagrange::GrdPhi::GrdPhi(Lagrange* owner, 
372
373
374
			   GeoIndex position, 
			   int positionIndex, 
			   int nodeIndex)
375
    : vertices(NULL)
376
377
378
379
380
381
382
383
384
385
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
386
    switch (owner->getDegree()) {
387
    case 0:
388
      switch (position) {
389
390
391
392
393
394
395
396
      case CENTER:
	func = grdPhi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
397
      switch (position) {
398
399
400
401
402
403
404
405
      case VERTEX:
	func = grdPhi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
406
      switch (position) {
407
408
409
410
411
412
413
      case VERTEX:
	func = grdPhi2v;
	break;
      case EDGE:
	func = grdPhi2e;
	break;
      case CENTER:
414
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
415
416
417
418
419
420
421
	func = grdPhi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
422
      switch (position) {
423
424
425
426
427
428
429
      case VERTEX:
	func = grdPhi3v;
	break;
      case EDGE:
	func = grdPhi3e;
	break;
      case FACE:
430
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
431
432
433
	func = grdPhi3f;
	break;
      case CENTER:
434
	switch (owner->getDim()) {
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
	case 1:
	  func = grdPhi3e;
	  break;
	case 2:
	  func = grdPhi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
451
      switch (position) {
452
453
454
455
      case VERTEX:
	func = grdPhi4v;
	break;
      case EDGE:
456
	if (nodeIndex == 1)
457
458
459
460
461
	  func = grdPhi4e1;
	else 
	  func = grdPhi4e0;
	break;
      case FACE:
462
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
463
464
465
	func = grdPhi4f;      
	break;
      case CENTER:
466
	switch (owner->getDim()) {
467
	case 1:
468
	  if (nodeIndex == 1)
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
	    func = grdPhi4e1;
	  else 
	    func = grdPhi4e0;
	  break;
	case 2:
	  func = grdPhi4f;
	  break;
	case 3:
	  func = grdPhi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

493
494
495
  Lagrange::GrdPhi::~GrdPhi() { 
    DELETE [] vertices; 
  }
496

497
  Lagrange::D2Phi::D2Phi(Lagrange* owner, 
498
499
500
			 GeoIndex position, 
			 int positionIndex, 
			 int nodeIndex)
501
    : vertices(NULL)
502
503
504
505
506
507
508
509
510
511
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
512
    switch (owner->getDegree()) {
513
    case 0:
514
      switch (position) {
515
516
517
518
519
520
521
522
      case CENTER:
	func = D2Phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
523
      switch (position) {
524
525
526
527
528
529
530
531
      case VERTEX:
	func = D2Phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
532
      switch (position) {
533
534
535
536
      case VERTEX:
	func = D2Phi2v;
	break;
      case EDGE:
537
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
538
539
540
	func = D2Phi2e;
	break;
      case CENTER:
541
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
542
543
544
545
546
547
548
	func = D2Phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
549
      switch (position) {
550
551
552
553
554
555
556
      case VERTEX:
	func = D2Phi3v;
	break;
      case EDGE:
	func = D2Phi3e;
	break;
      case FACE:
557
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
558
559
560
	func = D2Phi3f;
	break;
      case CENTER:
561
	switch (owner->getDim()) {
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
	case 1:
	  func = D2Phi3e;
	  break;
	case 2:
	  func = D2Phi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
578
      switch (position) {
579
580
581
582
583
584
585
586
587
588
      case VERTEX:
	func = D2Phi4v;
	break;
      case EDGE:
	if(nodeIndex == 1)
	  func = D2Phi4e1;
	else 
	  func = D2Phi4e0;
	break;
      case FACE:
589
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
590
591
592
	func = D2Phi4f;      
	break;
      case CENTER:
593
	switch (owner->getDim()) {
594
	case 1:
595
	  if (nodeIndex == 1)
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
	    func = D2Phi4e1;
	  else 
	    func = D2Phi4e0;
	  break;
	case 2:
	  func = D2Phi4f;
	  break;
	case 3:
	  func = D2Phi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

620
621
622
  Lagrange::D2Phi::~D2Phi() { 
    DELETE [] vertices; 
  }
623
624
625
626
627
628
629

  void Lagrange::createCoords(int* coordInd, 
			      int numCoords, 
			      int dimIndex, 
			      int rest, 
			      DimVec<double>* vec)
  {
630
631
632
    if (vec == NULL) {
      vec = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
    }
633

634
    if (dimIndex == numCoords-1) {
635
636
637
638
      (*vec)[coordInd[dimIndex]] = double(rest)/degree;
      DimVec<double>* newCoords = NEW DimVec<double>(*vec);
      bary->push_back(newCoords);
    } else {
639
      for (int i = rest - 1; i >= 1; i--) {
640
	(*vec)[coordInd[dimIndex]] = double(i) / degree;
641
	createCoords(coordInd, numCoords, dimIndex + 1, rest - i, vec);
642
643
644
645
646
647
648
      }
    }
  }

  void Lagrange::setBary()
  {
    bary  =  &baryDimDegree[dim][degree];
649
650
651

    if (static_cast<int>(bary->size()) == 0) {
      for (int i = 0; i <= dim; i++) { // for all positions
652
	int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
653
654
655
	for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
	  int *coordInd = GET_MEMORY(int, i + 1); // indices of relevant coords
	  for (int k = 0; k < i+1; k++) { 
656
657
658
	    coordInd[k] = Global::getReferenceElement(dim)->
	      getVertexOfPosition(INDEX_OF_DIM(i,dim), j, k);
	  }
659
660
661
	  createCoords(coordInd, i + 1, 0, degree);
	  FREE_MEMORY(coordInd, int, i + 1);
	  if(static_cast<int>(bary->size()) == nBasFcts) return;
662
663
664
665
666
667
668
669
	}
      }
    }
  }

  int Lagrange::getNumberOfDOFs(int dim, int degree)
  {
    int result = 0;
670
671
    for (int i = 0; i <= degree; i++) {
      result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
    }
    return result;
  }


  const DegreeOfFreedom *Lagrange::getDOFIndices(const Element *el, 
						 const DOFAdmin& admin,
						 DegreeOfFreedom *idof) const
  {
    WARNING("depricated: use getLocalIndices() instead\n");
    return getLocalIndices(el, &admin, idof);
  }

  int* Lagrange::orderOfPositionIndices(const Element* el, 
					GeoIndex position, 
					int positionIndex) const
  {
    static int sortedVertex = 0;
    static int sortedEdgeDeg2 = 0;
    static int sortedEdgeDeg3[2][2] = {{0, 1}, {1, 0}};
    static int sortedEdgeDeg4[2][3] = {{0, 1, 2}, {2, 1, 0}};
    static int sortedFaceDeg3 = 0;
    static int sortedFaceDeg4[7][3] = {{0,0,0}, {0,2,1}, {1,0,2}, {0,1,2},
				       {2,1,0}, {1,2,0}, {2,0,1}};
    static int sortedCenterDeg4 = 0;

    int dimOfPosition = DIM_OF_INDEX(position, dim);

    // vertex
701
    if (dimOfPosition == 0) {
702
703
704
705
      return &sortedVertex;
    }

    // edge
706
    if ((dimOfPosition == 1) && (degree == 2)) {
707
708
709
710
711
712
713
      return &sortedEdgeDeg2;
    }
      
    int vertex[3];
    int** dof = const_cast<int**>( el->getDOF());
    int verticesOfPosition = dimOfPosition + 1;

714
    for (int i = 0; i < verticesOfPosition; i++) {
715
716
717
718
      vertex[i] = Global::getReferenceElement(dim)->
	getVertexOfPosition(position, positionIndex, i);
    }

719
720
721
    if (dimOfPosition == 1) {
      if (degree == 3) {
	if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
722
723
724
725
726
	  return sortedEdgeDeg3[0];
	} else {
	  return sortedEdgeDeg3[1];
	}
      } else { // degree == 4
727
	if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
728
729
730
731
732
733
734
735
	  return sortedEdgeDeg4[0];
	} else {
	  return sortedEdgeDeg4[1];
	}
      }
    }

    // face
736
737
    if (dimOfPosition == 2) {
      if (degree == 3) {
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
	return &sortedFaceDeg3;
      } else {  // degree == 4!
	int no = 0;
	const Element *refElem = Global::getReferenceElement(dim);

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 0)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 1)][0]) 
	  no++;

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 1)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 2)][0])
	  no += 2;

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 2)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 0)][0])
	  no += 4;

	return sortedFaceDeg4[no];
      }
    }

    // center
760
761
    if ((dimOfPosition == 3) && (degree == 4)) {
      return &sortedCenterDeg4;
762
763
764
765
766
767
    } 

    ERROR_EXIT("should not be reached\n");
    return NULL;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
768
  void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
769
770
771
772
  {
    el_info->testFlag(Mesh::FILL_BOUND);

    // boundaries
773
    int index = 0;
774
775
    int offset = 0;
    BoundaryType boundaryType;
776
    for (int i = dim - 1; i > 0; i--)
777
      offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);
778

779
    for (int i = 0; i < dim; i++) {
780
781
      int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
      for (int j = offset; j < jto; j++) {
782
	boundaryType = el_info->getBoundary(j);
783
784
	int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
	for (int k = 0; k < kto; k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
785
	  bound[index++] = boundaryType;
786
787
	}
      }
788
      offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
789
790
791
    }

    // interior nodes in the center
792
    for (int i = 0; i < (*nDOF)[CENTER]; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
793
      bound[index++] = INTERIOR;
794
795
    }

796
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
797
798
799
800
801
802
803
804
  }


  const double* Lagrange::interpol(const ElInfo *el_info, 
				   int no, const int *b_no, 
				   AbstractFunction<double, WorldVector<double> > *f, 
				   double *vec)
  {
805
806
807
    FUNCNAME("Lagrange::interpol()");

    static double *inter = NULL;
808
809
810
811
    static int inter_size = 0;

    double *rvec =  NULL;

812
    if (vec) {
813
814
815
816
817
818
819
820
      rvec = vec;
    } else {
      if(inter) FREE_MEMORY(inter, double, nBasFcts);
      inter = GET_MEMORY(double, nBasFcts);
      inter_size = nBasFcts;
      rvec = inter;
    } 

821
    WorldVector<double> x;
822
823

    el_info->testFlag(Mesh::FILL_COORDS);
824
825
826
827
828
829
 
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	rvec[0] = 0.0;
	return(const_cast<const double *>( rvec));
830
      }
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
	
      for (int i = 0; i < no; i++) {
	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
	} else {
	  el_info->coordToWorld(*(*bary)[b_no[i]], &x);
	  rvec[i] = (*f)(x);
	}
      }
    } else {
      int vertices = Global::getGeo(VERTEX, dim);
      for (int i = 0; i < vertices; i++)
	rvec[i] = (*f)(el_info->getCoord(i));
      for (int i = vertices; i < nBasFcts; i++) {
	el_info->coordToWorld(*(*bary)[i], &x);
	rvec[i] = (*f)(x);
      }
    }  
849
850
851
852
853
854
855
856
857
858
    
    return(const_cast<const double *>( rvec));
  }

  const WorldVector<double>* 
  Lagrange::interpol(const ElInfo *el_info, int no, 
		     const int *b_no,
		     AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
		     WorldVector<double> *vec)
  {
859
860
    FUNCNAME("*Lagrange::interpol_d()");

861
    static WorldVector<double> *inter;
862
    WorldVector<double> *rvec = 
863
      vec ? vec : (inter?inter:inter= NEW WorldVector<double>[getNumber()]);
864
    WorldVector<double> x;
865
866
867
868
869

    el_info->testFlag(Mesh::FILL_COORDS);
  
    int vertices = Global::getGeo(VERTEX, dim);

870
871
872
873
874
875
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	for (int i = 0; i < dow; i++)
	  rvec[0][i] = 0.0;
	return(const_cast<const WorldVector<double> *>( rvec));
876
      }
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
      for (int i = 0; i < no; i++) {
	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
	} else {
	  el_info->coordToWorld(*(*bary)[b_no[i]], &x);
	  rvec[i] = (*f)(x);
	}
      }
    } else {
      for (int i = 0; i < vertices; i++)
	rvec[i] = (*f)(el_info->getCoord(i));
      for (int i = vertices; i < nBasFcts; i++) {
	el_info->coordToWorld(*(*bary)[i], &x);
	rvec[i] = (*f)(x);
      }
    }  
893
894
895
896
897
898
899
900
901
902
903
904
905
906
    
    return(const_cast<const WorldVector<double> *>( rvec));
  }


  const DegreeOfFreedom *Lagrange::getLocalIndices(const Element* el,
						   const DOFAdmin *admin,
						   DegreeOfFreedom *indices) const
  {
    static DegreeOfFreedom *localVec = NULL;
    static int localVecSize = 0;

    const DegreeOfFreedom **dof = el->getDOF();
    const int *indi;
907
    int nrDOFs, n0, node0, num = 0;
908
909
910
    GeoIndex posIndex;

    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
911
    
912
    if (indices) {
913
914
      result = indices;
    } else {
915
916
      if (localVec) 
	FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
917
918
919
920
      localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts);
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
921
    
922
    for (int pos = 0, j = 0; pos <= dim; pos++) {
923
924
925
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

Thomas Witkowski's avatar
Thomas Witkowski committed
926
927
928
929
930
931
      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

	for (int i = 0; i < num; node0++, i++) {
932
933
	  indi = orderOfPositionIndices(el, posIndex, i);

934
935
	  for (int k = 0; k < nrDOFs; k++)
	    result[j++] = dof[node0][n0 + indi[k]];
936
937
938
	}
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
939
    
940
941
942
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
943
944
945
946
947
948
949
950
951
952
953
  void Lagrange::getLocalIndicesVec(const Element* el,
				    const DOFAdmin *admin,
				    Vector<DegreeOfFreedom> *indices) const
  {
    if (indices->getSize() < nBasFcts) {
      indices->resize(nBasFcts);
    }

    getLocalIndices(el, admin, &((*indices)[0]));
  }

954
955
956
957
958
959
960
961
962
963
964
  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
    ERROR_EXIT("not yet\n");
  }

  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<double, WorldVector<double> >* f,
			     DOFVector<double>* fh)
  {
965
966
967
968
969
970
971
    FUNCNAME("Lagrange::l2ScpFctBas()");

    TraverseStack stack;
    double *wdetf_qp = NULL;
    double val, det;
    WorldVector<double> x;
    int iq, j;
972
973
974
    const DegreeOfFreedom  *dof;
    const Parametric *parametric;

975
976
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
    TEST_EXIT_DBG(fh->getFESpace())("no fe_space in DOF_REAL_VEC %s\n",
977
				fh->getName().c_str());
978
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)("wrong basis fcts for fh\n");
979
980

    Mesh* mesh = fh->getFESpace()->getMesh();
981
    int n_phi = nBasFcts;
982

983
984
985
    if (!q) {
      q = Quadrature::provideQuadrature(dim, 2 * degree - 2);
    }
986

987
    const FastQuadrature *quad_fast = FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
988

989
990
991
992
993
    if ((parametric = mesh->getParametric())) {
      ERROR_EXIT("not yet implemented\n");
    } else {
      wdetf_qp = GET_MEMORY(double, q->getNumPoints());
    }
994

995
    ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
996
997
998

    int numPoints = q->getNumPoints();
     
999
1000
1001
    while (el_info) {
      dof = getLocalIndices(el_info->getElement(), (fh->getFESpace()->getAdmin()), NULL);      
      det = el_info->getDet();
1002

1003
1004
1005
1006
1007
1008
1009
1010
1011
      for (iq = 0; iq < numPoints; iq++) {
	el_info->coordToWorld(q->getLambda(iq), &x);
	wdetf_qp[iq] = q->getWeight(iq)*det*((*f)(x));
      }
      
      for (j = 0; j < n_phi; j++) {
	for (val = iq = 0; iq < numPoints; iq++)
	  val += quad_fast->getPhi(iq,j)*wdetf_qp[iq];
	(*fh)[dof[j]] += val;
1012
1013
      }

1014
1015
1016
      el_info = stack.traverseNext(el_info);
    }

1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
    FREE_MEMORY(wdetf_qp, double, q->getNumPoints());

    return;
  }


  // ===== refineInter functions =====

  void  Lagrange::refineInter0(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct)
  {
1027
    FUNCNAME("Lagrange::refineInter0()");
1028

1029
1030
1031
1032
1033
1034
1035
1036
    if (n < 1) 
      return;

    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(CENTER);
    Element *el = list->getElement(0);
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);
    DegreeOfFreedom dof0 = el->getDOF(node, n0);           /* Parent center */
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0);  /*      newest vertex is center */
1037
1038
1039
1040
1041
1042
1043
1044

    (*drv)[dof_new] = (*drv)[dof0];
    dof_new = el->getChild(1)->getDOF(node, n0);  /*      newest vertex is center */
    (*drv)[dof_new] = (*drv)[dof0];
  }

  void  Lagrange::refineInter1(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct)
  {
1045
    FUNCNAME("Lagrange::refineInter1_1d()");
1046

1047
1048
    if (n < 1) 
      return;
1049

1050
1051
1052
1053
1054
1055
1056
    int dim = drv->getFESpace()->getMesh()->getDim();
    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
    Element *el = list->getElement(0);
    DegreeOfFreedom dof0 = el->getDOF(0, n0);           /* 1st endpoint of refinement edge */
    DegreeOfFreedom dof1 = el->getDOF(1, n0);           /* 2nd endpoint of refinement edge */
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0);  /*      newest vertex is DIM */
    (*drv)[dof_new] = 0.5 * ((*drv)[dof0] + (*drv)[dof1]);
1057
1058
1059
1060
  }

  void Lagrange::refineInter2_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct)
  {
1061
    FUNCNAME("Lagrange::refineInter2_1d()");
1062

1063
1064
1065
1066
1067
1068
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1069
  
1070
1071
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1072
1073
1074
1075
1076

    /****************************************************************************/
    /*  newest vertex of child[0] and child[1]                                  */
    /****************************************************************************/

1077
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+1, n0);  /*      newest vertex is DIM */
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
    (*drv)[cdof] = (*drv)[pdof[2]];

    node = drv->getFESpace()->getMesh()->getNode(CENTER);        
    n0 = admin->getNumberOfPreDOFs(CENTER);

    /****************************************************************************/
    /*  midpoint of edge on child[0] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(0)->getDOF(node, n0); 
    (*drv)[cdof] = 
1089
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1090
1091
1092
1093
1094
1095
1096

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node, n0); 
    (*drv)[cdof] = 
1097
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1098
1099
1100
1101
1102
1103
1104
    return;
  }

  void Lagrange::refineInter2_2d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
1105
    FUNCNAME("Lagrange::refineInter2_2d()");
1106

1107
1108
    if (n < 1) 
      return;
1109
1110

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1111
1112
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1113
  
1114
1115
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1116
1117
1118
1119
1120

    /****************************************************************************/
    /*  newest vertex of child[0] and child[1]                                  */
    /****************************************************************************/

1121
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
    (*drv)[cdof] = (*drv)[pdof[5]];

    node = drv->getFESpace()->getMesh()->getNode(EDGE);        
    n0 = admin->getNumberOfPreDOFs(EDGE);

    /****************************************************************************/
    /*  midpoint of edge on child[0] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(0)->getDOF(node, n0); 
    (*drv)[cdof] = 
1133
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
1134
1135
1136
1137
1138
1139
1140

    /****************************************************************************/
    /* node in the common edge of child[0] and child[1]                         */
    /****************************************************************************/

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
1141
1142
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
1143
1144
1145
1146
1147
1148
1149

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]]  + 0.75 * (*drv)[pdof[5]];

    if (n > 1) {
      /****************************************************************************/
      /*  adjust the value at the midpoint of the common edge of neigh's children */
      /****************************************************************************/
      el = list->getElement(1);
      pdof = basFct->getLocalIndices(el, admin, NULL);
      
      cdof = el->getChild(0)->getDOF(node+1, n0); 
      (*drv)[cdof] = 
	-0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
	+ 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
    }
1164
1165
1166
1167
1168
1169
  }

  void Lagrange::refineInter2_3d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
1170
    FUNCNAME("Lagrange::refineInter2_3d()");
1171

1172
1173
    if (n < 1) 
      return;
1174

1175
1176
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1177

1178
    DegreeOfFreedom pdof[10];
1179
1180
    basFct->getLocalIndices(el, admin, pdof);

1181
1182
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1183
1184
1185
1186
1187

    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/

1188
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1189
1190
1191

    (*drv)[cdof[3]] = ((*drv)[pdof[4]]);
    (*drv)[cdof[6]] = 
1192
1193
      (0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1194
    (*drv)[cdof[8]] = 
1195
1196
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
1197
    (*drv)[cdof[9]] = 
1198
1199
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
1200
1201
1202
1203
1204

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/
  
1205
    DegreeOfFreedom cdofi = el->getChild(1)->getDOF(node0 + 2, n0);
1206
    (*drv)[cdofi] = 
1207
1208
      (-0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1209
1210
1211
1212
1213

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
    for (int i = 1; i < n; i++) {
      el = list->getElement(i);
      basFct->getLocalIndices(el, admin, pdof);
      
      int lr_set = 0;
      if (list->getNeighbourElement(i, 0) && list->getNeighbourNr(i, 0) < i)
	lr_set = 1;
      
      if (list->getNeighbourElement(i, 1) && list->getNeighbourNr(i, 1) < i)
	lr_set += 2;
      
      TEST_EXIT(lr_set)("no values set on both neighbours\n");
      
      /****************************************************************************/
      /*  values on child[0]                                                      */
      /****************************************************************************/
      
      switch (lr_set) {
      case 1:
	cdofi = el->getChild(0)->getDOF(node0+4, n0);      
	(*drv)[cdofi] = 
	  (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
	   + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
	break;
      case 2:
	cdofi = el->getChild(0)->getDOF(node0+5, n0);      
	(*drv)[cdofi] = 
	  (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
	   + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
	
1244
      }
1245
    }
1246
1247
1248
1249
1250
1251
1252
    return;
  }

  void Lagrange::refineInter3_1d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
1253
    FUNCNAME("Lagrange::refineInter3_1d()");
1254

1255
1256
    if (n < 1) 
      return;
1257

1258
1259
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1260

1261
    DegreeOfFreedom *pdof = GET_MEMORY(DegreeOfFreedom, basFct->getNumber());
1262
1263
1264
1265
1266
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1267