ResidualEstimator.cc 21.6 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


22
23
24
25
26
27
#include "ResidualEstimator.h"
#include "Operator.h"
#include "DOFMatrix.h"
#include "DOFVector.h"
#include "Assembler.h"
#include "Traverse.h"
28
#include "Initfile.h"
29

30
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
31
#include <mpi.h>
32
#include "parallel/MeshDistributor.h"
Siqi Ling's avatar
Siqi Ling committed
33
#include "parallel/ParallelDebug.h"
34
35
#endif

36
37
using namespace std;

38
39
namespace AMDiS {

40
  ResidualEstimator::ResidualEstimator(std::string name, int r)
41
    : Estimator(name, r),
42
43
44
      C0(0.0),
      C1(0.0),
      C2(0.0),
45
      C3(0.0),
46
      jumpResidualOnly(false)
47
  {
48
49
    FUNCNAME("ResidualEstimator::ResidualEstimator()");

50
51
52
53
    Parameters::get(name + "->C0", C0);
    Parameters::get(name + "->C1", C1);
    Parameters::get(name + "->C2", C2);
    Parameters::get(name + "->C3", C3);
54
55
56
57
58

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

60
61
    if (C1 != 0.0 && C0 == 0.0 && C3 == 0.0)
      jumpResidualOnly = true;
62

63
    TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n");
64
65
  }

66

67
68
  void ResidualEstimator::init(double ts)
  {
69
    FUNCNAME_DBG("ResidualEstimator::init()");
70

71
    timestep = ts;
Thomas Witkowski's avatar
Thomas Witkowski committed
72
    nSystems = static_cast<int>(uh.size());
73

Thomas Witkowski's avatar
Thomas Witkowski committed
74
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
75
76

    dim = mesh->getDim();
77
78
    basFcts = new const BasisFunction*[nSystems];
    quadFast = new FastQuadrature*[nSystems];
79
80

    degree = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
81
    for (int system = 0; system < nSystems; system++) {
82
      basFcts[system] = uh[system]->getFeSpace()->getBasisFcts();
83
      degree = std::max(degree, basFcts[system]->getDegree());
84
85
86
87
    }
    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
88
    nPoints = quad->getNumPoints();
89
90

    Flag flag = INIT_PHI | INIT_GRD_PHI;
91
    if (degree > 2)
92
      flag |= INIT_D2_PHI;
93

94
    for (int system = 0; system < nSystems; system++)
95
96
97
98
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system],
							       *quad,
							       flag);

99
100
101
102
    uhEl.resize(nSystems);
    uhNeigh.resize(nSystems);
    if (timestep)
      uhOldEl.resize(nSystems);
103

Thomas Witkowski's avatar
Thomas Witkowski committed
104
    for (int system = 0; system < nSystems; system++) {
105
      uhEl[system].change_dim(basFcts[system]->getNumber());
106
      uhNeigh[system].change_dim(basFcts[system]->getNumber());
107
      if (timestep)
108
	uhOldEl[system].change_dim(basFcts[system]->getNumber());
109
110
    }

111
112
113
114
    if (timestep) {
      uhQP.change_dim(nPoints);
      uhOldQP.change_dim(nPoints);
    }
115
116

    riq.change_dim(nPoints);
117

118
    // clear error indicators and mark elements for jumpRes
119
120
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
121
    while (elInfo) {
122
123
124
      // SIMON: DEL LINE BELOW
      elInfo->getElement()->setEstimation(0.0, row);

125
126
127
128
129
130
131
132
      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;
133

134
    traverseFlag =
135
136
137
138
139
140
      Mesh::FILL_NEIGH      |
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND      |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_DET        |
141
      Mesh::CALL_LEAF_EL;
142
143
    neighInfo = mesh->createNewElInfo();

144
    // === Prepare date for computing jump residual. ===
145
146
147
148
149
150
151
152
153
154
    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);
155
156
157
158
159

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

160
	if (matrix[system] == NULL)
161
162
	  continue;

163
164
165
166
	for (std::vector<Operator*>::iterator it = matrix[system]->getOperators().begin();
	     it != matrix[system]->getOperators().end(); ++it)
	  secondOrderTerms[system] = secondOrderTerms[system] || (*it)->secondOrderTerms();
      }
167
    }
168

169
170
171
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    initParallel();
#endif
172
173
  }

