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

9
10
11
12
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "mpi.h"
#endif

13
14
namespace AMDiS {

15
  ResidualEstimator::ResidualEstimator(std::string name, int r) 
16
17
18
    : Estimator(name, r),
      C0(1.0), 
      C1(1.0), 
19
      C2(0.0), 
20
21
      C3(1.0),
      jumpResidualOnly(false)
22
  {
23
24
    FUNCNAME("ResidualEstimator::ResidualEstimator()");

25
26
27
28
29
30
31
32
33
    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;
34

35
36
37
38
    if (C1 != 0.0 && C0 == 0.0 && C3 == 0.0)
      jumpResidualOnly = true;
      

39
    TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n");
40
41
  }

42

43
44
45
  void ResidualEstimator::init(double ts)
  {
    FUNCNAME("ResidualEstimator::init()");
46

47
    timestep = ts;
Thomas Witkowski's avatar
Thomas Witkowski committed
48
    nSystems = static_cast<int>(uh.size());
49

Thomas Witkowski's avatar
Thomas Witkowski committed
50
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
51
52

    dim = mesh->getDim();
53
54
    basFcts = new const BasisFunction*[nSystems];
    quadFast = new FastQuadrature*[nSystems];
55
56

    degree = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
57
    for (int system = 0; system < nSystems; system++) {
58
      basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
59
      degree = std::max(degree, basFcts[system]->getDegree());
60
61
62
63
    }
    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
64
    nPoints = quad->getNumPoints();
65
66

    Flag flag = INIT_PHI | INIT_GRD_PHI;
67
68
    if (degree > 2)
      flag |= INIT_D2_PHI;    
69

70
    for (int system = 0; system < nSystems; system++)
71
72
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
							       *quad, 
73
							       flag);    
74
  
75
76
77
    uhEl = new double*[nSystems];
    uhNeigh = new double*[nSystems];
    uhOldEl = timestep ? new double*[nSystems] : NULL;
78

Thomas Witkowski's avatar
Thomas Witkowski committed
79
    for (int system = 0; system < nSystems; system++) {
80
81
      uhEl[system] = new double[basFcts[system]->getNumber()]; 
      uhNeigh[system] = new double[basFcts[system]->getNumber()];
82
      if (timestep)
83
	uhOldEl[system] = new double[basFcts[system]->getNumber()];
84
85
    }

86
87
88
    uhQP = timestep ? new double[nPoints] : NULL;
    uhOldQP = timestep ? new double[nPoints] : NULL;
    riq = new double[nPoints];
89
90
91
    grdUh_qp = NULL;
    D2uhqp = NULL;

92
    // clear error indicators and mark elements for jumpRes
93
94
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
    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;
114
115
    neighInfo = mesh->createNewElInfo();

116
    // === Prepare date for computing jump residual. ===
117
118
119
120
121
122
123
124
125
126
    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);
127
128
129
130
131

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

132
133
134
	if (matrix[system] == NULL)
	  continue;

135
136
137
138
	for (std::vector<Operator*>::iterator it = matrix[system]->getOperators().begin();
	     it != matrix[system]->getOperators().end(); ++it)
	  secondOrderTerms[system] = secondOrderTerms[system] || (*it)->secondOrderTerms();
      }
139
    }
140
141
  }

142

143
144
145
146
  void ResidualEstimator::exit(bool output)
  {
    FUNCNAME("ResidualEstimator::exit()");

147
148
149
150
151
152
153
154
155
156
157
158
#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

159
160
161
    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

Thomas Witkowski's avatar
Thomas Witkowski committed
162
    for (int system = 0; system < nSystems; system++) {
163
164
      delete [] uhEl[system];
      delete [] uhNeigh[system];
165
      if (timestep)
166
	delete [] uhOldEl[system];
167
168
    }

169
170
    delete [] uhEl;
    delete [] uhNeigh;
171
172

    if (timestep) {
173
174
175
      delete [] uhOldEl;
      delete [] uhQP;
      delete [] uhOldQP;
176
    } else {
177
178
      if (uhQP != NULL)
	delete [] uhQP;
179
180
181
182
    }

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

187
188
189
    delete [] riq;
    delete [] basFcts;
    delete [] quadFast;
190

191
192
193
194
    if (grdUh_qp != NULL)
      delete [] grdUh_qp;
    if (D2uhqp != NULL)
      delete [] D2uhqp;
195

196
    if (C1 && (dim > 1)) {
197
198
      delete lambdaNeigh;
      delete lambda;
199
200
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
201
    delete neighInfo;
202
203
  }

204

205
  void ResidualEstimator::estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo)
