ResidualEstimator.cc 19.1 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
181
182
183
  void ResidualEstimator::estimateElement(ElInfo *elInfo)
  {
    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;
Thomas Witkowski's avatar
Thomas Witkowski committed
188
    Element *el = elInfo->getElement();
189
    double det = elInfo->getDet();
190
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
191
    double est_el = el->getEstimation(row);
192
193
    double h2 = h2_from_det(det, dim);

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

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

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

      // init assemblers
      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
	   it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
	   ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
207
	(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
208
209
      }

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

      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);
	  
227
228
	  if (C3 && uhOldQP && system == std::max(row, 0)) {
	    val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
229
	    for (int iq = 0; iq < nPoints; iq++) {
230
231
232
233
234
235
236
237
238
	      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);
	  }
	}
      }
239
           
240
241
242
243
      if (C0) {  
	for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(); 
	     it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
	     ++it) {
244
	  if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
245
	    uhQP = GET_MEMORY(double, nPoints);
246
247
	    uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  }
248
	  if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
Thomas Witkowski's avatar
Thomas Witkowski committed
249
	    grdUh_qp = new WorldVector<double>[nPoints];
250
251
	    uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
	  }
252
	  if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) { 
Thomas Witkowski's avatar
Thomas Witkowski committed
253
	    D2uhqp = new WorldMatrix<double>[nPoints];
254
	    uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
255
256
257
258
	  }
	}
	
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
259
	  nPoints, 
260
261
262
263
264
265
266
267
268
	  uhQP,
	  grdUh_qp,
	  D2uhqp,
	  uhOldQP,
	  NULL,  // grdUhOldQP 
	  NULL,  // D2UhOldQP
	  matrix[system], 
	  fh[system],
	  quad,
269
	  riq);
270
271
272
273
      }     
    }

    // add integral over r square
274
    val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
275
    for (int iq = 0; iq < nPoints; iq++) {
276
277
      val += quad->getWeight(iq) * riq[iq] * riq[iq];
    }
278
279
   
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM) {
280
      val = C0 * h2 * h2 * det * val;
281
    } else {
282
      val = C0 * h2 * det * val;
283
    }
284
285
286
287
288

    est_el += val;

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

291
      for (int face = 0; face < neighbours_; face++) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
292
	Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
293
294
295
	if (neigh && neigh->getMark()) {      
	  int oppV = elInfo->getOppVertex(face);
	      
296
297
	  el->sortFaceIndices(face, &faceIndEl_);
	  neigh->sortFaceIndices(oppV, &faceIndNeigh_);
298
299
300
	    
	  neighInfo->setElement(const_cast<Element*>(neigh));
	  neighInfo->setFillFlag(Mesh::FILL_COORDS);
301
	      	
302
303
304
305
306
307
308
309
310
	  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) {
311
312
	    std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	    std::list<LeafDataPeriodic::PeriodicInfo>& infoList = 
313
		dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
314

315
	    for (it = infoList.begin(); it != infoList.end(); ++it) {
316
317
	      if (it->elementSide == face) {
		for (int i = 0; i < dim; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319
		  int i1 = faceIndEl_[i];
		  int i2 = faceIndNeigh_[i];
320

Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
		  int j = 0;
		  for (; j < dim; j++) {
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
		    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
343
344
345
	      int i1 = faceIndEl_[i];
	      int i2 = faceIndNeigh_[i];
	      for (int j = 0; j < dow; j++)
346
347
348
349
350
351
352
353
354
		neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
	    }
	  }
	      
	  Parametric *parametric = mesh->getParametric();
	  if (parametric) {
	    neighInfo = parametric->addParametricInfo(neighInfo);
	  }
	      
355
	  double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh_));
356
	      
357
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
358
	    jump_[iq].set(0.0);
359
360
361
	  }
	     

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

393
	    std::vector<double*>::iterator fac;
394
395
396
397
398

	    for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
		 ++it, ++fac) {
399
	      for (int iq = 0; iq < nPointsSurface_; iq++) {
400
		localJump_[iq].set(0.0);
401
402
	      }
		  
403
404
405
	      (*it)->weakEvalSecondOrder(nPointsSurface_,
					 grdUhEl_.getValArray(),
					 localJump_.getValArray());
406
407
	      double factor = *fac ? **fac : 1.0;
	      if (factor != 1.0) {
408
409
		for (int i = 0; i < nPointsSurface_; i++) {
		  localJump_[i] *= factor;
410
411
412
		}
	      }
		  
413
414
	      for (int i = 0; i < nPointsSurface_; i++) {
		jump_[i] += localJump_[i];
415
416
417
418
	      }
	    }				     
	  }
	      
419
	  val = 0.0;
420
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
421
	    val += surfaceQuad_->getWeight(iq) * (jump_[iq] * jump_[iq]);
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
	  }
	      
	  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;
	} 
      } 
       
440
      val = fh[std::max(row, 0)]->
441
	getBoundaryManager()->
442
	boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
      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);   
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
460
461
462
463
464
465
466
467
468
469
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



#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









649
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
650
	 int nPoints,
651
	 const double *uhIq,
652
653
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
654
	 const double *uhOldIq,
655
656
657
658
659
660
661
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
	 double *result)
  {
662
663
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
664
665
666
667
668
669

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

673
674
      if (factor) {
	if (D2UhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
675
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
676
677
678
	}

	if (grdUhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
681
682
683
	}
	
	if (uhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
684
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
685
686
687
688
689
690
691
692
693
694
	}
      }
    }
    
    // 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) {

695
696
      double factor = *fac ? **fac : 1.0;

697
698
699
      if (factor) {
	if ((*it)->getUhOld()) {
	  if (D2UhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
700
	    (*it)->evalSecondOrder(nPoints, 
701
702
703
704
				   uhOldIq, grdUhOldIq, D2UhOldIq, 
				   result, factor);
	  }
	  if (grdUhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
705
	    (*it)->evalFirstOrderGrdPsi(nPoints, 
706
707
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
708
	    (*it)->evalFirstOrderGrdPhi(nPoints, 
709
710
711
712
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
	  }
	  if (uhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
713
	    (*it)->evalZeroOrder(nPoints, 
714
715
716
717
				 uhOldIq, grdUhOldIq, D2UhOldIq, 
				 result, -factor);
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
	  double *fx = GET_MEMORY(double, nPoints);
	  for (int iq = 0; iq < nPoints; iq++) {
720
721
	    fx[iq] = 0.0;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
722
	  (*it)->getC(elInfo, nPoints, fx);
723

Thomas Witkowski's avatar
Thomas Witkowski committed
724
	  for (int iq = 0; iq < nPoints; iq++) {
725
726
	    result[iq] -= factor * fx[iq];
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
727
	  FREE_MEMORY(fx, double, nPoints);
728
729
730
731
732
733
734
	}
      }
    }    
  }


}