diff --git a/Examples/MFD-ES/main.cpp b/Examples/MFD-ES/main.cpp index 9ac265bfb441d84aa898f7a5412cc8e9420c5040..aed2c9c9a5d82fbe6c5b9987bc1b94a2ce9ed28c 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 de86a85bd8a817080af53dd9480be32af78f9c80..67fdd27116341f6e0c16c6e37d4fdf76016751f5 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 9bb16b5432a0e5d3992ee128e011bcf35f40670a..fcd1865739fb600465c4d4589d706cd097cef07c 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 f8ba64414a777ab511cdb9eff9746032931b6b1f..117e5289a654b170fc724e94989d047936cfb557 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++)