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

namespace AMDiS {

11
  ResidualEstimator::ResidualEstimator(std::string name, int r) 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
    : Estimator(name, r),
      C0(1.0), 
      C1(1.0), 
      C2(1.0), 
      C3(1.0)
  {
    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;
  }

  void ResidualEstimator::init(double ts)
  {
    FUNCNAME("ResidualEstimator::init()");
Thomas Witkowski's avatar
Thomas Witkowski committed
32
    
33
34
35
36
    timestep = ts;

    mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
37
38
    nSystems = static_cast<int>(uh.size());
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
39
40

    dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
    basFcts = GET_MEMORY(const BasisFunction*, nSystems);
    quadFast = GET_MEMORY(FastQuadrature*, nSystems);
43
44

    degree = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
45
    for (int system = 0; system < nSystems; system++) {
46
      basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
47
      degree = std::max(degree, basFcts[system]->getDegree());
48
49
50
51
52
    }

    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
53
    nPoints = quad->getNumPoints();
54
55
56
57
58
59

    Flag flag = INIT_PHI | INIT_GRD_PHI;
    if (degree > 2) {
      flag |= INIT_D2_PHI;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
60
    for (int system = 0; system < nSystems; system++) {
61
62
63
64
65
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
							       *quad, 
							       flag);
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
    uhEl = GET_MEMORY(double*, nSystems);
    uhNeigh = GET_MEMORY(double*, nSystems);
    uhOldEl = timestep ? GET_MEMORY(double*, nSystems) : NULL;
69

Thomas Witkowski's avatar
Thomas Witkowski committed
70
    for (int system = 0; system < nSystems; system++) {
71
      uhEl[system] = GET_MEMORY(double, basFcts[system]->getNumber()); 
72
      uhNeigh[system] = GET_MEMORY(double, basFcts[system]->getNumber());
73
74
75
76
      if (timestep)
	uhOldEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
    uhQP = timestep ? GET_MEMORY(double, nPoints) : NULL;
    uhOldQP = timestep ? GET_MEMORY(double, nPoints) : NULL;
79

Thomas Witkowski's avatar
Thomas Witkowski committed
80
    riq = GET_MEMORY(double, nPoints);
81

82
83
84
    grdUh_qp = NULL;
    D2uhqp = NULL;

85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    TraverseStack stack;
    ElInfo *elInfo = NULL;

    // clear error indicators and mark elements for jumpRes
    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    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
111
112
113
114
115
116
117
118
119
120

    neighInfo = mesh->createNewElInfo();

    // prepare date for computing jump residual
    if (C1 && (dim > 1)) {
      surfaceQuad_ = Quadrature::provideQuadrature(dim - 1, degree);
      nPointsSurface_ = surfaceQuad_->getNumPoints();
      grdUhEl_.resize(nPointsSurface_);
      grdUhNeigh_.resize(nPointsSurface_);
      jump_.resize(nPointsSurface_);
      localJump_.resize(nPointsSurface_);
      neighbours_ = Global::getGeo(NEIGH, dim);
121
122
      lambdaNeigh_ = NEW DimVec<WorldVector<double> >(dim, NO_INIT);
      lambda_ = NEW DimVec<double>(dim, NO_INIT);
123
    }
124
125
126
127
128
129
130
131
132
  }

  void ResidualEstimator::exit(bool output)
  {
    FUNCNAME("ResidualEstimator::exit()");

    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

Thomas Witkowski's avatar
Thomas Witkowski committed
133
    for (int system = 0; system < nSystems; system++) {
134
      FREE_MEMORY(uhEl[system], double, basFcts[system]->getNumber());
135
      FREE_MEMORY(uhNeigh[system], double, basFcts[system]->getNumber());
136
137
138
139
      if (timestep)
	FREE_MEMORY(uhOldEl[system], double, basFcts[system]->getNumber());    
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
    FREE_MEMORY(uhEl, double*, nSystems);
    FREE_MEMORY(uhNeigh, double*, nSystems);
142
143

    if (timestep) {
Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
146
      FREE_MEMORY(uhOldEl, double*, nSystems);
      FREE_MEMORY(uhQP, double, nPoints);
      FREE_MEMORY(uhOldQP, double, nPoints);
147
148
    } else {
      if (uhQP != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
149
	FREE_MEMORY(uhQP, double, nPoints);
150
      }
151
152
153
154
155
156
157
158
159
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
160
161
162
    FREE_MEMORY(riq, double, nPoints);
    FREE_MEMORY(basFcts, const BasisFunction*, nSystems);
    FREE_MEMORY(quadFast, FastQuadrature*, nSystems);
163
164

    if (grdUh_qp != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
165
      FREE_MEMORY(grdUh_qp, WorldVector<double>, nPoints);
166
167
    }
    if (D2uhqp != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
      FREE_MEMORY(D2uhqp, WorldMatrix<double>, nPoints);
169
    }
170

171
172
173
174
175
    if (C1 && (dim > 1)) {
      DELETE lambdaNeigh_;
      DELETE lambda_;
    }

176
    DELETE neighInfo;
177
178
179
  }

  void ResidualEstimator::estimateElement(ElInfo *elInfo)
180
  {    
181
182
    FUNCNAME("ResidualEstimator::estimateElement()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
185
    double val = 0.0;
186
    std::vector<Operator*>::iterator it;
187
    std::vector<double*>::iterator itfac;
Thomas Witkowski's avatar
Thomas Witkowski committed
188
    Element *el = elInfo->getElement();
189
    double det = elInfo->getDet();
190
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
191
    double est_el = el->getEstimation(row);
192
193
    double h2 = h2_from_det(det, dim);

Thomas Witkowski's avatar
Thomas Witkowski committed
194
    for (int iq = 0; iq < nPoints; iq++) {
195
196
197
      riq[iq] = 0.0;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
198
    for (int system = 0; system < nSystems; system++) {
199
200
201
202
203

      if (matrix[system] == NULL) 
	continue;

      // init assemblers
204
205
      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
	   itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
206
	   it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
207
	   ++it, ++itfac) {
208
	if (*itfac == NULL || **itfac != 0.0) {
209
210
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
	}
211
212
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
215
216
      if (C0) {
	for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
	     it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd(); 
	     ++it) {
217
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);	
Thomas Witkowski's avatar
Thomas Witkowski committed
218
	}
219
      }
Naumann, Andreas's avatar
Naumann, Andreas committed
220
	
221
      if (timestep && uhOld[system]) {
222
223
224
225
226
227
228
229
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
	uhOld[system]->getLocalVector(el, uhOldEl[system]);
  
	// ===== time and element residuals       
	if (C0 || C3) {   
	  uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
	  
230
231
	  if (C3 && uhOldQP && system == std::max(row, 0)) {
	    val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
232
	    for (int iq = 0; iq < nPoints; iq++) {
233
234
235
236
237
238
239
240
241
	      double tiq = (uhQP[iq] - uhOldQP[iq]);
	      val += quad->getWeight(iq) * tiq * tiq;
	    }
	    double v = C3 * det * val;
	    est_t_sum += v;
	    est_t_max = max(est_t_max, v);
	  }
	}
      }
242
           
243
      if (C0) {  
244
245
	for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
	     itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
246
	     it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
247
	     ++it, ++itfac) {
248
	  if (*itfac == NULL || **itfac != 0.0) {
249
250
251
252
253
254
255
256
257
258
259
260
	    if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
	      uhQP = GET_MEMORY(double, nPoints);
	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	    }
	    if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
	      grdUh_qp = new WorldVector<double>[nPoints];
	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
	    }
	    if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) { 
	      D2uhqp = new WorldMatrix<double>[nPoints];
	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
	    }
261
262
263
264
	  }
	}
	
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
265
	  nPoints, 
266
267
268
269
270
271
272
273
274
	  uhQP,
	  grdUh_qp,
	  D2uhqp,
	  uhOldQP,
	  NULL,  // grdUhOldQP 
	  NULL,  // D2UhOldQP
	  matrix[system], 
	  fh[system],
	  quad,
275
	  riq);
276
277
278
279
      }     
    }