206
  {    
207
208
    FUNCNAME("ResidualEstimator::estimateElement()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
211
212
    Element *el = elInfo->getElement();
    double est_el = el->getEstimation(row);
213
214
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator itfac;
215

216
    // === Init assemblers. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
217
    for (int system = 0; system < nSystems; system++) {
218
219
220
      if (matrix[system] == NULL) 
	continue;

221
222
223
224
225
      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)
226
227
228
229
230
231
232
	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)
233
234
235
	    (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo, 
								   dualElInfo->largeElInfo,
								   quad);
236
237
	  else
	    (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);	  
238
	}
239

240
      if (C0 > 0.0)
241
	for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it) {
242
	  if (dualElInfo)
243
244
245
	    (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo, 
								   dualElInfo->largeElInfo,
								   quad);
246
247
	  else
	    (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);	  
248
	}
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
    }


    // === 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;
    est_max = max(est_max, est_el);
  }


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

273
    TEST_EXIT(!dualElInfo)("Not yet implemented!\n");
274
275
276
277
278
279
280
281
282
283
284
285
286

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

    for (int iq = 0; iq < nPoints; iq++)
      riq[iq] = 0.0;

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

287
      if (timestep && uhOld[system]) {
288
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
289
	uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
290
  
291
292
293
	// === Compute time error. ===

	if (C0 > 0.0 || C3 > 0.0) {   
294
295
296
	  uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
	  
297
	  if (C3 > 0.0 && uhOldQP && system == std::max(row, 0)) {
298
	    double result = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
299
	    for (int iq = 0; iq < nPoints; iq++) {
300
	      double tiq = (uhQP[iq] - uhOldQP[iq]);
301
	      result += quad->getWeight(iq) * tiq * tiq;
302
	    }
303
	    double v = C3 * det * result;
304
305
306
307
308
	    est_t_sum += v;
	    est_t_max = max(est_t_max, v);
	  }
	}
      }
309
           
310
      // === Compute element residual. ===
311
312
313
314
      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
315
	for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
316
	     it != dofMat->getOperatorsEnd();  ++it, ++itfac) {
317
	  if (*itfac == NULL || **itfac != 0.0) {
318
	    if (uhQP == NULL && (*it)->zeroOrderTerms()) {
319
	      uhQP = new double[nPoints];
320
321
	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	    }
322
323
	    if (grdUh_qp == NULL && 
		((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
324
325
326
	      grdUh_qp = new WorldVector<double>[nPoints];
	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
	    }
327
	    if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) { 
328
329
330
	      D2uhqp = new WorldMatrix<double>[nPoints];
	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
	    }
331
332
333
	  }
	}
	
334
335
	// === Compute the element residual and store it in irq. ===

336
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
337
	  nPoints, 
338
339
340
341
342
343
	  uhQP,
	  grdUh_qp,
	  D2uhqp,
	  uhOldQP,
	  NULL,  // grdUhOldQP 
	  NULL,  // D2UhOldQP
Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
	  dofMat, 
	  dofVec,
346
	  quad,
347
	  riq);
348
349
350
351
      }     
    }

    // add integral over r square
352
    double result = 0.0;
353
    for (int iq = 0; iq < nPoints; iq++)
354
      result += quad->getWeight(iq) * riq[iq] * riq[iq];
355
   
356
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
357
      result = C0 * h2 * h2 * det * result;
