ResidualEstimator.cc 20.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
#include "ResidualEstimator.h"
#include "Operator.h"
#include "DOFMatrix.h"
#include "DOFVector.h"
#include "Assembler.h"
#include "Traverse.h"
19
#include "Initfile.h"
20

21
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
22
#include <mpi.h>
23
#include "parallel/MeshDistributor.h"
24
25
#endif

26
27
namespace AMDiS {

28
  ResidualEstimator::ResidualEstimator(std::string name, int r) 
29
    : Estimator(name, r),
30
31
      C0(0.0), 
      C1(0.0), 
32
      C2(0.0), 
33
      C3(0.0),
34
      jumpResidualOnly(false)
35
  {
36
37
    FUNCNAME("ResidualEstimator::ResidualEstimator()");

38
39
40
41
    Parameters::get(name + "->C0", C0);
    Parameters::get(name + "->C1", C1);
    Parameters::get(name + "->C2", C2);
    Parameters::get(name + "->C3", C3);
42
43
44
45
46

    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;
47

48
49
50
    if (C1 != 0.0 && C0 == 0.0 && C3 == 0.0)
      jumpResidualOnly = true;
      
51
    TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n");
52
53
  }

54

55
56
57
  void ResidualEstimator::init(double ts)
  {
    FUNCNAME("ResidualEstimator::init()");
58

59
    timestep = ts;
Thomas Witkowski's avatar
Thomas Witkowski committed
60
    nSystems = static_cast<int>(uh.size());
61

Thomas Witkowski's avatar
Thomas Witkowski committed
62
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
63
64

    dim = mesh->getDim();
65
66
    basFcts = new const BasisFunction*[nSystems];
    quadFast = new FastQuadrature*[nSystems];
67
68

    degree = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    for (int system = 0; system < nSystems; system++) {
70
      basFcts[system] = uh[system]->getFeSpace()->getBasisFcts();
71
      degree = std::max(degree, basFcts[system]->getDegree());
72
73
74
75
    }
    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
76
    nPoints = quad->getNumPoints();
77
78

    Flag flag = INIT_PHI | INIT_GRD_PHI;
79
80
    if (degree > 2)
      flag |= INIT_D2_PHI;    
81

82
    for (int system = 0; system < nSystems; system++)
83
84
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
							       *quad, 
85
							       flag);    
86
  
87
88
89
90
    uhEl.resize(nSystems);
    uhNeigh.resize(nSystems);
    if (timestep)
      uhOldEl.resize(nSystems);
91

Thomas Witkowski's avatar
Thomas Witkowski committed
92
    for (int system = 0; system < nSystems; system++) {
93
94
      uhEl[system].change_dim(basFcts[system]->getNumber()); 
      uhNeigh[system].change_dim(basFcts[system]->getNumber());
95
      if (timestep)
96
	uhOldEl[system].change_dim(basFcts[system]->getNumber());
97
98
    }

99
100
101
102
    if (timestep) {
      uhQP.change_dim(nPoints);
      uhOldQP.change_dim(nPoints);
    }
103
104

