Lagrange.cc 169 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
20
21
22
23
24
  std::vector<DimVec<double>* > Lagrange::baryDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  DimVec<int>* Lagrange::ndofDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  int Lagrange::nBasFctsDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  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];
25
  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
  {
    // 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()
  {
51
    for (int i = 0; i < static_cast<int>(bary->size()); i++)
52
      if ((*bary)[i]) {
53
	delete (*bary)[i];
54
55
	(*bary)[i] = NULL;
      }
56
57
58
59
  }

  Lagrange* Lagrange::getLagrange(int dim, int degree)
  {
60
    std::list<Lagrange*>::iterator it;
Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
    for (it = allBasFcts.begin(); it != allBasFcts.end(); it++)
      if (((*it)->dim == dim) && ((*it)->degree == degree))
63
64
	return (*it);

65
    Lagrange* newLagrange = new Lagrange(dim, degree);
66
67
68
69
    allBasFcts.push_back(newLagrange);
    return newLagrange;
  }

70
71
  void Lagrange::clear()
  {    
72
73
    for (std::list<Lagrange*>::iterator it = allBasFcts.begin(); 
	 it != allBasFcts.end(); it++)
74
      if (*it) {
75
	delete *it;
76
77
	*it = NULL;
      }     
78
79
  }

80
81
  void Lagrange::setFunctionPointer()
  {
82
    if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) {
83
      // for all positions
84
      for (int i = 0; i < dim+1; i++) {
85
	// no vertex dofs for degree 0 ?
86
87
	if (degree == 0 && i != dim) 
	  continue;
88
	// for all vertices/edges/...
89
	for (int j = 0; j < Global::getGeo(INDEX_OF_DIM(i, dim),dim); j++) {
90
	  // for all dofs at this position
91
	  for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
92
	    // basis function
93
	    phiDimDegree[dim][degree].push_back(new Phi(this, 
94
							INDEX_OF_DIM(i ,dim),j, k));
95
	    // gradients
96
	    grdPhiDimDegree[dim][degree].push_back(new GrdPhi(this, 
97
98
							      INDEX_OF_DIM(i, dim),
							      j, k));
99
	    // D2-Matrices
100
	    D2PhiDimDegree[dim][degree].push_back(new D2Phi(this, 
101
102
							    INDEX_OF_DIM(i, dim),
							    j, k));
103
104
105
106
	  }
	}
      }
    }
107
    phi = &phiDimDegree[dim][degree];
108
    grdPhi = &grdPhiDimDegree[dim][degree];
109
    d2Phi = &D2PhiDimDegree[dim][degree];
110

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    switch (degree) {
112
113
114
115
116
117
118
119
120
121
122
    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:
Thomas Witkowski's avatar
Thomas Witkowski committed
123
      switch (dim) {
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
      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:
Thomas Witkowski's avatar
Thomas Witkowski committed
143
      switch (dim) {
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
      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:
Thomas Witkowski's avatar
Thomas Witkowski committed
163
      switch (dim) {
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
      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()
  {
189
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
190
      ndofDimDegree[dim][degree] = new DimVec<int>(dim, DEFAULT_VALUE, 0);
191

Thomas Witkowski's avatar
Thomas Witkowski committed
192
      if (degree != 0)
193
	(*ndofDimDegree[dim][degree])[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
      else
	(*ndofDimDegree[dim][degree])[VERTEX] = 0;      
196

197
      for (int i = 1; i < dim + 1; i++) {
198
199
	nBasFcts = getNumberOfDOFs(i, degree);
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
200
	for (int j = 0; j < i; j++) {
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
	  (*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)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
221
222
    FUNCNAME("Lagrange::setVertices()");

223
    TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
224

225
226
    int dimOfPosition = DIM_OF_INDEX(position, dim);

Thomas Witkowski's avatar
Thomas Witkowski committed
227
    *vertices = new int[dimOfPosition + 1];
228

229
    if ((degree == 4) && (dimOfPosition==1)) {
230
231
232
233
234
235
236
237
238
      (*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)) {
239
240
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+dimOfPosition*nodeIndex) % (dimOfPosition + 1)] = 
241
242
243
244
245
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }    
    } else {
246
247
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+nodeIndex) % (dimOfPosition + 1)] = 
248
249
250
251
252
253
254
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }
    }
  }

255
  Lagrange::Phi::Phi(Lagrange* owner, 
256
257
258
		     GeoIndex position, 
		     int positionIndex, 
		     int nodeIndex)
259
    : vertices(NULL)
260
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
264
265
266
267
268
269
    FUNCNAME("Lagrange::Phi::Phi()");

    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);
