ResidualEstimator.cc 16.3 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)
  {
22
23
    FUNCNAME("ResidualEstimator::ResidualEstimator()");

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

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

37

38
39
40
  void ResidualEstimator::init(double ts)
  {
    FUNCNAME("ResidualEstimator::init()");
41

42
    timestep = ts;
Thomas Witkowski's avatar
Thomas Witkowski committed
43
    nSystems = static_cast<int>(uh.size());
44

Thomas Witkowski's avatar
Thomas Witkowski committed
45
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
46
47

    dim = mesh->getDim();
48
49
    basFcts = new const BasisFunction*[nSystems];
    quadFast = new FastQuadrature*[nSystems];
50
51

    degree = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
52
    for (int system = 0; system < nSystems; system++) {
53
      basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
54
      degree = std::max(degree, basFcts[system]->getDegree());
55
56
57
58
    }
    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
59
    nPoints = quad->getNumPoints();
60
61

    Flag flag = INIT_PHI | INIT_GRD_PHI;
62
63
    if (degree > 2)
      flag |= INIT_D2_PHI;    
64

65
    for (int system = 0; system < nSystems; system++)
66
67
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
							       *quad, 
68
							       flag);    
69
  
70
71
72
    uhEl = new double*[nSystems];
    uhNeigh = new double*[nSystems];
    uhOldEl = timestep ? new double*[nSystems] : NULL;
73

Thomas Witkowski's avatar
Thomas Witkowski committed
74
    for (int system = 0; system < nSystems; system++) {
75
76
      uhEl[system] = new double[basFcts[system]->getNumber()]; 
      uhNeigh[system] = new double[basFcts[system]->getNumber()];
77
      if (timestep)
78
	uhOldEl[system] = new double[basFcts[system]->getNumber()];
79
80
    }

81
82
83
    uhQP = timestep ? new double[nPoints] : NULL;
    uhOldQP = timestep ? new double[nPoints] : NULL;
    riq = new double[nPoints];
84
85
86
    grdUh_qp = NULL;
    D2uhqp = NULL;

87
    // clear error indicators and mark elements for jumpRes
88
89
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    while (elInfo) {
      elInfo->getElement()->setEstimation(0.0, row);
      elInfo->getElement()->setMark(1);
      elInfo = stack.traverseNext(elInfo);
    }

    est_sum = 0.0;
    est_max = 0.0;
    est_t_sum = 0.0;
    est_t_max = 0.0;

    traverseFlag = 
      Mesh::FILL_NEIGH      |
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND      |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_DET        |
      Mesh::CALL_LEAF_EL;
109
110
    neighInfo = mesh->createNewElInfo();

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

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

	for (std::vector<Operator*>::iterator it = matrix[system]->getOperators().begin();
	     it != matrix[system]->getOperators().end(); ++it)
	  secondOrderTerms[system] = secondOrderTerms[system] || (*it)->secondOrderTerms();
      }
131
    }
132
133
  }

134

135
136
137
138
  void ResidualEstimator::exit(bool output)
  {
    FUNCNAME("ResidualEstimator::exit()");

139
140
141
142
143
144
145
146
147
148
149
150
#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

151
152
153
    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

Thomas Witkowski's avatar
Thomas Witkowski committed
154
    for (int system = 0; system < nSystems; system++) {
155
156
      delete [] uhEl[system];
      delete [] uhNeigh[system];
157
      if (timestep)
158
	delete [] uhOldEl[system];
159
160
    }

161
162
    delete [] uhEl;
    delete [] uhNeigh;
163
164

    if (timestep) {
165
166
167
      delete [] uhOldEl;
      delete [] uhQP;
      delete [] uhOldQP;
168
    } else {
169
170
      if (uhQP != NULL)
	delete [] uhQP;
171
172
173
174
    }

    if (output) {
      MSG("estimate   = %.8e\n", est_sum);
175
      if (C3)
176
177
178
	MSG("time estimate   = %.8e\n", est_t_sum);
    }

179
180
181
    delete [] riq;
    delete [] basFcts;
    delete [] quadFast;
182

183
184
185
186
    if (grdUh_qp != NULL)
      delete [] grdUh_qp;
    if (D2uhqp != NULL)
      delete [] D2uhqp;
187

188
    if (C1 && (dim > 1)) {
189
190
      delete lambdaNeigh;
      delete lambda;
191
192
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
193
    delete neighInfo;
194
195
  }

196

197
  void ResidualEstimator::estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo)
