ResidualEstimator.cc 14.7 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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
    : 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()");

    timestep = ts;

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

    numSystems = static_cast<int>(uh.size());
    TEST_EXIT_DBG(numSystems > 0)("no system set\n");

    dim = mesh->getDim();
    basFcts = GET_MEMORY(const BasisFunction*, numSystems);
    quadFast = GET_MEMORY(FastQuadrature*, numSystems);

    degree = 0;
    for (int system = 0; system < numSystems; system++) {
      basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
47
      degree = std::max(degree, basFcts[system]->getDegree());
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
    }

    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
    numPoints = quad->getNumPoints();

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

    for (int system = 0; system < numSystems; system++) {
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
							       *quad, 
							       flag);
    }
  
    uhEl = GET_MEMORY(double*, numSystems);
67
    uhNeigh = GET_MEMORY(double*, numSystems);
68
69
70
71
    uhOldEl = timestep ? GET_MEMORY(double*, numSystems) : NULL;

    for (int system = 0; system < numSystems; system++) {
      uhEl[system] = GET_MEMORY(double, basFcts[system]->getNumber()); 
72
      uhNeigh[system] = GET_MEMORY(double, basFcts[system]->getNumber());
73
74
75
76
77
78
79
80
81
      if (timestep)
	uhOldEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
    }

    uhQP = timestep ? GET_MEMORY(double, numPoints) : NULL;
    uhOldQP = timestep ? GET_MEMORY(double, numPoints) : NULL;

    riq = GET_MEMORY(double, numPoints);

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
133
134
  }

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

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

    for (int system = 0; system < numSystems; system++) {
      FREE_MEMORY(uhEl[system], double, basFcts[system]->getNumber());
135
      FREE_MEMORY(uhNeigh[system], double, basFcts[system]->getNumber());
136
137
138
139
140
      if (timestep)
	FREE_MEMORY(uhOldEl[system], double, basFcts[system]->getNumber());    
    }

    FREE_MEMORY(uhEl, double*, numSystems);
141
    FREE_MEMORY(uhNeigh, double*, numSystems);
142
143

    if (timestep) {
144
      FREE_MEMORY(uhOldEl, double*, numSystems);
145
146
      FREE_MEMORY(uhQP, double, numPoints);
      FREE_MEMORY(uhOldQP, double, numPoints);
147
148
149
150
    } else {
      if (uhQP != NULL) {
	FREE_MEMORY(uhQP, double, numPoints);
      }
151
152
153
154
155
156
157
158
159
160
161
162
    }

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

    FREE_MEMORY(riq, double, numPoints);
    FREE_MEMORY(basFcts, const BasisFunction*, numSystems);
    FREE_MEMORY(quadFast, FastQuadrature*, numSystems);
163
164
165
166
167
168
169

    if (grdUh_qp != NULL) {
      FREE_MEMORY(grdUh_qp, WorldVector<double>, numPoints);
    }
    if (D2uhqp != NULL) {
      FREE_MEMORY(D2uhqp, WorldMatrix<double>, numPoints);
    }
170

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

176
    DELETE neighInfo;
177
178
179
180
181
182
183
184
185
186
187
188
  }

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

    double est_el = 0.0;
    double val = 0.0;
    Element *el, *neigh;

    TEST_EXIT_DBG(numSystems > 0)("no system set\n");

189
    std::vector<Operator*>::iterator it;
190
191
192
193

    el = elInfo->getElement();

    double det = elInfo->getDet();
194
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
195
196
197
198
199

    est_el = el->getEstimation(row);

    double h2 = h2_from_det(det, dim);

200
    for (int iq = 0; iq < numPoints; iq++) {
201
202
203
204
205
206
207
208
209
      riq[iq] = 0.0;
    }

    for (int system = 0; system < numSystems; system++) {

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

      // init assemblers
210
      std::vector<Operator*>::iterator it;
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232

      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
	   it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
	   ++it) {
	(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, quad);
      }

      for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
	   it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd(); 
	   ++it) {
	(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, quad);
      }

      if (timestep) {
	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);
	  
233
234
235
	  if (C3 && uhOldQP && system == std::max(row, 0)) {
	    val = 0.0;
	    for (int iq = 0; iq < numPoints; iq++) {
236
237
238
239
240
241
242
243
244
	      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);
	  }
	}
      }
245
           
246
247
248
249
      if (C0) {  
	for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(); 
	     it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
	     ++it) {
250
	  if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
251
252
253
	    uhQP = GET_MEMORY(double, numPoints);
	    uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  }
254
	  if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
255
256
257
	    grdUh_qp = new WorldVector<double>[numPoints];
	    uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
	  }
258
259
260
	  if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) { 
	    D2uhqp = new WorldMatrix<double>[numPoints];
	    uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
261
262
263
264
265
266
267
268
269
270
271
272
273
274
	  }
	}
	
	r(elInfo,
	  numPoints, 
	  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
281
    val = 0.0;
    for (int iq = 0; iq < numPoints; iq++) {
282
283
      val += quad->getWeight(iq) * riq[iq] * riq[iq];
    }
284
285
   
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM) {
286
      val = C0 * h2 * h2 * det * val;
287
    } else {
288
      val = C0 * h2 * det * val;
289
    }
290
291
292
293
294

    est_el += val;

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