270
271
  
    // set function pointer
272
    switch (owner->getDegree()) {
273
    case 0:
274
      switch (position) {
275
276
277
278
279
280
281
282
      case CENTER:
	func = phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
283
      switch (position) {
284
285
286
287
288
289
290
291
      case VERTEX:
	func = phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
292
      switch (position) {
293
294
295
296
      case VERTEX:
	func = phi2v;
	break;
      case EDGE:
297
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
298
299
300
	func = phi2e;
	break;
      case CENTER:
301
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
302
303
304
305
306
307
308
	func = phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
309
      switch (position) {
310
311
312
313
314
315
316
      case VERTEX:
	func = phi3v;
	break;
      case EDGE:
	func = phi3e;
	break;
      case FACE:
317
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
318
319
320
	func = phi3f;
	break;
      case CENTER:
321
	switch (owner->getDim()) {
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
	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:
338
      switch (position) {
339
340
341
342
      case VERTEX:
	func = phi4v;
	break;
      case EDGE:
343
	if (nodeIndex == 1)
344
345
346
347
348
	  func = phi4e1;
	else 
	  func = phi4e0;
	break;
      case FACE:
349
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
350
351
352
	func = phi4f;      
	break;
      case CENTER:
353
	switch (owner->getDim()) {
354
	case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
355
	  if (nodeIndex == 1)
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
	    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");
    }
  }

380
  Lagrange::Phi::~Phi() { 
381
    delete [] vertices; 
382
  }
383
384


385
  Lagrange::GrdPhi::GrdPhi(Lagrange* owner, 
386
387
388
			   GeoIndex position, 
			   int positionIndex, 
			   int nodeIndex)
389
    : vertices(NULL)
390
391
392
393
394
395
396
397
398
399
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
  Lagrange::GrdPhi::~GrdPhi() 
  { 
509
    delete [] vertices; 
510
  }
511

512
  Lagrange::D2Phi::D2Phi(Lagrange* owner, 
513
514
515
			 GeoIndex position, 
			 int positionIndex, 
			 int nodeIndex)
516
    : vertices(NULL)
517
518
519
520
521
522
523
524
525
526
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
527
    switch (owner->getDegree()) {
528
    case 0:
529
      switch (position) {
530
531
532
533
534
535
536
537
      case CENTER:
	func = D2Phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
538
      switch (position) {
539
540
541
542
543
544
545
546
      case VERTEX:
	func = D2Phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
547
      switch (position) {
548
549
550
551
      case VERTEX:
	func = D2Phi2v;
	break;
      case EDGE:
552
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
553
554
555
	func = D2Phi2e;
	break;
      case CENTER:
556
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
557
558
559
560
561
562
563
	func = D2Phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
564
      switch (position) {
565
566
567
568
569
570
571
      case VERTEX:
	func = D2Phi3v;
	break;
      case EDGE:
	func = D2Phi3e;
	break;
      case FACE:
572
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
573
574
575
	func = D2Phi3f;
	break;
      case CENTER:
576
	switch (owner->getDim()) {
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
	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:
593
      switch (position) {
594
595
596
597
      case VERTEX:
	func = D2Phi4v;
	break;
      case EDGE:
Thomas Witkowski's avatar
Thomas Witkowski committed
598
	if (nodeIndex == 1)
599
600
601
602
603
	  func = D2Phi4e1;
	else 
	  func = D2Phi4e0;
	break;
      case FACE:
604
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
605
606
607
	func = D2Phi4f;      
	break;
      case CENTER:
608
	switch (owner->getDim()) {
609
	case 1:
610
	  if (nodeIndex == 1)
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
	    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");
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
  Lagrange::D2Phi::~D2Phi() 
  { 
637
    delete [] vertices; 
638
  }
639

Thomas Witkowski's avatar
Thomas Witkowski committed
640
  void Lagrange::createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
641
642
			      DimVec<double>* vec)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
643
    if (vec == NULL)
644
      vec = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
645

646
647
    if (dimIndex == numCoords - 1) {
      (*vec)[coordInd[dimIndex]] = double(rest) / degree;
648
      DimVec<double>* newCoords = new DimVec<double>(*vec);
649
650
      bary->push_back(newCoords);
    } else {
651
      for (int i = rest - 1; i >= 1; i--) {
652
	(*vec)[coordInd[dimIndex]] = double(i) / degree;
653
	createCoords(coordInd, numCoords, dimIndex + 1, rest - i, vec);
654
655
656
657
658
659
      }
    }
  }

  void Lagrange::setBary()
  {
660
    bary = &baryDimDegree[dim][degree];
661
662
663

    if (static_cast<int>(bary->size()) == 0) {
      for (int i = 0; i <= dim; i++) { // for all positions
664
	int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
665
	for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
Thomas Witkowski's avatar
Thomas Witkowski committed
666
	  int *coordInd = new int[i + 1];      // indices of relevant coords
667
	  for (int k = 0; k < i + 1; k++) { 
668
	    coordInd[k] = Global::getReferenceElement(dim)->
Thomas Witkowski's avatar
Thomas Witkowski committed
669
	      getVertexOfPosition(INDEX_OF_DIM(i, dim), j, k);
670
	  }
671
	  createCoords(coordInd, i + 1, 0, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
672
673
674
	  delete [] coordInd;
	  if (static_cast<int>(bary->size()) == nBasFcts) 
	    return;
675
676
677
678
679
680
681
682
	}
      }
    }
  }

  int Lagrange::getNumberOfDOFs(int dim, int degree)
  {
    int result = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
683
    for (int i = 0; i <= degree; i++)
684
      result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));
Thomas Witkowski's avatar
Thomas Witkowski committed
685

686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
    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
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
    if (dimOfPosition == 0)
      return &sortedVertex;   
716
717

    // edge
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
    if ((dimOfPosition == 1) && (degree == 2))
      return &sortedEdgeDeg2;    
720
721
      
    int vertex[3];
722
    int** dof = const_cast<int**>(el->getDOF());
723
724
    int verticesOfPosition = dimOfPosition + 1;

Thomas Witkowski's avatar
Thomas Witkowski committed
725
    for (int i = 0; i < verticesOfPosition; i++)
726
      vertex[i] = Global::getReferenceElement(dim)->
Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
	getVertexOfPosition(position, positionIndex, i);   
    
729
730
    if (dimOfPosition == 1) {
      if (degree == 3) {
Thomas Witkowski's avatar
Thomas Witkowski committed
731
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
732
	  return sortedEdgeDeg3[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
733
734
	else
	  return sortedEdgeDeg3[1];	
735
      } else { // degree == 4
Thomas Witkowski's avatar
Thomas Witkowski committed
736
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
737
	  return sortedEdgeDeg4[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
738
739
	else
	  return sortedEdgeDeg4[1];	
740
741
742
743
      }
    }

    // face
744
745
    if (dimOfPosition == 2) {
      if (degree == 3) {
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
	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
Thomas Witkowski's avatar
Thomas Witkowski committed
768
769
    if ((dimOfPosition == 3) && (degree == 4))
      return &sortedCenterDeg4;    
770
771
772
773
774

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

Thomas Witkowski's avatar
Thomas Witkowski committed
775
  void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
776
777
778
779
  {
    el_info->testFlag(Mesh::FILL_BOUND);

    // boundaries
780
    int index = 0;
781
782
    int offset = 0;
    BoundaryType boundaryType;
783
    for (int i = dim - 1; i > 0; i--)
784
      offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);
785

786
    for (int i = 0; i < dim; i++) {
787
788
      int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
      for (int j = offset; j < jto; j++) {
789
	boundaryType = el_info->getBoundary(j);
790
791
	int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
	for (int k = 0; k < kto; k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
792
	  bound[index++] = boundaryType;
793
794
	}
      }
795
      offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
796
797
798
    }

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

803
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
804
805
806
807
808
809
810
811
  }


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

814
815
    static double* localVec = NULL;
    static int localVecSize = 0;
816

817
    double *rvec = NULL;
818

819
    if (vec) {
820
821
      rvec = vec;
    } else {
822
823
824
825
826
827
828
829
830
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new double[nBasFcts]; 
      }
      if (!localVec)
	localVec = new double[nBasFcts]; 

      localVecSize = nBasFcts;
      rvec = localVec;
831
832
    } 

833
    WorldVector<double> x;
834
835

    el_info->testFlag(Mesh::FILL_COORDS);
836
837
838
839
840
841
 
    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));
842
      }
843
844
845
846
847
	
      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 {
848
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
849
850
851
852
853
854
855
856
	  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++) {
857
	el_info->coordToWorld(*(*bary)[i], x);
858
859
860
	rvec[i] = (*f)(x);
      }
    }  
861
862
863
864
865
866
867
868
869
870
    
    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)
  {
871
872
    FUNCNAME("*Lagrange::interpol_d()");

873
    static WorldVector<double> *inter;
874
    WorldVector<double> *rvec = 
875
      vec ? vec : (inter ? inter : inter = new WorldVector<double>[getNumber()]);
876
    WorldVector<double> x;
877
878
879
880
881

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

882
883
884
885
886
887
    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));
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 {
893
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
894
895
896
897
898
899
900
	  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++) {
901
	el_info->coordToWorld(*(*bary)[i], x);
902
903
904
	rvec[i] = (*f)(x);
      }
    }  
905
906
907
908
909
910
911
912
913
914
915
916
917
    
    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();
Thomas Witkowski's avatar
Thomas Witkowski committed
918

919
    int nrDOFs, n0, node0, num = 0;
920
921
    GeoIndex posIndex;
    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
922

923
    if (indices) {
924
925
      result = indices;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
926
      if (localVec && nBasFcts > localVecSize) {
Thomas Witkowski's avatar
Thomas Witkowski committed
927
	delete [] localVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
928
929
930
931
932
	localVec = new DegreeOfFreedom[nBasFcts];
      }
      if (!localVec)
	localVec = new DegreeOfFreedom[nBasFcts];

933
934
935
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
936

937
    for (int pos = 0, j = 0; pos <= dim; pos++) {
938
939
940
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

Thomas Witkowski's avatar
Thomas Witkowski committed
941
942
943
944
945
946
      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

	for (int i = 0; i < num; node0++, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
947
	  const int *indi = orderOfPositionIndices(el, posIndex, i);
948

Thomas Witkowski's avatar
Thomas Witkowski committed
949
	  for (int k = 0; k < nrDOFs; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
950
	    result[j++] = dof[node0][n0 + indi[k]];  
951
952
953
	}
      }
    }
954

955
956
957
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
958
959
960
961
  void Lagrange::getLocalIndicesVec(const Element* el,
				    const DOFAdmin *admin,
				    Vector<DegreeOfFreedom> *indices) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
962
    if (indices->getSize() < nBasFcts)
Thomas Witkowski's avatar
Thomas Witkowski committed
963
964
965
966
967
      indices->resize(nBasFcts);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
  void Lagrange::getLocalDofPtrVec(const Element *el, 
				   const DOFAdmin *admin,
				   std::vector<const DegreeOfFreedom*>& vec) const
  {
    if (static_cast<int>(vec.size()) < nBasFcts)
      vec.resize(nBasFcts);

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

    for (int pos = 0, j = 0; pos <= dim; pos++) {
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

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

	  for (int k = 0; k < nrDOFs; k++)
	    vec[j++] = &(dof[node0][n0 + indi[k]]);
	}
      }
    }
  }

998
999
1000
1001
1002
1003
1004
  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
    ERROR_EXIT("not yet\n");
  }

1005
  void Lagrange::l2ScpFctBas(Quadrature *quad,
1006
1007
1008
			     AbstractFunction<double, WorldVector<double> >* f,
			     DOFVector<double>* fh)
  {
1009
1010
    FUNCNAME("Lagrange::l2ScpFctBas()");

1011
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1012
1013
1014
1015
    TEST_EXIT_DBG(fh->getFESpace())
      ("no fe_space in DOF_REAL_VEC %s\n", fh->getName().c_str());
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)
      ("wrong basis fcts for fh\n");
1016
1017
    TEST_EXIT_DBG(!fh->getFESpace()->getMesh()->getParametric())
      ("Not yet implemented!");
1018

1019
1020
    if (!quad)
      quad = Quadrature::provideQuadrature(dim, 2 * degree - 2);
1021

Thomas Witkowski's avatar
Thomas Witkowski committed
1022
    const FastQuadrature *quad_fast = 
1023
1024
1025
1026
1027
      FastQuadrature::provideFastQuadrature(this, *quad, INIT_PHI);
    double *wdetf_qp = new double[quad->getNumPoints()];
    int nPoints = quad->getNumPoints();
    DOFAdmin *admin = fh->getFESpace()->getAdmin();
    WorldVector<double> x;
1028

1029
1030
    TraverseStack stack;
    ElInfo *el_info = stack.traverseFirst(fh->getFESpace()->getMesh(), -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
1031
					  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
1032
    while (el_info) {
1033
1034
      const DegreeOfFreedom  *dof = getLocalIndices(el_info->getElement(), admin, NULL);
      double det = el_info->getDet();
1035

1036
1037
1038
      for (int iq = 0; iq < nPoints; iq++) {
	el_info->coordToWorld(quad->getLambda(iq), x);
	wdetf_qp[iq] = quad->getWeight(iq) * det * ((*f)(x));
1039
1040
      }
      
1041
1042
1043
1044
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
	for (int iq = 0; iq < nPoints; iq++)
	  val += quad_fast->getPhi(iq, j) * wdetf_qp[iq];
Thomas Witkowski's avatar
Thomas Witkowski committed
1045

1046
	(*fh)[dof[j]] += val;
1047
1048
      }

1049
1050
1051
      el_info = stack.traverseNext(el_info);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1052
    delete [] wdetf_qp;
1053
1054
1055
1056
1057
1058
1059
  }


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

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

1062
1063
1064
1065
1066
1067
    if (n < 1) 
      return;

    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(CENTER);
    Element *el = list->getElement(0);
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);
1068
1069
1070
1071
    // Parent center
    DegreeOfFreedom dof0 = el->getDOF(node, n0);           
    // Newest vertex is center 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0);  
1072
1073

    (*drv)[dof_new] = (*drv)[dof0];
1074
1075
    // Newest vertex is center 
    dof_new = el->getChild(1)->getDOF(node, n0);  
1076
1077
1078
    (*drv)[dof_new] = (*drv)[dof0];
  }

1079
1080
  void Lagrange::refineInter1(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
			      BasisFunction* basFct)
1081
  {
1082
    FUNCNAME("Lagrange::refineInter1_1d()");
1083

1084
1085
    if (n < 1) 
      return;
1086

1087
1088
1089
    int dim = drv->getFESpace()->getMesh()->getDim();
    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
    Element *el = list->getElement(0);
1090
1091
1092
1093
1094
1095
    // 1st endpoint of refinement edge
    DegreeOfFreedom dof0 = el->getDOF(0, n0);
    // 2nd endpoint of refinement edge          
    DegreeOfFreedom dof1 = el->getDOF(1, n0);
    // newest vertex is DIM 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0); 
1096
    (*drv)[dof_new] = 0.5 * ((*drv)[dof0] + (*drv)[dof1]);
1097
1098
  }

1099
1100
  void Lagrange::refineInter2_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
				 BasisFunction* basFct)
1101
  {
1102
    FUNCNAME("Lagrange::refineInter2_1d()");
1103

1104
1105
1106
1107
1108
1109
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1110
  
1111
1112
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1113
1114
1115
1116
1117

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

1118
1119
    // newest vertex is DIM
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node + 1, n0);
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
    (*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] = 
1131
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1132
1133
1134
1135
1136
1137
1138

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node, n0); 
    (*drv)[cdof] = 
