PollutionError.cc 4.28 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
#include "PollutionError.h"
#include "mpi.h"

namespace AMDiS
{

  PollutionErrorKonvex::BoundaryTraverseStack::BoundaryTraverseStack(Mesh *mesh, 
								     PartitionStatus ps_)
  {
    TraverseStack myStack;
    ElInfo *swap = 
      myStack.traverseFirst(mesh,-1, 
			    Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_NEIGH);
    ps = ps_;

    WorldVector<double>* centerWorld;
    DimVec<double> center = DimVec<double>(2, DEFAULT_VALUE, 1.0 / 3.0);
    while (swap != NULL) {
      if (proofConditions(swap)) {
	// if conditions are true, push the midpoint to the end of the list. a better 
	// algorithm should use the boundary edge
	centerWorld = new WorldVector<double>();
	swap->coordToWorld(center, *centerWorld);
	push_back(centerWorld);
      }
      swap = myStack.traverseNext(swap);
    }
  }
	
  bool AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::proofConditions(ElInfo *elInfo) 
  {
    if (elInfo == NULL) 
      return false;

    Element *curEl = elInfo->getElement();
    PartitionElementData *ped = 
      static_cast<PartitionElementData*>(curEl->getElementData(PARTITION_ED));

    TEST_EXIT_DBG(ped != NULL)("No PartitionElementData found");    

    if (ped != NULL && ped->getPartitionStatus() == OVERLAP) {
      //check neighbour
      int nrNeigh = curEl->getGeo(NEIGH);
      for (int i = 0; i < nrNeigh; i++) {
	Element* curNeigh = elInfo->getNeighbour(i);
	if (curNeigh != NULL) { 		    
	  ped = 
	    static_cast<PartitionElementData*>(curNeigh->getElementData(PARTITION_ED));
	  if (ped != NULL && ped->getPartitionStatus() == ps)
	    return true;
	}
      }
    }
    return false;
  }

  PollutionErrorKonvex::PollutionErrorKonvex(std::string name, int r)
    : ResidualEstimator(name,r),
      center(2,DEFAULT_VALUE, 1.0 / 3.0)
  {
    FUNCNAME("PollutionErrorKonvexEstimator::PollutionErrorKonvexEstimator()");
    GET_PARAMETER(0, name + "->C1p", "%f", &C1p);
    GET_PARAMETER(0, name + "->C2p", "%f", &C2p);   
  }

  void PollutionErrorKonvex::init(double timestep)  
  {
    FUNCNAME("PollutionErrorKonvex::init");
    ResidualEstimator::init(timestep);

    // set the distance (ie. width of the overlapped area)
    // approximate the overlapped area as rectangle. then set d=max(edge1, edge2)
    // retrieve all elements at the boundary overlapped -> out
    
    bts = new BoundaryTraverseStack(mesh, IN);	
  }

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

    // switch between inner and outer
    // inner part --> use H_1Norm
    // outer part --> use L_2Norm

    PartitionElementData* ped = 
      static_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));
    TEST_EXIT_DBG(ped != NULL)("No partition Elementdata found");

    PartitionStatus ps = ped->getPartitionStatus();
    TEST_EXIT_DBG(ps != UNDEFINED)("Elementstatus undefined");

    // save old values
    double C0 = ResidualEstimator::C0;
    double C1 = ResidualEstimator::C1;

    if (ps == OUT) {
      //outer part
      ResidualEstimator::norm = L2_NORM;		
      ResidualEstimator::C0 *= C2p;

      //compute distance
      double d = getDistance(elInfo);
      ResidualEstimator::C1 *= C2p / d;
      
    } else {
      //inner part
      ResidualEstimator::norm = H1_NORM;
      ResidualEstimator::C0 *= C1p;
      ResidualEstimator::C1 *= C1p;
    }

    ResidualEstimator::estimateElement(elInfo);
    //set old values
    ResidualEstimator::C0 = C0;
    ResidualEstimator::C1 = C1;
  }

  double PollutionErrorKonvex::getDistance(ElInfo* element)  
  {
    if (bts==NULL)
      MSG("bts is null");
    TEST_EXIT_DBG(bts!=NULL)("no bts");
    double distance = 0.0;
    double curDistance = 1.0;
    //get center in world coord
    WorldVector<double> centerWorld;
    element->coordToWorld(center,centerWorld);

    //start traversal over all Midpoints at the boundary overlap <-> in      
    for (BoundaryTraverseStack::iterator myIt = bts->begin(); 
	 myIt != bts->end(); myIt++) {
      curDistance = 0;	
      for (int i = 0;i < dim;i++)
	curDistance += ((*(*myIt))[i]-centerWorld[i]) * ((*(*myIt))[i] - centerWorld[i]);

      distance = max(curDistance, distance);
    }
    return distance;
  }
  
  void PollutionErrorKonvex::exit(bool output = true) 
  {
    ResidualEstimator::exit(output);
    delete bts;
  }
}