CFE_NormAndErrorFcts.cc 30.2 KB
Newer Older
1
#include <vector>
2
#include "CFE_NormAndErrorFcts.h"
3
4
5
6
#include "Mesh.h"
#include "Traverse.h"
#include "SubElInfo.h"

7
namespace AMDiS {
8

9
10
11
12
13
14
15
16
17
18
19
  double CFE_NormAndErrorFcts::L2_err_abs = 0.0;
  double CFE_NormAndErrorFcts::L2_u_norm = 0.0;
  double CFE_NormAndErrorFcts::H1_err_abs = 0.0;
  double CFE_NormAndErrorFcts::H1_u_norm = 0.0;

  double
  ElementL1Norm_Analyt::calcElNorm(ElInfo *elInfo, 
				   const double &det,
				   const double &fac)
  {
    double val = 0.0;
20
    WorldVector<double> worldCoordsAtQP;
21
22

    for (int iq = 0; iq < nQPts; ++iq) {
23
24
      elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
      val += q->getWeight(iq) * fabs((*f)(worldCoordsAtQP));
25
26
    }
    double nrm = det * val;
27

28
    return nrm;
29
30
  }

31
32
33
34
35
36
  double
  ElementL2Norm_Analyt::calcElNorm(ElInfo *elInfo, 
				   const double &det,
				   const double &fac)
  {
    double val = 0.0;
37
    WorldVector<double> worldCoordsAtQP;
38
39

    for (int iq = 0; iq < nQPts; ++iq) {
40
41
      elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
      val += q->getWeight(iq) * sqr((*f)(worldCoordsAtQP));
42
43
44
45
46
    }
    double nrm = det * val;

    return nrm;
  }
47

48
49
50
51
52
53
  double
  ElementH1Norm_Analyt::calcElNorm(ElInfo *elInfo, 
				   const double &det,
				   const double &fac)
  {
    double val = 0.0;
54
55
    double norm_grd2 = 0.0;
    WorldVector<double> worldCoordsAtQP;
56
57

    for (int iq = 0; iq < nQPts; ++iq) {
58
      elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
59
    
60
61
      norm_grd2 = 0.0;
      for (int j = 0; j < dim; j++)
62
	norm_grd2 += sqr(((*grd)(worldCoordsAtQP))[j]);
63
    
64
65
66
67
68
      val += q->getWeight(iq) * norm_grd2;
    }
    double nrm = det * val;

    return nrm;
69
70
  }

71
72
73
74

  double ElementL1Norm_DOF::calcElNorm(ElInfo *elInfo, 
				       const double &det,
				       const double &fac)
75
76
  {
    double val = 0.0;
77
78
    mtl::dense_vector<double> dofAtQPs(q->getNumPoints());
    dofVec->getVecAtQPs(elInfo, q, NULL, dofAtQPs); 
79

80
    for (int iq = 0; iq < nQPts; ++iq)
81
      val += q->getWeight(iq) * fabs(dofAtQPs[iq]);
82
83
    
    return det * val;
84
85
  }

86
87
88
89

  double ElementL2Norm_DOF::calcElNorm(ElInfo *elInfo, 
				       const double &det,
				       const double &fac)
90
91
  {
    double val = 0.0;
92
93
    mtl::dense_vector<double> dofAtQPs(q->getNumPoints());
    dofVec->getVecAtQPs(elInfo, q, NULL, dofAtQPs); 
94

95
    for (int iq = 0; iq < nQPts; ++iq)
96
      val += q->getWeight(iq) * sqr(dofAtQPs[iq]);
97
98
    
    return det * val;
99
100
  }

101
102
103
104
105
106
107
108
109
110
111
112
113
  double
  ElementH1Norm_DOF::calcElNorm(ElInfo *elInfo, 
				const double &det,
				const double &fac)
  {
    double val = 0.0;
    double norm_grd2;
    const WorldVector<double> *grdDofAtQPs = dofVec->getGrdAtQPs(elInfo, 
								 q, 
								 NULL, 
								 NULL);

    for (int iq = 0; iq < nQPts; ++iq) {
114
    
115
116
117
      norm_grd2 = 0.0;
      for (int j = 0; j < dim; ++j)
	norm_grd2 += sqr(grdDofAtQPs[iq][j]);
118
    
119
120
121
      val += q->getWeight(iq) * norm_grd2;
    }
    double nrm = det * val;
122

123
124
    return nrm;
  }
125

126
127
128
129