1139
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1140
1141
1142
1143
1144
1145
  }

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

1148
1149
    if (n < 1) 
      return;
1150
1151

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1152
1153
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1154
  
1155
1156
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1157
1158
1159
1160
1161

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

1162
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
    (*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] = 
1174
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
1175
1176
1177
1178
1179
1180
1181

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

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
1182
1183
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
1184
1185
1186
1187
1188
1189
1190

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
      -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]]);
    }
1205
1206
1207
1208
1209
1210
  }

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

1213
1214
    if (n < 1) 
      return;
1215

1216
1217
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1218

1219
    DegreeOfFreedom pdof[10];
1220
1221
    basFct->getLocalIndices(el, admin, pdof);

1222
1223
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1224
1225
1226
1227
1228

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

1229
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1230
1231
1232

    (*drv)[cdof[3]] = ((*drv)[pdof[4]]);
    (*drv)[cdof[6]] = 
1233
1234
      (0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1235
    (*drv)[cdof[8]] = 
1236
1237
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
1238
    (*drv)[cdof[9]] = 
1239
1240
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
1241
1242
1243
1244
1245

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/
  
1246
    DegreeOfFreedom cdofi = el->getChild(1)->getDOF(node0 + 2, n0);
1247
    (*drv)[cdofi] = 
1248
1249
      (-0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1250
1251
1252
1253
1254

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
    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]]));
	
