Estimator.cc 2.73 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
#include "Estimator.h"
#include "Traverse.h"
15
#include "Initfile.h"
16
#include "DualTraverse.h"
17
18
19

namespace AMDiS {

20
  Estimator::Estimator(std::string name_, int r) 
21
22
    : name(name_),
      norm(NO_NORM),
23
24
25
26
      row(r),
      mesh(NULL),
      auxMesh(NULL),
      traverseInfo(0)
27
  {
28
29
    FUNCNAME("Estimator::Estimator()");

30
	Parameters::get(name + "->error norm", norm);
31
32
  }

33

34
35
36
  double Estimator::estimate(double ts)
  {
    FUNCNAME("Estimator::estimate()");
37
38

    bool dualTraverse = false;
39
40

    /*
41
42
43
44
45
46
47
48
49
50
51
    for (unsigned int i = 0; i < matrix.size(); i++) {
      TEST_EXIT(traverseInfo.getStatus(row, i) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX)
	("Not yet implemented!\n");

      if (traverseInfo.getStatus(row, i) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX ||
	  traverseInfo.getStatus(row, i) == SingleComponentInfo::DIF_SPACES_NO_AUX ||
	  traverseInfo.getStatus(row, i) == SingleComponentInfo::DIF_SPACES_WITH_AUX)
	dualTraverse = true;
    }

    if (!dualTraverse) {
52
      mesh = uh[row == -1 ? 0 : row]->getFeSpace()->getMesh();
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
      auxMesh = NULL;
    } else {
      const FiniteElemSpace *mainFeSpace = traverseInfo.getRowFeSpace(row);
      const FiniteElemSpace *auxFeSpace = traverseInfo.getNonRowFeSpace(row);

      TEST_EXIT(mainFeSpace)("No main FE space!\n");
      TEST_EXIT(auxFeSpace)("No aux FE space!\n"); 

      mesh = mainFeSpace->getMesh();
      auxMesh = auxFeSpace->getMesh();

      TEST_EXIT_DBG(mainFeSpace->getBasisFcts()->getDegree() ==
		    auxFeSpace->getBasisFcts()->getDegree())
	("Mh, do you really want to do this? Think about it ...\n");
    }
68
69
    */

70
    mesh = uh[row == -1 ? 0 : row]->getFeSpace()->getMesh();
71
72
    auxMesh = NULL;

73
    init(ts);
74

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
    if (!dualTraverse)
      singleMeshTraverse();
    else
      dualMeshTraverse();

    exit();

    return est_sum;
  }


  void Estimator::singleMeshTraverse()
  {
    FUNCNAME("Estimator::singleMeshTraverse()");

90
91
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
92
    while (elInfo) {
93
94
      estimateElement(elInfo);
      elInfo = stack.traverseNext(elInfo);
95
    }  
96
  }
97
98


99
100
101
102
103
104
105
106
107
108
109
110
111
112
  void Estimator::dualMeshTraverse()
  {
    FUNCNAME("Estimator::dualMeshTraverse()");

    DualTraverse dualTraverse;
    DualElInfo dualElInfo;

    bool cont = dualTraverse.traverseFirst(mesh, auxMesh, -1, -1, 
					   traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
      estimateElement(dualElInfo.rowElInfo, &dualElInfo);      
      cont = dualTraverse.traverseNext(dualElInfo);
    }
113
114
  }
}