  double ElementL2Err::calcElNorm(ElInfo *elInfo, 
				  const double &det, 
				  const double &fac)
130
131
132
  {
    double val = 0.0;
    double val_nrm = 0.0;
133
134
    mtl::dense_vector<double> uhAtQPs(q->getNumPoints());
    uh->getVecAtQPs(elInfo, q, NULL, uhAtQPs); 
135
    WorldVector<double> worldCoordsAtQP;
136
137

    for (int iq = 0; iq < nQPts; ++iq) {
138
139
      elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
      val += q->getWeight(iq) * sqr((*u)(worldCoordsAtQP) - uhAtQPs[iq]);
140
    
141
      if (relErr) 
142
	val_nrm += q->getWeight(iq) * sqr((*u)(worldCoordsAtQP));
143
    }
144

145
    double err = det * val;
146

147
148
    if (relErr)
      nrmU += fac * det * val_nrm;
149

150
151
    return err;
  }
152

153
154
155
156

  double ElementH1Err::calcElNorm(ElInfo *elInfo, 
				  const double &det,
				  const double &fac)
157
158
159
160
161
162
163
164
165
  {
    double val = 0.0;
    double val_nrm = 0.0;
    double norm_err_grd2;
    double norm_grd2;
    const WorldVector<double> *grdUhAtQPs = uh->getGrdAtQPs(elInfo,
							    q,
							    NULL,
							    NULL); 
166
    WorldVector<double> worldCoordsAtQP;
167
168

    for (int iq = 0; iq < nQPts; ++iq) {
169
      elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
170
171
172
173

      norm_err_grd2 = 0.0;
      for (int j = 0; j < dim; ++j)
	norm_err_grd2 += 
174
	  sqr(((*grdu)(worldCoordsAtQP))[j] - grdUhAtQPs[iq][j]);
175
    
176
      val += q->getWeight(iq) * norm_err_grd2;
177
    
178
179
180
      if (relErr) {
	norm_grd2 = 0.0;
	for (int j = 0; j < dim; ++j)
181
	  norm_grd2 += sqr(((*grdu)(worldCoordsAtQP))[j]);
182
      
183
184
	val_nrm += q->getWeight(iq) * norm_grd2;
      }
185
    }
186
    double err = det * val;
187
  
188
189
    if (relErr) 
      nrmGrdU += fac * det * val_nrm;
190
  
191
    return err;
192
193
  }

194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
  double
  CFE_NormAndErrorFcts::Norm_IntNoBound(ElementNorm *elNorm,
					ElementLevelSet *elLS,
					Flag fillFlag,
					int deg, 
					Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::Norm_IntNoBound()");

    int dim = elLS->getDim();
    Mesh *mesh = elLS->getMesh();
    double nrm = 0.0;
    int elStatus;

    // ===== Get quadratures. =====
    if (!q) {
      q = Quadrature::provideQuadrature(dim, deg);
    }
    elNorm->setQuadrature(q);
213
214


215
216
    // ===== Traverse mesh and calculate integral on each element. =====
    TraverseStack stack;
217

218
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
219

220
    while(elInfo) {
221

222
223
      // Check whether current element is cut by the zero level set.
      elStatus = elLS->createElementLevelSet(elInfo);
224

225
      if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
226

227
228
229
	// -------------------------------------------------------------------
	//  Element is in the domain with negative level set function values.
	// -------------------------------------------------------------------
230

231
	elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
232

233
234
	nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
      }
235

236
      elInfo = stack.traverseNext(elInfo);
237

238
    }  // end of: mesh traverse
239

240
    return nrm;
241
242
  }