1285
      }
1286
    }
1287
1288
1289
1290
1291
1292
  }

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

1295
1296
    if (n < 1) 
      return;
1297

1298
1299
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1300

Thomas Witkowski's avatar
Thomas Witkowski committed
1301
    DegreeOfFreedom *pdof = new DegreeOfFreedom[basFct->getNumber()];
1302
1303
1304
1305
1306
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1307
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1308
1309

    (*drv)[cdof[2]] = 
1310
1311
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1312
    (*drv)[cdof[3]] = 
1313
1314
1315
1316
      (0.3125 * ((*drv)[pdof[0]] - (*drv)[pdof[8]]) + 0.0625 * (*drv)[pdof[1]]
       + 0.9375 * (*drv)[pdof[7]]);
    (*drv)[cdof[4]] = (*drv)[pdof[7]];
    (*drv)[cdof[5]] = (*drv)[pdof[9]];
1317
    (*drv)[cdof[6]] =  
1318
1319
1320
1321
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1322
    (*drv)[cdof[9]] = 
1323
1324
1325
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.125 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1326
1327
1328
1329
1330
1331
1332
1333
1334
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

    cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
1335
1336
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1337
    (*drv)[cdof[9]] =  
1338
1339
1340
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1341
1342

    if (n <= 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1343
      delete [] pdof;
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
      return;
    }
    /****************************************************************************/
    /*  adjust the value on the neihgbour                                       */
    /****************************************************************************/

    el = list->getElement(1);
    basFct->getLocalIndices(el, admin, pdof);

    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/
    cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
  
    (*drv)[cdof[5]] =  (*drv)[pdof[9]];
    (*drv)[cdof[6]] = 