    riq.change_dim(nPoints);
105

106
    // clear error indicators and mark elements for jumpRes
107
108
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
109
    while (elInfo) {
110
111
112
      // SIMON: DEL LINE BELOW
      elInfo->getElement()->setEstimation(0.0, row);

113
114
115
116
117
118
119
120
      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;
121

122
    traverseFlag = 
123
124
125
126
127
128
      Mesh::FILL_NEIGH      |
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND      |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_DET        |
129
      Mesh::CALL_LEAF_EL;
130
131
    neighInfo = mesh->createNewElInfo();

132
    // === Prepare date for computing jump residual. ===
133
134
135
136
137
138
139
140
141
142
    if (C1 > 0.0 && dim > 1) {
      surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
      nPointsSurface = surfaceQuad->getNumPoints();
      grdUhEl.resize(nPointsSurface);
      grdUhNeigh.resize(nPointsSurface);
      jump.resize(nPointsSurface);
      localJump.resize(nPointsSurface);
      nNeighbours = Global::getGeo(NEIGH, dim);
      lambdaNeigh = new DimVec<WorldVector<double> >(dim, NO_INIT);
      lambda = new DimVec<double>(dim, NO_INIT);
143
144
145
146
147

      secondOrderTerms.resize(nSystems);
      for (int system = 0; system < nSystems; system++) {
	secondOrderTerms[system] = false;

148
149
150
	if (matrix[system] == NULL)
	  continue;

151
152
153
154
	for (std::vector<Operator*>::iterator it = matrix[system]->getOperators().begin();
	     it != matrix[system]->getOperators().end(); ++it)
	  secondOrderTerms[system] = secondOrderTerms[system] || (*it)->secondOrderTerms();
      }
155
    }
156
157
158
159

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    initParallel();
#endif
160
161
  }

162

163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  void ResidualEstimator::initParallel()
  {
    FUNCNAME("ResidualEstimator::initParallel()");

    if (C1 == 0.0)
      return;

    DOFVector<WorldVector<double> > coords(uh[0]->getFeSpace(), "tmp");
    mesh->getDofIndexCoords(uh[0]->getFeSpace(), coords);

    InteriorBoundary &intBoundary = 
      MeshDistributor::globalMeshDistributor->getIntBoundary();

    ElInfo *elInfo = mesh->createNewElInfo();
    elInfo->setFillFlag(Mesh::FILL_COORDS);

    StdMpi<vector<double> > stdMpiDet(MeshDistributor::globalMeshDistributor->getMpiComm());
    StdMpi<vector<vector<WorldVector<double> > > > stdMpiGrdUh(MeshDistributor::globalMeshDistributor->getMpiComm());

183
184
185
186
187
188
189
190
    RankToBoundMap allBounds = intBoundary.getOther();
    allBounds.insert(intBoundary.getOwn().begin(), intBoundary.getOwn().end());

    for (RankToBoundMap::iterator it = allBounds.begin();
	 it != allBounds.end(); ++it) {

      vector<BoundaryObject> subBound;

191
192
193
194
195
196
      for (unsigned int i = 0; i < it->second.size(); i++) {
	BoundaryObject &bObj = it->second[i].rankObj;
	if (bObj.subObj == VERTEX)
	  continue;

	bObj.el->getSubBoundary(bObj, subBound);
197
      }
198

199
200
      if (subBound.size() == 0)
	continue;
201

202
      WorldVector<int> faceIndices;
203

204
205
206
207
208
209
210
211
212
213
214
      for (unsigned int i = 0; i < subBound.size(); i++) {
	Element *el = subBound[i].el;	  
	int oppV = subBound[i].ithObj;
	
	elInfo->setElement(el);
	el->sortFaceIndices(oppV, faceIndices);	
	for (int k = 0; k <= dim; k++)
	  elInfo->getCoord(k) = coords[el->getDof(k, 0)];
	
	double detNeigh = abs(elInfo->calcGrdLambda(*lambdaNeigh));
	stdMpiDet.getSendData(it->first).push_back(detNeigh);
215

216
217
218
219
220
221
222
223
224
225
226
227
228
	for (int system = 0; system < nSystems; system++) {
	  if (matrix[system] == NULL || secondOrderTerms[system] == false)
	    continue;
	  
	  uh[system]->getLocalVector(el, uhNeigh[system]);
	  
	  for (int iq = 0; iq < nPointsSurface; iq++) {
	    
	    (*lambda)[oppV] = 0.0;
	    for (int k = 0; k < dim; k++)
	      (*lambda)[faceIndices[k]] = surfaceQuad->getLambda(iq, k);
	    
	    basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
229
	  }
230
231
	  
	  stdMpiGrdUh.getSendData(it->first).push_back(grdUhNeigh);	  
232
233
	}
      }
234
235
236

      stdMpiDet.recv(it->first);
      stdMpiGrdUh.recv(it->first);      
237
238
239
240
241
    }

    stdMpiDet.updateSendDataSize();
    stdMpiGrdUh.updateSendDataSize();

242
#if 0
243
244
245
246
247
248
249
250
251
252
253
254
255
256
    RankToBoundMap& ownBound = intBoundary.getOwn();

    for (RankToBoundMap::iterator it = ownBound.begin();
	 it != ownBound.end(); ++it) {

      for (unsigned int i = 0; i < it->second.size(); i++) {
	BoundaryObject &bObj = it->second[i].rankObj;
	if (bObj.subObj != VERTEX) {
	  stdMpiDet.recv(it->first);
	  stdMpiGrdUh.recv(it->first);
	  break;
	}
      }
    }
257
#endif
258
259
260
261
262

    stdMpiDet.startCommunication();
    stdMpiGrdUh.startCommunication();


263
264
    for (RankToBoundMap::iterator it = allBounds.begin();
	 it != allBounds.end(); ++it) {
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
      vector<BoundaryObject> subBound;
	
      for (unsigned int i = 0; i < it->second.size(); i++) {
	BoundaryObject &bObj = it->second[i].rankObj;
	if (bObj.subObj == VERTEX)
	  continue;

	bObj.el->getSubBoundary(bObj, subBound);
      }
      
      if (subBound.size() == 0)
	continue;

      TEST_EXIT_DBG(subBound.size() == stdMpiDet.getRecvData(it->first).size())
	("Should not happen: %d %d from rank %d\n", subBound.size(), stdMpiDet.getRecvData(it->first).size(), it->first);
      
      TEST_EXIT_DBG(subBound.size() == stdMpiGrdUh.getRecvData(it->first).size())
	("Should not happen!\n");      
283
284
285
286
287

      for (unsigned int i = 0; i < subBound.size(); i++) {
	elBoundDet[subBound[i]] = stdMpiDet.getRecvData(it->first)[i];
	elBoundGrdUhNeigh[subBound[i]] = stdMpiGrdUh.getRecvData(it->first)[i];
      }
288
289
290
291
292
293
294
    }    

    MSG("INIT PARALLEL FINISCHED!\n");
  }
#endif


295
296
297
298
  void ResidualEstimator::exit(bool output)
  {
    FUNCNAME("ResidualEstimator::exit()");

299
300
301
302
303
304
305
306
307
308
309
310
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double send_est_sum = est_sum;
    double send_est_max = est_max;
    double send_est_t_sum = est_t_sum;
    double send_est_t_max = est_t_max;

    MPI::COMM_WORLD.Allreduce(&send_est_sum, &est_sum, 1, MPI_DOUBLE, MPI_SUM);
    MPI::COMM_WORLD.Allreduce(&send_est_max, &est_max, 1, MPI_DOUBLE, MPI_MAX);
    MPI::COMM_WORLD.Allreduce(&send_est_t_sum, &est_t_sum, 1, MPI_DOUBLE, MPI_SUM);
    MPI::COMM_WORLD.Allreduce(&send_est_t_max, &est_t_max, 1, MPI_DOUBLE, MPI_MAX);
#endif

311
312
313
314
    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

    if (output) {
315
      MSG("estimate for component %d = %.8e\n", row, est_sum);
316
      if (C3)
317
	MSG("time estimate for component %d = %.8e\n", row, est_t_sum);
318
319
    }

320
321
    delete [] basFcts;
    delete [] quadFast;
322

323
    if (C1 && (dim > 1)) {
324
325
      delete lambdaNeigh;
      delete lambda;
326
327
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
328
    delete neighInfo;
329
330
  }

331

332
  void ResidualEstimator::estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo)
333
  {    
334
335
    FUNCNAME("ResidualEstimator::estimateElement()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    Element *el = elInfo->getElement();
339
340
    double est_el = el->getEstimation(row);
    // SIMON    double est_el = 0.0;
341
342
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator itfac;
343

344
    // === Init assemblers. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
345
    for (int system = 0; system < nSystems; system++) {
346
347
348
      if (matrix[system] == NULL) 
	continue;

349
350
351
352
353
      DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[system]);
      DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[system]);

      for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
	   it != dofMat->getOperatorsEnd(); ++it, ++itfac)
354
355
356
357
358
359
360
	if (*itfac == NULL || **itfac != 0.0) {	  
	  // If the estimator must only compute the jump residual but there are no
	  // second order terms in the operator, it can be skipped.
	  if (jumpResidualOnly && (*it)->secondOrderTerms() == false)
	    continue;
	  
	  if (dualElInfo)
361
362
363
	    (*it)->getAssembler()->initElement(dualElInfo->smallElInfo, 
					       dualElInfo->largeElInfo,
					       quad);
364
	  else
365
	    (*it)->getAssembler()->initElement(elInfo, NULL, quad);	  
366
	}
367

368
      if (C0 > 0.0)
369
	for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it) {
370
	  if (dualElInfo)
371
372
373
	    (*it)->getAssembler()->initElement(dualElInfo->smallElInfo, 
					       dualElInfo->largeElInfo,
					       quad);
374
	  else
375
	    (*it)->getAssembler()->initElement(elInfo, NULL, quad);	  
376
	}
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
    }


    // === Compute element residuals and time error estimation. ===
    if (C0 || C3)
      est_el += computeElementResidual(elInfo, dualElInfo);

    // === Compute jump residuals. ===
    if (C1 && dim > 1)
      est_el += computeJumpResidual(elInfo, dualElInfo);

    // === Update global residual variables. ===
    el->setEstimation(est_el, row);
    el->setMark(0);
    est_sum += est_el;