174

175
176
177
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  void ResidualEstimator::initParallel()
  {
178
    FUNCNAME_DBG("ResidualEstimator::initParallel()");
179
180
181
182

    if (C1 == 0.0)
      return;

183
184
185
186
    size_t comp = 0;
    for(; comp < uh.size(); comp++)
      if(uh[comp]->getFeSpace()->getMesh() == mesh)
	break;
187

188
189
    TEST_EXIT_DBG(comp < uh.size())("No mesh in the solution is found?.\n");
    TEST_EXIT_DBG(*(uh[comp]->getFeSpace()->getAdmin()) == mesh->getDofAdmin(0))("The admin is not the first one in the mesh.\n");
190

191
    DOFVector<WorldVector<double> > coords(uh[comp]->getFeSpace(), "tmp");
Thomas Witkowski's avatar
Thomas Witkowski committed
192
    mesh->getDofIndexCoords(coords);
193

194
    Parallel::InteriorBoundary &intBoundary =
195
      Parallel::MeshDistributor::globalMeshDistributor->getIntBoundary(0);
196

197
#ifndef NDEBUG
Siqi Ling's avatar
Siqi Ling committed
198
199
200
    // Make sure interior boundary is correct
    Parallel::ParallelDebug::testInteriorBoundary(*(Parallel::MeshDistributor::globalMeshDistributor));
#endif
201
202
203
204

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

205
206
    Parallel::StdMpi<vector<double> > stdMpiDet(Parallel::MeshDistributor::globalMeshDistributor->getMpiComm(0));
    Parallel::StdMpi<vector<vector<WorldVector<double> > > > stdMpiGrdUh(Parallel::MeshDistributor::globalMeshDistributor->getMpiComm(0));
207

208
    Parallel::RankToBoundMap allBounds = intBoundary.getOther();
209
210
    allBounds.insert(intBoundary.getOwn().begin(), intBoundary.getOwn().end());

211
    for (Parallel::RankToBoundMap::iterator it = allBounds.begin();
212
213
214
215
	 it != allBounds.end(); ++it) {

      vector<BoundaryObject> subBound;

216
217
218
219
220
      for (unsigned int i = 0; i < it->second.size(); i++) {
	BoundaryObject &bObj = it->second[i].rankObj;
	if (bObj.subObj == VERTEX)
	  continue;

221
222
	bObj.el = intBoundary.getElementPtr(bObj.elIndex, mesh);
	TEST_EXIT_DBG(bObj.el)("No element found in the map.\n");
223
	bObj.el->getSubBoundary(bObj, subBound);
224
      }
225

226
227
      if (subBound.size() == 0)
	continue;
228

229
      WorldVector<int> faceIndices;
230

231
      for (unsigned int i = 0; i < subBound.size(); i++) {
232
	Element *el = subBound[i].el;
233
	int oppV = subBound[i].ithObj;
234

235
	elInfo->setElement(el);
236
	el->sortFaceIndices(oppV, faceIndices);
237
238
	for (int k = 0; k <= dim; k++)
	  elInfo->getCoord(k) = coords[el->getDof(k, 0)];
239

240
241
	double detNeigh = abs(elInfo->calcGrdLambda(*lambdaNeigh));
	stdMpiDet.getSendData(it->first).push_back(detNeigh);
242
243


244
	for (int system = 0; system < nSystems; system++) {
245
	  if (matrix[system] == NULL || secondOrderTerms[system] == false)
246
	    continue;
247

248
	  uh[system]->getLocalVector(el, uhNeigh[system]);
249

250
	  for (int iq = 0; iq < nPointsSurface; iq++) {
251

252
253
254
	    (*lambda)[oppV] = 0.0;
	    for (int k = 0; k < dim; k++)
	      (*lambda)[faceIndices[k]] = surfaceQuad->getLambda(iq, k);
255

256
	    basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
257
	  }
258
259

	  stdMpiGrdUh.getSendData(it->first).push_back(grdUhNeigh);
260
261
	}
      }
262
263

      stdMpiDet.recv(it->first);
264
      stdMpiGrdUh.recv(it->first);
265
266
267
268
269
270
271
    }

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

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

273
    for (Parallel::RankToBoundMap::iterator it = allBounds.begin();
274
	 it != allBounds.end(); ++it) {
275
      vector<BoundaryObject> subBound;
276

277
278
279
280
281
      for (unsigned int i = 0; i < it->second.size(); i++) {
	BoundaryObject &bObj = it->second[i].rankObj;
	if (bObj.subObj == VERTEX)
	  continue;

282
283
	bObj.el = intBoundary.getElementPtr(bObj.elIndex, mesh);
	TEST_EXIT(bObj.el)("No element found in the map.\n");
284
285
	bObj.el->getSubBoundary(bObj, subBound);
      }
286

287
288
      if (subBound.size() == 0)
	continue;
289

290
      // Highly possible mesh is not correct, for example: hanging node, if the following error occurs
Siqi Ling's avatar
Siqi Ling committed
291
      // Highly possible mesh is not correct, for example: hanging node
292
      TEST_EXIT_DBG(subBound.size() == stdMpiDet.getRecvData(it->first).size())
293
294
	("Should not happen: rank %d from rank %d should recv %d (recv %d)\n",
	   MPI::COMM_WORLD.Get_rank(), it->first, subBound.size(), stdMpiDet.getRecvData(it->first).size());
295

296
      TEST_EXIT_DBG(subBound.size() == stdMpiGrdUh.getRecvData(it->first).size())
297
	("Should not happen!\n");
298
299
300
301
302

      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];
      }
