ResidualEstimator.cc 19.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#include "ResidualEstimator.h"
#include "Operator.h"
#include "DOFMatrix.h"
#include "DOFVector.h"
#include "Assembler.h"
#include "Traverse.h"
#include "Parameters.h"

namespace AMDiS {

11
  ResidualEstimator::ResidualEstimator(std::string name, int r) 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
    : Estimator(name, r),
      C0(1.0), 
      C1(1.0), 
      C2(1.0), 
      C3(1.0)
  {
    GET_PARAMETER(0, name + "->C0", "%f", &C0);
    GET_PARAMETER(0, name + "->C1", "%f", &C1);
    GET_PARAMETER(0, name + "->C2", "%f", &C2);
    GET_PARAMETER(0, name + "->C3", "%f", &C3);

    C0 = C0 > 1.e-25 ? sqr(C0) : 0.0;
    C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
    C2 = C2 > 1.e-25 ? sqr(C2) : 0.0;
    C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;
  }

  void ResidualEstimator::init(double ts)
  {
    FUNCNAME("ResidualEstimator::init()");
Thomas Witkowski's avatar
Thomas Witkowski committed
32
    
33
34
35
36
    timestep = ts;

    mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
37
38
    nSystems = static_cast<int>(uh.size());
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
39
40

    dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
    basFcts = GET_MEMORY(const BasisFunction*, nSystems);
    quadFast = GET_MEMORY(FastQuadrature*, nSystems);
43
44

    degree = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
45
    for (int system = 0; system < nSystems; system++) {
46
      basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
47
      degree = std::max(degree, basFcts[system]->getDegree());
48
49
50
51
52
    }

    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
53
    nPoints = quad->getNumPoints();
54
55
56
57
58
59

    Flag flag = INIT_PHI | INIT_GRD_PHI;
    if (degree > 2) {
      flag |= INIT_D2_PHI;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
60
    for (int system = 0; system < nSystems; system++) {
61
62
63
64
65
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
							       *quad, 
							       flag);
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
    uhEl = GET_MEMORY(double*, nSystems);
    uhNeigh = GET_MEMORY(double*, nSystems);
    uhOldEl = timestep ? GET_MEMORY(double*, nSystems) : NULL;
69

Thomas Witkowski's avatar
Thomas Witkowski committed
70
    for (int system = 0; system < nSystems; system++) {
71
      uhEl[system] = GET_MEMORY(double, basFcts[system]->getNumber()); 
72
      uhNeigh[system] = GET_MEMORY(double, basFcts[system]->getNumber());
73
74
75
76
      if (timestep)
	uhOldEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
    uhQP = timestep ? GET_MEMORY(double, nPoints) : NULL;
    uhOldQP = timestep ? GET_MEMORY(double, nPoints) : NULL;
79

Thomas Witkowski's avatar
Thomas Witkowski committed
80
    riq = GET_MEMORY(double, nPoints);
81

82
83
84
    grdUh_qp = NULL;
    D2uhqp = NULL;

85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    TraverseStack stack;
    ElInfo *elInfo = NULL;

    // clear error indicators and mark elements for jumpRes
    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      elInfo->getElement()->setEstimation(0.0, row);
      elInfo->getElement()->setMark(1);
      elInfo = stack.traverseNext(elInfo);
    }

    est_sum = 0.0;
    est_max = 0.0;
    est_t_sum = 0.0;
    est_t_max = 0.0;

    traverseFlag = 
      Mesh::FILL_NEIGH      |
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND      |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_DET        |
      Mesh::CALL_LEAF_EL;
109
110
111
112
113
114
115
116
117
118
119
120

    neighInfo = mesh->createNewElInfo();

    // prepare date for computing jump residual
    if (C1 && (dim > 1)) {
      surfaceQuad_ = Quadrature::provideQuadrature(dim - 1, degree);
      nPointsSurface_ = surfaceQuad_->getNumPoints();
      grdUhEl_.resize(nPointsSurface_);
      grdUhNeigh_.resize(nPointsSurface_);
      jump_.resize(nPointsSurface_);
      localJump_.resize(nPointsSurface_);
      neighbours_ = Global::getGeo(NEIGH, dim);
121
122
      lambdaNeigh_ = NEW DimVec<WorldVector<double> >(dim, NO_INIT);
      lambda_ = NEW DimVec<double>(dim, NO_INIT);
123
    }
124
125
126
127
128
129
130
131
132
  }

  void ResidualEstimator::exit(bool output)
  {
    FUNCNAME("ResidualEstimator::exit()");

    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

Thomas Witkowski's avatar
Thomas Witkowski committed
133
    for (int system = 0; system < nSystems; system++) {
134
      FREE_MEMORY(uhEl[system], double, basFcts[system]->getNumber());
135
      FREE_MEMORY(uhNeigh[system], double, basFcts[system]->getNumber());
136
137
138
139
      if (timestep)
	FREE_MEMORY(uhOldEl[system], double, basFcts[system]->getNumber());    
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
    FREE_MEMORY(uhEl, double*, nSystems);
    FREE_MEMORY(uhNeigh, double*, nSystems);
142
143

    if (timestep) {
Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
146
      FREE_MEMORY(uhOldEl, double*, nSystems);
      FREE_MEMORY(uhQP, double, nPoints);
      FREE_MEMORY(uhOldQP, double, nPoints);
147
148
    } else {
      if (uhQP != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
149
	FREE_MEMORY(uhQP, double, nPoints);
150
      }
151
152
153
154
155
156
157
158
159
    }

    if (output) {
      MSG("estimate   = %.8e\n", est_sum);
      if (C3) {
	MSG("time estimate   = %.8e\n", est_t_sum);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
160
161
162
    FREE_MEMORY(riq, double, nPoints);
    FREE_MEMORY(basFcts, const BasisFunction*, nSystems);
    FREE_MEMORY(quadFast, FastQuadrature*, nSystems);
163
164

    if (grdUh_qp != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
165
      FREE_MEMORY(grdUh_qp, WorldVector<double>, nPoints);
166
167
    }
    if (D2uhqp != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
      FREE_MEMORY(D2uhqp, WorldMatrix<double>, nPoints);
169
    }
170

171
172
173
174
175
    if (C1 && (dim > 1)) {
      DELETE lambdaNeigh_;
      DELETE lambda_;
    }

176
    DELETE neighInfo;
177
178
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
179
#if 1
180
  void ResidualEstimator::estimateElement(ElInfo *elInfo)
181
  {    
182
183
    FUNCNAME("ResidualEstimator::estimateElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
184
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
185

Thomas Witkowski's avatar
Thomas Witkowski committed
186
    double val = 0.0;
187
    std::vector<Operator*>::iterator it;
188
    std::vector<double*>::iterator itfac;
Thomas Witkowski's avatar
Thomas Witkowski committed
189
    Element *el = elInfo->getElement();
190
    double det = elInfo->getDet();
191
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
192
    double est_el = el->getEstimation(row);
193
194
    double h2 = h2_from_det(det, dim);

Thomas Witkowski's avatar
Thomas Witkowski committed
195
    for (int iq = 0; iq < nPoints; iq++) {
196
197
198
      riq[iq] = 0.0;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
199
    for (int system = 0; system < nSystems; system++) {
200
201
202
203
204

      if (matrix[system] == NULL) 
	continue;

      // init assemblers
205
206
      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
	   itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
207
	   it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
208
	   ++it, ++itfac) {
209
	if (*itfac == NULL || **itfac != 0.0) {
210
211
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
	}
212
213
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
216
217
      if (C0) {
	for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
	     it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd(); 
	     ++it) {
218
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);	
Thomas Witkowski's avatar
Thomas Witkowski committed
219
	}
220
221
222
223
224
225
226
227
228
229
230
      }

      if (timestep) {
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
	uhOld[system]->getLocalVector(el, uhOldEl[system]);
  
	// ===== time and element residuals       
	if (C0 || C3) {   
	  uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
	  
231
232
	  if (C3 && uhOldQP && system == std::max(row, 0)) {
	    val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
233
	    for (int iq = 0; iq < nPoints; iq++) {
234
235
236
237
238
239
240
241
242
	      double tiq = (uhQP[iq] - uhOldQP[iq]);
	      val += quad->getWeight(iq) * tiq * tiq;
	    }
	    double v = C3 * det * val;
	    est_t_sum += v;
	    est_t_max = max(est_t_max, v);
	  }
	}
      }
243
           
244
      if (C0) {  
245
246
	for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
	     itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
247
	     it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
248
	     ++it, ++itfac) {
249
	  if (*itfac == NULL || **itfac != 0.0) {
250
251
252
253
254
255
256
257
258
259
260
261
	    if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
	      uhQP = GET_MEMORY(double, nPoints);
	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	    }
	    if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
	      grdUh_qp = new WorldVector<double>[nPoints];
	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
	    }
	    if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) { 
	      D2uhqp = new WorldMatrix<double>[nPoints];
	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
	    }
262
263
264
265
	  }
	}
	
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
266
	  nPoints, 
267
268
269
270
271
272
273
274
275
	  uhQP,
	  grdUh_qp,
	  D2uhqp,
	  uhOldQP,
	  NULL,  // grdUhOldQP 
	  NULL,  // D2UhOldQP
	  matrix[system], 
	  fh[system],
	  quad,
276
	  riq);
277
278
279
280
      }     
    }

    // add integral over r square
281
    val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    for (int iq = 0; iq < nPoints; iq++) {
283
284
      val += quad->getWeight(iq) * riq[iq] * riq[iq];
    }
285
286
   
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM) {
287
      val = C0 * h2 * h2 * det * val;
288
    } else {
289
      val = C0 * h2 * det * val;
290
    }
291
292
293
294
295

