Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

Lagrange.cc 169 KB
Newer Older
1
2
3
4
#include <stdio.h>
#include <algorithm>
#include <list>
#include "boost/lexical_cast.hpp"
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"

namespace AMDiS {

19
20
  using boost::lexical_cast;

21
22
23
24
25
26
  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];
27
  std::list<Lagrange*> Lagrange::allBasFcts;
28
29
30


  Lagrange::Lagrange(int dim, int degree)
31
    : BasisFunction(std::string("Lagrange"), dim, degree)
32
33
  {
    // set name
34
    name += lexical_cast<std::string>(dim) + " " + lexical_cast<std::string>(degree);
35
36
37
38
39
40
41
42
43
44
45
46
47

    // set nDOF
    setNDOF();

    // set barycentric coordinates
    setBary();

    // set function pointer
    setFunctionPointer();
  }

  Lagrange::~Lagrange()
  {
48
    for (int i = 0; i < static_cast<int>(bary->size()); i++)
49
      if ((*bary)[i]) {
50
	delete (*bary)[i];
51
52
	(*bary)[i] = NULL;
      }
53
54
55
56
  }

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

62
    Lagrange* newLagrange = new Lagrange(dim, degree);
63
64
65
66
    allBasFcts.push_back(newLagrange);
    return newLagrange;
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
108
    switch (degree) {
109
110
111
112
113
114
115
116
117
118
119
    case 0:
      refineInter_fct = refineInter0;
      coarseRestr_fct = coarseRestr0;
      coarseInter_fct = coarseInter0;  // not yet implemented
      break;
    case 1:
      refineInter_fct = refineInter1;
      coarseRestr_fct = coarseRestr1;
      coarseInter_fct = NULL;  // not yet implemented
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      switch (dim) {
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
      case 1:
	refineInter_fct = refineInter2_1d;
	coarseRestr_fct = NULL; // not yet implemented
	coarseInter_fct = coarseInter2_1d;
	break;
      case 2: 
	refineInter_fct = refineInter2_2d;
	coarseRestr_fct = coarseRestr2_2d;
	coarseInter_fct = coarseInter2_2d;
	break;
      case 3:
	refineInter_fct = refineInter2_3d;
	coarseRestr_fct = coarseRestr2_3d;
	coarseInter_fct = coarseInter2_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
140
      switch (dim) {
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
      case 1:
	refineInter_fct = refineInter3_1d;
	coarseRestr_fct = coarseRestr3_1d;
	coarseInter_fct = coarseInter3_1d;
	break;
      case 2: 
	refineInter_fct = refineInter3_2d;
	coarseRestr_fct = coarseRestr3_2d;
	coarseInter_fct = coarseInter3_2d;
	break;
      case 3:
	refineInter_fct = refineInter3_3d;
	coarseRestr_fct = coarseRestr3_3d;
	coarseInter_fct = coarseInter3_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 4:
Thomas Witkowski's avatar
Thomas Witkowski committed
160
      switch (dim) {
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
      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()
  {
186
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
187
      ndofDimDegree[dim][degree] = new DimVec<int>(dim, DEFAULT_VALUE, 0);
188

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

194
      for (int i = 1; i < dim + 1; i++) {
195
196
	nBasFcts = getNumberOfDOFs(i, degree);
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
197
	for (int j = 0; j < i; j++) {
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
	  (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] -= 
	    Global::getGeo(INDEX_OF_DIM(j, dim), i) * 
	    (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(j, dim)];
	}
      }
      nBasFctsDimDegree[dim][degree] = nBasFcts;
    }
    nBasFcts = nBasFctsDimDegree[dim][degree];
    nDOF = ndofDimDegree[dim][degree];
  } 

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

  void Lagrange::setVertices(int dim, int degree, 
			     GeoIndex position, int positionIndex, int nodeIndex, 
			     int** vertices)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
218
219
    FUNCNAME("Lagrange::setVertices()");

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

222
223
    int dimOfPosition = DIM_OF_INDEX(position, dim);

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

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

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

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

377
  Lagrange::Phi::~Phi() { 
378
    delete [] vertices; 
379
  }
380
381


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
  Lagrange::GrdPhi::~GrdPhi() 
  { 
506
    delete [] vertices; 
507
  }
508

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
  Lagrange::D2Phi::~D2Phi() 
  { 
634
    delete [] vertices; 
635
  }
636

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

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

  void Lagrange::setBary()
  {
657
    bary = &baryDimDegree[dim][degree];
658
659
660

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

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

683
684
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
    return result;
  }


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

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

    int dimOfPosition = DIM_OF_INDEX(position, dim);

    // vertex
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
    if (dimOfPosition == 0)
      return &sortedVertex;   
713
714

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

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

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

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

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

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

	return sortedFaceDeg4[no];
      }
    }

    // center
Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
    if ((dimOfPosition == 3) && (degree == 4))
      return &sortedCenterDeg4;    
767
768
769
770
771

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

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

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

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

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

800
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
801
802
803
804
805
806
807
808
  }


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

