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

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

namespace AMDiS {

19
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
51
  {
    // 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()
  {
52
53
54
55
56
57
    for (int i = 0; i < static_cast<int>(bary->size()); i++) {
      if ((*bary)[i]) {
	DELETE (*bary)[i];
	(*bary)[i] = NULL;
      }
    }
58
59
60
61
  }

  Lagrange* Lagrange::getLagrange(int dim, int degree)
  {
62
    std::list<Lagrange*>::iterator it;
63
64
    for (it = allBasFcts.begin(); it != allBasFcts.end(); it++) {
      if (((*it)->dim == dim) && ((*it)->degree == degree)) {
65
66
67
68
69
70
71
72
73
	return (*it);
      } 
    }

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

74
75
76
77
78
79
80
81
  void Lagrange::clear()
  {    
    std::list<Lagrange*>::iterator it;
    for (it = allBasFcts.begin(); it != allBasFcts.end(); it++) {
      DELETE *it;
    }
  }

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

    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()
  {
191
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
192
193
      ndofDimDegree[dim][degree] = NEW DimVec<int>(dim, DEFAULT_VALUE, 0);

194
195
196
197
198
      if (degree != 0) {
	(*ndofDimDegree[dim][degree])[VERTEX] = 1;
      } else {
	(*ndofDimDegree[dim][degree])[VERTEX] = 0;
      }
199

200
      for (int i = 1; i < dim + 1; i++) {
201
202
	nBasFcts = getNumberOfDOFs(i, degree);
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
203
	for (int j = 0; j < i; j++) {
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
	  (*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)
  {
224
    TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");
225
226
227
228
    int dimOfPosition = DIM_OF_INDEX(position, dim);

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

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

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


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

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

506
507
508
  Lagrange::GrdPhi::~GrdPhi() { 
    DELETE [] vertices; 
  }
509

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

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

633
634
635
  Lagrange::D2Phi::~D2Phi() { 
    DELETE [] vertices; 
  }
636
637
638
639
640
641
642

  void Lagrange::createCoords(int* coordInd, 
			      int numCoords, 
			      int dimIndex, 
			      int rest, 
			      DimVec<double>* vec)
  {
643
644
645
    if (vec == NULL) {
      vec = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
    }
646

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

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

    if (static_cast<int>(bary->size()) == 0) {
      for (int i = 0; i <= dim; i++) { // for all positions
665
	int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
666
667
	for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
	  int *coordInd = GET_MEMORY(int, i + 1); // indices of relevant coords
668
	  for (int k = 0; k < i + 1; k++) { 
669
670
671
	    coordInd[k] = Global::getReferenceElement(dim)->
	      getVertexOfPosition(INDEX_OF_DIM(i,dim), j, k);
	  }
672
673
674
	  createCoords(coordInd, i + 1, 0, degree);
	  FREE_MEMORY(coordInd, int, i + 1);
	  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;
683
684
    for (int i = 0; i <= degree; i++) {
      result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));
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
714
    if (dimOfPosition == 0) {
715
716
717
718
      return &sortedVertex;
    }

    // edge
719
    if ((dimOfPosition == 1) && (degree == 2)) {
720
721
722
723
      return &sortedEdgeDeg2;
    }
      
    int vertex[3];
724
    int** dof = const_cast<int**>(el->getDOF());
725
726
    int verticesOfPosition = dimOfPosition + 1;

727
    for (int i = 0; i < verticesOfPosition; i++) {
728
729
730
731
      vertex[i] = Global::getReferenceElement(dim)->
	getVertexOfPosition(position, positionIndex, i);
    }

732
733
734
    if (dimOfPosition == 1) {
      if (degree == 3) {
	if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
735
736
737
738
739
	  return sortedEdgeDeg3[0];
	} else {
	  return sortedEdgeDeg3[1];
	}
      } else { // degree == 4
740
	if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
741
742
743
744
745
746
747
748
	  return sortedEdgeDeg4[0];
	} else {
	  return sortedEdgeDeg4[1];
	}
      }
    }

    // face
749
750
    if (dimOfPosition == 2) {
      if (degree == 3) {
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
	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
773
774
    if ((dimOfPosition == 3) && (degree == 4)) {
      return &sortedCenterDeg4;
775
776
777
778
779
780
    } 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
781
  void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
782
783
784
785
  {
    el_info->testFlag(Mesh::FILL_BOUND);

    // boundaries
786
    int index = 0;
787
788
    int offset = 0;
    BoundaryType boundaryType;
789
    for (int i = dim - 1; i > 0; i--)
790
      offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);
791

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

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

809
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
810
811
812
813
814
815
816
817
  }


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

    static double *inter = NULL;
821
822
823
824
    static int inter_size = 0;

    double *rvec =  NULL;

825
    if (vec) {
826
827
828
829
830
831
832
833
      rvec = vec;
    } else {
      if(inter) FREE_MEMORY(inter, double, nBasFcts);
      inter = GET_MEMORY(double, nBasFcts);
      inter_size = nBasFcts;
      rvec = inter;
    } 

834
    WorldVector<double> x;
835
836

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

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

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

883
884
885
886
887
888
    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));
889
      }
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
      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);
      }
    }  
906
907
908
909
910
911
912
913
914
915
916
917
918
919
    
    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;
920
    int nrDOFs, n0, node0, num = 0;
921
922
923
    GeoIndex posIndex;

    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
924
    
925
    if (indices) {
926
927
      result = indices;
    } else {
928
929
      if (localVec) 
	FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
930
931
932
933
      localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts);
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
934
    
935
    for (int pos = 0, j = 0; pos <= dim; pos++) {
936
937
938
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

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

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

947
	  for (int k = 0; k < nrDOFs; k++) {
948
	    result[j++] = dof[node0][n0 + indi[k]];
949
	  }
950
951
952
	}
      }
    }