    est_el += val;

    // ===== jump residuals 
    if (C1 && (dim > 1)) {
296
      int dow = Global::getGeo(WORLD);
297

298
      for (int face = 0; face < neighbours_; face++) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
299
	Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
300
301
302
	if (neigh && neigh->getMark()) {      
	  int oppV = elInfo->getOppVertex(face);
	      
303
304
	  el->sortFaceIndices(face, &faceIndEl_);
	  neigh->sortFaceIndices(oppV, &faceIndNeigh_);
305
306
307
	    
	  neighInfo->setElement(const_cast<Element*>(neigh));
	  neighInfo->setFillFlag(Mesh::FILL_COORDS);
308
	      	
309
310
311
312
313
314
315
316
317
	  for (int i = 0; i < dow; i++)
	    neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
		
	  // periodic leaf data ?
	  ElementData *ldp = el->getElementData()->getElementData(PERIODIC);

	  bool periodicCoords = false;

	  if (ldp) {
318
319
	    std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	    std::list<LeafDataPeriodic::PeriodicInfo>& infoList = 
320
		dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
321

322
	    for (it = infoList.begin(); it != infoList.end(); ++it) {
323
324
	      if (it->elementSide == face) {
		for (int i = 0; i < dim; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
		  int i1 = faceIndEl_[i];
		  int i2 = faceIndNeigh_[i];
327

Thomas Witkowski's avatar
Thomas Witkowski committed
328
329
		  int j = 0;
		  for (; j < dim; j++) {
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
		    if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, 
								   dim),
						      face,
						      j)) {
		      break;
		    }
		  }

		  TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
		      
		  neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
		}
		periodicCoords = true;
		break;
	      }
	    }
	  }
      
	  if (!periodicCoords) {
	    for (int i = 0; i < dim; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
352
	      int i1 = faceIndEl_[i];
	      int i2 = faceIndNeigh_[i];
	      for (int j = 0; j < dow; j++)
353
354
355
356
357
358
359
360
361
		neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
	    }
	  }
	      
	  Parametric *parametric = mesh->getParametric();
	  if (parametric) {
	    neighInfo = parametric->addParametricInfo(neighInfo);
	  }
	      