303
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
304
305

    delete elInfo;
306
307
308
309
  }
#endif


310
311
312
313
  void ResidualEstimator::exit(bool output)
  {
    FUNCNAME("ResidualEstimator::exit()");

314
315
316
317
318
319
320
321
322
323
324
325
#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

326
327
328
329
    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

    if (output) {
330
      MSG("estimate for component %d = %.8e\n", row, est_sum);
331
      if (C3)
332
	MSG("time estimate for component %d = %.8e\n", row, est_t_sum);
333
334
    }

335
336
    delete [] basFcts;
    delete [] quadFast;
337

338
    if (C1 && (dim > 1)) {
339
340
      delete lambdaNeigh;
      delete lambda;
341
342
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
343
    delete neighInfo;
344
345
  }

346

347
  void ResidualEstimator::estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo)
348
  {
349
    FUNCNAME_DBG("ResidualEstimator::estimateElement()");
350

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

Thomas Witkowski's avatar
Thomas Witkowski committed
353
    Element *el = elInfo->getElement();
354
355
    double est_el = el->getEstimation(row);
    // SIMON    double est_el = 0.0;
356
357
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator itfac;
358

359
    // === Init assemblers. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
360
    for (int system = 0; system < nSystems; system++) {
361
      if (matrix[system] == NULL)
362
363
	continue;

364
365
366
367
368
      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)
369
	if (*itfac == NULL || **itfac != 0.0) {
370
371
372
373
	  // 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;
374

375
	  if (dualElInfo)
376
	    (*it)->getAssembler()->initElement(dualElInfo->smallElInfo,
377
378
					       dualElInfo->largeElInfo,
					       quad);
379
	  else
380
	    (*it)->getAssembler()->initElement(elInfo, NULL, quad);
381
	}
382

383
      if (C0 > 0.0)
384
	for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it) {
385
	  if (dualElInfo)
386
	    (*it)->getAssembler()->initElement(dualElInfo->smallElInfo,
387
388
					       dualElInfo->largeElInfo,
					       quad);
389
	  else
390
	    (*it)->getAssembler()->initElement(elInfo, NULL, quad);
391
	}
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
    }


    // === 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;
407
    est_max = std::max(est_max, est_el);
408
409
410
  }


411
  double ResidualEstimator::computeElementResidual(ElInfo *elInfo,
412
413
414
415
						   DualElInfo *dualElInfo)
  {
    FUNCNAME("ResidualEstimator::computeElementResidual()");

416
    TEST_EXIT(!dualElInfo)("Not yet implemented!\n");
417
418
419
420
421

    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator itfac;
    double det = elInfo->getDet();
    double h2 = h2_from_det(det, dim);
422
    riq = 0.0;
423
424

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

428
      if (timestep && uhOld[system]) {
429
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
430
	uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
431

432
433
	// === Compute time error. ===

434
	if (C0 > 0.0 || C3 > 0.0) {
435
436
	  uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
437

438
	  if (C3 > 0.0 && system == std::max(row, 0)) {
439
	    double result = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
440
	    for (int iq = 0; iq < nPoints; iq++) {
441
	      double tiq = (uhQP[iq] - uhOldQP[iq]);
442
	      result += quad->getWeight(iq) * tiq * tiq;
443
	    }
444
	    double v = C3 * det * result;
445
	    est_t_sum += v;
446
	    est_t_max = std::max(est_t_max, v);
447
448
449
	  }
	}
      }
450

451
      // === Compute element residual. ===
452
453
454
      if (C0 > 0.0) {
	DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[system]);
	DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[system]);
