-
Notifications
You must be signed in to change notification settings - Fork 0
/
filterstruct.h
109 lines (100 loc) · 4.66 KB
/
filterstruct.h
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
#include <vector>
#include <utility>
#include "pzreal.h"
#include "pzvec_extras.h"
struct FilterStruct {
std::vector<std::pair<int, REAL>> fneighIndexHf;
REAL fsumHf;
int findex;
FilterStruct(int index) : findex(index), fsumHf(0.) {
}
FilterStruct() : findex(-1), fsumHf(0.) {
}
REAL ComputeFiltereddcdxe(TPZVec<REAL> &dcdxf, TPZVec<REAL>& xf, TPZCompMesh* cmesh) {
// TPZGeoMesh* gmesh = cmesh->Reference();
REAL xe = xf[findex];
REAL sum = 0.;
if (fneighIndexHf.size() == 0) {
return dcdxf[findex];
}
for (const auto& neighIndexHf : fneighIndexHf) {
int64_t neighIndex = neighIndexHf.first;
REAL Hfneigh = neighIndexHf.second;
if (Hfneigh < 0) DebugStop(); // what happened? It is defined as rmin - distcenters and should only be here if positeve
TPZCompEl* cel = cmesh->Element(neighIndex);
if(!cel) continue;
const int64_t celindex = cel->Index(); // geoel and compel indexes are not the same
REAL xneigh = xf[celindex];
REAL dcdxneigh = dcdxf[celindex];
sum += Hfneigh * xneigh * dcdxneigh;
}
sum /= (fsumHf*xe);
return sum;
}
void Print(std::ostream &out) {
out << "Filter index " << findex << " sumHf " << fsumHf << " neighIndexHf " << " nneigh " << fneighIndexHf.size() << std::endl;
for (const auto& neighIndexHf : fneighIndexHf) {
out << neighIndexHf.first << " " << neighIndexHf.second << "\n";
}
out << std::endl;
}
void ComputeNeighIndexHf(TPZGeoMesh* gmesh, REAL rmin) {
TPZCompMesh* cmesh = gmesh->Reference();
TPZCompEl* cel = cmesh->Element(findex);
TPZGeoEl* geo = cel->Reference();
TPZGeoElSide geoside(geo);
TPZManVector<REAL,3> rootcenter(3,0.);
geoside.CenterX(rootcenter);
TPZManVector<REAL,3> centerneigh(3,0.);
int ncornernodes = geo->NCornerNodes();
std::set<int64_t> tocheck, included, checked;
tocheck.insert(findex);
while (tocheck.size()) {
TPZCompEl* celcheck = cmesh->Element(*tocheck.begin());
TPZGeoEl* gel = celcheck->Reference();
tocheck.erase(tocheck.begin());
TPZGeoElSide gelside(gel);
gelside.CenterX(centerneigh);
REAL d = dist(centerneigh,rootcenter);
REAL Hf = rmin - d;
// std::cout << "Hf " << Hf << " d " << d << " rmin " << rmin << std::endl;
if (Hf > 0) {
included.insert(celcheck->Index());
fneighIndexHf.push_back(std::make_pair(celcheck->Index(),Hf));
fsumHf += Hf;
for (int iside = 0; iside < gel->NCornerNodes(); iside++) {
TPZGeoElSide nodeside(gel,iside);
TPZGeoElSide neighbour = nodeside.Neighbour();
if (neighbour.Element()->HasSubElement()) {
TPZStack<TPZGeoEl*> subels;
// Get all the youngest children and insert them in the tocheck set
neighbour.Element()->YoungestChildren(subels);
for (auto subel : subels) {
if (subel->Dimension() == gmesh->Dimension()) {
int64_t neighindex = subel->Reference()->Index();
if (checked.find(neighindex) == checked.end() && included.find(neighindex) == included.end()) {
tocheck.insert(neighindex);
}
}
}
continue;
}
while(neighbour != nodeside) {
if (neighbour.Element()->Dimension() == gmesh->Dimension() && neighbour.Element()->Reference()) {
int64_t neighindex = neighbour.Element()->Reference()->Index();
// std::cout << "Verifying neighbour " << neighindex << std::endl;
if (checked.find(neighindex) == checked.end() && included.find(neighindex) == included.end()) {
// std::cout << "Inserting neighbour " << neighindex << std::endl;
tocheck.insert(neighindex);
}
}
neighbour = neighbour.Neighbour();
}
}
}
else {
checked.insert(gel->Reference()->Index());
}
}
}
};