362
	  double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh_));
363
	      
364
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
365
	    jump_[iq].set(0.0);
366
367
368
	  }
	     

Thomas Witkowski's avatar
Thomas Witkowski committed
369
	  for (int system = 0; system < nSystems; system++) {	
370
371
372
	    if (matrix[system] == NULL) 
	      continue;
	      
373
374
	    uh[system]->getLocalVector(el, uhEl[system]);	
	    uh[system]->getLocalVector(neigh, uhNeigh[system]);
375
			
376
	    for (int iq = 0; iq < nPointsSurface_; iq++) {
377
	      (*lambda_)[face] = 0.0;
378
	      for (int i = 0; i < dim; i++) {
379
		(*lambda_)[faceIndEl_[i]] = surfaceQuad_->getLambda(iq, i);
380
381
	      }
		  
382
383
	      basFcts[system]->evalGrdUh(*lambda_, 
					 grdLambda, 
384
					 uhEl[system], 
385
					 &grdUhEl_[iq]);
386
		  
387
	      (*lambda_)[oppV] = 0.0;
388
	      for (int i = 0; i < dim; i++) {
389
		(*lambda_)[faceIndNeigh_[i]] = surfaceQuad_->getLambda(iq, i);
390
391
	      }
		  
392
393
	      basFcts[system]->evalGrdUh(*lambda_, 
					 *lambdaNeigh_, 
394
395
					 uhNeigh[system], 
					 &grdUhNeigh_[iq]);
396
		  
397
	      grdUhEl_[iq] -= grdUhNeigh_[iq];
398
399
	    }				

400
	    std::vector<double*>::iterator fac;
401
402
403
404
405

	    for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
		 ++it, ++fac) {
406

407
	      if (*fac == NULL || **fac != 0.0) {
408
409
410
411
412
413
414
415
416
417
418
419
420
421
		for (int iq = 0; iq < nPointsSurface_; iq++) {
		  localJump_[iq].set(0.0);
		}
		
		(*it)->weakEvalSecondOrder(nPointsSurface_,
					   grdUhEl_.getValArray(),
					   localJump_.getValArray());
		double factor = *fac ? **fac : 1.0;
		if (factor != 1.0) {
		  for (int i = 0; i < nPointsSurface_; i++) {
		    localJump_[i] *= factor;
		  }
		}
		
422
		for (int i = 0; i < nPointsSurface_; i++) {
423
		  jump_[i] += localJump_[i];
424
		}
425
426
	      }		
	    }