    // add integral over r square
280
    val = 0.0;
281
    for (int iq = 0; iq < nPoints; iq++)
282
      val += quad->getWeight(iq) * riq[iq] * riq[iq];
283
   
284
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
285
      val = C0 * h2 * h2 * det * val;
286
    else
287
      val = C0 * h2 * det * val;
288
    
289
290
291
292
    est_el += val;

    // ===== jump residuals 
    if (C1 && (dim > 1)) {
293
      int dow = Global::getGeo(WORLD);
294

295
      for (int face = 0; face < neighbours_; face++) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
296
	Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
297
298
299
	if (neigh && neigh->getMark()) {      
	  int oppV = elInfo->getOppVertex(face);
	      
300
301
	  el->sortFaceIndices(face, &faceIndEl_);
	  neigh->sortFaceIndices(oppV, &faceIndNeigh_);
302
303
304
	    
	  neighInfo->setElement(const_cast<Element*>(neigh));
	  neighInfo->setFillFlag(Mesh::FILL_COORDS);
305
	      	
306
307
308
309
310
311
312
313
314
	  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) {
315
316
	    std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	    std::list<LeafDataPeriodic::PeriodicInfo>& infoList = 
317
		dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
318

319
	    for (it = infoList.begin(); it != infoList.end(); ++it) {
320
321
	      if (it->elementSide == face) {
		for (int i = 0; i < dim; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
		  int i1 = faceIndEl_[i];
		  int i2 = faceIndNeigh_[i];
324

Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
		  int j = 0;
		  for (; j < dim; j++) {
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
		    if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, 
								   dim),
						      face,
						      j)) {
		      break;
		    }
		  }

		  TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
		      
		  neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
		}
		periodicCoords = true;
		break;
	      }
	    }
	  }
      
	  if (!periodicCoords) {
	    for (int i = 0; i < dim; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
	      int i1 = faceIndEl_[i];
	      int i2 = faceIndNeigh_[i];
	      for (int j = 0; j < dow; j++)
350
351
352
353
354
355
356
357
358
		neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
	    }
	  }
	      
	  Parametric *parametric = mesh->getParametric();
	  if (parametric) {
	    neighInfo = parametric->addParametricInfo(neighInfo);
	  }
	      
