Commit 4457ab38 authored by Kirill Terekhov's avatar Kirill Terekhov

MFD-ES works in parallel

Restricted communication of geometrical parameters when exchanging
parts of the mesh between processors.

Added function MarkNormalOrientation to detect situation when ghost
faces are oriented opposite to shared face on owner processor.

Corrected MFD-ES to account for orientation at ghost faces and to flip
sign of normal velocity component on ghost faces when it is known to
differ on master face on owner processor.
parent a71f2eaf
......@@ -73,7 +73,7 @@ real refP(real xyz[3], real nu, real t)
}
// (done)todo: adapt for parallel use
void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real visc, real T)
void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real visc, real T, MarkerType mrk)
{
{
TagReal tag_refP = m.CreateTag("refP" ,DATA_REAL,CELL|FACE,NONE,1);
......@@ -83,6 +83,7 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
for(integer it = 0; it < m.LastLocalID(etype); ++it) if( m.isValidElement(etype,it) )
{
Element e = m.ElementByLocalID(etype,it);
if( e.GetStatus() == Element::Ghost ) continue;
real x[3],p;
e.Centroid(x);
p = refP(x,visc,T);
......@@ -100,6 +101,7 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
for(integer it = 0; it < m.LastLocalID(etype); ++it) if( m.isValidElement(etype,it) )
{
Element e = m.ElementByLocalID(etype,it);
if( e.GetStatus() == Element::Ghost ) continue;
real x[3];
e.Centroid(x);
tag_refP[e] = refP(x,visc,T) - p_ref_min;
......@@ -112,6 +114,7 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
TagReal tag_E = m.CreateTag("ERROR_P" ,DATA_REAL,CELL|FACE,NONE,1);
TagRealArray tag_EU = m.CreateTag("ERROR_U" ,DATA_REAL,CELL,NONE,3);
TagRealArray tag_U = m.CreateTag("U" ,DATA_REAL,CELL,NONE,3);
TagRealArray tag_refU = m.CreateTag("refU" ,DATA_REAL,CELL,NONE,3);
TagReal tag_EnU = m.CreateTag("ERROR_nU",DATA_REAL,FACE,NONE,1);
real C, L2, volume, Cu[3], L2u[3], Cnu, L2nu;
C = L2 = volume = 0.0;
......@@ -120,7 +123,7 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
{
Cell c = m.CellByLocalID(q);
if( c.GetStatus() == Element::Ghost ) continue;
real x[3], z[3];
real x[3], z[3], n[3];
c.Centroid(x);
real err = tag_P[c] - p_min - tag_refP[c];
real vol = c.Volume();
......@@ -129,15 +132,25 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
volume += vol;
tag_E[c] = err;
real V[3] = {0,0,0}; //restored velocity
real V[3],Q[3]; //restored velocity
V[0] = V[1] = V[2] = 0;
Q[0] = Q[1] = Q[2] = 0;
real a,s,v;
ElementArray<Face> cfaces = c.getFaces();
for(integer kt = 0; kt < cfaces.size(); ++kt)
{
//cfaces[kt].FixNormalOrientation();
//if( !cfaces[kt].FixNormalOrientation() ) std::cout << "Orientation is bad" << std::endl;
cfaces[kt].Centroid(z);
cfaces[kt].UnitNormal(n);
a = cfaces[kt].Area();
s = cfaces[kt].FaceOrientedOutside(c) ? 1.0 : -1.0;
s*= cfaces[kt].GetMarker(mrk) ? -1.0 : 1.0;
v = tag_nU[cfaces[kt]]*s*a;
//v = n[0]*refU(z,visc,T).GetValue()
// + n[1]*refV(z,visc,T).GetValue()
// + n[2]*refW(z,visc,T).GetValue();
//v*= a*s;
V[0] += v*(z[0]-x[0]);
V[1] += v*(z[1]-x[1]);
V[2] += v*(z[2]-x[2]);
......@@ -150,23 +163,27 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
tag_U[c][1] = V[1];
tag_U[c][2] = V[2];
err = V[0] - refU(x,visc,T).GetValue();
tag_refU[c][0] = refU(x,visc,T).GetValue();
tag_refU[c][1] = refV(x,visc,T).GetValue();
tag_refU[c][2] = refW(x,visc,T).GetValue();
err = V[0] - tag_refU[c][0];
if( Cu[0] < fabs(err) ) Cu[0] = fabs(err);
L2u[0] += err*err*vol;
tag_EU[c][0] = err;
err = V[1] - refV(x,visc,T).GetValue();
err = V[1] - tag_refU[c][1];
if( Cu[1] < fabs(err) ) Cu[1] = fabs(err);
L2u[1] += err*err*vol;
tag_EU[c][1] = err;
err = V[2] - refW(x,visc,T).GetValue();
err = V[2] - tag_refU[c][2];
if( Cu[2] < fabs(err) ) Cu[2] = fabs(err);
L2u[2] += err*err*vol;
tag_EU[c][2] = err;
}
m.ExchangeData(tag_U,CELL,0);
m.ExchangeData(tag_EU,CELL,0);
m.ExchangeData(tag_refU,CELL,0);
volume = m.Integrate(volume);
L2 = sqrt(m.Integrate(L2)/volume);
m.Integrate(L2u,3);
......@@ -182,6 +199,7 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
std::cout << "Reconstructed V error on cells, C-norm " << Cu[1] << " L2-norm " << L2u[1] << std::endl;
std::cout << "Reconstructed W error on cells, C-norm " << Cu[2] << " L2-norm " << L2u[2] << std::endl;
}
C = L2 = volume = 0.0;
Cnu = L2nu = 0.0;
for( int q = 0; q < m.FaceLastLocalID(); ++q ) if( m.isValidFace(q) )
......@@ -211,7 +229,7 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
V[0] = refU(x,visc,T).GetValue();
V[1] = refV(x,visc,T).GetValue();
V[2] = refW(x,visc,T).GetValue();
nV = n[0]*V[0]+n[1]*V[1]+n[2]*V[2];
nV = (n[0]*V[0]+n[1]*V[1]+n[2]*V[2])*(f.GetMarker(mrk)?-1:1);
err = tag_nU[f] - nV;
if( Cnu < fabs(err) ) Cnu = fabs(err);
......@@ -241,6 +259,7 @@ void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real
{
a = cfaces[kt].Area();
s = cfaces[kt].FaceOrientedOutside(c)?1.0:-1.0;
s*= cfaces[kt].GetMarker(mrk) ? -1.0 : 1.0;
div += tag_nU[cfaces[kt]]*a*s;
}
if( div > div_max ) div_max = div;
......@@ -332,7 +351,9 @@ int main(int argc,char ** argv)
for(int k = 0; k < 3; ++k)
c[k] = 2.0*(c[k]-xmin[k])/(xmax[k]-xmin[k])-1.0;
}
{ // prepare geometrical data on the mesh
Mesh::GeomParam table;
table[CENTROID] = CELL | FACE; //Compute averaged center of mass
......@@ -340,9 +361,11 @@ int main(int argc,char ** argv)
table[ORIENTATION] = FACE; //Check and fix normal orientation
table[MEASURE] = CELL | FACE; //Compute volumes and areas
//table[BARYCENTER] = CELL | FACE; //Compute volumetric center of mass
m.RemoveGeometricData(table);
m.PrepareGeometricData(table); //Ask to precompute the data
}
// data tags for
TagReal tag_nU; // Normal component of velocity
TagReal tag_nUold; // on previous step
......@@ -352,11 +375,19 @@ int main(int argc,char ** argv)
TagRealArray tag_W; // Gradient matrix
TagInteger tag_i; // row number for current cell in back cell matrix
if( m.GetProcessorsNumber() > 1 ) //skip for one processor job
{ // Exchange ghost cells
m.ExchangeGhost(1,FACE); // Produce layer of ghost cells
}
MarkerType nrm_ornt = m.CreateMarker();
MarkerType boundary = m.CreateMarker();
m.MarkNormalOrientation(nrm_ornt);
m.MarkBoundaryFaces(boundary);
{ //initialize data
tag_P = m.CreateTag("PRESSURE",DATA_REAL,CELL|FACE,NONE,1); // Create a new tag for the pressure
tag_W = m.CreateTag("W",DATA_REAL,CELL,NONE);
......@@ -416,11 +447,12 @@ int main(int argc,char ** argv)
u[0] = refU(x,visc,0).GetValue();
u[1] = refV(x,visc,0).GetValue();
u[2] = refW(x,visc,0).GetValue();
tag_nUold[f] = tag_nU[f] = u[0]*n[0]+u[1]*n[1]+u[2]*n[2];
tag_nUold[f] = tag_nU[f] = (u[0]*n[0]+u[1]*n[1]+u[2]*n[2])*(f.GetMarker(nrm_ornt)?-1:1);
//pressure
tag_P[f] = refP(x,visc,0);
}
//initialize pressure
#if defined(USE_OMP)
#pragma omp parallel for
......@@ -432,12 +464,17 @@ int main(int argc,char ** argv)
c.Centroid(x);
tag_P[c] = refP(x,visc,0);
}
//CheckResidual(m,tag_P,tag_nU,visc,0,nrm_ornt);
m.ExchangeData(tag_nU,FACE,0);
m.ExchangeData(tag_nUold,FACE,0);
//CheckResidual(m,tag_P,tag_nU,visc,0,nrm_ornt);
} //end of initialize data
MarkerType boundary = m.CreateMarker();
m.MarkBoundaryFaces(boundary);
if( m.GetProcessorRank() == 0 ) std::cout << "Initialization done" << std::endl;
{ //Main loop for problem solution
......@@ -510,7 +547,7 @@ int main(int argc,char ** argv)
if( report_residual )
{
if( m.GetProcessorRank() == 0 ) std::cout << "Before velocity update " << std::endl;
CheckResidual(m,tag_P,tag_nU,visc,t);
CheckResidual(m,tag_P,tag_nU,visc,t,nrm_ornt);
}
......@@ -520,7 +557,9 @@ int main(int argc,char ** argv)
for(integer it = 0; it < m.FaceLastLocalID(); ++it) if( m.isValidFace(it) )
{
Face f = m.FaceByLocalID(it);
real x[3], n[3], V[3], nV, U, Uold;
if( f.GetStatus() == Element::Ghost ) continue;
real x[3], n[3], V[3], nV, U, Uold, s;
s = f.GetMarker(nrm_ornt) ? -1 : 1;
f.UnitNormal(n);
f.Centroid(x);
U = tag_nU[f];
......@@ -531,7 +570,7 @@ int main(int argc,char ** argv)
V[0] = refU(x,visc,t+dT).GetValue();
V[1] = refV(x,visc,t+dT).GetValue();
V[2] = refW(x,visc,t+dT).GetValue();
nV = n[0]*V[0]+n[1]*V[1]+n[2]*V[2];
nV = (n[0]*V[0]+n[1]*V[1]+n[2]*V[2])*s;
tag_nU[f] = nV;
tag_nUold[f] = U;
}
......@@ -564,19 +603,19 @@ int main(int argc,char ** argv)
V[0] = u.GetValue()*u.GetRow()[0] + v.GetValue()*u.GetRow()[1] + w.GetValue()*u.GetRow()[2];
V[1] = u.GetValue()*v.GetRow()[0] + v.GetValue()*v.GetRow()[1] + w.GetValue()*v.GetRow()[2];
V[2] = u.GetValue()*w.GetRow()[0] + v.GetValue()*w.GetRow()[1] + w.GetValue()*w.GetRow()[2];
nV = n[0]*V[0]+n[1]*V[1]+n[2]*V[2];
nV = (n[0]*V[0]+n[1]*V[1]+n[2]*V[2])*s;
tag_nU[f] = -(beta*U + gamma*Uold + (nV + ngradP)*dT)/alpha;
tag_nUold[f] = U;
}
}
m.ExchangeData(tag_nU,FACE,0);
m.ExchangeData(tag_nUold,FACE,0);
if( report_residual )
{
if( m.GetProcessorRank() == 0 ) std::cout << "After velocity update " << std::endl;
CheckResidual(m,tag_P,tag_nU,visc,t+dT);
CheckResidual(m,tag_P,tag_nU,visc,t+dT,nrm_ornt);
}
R.Clear(); //clean up the residual
......@@ -630,6 +669,7 @@ int main(int argc,char ** argv)
{
a = cfaces[kt].Area();
s = cfaces[kt].FaceOrientedOutside(c)?1.0:-1.0;
s*= cfaces[kt].GetMarker(nrm_ornt) ? -1 : 1;
div += tag_nU[cfaces[kt]]*a*s;
}
tag_F[c] = div/dT;
......@@ -680,6 +720,7 @@ int main(int argc,char ** argv)
S.SetParameter("absolute_tolerance", "1.0e-12");
S.SetParameter("drop_tolerance", "5.0e-2");
S.SetParameter("reuse_tolerance", "1.0e-3");
S.SetParameter("schwartz_overlap", "2");
S.SetMatrix(R.GetJacobian());;
if( S.Solve(R.GetResidual(),Update) )
{
......@@ -694,6 +735,7 @@ int main(int argc,char ** argv)
for(integer it = 0; it < m.LastLocalID(etype); ++it) if( m.isValidElement(etype,it) )
{
Element e = m.ElementByLocalID(etype,it);
if( e.GetStatus() == Element::Ghost ) continue;
tag_Q[e] -= Update[Q.Index(e)];
//if( tag_Q[e] < Q_min ) Q_min = tag_Q[e];
//if( tag_Q[e] > Q_max ) Q_max = tag_Q[e];
......@@ -756,7 +798,7 @@ int main(int argc,char ** argv)
if( report_residual )
{
if( m.GetProcessorRank() == 0 ) std::cout << "After projection" << std::endl;
CheckResidual(m,tag_P,tag_nU,visc,t+dT);
CheckResidual(m,tag_P,tag_nU,visc,t+dT,nrm_ornt);
}
if( write_steps )
......@@ -787,7 +829,7 @@ int main(int argc,char ** argv)
}
CheckResidual(m,tag_P,tag_nU,visc,T);
CheckResidual(m,tag_P,tag_nU,visc,T,nrm_ornt);
if( m.GetProcessorsNumber() == 1 )
......@@ -797,7 +839,7 @@ int main(int argc,char ** argv)
m.ReleaseMarker(boundary,FACE);
m.ReleaseMarker(nrm_ornt,FACE);
}
else
{
......
......@@ -3080,6 +3080,7 @@ namespace INMOST
private:
void RestoreGeometricTags();
void RepairGeometricTags();
public:
bool HideGeometricData (GeometricData type, ElementType mask) {return remember[type][ElementNum(mask)-1] = false;}
bool ShowGeometricData (GeometricData type, ElementType mask) {return remember[type][ElementNum(mask)-1] = true;}
public:
......@@ -3104,6 +3105,16 @@ namespace INMOST
/// Sets marker for all the faces that have only one neighbouring cell, works correctly in parallel environment.
/// @param boundary_marker Non-private marker that will indicate boundary faces.
void MarkBoundaryFaces(MarkerType boundary_marker);
/// This function should be used to detect normal inversion on ghost interfaces
/// with respect to normal orientation on owner of the interface.
/// Due to automatic control over normal orientation in the grid, it may
/// happen that some ghost faces have different orientation rather then
/// face on owner processor.
/// It may happen that some data depends on normal orientation, then
/// one should be aware on a local processor orientation may be different
/// from owner value, then the data may have incorrect sign.
/// @param mrk Non-private marker that will indicate inverted normals
void MarkNormalOrientation(MarkerType mrk);
//implemented in modify.cpp
private:
MarkerType hide_element, new_element, temp_hide_element;
......
......@@ -3,6 +3,11 @@
#include <deque>
using namespace std;
const std::string normal_name = "PROTECTED_GEOM_UTIL_NORMAL";
const std::string measure_name = "PROTECTED_GEOM_UTIL_MEASURE";
const std::string centroid_name = "PROTECTED_GEOM_UTIL_CENTROID";
const std::string barycenter_name = "PROTECTED_GEOM_UTIL_BARYCENTER";
namespace INMOST
{
typedef struct orient_face_t
......@@ -548,10 +553,10 @@ namespace INMOST
{
switch(gtype)
{
case MEASURE: measure_tag = GetTag("GEOM_UTIL_MEASURE"); break;
case CENTROID: centroid_tag = GetTag("GEOM_UTIL_CENTROID"); break;
case BARYCENTER: barycenter_tag = GetTag("GEOM_UTIL_BARYCENTER"); break;
case NORMAL: normal_tag = GetTag("GEOM_UTIL_NORMAL"); break;
case MEASURE: measure_tag = GetTag(measure_name); break;
case CENTROID: centroid_tag = GetTag(centroid_name); break;
case BARYCENTER: barycenter_tag = GetTag(barycenter_name); break;
case NORMAL: normal_tag = GetTag(normal_name); break;
}
}
}
......@@ -559,30 +564,30 @@ namespace INMOST
void Mesh::RepairGeometricTags()
{
if( HaveTag("GEOM_UTIL_MEASURE") )
if( HaveTag(measure_name) )
{
measure_tag = GetTag("GEOM_UTIL_MEASURE");
measure_tag = GetTag(measure_name);
for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype))
if( measure_tag.isDefined(etype) && !HaveGeometricData(MEASURE,etype) )
ShowGeometricData(MEASURE,etype);
}
if( HaveTag("GEOM_UTIL_CENTROID") )
if( HaveTag(centroid_name) )
{
centroid_tag = GetTag("GEOM_UTIL_CENTROID");
centroid_tag = GetTag(centroid_name);
for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype))
if( centroid_tag.isDefined(etype) && !HaveGeometricData(CENTROID,etype) )
ShowGeometricData(CENTROID,etype);
}
if( HaveTag("GEOM_UTIL_BARYCENTER") )
if( HaveTag(barycenter_name) )
{
barycenter_tag = GetTag("GEOM_UTIL_BARYCENTER");
barycenter_tag = GetTag(barycenter_name);
for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype))
if( barycenter_tag.isDefined(etype) && !HaveGeometricData(BARYCENTER,etype) )
ShowGeometricData(BARYCENTER,etype);
}
if( HaveTag("GEOM_UTIL_NORMAL") )
if( HaveTag(normal_name) )
{
normal_tag = GetTag("GEOM_UTIL_NORMAL");
normal_tag = GetTag(normal_name);
for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype))
if( normal_tag.isDefined(etype) && !HaveGeometricData(NORMAL,etype) )
ShowGeometricData(NORMAL,etype);
......@@ -619,7 +624,7 @@ namespace INMOST
{
if( (mask & etype) && !HaveGeometricData(MEASURE,etype))
{
measure_tag = CreateTag("GEOM_UTIL_MEASURE",DATA_REAL,etype,NONE,1);
measure_tag = CreateTag(measure_name,DATA_REAL,etype,NONE,1);
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -639,7 +644,7 @@ namespace INMOST
{
if( (mask & etype) && !HaveGeometricData(CENTROID,etype))
{
centroid_tag = CreateTag("GEOM_UTIL_CENTROID",DATA_REAL,etype,NONE,GetDimensions());
centroid_tag = CreateTag(centroid_name,DATA_REAL,etype,NONE,GetDimensions());
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -659,7 +664,7 @@ namespace INMOST
{
if( (mask & etype) && !HaveGeometricData(BARYCENTER,etype))
{
barycenter_tag = CreateTag("GEOM_UTIL_BARYCENTER",DATA_REAL,etype,NONE,GetDimensions());
barycenter_tag = CreateTag(barycenter_name,DATA_REAL,etype,NONE,GetDimensions());
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -679,7 +684,7 @@ namespace INMOST
{
if( (mask & etype) && !HaveGeometricData(NORMAL,etype))
{
normal_tag = CreateTag("GEOM_UTIL_NORMAL",DATA_REAL,etype,NONE,GetDimensions());
normal_tag = CreateTag(normal_name,DATA_REAL,etype,NONE,GetDimensions());
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -1027,57 +1032,74 @@ namespace INMOST
}
break;
case NORMAL:
if( HaveGeometricData(NORMAL,etype) )
memcpy(ret,MGetDenseLink(e,normal_tag),sizeof(real)*mdim);
else
{
memset(ret,0,sizeof(real)*mdim);
if( edim == 2 )//&& mdim == 3)
// real sret[3];
// bool cmp = false;
if( HaveGeometricData(NORMAL,etype) )
{
ElementArray<Node> nodes = Element(this,e)->getNodes();
Storage::real_array x0 = nodes[0].Coords(), a = x0, b;
for(ElementArray<Node>::size_type i = 0; i < nodes.size(); i++)
{
b = nodes[(i+1)%nodes.size()].Coords();
ret[0] += (a[1]-x0[1])*(b[2]-x0[2]) - (a[2]-x0[2])*(b[1]-x0[1]);
ret[1] += (a[2]-x0[2])*(b[0]-x0[0]) - (a[0]-x0[0])*(b[2]-x0[2]);
ret[2] += (a[0]-x0[0])*(b[1]-x0[1]) - (a[1]-x0[1])*(b[0]-x0[0]);
a.swap(b);
}
/*
for(unsigned i = 0; i < nodes.size(); i++)
{
Storage::real_array a = nodes[i].Coords();
Storage::real_array b = nodes[(i+1)%nodes.size()].Coords();
ret[0] += a[1]*b[2] - a[2]*b[1];
ret[1] += a[2]*b[0] - a[0]*b[2];
ret[2] += a[0]*b[1] - a[1]*b[0];
}
*/
ret[0] *= 0.5;
ret[1] *= 0.5;
ret[2] *= 0.5;
memcpy(ret,MGetDenseLink(e,normal_tag),sizeof(real)*mdim);
// cmp = true;
}
else if( edim == 1 )//&& mdim == 2 )
else
{
ElementArray<Node> nodes = Element(this,e)->getNodes();
if( nodes.size() > 1 )
memset(ret,0,sizeof(real)*mdim);
if( edim == 2 )//&& mdim == 3)
{
ElementArray<Node> nodes = Element(this,e)->getNodes();
Storage::real_array x0 = nodes[0].Coords(), a = x0, b;
for(ElementArray<Node>::size_type i = 0; i < nodes.size(); i++)
{
b = nodes[(i+1)%nodes.size()].Coords();
ret[0] += (a[1]-x0[1])*(b[2]-x0[2]) - (a[2]-x0[2])*(b[1]-x0[1]);
ret[1] += (a[2]-x0[2])*(b[0]-x0[0]) - (a[0]-x0[0])*(b[2]-x0[2]);
ret[2] += (a[0]-x0[0])*(b[1]-x0[1]) - (a[1]-x0[1])*(b[0]-x0[0]);
a.swap(b);
}
/*
for(unsigned i = 0; i < nodes.size(); i++)
{
Storage::real_array a = nodes[i].Coords();
Storage::real_array b = nodes[(i+1)%nodes.size()].Coords();
ret[0] += a[1]*b[2] - a[2]*b[1];
ret[1] += a[2]*b[0] - a[0]*b[2];
ret[2] += a[0]*b[1] - a[1]*b[0];
}
*/
ret[0] *= 0.5;
ret[1] *= 0.5;
ret[2] *= 0.5;
}
else if( edim == 1 )//&& mdim == 2 )
{
Storage::real_array a = nodes[0].Coords();
Storage::real_array b = nodes[1].Coords();
ret[0] = b[1] - a[1];
ret[1] = a[0] - b[0];
Storage::real l = ::sqrt(ret[0]*ret[0]+ret[1]*ret[1]);
if( l )
ElementArray<Node> nodes = Element(this,e)->getNodes();
if( nodes.size() > 1 )
{
ret[0] /= l;
ret[1] /= l;
Storage::real_array a = nodes[0].Coords();
Storage::real_array b = nodes[1].Coords();
ret[0] = b[1] - a[1];
ret[1] = a[0] - b[0];
Storage::real l = ::sqrt(ret[0]*ret[0]+ret[1]*ret[1]);
if( l )
{
ret[0] /= l;
ret[1] /= l;
}
l = ::sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]));
ret[0] *= l;
ret[1] *= l;
}
l = ::sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]));
ret[0] *= l;
ret[1] *= l;
}
// real err = 0;
// for(int k = 0; k < 3; ++k)
// err += (ret[0]-sret[0])*(ret[0]-sret[0]);
// err = sqrt(err);
// if( err > 1.0e-8 && cmp )
// {
// std::cout << "face " << GetHandleID(e) << " rank " << GetProcessorRank() << std::endl;
// std::cout << "cn " << ret[0] << " " << ret[1] << " " << ret[2] << std::endl;
// std::cout << "sn " << sret[0] << " " << sret[1] << " " << sret[2] << std::endl;
// }
}
}
break;
......@@ -1261,6 +1283,7 @@ namespace INMOST
c1->Centroid(ccnt);
for(unsigned j = 0; j < data.size(); j++)
{
//std::cout << (data[j].GetPrivateMarker(rev) ? 0:1);
//compute normal to face
data[j].Centroid(cnt);
data[j].Normal(nrm);
......@@ -1268,6 +1291,7 @@ namespace INMOST
cnt[r] = cnt[r]-ccnt[r];
measure += (data[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*vec_dot_product(cnt,nrm,3);
}
//std::cout << "cur" << (cur->GetPrivateMarker(rev) ? 0:1) << " " << measure << " ";
//bool have_rev = cur->GetPrivateMarker(rev);
data.RemPrivateMarker(rev);
mesh->ReleasePrivateMarker(rev);
......
......@@ -3216,6 +3216,9 @@ namespace INMOST
int position = 0;
elements_by_type selems;
MarkerType unpack_tags_mrk = CreateMarker();
//MarkerType swap_normal = 0;
//if( HaveGeometricData(NORMAL,FACE) )
// swap_normal = CreateMarker();
//TODO 46 old
//elements_by_type unpack_tags;
std::vector<HandleType> old_nodes(NumberOfNodes());
......@@ -3478,15 +3481,40 @@ namespace INMOST
}
else
{
//current solution - no communication of geometric data
//here we should check tat remote face orientation match
//existing face orientation, otherwise we would need
//to swap normal if it is precomputed and comes from remote
//processor the worong way
if( IntegerDF(new_face,tag_owner) != mpirank )
{
//TODO 46 old
//unpack_tags[2].push_back(new_face);
SetMarker(new_face,unpack_tags_mrk);
++marked_for_data;
// if( swap_normal )//&& GetMarker(new_face,unpack_tags_mrk) )
// {
// const Element::adj_type & lc = LowConn(new_face);
// assert(lc.size() > 1 );
// assert(lc.size() == f_edges.size());
// if(lc.size() != f_edges.size() )
// std::cout << "This should not happen" << std::endl;
// for(ElementArray<Edge>::size_type q = 0; q < f_edges.size(); ++q)
// {
// if( f_edges.at(q) == lc[0] )
// {
// if( f_edges.at((q+1)%f_edges.size()) != lc[1] )
// SetMarker(new_face,swap_normal);
// break;
// }
// }
// }
}
++found;
}
selems[2].push_back(new_face);
shift+= low_conn_size[i];
}
......@@ -3636,6 +3664,24 @@ namespace INMOST
for(integer k = ElementNum(NODE); k <= ElementNum(CELL); ++k)
if( !selems[k].empty() ) RemMarkerArray(&selems[k][0],static_cast<enumerator>(selems[k].size()),unpack_tags_mrk);
ReleaseMarker(unpack_tags_mrk);
// if( swap_normal )
// {
// for(enumerator q = 0; q < selems[2].size(); ++q)
// {
// if( GetMarker(selems[2][q],swap_normal) )
// {
// real_array nrm = RealArrayDF(selems[2][q],GetGeometricTag(NORMAL));
// std::cout << "Swap normal face " << GetHandleID(selems[2][q]) << " rank " << GetProcessorRank() << std::endl;
// for(real_array::size_type it = 0; it < nrm.size(); ++it)
// nrm[it] = -nrm[it];
// RemMarker(selems[2][q],swap_normal);
// }
// else
// std::cout << "No swap normal face " << GetHandleID(selems[2][q]) << " rank " << GetProcessorRank() << std::endl;
//
// }
// ReleaseMarker(swap_normal);
// }
time = Timer() - time;
REPORT_STR("unpack tag data");
REPORT_VAL("time", time);
......@@ -4767,7 +4813,49 @@ namespace INMOST
EXIT_FUNC();
}
void Mesh::BeginSequentialCode()
void UnpackMarkNormalOrientation(const Tag & tag,
const Element & element,
const INMOST_DATA_BULK_TYPE * data,
INMOST_DATA_ENUM_TYPE size)
{
Storage::real_array local_nrm = element.RealArrayDF(tag);
const Storage::real * remote_nrm = (const Storage::real *) data;
assert(size == local_nrm.size());
Storage::real dot = 0;
for(Storage::enumerator k = 0; k < local_nrm.size(); ++k)
dot += local_nrm[k]*remote_nrm[k];
local_nrm[0] = dot;
}
void Mesh::MarkNormalOrientation(MarkerType mrk)
{
if( m_state == Serial ) return;
ENTER_FUNC();
#if defined(USE_MPI)
#if !defined(USE_PARALLEL_STORAGE)
parallel_storage ghost_elements, shared_elements;
GatherParallelStorage(ghost_elements,shared_elements,mask);
#endif //USE_PARALLEL_STORAGE
TagRealArray tag_nrm = CreateTag("TEMPORARY_NORMAL",DATA_REAL,FACE,NONE,3);
for(iteratorFace it = BeginFace(); it != EndFace(); ++it)
if( it->GetStatus() & (Element::Ghost | Element::Shared) )
it->UnitNormal(tag_nrm[it->self()].data());
exchange_data storage;
ExchangeDataInnerBegin(tag_set(1,tag_nrm),shared_elements,ghost_elements,FACE,0,storage);
ExchangeDataInnerEnd(tag_set(1,tag_nrm),shared_elements,ghost_elements,FACE,0,UnpackMarkNormalOrientation,storage);
for(iteratorFace it = BeginFace(); it != EndFace(); ++it)
if( it->GetStatus() == Element::Ghost && tag_nrm[it->self()][0] < 0 )
it->SetMarker(mrk);
DeleteTag(tag_nrm);
#else //USE_MPI
(void) mrk;
#endif //USE_MPI
EXIT_FUNC();
}
void Mesh::BeginSequentialCode()
{
#if defined(USE_MPI)
for(int i = 0; i < GetProcessorRank(); i++)
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Pr