MeshKit
1.0
|
00001 #include "EquipotentialSmooth.hpp" 00002 #include <algorithm> 00003 00004 00005 //---------------------------------------------------------------------------// 00006 // construction function for EquipotentialSmooth class 00007 EquipotentialSmooth::EquipotentialSmooth() 00008 { 00009 00010 00011 } 00012 00013 00014 00015 //---------------------------------------------------------------------------// 00016 // deconstruction function for EquipotentialSmooth class 00017 EquipotentialSmooth::~EquipotentialSmooth() 00018 { 00019 00020 } 00021 00022 //setup the data for EquipotentialSmooth 00023 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) 00024 { 00025 //pass the node coordinates to the local variables 00026 coordinates.resize(coords.size()); 00027 for (unsigned int i = 0; i < coords.size(); i++){ 00028 coordinates[i].resize(coords[i].size()); 00029 for (unsigned int j = 0; j < coords[i].size(); j++) 00030 coordinates[i][j] = coords[i][j]; 00031 } 00032 00033 //pass the connectivity information to the local variables 00034 adjElements.resize(AdjElements.size()); 00035 for (unsigned int i = 0; i < AdjElements.size(); i++){ 00036 for (std::set<int>::iterator it = AdjElements[i].begin(); it != AdjElements[i].end(); it++) 00037 adjElements[i].insert(*it); 00038 } 00039 00040 quads.resize(Quads.size()); 00041 for (unsigned int i = 0; i < Quads.size(); i++){ 00042 quads[i].resize(Quads[i].size()); 00043 for (unsigned int j = 0; j < Quads[i].size(); j++) 00044 quads[i][j] = Quads[i][j]; 00045 } 00046 00047 //pass the boundary info to the local variables 00048 isBoundary.resize(isBnd.size()); 00049 for (unsigned int i = 0; i < isBnd.size(); i++) 00050 isBoundary[i] = isBnd[i]; 00051 00052 //pass the connectivity to the local variables 00053 connect.resize(conn.size()); 00054 for (unsigned int i = 0; i < conn.size(); i++){ 00055 connect[i].resize(conn[i].size()); 00056 for (unsigned int j = 0; j < conn[i].size(); j++) 00057 connect[i][j] = conn[i][j]; 00058 } 00059 //done 00060 } 00061 00062 //return the results 00063 void EquipotentialSmooth::GetCoords(std::vector<std::vector<double> > &coords) 00064 { 00065 //return the results to the user 00066 coords.resize(coordinates.size()); 00067 for (unsigned int i = 0; i < coords.size(); i++){ 00068 if (!isBoundary[i]){ 00069 coords[i].resize(coordinates[i].size()); 00070 for (unsigned int j = 0; j < coordinates[i].size(); j++) 00071 coords[i][j] = coordinates[i][j]; 00072 } 00073 } 00074 } 00075 00076 //Execute function 00077 void EquipotentialSmooth::Execute() 00078 { 00079 /* Algorithm --- Iterative smoothing */ 00080 /* p'(i) = p(i) + sum{w(n)*(p(n)-p(i)), n = 1...N} elements 1...N are the neighbors of node i */ 00081 /* w(1) = -belta/2 w(2) = alpha w(3) = belta/2 w(4) = gamma */ 00082 /* w(5) = -belta/2 w(6) = alpha w(7) = belta/2 w(8) = gamma */ 00083 /* where alpha = xp^2 + yp^2 + zp^2 */ 00084 /* belta = xp*xq + yp*yq + zp*zq */ 00085 /* gamma = xq^2 + yq^2 + zq^2 */ 00086 /* xp = (x2-x6)/2 yp = (y2-y6)/2 zp = (z2-z6)/2 */ 00087 /* xq = (x8-x4)/2 yq = (y8-y4)/2 zq = (z8-z4)/2 */ 00088 00089 // 3------------2------------1 00090 // | | | 00091 // | | | 00092 // | | | 00093 // 4------------i------------8 00094 // | | | 00095 // | | | 00096 // | | | 00097 // 5------------6------------7 00098 00099 double epilson = 0.01; 00100 double error = 0.0; 00101 int step = 0; 00102 while(true){ 00103 error = 0.0; 00104 //loop over all the nodes 00105 for (unsigned int i = 0; i < coordinates.size(); i++){ 00106 if (!(isBoundary[i])){ 00107 double weight[9]; 00108 double alpha = 0.0, belta = 0.0, gamma = 0.0; 00109 double p[3] = {0.0, 0.0, 0.0}; 00110 double q[3] = {0.0, 0.0, 0.0}; 00111 for (int j = 0; j < 3; j++){ 00112 p[j] = (coordinates[connect[i][2]][j] - coordinates[connect[i][6]][j])/2.0; 00113 q[j] = (coordinates[connect[i][4]][j] - coordinates[connect[i][8]][j])/2.0; 00114 } 00115 alpha = pow(p[0],2) + pow(p[1],2) + pow(p[2],2); 00116 gamma = pow(q[0],2) + pow(q[1],2) + pow(q[2],2); 00117 belta = p[0]*q[0] + p[1]*q[1] + p[2]*q[2]; 00118 00119 weight[1] = belta/2.0; 00120 weight[2] = gamma; 00121 weight[3] = -belta/2.0; 00122 weight[4] = alpha; 00123 weight[5] = belta/2.0; 00124 weight[6] = gamma; 00125 weight[7] = -belta/2.0; 00126 weight[8] = alpha; 00127 00128 double sum_dis[3] = {0.0, 0.0, 0.0}; 00129 for (int j = 0; j < 3; j++){ 00130 for (int k = 1; k < 9; k++){ 00131 sum_dis[j]+=weight[k]*coordinates[connect[i][k]][j]/(2.0*(alpha+gamma)); 00132 } 00133 } 00134 double e = fabs(sum_dis[0] - coordinates[i][0]) + fabs(sum_dis[1] - coordinates[i][1]) + fabs(sum_dis[2]- coordinates[i][2]); 00135 for (int j = 0; j < 3; j++) 00136 coordinates[i][j] = sum_dis[j]; 00137 if(e > error) error = e; 00138 } 00139 } 00140 step++; 00141 00142 std::cout << "EquipotentialSmooth smoothing step = " << step << "\tError = " << error << std::endl; 00143 if (error < epilson) break; 00144 } 00145 } 00146