358
    else
359
      result = C0 * h2 * det * result;
360
    
361
362
    return result;
  }
363

364

365
366
367
368
369
370
371
372
373
374
375
  double ResidualEstimator::computeJumpResidual(ElInfo *elInfo, 
						DualElInfo *dualElInfo)
  {
    FUNCNAME("ResidualEstimator::computeJumpResidual()");

    double result = 0.0;
    int dow = Global::getGeo(WORLD);
    Element *el = elInfo->getElement();
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    double det = elInfo->getDet();
    double h2 = h2_from_det(det, dim);
376

377
378
    for (int face = 0; face < nNeighbours; face++) {  
      Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
379

380
381
382
383
      if (dualElInfo) {
	int smallFace = DualTraverse::getFace(dualElInfo, face);
      }

384
385
      if (!(neigh && neigh->getMark()))
	continue;
386

387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
      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];
      
      // 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) {
407
	    for (int i = 0; i < dim; i++) {
408
409
	      int i1 = faceIndEl[i];
	      int i2 = faceIndNeigh[i];
410
	      
411
412
413
414
	      int j = 0;
	      for (; j < dim; j++)
		if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), face, j))
		  break;
415
	      
416
	      TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
417
	      
418
	      neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
419
	    }
420
421
	    periodicCoords = true;
	    break;
422
	  }
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
	}
      }  // if (ldp)
      
      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];
	}
      }
	
      Parametric *parametric = mesh->getParametric();
      if (parametric)
	neighInfo = parametric->addParametricInfo(neighInfo);	  
438

439
440
441
442
443
444
445
446
      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 || secondOrderTerms[system] == false) 
	  continue;
447
	      
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
	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];
	} // 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);
	    
477
	    (*it)->weakEvalSecondOrder(grdUhEl, localJump);
478

479
480
481
482
483
484
485
486
487
488
489
490
491
492
	    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
493

494
495
      double d = 0.5 * (det + detNeigh);
   
496
      if (norm == NO_NORM || norm == L2_NORM)
497
	val *= C1 * h2_from_det(d, dim) * d;
498
      else
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
	val *= C1 * d;
      
      if (parametric)
	neighInfo = parametric->removeParametricInfo(neighInfo);
      
      neigh->setEstimation(neigh->getEstimation(row) + val, row);
      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;
515

516
    return result;
517
518
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
519

520
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
521
	 int nPoints,
522
	 const double *uhIq,
523
524
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
525
	 const double *uhOldIq,
526
527
528
529
530
531
532
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
	 double *result)
  {
533
534
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
535
536

    // lhs
Thomas Witkowski's avatar
Thomas Witkowski committed
537
538
    for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin(); 
	 it != A->getOperatorsEnd(); ++it, ++fac) {
539
540
541
     
      double factor = *fac ? **fac : 1.0;

542
      if (factor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
543
544
	if (D2UhIq)
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);	
545
546

	if (grdUhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
547
548
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
549
550
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
	if (uhIq)
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);	
553
554
555
556
      }
    }
    
    // rhs
Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
    for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin(); 
	 it != fh->getOperatorsEnd(); ++it, ++fac) {
559

560
561
      double factor = *fac ? **fac : 1.0;

562
563
      if (factor) {
	if ((*it)->getUhOld()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
564
565
566
	  if (D2UhOldIq)
	    (*it)->evalSecondOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, factor);
	  
567
	  if (grdUhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
568
569
	    (*it)->evalFirstOrderGrdPsi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
	    (*it)->evalFirstOrderGrdPhi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
570
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
571
572
573

	  if (uhOldIq)
	    (*it)->evalZeroOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);	  
574
	} else {
575
	  std::vector<double> fx(nPoints, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
576
	  (*it)->getC(elInfo, nPoints, fx);
577

578
	  for (int iq = 0; iq < nPoints; iq++)
579
580
581
582
583
584
585
586
	    result[iq] -= factor * fx[iq];
	}
      }
    }    
  }


}