427
428
	  }
	      
429
	  val = 0.0;
430
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
431
	    val += surfaceQuad_->getWeight(iq) * (jump_[iq] * jump_[iq]);
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
	  }
	      
	  double d = 0.5 * (det + detNeigh);

	  if (norm == NO_NORM || norm == L2_NORM)
	    val *= C1 * h2_from_det(d, dim) * d;
	  else
	    val *= C1 * d;
	      
	  if (parametric) {
	    neighInfo = parametric->removeParametricInfo(neighInfo);
	  }

	  neigh->setEstimation(neigh->getEstimation(row) + val, row);
	  est_el += val;
	} 
      } 
       
450
      val = fh[std::max(row, 0)]->
451
	getBoundaryManager()->
452
	boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
453
454
455
456
457
458
459
460
461
462
463
464
465
466
      if (norm == NO_NORM || norm == L2_NORM)
	val *= C1 * h2;
      else
	val *= C1;
	
      est_el += val;
    } 
  

    el->setEstimation(est_el, row);

    est_sum += est_el;
    est_max = max(est_max, est_el);

467
    elInfo->getElement()->setMark(0);  
468
469
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
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
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
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
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658



#endif


#if 0


  void ResidualEstimator::estimateElement(ElInfo *elInfo)
  {
    FUNCNAME("ResidualEstimator::estimateElement()");

    TEST_EXIT(!C0 && !C2 && !C3)("NOT YET!\n");

    Element *neigh;
    std::vector<Operator*>::iterator it;
    Element *el = elInfo->getElement();
    double det = elInfo->getDet();
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    double est_el = el->getEstimation(row);

    for (int system = 0; system < nSystems; system++) {

      if (matrix[system] == NULL) 
	continue;

      // init assemblers
      std::vector<Operator*>::iterator it;

      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
	   it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
	   ++it) {
	(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
      }

      if (timestep) {
	uhOld[system]->getLocalVector(el, uhOldEl[system]);
      }          
    }

    double val = 0.0;
    double h2 = h2_from_det(det, dim);
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM) {
      val = C0 * h2 * h2 * det * val;
    } else {
      val = C0 * h2 * det * val;
    }

    est_el += val;

    // ===== jump residuals 
    if (C1 && (dim > 1)) {
      int dow = Global::getGeo(WORLD);

      for (int face = 0; face < neighbours_; face++) {  
	neigh = const_cast<Element*>(elInfo->getNeighbour(face));
	if (neigh && neigh->getMark()) {      
	  int oppV = elInfo->getOppVertex(face);
	      
	  el->sortFaceIndices(face, &faceIndEl_);
	  neigh->sortFaceIndices(oppV, &faceIndNeigh_);
	    
	  neighInfo->setElement(const_cast<Element*>(neigh));
	  neighInfo->setFillFlag(Mesh::FILL_COORDS);
	      	
	  for (int i = 0; i < dow; i++)
	    neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
		
	  bool periodicCoords = false;     
	  if (!periodicCoords) {
	    for (int i = 0; i < dim; i++) {
	      int i1 = faceIndEl_[i];
	      int i2 = faceIndNeigh_[i];
	      for (int j = 0; j < dow; j++)
		neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
	    }
	  }
	      
	  double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh_));
	      
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
	    jump_[iq].set(0.0);
	  }	     

	  for (int system = 0; system < nSystems; system++) {	
	    if (matrix[system] == NULL) 
	      continue;
	      
	    uh[system]->getLocalVector(el, uhEl[system]);	
	    uh[system]->getLocalVector(neigh, uhNeigh[system]);
			
	    for (int iq = 0; iq < nPointsSurface_; iq++) {
	      (*lambda_)[face] = 0.0;
	      for (int i = 0; i < dim; i++) {
		(*lambda_)[faceIndEl_[i]] = surfaceQuad_->getLambda(iq, i);
	      }
		  
	      basFcts[system]->evalGrdUh(*lambda_, 
					 grdLambda, 
					 uhEl[system], 
					 &grdUhEl_[iq]);
		  
	      (*lambda_)[oppV] = 0.0;
	      for (int i = 0; i < dim; i++) {
		(*lambda_)[faceIndNeigh_[i]] = surfaceQuad_->getLambda(iq, i);
	      }
		  
	      basFcts[system]->evalGrdUh(*lambda_, 
					 *lambdaNeigh_, 
					 uhNeigh[system], 
					 &grdUhNeigh_[iq]);
		  
	      grdUhEl_[iq] -= grdUhNeigh_[iq];
	    }				

	    std::vector<double*>::iterator fac;

	    for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
		 ++it, ++fac) {
	      for (int iq = 0; iq < nPointsSurface_; iq++) {
		localJump_[iq].set(0.0);
	      }
		  
	      (*it)->weakEvalSecondOrder(nPointsSurface_,
					 grdUhEl_.getValArray(),
					 localJump_.getValArray());
	      double factor = *fac ? **fac : 1.0;
	      if (factor != 1.0) {
		for (int i = 0; i < nPointsSurface_; i++) {
		  localJump_[i] *= factor;
		}
	      }
		  
	      for (int i = 0; i < nPointsSurface_; i++) {
		jump_[i] += localJump_[i];
	      }
	    }				     
	  }
	      
	  val = 0.0;
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
	    val += surfaceQuad_->getWeight(iq) * (jump_[iq] * jump_[iq]);
	  }
	      
	  double d = 0.5 * (det + detNeigh);

	  if (norm == NO_NORM || norm == L2_NORM)
	    val *= C1 * h2_from_det(d, dim) * d;
	  else
	    val *= C1 * d;
	      
	  neigh->setEstimation(neigh->getEstimation(row) + val, row);
	  est_el += val;
	} 
      } 
       
      val = fh[std::max(row, 0)]->
	getBoundaryManager()->
	boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
      if (norm == NO_NORM || norm == L2_NORM)
	val *= C1 * h2;
      else
	val *= C1;
	
      est_el += val;
    } 
  

    el->setEstimation(est_el, row);

    est_sum += est_el;
    est_max = max(est_max, est_el);

    elInfo->getElement()->setMark(0);   
  }

