Commit e639c57d authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Fixes to geometrical algorithms for polyhedral mesh; update crevice N-S...

Fixes to geometrical algorithms for polyhedral mesh; update crevice N-S example with parabolic velocity profile
parent c6bf5454
...@@ -119,7 +119,7 @@ int main(int argc, char ** argv) ...@@ -119,7 +119,7 @@ int main(int argc, char ** argv)
bc[*it][1] = 0; bc[*it][1] = 0;
bc[*it][2] = 1; bc[*it][2] = 1;
bc[*it][3] = 0; bc[*it][3] = 0;
bc[*it][4] = 17.3; bc[*it][4] = 17.3 * c[1]*(3.0-c[1])/9.0 * 4.0 * (3.0 / 2.0); //parabolic profile with average of 17.3
bc[*it][5] = 0; bc[*it][5] = 0;
bc[*it][6] = 0; bc[*it][6] = 0;
} }
......
...@@ -27,6 +27,12 @@ namespace INMOST ...@@ -27,6 +27,12 @@ namespace INMOST
vecout[i] = vecin1[i] - vecin2[i]; vecout[i] = vecin1[i] - vecin2[i];
} }
__INLINE static void vec_diff(const Storage::real_array & vecin1, const Storage::real * vecin2, Storage::real * vecout, unsigned int size)
{
for(unsigned int i = 0; i < size; i++)
vecout[i] = vecin1[i] - vecin2[i];
}
__INLINE static void vec_diff(const Storage::real_array & vecin1,const Storage::real_array & vecin2, Storage::real * vecout, unsigned int size) __INLINE static void vec_diff(const Storage::real_array & vecin1,const Storage::real_array & vecin2, Storage::real * vecout, unsigned int size)
{ {
for(unsigned int i = 0; i < size; i++) for(unsigned int i = 0; i < size; i++)
...@@ -1073,34 +1079,30 @@ namespace INMOST ...@@ -1073,34 +1079,30 @@ namespace INMOST
if( nodes.size() > 2 ) if( nodes.size() > 2 )
{ {
*ret = 0; *ret = 0;
real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0},n0[3] = {0,0,0}, ss; real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0},n0[3] = {0,0,0}, cxf[3] = {0,0,0}, ss, at;
//real_array v0 = nodes[0].Coords(); real_array v1, v2;
//real_array v1 = nodes[1].Coords(); Element(this,e)->Centroid(cxf);
//real_array v2 = nodes[2].Coords(); for(int i = 0; i < (int)nodes.size(); i++)
//vec_diff(v1,v0,l1,mdim);
//vec_diff(v2,v0,l2,mdim);
//vec_cross_product(l1,l2,n0);
real_array v0 = nodes[0].Coords(), v1, v2;
for(int i = 1; i < (int)nodes.size()-1; i++)
{ {
v1 = nodes[i].Coords(); v1 = nodes[i].Coords();
v2 = nodes[i+1].Coords(); v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
n0[q] += nt[q]*0.5; n0[q] += nt[q]*0.5;
} }
for(int i = 1; i < (int)nodes.size()-1; i++) for(int i = 0; i < (int)nodes.size(); i++)
{ {
v1 = nodes[i].Coords(); v1 = nodes[i].Coords();
v2 = nodes[i+1].Coords(); v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
ss = vec_dot_product(n0,nt,3); ss = vec_dot_product(n0,nt,3);
if( ss ) ss /= fabs(ss); if( ss ) ss /= fabs(ss);
*ret += sqrt(vec_dot_product(nt,nt,3))*0.5*ss; at = sqrt(vec_dot_product(nt,nt,3))*ss;
*ret += at;
} }
if( *ret != *ret ) std::cout << "area is nan" << std::endl; if( *ret != *ret ) std::cout << "area is nan" << std::endl;
*ret = fabs(*ret); *ret = fabs(*ret);
...@@ -1124,8 +1126,8 @@ namespace INMOST ...@@ -1124,8 +1126,8 @@ namespace INMOST
real vol = 0, a, at; real vol = 0, a, at;
real x[3] = {0,0,0}, n[3] = {0,0,0}, n0[3] = {0,0,0}, s, ss; real x[3] = {0,0,0}, n[3] = {0,0,0}, n0[3] = {0,0,0}, s, ss;
real l1[3] = {0,0,0}, l2[3] = {0,0,0}; real l1[3] = {0,0,0}, l2[3] = {0,0,0};
real nt[3] = {0,0,0}, cx[3] = {0,0,0}; real nt[3] = {0,0,0}, cxf[3] = {0,0,0};
me.Centroid(cx); //me.Centroid(cx);
//if( print ) std::cout << "cx: " << cx[0] << " " << cx[1] << " " << cx[2] << std::endl; //if( print ) std::cout << "cx: " << cx[0] << " " << cx[1] << " " << cx[2] << std::endl;
for(unsigned j = 0; j < faces.size(); j++) for(unsigned j = 0; j < faces.size(); j++)
{ {
...@@ -1138,6 +1140,7 @@ namespace INMOST ...@@ -1138,6 +1140,7 @@ namespace INMOST
x[0] = x[1] = x[2] = 0; x[0] = x[1] = x[2] = 0;
n[0] = n[1] = n[2] = 0; n[0] = n[1] = n[2] = 0;
n0[0] = n0[1] = n0[2] = 0; n0[0] = n0[1] = n0[2] = 0;
cxf[0] = cxf[1] = cxf[2] = 0;
a = 0; a = 0;
//real_array v0 = nodes[0].Coords(); //real_array v0 = nodes[0].Coords();
//real_array v1 = nodes[1].Coords(); //real_array v1 = nodes[1].Coords();
...@@ -1145,27 +1148,29 @@ namespace INMOST ...@@ -1145,27 +1148,29 @@ namespace INMOST
//vec_diff(v1,v0,l1,mdim); //vec_diff(v1,v0,l1,mdim);
//vec_diff(v2,v0,l2,mdim); //vec_diff(v2,v0,l2,mdim);
//vec_cross_product(l1,l2,n0); //vec_cross_product(l1,l2,n0);
real_array v0 = nodes[0].Coords(), v1, v2; real_array /*v0 = nodes[0].Coords(),*/ v1, v2;
for(int i = 1; i < (int)nodes.size()-1; i++) faces[j].Centroid(cxf);
for(int i = 0; i < (int)nodes.size(); i++)
{ {
v1 = nodes[i].Coords(); v1 = nodes[i].Coords();
v2 = nodes[i+1].Coords(); v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
//if( print ) std::cout << "nt: " << nt[0] << " " << nt[1] << " " << nt[2] << std::endl; //if( print ) std::cout << "nt: " << nt[0] << " " << nt[1] << " " << nt[2] << std::endl;
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
n0[q] += nt[q]*0.5; n0[q] += nt[q]*0.5;
} }
//if( print ) std::cout << "n0: " << n0[0] << " " << n0[1] << " " << n0[2] << std::endl; //if( print ) std::cout << "n0: " << n0[0] << " " << n0[1] << " " << n0[2] << std::endl;
for(int k = 1; k < (int)nodes.size()-1; k++) for(int k = 0; k < (int)nodes.size(); k++)
{ {
v1 = nodes[k+0].Coords(); v1 = nodes[k].Coords();
v2 = nodes[k+1].Coords(); v2 = nodes[(k+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
for(int q = 0; q < 3; ++q) nt[q] *= 0.5; for(int q = 0; q < 3; ++q)
nt[q] *= 0.5;
ss = vec_dot_product(n0,nt,3); ss = vec_dot_product(n0,nt,3);
if( ss ) ss /= fabs(ss); if( ss ) ss /= fabs(ss);
at = sqrt(vec_dot_product(nt,nt,3))*ss; at = sqrt(vec_dot_product(nt,nt,3))*ss;
...@@ -1173,14 +1178,22 @@ namespace INMOST ...@@ -1173,14 +1178,22 @@ namespace INMOST
//same as faces[j].Normal(n) //same as faces[j].Normal(n)
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
n[q] += nt[q]; n[q] += nt[q];
//same as faces[j].Centroid(x) //same as faces[j].Barycenter(x)
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
x[q] += at*((v0[q]-cx[q])+(v1[q]-cx[q])+(v2[q]-cx[q]))/3.0; x[q] += at*(/*(v0[q]-cx[q])+*/(v1[q]-cxf[q])+(v2[q]-cxf[q]))/3.0;
a += at; a += at;
} }
//if( print ) std::cout << "n: " << n[0] << " " << n[1] << " " << n[2] << std::endl; //if( print ) std::cout << "n: " << n[0] << " " << n[1] << " " << n[2] << std::endl;
//if( print ) std::cout << "x: " << x[0] << " " << x[1] << " " << x[2] << " a " << a << std::endl; //if( print ) std::cout << "x: " << x[0] << " " << x[1] << " " << x[2] << " a " << a << std::endl;
for(int q = 0; q < 3; ++q) x[q] /= a; if( a )
{
//~ std::cout << "Correction: " << x[0]/a << " " << x[1]/a << " " << x[2]/a << " area " << a << std::endl;
for(int q = 0; q < 3; ++q)
x[q] = x[q]/a + cxf[q];
//~ std::cout << "Centroid: " << cxf[0] << " " << cxf[1] << " " << cxf[2] << std::endl;
//~ std::cout << "Barycenter: " << x[0] << " " << x[1] << " " << x[2] << std::endl;
}
else for(int q = 0; q < 3; ++q) x[q] = cxf[q];
vol += s*vec_dot_product(x,n,3); vol += s*vec_dot_product(x,n,3);
} }
//if( print ) std::cout << "vol: " << vol << std::endl; //if( print ) std::cout << "vol: " << vol << std::endl;
...@@ -1256,7 +1269,7 @@ namespace INMOST ...@@ -1256,7 +1269,7 @@ namespace INMOST
{ {
*ret = 0; *ret = 0;
real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0}; real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0};
real c[3] = {0,0,0}, n0[3] = {0,0,0}, ss, cx[3] = {0,0,0}; real c[3] = {0,0,0}, n0[3] = {0,0,0}, ss, cxf[3] = {0,0,0};
//real_array v0 = nodes[0].Coords(); //real_array v0 = nodes[0].Coords();
//real_array v1 = nodes[1].Coords(); //real_array v1 = nodes[1].Coords();
//real_array v2 = nodes[2].Coords(); //real_array v2 = nodes[2].Coords();
...@@ -1264,33 +1277,37 @@ namespace INMOST ...@@ -1264,33 +1277,37 @@ namespace INMOST
//vec_diff(v2,v0,l2,mdim); //vec_diff(v2,v0,l2,mdim);
//vec_cross_product(l1,l2,n0); //vec_cross_product(l1,l2,n0);
real_array v0 = nodes[0].Coords(), v1, v2; real_array v0 = nodes[0].Coords(), v1, v2;
for(int i = 1; i < (int)nodes.size()-1; i++) Element(this,e).Centroid(cxf);
for(int i = 0; i < (int)nodes.size(); i++)
{ {
v1 = nodes[i].Coords(); v1 = nodes[i].Coords();
v2 = nodes[i+1].Coords(); v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
n0[q] += nt[q]*0.5; n0[q] += nt[q]*0.5;
} }
real a = 0, at; real a = 0, at;
Element(this,e).Centroid(cx); for(int i = 0; i < (int)nodes.size(); i++)
for(int i = 1; i < (int)nodes.size()-1; i++)
{ {
real_array v1 = nodes[i].Coords(); real_array v1 = nodes[i].Coords();
real_array v2 = nodes[i+1].Coords(); real_array v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
ss = vec_dot_product(n0,nt,3); ss = vec_dot_product(n0,nt,3);
if( ss ) ss /= fabs(ss); if( ss ) ss /= fabs(ss);
at = sqrt(vec_dot_product(nt,nt,3))*0.5*ss; at = sqrt(vec_dot_product(nt,nt,3))*0.5*ss;
for(int q = 0; q < mdim; ++q) for(int q = 0; q < mdim; ++q)
c[q] += at*((v0[q]-cx[q])+(v1[q]-cx[q])+(v2[q]-cx[q]))/3.0; c[q] += at*(/*(v0[q]-cxf[q])+*/(v1[q]-cxf[q])+(v2[q]-cxf[q]))/3.0;
a += at; a += at;
} }
for(int q = 0; q < mdim; ++q) ret[q] = c[q]/a+cx[q]; if( a )
{
for(int q = 0; q < mdim; ++q)
ret[q] = c[q]/a+cxf[q];
} else for(int q = 0; q < mdim; ++q) ret[q] = cxf[q];
} }
//std::cout << ret[0] << " " << ret[1] << " " << ret[2] << std::endl; //std::cout << ret[0] << " " << ret[1] << " " << ret[2] << std::endl;
} }
...@@ -1310,7 +1327,7 @@ namespace INMOST ...@@ -1310,7 +1327,7 @@ namespace INMOST
real c[3] = {0,0,0}, n[3] = {0,0,0}; real c[3] = {0,0,0}, n[3] = {0,0,0};
real n0[3] = {0,0,0}, ss; real n0[3] = {0,0,0}, ss;
real l1[3] = {0,0,0}, l2[3] = {0,0,0}; real l1[3] = {0,0,0}, l2[3] = {0,0,0};
real cx[3] = {0,0,0}; real cx[3] = {0,0,0}, cxf[3] = {0,0,0};
me.Centroid(cx); me.Centroid(cx);
for(unsigned j = 0; j < faces.size(); j++) for(unsigned j = 0; j < faces.size(); j++)
{ {
...@@ -1329,24 +1346,26 @@ namespace INMOST ...@@ -1329,24 +1346,26 @@ namespace INMOST
n0[0] = n0[1] = n0[2] = 0; n0[0] = n0[1] = n0[2] = 0;
x[0] = x[1] = x[2] = 0; x[0] = x[1] = x[2] = 0;
n[0] = n[1] = n[2] = 0; n[0] = n[1] = n[2] = 0;
cxf[0] = cxf[1] = cxf[2] = 0;
a = 0; a = 0;
real_array v0 = nodes[0].Coords(), v1, v2; real_array /*v0 = nodes[0].Coords(),*/ v1, v2;
for(int i = 1; i < (int)nodes.size()-1; i++) faces[j].Centroid(cxf);
for(int i = 0; i < (int)nodes.size(); i++)
{ {
v1 = nodes[i].Coords(); v1 = nodes[i].Coords();
v2 = nodes[i+1].Coords(); v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
n0[q] += nt[q]*0.5; n0[q] += nt[q]*0.5;
} }
for(int i = 1; i < (int)nodes.size()-1; i++) for(int i = 0; i < (int)nodes.size(); i++)
{ {
real_array v1 = nodes[i].Coords(); real_array v1 = nodes[i].Coords();
real_array v2 = nodes[i+1].Coords(); real_array v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
nt[q] *= 0.5; nt[q] *= 0.5;
...@@ -1358,13 +1377,17 @@ namespace INMOST ...@@ -1358,13 +1377,17 @@ namespace INMOST
n[q] += nt[q]; n[q] += nt[q];
//same as faces[j].Centroid(x) //same as faces[j].Centroid(x)
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
x[q] += at*((v0[q]-cx[q])+(v1[q]-cx[q])+(v2[q]-cx[q]))/3.0; x[q] += at*(/*(v0[q]-cxf[q])+*/(v1[q]-cxf[q])+(v2[q]-cxf[q]))/3.0;
a += at; a += at;
//second-order midpoint formula //second-order midpoint formula
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
c[q] += s*nt[q]*(pow((v0[q]-cx[q])+(v1[q]-cx[q]),2)+pow((v0[q]-cx[q])+(v2[q]-cx[q]),2)+pow((v1[q]-cx[q])+(v2[q]-cx[q]),2))/24.0; c[q] += s*nt[q]*(pow((/*v0*/cxf[q]-cx[q])+(v1[q]-cx[q]),2)+pow((/*v0*/cxf[q]-cx[q])+(v2[q]-cx[q]),2)+pow((v1[q]-cx[q])+(v2[q]-cx[q]),2))/24.0;
} }
for(int q = 0; q < 3; ++q) x[q] = x[q]/a; if( a )
{
for(int q = 0; q < 3; ++q)
x[q] = x[q]/a + cxf[q];
} else for(int q = 0; q < 3; ++q) x[q] = cxf[q];
vol += s*vec_dot_product(x,n,3); vol += s*vec_dot_product(x,n,3);
} }
if( ornt ) if( ornt )
...@@ -1402,14 +1425,15 @@ namespace INMOST ...@@ -1402,14 +1425,15 @@ namespace INMOST
{ {
ElementArray<Node> nodes = Element(this,e)->getNodes(); ElementArray<Node> nodes = Element(this,e)->getNodes();
real n[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0}; real n[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0};
real nt[3] = {0,0,0}; real nt[3] = {0,0,0}, cxf[3] = {0,0,0};
real_array v0 = nodes[0].Coords(); real_array v0 = nodes[0].Coords();
for(int i = 1; i < (int)nodes.size()-1; i++) Element(this,e)->Centroid(cxf);
for(int i = 0; i < (int)nodes.size(); i++)
{ {
real_array v1 = nodes[i].Coords(); real_array v1 = nodes[i].Coords();
real_array v2 = nodes[i+1].Coords(); real_array v2 = nodes[(i+1)%nodes.size()].Coords();
vec_diff(v1,v0,l1,mdim); vec_diff(v1,cxf,l1,mdim);
vec_diff(v2,v0,l2,mdim); vec_diff(v2,cxf,l2,mdim);
vec_cross_product(l1,l2,nt); vec_cross_product(l1,l2,nt);
for(int q = 0; q < 3; ++q) for(int q = 0; q < 3; ++q)
n[q] += nt[q]*0.5; n[q] += nt[q]*0.5;
...@@ -1626,7 +1650,7 @@ namespace INMOST ...@@ -1626,7 +1650,7 @@ namespace INMOST
{ {
//std::cout << (data[j].GetPrivateMarker(rev) ? 0:1); //std::cout << (data[j].GetPrivateMarker(rev) ? 0:1);
//compute normal to face //compute normal to face
data[j].Centroid(cnt); data[j].Barycenter(cnt);
data[j].Normal(nrm); data[j].Normal(nrm);
measure += (data[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*vec_dot_product(cnt,nrm,3); measure += (data[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*vec_dot_product(cnt,nrm,3);
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment