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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195 | #include "IsoLaplace.hpp"
#include <algorithm>
//---------------------------------------------------------------------------//
// construction function for IsoLaplace class
IsoLaplace::IsoLaplace()
{
}
//---------------------------------------------------------------------------//
// deconstruction function for IsoLaplace class
IsoLaplace::~IsoLaplace()
{
}
//setup the data for IsoLaplace
void IsoLaplace::SetupData(std::vector<std::set<int> > AdjElements, std::vector<std::vector<int> > Quads, std::vector<std::vector<double> > coords, std::vector<bool> isBnd, std::vector<double> w)
{
//pass the node coordinates to the local variables
coordinates.resize(coords.size());
for (unsigned int i = 0; i < coords.size(); i++){
coordinates[i].resize(coords[i].size());
for (unsigned int j = 0; j < coords[i].size(); j++)
coordinates[i][j] = coords[i][j];
}
//pass the connectivity information to the local variables
adjElements.resize(AdjElements.size());
for (unsigned int i = 0; i < AdjElements.size(); i++){
for (std::set<int>::iterator it = AdjElements[i].begin(); it != AdjElements[i].end(); it++)<--- Prefer prefix ++/-- operators for non-primitive types.
adjElements[i].insert(*it);
}
quads.resize(Quads.size());
for (unsigned int i = 0; i < Quads.size(); i++){
quads[i].resize(Quads[i].size());
for (unsigned int j = 0; j < Quads[i].size(); j++)
quads[i][j] = Quads[i][j];
}
//pass the boundary info to the local variables
isBoundary.resize(isBnd.size());
for (unsigned int i = 0; i < isBnd.size(); i++)
isBoundary[i] = isBnd[i];
//pass the element's weight to the local variables
weight.resize(w.size());
for (unsigned int i = 0; i < w.size(); i++)
weight[i] = w[i];
//done
}
//return the results
void IsoLaplace::GetCoords(std::vector<std::vector<double> > &coords)
{
//return the results to the user
coords.resize(coordinates.size());
for (unsigned int i = 0; i < coords.size(); i++){
if (!isBoundary[i]){
coords[i].resize(coordinates[i].size());
for (unsigned int j = 0; j < coordinates[i].size(); j++)
coords[i][j] = coordinates[i][j];
}
}
}
//Execute function
void IsoLaplace::Execute()
{
//test the input
for (unsigned int i = 0; i < coordinates.size(); i++){
if (!isBoundary[i]){
std::cout << "IsoLaplace index = " << i << "\t[" << coordinates[i][0] << ", " << coordinates[i][1] << ", " << coordinates[i][2] << "]\n";
std::cout << "\n\n";
}
}
std::cout << "test the adjacent element info\n";
for (unsigned int i = 0; i < adjElements.size(); i++){
std::set<int>::iterator it = adjElements[i].begin();
std::cout << "node index = " << i << "\t";
for (; it != adjElements[i].end(); it++)<--- Prefer prefix ++/-- operators for non-primitive types.
std::cout << *it << "\t";
std::cout << "\n";
}
std::cout << "test the quad info\n";
for (unsigned int i = 0; i < quads.size(); i++){
std::cout << "quad index = " << i << "\nthe adjacent nodes are ";
for (unsigned int j = 0; j < quads[i].size(); j++){
std::cout << quads[i][j] << "\t";
}
std::cout << std::endl;
}
std::cout << "test the weights\n";
for (unsigned int i = 0; i < weight.size(); i++){
std::cout << "index = " << i << "\tweight is " << weight[i] << std::endl;
}
double epilson = 0.01;
double error = 0.0;<--- The scope of the variable 'error' can be reduced.
int step = 0;
double r = 0.1;
while(true){
error = 0.0;
//loop over all the nodes
for (unsigned int i = 0; i < coordinates.size(); i++){
if (!(isBoundary[i])){
double sum_neighbors[3] = {0.0, 0.0, 0.0};
double sum_opposite[3] = {0.0, 0.0, 0.0};
double sum_weight = 0.0;
//loop over the adjacent quads
for (std::set<int>::iterator it = adjElements[i].begin(); it != adjElements[i].end(); it++){<--- Prefer prefix ++/-- operators for non-primitive types.
int stNodeIndex = -1;
if (int(i) == quads[*it][0])
stNodeIndex = 0;
else if (int(i) == quads[*it][1])
stNodeIndex = 1;
else if (int(i) == quads[*it][2])
stNodeIndex = 2;
else
stNodeIndex = 3;
//sum up the neighbor nodes and opposite nodes
for (int k = 0; k < 3; k++){
sum_neighbors[k] += weight[*it]*coordinates[quads[*it][(stNodeIndex+1)%4]][k];
sum_neighbors[k] += weight[*it]*coordinates[quads[*it][(stNodeIndex+3)%4]][k];
sum_opposite[k] += weight[*it]*r*coordinates[quads[*it][(stNodeIndex+2)%4]][k];
}
sum_weight += weight[*it];
}
double tmp[3];
double e = 0.0;
//average the sum_neighbors and sum_opposite
for (int k = 0; k < 3; k++){
sum_neighbors[k] = sum_neighbors[k]/sum_weight;
sum_opposite[k] = sum_opposite[k]/sum_weight;
}
//update the node position and calculate the error
for (int j = 0; j < 3; j++){
tmp[j] = (sum_neighbors[j] - sum_opposite[j])/(2-r);
e += fabs(tmp[j] - coordinates[i][j]);
coordinates[i][j] = tmp[j];
}
if(e > error)
error = e;
}
}
step++;
//UpdateWeight();
std::cout << "IsoLaplace smoothing step = " << step << "\tError = " << error << std::endl;
if (error < epilson) break;
}
for (unsigned int i = 0; i < coordinates.size(); i++){
if (!isBoundary[i]){
std::cout << "IsoLaplace index = " << i << "\t[" << coordinates[i][0] << ", " << coordinates[i][1] << ", " << coordinates[i][2] << "]\n";
std::cout << "\n\n";
}
}
std::cout << std::endl;
}
void IsoLaplace::UpdateWeight(){<--- The function 'UpdateWeight' is never used.
for (unsigned int i = 0; i < quads.size(); i++){
//calculate the area
double a, b, c, s;
weight[i] = 0.0;
//calculate one half triangle
a = sqrt(pow(coordinates[quads[i][0]][0] - coordinates[quads[i][1]][0], 2) + pow(coordinates[quads[i][0]][1] - coordinates[quads[i][1]][1],2) + pow(coordinates[quads[i][0]][2] - coordinates[quads[i][1]][2], 2));
b = sqrt(pow(coordinates[quads[i][0]][0] - coordinates[quads[i][2]][0], 2) + pow(coordinates[quads[i][0]][1] - coordinates[quads[i][2]][1],2) + pow(coordinates[quads[i][0]][2] - coordinates[quads[i][2]][2], 2));
c = sqrt(pow(coordinates[quads[i][1]][0] - coordinates[quads[i][2]][0], 2) + pow(coordinates[quads[i][1]][1] - coordinates[quads[i][2]][1],2) + pow(coordinates[quads[i][1]][2] - coordinates[quads[i][2]][2], 2));
s = 0.5*(a + b + c);
weight[i] += sqrt(fabs(s*(s-a)*(s-b)*(s-c)));
//calculate the other half triangle
a = sqrt(pow(coordinates[quads[i][0]][0] - coordinates[quads[i][3]][0], 2) + pow(coordinates[quads[i][0]][1] - coordinates[quads[i][3]][1],2) + pow(coordinates[quads[i][0]][2] - coordinates[quads[i][3]][2], 2));
b = sqrt(pow(coordinates[quads[i][3]][0] - coordinates[quads[i][2]][0], 2) + pow(coordinates[quads[i][3]][1] - coordinates[quads[i][2]][1],2) + pow(coordinates[quads[i][3]][2] - coordinates[quads[i][2]][2], 2));
c = sqrt(pow(coordinates[quads[i][1]][0] - coordinates[quads[i][2]][0], 2) + pow(coordinates[quads[i][1]][1] - coordinates[quads[i][2]][1],2) + pow(coordinates[quads[i][1]][2] - coordinates[quads[i][2]][2], 2));
s = 0.5*(a + b + c);
weight[i] += sqrt(fabs(s*(s-a)*(s-b)*(s-c)));
}
}
|