Lagrange.cc 179 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
62
63
64
65
66
67
68
69
70
#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 {

  ::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];
  ::std::list<Lagrange*> Lagrange::allBasFcts;


  Lagrange::Lagrange(int dim, int degree)
    : BasisFunction(::std::string("Lagrange"), dim, degree)
  {
    // 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)
  {
    ::std::list<Lagrange*>::iterator it;

    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) {
      if (degree == 4) {
762
763
764
765
766
767
768
769
770
771
772
773
	return &sortedCenterDeg4;
      }
    } 

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

  const BoundaryType* Lagrange::getBound(const ElInfo* el_info, 
					 BoundaryType* bound) const
  {
    static BoundaryType *bound_vec = NULL;
774
775
    static int bound_vec_size = 0;
    BoundaryType *result;
776

777
778
779
780
781
782
783
784
785
786
    if (bound) {
      result = bound;
    } else {
      // realloc memory if neccessary
      if (bound_vec_size < nBasFcts) {
	if (bound_vec) 
	  FREE_MEMORY(bound_vec, BoundaryType, bound_vec_size);
	bound_vec = GET_MEMORY(BoundaryType, nBasFcts);
	bound_vec_size = nBasFcts;
      }
787

788
      result = bound_vec;
789
790
791
792
793
    }

    el_info->testFlag(Mesh::FILL_BOUND);

    DimVec<int> parts(dim, NO_INIT);
794
    for (int i = 0; i < dim + 1; i++)
795
      parts[i] = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
796
797
798
799
800
801

    int index = 0;

    // boundaries
    int offset = 0;
    BoundaryType boundaryType;
802
    for (int i = dim - 1; i > 0; i--)
803
804
      offset += parts[i];

805
806
    for (int i = 0; i < dim; i++) {
      for (int j = offset; j < offset + parts[i]; j++) {
807
	boundaryType = el_info->getBoundary(j);
808
	for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
809
810
811
812
813
814
815
	  result[index++] = boundaryType;
	}
      }
      offset -= parts[i+1];
    }

    // interior nodes in the center
816
    for (int i = 0; i < (*nDOF)[CENTER]; i++) {
817
818
819
      result[index++] = INTERIOR;
    }

820
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
821
822
823
824
825
826
827
828
829
830

    return result;
  }


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

    static double *inter = NULL;
834
835
836
837
    static int inter_size = 0;

    double *rvec =  NULL;

838
    if (vec) {
839
840
841
842
843
844
845
846
      rvec = vec;
    } else {
      if(inter) FREE_MEMORY(inter, double, nBasFcts);
      inter = GET_MEMORY(double, nBasFcts);
      inter_size = nBasFcts;
      rvec = inter;
    } 

847
    WorldVector<double> x;
848
849

    el_info->testFlag(Mesh::FILL_COORDS);
850
851
852
853
854
855
 
    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));
856
      }
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
	
      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);
      }
    }  
875
876
877
878
879
880
881
882
883
884
    
    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)
  {
885
886
    FUNCNAME("*Lagrange::interpol_d()");

887
    static WorldVector<double> *inter;
888
    WorldVector<double> *rvec = 
889
      vec ? vec : (inter?inter:inter= NEW WorldVector<double>[getNumber()]);
890
    WorldVector<double> x;
891
892
893
894
895

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

896
897
898
899
900
901
    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));
902
      }
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
      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);
      }
    }  
919
920
921
922
923
924
925
926
927
928
929
930
931
932
    
    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;
933
    int nrDOFs, n0, node0, num = 0;
934
935
936
    GeoIndex posIndex;

    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
937
    
938
    if (indices) {
939
940
      result = indices;
    } else {
941
942
      if (localVec) 
	FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
943
944
945
946
      localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts);
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
947
    
948
    for (int pos = 0, j = 0; pos <= dim; pos++) {
949
950
951
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

Thomas Witkowski's avatar
Thomas Witkowski committed
952
953
954
955
956
957
      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

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

960
961
	  for (int k = 0; k < nrDOFs; k++)
	    result[j++] = dof[node0][n0 + indi[k]];
962
963
964
	}
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
965
    
966
967
968
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
969
970
971
972
973
974
975
976
977
978
979
  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]));
  }

980
981
982
983
984
985
986
987
988
989
990
  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)
  {
991
992
993
994
995
996
997
    FUNCNAME("Lagrange::l2ScpFctBas()");

    TraverseStack stack;
    double *wdetf_qp = NULL;
    double val, det;
    WorldVector<double> x;
    int iq, j;
998
999
1000
    const DegreeOfFreedom  *dof;
    const Parametric *parametric;

1001
1002
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
    TEST_EXIT_DBG(fh->getFESpace())("no fe_space in DOF_REAL_VEC %s\n",
1003
				fh->getName().c_str());
1004
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)("wrong basis fcts for fh\n");
1005
1006

    Mesh* mesh = fh->getFESpace()->getMesh();