243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
  double
  CFE_NormAndErrorFcts::Norm_IntBound(ElementNorm *elNorm,
				      ElementLevelSet *elLS,
				      Flag fillFlag,
				      int deg, 
				      Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::Norm_IntBound()");

    int dim = elLS->getDim();
    Mesh *mesh = elLS->getMesh();
    double nrm = 0.0;
    double el_norm;
    VectorOfFixVecs<DimVec<double> > *intersecPts;
    int numIntersecPts;
    SubPolytope *subPolytope;
    ScalableQuadrature *scalQuad;
    int nScalQPts;
    int elStatus;

    // ===== Get quadratures. =====
    if (!q) {
      q = Quadrature::provideQuadrature(dim, deg);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
267
    scalQuad = new ScalableQuadrature(q);
268
    nScalQPts = scalQuad->getNumPoints();
269
270


271
272
    // ===== Traverse mesh and calculate integral on each element. =====
    TraverseStack stack;
273

274
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
275

276
277
    while(elInfo) {
      el_norm = 0.0;
278

279
280
      // Check whether current element is cut by the zero level set.
      elStatus = elLS->createElementLevelSet(elInfo);
281

282
      if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
283

284
285
286
	// -------------------------------------------------------------------
	//  Element is cut by the zero level set.
	// -------------------------------------------------------------------
287

288
289
290
	// Calculate norm on subpolyope.
	intersecPts = elLS->getElIntersecPoints();
	numIntersecPts = elLS->getNumElIntersecPoints();
291

292
293
294
295
296
297
298
299
300
301
302
303
	// -----------------------------------------------------------------
	//  Subelement may be inside the domain with negative level set
	//  function value as well as inside the domain with positive
	//  function value.
	//
	//  Whether a subelement is in the domain with negative or positive
	//  level set function values is checked by the level set function
	//  value of the first vertex of the subelement. (The subelements 
	//  are created in such a way that this vertex always is a vertex 
	//  of the element and not an intersection point. Thus the level set 
	//  function value of this vertex really is unequal to zero.)

Thomas Witkowski's avatar
Thomas Witkowski committed
304
	subPolytope = new SubPolytope(elInfo, 
305
306
307
				      intersecPts, 
				      numIntersecPts, 
				      0);
308

309
310
	elLS->setLevelSetDomain(
				elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
311

312
313
314
315
	el_norm = calcSubPolNorm(elInfo,
				 subPolytope,
				 elNorm,
				 scalQuad);
316

317
318
319
320
321
	// Calculate integral on the other element part.
	if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR) 
	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
	else 
	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
322

323
324
	elNorm->setQuadrature(q);
	el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
325

326
327
328
329
330
	el_norm -= calcSubPolNorm(elInfo,
				  subPolytope,
				  elNorm,
				  scalQuad,
				  -1.0);
331

332
	nrm += el_norm;
333
    
334
	// Free data.
Thomas Witkowski's avatar
Thomas Witkowski committed
335
	delete subPolytope;
336
337
      }
      else if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
338

339
340
341
	// -------------------------------------------------------------------
	//  Element is in the domain with negative level set function values.
	// -------------------------------------------------------------------
342

343
	elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
344

345
346
347
	elNorm->setQuadrature(q);
	nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
      }
348

349
      elInfo = stack.traverseNext(elInfo);
350

351
    }  // end of: mesh traverse
352

Thomas Witkowski's avatar
Thomas Witkowski committed
353
    delete scalQuad;
354

355
    return nrm;
356
357
  }