1360
1361
1362
1363
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1364
    (*drv)[cdof[9]] = 
1365
1366
1367
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.12500 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1368
1369
1370
1371
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

1372
1373
1374
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1375
1376

    (*drv)[dof9] =  
1377
1378
1379
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) +  0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75  *(*drv)[pdof[9]]);
1380

Thomas Witkowski's avatar
Thomas Witkowski committed
1381
    delete [] pdof;
1382
1383
1384
1385
1386
1387
1388
  }

  void Lagrange::refineInter3_2d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, 
				 BasisFunction* basFct)
  {
1389
    FUNCNAME("Lagrange::refineInter3_2d()");
1390

1391
1392
1393
1394
1395
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1396

1397
    DegreeOfFreedom pdof[10];
1398
1399
1400
1401
1402
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1403
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1404
1405

    (*drv)[cdof[2]] =
1406
1407
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1408
    (*drv)[cdof[3]] = 
1409
1410
1411
1412
      (0.3125 * ((*drv)[pdof[0]] - (*drv)[pdof[8]]) + 0.0625 * (*drv)[pdof[1]]
       + 0.9375 * (*drv)[pdof[7]]);
    (*drv)[cdof[4]] = (*drv)[pdof[7]];
    (*drv)[cdof[5]] = (*drv)[pdof[9]];
1413
    (*drv)[cdof[6]] =  
1414
1415
1416
1417
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1418
    (*drv)[cdof[9]] = 
1419
1420
1421
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.125 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1422
1423
1424
1425
1426
1427
1428

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

    cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);

1429
    (*drv)[cdof[5]] = (*drv)[pdof[8]];
1430
    (*drv)[cdof[6]] =  
1431
1432
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));