Branch data Line data Source code
1 : : #include "EquipotentialSmooth.hpp"
2 : : #include <algorithm>
3 : :
4 : :
5 : : //---------------------------------------------------------------------------//
6 : : // construction function for EquipotentialSmooth class
7 [ # # ][ # # ]: 0 : EquipotentialSmooth::EquipotentialSmooth()
[ # # ][ # # ]
8 : : {
9 : :
10 : :
11 : 0 : }
12 : :
13 : :
14 : :
15 : : //---------------------------------------------------------------------------//
16 : : // deconstruction function for EquipotentialSmooth class
17 : 0 : EquipotentialSmooth::~EquipotentialSmooth()
18 : : {
19 : :
20 : 0 : }
21 : :
22 : : //setup the data for EquipotentialSmooth
23 : 0 : void EquipotentialSmooth::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<std::vector<int> > conn)
24 : : {
25 : : //pass the node coordinates to the local variables
26 : 0 : coordinates.resize(coords.size());
27 [ # # ]: 0 : for (unsigned int i = 0; i < coords.size(); i++){
28 : 0 : coordinates[i].resize(coords[i].size());
29 [ # # ]: 0 : for (unsigned int j = 0; j < coords[i].size(); j++)
30 : 0 : coordinates[i][j] = coords[i][j];
31 : : }
32 : :
33 : : //pass the connectivity information to the local variables
34 : 0 : adjElements.resize(AdjElements.size());
35 [ # # ]: 0 : for (unsigned int i = 0; i < AdjElements.size(); i++){
36 [ # # ][ # # ]: 0 : for (std::set<int>::iterator it = AdjElements[i].begin(); it != AdjElements[i].end(); it++)
[ # # ][ # # ]
[ # # ]
37 [ # # ][ # # ]: 0 : adjElements[i].insert(*it);
[ # # ]
38 : : }
39 : :
40 : 0 : quads.resize(Quads.size());
41 [ # # ]: 0 : for (unsigned int i = 0; i < Quads.size(); i++){
42 : 0 : quads[i].resize(Quads[i].size());
43 [ # # ]: 0 : for (unsigned int j = 0; j < Quads[i].size(); j++)
44 : 0 : quads[i][j] = Quads[i][j];
45 : : }
46 : :
47 : : //pass the boundary info to the local variables
48 : 0 : isBoundary.resize(isBnd.size());
49 [ # # ]: 0 : for (unsigned int i = 0; i < isBnd.size(); i++)
50 [ # # ]: 0 : isBoundary[i] = isBnd[i];
51 : :
52 : : //pass the connectivity to the local variables
53 : 0 : connect.resize(conn.size());
54 [ # # ]: 0 : for (unsigned int i = 0; i < conn.size(); i++){
55 : 0 : connect[i].resize(conn[i].size());
56 [ # # ]: 0 : for (unsigned int j = 0; j < conn[i].size(); j++)
57 : 0 : connect[i][j] = conn[i][j];
58 : : }
59 : : //done
60 : 0 : }
61 : :
62 : : //return the results
63 : 0 : void EquipotentialSmooth::GetCoords(std::vector<std::vector<double> > &coords)
64 : : {
65 : : //return the results to the user
66 : 0 : coords.resize(coordinates.size());
67 [ # # ]: 0 : for (unsigned int i = 0; i < coords.size(); i++){
68 [ # # ]: 0 : if (!isBoundary[i]){
69 : 0 : coords[i].resize(coordinates[i].size());
70 [ # # ]: 0 : for (unsigned int j = 0; j < coordinates[i].size(); j++)
71 : 0 : coords[i][j] = coordinates[i][j];
72 : : }
73 : : }
74 : 0 : }
75 : :
76 : : //Execute function
77 : 0 : void EquipotentialSmooth::Execute()
78 : : {
79 : : /* Algorithm --- Iterative smoothing */
80 : : /* p'(i) = p(i) + sum{w(n)*(p(n)-p(i)), n = 1...N} elements 1...N are the neighbors of node i */
81 : : /* w(1) = -belta/2 w(2) = alpha w(3) = belta/2 w(4) = gamma */
82 : : /* w(5) = -belta/2 w(6) = alpha w(7) = belta/2 w(8) = gamma */
83 : : /* where alpha = xp^2 + yp^2 + zp^2 */
84 : : /* belta = xp*xq + yp*yq + zp*zq */
85 : : /* gamma = xq^2 + yq^2 + zq^2 */
86 : : /* xp = (x2-x6)/2 yp = (y2-y6)/2 zp = (z2-z6)/2 */
87 : : /* xq = (x8-x4)/2 yq = (y8-y4)/2 zq = (z8-z4)/2 */
88 : :
89 : : // 3------------2------------1
90 : : // | | |
91 : : // | | |
92 : : // | | |
93 : : // 4------------i------------8
94 : : // | | |
95 : : // | | |
96 : : // | | |
97 : : // 5------------6------------7
98 : :
99 : 0 : double epilson = 0.01;
100 : 0 : double error = 0.0;
101 : 0 : int step = 0;
102 : : while(true){
103 : 0 : error = 0.0;
104 : : //loop over all the nodes
105 [ # # ]: 0 : for (unsigned int i = 0; i < coordinates.size(); i++){
106 [ # # ]: 0 : if (!(isBoundary[i])){
107 : : double weight[9];
108 : 0 : double alpha = 0.0, belta = 0.0, gamma = 0.0;
109 : 0 : double p[3] = {0.0, 0.0, 0.0};
110 : 0 : double q[3] = {0.0, 0.0, 0.0};
111 [ # # ]: 0 : for (int j = 0; j < 3; j++){
112 [ # # ][ # # ]: 0 : p[j] = (coordinates[connect[i][2]][j] - coordinates[connect[i][6]][j])/2.0;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
113 [ # # ][ # # ]: 0 : q[j] = (coordinates[connect[i][4]][j] - coordinates[connect[i][8]][j])/2.0;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
114 : : }
115 [ # # ][ # # ]: 0 : alpha = pow(p[0],2) + pow(p[1],2) + pow(p[2],2);
[ # # ]
116 [ # # ][ # # ]: 0 : gamma = pow(q[0],2) + pow(q[1],2) + pow(q[2],2);
[ # # ]
117 : 0 : belta = p[0]*q[0] + p[1]*q[1] + p[2]*q[2];
118 : :
119 : 0 : weight[1] = belta/2.0;
120 : 0 : weight[2] = gamma;
121 : 0 : weight[3] = -belta/2.0;
122 : 0 : weight[4] = alpha;
123 : 0 : weight[5] = belta/2.0;
124 : 0 : weight[6] = gamma;
125 : 0 : weight[7] = -belta/2.0;
126 : 0 : weight[8] = alpha;
127 : :
128 : 0 : double sum_dis[3] = {0.0, 0.0, 0.0};
129 [ # # ]: 0 : for (int j = 0; j < 3; j++){
130 [ # # ]: 0 : for (int k = 1; k < 9; k++){
131 [ # # ][ # # ]: 0 : sum_dis[j]+=weight[k]*coordinates[connect[i][k]][j]/(2.0*(alpha+gamma));
[ # # ][ # # ]
132 : : }
133 : : }
134 [ # # ][ # # ]: 0 : double e = fabs(sum_dis[0] - coordinates[i][0]) + fabs(sum_dis[1] - coordinates[i][1]) + fabs(sum_dis[2]- coordinates[i][2]);
[ # # ][ # # ]
[ # # ][ # # ]
135 [ # # ]: 0 : for (int j = 0; j < 3; j++)
136 [ # # ][ # # ]: 0 : coordinates[i][j] = sum_dis[j];
137 [ # # ]: 0 : if(e > error) error = e;
138 : : }
139 : : }
140 : 0 : step++;
141 : :
142 : 0 : std::cout << "EquipotentialSmooth smoothing step = " << step << "\tError = " << error << std::endl;
143 [ # # ]: 0 : if (error < epilson) break;
144 : 0 : }
145 [ + - ][ + - ]: 1872 : }
146 : :
|
