ResidualEstimator.cc 15.9 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
	if (*itfac == NULL || **itfac != 0.0)
219
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
220

221
      if (C0 > 0.0)
222
	for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it)
223
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);	
224
225
226
227
228
229
230
231
232
233
234
235
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
    }


    // === 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;

262
      if (timestep && uhOld[system]) {
263
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
264
	uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
265
  
266
267
268
	// === Compute time error. ===

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

311
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
312
	  nPoints, 
313
314
315
316
317
318
	  uhQP,
	  grdUh_qp,
	  D2uhqp,
	  uhOldQP,
	  NULL,  // grdUhOldQP 
	  NULL,  // D2UhOldQP
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
	  dofMat, 
	  dofVec,
321
	  quad,
322
	  riq);
323
324
325
326
      }     
    }

    // add integral over r square
327
    double result = 0.0;
328
    for (int iq = 0; iq < nPoints; iq++)
329
      result += quad->getWeight(iq) * riq[iq] * riq[iq];
330
   
331
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
332
      result = C0 * h2 * h2 * det * result;
333
    else
334
      result = C0 * h2 * det * result;
335
    
336
337
    return result;
  }
338

339

340
341
342
343
344
345
346
347
348
349
350
  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);
351

352
353
    for (int face = 0; face < nNeighbours; face++) {  
      Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
354

355
356
      if (!(neigh && neigh->getMark()))
	continue;
357

358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
      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) {
378
	    for (int i = 0; i < dim; i++) {
379
380
	      int i1 = faceIndEl[i];
	      int i2 = faceIndNeigh[i];
381
	      
382
383
384
385
	      int j = 0;
	      for (; j < dim; j++)
		if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), face, j))
		  break;
386
	      
387
	      TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
388
	      
389
	      neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
390
	    }
391
392
	    periodicCoords = true;
	    break;
393
	  }
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
	}
      }  // 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);	  
409

410
411
412
413
414
415
416
417
      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;
418
	      
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
	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);
	    
	    (*it)->weakEvalSecondOrder(nPointsSurface, &(grdUhEl[0]), &(localJump[0]));
449

450
451
452
453
454
455
456
457
458
459
460
461
462
463
	    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
464

465
466
      double d = 0.5 * (det + detNeigh);
   
467
      if (norm == NO_NORM || norm == L2_NORM)
468
	val *= C1 * h2_from_det(d, dim) * d;
469
      else
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
	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;
486

487
    return result;
488
489
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
490

491
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
492
	 int nPoints,
493
	 const double *uhIq,
494
495
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
496
	 const double *uhOldIq,
497
498
499
500
501
502
503
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
	 double *result)
  {
504
505
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
506
507

    // lhs
Thomas Witkowski's avatar
Thomas Witkowski committed
508
509
    for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin(); 
	 it != A->getOperatorsEnd(); ++it, ++fac) {
510
511
512
     
      double factor = *fac ? **fac : 1.0;

513
      if (factor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
514
515
	if (D2UhIq)
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);	
516
517

	if (grdUhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
520
521
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
522
523
	if (uhIq)
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);	
524
525
526
527
      }
    }
    
    // rhs
Thomas Witkowski's avatar
Thomas Witkowski committed
528
529
    for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin(); 
	 it != fh->getOperatorsEnd(); ++it, ++fac) {
530

531
532
      double factor = *fac ? **fac : 1.0;

533
534
      if (factor) {
	if ((*it)->getUhOld()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
535
536
537
	  if (D2UhOldIq)
	    (*it)->evalSecondOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, factor);
	  
538
	  if (grdUhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
539
540
	    (*it)->evalFirstOrderGrdPsi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
	    (*it)->evalFirstOrderGrdPhi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
541
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
542
543
544

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

549
	  for (int iq = 0; iq < nPoints; iq++)
550
551
552
553
554
555
556
557
	    result[iq] -= factor * fx[iq];
	}
      }
    }    
  }


}