392
    est_max = std::max(est_max, est_el);
393
394
395
396
397
398
399
400
  }


  double ResidualEstimator::computeElementResidual(ElInfo *elInfo, 
						   DualElInfo *dualElInfo)
  {
    FUNCNAME("ResidualEstimator::computeElementResidual()");

401
    TEST_EXIT(!dualElInfo)("Not yet implemented!\n");
402
403
404
405
406

    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator itfac;
    double det = elInfo->getDet();
    double h2 = h2_from_det(det, dim);
407
    riq = 0.0;
408
409
410
411
412

    for (int system = 0; system < nSystems; system++) {
      if (matrix[system] == NULL) 
	continue;

413
      if (timestep && uhOld[system]) {
414
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
415
	uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
416
  
417
418
419
	// === Compute time error. ===

	if (C0 > 0.0 || C3 > 0.0) {   
420
421
422
	  uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
	  
423
	  if (C3 > 0.0 && system == std::max(row, 0)) {
424
	    double result = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
425
	    for (int iq = 0; iq < nPoints; iq++) {
426
	      double tiq = (uhQP[iq] - uhOldQP[iq]);
427
	      result += quad->getWeight(iq) * tiq * tiq;
428
	    }
429
	    double v = C3 * det * result;
430
	    est_t_sum += v;
431
	    est_t_max = std::max(est_t_max, v);
432
433
434
	  }
	}
      }
435
           
436
      // === Compute element residual. ===
437
438
439
440
      if (C0 > 0.0) {
	DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[system]);
	DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[system]);
  
Thomas Witkowski's avatar
Thomas Witkowski committed
441
	for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
442
	     it != dofMat->getOperatorsEnd();  ++it, ++itfac) {
443
	  if (*itfac == NULL || **itfac != 0.0) {
444
	    if ((*it)->zeroOrderTerms()) {
445
	      uhQP.change_dim(nPoints);
446
447
	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	    }
448
	    if ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi()) {
449
450
	      grdUhQp.change_dim(nPoints);
	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUhQp);
451
	    }