953

954
955
956
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
957
958
959
960
961
962
963
964
965
966
967
  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]));
  }

968
969
970
971
972
973
974
975
976
977
978
  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)
  {
979
980
981
982
983
984
985
    FUNCNAME("Lagrange::l2ScpFctBas()");

    TraverseStack stack;
    double *wdetf_qp = NULL;
    double val, det;
    WorldVector<double> x;
    int iq, j;
986
987
988
    const DegreeOfFreedom  *dof;
    const Parametric *parametric;

989
990
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
    TEST_EXIT_DBG(fh->getFESpace())("no fe_space in DOF_REAL_VEC %s\n",
991
				fh->getName().c_str());
992
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)("wrong basis fcts for fh\n");
993
994

    Mesh* mesh = fh->getFESpace()->getMesh();
995
    int n_phi = nBasFcts;
996

997
998
999
    if (!q) {
      q = Quadrature::provideQuadrature(dim, 2 * degree - 2);
    }
1000

1001
    const FastQuadrature *quad_fast = FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
1002

1003
1004
1005
1006
1007
    if ((parametric = mesh->getParametric())) {
      ERROR_EXIT("not yet implemented\n");
    } else {
      wdetf_qp = GET_MEMORY(double, q->getNumPoints());
    }
1008

1009
    ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
1010
1011
1012

    int numPoints = q->getNumPoints();
     
1013
1014
1015
    while (el_info) {
      dof = getLocalIndices(el_info->getElement(), (fh->getFESpace()->getAdmin()), NULL);      
      det = el_info->getDet();
1016

1017
1018
1019
1020
1021
1022
1023
1024
1025
      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;
1026
1027
      }

1028
1029
1030
      el_info = stack.traverseNext(el_info);
    }

1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
    FREE_MEMORY(wdetf_qp, double, q->getNumPoints());

    return;
  }


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

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

1043
1044
1045
1046
1047
1048
1049
1050
    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 */
1051
1052
1053
1054
1055
1056
1057
1058

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

1061
1062
    if (n < 1) 
      return;
1063

1064
1065
1066
1067
1068
1069
1070
    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]);
1071
1072
1073
1074
  }

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

1077
1078
1079
1080
1081
1082
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1083
  
1084
1085
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1086
1087
1088
1089
1090

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