358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
  double
  CFE_NormAndErrorFcts::Norm_Int(ElementNorm *elNorm,
				 ElementLevelSet *elLS,
				 Flag fillFlag,
				 int deg, 
				 Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::Norm_Int()");

    int dim = elLS->getDim();
    Mesh *mesh = elLS->getMesh();
    double nrm = 0.0;
    double el_norm;
    VectorOfFixVecs<DimVec<double> > *intersecPts;
    int numIntersecPts;
    int vertex_interior;
    SubPolytope *subPolytope;
    ScalableQuadrature *scalQuad;
    int nScalQPts;
    int elStatus;

    // ===== Get quadratures. =====
    if (!q) {
      q = Quadrature::provideQuadrature(dim, deg);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
383
    scalQuad = new ScalableQuadrature(q);
384
    nScalQPts = scalQuad->getNumPoints();
385
386


387
388
    // ===== Traverse mesh and calculate integral on each element. =====
    TraverseStack stack;
389

390
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
391

392
393
    while(elInfo) {
      el_norm = 0.0;
394

395
396
      // Check whether current element is cut by the zero level set.
      elStatus = elLS->createElementLevelSet(elInfo);
397

398
      if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
399

400
401
402
	// -------------------------------------------------------------------
	//  Element is cut by the zero level set.
	// -------------------------------------------------------------------
403

404
405
406
	// Create subelements.
	intersecPts = elLS->getElIntersecPoints();
	numIntersecPts = elLS->getNumElIntersecPoints();
407

408
409
410
411
412
	if (dim == 1 || (dim == 3 && numIntersecPts == 4)) {

	  // -----------------------------------------------------------------
	  //  Subelement(s) are inside the domain with negative level set
	  //  function value.
413

414
415
416
417
418
419
	  // Get vertex with negative level set function value.
	  for (int i=0; i<=dim; ++i) {
	    if (elLS->getElVertLevelSetVec(i) < 0) {
	      vertex_interior = i;
	      break;
	    }
420
	  }
421

Thomas Witkowski's avatar
Thomas Witkowski committed
422
	  subPolytope = new SubPolytope(elInfo, 
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
					intersecPts, 
					numIntersecPts, 
					vertex_interior);

	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
	}
	else {

	  // -----------------------------------------------------------------
	  //  Subelement may be inside the domain with negative level set
	  //  function value as well as inside the domain with positive
	  //  function value.
	  //
	  //  Whether a subelement is in the domain with negative or positive
	  //  level set function values is checked by the level set function
	  //  value of the first vertex of the subelement. (The subelements 
	  //  are created in such a way that this vertex always is a vertex 
	  //  of the element and not an intersection point. Thus the level set 
	  //  function value of this vertex really is unequal to zero.)

Thomas Witkowski's avatar
Thomas Witkowski committed
443
	  subPolytope = new SubPolytope(elInfo, 
444
445
446
447
448
449
					intersecPts, 
					numIntersecPts, 
					0);

	  elLS->setLevelSetDomain(
				  elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
450
451
	}

452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
	// Calculate norm on subpolytope.
	if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_INTERIOR)
	  el_norm = calcSubPolNorm(elInfo,
				   subPolytope,
				   elNorm,
				   scalQuad);
	else 
	  el_norm = calcSubPolNorm(elInfo,
				   subPolytope,
				   elNorm,
				   scalQuad,
				   -1.0);

	// -------------------------------------------------------------------
	// In case the subelement is in the domain with positive level set
	// function values:
	// Calculate the integral on the element part with negative
	// level set function values by substracting the integral on the
	// subelement from the integral on the complete element.
	if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR) {

	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);

	  elNorm->setQuadrature(q);
	  el_norm *= -1.0;
	  el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
	}

	// Free data.
Thomas Witkowski's avatar
Thomas Witkowski committed
481
	delete subPolytope;
482
483
484
485
486
487
      }
      else if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {

	// -------------------------------------------------------------------
	//  Element is in the domain with negative level set function values.
	// -------------------------------------------------------------------
488
489

	elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
490
491
492

	elNorm->setQuadrature(q);
	el_norm = elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
493
      }
494
495
496
497
498
499

      nrm += el_norm;
      elInfo = stack.traverseNext(elInfo);

    }  // end of: mesh traverse