359
	  double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh_));
360
	      
361
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
362
	    jump_[iq].set(0.0);
363
364
365
	  }
	     

Thomas Witkowski's avatar
Thomas Witkowski committed
366
	  for (int system = 0; system < nSystems; system++) {	
367
368
369
	    if (matrix[system] == NULL) 
	      continue;
	      
370
371
	    uh[system]->getLocalVector(el, uhEl[system]);	
	    uh[system]->getLocalVector(neigh, uhNeigh[system]);
372
			
373
	    for (int iq = 0; iq < nPointsSurface_; iq++) {
374
	      (*lambda_)[face] = 0.0;
375
	      for (int i = 0; i < dim; i++) {
376
		(*lambda_)[faceIndEl_[i]] = surfaceQuad_->getLambda(iq, i);
377
378
	      }
		  
379
380
	      basFcts[system]->evalGrdUh(*lambda_, 
					 grdLambda, 
381
					 uhEl[system], 
382
					 &grdUhEl_[iq]);
383
		  
384
	      (*lambda_)[oppV] = 0.0;
385
	      for (int i = 0; i < dim; i++) {
386
		(*lambda_)[faceIndNeigh_[i]] = surfaceQuad_->getLambda(iq, i);
387
388
	      }
		  
389
390
	      basFcts[system]->evalGrdUh(*lambda_, 
					 *lambdaNeigh_, 
391
392
					 uhNeigh[system], 
					 &grdUhNeigh_[iq]);
393
		  
394
	      grdUhEl_[iq] -= grdUhNeigh_[iq];
395
396
	    }				

397
	    std::vector<double*>::iterator fac;
398
399
400
401
402

	    for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
		 ++it, ++fac) {
403

404
	      if (*fac == NULL || **fac != 0.0) {
405
406
407
408
409
410
411
412
413
414
415
416
417
418
		for (int iq = 0; iq < nPointsSurface_; iq++) {
		  localJump_[iq].set(0.0);
		}
		
		(*it)->weakEvalSecondOrder(nPointsSurface_,
					   grdUhEl_.getValArray(),
					   localJump_.getValArray());
		double factor = *fac ? **fac : 1.0;
		if (factor != 1.0) {
		  for (int i = 0; i < nPointsSurface_; i++) {
		    localJump_[i] *= factor;
		  }
		}
		
419
		for (int i = 0; i < nPointsSurface_; i++) {
420
		  jump_[i] += localJump_[i];
421
		}
422
423
	      }		
	    }