455

Thomas Witkowski's avatar
Thomas Witkowski committed
456
	for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
457
	     it != dofMat->getOperatorsEnd();  ++it, ++itfac) {
458
	  if (*itfac == NULL || **itfac != 0.0) {
459
	    if ((*it)->zeroOrderTerms()) {
460
	      uhQP.change_dim(nPoints);
461
	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
462
	    }
463
	    if ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi()) {
464
	      grdUhQp.change_dim(nPoints);
465
	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUhQp);
466
	    }
467
	    if (degree > 2 && (*it)->secondOrderTerms()) {
468
	      D2UhQp.change_dim(nPoints);
469
	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2UhQp);
470
	    }
471
472
	  }
	}
473

474
475
	// === Compute the element residual and store it in irq. ===

476
	r(elInfo,
477
	  nPoints,
478
	  uhQP,
479
480
	  grdUhQp,
	  D2UhQp,
481
	  uhOldQP,
482
483
	  grdUhOldQp,  // not used
	  D2UhOldQp,  // not used
484
	  dofMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
485
	  dofVec,
486
	  quad,
487
	  riq);
488
      }
489
490
491
    }

    // add integral over r square
492
    double result = 0.0;
493
    for (int iq = 0; iq < nPoints; iq++)
494
      result += quad->getWeight(iq) * riq[iq] * riq[iq];
495

496
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
497
      result = C0 * h2 * h2 * det * result;
498
    else
499
      result = C0 * h2 * det * result;
500

501
502
    return result;
  }
503

504

505
  double ResidualEstimator::computeJumpResidual(ElInfo *elInfo,
506
507
						DualElInfo *dualElInfo)
  {
508
    FUNCNAME_DBG("ResidualEstimator::computeJumpResidual()");
509
510
511
512
513
514

    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);
515
    BoundaryObject blub;
516

517
    for (int face = 0; face < nNeighbours; face++) {
518
      Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
519

520
521
522
      bool parallelMode = false;

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
523
      if (neigh == NULL) {
524
	BoundaryObject testObj(el, elInfo->getType(), EDGE, face);
525

526
527
528
529
530
531
532
533
534
535
536
	if (elBoundDet.count(testObj)) {
	  parallelMode = true;
	  blub = testObj;
	} else {
	  continue;
	}
      }

      if (!parallelMode && !(neigh && neigh->getMark()))
	continue;
#else
537
538
      if (!(neigh && neigh->getMark()))
	continue;
539
#endif
540

541
      int oppV = elInfo->getOppVertex(face);
542
      el->sortFaceIndices(face, faceIndEl);
543
544

      if (!parallelMode) {
545
	neigh->sortFaceIndices(oppV, faceIndNeigh);
546
	neighInfo->setElement(const_cast<Element*>(neigh));
547
	neighInfo->setFillFlag(Mesh::FILL_COORDS);
548
549
	neighInfo->getCoord(oppV) = elInfo->getOppCoord(face);
      }
550

551
      // periodic leaf data ?
552
      ElementData *ldp = el->getElementData()->getElementData(PERIODIC);
553
      bool periodicCoords = false;
554

555
556
557
      if (ldp) {
	typedef std::list<LeafDataPeriodic::PeriodicInfo> PerInfList;
	PerInfList& infoList = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
558

559
560
	for (PerInfList::iterator it = infoList.begin(); it != infoList.end(); ++it) {
	  if (it->elementSide == face) {
561
	    for (int i = 0; i < dim; i++) {
562
563
	      int i1 = faceIndEl[i];
	      int i2 = faceIndNeigh[i];
564

565
566
567
568
	      int j = 0;
	      for (; j < dim; j++)
		if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), face, j))
		  break;
569

570
	      TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
571

572
	      neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
573
	    }
574
575
	    periodicCoords = true;
	    break;
576
	  }