Thomas Witkowski's avatar
Thomas Witkowski committed
500
    delete scalQuad;
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528

    return nrm;
  }

  double
  CFE_NormAndErrorFcts::Norm_Bound(ElementNorm *elNorm,
				   ElementLevelSet *elLS,
				   Flag fillFlag,
				   int deg, 
				   Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::Norm_Bound()");

    int dim = elLS->getDim();
    Mesh *mesh = elLS->getMesh();
    double nrm = 0.0;
    double el_norm;
    VectorOfFixVecs<DimVec<double> > *intersecPts;
    int numIntersecPts;
    SubPolytope *subPolytope;
    ScalableQuadrature *scalQuad;
    int nScalQPts;
    int elStatus;

    // ===== Get quadratures. =====
    if (!q) {
      q = Quadrature::provideQuadrature(dim, deg);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
529
    scalQuad = new ScalableQuadrature(q);
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
    nScalQPts = scalQuad->getNumPoints();


    // ===== Traverse mesh and calculate integral on each element. =====
    TraverseStack stack;

    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

    while(elInfo) {
      el_norm = 0.0;

      // Check whether current element is cut by the zero level set.
      elStatus = elLS->createElementLevelSet(elInfo);

      if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {

	// -------------------------------------------------------------------
	//  Element is cut by the zero level set.
	// -------------------------------------------------------------------

	// Calculate norm on subpolyope.
	intersecPts = elLS->getElIntersecPoints();
	numIntersecPts = elLS->getNumElIntersecPoints();
553
554
555
556
557
558
559
560
561
562
563
564
565

	// -----------------------------------------------------------------
	//  Subelement may be inside the domain with negative level set
	//  function value as well as inside the domain with positive
	//  function value.
	//
	//  Whether a subelement is in the domain with negative or positive
	//  level set function values is checked by the level set function
	//  value of the first vertex of the subelement. (The subelements 
	//  are created in such a way that this vertex always is a vertex 
	//  of the element and not an intersection point. Thus the level set 
	//  function value of this vertex really is unequal to zero.)

Thomas Witkowski's avatar
Thomas Witkowski committed
566
	subPolytope = new SubPolytope(elInfo, 
567
568
569
570
571
				      intersecPts, 
				      numIntersecPts, 
				      0);

	elLS->setLevelSetDomain(
572
				elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
573
574
575
576
577
578

	el_norm = calcSubPolNorm(elInfo,
				 subPolytope,
				 elNorm,
				 scalQuad);

579
580
581
582
583
	// Calculate integral on the other element part.
	if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR) 
	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
	else 
	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
584
585
586
587

	elNorm->setQuadrature(q);
	el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));

588
589
590
591
592
	el_norm -= calcSubPolNorm(elInfo,
				  subPolytope,
				  elNorm,
				  scalQuad,
				  -1.0);
593

594
595
596
	nrm += el_norm;
    
	// Free data.
Thomas Witkowski's avatar
Thomas Witkowski committed
597
	delete subPolytope;
598
      }
599

600
      elInfo = stack.traverseNext(elInfo);
601

602
    }  // end of: mesh traverse
603

Thomas Witkowski's avatar
Thomas Witkowski committed
604
    delete scalQuad;
605

606
    return nrm;
607
608
  }