#endif









659
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
660
	 int nPoints,
661
	 const double *uhIq,
662
663
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
664
	 const double *uhOldIq,
665
666
667
668
669
670
671
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
	 double *result)
  {
672
673
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
674
675
676
677
678
679

    // lhs
    for (it = const_cast<DOFMatrix*>(A)->getOperatorsBegin(),
	   fac = const_cast<DOFMatrix*>(A)->getOperatorEstFactorBegin(); 
	 it != const_cast<DOFMatrix*>(A)->getOperatorsEnd(); 
	 ++it, ++fac) {
680
681
682
     
      double factor = *fac ? **fac : 1.0;

683
684
      if (factor) {
	if (D2UhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
685
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
686
687
688
	}

	if (grdUhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
689
690
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
691
692
693
	}
	
	if (uhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
694
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
695
696
697
698
699
700
701
702
703
704
	}
      }
    }
    
    // rhs
    for (it = const_cast<DOFVector<double>*>(fh)->getOperatorsBegin(),
	 fac = const_cast<DOFVector<double>*>(fh)->getOperatorEstFactorBegin(); 
	 it != const_cast<DOFVector<double>*>(fh)->getOperatorsEnd(); 
	 ++it, ++fac) {

705
706
      double factor = *fac ? **fac : 1.0;

707
708
709
      if (factor) {
	if ((*it)->getUhOld()) {
	  if (D2UhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
710
	    (*it)->evalSecondOrder(nPoints, 
711
712
713
714
				   uhOldIq, grdUhOldIq, D2UhOldIq, 
				   result, factor);
	  }
	  if (grdUhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
715
	    (*it)->evalFirstOrderGrdPsi(nPoints, 
716
717
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
718
	    (*it)->evalFirstOrderGrdPhi(nPoints, 
719
720
721
722
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
	  }
	  if (uhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
723
	    (*it)->evalZeroOrder(nPoints, 
724
725
726
727
				 uhOldIq, grdUhOldIq, D2UhOldIq, 
				 result, -factor);
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
728
729
	  double *fx = GET_MEMORY(double, nPoints);
	  for (int iq = 0; iq < nPoints; iq++) {
730
731
	    fx[iq] = 0.0;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	  (*it)->getC(elInfo, nPoints, fx);
733

Thomas Witkowski's avatar
Thomas Witkowski committed
734
	  for (int iq = 0; iq < nPoints; iq++) {
735
736
	    result[iq] -= factor * fx[iq];
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
737
	  FREE_MEMORY(fx, double, nPoints);
738
739
740
741
742
743
744
	}
      }
    }    
  }


}