ResidualEstimator.cc 16.8 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
    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;
105
    
106
    traverseFlag = 
107
108
109
110
111
112
113
      Mesh::FILL_NEIGH         |
      Mesh::FILL_COORDS        |
      Mesh::FILL_OPP_COORDS    |
      Mesh::FILL_BOUND         |
      Mesh::FILL_GRD_LAMBDA    |
      Mesh::FILL_DET           |
      Mesh::FILL_LOCAL_INDICES |
114
      Mesh::CALL_LEAF_EL;
115
116
    neighInfo = mesh->createNewElInfo();

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

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

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

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

143

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

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

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

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

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

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

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

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

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

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

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

205

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

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

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

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

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

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


    // === 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()");

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

    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;

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

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

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

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

365

366
367
368
369
370
371
372
373
374
375
376
  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);
377

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

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

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

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

440
441
442
443
444
445
446
447
      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;
448
	      
449
	uh[system]->getLocalVector(elInfo, uhEl[system]);	
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
477
	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);
	    
478
	    (*it)->weakEvalSecondOrder(grdUhEl, localJump);
479

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

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

517
    return result;
518
519
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
520

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

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

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

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

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

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

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

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


}