609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
  double
  CFE_NormAndErrorFcts::Norm_Complete(ElementNorm *elNorm,
				      ElementLevelSet *elLS,
				      Flag fillFlag,
				      int deg, 
				      Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::Norm_Complete()");

    int dim = elLS->getDim();
    Mesh *mesh = elLS->getMesh();
    double nrm = 0.0;
    double el_norm;
    VectorOfFixVecs<DimVec<double> > *intersecPts;
    int numIntersecPts;
    SubPolytope *subPolytope;
    ScalableQuadrature *scalQuad;
    int nScalQPts;
    int elStatus;

    // ===== Get quadratures. =====
    if (!q) {
      q = Quadrature::provideQuadrature(dim, deg);
632
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
633
    scalQuad = new ScalableQuadrature(q);
634
    nScalQPts = scalQuad->getNumPoints();
635
636


637
638
    // ===== Traverse mesh and calculate integral on each element. =====
    TraverseStack stack;
639

640
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
641

642
643
    while(elInfo) {
      el_norm = 0.0;
644

645
646
      // Check whether current element is cut by the zero level set.
      elStatus = elLS->createElementLevelSet(elInfo);
647

648
      if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
649

650
651
652
	// -------------------------------------------------------------------
	//  Element is cut by the zero level set.
	// -------------------------------------------------------------------
653

654
655
656
	// Calculate norm on subpolyope.
	intersecPts = elLS->getElIntersecPoints();
	numIntersecPts = elLS->getNumElIntersecPoints();
657

658
659
660
661
662
663
664
665
666
667
668
	// -----------------------------------------------------------------
	//  Subelement may be inside the domain with negative level set
	//  function value as well as inside the domain with positive
	//  function value.
	//
	//  Whether a subelement is in the domain with negative or positive
	//  level set function values is checked by the level set function
	//  value of the first vertex of the subelement. (The subelements 
	//  are created in such a way that this vertex always is a vertex 
	//  of the element and not an intersection point. Thus the level set 
	//  function value of this vertex really is unequal to zero.)
669

Thomas Witkowski's avatar
Thomas Witkowski committed
670
	subPolytope = new SubPolytope(elInfo, 
671
672
673
				      intersecPts, 
				      numIntersecPts, 
				      0);
674

675
676
	elLS->setLevelSetDomain(
				elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
677

678
679
680
681
	el_norm = calcSubPolNorm(elInfo,
				 subPolytope,
				 elNorm,
				 scalQuad);
682

683
684
685
686
687
	// Calculate integral on the other element part.
	if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR) 
	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
	else 
	  elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
688

689
690
	elNorm->setQuadrature(q);
	el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
691

692
693
694
695
696
	el_norm -= calcSubPolNorm(elInfo,
				  subPolytope,
				  elNorm,
				  scalQuad,
				  -1.0);
697

698
	nrm += el_norm;
699
    
700
	// Free data.
Thomas Witkowski's avatar
Thomas Witkowski committed
701
	delete subPolytope;
702
703
      }
      else {
704

705
706
707
708
	// -------------------------------------------------------------------
	//  Element is either completely in the domain with negative 
	//  or positive level set function values.
	// -------------------------------------------------------------------
709

710
	elLS->setLevelSetDomain(elStatus);
711

712
713
714
	elNorm->setQuadrature(q);
	nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
      }
715

716
      elInfo = stack.traverseNext(elInfo);
717

718
    }  // end of: mesh traverse
719

Thomas Witkowski's avatar
Thomas Witkowski committed
720
    delete scalQuad;
721

722
    return nrm;
723
724
  }