198
  {    
199
200
    FUNCNAME("ResidualEstimator::estimateElement()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
    Element *el = elInfo->getElement();
    double est_el = el->getEstimation(row);
205
206
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator itfac;
207

208
    // === Init assemblers. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
209
    for (int system = 0; system < nSystems; system++) {
210
211
212
      if (matrix[system] == NULL) 
	continue;

213
214
215
216
217
      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)
218
219
220
221
222
223
224
225
	if (*itfac == NULL || **itfac != 0.0) {
	  if (dualElInfo) {
	    (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo, 
								   dualElInfo->largeElInfo,
								   quad);
	  } else
	    (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
	}
226

227
      if (C0 > 0.0)
228
229
230
231
232
233
234
235
	for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it) {
	  if (dualElInfo) {
	    (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo, 
								   dualElInfo->largeElInfo,
								   quad);
	  } else
	    (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
	}
236
237
238
239
240
241
242
243
244
245
246
247
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
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()");

    TEST_EXIT(dualElInfo)("Not yet implemented!\n");

    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;

274
      if (timestep && uhOld[system]) {
275
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
276
	uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
277
  
278
279
280
	// === Compute time error. ===

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

323
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
324
	  nPoints, 
325
326
327
328
329
330
	  uhQP,
	  grdUh_qp,
	  D2uhqp,
	  uhOldQP,
	  NULL,  // grdUhOldQP 
	  NULL,  // D2UhOldQP
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
	  dofMat, 
	  dofVec,
333
	  quad,
334
	  riq);
335
336
337
338
      }     
    }

    // add integral over r square
339
    double result = 0.0;
340
    for (int iq = 0; iq < nPoints; iq++)
341
      result += quad->getWeight(iq) * riq[iq] * riq[iq];
342
   
343
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
344
      result = C0 * h2 * h2 * det * result;
345
    else
346
      result = C0 * h2 * det * result;
347
    
348
349
    return result;
  }
350

351

352
353
354
355
356
357
358
359
360
361
362
  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);
363

364
365
    for (int face = 0; face < nNeighbours; face++) {  
      Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
366

367
368
369
370
      if (dualElInfo) {
	int smallFace = DualTraverse::getFace(dualElInfo, face);
      }

371
372
      if (!(neigh && neigh->getMark()))
	continue;
373

374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
      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) {
394
	    for (int i = 0; i < dim; i++) {
395
396
	      int i1 = faceIndEl[i];
	      int i2 = faceIndNeigh[i];
397
	      
398
399
400
401
	      int j = 0;
	      for (; j < dim; j++)
		if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), face, j))
		  break;
402
	      
403
	      TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
404
	      
405
	      neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
406
	    }
407
408
	    periodicCoords = true;
	    break;
409
	  }
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
	}
      }  // 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);	  
425

426
427
428
429
430
431
432
433
      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;
434
	      
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
	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);
	    
464
	    (*it)->weakEvalSecondOrder(grdUhEl, localJump);
465

466
467
468
469
470
471
472
473
474
475
476
477
478
479
	    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
480

481
482
      double d = 0.5 * (det + detNeigh);
   
483
      if (norm == NO_NORM || norm == L2_NORM)
484
	val *= C1 * h2_from_det(d, dim) * d;
485
      else
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
	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;
502

503
    return result;
504
505
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
506

507
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	 int nPoints,
509
	 const double *uhIq,
510
511
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
512
	 const double *uhOldIq,
513
514
515
516
517
518
519
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
	 double *result)
  {
520
521
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
522
523

    // lhs
Thomas Witkowski's avatar
Thomas Witkowski committed
524
525
    for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin(); 
	 it != A->getOperatorsEnd(); ++it, ++fac) {
526
527
528
     
      double factor = *fac ? **fac : 1.0;

529
      if (factor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
530
531
	if (D2UhIq)
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);	
532
533

	if (grdUhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
536
537
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
538
539
	if (uhIq)
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);	
540
541
542
543
      }
    }
    
    // rhs
Thomas Witkowski's avatar
Thomas Witkowski committed
544
545
    for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin(); 
	 it != fh->getOperatorsEnd(); ++it, ++fac) {
546

547
548
      double factor = *fac ? **fac : 1.0;

549
550
      if (factor) {
	if ((*it)->getUhOld()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
553
	  if (D2UhOldIq)
	    (*it)->evalSecondOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, factor);
	  
554
	  if (grdUhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
555
556
	    (*it)->evalFirstOrderGrdPsi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
	    (*it)->evalFirstOrderGrdPhi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
557
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
560

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

565
	  for (int iq = 0; iq < nPoints; iq++)
566
567
568
569
570
571
572
573
	    result[iq] -= factor * fx[iq];
	}
      }
    }    
  }


}