577
578
	}
      }  // if (ldp)
579

580
      if (!periodicCoords && !parallelMode) {
581
582
583
	for (int i = 0; i < dim; i++) {
	  int i1 = faceIndEl[i];
	  int i2 = faceIndNeigh[i];
584
	  neighInfo->getCoord(i2) = elInfo->getCoord(i1);
585
586
	}
      }
587

588
      double detNeigh = 0.0;
589
      Parametric *parametric = mesh->getParametric();
590
591
      if (!parallelMode) {
	if (parametric)
592
	  neighInfo = parametric->addParametricInfo(neighInfo);
593
594
595
596

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

598
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
599
	detNeigh = elBoundDet[blub];
600
#endif
601
      }
602
603


604
      for (int iq = 0; iq < nPointsSurface; iq++)
605
606
	jump[iq].set(0.0);

607
      for (int system = 0; system < nSystems; system++) {
608
	if (matrix[system] == NULL || secondOrderTerms[system] == false)
609
	  continue;
610
611

	uh[system]->getLocalVector(el, uhEl[system]);
612
613
	if (!parallelMode)
	  uh[system]->getLocalVector(neigh, uhNeigh[system]);
614

615
616
617
618
	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);
619

620
	  basFcts[system]->evalGrdUh(*lambda, grdLambda, uhEl[system], grdUhEl[iq]);
621
622
623
624
625
626
627
628

	  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 {
629
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
630
	    grdUhNeigh[iq] = elBoundGrdUhNeigh[blub][iq];
631
#endif
632
633
	  }

634
635


636
	  grdUhEl[iq] -= grdUhNeigh[iq];
637
638
	} // for iq

639
640
641
	std::vector<double*>::iterator fac;
	std::vector<Operator*>::iterator it;
	DOFMatrix *mat = const_cast<DOFMatrix*>(matrix[system]);
642
        for (it = mat->getOperatorsBegin(), fac = mat->getOperatorEstFactorBegin();
643
	     it != mat->getOperatorsEnd(); ++it, ++fac) {
644

645
	  if (*fac == NULL || **fac != 0.0) {
646
647
	    for (int iq = 0; iq < nPointsSurface; iq++)
	      localJump[iq].set(0.0);
648

649
	    (*it)->weakEvalSecondOrder(grdUhEl, localJump);
650

651
652
653
654
	    double factor = *fac ? **fac : 1.0;
	    if (factor != 1.0)
	      for (int i = 0; i < nPointsSurface; i++)
		localJump[i] *= factor;
655

656
657
	    for (int i = 0; i < nPointsSurface; i++)
	      jump[i] += localJump[i];
658
	  }
659
660
	} // for (it = ...
      } // for system
661

662
663
664
      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
665

666
      double d = 0.5 * (det + detNeigh);
667

668
      if (norm == NO_NORM || norm == L2_NORM)
669
	val *= C1 * h2_from_det(d, dim) * d;
670
      else
671
	val *= C1 * d;
672
673
674

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

679
680
      result += val;
    } // for face
681

682
    double val = fh[std::max(row, 0)]->getBoundaryManager()->
683
      boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
684
685
686
    if (norm == NO_NORM || norm == L2_NORM)
      val *= C1 * h2;
    else
687
      val *= C1;
688
    result += val;
689

690
    return result;
691
692
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
693

694
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
695
	 int nPoints,
696
697
698
699
700
701
	 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,
702
	 DOFMatrix *A,
703
704
	 DOFVector<double> *fh,
	 Quadrature *quad,
705
	 ElementVector& result)
706
  {
707
708
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
709
710

    // lhs
711
    for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin();
Thomas Witkowski's avatar
Thomas Witkowski committed
712
	 it != A->getOperatorsEnd(); ++it, ++fac) {
713

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

716
      if (factor) {
717
	if (num_rows(D2UhIq) > 0)
718
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
719
720
721
722

	if (num_rows(grdUhIq) > 0) {
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
723
	}
724

725
	if (num_rows(uhIq) > 0)
726
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
727
728
      }
    }
729

730
    // rhs
731
    for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin();
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	 it != fh->getOperatorsEnd(); ++it, ++fac) {
733

734
735
      double factor = *fac ? **fac : 1.0;

736
737
      if (factor) {
	if ((*it)->getUhOld()) {