1007
    int n_phi = nBasFcts;
1008

1009
1010
1011
    if (!q) {
      q = Quadrature::provideQuadrature(dim, 2 * degree - 2);
    }
1012

1013
    const FastQuadrature *quad_fast = FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
1014

1015
1016
1017
1018
1019
    if ((parametric = mesh->getParametric())) {
      ERROR_EXIT("not yet implemented\n");
    } else {
      wdetf_qp = GET_MEMORY(double, q->getNumPoints());
    }
1020

1021
    ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
1022
1023
1024

    int numPoints = q->getNumPoints();
     
1025
1026
1027
    while (el_info) {
      dof = getLocalIndices(el_info->getElement(), (fh->getFESpace()->getAdmin()), NULL);      
      det = el_info->getDet();
1028

1029
1030
1031
1032
1033
1034
1035
1036
1037
      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;
1038
1039
      }

1040
1041
1042
      el_info = stack.traverseNext(el_info);
    }

1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
    FREE_MEMORY(wdetf_qp, double, q->getNumPoints());

    return;
  }


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

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

1055
1056
1057
1058
1059
1060
1061
1062
    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 */
1063
1064
1065
1066
1067
1068
1069
1070

    (*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)
  {
1071
    FUNCNAME("Lagrange::refineInter1_1d()");
1072

1073
1074
    if (n < 1) 
      return;
1075

1076
1077
1078
1079
1080
1081
1082
    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]);
1083
1084
1085
1086
  }

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

1089
1090
1091
1092
1093
1094
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1095
  
1096
1097
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1098
1099
1100
1101
1102

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

1103
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+1, n0);  /*      newest vertex is DIM */
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
    (*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] = 
1115
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1116
1117
1118
1119
1120
1121
1122

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node, n0); 
    (*drv)[cdof] = 
1123
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1124
1125
1126
1127
1128
1129
1130
    return;
  }

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

1133
1134
    if (n < 1) 
      return;
1135
1136

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1137
1138
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1139
  
1140
1141
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1142
1143
1144
1145
1146

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

1147
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
    (*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] = 
1159
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
1160
1161
1162
1163
1164
1165
1166

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

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
1167
1168
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
1169
1170
1171
1172
1173
1174
1175

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
      -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]]);
    }
1190
1191
1192
1193
1194
1195
  }

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

1198
1199
    if (n < 1) 
      return;
1200

1201
1202
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1203

1204
    DegreeOfFreedom pdof[10];
1205
1206
    basFct->getLocalIndices(el, admin, pdof);

1207
1208
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1209
1210
1211
1212
1213

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

1214
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1215
1216
1217

    (*drv)[cdof[3]] = ((*drv)[pdof[4]]);
    (*drv)[cdof[6]] = 
1218
1219
      (0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1220
    (*drv)[cdof[8]] = 
1221
1222
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
1223
    (*drv)[cdof[9]] = 
1224
1225
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
1226
1227
1228
1229
1230

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/
  
1231
    DegreeOfFreedom cdofi = el->getChild(1)->getDOF(node0 + 2, n0);
1232
    (*drv)[cdofi] = 
1233
1234
      (-0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1235
1236
1237
1238
1239

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
    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]]));
	
1270
      }
1271
    }
1272
1273
1274
1275
1276
1277
1278
    return;
  }

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

1281
1282
    if (n < 1) 
      return;
1283

1284
1285
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1286

1287
    DegreeOfFreedom *pdof = GET_MEMORY(DegreeOfFreedom, basFct->getNumber());
1288
1289
1290
1291
1292
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1293
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1294
1295

    (*drv)[cdof[2]] = 
1296
1297
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1298
    (*drv)[cdof[3]] = 
1299
1300
1301
1302
      (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]];
1303
    (*drv)[cdof[6]] =  
1304
1305
1306
1307
      (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]]));
1308
    (*drv)[cdof[9]] = 
1309
1310
1311
      (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]]);
1312
1313
1314
1315
1316
1317
1318
1319
1320
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

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

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
1321
1322
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1323
    (*drv)[cdof[9]] =  
1324
1325
1326
      (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]]);
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345

    if (n <= 1) {
      FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber());
      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]] = 
1346
1347
1348
1349
      (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]]));
1350
    (*drv)[cdof[9]] = 
1351
1352
1353
      (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]]);
1354
1355
1356
1357
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

1358
1359
1360
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1361
1362

    (*drv)[dof9] =  
1363
1364
1365
      (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]]);
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375

    FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber());
    return;
  }

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

1378
1379
1380
1381
1382
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1383

1384
    DegreeOfFreedom pdof[10];