424
425
	  }
	      
426
	  val = 0.0;
427
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
428
	    val += surfaceQuad_->getWeight(iq) * (jump_[iq] * jump_[iq]);
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
	  }
	      
	  double d = 0.5 * (det + detNeigh);

	  if (norm == NO_NORM || norm == L2_NORM)
	    val *= C1 * h2_from_det(d, dim) * d;
	  else
	    val *= C1 * d;
	      
	  if (parametric) {
	    neighInfo = parametric->removeParametricInfo(neighInfo);
	  }

	  neigh->setEstimation(neigh->getEstimation(row) + val, row);
	  est_el += val;
	} 
      } 
       
447
      val = fh[std::max(row, 0)]->
448
	getBoundaryManager()->
449
	boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
450
451
452
453
454
455
456
457
458
459
460
461
462
463
      if (norm == NO_NORM || norm == L2_NORM)
	val *= C1 * h2;
      else
	val *= C1;
	
      est_el += val;
    } 
  

    el->setEstimation(est_el, row);

    est_sum += est_el;
    est_max = max(est_max, est_el);

464
    elInfo->getElement()->setMark(0);  
465
466
  }

467
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
468
	 int nPoints,
469
	 const double *uhIq,
470
471
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
472
	 const double *uhOldIq,
473
474
475
476
477
478
479
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
	 double *result)
  {
480
481
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
482
483
484
485
486
487

    // lhs
    for (it = const_cast<DOFMatrix*>(A)->getOperatorsBegin(),
	   fac = const_cast<DOFMatrix*>(A)->getOperatorEstFactorBegin(); 
	 it != const_cast<DOFMatrix*>(A)->getOperatorsEnd(); 
	 ++it, ++fac) {
488
489
490
     
      double factor = *fac ? **fac : 1.0;

491
492
      if (factor) {
	if (D2UhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
493
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
494
495
496
	}

	if (grdUhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
497
498
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
499
500
501
	}
	
	if (uhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
502
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
503
504
505
506
507
508
509
510
511
512
	}
      }
    }
    
    // rhs
    for (it = const_cast<DOFVector<double>*>(fh)->getOperatorsBegin(),
	 fac = const_cast<DOFVector<double>*>(fh)->getOperatorEstFactorBegin(); 
	 it != const_cast<DOFVector<double>*>(fh)->getOperatorsEnd(); 
	 ++it, ++fac) {

513
514
      double factor = *fac ? **fac : 1.0;

515
516
517
      if (factor) {
	if ((*it)->getUhOld()) {
	  if (D2UhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
518
	    (*it)->evalSecondOrder(nPoints, 
519
520
521
522
				   uhOldIq, grdUhOldIq, D2UhOldIq, 
				   result, factor);
	  }
	  if (grdUhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
523
	    (*it)->evalFirstOrderGrdPsi(nPoints, 
524
525
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
526
	    (*it)->evalFirstOrderGrdPhi(nPoints, 
527
528
529
530
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
	  }
	  if (uhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
531
	    (*it)->evalZeroOrder(nPoints, 
532
533
534
535
				 uhOldIq, grdUhOldIq, D2UhOldIq, 
				 result, -factor);
	  }
	} else {
536
	  std::vector<double> fx(nPoints, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
537
	  (*it)->getC(elInfo, nPoints, fx);
538

539
	  for (int iq = 0; iq < nPoints; iq++)
540
541
542
543
544
545
546
547
	    result[iq] -= factor * fx[iq];
	}
      }
    }    
  }


}