452
	    if (degree > 2 && (*it)->secondOrderTerms()) {
453
454
	      D2UhQp.change_dim(nPoints);
	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2UhQp);
455
	    }
456
457
458
	  }
	}
	
459
460
	// === Compute the element residual and store it in irq. ===

461
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
462
	  nPoints, 
463
	  uhQP,
464
465
	  grdUhQp,
	  D2UhQp,
466
	  uhOldQP,
467
468
	  grdUhOldQp,  // not used
	  D2UhOldQp,  // not used
Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
	  dofMat, 
	  dofVec,
471
	  quad,
472
	  riq);
473
474
475
476
      }     
    }

    // add integral over r square
477
    double result = 0.0;
478
    for (int iq = 0; iq < nPoints; iq++)
479
      result += quad->getWeight(iq) * riq[iq] * riq[iq];
480
   
481
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
482
      result = C0 * h2 * h2 * det * result;
483
    else
484
      result = C0 * h2 * det * result;
485
    
486
487
    return result;
  }
488

489

490
491
492
493
494
495
496
497
498
499
  double ResidualEstimator::computeJumpResidual(ElInfo *elInfo, 
						DualElInfo *dualElInfo)
  {
    FUNCNAME("ResidualEstimator::computeJumpResidual()");

    double result = 0.0;
    Element *el = elInfo->getElement();
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    double det = elInfo->getDet();
    double h2 = h2_from_det(det, dim);
500
    BoundaryObject blub;
501

502
503
    for (int face = 0; face < nNeighbours; face++) {  
      Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
504

505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
      bool parallelMode = false;

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      if (neigh == NULL) {
	BoundaryObject testObj(el, elInfo->getType(), EDGE, face);
	
	if (elBoundDet.count(testObj)) {
	  parallelMode = true;
	  blub = testObj;
	} else {
	  continue;
	}
      }

      if (!parallelMode && !(neigh && neigh->getMark()))
	continue;
#else
522
523
      if (!(neigh && neigh->getMark()))
	continue;
524
#endif
525

526
      int oppV = elInfo->getOppVertex(face);	
527
      el->sortFaceIndices(face, faceIndEl);
528
529
530
531
532
533
534

      if (!parallelMode) {
	neigh->sortFaceIndices(oppV, faceIndNeigh);	
	neighInfo->setElement(const_cast<Element*>(neigh));
	neighInfo->setFillFlag(Mesh::FILL_COORDS);     
	neighInfo->getCoord(oppV) = elInfo->getOppCoord(face);
      }
535
536
537
538
539
540
541
542
543
544
545
      
      // periodic leaf data ?
      ElementData *ldp = el->getElementData()->getElementData(PERIODIC);	
      bool periodicCoords = false;
      
      if (ldp) {
	typedef std::list<LeafDataPeriodic::PeriodicInfo> PerInfList;
	PerInfList& infoList = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
	
	for (PerInfList::iterator it = infoList.begin(); it != infoList.end(); ++it) {
	  if (it->elementSide == face) {
546
	    for (int i = 0; i < dim; i++) {
547
548
	      int i1 = faceIndEl[i];
	      int i2 = faceIndNeigh[i];
549
	      
550
551
552
553
	      int j = 0;
	      for (; j < dim; j++)
		if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), face, j))
		  break;
554
	      
555
	      TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
556
	      
557
	      neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
558
	    }
559
560
	    periodicCoords = true;
	    break;
561
	  }
562
563
564
	}
      }  // if (ldp)
      