811
812
    static double* localVec = NULL;
    static int localVecSize = 0;
813

814
    double *rvec = NULL;
815

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

      localVecSize = nBasFcts;
      rvec = localVec;
828
829
    } 

830
    WorldVector<double> x;
831
832

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

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

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

879
880
881
882
883
884
    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));
885
      }
886
887
888
889
      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 {
890
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
891
892
893
894
895
896
897
	  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++) {
898
	el_info->coordToWorld(*(*bary)[i], x);
899
900
901
	rvec[i] = (*f)(x);
      }
    }  
902
903
904
905
906
907
908
909
910
911
912
913
914
    
    return(const_cast<const WorldVector<double> *>( rvec));
  }


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

    const DegreeOfFreedom **dof = el->getDOF();
Thomas Witkowski's avatar
Thomas Witkowski committed
915

916
    int nrDOFs, n0, node0, num = 0;
917
918
    GeoIndex posIndex;
    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
919

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

930
931
932
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
933

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

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

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

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

952
953
954
    return result;
  }

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

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

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

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

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

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

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

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

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

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

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

1016
1017
    if (!quad)
      quad = Quadrature::provideQuadrature(dim, 2 * degree - 2);
1018

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

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

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

1043
	(*fh)[dof[j]] += val;
1044
1045
      }

1046
1047
1048
      el_info = stack.traverseNext(el_info);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1049
    delete [] wdetf_qp;
1050
1051
1052
1053
1054
1055
1056
  }


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

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

1059
1060
1061
1062
1063
1064
    if (n < 1) 
      return;

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

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

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

1081
1082
    if (n < 1) 
      return;
1083

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

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

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

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

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

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

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

1145
1146
    if (n < 1) 
      return;
1147
1148

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

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

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

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

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

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

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

1210
1211
    if (n < 1) 
      return;
1212

1213
1214
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1215

1216
    DegreeOfFreedom pdof[10];
1217
1218
    basFct->getLocalIndices(el, admin, pdof);

1219
1220
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1221
1222
1223
1224
1225

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

1226
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1227
1228
1229

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

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

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

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

1292
1293
    if (n < 1) 
      return;
1294

1295
1296
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1297

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

    (*drv)[cdof[2]] = 
1307
1308
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1309
    (*drv)[cdof[3]] = 
1310
1311
1312
1313
      (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]];
1314
    (*drv)[cdof[6]] =  
1315
1316
1317
1318
      (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]]));
1319
    (*drv)[cdof[9]] = 
1320
1321
1322
      (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]]);
1323
1324
1325
1326
1327
1328
1329
1330
1331
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

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

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
1332
1333
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1334
    (*drv)[cdof[9]] =  
1335
1336
1337
      (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]]);
1338
1339

    if (n <= 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1340
      delete [] pdof;
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
      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]] = 
1357
1358
1359
1360
      (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]]));
1361
    (*drv)[cdof[9]] = 
1362
1363
1364
      (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]]);
1365
1366
1367
1368
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

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

    (*drv)[dof9] =  
1374
1375
1376
      (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]]);
1377

Thomas Witkowski's avatar
Thomas Witkowski committed
1378
    delete [] pdof;
1379
1380
1381
1382
1383
1384
1385
  }

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

1388