1091
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+1, n0);  /*      newest vertex is DIM */
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
    (*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] = 
1103
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1104
1105
1106
1107
1108
1109
1110

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node, n0); 
    (*drv)[cdof] = 
1111
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1112
1113
1114
1115
1116
1117
1118
    return;
  }

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

1121
1122
    if (n < 1) 
      return;
1123
1124

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1125
1126
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1127
  
1128
1129
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1130
1131
1132
1133
1134

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

1135
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
    (*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] = 
1147
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
1148
1149
1150
1151
1152
1153
1154

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

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
1155
1156
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
1157
1158
1159
1160
1161
1162
1163

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
      -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]]);
    }
1178
1179
1180
1181
1182
1183
  }

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

1186
1187
    if (n < 1) 
      return;
1188

1189
1190
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1191

1192
    DegreeOfFreedom pdof[10];
1193
1194
    basFct->getLocalIndices(el, admin, pdof);

1195
1196
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1197
1198
1199
1200
1201

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

1202
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1203
1204
1205

    (*drv)[cdof[3]] = ((*drv)[pdof[4]]);
    (*drv)[cdof[6]] = 
1206
1207
      (0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1208
    (*drv)[cdof[8]] = 
1209
1210
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
1211
    (*drv)[cdof[9]] = 
1212
1213
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
1214
1215
1216
1217
1218

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/
  
1219
    DegreeOfFreedom cdofi = el->getChild(1)->getDOF(node0 + 2, n0);
1220
    (*drv)[cdofi] = 
1221
1222
      (-0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1223
1224
1225
1226
1227

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
    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]]));
	
1258
      }
1259
    }
1260
1261
1262
1263
1264
1265
1266
    return;
  }

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

1269
1270
    if (n < 1) 
      return;
1271

1272
1273
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1274

1275
    DegreeOfFreedom *pdof = GET_MEMORY(DegreeOfFreedom, basFct->getNumber());
1276
1277
1278
1279
1280
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1281
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1282
1283

    (*drv)[cdof[2]] = 
1284
1285
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1286
    (*drv)[cdof[3]] = 
1287
1288
1289
1290
      (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]];
1291
    (*drv)[cdof[6]] =  
1292
1293
1294
1295
      (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]]));
1296
    (*drv)[cdof[9]] = 
1297
1298
1299
      (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]]);
1300
1301
1302
1303
1304
1305
1306
1307
1308
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

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

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
1309
1310
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1311
    (*drv)[cdof[9]] =  
1312
1313
1314
      (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]]);
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333

    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]] = 
1334
1335
1336
1337
      (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]]));
1338
    (*drv)[cdof[9]] = 
1339
1340
1341
      (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]]);
1342
1343
1344
1345
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

1346
1347
1348
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1349
1350

    (*drv)[dof9] =  
1351
1352
1353
      (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]]);
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363

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

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

1366
1367
1368
1369
1370
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1371

1372
    DegreeOfFreedom pdof[10];
1373
1374
1375
1376
1377
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1378
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1379
1380

    (*drv)[cdof[2]] =
1381
1382
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1383
    (*drv)[cdof[3]] = 
1384
1385
1386
1387
      (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]];
1388
    (*drv)[cdof[6]] =  
1389
1390
1391
1392
      (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]]));
1393
    (*drv)[cdof[9]] = 
1394
1395
1396
      (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]]);
1397
1398
1399
1400
1401
1402
1403

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

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

1404
    (*drv)[cdof[5]] = (*drv)[pdof[8]];
1405
    (*drv)[cdof[6]] =  
1406
1407
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1408
    (*drv)[cdof[9]] =  
1409
1410
1411
      (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]]);
1412

1413
1414
    if (n <= 1)
      return;
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429

    /****************************************************************************/
    /*  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]] =  
1430
1431
1432
1433
      (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]]));
1434
    (*drv)[cdof[9]] = 
1435
1436
1437
      (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]]);
1438
1439
1440
1441
    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/

1442
1443
1444
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1445
1446

    (*drv)[dof9] =  
1447
1448
1449
      (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]]);
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
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