565
      if (!periodicCoords && !parallelMode) {
566
567
568
	for (int i = 0; i < dim; i++) {
	  int i1 = faceIndEl[i];
	  int i2 = faceIndNeigh[i];
569
	  neighInfo->getCoord(i2) = elInfo->getCoord(i1);
570
571
572
	}
      }
	
573
      double detNeigh = 0.0;
574
      Parametric *parametric = mesh->getParametric();
575
576
577
578
579
580
581
      if (!parallelMode) {
	if (parametric)
	  neighInfo = parametric->addParametricInfo(neighInfo);	  

	detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
      } else {
	TEST_EXIT_DBG(!parametric)("No yet implemented!\n");
582

583
584
585
	detNeigh = elBoundDet[blub];
      }
      
586
587
588
589
590
591
592
           
      for (int iq = 0; iq < nPointsSurface; iq++)
	jump[iq].set(0.0);     
      
      for (int system = 0; system < nSystems; system++) {
	if (matrix[system] == NULL || secondOrderTerms[system] == false) 
	  continue;
593
	      
594
	uh[system]->getLocalVector(el, uhEl[system]);	
595
596
	if (!parallelMode)
	  uh[system]->getLocalVector(neigh, uhNeigh[system]);
597
598
599
600
601
602
	  
	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);
	  
603
	  basFcts[system]->evalGrdUh(*lambda, grdLambda, uhEl[system], grdUhEl[iq]);