297
      for (int face = 0; face < neighbours_; face++) {  
298
299
300
301
	neigh = const_cast<Element*>(elInfo->getNeighbour(face));
	if (neigh && neigh->getMark()) {      
	  int oppV = elInfo->getOppVertex(face);
	      
302
303
	  el->sortFaceIndices(face, &faceIndEl_);
	  neigh->sortFaceIndices(oppV, &faceIndNeigh_);
304
305
306
	    
	  neighInfo->setElement(const_cast<Element*>(neigh));
	  neighInfo->setFillFlag(Mesh::FILL_COORDS);
307
	      	
308
309
310
311
312
313
314
315
316
317
318
	  int j, i1, i2;

	  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) {
319
320
	    std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	    std::list<LeafDataPeriodic::PeriodicInfo>& infoList = 
321
		dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
322

323
	    for (it = infoList.begin(); it != infoList.end(); ++it) {
324
325
	      if (it->elementSide == face) {
		for (int i = 0; i < dim; i++) {
326
327
		  i1 = faceIndEl_[i];
		  i2 = faceIndNeigh_[i];
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349

		  for (j = 0; j < dim; j++) {
		    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++) {
350
351
	      i1 = faceIndEl_[i];
	      i2 = faceIndNeigh_[i];
352
353
354
355
356
357
358
359
360
361
	      for (j = 0; j < dow; j++)
		neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
	    }
	  }
	      
	  Parametric *parametric = mesh->getParametric();
	  if (parametric) {
	    neighInfo = parametric->addParametricInfo(neighInfo);
	  }
	      
362
	  double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh_));
363
	      
364
	  for (int iq = 0; iq < nPointsSurface_; iq++) {
365
	    jump_[iq].set(0.0);
366
367
368
369
370
371
372
	  }
	     

	  for (int system = 0; system < numSystems; system++) {	
	    if (matrix[system] == NULL) 
	      continue;
	      
373
374
	    uh[system]->getLocalVector(el, uhEl[system]);	
	    uh[system]->getLocalVector(neigh, uhNeigh[system]);
375
			
376
	    for (int iq = 0; iq < nPointsSurface_; iq++) {
377
	      (*lambda_)[face] = 0.0;
378
	      for (int i = 0; i < dim; i++) {
379
		(*lambda_)[faceIndEl_[i]] = surfaceQuad_->getLambda(iq, i);
380
381
	      }
		  
382
383
	      basFcts[system]->evalGrdUh(*lambda_, 
					 grdLambda, 
384
					 uhEl[system], 
385
					 &grdUhEl_[iq]);
386
		  
387
	      (*lambda_)[oppV] = 0.0;
388
	      for (int i = 0; i < dim; i++) {
389
		(*lambda_)[faceIndNeigh_[i]] = surfaceQuad_->getLambda(iq, i);
390
391
	      }
		  
392
393
	      basFcts[system]->evalGrdUh(*lambda_, 
					 *lambdaNeigh_, 
394
395
					 uhNeigh[system], 
					 &grdUhNeigh_[iq]);
396
		  
397
	      grdUhEl_[iq] -= grdUhNeigh_[iq];
398
399
	    }				

400
	    std::vector<double*>::iterator fac;
401
402
403
404
405

	    for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
		 ++it, ++fac) {
406
	      for (int iq = 0; iq < nPointsSurface_; iq++) {
407
		localJump_[iq].set(0.0);
408
409
	      }
		  
410
411
412
	      (*it)->weakEvalSecondOrder(nPointsSurface_,
					 grdUhEl_.getValArray(),
					 localJump_.getValArray());
413
414
	      double factor = *fac ? **fac : 1.0;
	      if (factor != 1.0) {
415
416
		for (int i = 0; i < nPointsSurface_; i++) {
		  localJump_[i] *= factor;
417
418
419
		}
	      }
		  
420
421
	      for (int i = 0; i < nPointsSurface_; i++) {
		jump_[i] += localJump_[i];
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
464
465
466
      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);

    elInfo->getElement()->setMark(0);   
  }

467
468
469
  void r(const ElInfo *elInfo,
	 int numPoints,
	 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
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
      if (factor) {
	if (D2UhIq) {
	  (*it)->evalSecondOrder(numPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
	}

	if (grdUhIq) {
	  (*it)->evalFirstOrderGrdPsi(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	}
	
	if (uhIq) {
	  (*it)->evalZeroOrder(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	}
      }
    }
    
    // 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
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
      if (factor) {
	if ((*it)->getUhOld()) {
	  if (D2UhOldIq) {
	    (*it)->evalSecondOrder(numPoints, 
				   uhOldIq, grdUhOldIq, D2UhOldIq, 
				   result, factor);
	  }
	  if (grdUhOldIq) {
	    (*it)->evalFirstOrderGrdPsi(numPoints, 
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
	    (*it)->evalFirstOrderGrdPhi(numPoints, 
					uhOldIq, grdUhOldIq, D2UhOldIq, 
					result, -factor);
	  }
	  if (uhOldIq) {
	    (*it)->evalZeroOrder(numPoints, 
				 uhOldIq, grdUhOldIq, D2UhOldIq, 
				 result, -factor);
	  }
	} else {
	  double *fx = GET_MEMORY(double, numPoints);
	  for (int iq = 0; iq < numPoints; iq++) {
	    fx[iq] = 0.0;
	  }
	  (*it)->getC(elInfo, numPoints, fx);

	  for (int iq = 0; iq < numPoints; iq++) {
	    result[iq] -= factor * fx[iq];
	  }
	  FREE_MEMORY(fx, double, numPoints);
	}
      }
    }    
  }


}