From 4457ab38e8ad55a205a547fe441d13def2a41757 Mon Sep 17 00:00:00 2001 From: Kirill Terekhov Date: Fri, 14 Apr 2017 17:14:46 -0700 Subject: [PATCH] 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. --- Examples/MFD-ES/main.cpp | 90 ++++++++++++++++------ Source/Headers/inmost_mesh.h | 11 +++ Source/Mesh/geometry.cpp | 142 ++++++++++++++++++++--------------- Source/Mesh/parallel.cpp | 90 +++++++++++++++++++++- 4 files changed, 249 insertions(+), 84 deletions(-) diff --git a/Examples/MFD-ES/main.cpp b/Examples/MFD-ES/main.cpp index 9ac265b..aed2c9c 100644 --- a/Examples/MFD-ES/main.cpp +++ b/Examples/MFD-ES/main.cpp @@ -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 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 { diff --git a/Source/Headers/inmost_mesh.h b/Source/Headers/inmost_mesh.h index de86a85..67fdd27 100644 --- a/Source/Headers/inmost_mesh.h +++ b/Source/Headers/inmost_mesh.h @@ -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; diff --git a/Source/Mesh/geometry.cpp b/Source/Mesh/geometry.cpp index 9bb16b5..fcd1865 100644 --- a/Source/Mesh/geometry.cpp +++ b/Source/Mesh/geometry.cpp @@ -3,6 +3,11 @@ #include 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 nodes = Element(this,e)->getNodes(); - - Storage::real_array x0 = nodes[0].Coords(), a = x0, b; - for(ElementArray::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 nodes = Element(this,e)->getNodes(); - if( nodes.size() > 1 ) + memset(ret,0,sizeof(real)*mdim); + if( edim == 2 )//&& mdim == 3) + { + ElementArray nodes = Element(this,e)->getNodes(); + + Storage::real_array x0 = nodes[0].Coords(), a = x0, b; + for(ElementArray::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 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); diff --git a/Source/Mesh/parallel.cpp b/Source/Mesh/parallel.cpp index f8ba644..117e528 100644 --- a/Source/Mesh/parallel.cpp +++ b/Source/Mesh/parallel.cpp @@ -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 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::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(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++) -- 2.26.2