604
605
606
607
608
609
610
611
612
613

	  if (!parallelMode) {
	    (*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]);
	  } else {
	    grdUhNeigh[iq] = elBoundGrdUhNeigh[blub][iq];
	  }
614
	  
615

616
617
618
619
620
621
622
623
624
625
626
627
628
629
	  
	  grdUhEl[iq] -= grdUhNeigh[iq];
	} // for iq				
	
	std::vector<double*>::iterator fac;
	std::vector<Operator*>::iterator it;
	DOFMatrix *mat = const_cast<DOFMatrix*>(matrix[system]);
        for (it = mat->getOperatorsBegin(), fac = mat->getOperatorEstFactorBegin(); 
	     it != mat->getOperatorsEnd(); ++it, ++fac) {
	
	  if (*fac == NULL || **fac != 0.0) {
	    for (int iq = 0; iq < nPointsSurface; iq++)
	      localJump[iq].set(0.0);
	    
630
	    (*it)->weakEvalSecondOrder(grdUhEl, localJump);
631

632
633
634
635
636
637
638
639
640
641
642
643
644
645
	    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];
	  }		
	} // for (it = ...
      } // for system
    
      double val = 0.0;
      for (int iq = 0; iq < nPointsSurface; iq++)
	val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
Thomas Witkowski's avatar
Thomas Witkowski committed
646

647
648
      double d = 0.5 * (det + detNeigh);
   
649
      if (norm == NO_NORM || norm == L2_NORM)
650
	val *= C1 * h2_from_det(d, dim) * d;
651
      else
652
	val *= C1 * d;
653
654
655
656
657
658
659

      if (!parallelMode) {
	if (parametric)
	  neighInfo = parametric->removeParametricInfo(neighInfo);      
	neigh->setEstimation(neigh->getEstimation(row) + val, row);
      }

660
661
662
663
664
665
666
667
668
669
      result += val;
    } // for face
    
    double 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;   
    result += val;
670

671
    return result;
672
673
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
674

675
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
676
	 int nPoints,
677
678
679
680
681
682
	 const mtl::dense_vector<double>& uhIq,
	 const mtl::dense_vector<WorldVector<double> > &grdUhIq,
	 const mtl::dense_vector<WorldMatrix<double> > &D2UhIq,
	 const mtl::dense_vector<double>& uhOldIq,
	 const mtl::dense_vector<WorldVector<double> > &grdUhOldIq,
	 const mtl::dense_vector<WorldMatrix<double> > &D2UhOldIq,
683
684
685
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
686
	 ElementVector& result)
687
  {
688
689
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
690
691

    // lhs
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
    for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin(); 
	 it != A->getOperatorsEnd(); ++it, ++fac) {
694
695
696
     
      double factor = *fac ? **fac : 1.0;

697
      if (factor) {
698
699
700
701
702
703
	if (num_rows(D2UhIq) > 0)
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);	

	if (num_rows(grdUhIq) > 0) {
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
704
705
	}
	
706
	if (num_rows(uhIq) > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
707
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);	
708
709
710
711
      }
    }
    
    // rhs
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
    for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin(); 
	 it != fh->getOperatorsEnd(); ++it, ++fac) {
714

715
716
      double factor = *fac ? **fac : 1.0;

717
718
      if (factor) {
	if ((*it)->getUhOld()) {
719
	  if (num_rows(D2UhOldIq) > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
720
721
	    (*it)->evalSecondOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, factor);
	  
722
	  if (num_rows(grdUhOldIq) > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
723
724
	    (*it)->evalFirstOrderGrdPsi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
	    (*it)->evalFirstOrderGrdPhi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
725
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
726

727
	  if (num_rows(uhOldIq) > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
728
	    (*it)->evalZeroOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);	  
729
	} else {
730
	  ElementVector fx(nPoints, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
731
	  (*it)->getC(elInfo, nPoints, fx);
732

733
	  for (int iq = 0; iq < nPoints; iq++)
734
735
736
737
738
739
740
741
	    result[iq] -= factor * fx[iq];
	}
      }
    }    
  }


}