1385
1386
1387
1388
1389
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1390
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1391
1392

    (*drv)[cdof[2]] =
1393
1394
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1395
    (*drv)[cdof[3]] = 
1396
1397
1398
1399
      (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]];
1400
    (*drv)[cdof[6]] =  
1401
1402
1403
1404
      (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]]));
1405
    (*drv)[cdof[9]] = 
1406
1407
1408
      (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]]);
1409
1410
1411
1412
1413
1414
1415

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

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

1416
    (*drv)[cdof[5]] = (*drv)[pdof[8]];
1417
    (*drv)[cdof[6]] =  
1418
1419
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1420
    (*drv)[cdof[9]] =  
1421
1422
1423
      (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]]);
1424

1425
1426
    if (n <= 1)
      return;
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441

    /****************************************************************************/
    /*  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]] =  
1442
1443
1444
1445
      (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]]));
1446
    (*drv)[cdof[9]] = 
1447
1448
1449
      (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]]);
1450
1451
1452
1453
    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/

1454
1455
1456
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1457
1458

    (*drv)[dof9] =  
1459
1460
1461
      (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]]);
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586

    return;
  }

  void Lagrange::refineInter3_3d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter3_3d");
    Element *el;
    const DegreeOfFreedom *cd;
    DegreeOfFreedom pd[20], cdi;
    int       i, typ, lr_set, node0, n0;
    const DOFAdmin *admin;

    if (n < 1) return;
    el = list->getElement(0);
    typ = list->getType(0);

    admin = drv->getFESpace()->getAdmin();

    basFct->getLocalIndices(el, admin, pd);

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

    cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);

    (*drv)[cd[3]] = 
      (0.0625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) 
       + 0.5625*((*drv)[pd[4]] + (*drv)[pd[5]]));
    (*drv)[cd[8]] = 
      (0.3125*((*drv)[pd[0]] - (*drv)[pd[5]]) + 0.0625*(*drv)[pd[1]]
       + 0.937500*(*drv)[pd[4]]);
    (*drv)[cd[9]] = (*drv)[pd[4]];
    (*drv)[cd[12]] = 
      (0.0625*((*drv)[pd[0]] + (*drv)[pd[1]] - (*drv)[pd[4]] - (*drv)[pd[5]])
       + 0.25*(-(*drv)[pd[6]] - (*drv)[pd[10]])
       + 0.5*((*drv)[pd[7]] + (*drv)[pd[11]] + (*drv)[pd[19]]));
    (*drv)[cd[13]] = (*drv)[pd[19]];
    (*drv)[cd[14]] = 
      (0.0625*((*drv)[pd[0]] + (*drv)[pd[1]] - (*drv)[pd[4]] - (*drv)[pd[5]])
       + 0.25*(-(*drv)[pd[8]] - (*drv)[pd[12]]) 
       + 0.5*((*drv)[pd[9]] + (*drv)[pd[13]] + (*drv)[pd[18]]));
    (*drv)[cd[15]] = ((*drv)[pd[18]]);
    (*drv)[cd[16]] = 
      (0.0625*((*drv)[pd[0]] + (*drv)[pd[1]] - (*drv)[pd[4]] - (*drv)[pd[5]])
       + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]])
       + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) 
       + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]]));
    (*drv)[cd[17]] = 
      (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) 
       + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[8]]
       - 0.125*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
    (*drv)[cd[18]] = 
      (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) 
       + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[6]]
       - 0.125*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

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

    if (typ == 0)
      {
	/****************************************************************************/
	/*  parent of el_type 0                                                     */
	/****************************************************************************/

	(*drv)[cd[8]] = 
	  (0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
	   + 0.9375*(*drv)[pd[5]]);
	(*drv)[cd[9]] = (*drv)[pd[5]];
	(*drv)[cd[17]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) 
	   + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]])
	   - 0.125*(*drv)[pd[6]] + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
	(*drv)[cd[18]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) 
	   + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]] 
	   + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
      }
    else
      {
	/****************************************************************************/
	/*  parent of el_type 0                                                     */
	/****************************************************************************/

	(*drv)[cd[8]] = 
	  (0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
	   + 0.9375*(*drv)[pd[5]]);
	(*drv)[cd[9]] = (*drv)[pd[5]];
	(*drv)[cd[17]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) 
	   + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
	   + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
	(*drv)[cd[18]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) + 
	   0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
	   + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
      }

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/

    node0 = drv->getFESpace()->getMesh()->getNode(FACE);
    n0 = admin->getNumberOfPreDOFs(FACE);
  
    for (i = 1; i < n; i++)
      {
	el = list->getElement(i);
	typ = list->getType(i);
	basFct->getLocalIndices(el, admin, pd);

	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;

1587
	TEST_EXIT_DBG(lr_set)("no values set on both neighbours\n");
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629