725
726
727
728
729
730
731
732
733
734
  double 
  CFE_NormAndErrorFcts::L1Norm_Analyt(
				      AbstractFunction<double, WorldVector<double> > *f, 
				      ElementLevelSet *elLS,
				      int domainFlag, 
				      int deg, 
				      Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::L1Norm_Analyt");

Thomas Witkowski's avatar
Thomas Witkowski committed
735
    ElementL1Norm_Analyt *elNorm = new ElementL1Norm_Analyt(q, f);
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
    int dim = elLS->getDim();

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET;

    double nrm = 0.0;
    switch(domainFlag) {
    case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
760

Thomas Witkowski's avatar
Thomas Witkowski committed
761
    delete elNorm;
762

763
    return nrm;  
764
765
  }

766
767
768
769
770
771
772
773
774
775
  double 
  CFE_NormAndErrorFcts::L2NormSquare_Analyt(
					    AbstractFunction<double, WorldVector<double> > *f, 
					    ElementLevelSet *elLS,
					    int domainFlag, 
					    int deg, 
					    Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::L2NormSquare_Analyt");

Thomas Witkowski's avatar
Thomas Witkowski committed
776
    ElementL2Norm_Analyt *elNorm = new ElementL2Norm_Analyt(q, f);
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
    int dim = elLS->getDim();

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET;

    double nrm = 0.0;
    switch(domainFlag) {
    case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
801

Thomas Witkowski's avatar
Thomas Witkowski committed
802
    delete elNorm;
803

804
805
    return nrm;  
  }
806

807
808
809
810
811
812
813
814
815
  double 
  CFE_NormAndErrorFcts::L2Norm_Analyt(
				      AbstractFunction<double, WorldVector<double> > *f, 
				      ElementLevelSet *elLS,
				      int domainFlag, 
				      int deg, 
				      Quadrature *q)
  {
    return sqrt(L2NormSquare_Analyt(f, elLS, domainFlag, deg, q));
816
817
  }

818
819
820
821
822
823
824
825
826
827
828
  double 
  CFE_NormAndErrorFcts::H1NormSquare_Analyt(
					    AbstractFunction<WorldVector<double>, WorldVector<double> > *grd, 
					    ElementLevelSet *elLS,
					    int domainFlag, 
					    int deg, 
					    Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::H1NormSquare_Analyt");

    int dim = elLS->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
829
    ElementH1Norm_Analyt *elNorm = new ElementH1Norm_Analyt(q, grd, dim);
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET;

    double nrm = 0.0;
    switch(domainFlag) {
    case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
853

Thomas Witkowski's avatar
Thomas Witkowski committed
854
    delete elNorm;
855

856
857
    return nrm;  
  }
858

859
860
861
862
863
864
865
866
867
  double 
  CFE_NormAndErrorFcts::H1Norm_Analyt(
				      AbstractFunction<WorldVector<double>, WorldVector<double> > *grd, 
				      ElementLevelSet *elLS,
				      int domainFlag, 
				      int deg, 
				      Quadrature *q)
  {
    return sqrt(H1NormSquare_Analyt(grd, elLS, domainFlag, deg, q));
868
869
  }

870
871
872
873
874
875
876
877
878
  double 
  CFE_NormAndErrorFcts::L1Norm_DOF(DOFVector<double> *dof, 
				   ElementLevelSet *elLS,
				   int domainFlag, 
				   int deg, 
				   Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::L1Norm_DOF");

Thomas Witkowski's avatar
Thomas Witkowski committed
879
    ElementL1Norm_DOF *elNorm = new ElementL1Norm_DOF(q, dof);
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
    int dim = elLS->getDim();

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET;

    double nrm = 0.0;
    switch(domainFlag) {
    case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
904

Thomas Witkowski's avatar
Thomas Witkowski committed
905
    delete elNorm;
906

907
    return nrm;  
908
909
  }

910
911
912
913
914
915
916
917
918
  double 
  CFE_NormAndErrorFcts::L2NormSquare_DOF(DOFVector<double> *dof, 
					 ElementLevelSet *elLS,
					 int domainFlag, 
					 int deg, 
					 Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::L2NormSquare_DOF");

Thomas Witkowski's avatar
Thomas Witkowski committed
919
    ElementL2Norm_DOF *elNorm = new ElementL2Norm_DOF(q, dof);
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
    int dim = elLS->getDim();

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET;

    double nrm = 0.0;
    switch(domainFlag) {
    case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
944

Thomas Witkowski's avatar
Thomas Witkowski committed
945
    delete elNorm;
946

947
948
    return nrm;  
  }
949

950
951
952
953
954
955
956
957
  double 
  CFE_NormAndErrorFcts::L2Norm_DOF(DOFVector<double> *dof, 
				   ElementLevelSet *elLS,
				   int domainFlag, 
				   int deg, 
				   Quadrature *q)
  {
    return sqrt(L2NormSquare_DOF(dof, elLS, domainFlag, deg, q));
958
959
  }

960
961
962
963
964
965
966
967
968
969
  double 
  CFE_NormAndErrorFcts::H1NormSquare_DOF(DOFVector<double> *dof, 
					 ElementLevelSet *elLS,
					 int domainFlag, 
					 int deg, 
					 Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::H1NormSquare_DOF");

    int dim = elLS->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
970
    ElementH1Norm_DOF *elNorm = new ElementH1Norm_DOF(q, dof, dim);
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET | 
      Mesh::FILL_GRD_LAMBDA;

    double nrm = 0.0;
    switch(domainFlag) {
    case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
995

Thomas Witkowski's avatar
Thomas Witkowski committed
996
    delete elNorm;
997

998
999
    return nrm;  
  }
1000

1001
1002
1003
1004
1005
1006
1007
1008
  double 
  CFE_NormAndErrorFcts::H1Norm_DOF(DOFVector<double> *dof, 
				   ElementLevelSet *elLS,
				   int domainFlag, 
				   int deg, 
				   Quadrature *q)
  {
    return sqrt(H1NormSquare_DOF(dof, elLS, domainFlag, deg, q));
1009
1010
  }

1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
  double 
  CFE_NormAndErrorFcts::L2Err(
			      AbstractFunction<double, WorldVector<double> > *u,
			      DOFVector<double> *uh,
			      ElementLevelSet *elLS,
			      int domainFlag,
			      int relErr,
			      int deg,
			      Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::L2Err()");

Thomas Witkowski's avatar
Thomas Witkowski committed
1023
    ElementL2Err *elNorm = new ElementL2Err(q, u, uh, relErr);
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
    int dim = elLS->getDim();

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET;

    double err = 0.0;
    switch(domainFlag) {
    case -3: err = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: err = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: err = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: err = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: err = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
1048

1049
1050
    L2_err_abs = sqrt(err);
    L2_u_norm = sqrt(elNorm->getNormU());
1051

1052
1053
1054
1055
    if (relErr)
      err = L2_err_abs / (L2_u_norm + 1.e-15);
    else 
      err = L2_err_abs;
1056

Thomas Witkowski's avatar
Thomas Witkowski committed
1057
    delete elNorm;
1058

1059
    return err;  
1060
1061
  }

1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
  double 
  CFE_NormAndErrorFcts::H1Err(
			      AbstractFunction<WorldVector<double>, WorldVector<double> > *u,
			      DOFVector<double> *uh,
			      ElementLevelSet *elLS,
			      int domainFlag,
			      int relErr,
			      int deg,
			      Quadrature *q)
  {
    FUNCNAME("CFE_NormAndErrorFcts::H1Err()");

    int dim = elLS->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
1075
    ElementH1Err *elNorm = new ElementH1Err(q, u, uh, relErr, dim);
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099

    TEST_EXIT(dim == Global::getGeo(WORLD))
      ("doesn't work for dimension of problem != dimension of world!\n");

    Flag fillFlag = Mesh::CALL_LEAF_EL | 
      Mesh::FILL_COORDS |
      Mesh::FILL_DET | 
      Mesh::FILL_GRD_LAMBDA;

    double err = 0.0;
    switch(domainFlag) {
    case -3: err = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -2: err = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
      break;
    case -1: err = Norm_Int(elNorm, elLS, fillFlag, deg, q);
      break;
    case 0: err = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
      break;
    case 1: err = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
      break;
    default: ERROR_EXIT("illegal flag !\n");
      break;
    }
1100

1101
1102
    H1_err_abs = sqrt(err);
    H1_u_norm = sqrt(elNorm->getNormGrdU());
1103

1104
1105
1106
1107
    if (relErr)
      err = H1_err_abs / (H1_u_norm + 1.e-15);
    else 
      err = H1_err_abs;
1108

Thomas Witkowski's avatar
Thomas Witkowski committed
1109
    delete elNorm;
1110

1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
    return err;  
  }

  double
  CFE_NormAndErrorFcts::calcSubPolNorm(ElInfo *elInfo,
				       SubPolytope *subPolytope,
				       ElementNorm *elNorm,
				       ScalableQuadrature *scalQuad,
				       const double &subPolFac)
  {
    double nrm = 0.0;

    for (std::vector<SubElInfo *>::iterator it = 
	   subPolytope->getSubElementsBegin(); 
	 it != subPolytope->getSubElementsEnd();
	 it++) {

      scalQuad->scaleQuadrature(**it);
      elNorm->setQuadrature(scalQuad);
1130
    
1131
1132
1133
1134
1135
1136
      nrm += elNorm->calcElNorm(elInfo, 
				fabs((*it)->getDet()),
				subPolFac);
    }

    return nrm;
1137
1138
1139
  }

}