diff --git a/configure.ac b/configure.ac index bd821ba0..6453d4bb 100644 --- a/configure.ac +++ b/configure.ac @@ -1,4 +1,4 @@ -AC_INIT([netgen],[4.9.7-dev],[],[]) +AC_INIT([netgen],[4.9.8-dev],[],[]) AM_INIT_AUTOMAKE([-Wall -Werror foreign]) AC_CONFIG_MACRO_DIR([m4]) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index 0506aabb..b643c1a5 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -990,14 +990,11 @@ namespace netgen for (int k = 0; k < tas.GetNP(); k++) { tams -> AddPoint (tas.GetPoint(k)); - Vec<3> n = GetSurface(j) -> GetNormalVector (tas.GetPoint(k)); + Vec<3> n = GetSurface(j) -> GetNormalVector (tas.GetPoint(k)); n.Normalize(); if (GetSurface(j)->Inverse()) n *= -1; tams -> AddNormal (n); - //(*testout) << "point " << tas.GetPoint(k) << " normal " << n << endl; - //cout << "added point, normal=" << n << endl; } - BoxSphere<3> surfbox; @@ -1033,9 +1030,12 @@ namespace netgen tri[0] + oldnp, tri[1] + oldnp, tri[2] + oldnp); + + // tams -> AddTriangle (tria); RefineTriangleApprox (locsol, j, box, detail, - tria, *tams, iset); + tria, *tams, iset, 1); + delete locsol; } } @@ -1057,9 +1057,11 @@ namespace netgen double detail, const TATriangle & tria, TriangleApproximation & tams, - IndexSet & iset) + IndexSet & iset, + int level) { - + if (level > 10) return; + //tams.AddTriangle (tria); //(*testout) << "tria " << tams.GetPoint(tria[0]) << " - " << tams.GetPoint(tria[1]) << " - " << tams.GetPoint(tria[2]) // << " ( " << tria[0] << " " << tria[1] << " " << tria[2] << ")" < n; GetSurface(surfind)->Project (newp); + n = GetSurface(surfind)->GetNormalVector (newp); pinds[between[i][2]] = tams.AddPoint (newp); @@ -1377,7 +1380,7 @@ namespace netgen nbox.Set (tams.GetPoint (ntri[0])); nbox.Add (tams.GetPoint (ntri[1])); nbox.Add (tams.GetPoint (ntri[2])); - nbox.Increase (1e-6); + nbox.Increase (1e-8); nbox.CalcDiamCenter(); Solid * nsol = locsol -> GetReducedSolid (nbox); @@ -1385,7 +1388,7 @@ namespace netgen if (nsol) { RefineTriangleApprox (nsol, surfind, nbox, - detail, ntri, tams, iset); + detail, ntri, tams, iset, level+1); delete nsol; } diff --git a/libsrc/csg/csgeom.hpp b/libsrc/csg/csgeom.hpp index 9e22c76b..a96672b0 100644 --- a/libsrc/csg/csgeom.hpp +++ b/libsrc/csg/csgeom.hpp @@ -274,7 +274,8 @@ public: double detail, const TATriangle & tria, TriangleApproximation & tams, - IndexSet & iset); + IndexSet & iset, + int level); const Box<3> & BoundingBox () const { return boundingbox; } diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index cb24a9ba..adf52261 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -529,7 +529,7 @@ namespace netgen { if (osedges.Get(seg.edgenr)) { - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort (); if (osedgesht.Used (i2)) osedgesht.Set (i2, 2); @@ -592,26 +592,26 @@ namespace netgen const Segment & seg = mesh[si]; if (seg.seginfo && seg.edgenr >= 1 && seg.edgenr <= cntedge) { - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort (); if (osedgesht.Used (i2) && osedgesht.Get (i2) == 2 && osedges.Elem(seg.edgenr) == -1) { - Point<3> newp = Center (mesh[PointIndex(seg.p1)], - mesh[PointIndex(seg.p2)]); + Point<3> newp = Center (mesh[PointIndex(seg[0])], + mesh[PointIndex(seg[1])]); ProjectToEdge (geometry.GetSurface(seg.surfnr1), geometry.GetSurface(seg.surfnr2), newp); osedges.Elem(seg.edgenr) = - mesh.AddPoint (newp, mesh[PointIndex(seg.p1)].GetLayer(), EDGEPOINT); + mesh.AddPoint (newp, mesh[PointIndex(seg[0])].GetLayer(), EDGEPOINT); meshpoint_tree -> Insert (newp, osedges.Elem(seg.edgenr)); } } } - + for (int i = 1; i <= nseg; i++) { @@ -621,8 +621,8 @@ namespace netgen if (osedges.Get(seg.edgenr) != -1) { Segment newseg = seg; - newseg.p1 = osedges.Get(seg.edgenr); - seg.p2 = osedges.Get(seg.edgenr); + newseg[0] = osedges.Get(seg.edgenr); + seg[1] = osedges.Get(seg.edgenr); mesh.AddSegment (newseg); } } @@ -1329,13 +1329,13 @@ namespace netgen { if (refedgesinv.Get(k)) { - seg.p1 = lastpi; - seg.p2 = thispi; + seg[0] = lastpi; + seg[1] = thispi; } else { - seg.p1 = thispi; - seg.p2 = lastpi; + seg[0] = thispi; + seg[1] = lastpi; } seg.si = refedges.Get(k).si; seg.domin = refedges.Get(k).domin; @@ -1496,13 +1496,13 @@ namespace netgen { if (refedgesinv.Get(k)) { - seg.p1 = pi1; - seg.p2 = pi2; + seg[0] = pi1; + seg[1] = pi2; } else { - seg.p1 = pi2; - seg.p2 = pi1; + seg[0] = pi2; + seg[1] = pi1; } seg.si = refedges.Get(k).si; @@ -1515,7 +1515,7 @@ namespace netgen seg.seginfo = 0; if (k == 1) seg.seginfo = (refedgesinv.Get(k)) ? 2 : 1; mesh.AddSegment (seg); - // (*testout) << "add seg " << seg.p1 << "-" << seg.p2 << endl; + // (*testout) << "add seg " << seg[0] << "-" << seg[1] << endl; } } @@ -1604,8 +1604,8 @@ namespace netgen if (oldseg.seginfo == 0) continue; - int pi1 = oldseg.p1; - int pi2 = oldseg.p2; + int pi1 = oldseg[0]; + int pi2 = oldseg[1]; int npi1 = geometry.identifications.Get(copyedgeidentification) -> GetIdentifiedPoint (mesh, pi1); @@ -1628,13 +1628,13 @@ namespace netgen if (inv) { - seg.p1 = npi1; - seg.p2 = npi2; + seg[0] = npi1; + seg[1] = npi2; } else { - seg.p1 = npi2; - seg.p2 = npi1; + seg[0] = npi2; + seg[1] = npi1; } seg.si = refedges.Get(k).si; seg.domin = refedges.Get(k).domin; @@ -1646,12 +1646,12 @@ namespace netgen seg.seginfo = 0; if (k == 1) seg.seginfo = refedgesinv.Get(k) ? 2 : 1; mesh.AddSegment (seg); - // (*testout) << "copy seg " << seg.p1 << "-" << seg.p2 << endl; + // (*testout) << "copy seg " << seg[0] << "-" << seg[1] << endl; #ifdef DEVELOP (*testout) << "copy seg, face = " << seg.si << ": " << " inv = " << inv << ", refinv = " << refedgesinv.Get(k) - << mesh.Point(seg.p1) << ", " << mesh.Point(seg.p2) << endl; + << mesh.Point(seg[0]) << ", " << mesh.Point(seg[1]) << endl; #endif } @@ -1792,10 +1792,10 @@ namespace netgen { mesh.AddPoint (p1, layer, EDGEPOINT); mesh.AddPoint (p2, layer, EDGEPOINT); - seg1.p1 = mesh.GetNP()-1; - seg1.p2 = mesh.GetNP(); - seg2.p2 = mesh.GetNP()-1; - seg2.p1 = mesh.GetNP(); + seg1[0] = mesh.GetNP()-1; + seg1[1] = mesh.GetNP(); + seg2[1] = mesh.GetNP()-1; + seg2[0] = mesh.GetNP(); seg1.geominfo[0].trignum = 1; seg1.geominfo[1].trignum = 1; seg2.geominfo[0].trignum = 1; diff --git a/libsrc/csg/extrusion.cpp b/libsrc/csg/extrusion.cpp index c9e1bbf2..1e0e2a18 100644 --- a/libsrc/csg/extrusion.cpp +++ b/libsrc/csg/extrusion.cpp @@ -7,6 +7,10 @@ namespace netgen { + Array > project1, project2; + + + void ExtrusionFace :: Init(void) { p0.SetSize(path->GetNSplines()); @@ -40,7 +44,6 @@ namespace netgen profile->GetCoeff(profile_spline_coeff); latest_point3d = -1.111e30; - } @@ -73,31 +76,24 @@ namespace netgen { profile = new LineSeg<2>(GeomPoint<2>(p[0],1), GeomPoint<2>(p[1],1)); - //(*testout) << "appending LineSeg<2> " << p[0] - // << " to " << p[1] << endl; } else if(ptype == 3) { profile = new SplineSeg3<2>(GeomPoint<2>(p[0],1), GeomPoint<2>(p[1],1), GeomPoint<2>(p[2],1)); - //(*testout) << "appending SplineSeg<3> " - // << p[0] << " -- " << p[1] << " -- " << p[2] << endl; } path = new SplineGeometry<3>; pos = const_cast< SplineGeometry<3> *>(path)->Load(raw_data,pos); - for(int i=0; i<3; i++) + for(int i = 0; i < 3; i++) { glob_z_direction(i) = raw_data[pos]; pos++; } - - //(*testout) << "read glob_z_direction " << glob_z_direction << endl; Init(); - } ExtrusionFace :: ~ExtrusionFace() @@ -128,11 +124,34 @@ namespace netgen v2.Normalize(); } - + +#define NEWJS void ExtrusionFace :: CalcProj(const Point<3> & point3d, Point<2> & point2d, int & seg, double & t) const { - if(Dist2(point3d,latest_point3d) < 1e-25*Dist2(path->GetSpline(0).StartPI(),path->GetSpline(0).EndPI())) +#ifdef NEWJS + // JS Version + double mindist = 0; + for (int i = 0; i < path->GetNSplines(); i++) + { + Point<2> hpoint2d; + double ht = CalcProj (point3d, hpoint2d, i); + double hdist = Dist2(point3d, p0[i]); + if (i == 0 || hdist < mindist) + { + seg = i; + t = ht; + mindist = hdist; + point2d = hpoint2d; + latest_seg = i; + } + } + return; +#endif + +#ifdef OLDWM + if (Dist2 (point3d, latest_point3d) < + 1e-25 * Dist2(path->GetSpline(0).StartPI(), path->GetSpline(0).EndPI())) { point2d = latest_point2d; seg = latest_seg; @@ -142,18 +161,15 @@ namespace netgen latest_point3d = point3d; - double cutdist(-1); - - + double cutdist = -1; Array mindist(path->GetNSplines()); - for(int i=0; iGetNSplines(); i++) + for(int i = 0; i < path->GetNSplines(); i++) { - - double auxcut(-1); - double auxmin(-1); + double auxcut = -1; + double auxmin = -1; if(spline3_path[i]) { @@ -162,9 +178,9 @@ namespace netgen Point<3> tanp(spline3_path[i]->TangentPoint()); double da,db,dc; - double l; Vec<3> dir = endp-startp; - l = dir.Length(); dir *= 1./l; + double l = dir.Length(); + dir *= 1./l; Vec<3> topoint = point3d - startp; double s = topoint * dir; if(s<=0) @@ -267,6 +283,7 @@ namespace netgen double thist = CalcProj(point3d,testpoint2d,i); testpoint3d = p0[i] + testpoint2d(0)*x_dir[i] + testpoint2d(1)*loc_z_dir[i]; double d = Dist2(point3d,testpoint3d); + //(*testout) << "(d="< & point3d, Point<2> & point2d, - const int seg) const + int seg) const { - double t(-1); + double t = -1; if(line_path[seg]) { @@ -306,7 +323,8 @@ namespace netgen { spline3_path[seg]->Project(point3d,p0[seg],t); - y_dir[seg] = spline3_path[seg]->GetTangent(t); y_dir[seg].Normalize(); + y_dir[seg] = spline3_path[seg]->GetTangent(t); + y_dir[seg].Normalize(); loc_z_dir[seg] = z_dir[seg]; Orthogonalize(y_dir[seg],loc_z_dir[seg]); x_dir[seg] = Cross(y_dir[seg],loc_z_dir[seg]); @@ -482,7 +500,7 @@ namespace netgen profile->Project(p2d,p2d,profile_par); p = p0[seg] + p2d(0)*x_dir[seg] + p2d(1)*loc_z_dir[seg]; - + Vec<2> tangent2d = profile->GetTangent(profile_par); profile_tangent = tangent2d(0)*x_dir[seg] + tangent2d(1)*y_dir[seg]; } @@ -505,6 +523,7 @@ namespace netgen return p0[0] + locpoint(0)*x_dir[0] + locpoint(1)*loc_z_dir[0]; } + bool ExtrusionFace :: BoxIntersectsFace(const Box<3> & box) const { @@ -515,7 +534,7 @@ namespace netgen //(*testout) << "box.Center() " << box.Center() << " projected " << center << " diam " << box.Diam() // << " dist " << Dist(box.Center(),center) << endl; - return (Dist(box.Center(),center) < 0.5*box.Diam()); + return (Dist(box.Center(),center) < 1*box.Diam()); } @@ -661,16 +680,12 @@ namespace netgen double facets) const { int n = int(facets) + 1; - - int i,j,k; - - - int nump = 0; - for(k=0; kGetNSplines(); k++) + + for(int k = 0; k < path -> GetNSplines(); k++) { - for(i=0; i<=n; i++) + for(int i = 0; i <= n; i++) { - Point<3> origin = path->GetSpline(k).GetPoint(double(i)/double(n)); + Point<3> origin = path -> GetSpline(k).GetPoint(double(i)/double(n)); if(!line_path[k]) { y_dir[k] = path->GetSpline(k).GetTangent(double(i)/double(n)); @@ -681,23 +696,22 @@ namespace netgen if(!line_path[k]) x_dir[k] = Cross(y_dir[k],loc_z_dir[k]); - for(j=0; j<=n; j++) + for(int j = 0; j <= n; j++) { Point<2> locp = profile->GetPoint(double(j)/double(n)); tas.AddPoint(origin + locp(0)*x_dir[k] + locp(1)*loc_z_dir[k]); - nump++; } } } - for(k=0; kGetNSplines(); k++) - for(i=0; iGetNSplines(); k++) + for(int i = 0; i < n; i++) + for(int j = 0; j < n; j++) { int pi = k*(n+1)*(n+1) + (n+1)*i +j; - tas.AddTriangle( TATriangle (0, pi,pi+1,pi+n+1)); - tas.AddTriangle( TATriangle (0, pi+1,pi+n+1,pi+n+2)); + tas.AddTriangle( TATriangle (0, pi,pi+1,pi+n+1) ); + tas.AddTriangle( TATriangle (0, pi+1,pi+n+1,pi+n+2) ); } } @@ -926,13 +940,13 @@ namespace netgen void Extrusion :: Reduce (const BoxSphere<3> & box) { - for(int i=0; iBoxIntersectsFace(box); } void Extrusion :: UnReduce () { - for(int i=0; i & point3d, Point<2> & point2d, - const int seg) const; + int seg) const; void CalcProj(const Point<3> & point3d, Point<2> & point2d, int & seg, double & t) const; diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index bf6ce7db..ff93c2ce 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -124,8 +124,8 @@ namespace netgen if (mesh[si].seginfo) { Box<3> hbox; - hbox.Set (mesh[mesh[si].p1]); - hbox.Add (mesh[mesh[si].p2]); + hbox.Set (mesh[mesh[si][0]]); + hbox.Add (mesh[mesh[si][1]]); segtree.Insert (hbox.PMin(), hbox.PMax(), si); } } @@ -137,8 +137,8 @@ namespace netgen if (!mesh[si].seginfo) continue; Box<3> hbox; - hbox.Set (mesh[mesh[si].p1]); - hbox.Add (mesh[mesh[si].p2]); + hbox.Set (mesh[mesh[si][0]]); + hbox.Add (mesh[mesh[si][1]]); hbox.Increase (1e-6); segtree.GetIntersecting (hbox.PMin(), hbox.PMax(), loc); @@ -148,12 +148,12 @@ namespace netgen SegmentIndex sj = loc[j]; if (sj >= si) continue; if (!mesh[si].seginfo || !mesh[sj].seginfo) continue; - if (mesh[mesh[si].p1].GetLayer() != mesh[mesh[sj].p2].GetLayer()) continue; + if (mesh[mesh[si][0]].GetLayer() != mesh[mesh[sj][1]].GetLayer()) continue; - Point<3> pi1 = mesh[mesh[si].p1]; - Point<3> pi2 = mesh[mesh[si].p2]; - Point<3> pj1 = mesh[mesh[sj].p1]; - Point<3> pj2 = mesh[mesh[sj].p2]; + Point<3> pi1 = mesh[mesh[si][0]]; + Point<3> pi2 = mesh[mesh[si][1]]; + Point<3> pj1 = mesh[mesh[sj][0]]; + Point<3> pj2 = mesh[mesh[sj][1]]; Vec<3> vi = pi2 - pi1; Vec<3> vj = pj2 - pj1; @@ -201,7 +201,7 @@ namespace netgen cout << "Intersection at " << ip << endl; geom.AddUserPoint (ip); - spoints.Append (MeshPoint (ip, mesh[mesh[si].p1].GetLayer())); + spoints.Append (MeshPoint (ip, mesh[mesh[si][0]].GetLayer())); mesh.AddPoint (ip); (*testout) << "found intersection at " << ip << endl; @@ -447,8 +447,8 @@ namespace netgen { PointGeomInfo gi; gi.trignum = k; - meshing.AddBoundaryElement (segments[si].p1 + 1 - PointIndex::BASE, - segments[si].p2 + 1 - PointIndex::BASE, + meshing.AddBoundaryElement (segments[si][0] + 1 - PointIndex::BASE, + segments[si][1] + 1 - PointIndex::BASE, gi, gi); } diff --git a/libsrc/csg/identify.cpp b/libsrc/csg/identify.cpp index 2e42117f..27478160 100644 --- a/libsrc/csg/identify.cpp +++ b/libsrc/csg/identify.cpp @@ -380,19 +380,19 @@ void PeriodicIdentification :: IdentifyFaces (class Mesh & mesh) if (seg2.si != fi2) continue; - // (*testout) << "seg1 = " << seg1.p1 << "-" << seg1.p2 << ", seg2 = " << seg2.p1 << "-" << seg2.p2; + // (*testout) << "seg1 = " << seg1[0] << "-" << seg1[1] << ", seg2 = " << seg2[0] << "-" << seg2[1]; if (side == 1) { - if (mesh.GetIdentifications().Get (seg1.p1, seg2.p1) && - mesh.GetIdentifications().Get (seg1.p2, seg2.p2)) + if (mesh.GetIdentifications().Get (seg1[0], seg2[0]) && + mesh.GetIdentifications().Get (seg1[1], seg2[1])) { foundother = 1; break; } - if (mesh.GetIdentifications().Get (seg1.p1, seg2.p2) && - mesh.GetIdentifications().Get (seg1.p2, seg2.p1)) + if (mesh.GetIdentifications().Get (seg1[0], seg2[1]) && + mesh.GetIdentifications().Get (seg1[1], seg2[0])) { foundother = 1; break; @@ -400,15 +400,15 @@ void PeriodicIdentification :: IdentifyFaces (class Mesh & mesh) } else { - if (mesh.GetIdentifications().Get (seg2.p1, seg1.p1) && - mesh.GetIdentifications().Get (seg2.p2, seg1.p2)) + if (mesh.GetIdentifications().Get (seg2[0], seg1[0]) && + mesh.GetIdentifications().Get (seg2[1], seg1[1])) { foundother = 1; break; } - if (mesh.GetIdentifications().Get (seg2.p1, seg1.p2) && - mesh.GetIdentifications().Get (seg2.p2, seg1.p1)) + if (mesh.GetIdentifications().Get (seg2[0], seg1[1]) && + mesh.GetIdentifications().Get (seg2[1], seg1[0])) { foundother = 1; break; @@ -1146,15 +1146,15 @@ void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh) if (side == 1) { - if (mesh.GetIdentifications().Get (seg1.p1, seg2.p1) && - mesh.GetIdentifications().Get (seg1.p2, seg2.p2)) + if (mesh.GetIdentifications().Get (seg1[0], seg2[0]) && + mesh.GetIdentifications().Get (seg1[1], seg2[1])) { foundother = 1; break; } - if (mesh.GetIdentifications().Get (seg1.p1, seg2.p2) && - mesh.GetIdentifications().Get (seg1.p2, seg2.p1)) + if (mesh.GetIdentifications().Get (seg1[0], seg2[1]) && + mesh.GetIdentifications().Get (seg1[1], seg2[0])) { foundother = 1; break; @@ -1162,15 +1162,15 @@ void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh) } else { - if (mesh.GetIdentifications().Get (seg2.p1, seg1.p1) && - mesh.GetIdentifications().Get (seg2.p2, seg1.p2)) + if (mesh.GetIdentifications().Get (seg2[0], seg1[0]) && + mesh.GetIdentifications().Get (seg2[1], seg1[1])) { foundother = 1; break; } - if (mesh.GetIdentifications().Get (seg2.p1, seg1.p2) && - mesh.GetIdentifications().Get (seg2.p2, seg1.p1)) + if (mesh.GetIdentifications().Get (seg2[0], seg1[1]) && + mesh.GetIdentifications().Get (seg2[1], seg1[0])) { foundother = 1; break; @@ -1229,7 +1229,7 @@ BuildSurfaceElements (Array & segs, for (int i1 = 0; i1 < segs.Size(); i1++) { const Segment & s1 = segs[i1]; - if (identmap[s1.p1] && identmap[s1.p2]) + if (identmap[s1[0]] && identmap[s1[1]]) for (int i2 = 0; i2 < i1; i2++) { const Segment & s2 = segs[i2]; @@ -1241,12 +1241,12 @@ BuildSurfaceElements (Array & segs, s2.domout == dom_nr))) continue; - if ((mesh.GetIdentifications().Get (s1.p1, s2.p2, nr) && - mesh.GetIdentifications().Get (s1.p2, s2.p1, nr)) || - (mesh.GetIdentifications().Get (s2.p1, s1.p2, nr) && - mesh.GetIdentifications().Get (s2.p2, s1.p1, nr))) + if ((mesh.GetIdentifications().Get (s1[0], s2[1], nr) && + mesh.GetIdentifications().Get (s1[1], s2[0], nr)) || + (mesh.GetIdentifications().Get (s2[0], s1[1], nr) && + mesh.GetIdentifications().Get (s2[1], s1[0], nr))) { - Element2d el(s1.p1, s1.p2, s2.p1, s2.p2); + Element2d el(s1[0], s1[1], s2[0], s2[1]); Vec<3> n = Cross (mesh[el[1]] - mesh[el[0]], mesh[el[3]] - mesh[el[0]]); @@ -1629,14 +1629,14 @@ BuildSurfaceElements (Array & segs, { const Segment & s1 = segs.Get(i1); const Segment & s2 = segs.Get(i2); - if (mesh.GetIdentifications().Get (s1.p1, s2.p2) && - mesh.GetIdentifications().Get (s1.p2, s2.p1)) + if (mesh.GetIdentifications().Get (s1[0], s2[1]) && + mesh.GetIdentifications().Get (s1[1], s2[0])) { Element2d el(QUAD); - el.PNum(1) = s1.p1; - el.PNum(2) = s1.p2; - el.PNum(3) = s2.p2; - el.PNum(4) = s2.p1; + el.PNum(1) = s1[0]; + el.PNum(2) = s1[1]; + el.PNum(3) = s2[1]; + el.PNum(4) = s2[0]; Vec<3> n = Cross (Point<3> (mesh.Point(el.PNum(2)))- Point<3> (mesh.Point(el.PNum(1))), diff --git a/libsrc/csg/singularref.cpp b/libsrc/csg/singularref.cpp index 8d65f2a1..9a37a146 100644 --- a/libsrc/csg/singularref.cpp +++ b/libsrc/csg/singularref.cpp @@ -53,7 +53,7 @@ void SingularEdge :: FindPointsOnEdge (class Mesh & mesh) for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++) { - INDEX_2 i2 (mesh[si].p1, mesh[si].p2); + INDEX_2 i2 (mesh[si][0], mesh[si][1]); /* bool onedge = 1; diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 8c606eac..63679c65 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -117,9 +117,9 @@ namespace netgen int p1 = -1, p2 = -2; if ( (*mesh)[si].domin == domnr) - { p1 = (*mesh)[si].p1; p2 = (*mesh)[si].p2; } + { p1 = (*mesh)[si][0]; p2 = (*mesh)[si][1]; } if ( (*mesh)[si].domout == domnr) - { p1 = (*mesh)[si].p2; p2 = (*mesh)[si].p1; } + { p1 = (*mesh)[si][1]; p2 = (*mesh)[si][0]; } if (p1 == -1) continue; @@ -207,11 +207,11 @@ namespace netgen for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++) { if ( (*mesh)[si].domin == domnr) - meshing.AddBoundaryElement ( (*mesh)[si].p1 + 1 - PointIndex::BASE, - (*mesh)[si].p2 + 1 - PointIndex::BASE, gi, gi); + meshing.AddBoundaryElement ( (*mesh)[si][0] + 1 - PointIndex::BASE, + (*mesh)[si][1] + 1 - PointIndex::BASE, gi, gi); if ( (*mesh)[si].domout == domnr) - meshing.AddBoundaryElement ( (*mesh)[si].p2 + 1 - PointIndex::BASE, - (*mesh)[si].p1 + 1 - PointIndex::BASE, gi, gi); + meshing.AddBoundaryElement ( (*mesh)[si][1] + 1 - PointIndex::BASE, + (*mesh)[si][0] + 1 - PointIndex::BASE, gi, gi); } diff --git a/libsrc/geom2d/spline.cpp b/libsrc/geom2d/spline.cpp index 51b43499..a34b6ea4 100644 --- a/libsrc/geom2d/spline.cpp +++ b/libsrc/geom2d/spline.cpp @@ -13,273 +13,197 @@ namespace netgen { #include "spline.hpp" - /* - template<> void SplineSeg3<2> :: Project (const Point<2> point, Point<2> & point_on_curve, double & t) const + + // just for testing (JS) + template + void ProjectTrivial (const SplineSeg3 & seg, + const Point point, Point & point_on_curve, double & t) { - double t_old = 0; + double mindist = -1; + for (int i = 0; i <= 1000; i++) + { + double ht = double(i)/1000; + Point p = seg.GetPoint(ht); + double dist = Dist2 (p, point); + if (i == 0 || dist < mindist) + { + mindist = dist; + t = ht; + } + } + point_on_curve = seg.GetPoint(t); + } + + + +template +void SplineSeg3 :: Project (const Point point, Point & point_on_curve, double & t) const +{ + double t_old = -1; + + if(proj_latest_t > 0. && proj_latest_t < 1.) + t = proj_latest_t; + else t = 0.5; + + Point phi; + Vec phip,phipp,phimp; - Point<2> phi; - Vec<2> phip,phipp,phimp; + int i=0; + + while(t > -0.5 && t < 1.5 && i<20 && fabs(t-t_old) > 1e-15 ) + { + GetDerivatives(t,phi,phip,phipp); + + t_old = t; + + phimp = phi-point; + + //t = min2(max2(t-(phip*phimp)/(phipp*phimp + phip*phip),0.),1.); + t -= (phip*phimp)/(phipp*phimp + phip*phip); + + i++; + } - int i=0; - - while(fabs(t-t_old) > 1e-8 && i<10) - { - GetDerivatives(t,phi,phip,phipp); - - t_old = t; - - phimp = phi-point; - - t = min2(max2(t-(phip*phimp)/(phipp*phimp + phip*phip),0.),1.); - - i++; - } - - if(i<10) - { - point_on_curve = GetPoint(t); - - double dist = Dist(point,point_on_curve); - - phi = GetPoint(0); - double auxdist = Dist(phi,point); - if(auxdist < dist) - { - t = 0.; - point_on_curve = phi; - dist = auxdist; - } - phi = GetPoint(1); - auxdist = Dist(phi,point); - if(auxdist < dist) - { - t = 1.; - point_on_curve = phi; - dist = auxdist; - } - - } - else - { - double t0 = 0; - double t1 = 0.5; - double t2 = 1.; - - double d0,d1,d2; - - - //(*testout) << "2d newtonersatz" << endl; - while(t2-t0 > 1e-8) - { - - phi = GetPoint(t0); d0 = Dist(phi,point); - phi = GetPoint(t1); d1 = Dist(phi,point); - phi = GetPoint(t2); d2 = Dist(phi,point); - - double a = (2.*d0 - 4.*d1 +2.*d2)/pow(t2-t0,2); - - if(a <= 0) - { - if(d0 < d2) - t2 -= 0.3*(t2-t0); - else - t0 += 0.3*(t2-t0); - - t1 = 0.5*(t2+t0); - } - else - { - double b = (d1-d0-a*(t1*t1-t0*t0))/(t1-t0); - - double auxt1 = -0.5*b/a; - - if(auxt1 < t0) - { - t2 -= 0.4*(t2-t0); - t0 = max2(0.,t0-0.1*(t2-t0)); - } - else if (auxt1 > t2) - { - t0 += 0.4*(t2-t0); - t2 = min2(1.,t2+0.1*(t2-t0)); - } - else - { - t1 = auxt1; - auxt1 = 0.25*(t2-t0); - t0 = max2(0.,t1-auxt1); - t2 = min2(1.,t1+auxt1); - } - - t1 = 0.5*(t2+t0); - } - - } - - - phi = GetPoint(t0); d0 = Dist(phi,point); - phi = GetPoint(t1); d1 = Dist(phi,point); - phi = GetPoint(t2); d2 = Dist(phi,point); - - double mind = d0; - t = t0; - if(d1 < mind) - { - t = t1; - mind = d1; - } - if(d2 < mind) - { - t = t2; - mind = d2; - } - - point_on_curve = GetPoint(t); - - } - } - - template<> void SplineSeg3<3> :: Project (const Point<3> point, Point<3> & point_on_curve, double & t) const - { - double t_old = -1; - - if(proj_latest_t > 0. && proj_latest_t < 1.) - t = proj_latest_t; - else - t = 0.5; - - Point<3> phi; - Vec<3> phip,phipp,phimp; - - int i=0; - - while(fabs(t-t_old) > 1e-8 && t > -0.5 && t < 1.5 && i<10) - { - GetDerivatives(t,phi,phip,phipp); - - t_old = t; - - phimp = phi-point; - - //t = min2(max2(t-(phip*phimp)/(phipp*phimp + phip*phip),0.),1.); - t -= (phip*phimp)/(phipp*phimp + phip*phip); - - i++; - } - - //if(i<10 && t > 0. && t < 1.) - if(i<10 && t > -0.4 && t < 1.4) - { - if(t < 0) + //if(i<10 && t > 0. && t < 1.) + if(i<20 && t > -0.4 && t < 1.4) + { + if(t < 0) + { t = 0.; - if(t > 1) + } + if(t > 1) + { t = 1.; + } - point_on_curve = GetPoint(t); + point_on_curve = SplineSeg3::GetPoint(t); - double dist = Dist(point,point_on_curve); + double dist = Dist(point,point_on_curve); - phi = GetPoint(0); - double auxdist = Dist(phi,point); - if(auxdist < dist) - { - t = 0.; - point_on_curve = phi; - dist = auxdist; - } - phi = GetPoint(1); - auxdist = Dist(phi,point); - if(auxdist < dist) - { - t = 1.; - point_on_curve = phi; - dist = auxdist; - } - } - else - { - double t0 = 0; - double t1 = 0.5; - double t2 = 1.; + phi = SplineSeg3 ::GetPoint(0); + double auxdist = Dist(phi,point); + if(auxdist < dist) + { + t = 0.; + point_on_curve = phi; + dist = auxdist; + } + phi = SplineSeg3 ::GetPoint(1); + auxdist = Dist(phi,point); + if(auxdist < dist) + { + t = 1.; + point_on_curve = phi; + dist = auxdist; + } + } + else + { + double t0 = 0; + double t1 = 0.5; + double t2 = 1.; - double d0,d1,d2; + double d0,d1,d2; - //(*testout) << "newtonersatz" << endl; - while(t2-t0 > 1e-8) - { + //(*testout) << "newtonersatz" << endl; + while(t2-t0 > 1e-8) + { - phi = GetPoint(t0); d0 = Dist(phi,point); - phi = GetPoint(t1); d1 = Dist(phi,point); - phi = GetPoint(t2); d2 = Dist(phi,point); + phi = SplineSeg3 ::GetPoint(t0); d0 = Dist(phi,point); + phi = SplineSeg3 ::GetPoint(t1); d1 = Dist(phi,point); + phi = SplineSeg3 ::GetPoint(t2); d2 = Dist(phi,point); - double a = (2.*d0 - 4.*d1 +2.*d2)/pow(t2-t0,2); + double a = (2.*d0 - 4.*d1 +2.*d2)/pow(t2-t0,2); - if(a <= 0) - { - if(d0 < d2) - t2 -= 0.3*(t2-t0); - else - t0 += 0.3*(t2-t0); + if(a <= 0) + { + if(d0 < d2) + t2 -= 0.3*(t2-t0); + else + t0 += 0.3*(t2-t0); - t1 = 0.5*(t2+t0); - } - else - { - double b = (d1-d0-a*(t1*t1-t0*t0))/(t1-t0); + t1 = 0.5*(t2+t0); + } + else + { + double b = (d1-d0-a*(t1*t1-t0*t0))/(t1-t0); - double auxt1 = -0.5*b/a; + double auxt1 = -0.5*b/a; - if(auxt1 < t0) - { - t2 -= 0.4*(t2-t0); - t0 = max2(0.,t0-0.1*(t2-t0)); - } - else if (auxt1 > t2) - { - t0 += 0.4*(t2-t0); - t2 = min2(1.,t2+0.1*(t2-t0)); - } - else - { - t1 = auxt1; - auxt1 = 0.25*(t2-t0); - t0 = max2(0.,t1-auxt1); - t2 = min2(1.,t1+auxt1); - } + if(auxt1 < t0) + { + t2 -= 0.4*(t2-t0); + t0 = max2(0.,t0-0.1*(t2-t0)); + } + else if (auxt1 > t2) + { + t0 += 0.4*(t2-t0); + t2 = min2(1.,t2+0.1*(t2-t0)); + } + else + { + t1 = auxt1; + auxt1 = 0.25*(t2-t0); + t0 = max2(0.,t1-auxt1); + t2 = min2(1.,t1+auxt1); + } - t1 = 0.5*(t2+t0); - } + t1 = 0.5*(t2+t0); + } - } + } - phi = GetPoint(t0); d0 = Dist(phi,point); - phi = GetPoint(t1); d1 = Dist(phi,point); - phi = GetPoint(t2); d2 = Dist(phi,point); + phi = SplineSeg3 ::GetPoint(t0); d0 = Dist(phi,point); + phi = SplineSeg3 ::GetPoint(t1); d1 = Dist(phi,point); + phi = SplineSeg3 ::GetPoint(t2); d2 = Dist(phi,point); - double mind = d0; - t = t0; - if(d1 < mind) - { - t = t1; - mind = d1; - } - if(d2 < mind) - { - t = t2; - mind = d2; - } + double mind = d0; + t = t0; + if(d1 < mind) + { + t = t1; + mind = d1; + } + if(d2 < mind) + { + t = t2; + mind = d2; + } - point_on_curve = GetPoint(t); - } - //(*testout) << " latest_t " << proj_latest_t << " t " << t << endl; + point_on_curve = SplineSeg3 ::GetPoint(t); + } + //(*testout) << " latest_t " << proj_latest_t << " t " << t << endl; - proj_latest_t = t; - } + proj_latest_t = t; + + /* + // test it by trivial sampling + double ht; + Point hp; + ProjectTrivial (*this, point, hp, ht); + if (fabs (t-ht) > 1e-3) + { + // if (Dist2 (point, hp) < Dist2 (point, point_on_curve)) + cout << "project is wrong" << endl; + cout << "t = " << t << ", ht = " << ht << endl; + cout << "dist org = " << Dist(point, point_on_curve) << endl; + cout << "dist trivial = " << Dist(point, hp) << endl; + } */ +} + + + template class SplineSeg3<2>; + template class SplineSeg3<3>; + + + + + + void CalcPartition (double l, double h, double r1, double r2, double ra, double elto0, Array & points) diff --git a/libsrc/geom2d/spline.hpp b/libsrc/geom2d/spline.hpp index 264a3234..de5835e0 100644 --- a/libsrc/geom2d/spline.hpp +++ b/libsrc/geom2d/spline.hpp @@ -21,24 +21,23 @@ template < int D > class GeomPoint : public Point { public: - /// refinement to point + /// refinement factor at point double refatpoint; + + /// hp-refinement bool hpref; - GeomPoint () - { ; } + /// + GeomPoint () { ; } /// - GeomPoint (double ax, double ay, double aref = 1, bool ahpref=false) - : Point (ax, ay), refatpoint(aref), hpref(ahpref) { ; } - GeomPoint (double ax, double ay, double az, double aref, bool ahpref=false) - : Point (ax, ay, az), refatpoint(aref), hpref(ahpref) { ; } GeomPoint (const Point & ap, double aref = 1, bool ahpref=false) : Point(ap), refatpoint(aref), hpref(ahpref) { ; } }; + /// base class for 2d - segment template < int D > class SplineSeg @@ -56,12 +55,14 @@ public: int copyfrom; /// perfrom anisotropic refinement (hp-refinement) to edge bool hpref_left; + /// perfrom anisotropic refinement (hp-refinement) to edge bool hpref_right; + /// calculates length of curve virtual double Length () const; /// returns point at curve, 0 <= t <= 1 virtual Point GetPoint (double t) const = 0; - /// returns a (not necessarily uniform) tangent vector for 0 <= t <= 1 + /// returns a (not necessarily unit-length) tangent vector for 0 <= t <= 1 virtual Vec GetTangent (const double t) const { cerr << "GetTangent not implemented for spline base-class" << endl; Vec dummy; return dummy;} virtual void GetDerivatives (const double t, @@ -112,14 +113,13 @@ class LineSeg : public SplineSeg { /// GeomPoint p1, p2; - //const GeomPoint &p1, &p2; public: /// LineSeg (const GeomPoint & ap1, const GeomPoint & ap2); /// virtual double Length () const; /// - virtual Point GetPoint (double t) const; + inline virtual Point GetPoint (double t) const; /// virtual Vec GetTangent (const double t) const; @@ -154,7 +154,6 @@ class SplineSeg3 : public SplineSeg { /// GeomPoint p1, p2, p3; - //const GeomPoint &p1, &p2, &p3; mutable double proj_latest_t; public: @@ -163,7 +162,7 @@ public: const GeomPoint & ap2, const GeomPoint & ap3); /// - virtual Point GetPoint (double t) const; + inline virtual Point GetPoint (double t) const; /// virtual Vec GetTangent (const double t) const; @@ -376,8 +375,8 @@ void SplineSeg :: Partition (double h, double elto0, Segment seg; seg.edgenr = segnr; seg.si = bc; // segnr; - seg.p1 = pi1; - seg.p2 = pi2; + seg[0] = pi1; + seg[1] = pi2; seg.domin = leftdom; seg.domout = rightdom; seg.epgeominfo[0].edgenr = segnr; @@ -409,6 +408,7 @@ void SplineSeg :: GetPoints (int n, Array > & points) points[i] = GetPoint(double(i) / (n-1)); } + template void SplineSeg :: PrintCoeff (ostream & ost) const { @@ -438,7 +438,7 @@ LineSeg :: LineSeg (const GeomPoint & ap1, template -Point LineSeg :: GetPoint (double t) const +inline Point LineSeg :: GetPoint (double t) const { return p1 + t * (p2 - p1); } @@ -529,150 +529,6 @@ void LineSeg :: GetRawData (Array & data) const } -template -void SplineSeg3 :: Project (const Point point, Point & point_on_curve, double & t) const -{ - double t_old = -1; - - if(proj_latest_t > 0. && proj_latest_t < 1.) - t = proj_latest_t; - else - t = 0.5; - - Point phi; - Vec phip,phipp,phimp; - - int i=0; - - while(t > -0.5 && t < 1.5 && i<20 && fabs(t-t_old) > 1e-15 ) - { - GetDerivatives(t,phi,phip,phipp); - - t_old = t; - - phimp = phi-point; - - //t = min2(max2(t-(phip*phimp)/(phipp*phimp + phip*phip),0.),1.); - t -= (phip*phimp)/(phipp*phimp + phip*phip); - - i++; - } - - //if(i<10 && t > 0. && t < 1.) - if(i<20 && t > -0.4 && t < 1.4) - { - if(t < 0) - { - t = 0.; - } - if(t > 1) - { - t = 1.; - } - - point_on_curve = GetPoint(t); - - double dist = Dist(point,point_on_curve); - - phi = GetPoint(0); - double auxdist = Dist(phi,point); - if(auxdist < dist) - { - t = 0.; - point_on_curve = phi; - dist = auxdist; - } - phi = GetPoint(1); - auxdist = Dist(phi,point); - if(auxdist < dist) - { - t = 1.; - point_on_curve = phi; - dist = auxdist; - } - } - else - { - double t0 = 0; - double t1 = 0.5; - double t2 = 1.; - - double d0,d1,d2; - - - //(*testout) << "newtonersatz" << endl; - while(t2-t0 > 1e-8) - { - - phi = GetPoint(t0); d0 = Dist(phi,point); - phi = GetPoint(t1); d1 = Dist(phi,point); - phi = GetPoint(t2); d2 = Dist(phi,point); - - double a = (2.*d0 - 4.*d1 +2.*d2)/pow(t2-t0,2); - - if(a <= 0) - { - if(d0 < d2) - t2 -= 0.3*(t2-t0); - else - t0 += 0.3*(t2-t0); - - t1 = 0.5*(t2+t0); - } - else - { - double b = (d1-d0-a*(t1*t1-t0*t0))/(t1-t0); - - double auxt1 = -0.5*b/a; - - if(auxt1 < t0) - { - t2 -= 0.4*(t2-t0); - t0 = max2(0.,t0-0.1*(t2-t0)); - } - else if (auxt1 > t2) - { - t0 += 0.4*(t2-t0); - t2 = min2(1.,t2+0.1*(t2-t0)); - } - else - { - t1 = auxt1; - auxt1 = 0.25*(t2-t0); - t0 = max2(0.,t1-auxt1); - t2 = min2(1.,t1+auxt1); - } - - t1 = 0.5*(t2+t0); - } - - } - - - phi = GetPoint(t0); d0 = Dist(phi,point); - phi = GetPoint(t1); d1 = Dist(phi,point); - phi = GetPoint(t2); d2 = Dist(phi,point); - - double mind = d0; - t = t0; - if(d1 < mind) - { - t = t1; - mind = d1; - } - if(d2 < mind) - { - t = t2; - mind = d2; - } - - point_on_curve = GetPoint(t); - } - //(*testout) << " latest_t " << proj_latest_t << " t " << t << endl; - - proj_latest_t = t; -} - @@ -686,7 +542,7 @@ SplineSeg3 :: SplineSeg3 (const GeomPoint & ap1, } template -Point SplineSeg3 :: GetPoint (double t) const +inline Point SplineSeg3 :: GetPoint (double t) const { double x, y, w; double b1, b2, b3; diff --git a/libsrc/geom2d/splinegeometry.cpp b/libsrc/geom2d/splinegeometry.cpp index 7d3a7989..a9052956 100644 --- a/libsrc/geom2d/splinegeometry.cpp +++ b/libsrc/geom2d/splinegeometry.cpp @@ -1035,11 +1035,11 @@ void SplineGeometry :: CopyEdgeMesh (int from, int to, Mesh & mesh, Point3dTr const Segment & seg = mesh.LineSegment(i); if (seg.edgenr == from) { - mappoints.Elem(seg.p1) = 1; - param.Elem(seg.p1) = seg.epgeominfo[0].dist; + mappoints.Elem(seg[0]) = 1; + param.Elem(seg[0]) = seg.epgeominfo[0].dist; - mappoints.Elem(seg.p2) = 1; - param.Elem(seg.p2) = seg.epgeominfo[1].dist; + mappoints.Elem(seg[1]) = 1; + param.Elem(seg[1]) = seg.epgeominfo[1].dist; } } @@ -1087,15 +1087,15 @@ void SplineGeometry :: CopyEdgeMesh (int from, int to, Mesh & mesh, Point3dTr Segment nseg; nseg.edgenr = to; nseg.si = splines.Get(to)->bc; - nseg.p1 = mappoints.Get(seg.p1); - nseg.p2 = mappoints.Get(seg.p2); + nseg[0] = mappoints.Get(seg[0]); + nseg[1] = mappoints.Get(seg[1]); nseg.domin = splines.Get(to)->leftdom; nseg.domout = splines.Get(to)->rightdom; nseg.epgeominfo[0].edgenr = to; - nseg.epgeominfo[0].dist = param.Get(seg.p1); + nseg.epgeominfo[0].dist = param.Get(seg[0]); nseg.epgeominfo[1].edgenr = to; - nseg.epgeominfo[1].dist = param.Get(seg.p2); + nseg.epgeominfo[1].dist = param.Get(seg[1]); mesh.AddSegment (nseg); } } diff --git a/libsrc/include/Makefile.am b/libsrc/include/Makefile.am index c9d97fe5..f1f6f9f6 100644 --- a/libsrc/include/Makefile.am +++ b/libsrc/include/Makefile.am @@ -2,7 +2,7 @@ noinst_HEADERS = acisgeom.hpp gprim.hpp meshing.hpp occgeom.hpp \ visual.hpp csg.hpp incvis.hpp myadt.hpp opti.hpp geometry2d.hpp \ linalg.hpp mydefs.hpp parallel.hpp stlgeom.hpp mystdlib.h -include_HEADERS = nginterface.h parallelinterface.hpp +include_HEADERS = nginterface.h nginterface_v2.hpp parallelinterface.hpp AM_CPPFLAGS = METASOURCES = AUTO diff --git a/libsrc/include/mydefs.hpp b/libsrc/include/mydefs.hpp index ea325728..fb01297b 100644 --- a/libsrc/include/mydefs.hpp +++ b/libsrc/include/mydefs.hpp @@ -20,7 +20,7 @@ // in the configure/make phases, with the // right version number #ifdef WIN32 -#define PACKAGE_VERSION "4.9.7-dev" +#define PACKAGE_VERSION "4.9.7" #endif diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp new file mode 100644 index 00000000..db849bf6 --- /dev/null +++ b/libsrc/include/nginterface_v2.hpp @@ -0,0 +1,29 @@ +class Ng_Element +{ +public: + NG_ELEMENT_TYPE type; + int npoints; + int nv; + int * points; + + NG_ELEMENT_TYPE GetType() const + { + return type; + } + + int GetNP() const + { + return npoints; + } + + int operator[] (int i) const + { + return points[i]-1; + } +}; + + + +template int Ng_GetNElements (); +template Ng_Element Ng_GetElement (int nr); + diff --git a/libsrc/interface/read_fnf_mesh.cpp b/libsrc/interface/read_fnf_mesh.cpp index 8e3d35d5..7b0586c9 100644 --- a/libsrc/interface/read_fnf_mesh.cpp +++ b/libsrc/interface/read_fnf_mesh.cpp @@ -298,9 +298,9 @@ namespace netgen for (int j = 0; j+2 < enums.Size(); j+=2) { Segment seg; - seg.p1 = enums[j]; - seg.p2 = enums[j+2]; - seg.pmid = enums[j+1]; + seg[0] = enums[j]; + seg[1] = enums[j+2]; + seg[2] = enums[j+1]; seg.edgenr = nr; mesh.AddSegment (seg); } diff --git a/libsrc/interface/readtetmesh.cpp b/libsrc/interface/readtetmesh.cpp index 0a8a7bca..524b0bcc 100644 --- a/libsrc/interface/readtetmesh.cpp +++ b/libsrc/interface/readtetmesh.cpp @@ -758,8 +758,8 @@ namespace netgen if((atof(version.c_str()) <= 1.999999 && (*segmentdata[i])[2] > 0) || (atof(version.c_str()) > 1.999999 && (*segmentdata[i])[2] > 0 && (*segmentdata[i])[2] < minId2D)) { - seg.p1 = (*segmentdata[i])[0]; - seg.p2 = (*segmentdata[i])[1]; + seg[0] = (*segmentdata[i])[0]; + seg[1] = (*segmentdata[i])[1]; seg.edgenr = (*segmentdata[i])[2]; seg.epgeominfo[0].edgenr = (*segmentdata[i])[2]; seg.epgeominfo[1].edgenr = (*segmentdata[i])[2]; @@ -770,8 +770,8 @@ namespace netgen seg.geominfo[1].trignum = (*segmentdata[i])[5]; mesh.AddSegment(seg); - seg.p1 = (*segmentdata[i])[1]; - seg.p2 = (*segmentdata[i])[0]; + seg[0] = (*segmentdata[i])[1]; + seg[1] = (*segmentdata[i])[0]; seg.si = (*segmentdata[i])[4]-minId2D+1; seg.surfnr1 = -1;//(*segmentdata[i])[3]; seg.surfnr2 = -1;//(*segmentdata[i])[4]; diff --git a/libsrc/interface/writetet.cpp b/libsrc/interface/writetet.cpp index 93c7074f..221f0301 100644 --- a/libsrc/interface/writetet.cpp +++ b/libsrc/interface/writetet.cpp @@ -107,7 +107,7 @@ namespace netgen for(SegmentIndex si = 0; si < mesh.GetNSeg(); si++) { const Segment & seg = mesh[si]; - INDEX_2 i2(seg.p1,seg.p2); + INDEX_2 i2(seg[0],seg[1]); i2.Sort(); if(edgenumbers.Used(i2)) continue; @@ -118,10 +118,10 @@ namespace netgen edge_ids.Append(seg.edgenr); - if(point_ids[seg.p1] == -1) - point_ids[seg.p1] = (version >= 2) ? seg.edgenr : 0; - if(point_ids[seg.p2] == -1) - point_ids[seg.p2] = (version >= 2) ? seg.edgenr : 0; + if(point_ids[seg[0]] == -1) + point_ids[seg[0]] = (version >= 2) ? seg.edgenr : 0; + if(point_ids[seg[1]] == -1) + point_ids[seg[1]] = (version >= 2) ? seg.edgenr : 0; } for(SurfaceElementIndex si = 0; si < mesh.GetNSE(); si++) diff --git a/libsrc/interface/writeuser.cpp b/libsrc/interface/writeuser.cpp index 4249179b..148f90a3 100644 --- a/libsrc/interface/writeuser.cpp +++ b/libsrc/interface/writeuser.cpp @@ -218,10 +218,10 @@ void WriteNeutralFormat (const Mesh & mesh, outfile << " "; outfile.width(8); - outfile << seg.p1; + outfile << seg[0]; outfile << " "; outfile.width(8); - outfile << seg.p2; + outfile << seg[1]; outfile << "\n"; } diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 1e981761..712a970d 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -3088,8 +3088,8 @@ namespace netgen for (i = 1; i <= mesh.GetNSeg(); i++) { const Segment & seg = mesh.LineSegment(i); - singv.Set (seg.p1); - singv.Set (seg.p2); + singv.Set (seg[0]); + singv.Set (seg[1]); } /* for ( i=1; i<= mesh.GetNSE(); i++) @@ -3110,7 +3110,7 @@ namespace netgen const Segment & seg = mesh.LineSegment(i); for (j = 0; j < 2; j++) { - int pi = (j == 0) ? seg.p1 : seg.p2; + int pi = (j == 0) ? seg[0] : seg[1]; if (bndind.Elem(pi) == 0) bndind.Elem(pi) = seg.edgenr; else if (bndind.Elem(pi) != seg.edgenr) @@ -3478,7 +3478,7 @@ namespace netgen for (i = 1; i <= nseg; i++) { Segment & seg = mesh.LineSegment (i); - INDEX_2 edge(seg.p1, seg.p2); + INDEX_2 edge(seg[0], seg[1]); edge.Sort(); if (cutedges.Used (edge)) { @@ -3488,15 +3488,15 @@ namespace netgen int newpi = cutedges.Get(edge); - nseg1.p2 = newpi; - nseg2.p1 = newpi; + nseg1[1] = newpi; + nseg2[0] = newpi; EdgePointGeomInfo newepgi; // // cerr << "move edgepoint " << newpi << " from " << mesh.Point(newpi); - PointBetween (mesh.Point (seg.p1), mesh.Point (seg.p2), + PointBetween (mesh.Point (seg[0]), mesh.Point (seg[1]), 0.5, seg.surfnr1, seg.surfnr2, seg.epgeominfo[0], seg.epgeominfo[1], mesh.Point (newpi), newepgi); diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 359b0cab..73dd7860 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -29,8 +29,8 @@ void InsertVirtualBoundaryLayer (Mesh & mesh) cout << "snr = " << snr << endl; if (snr == surfid) { - bndnodes.Set (mesh.LineSegment(i).p1); - bndnodes.Set (mesh.LineSegment(i).p2); + bndnodes.Set (mesh.LineSegment(i)[0]); + bndnodes.Set (mesh.LineSegment(i)[1]); } } for (i = 1; i <= mesh.GetNSeg(); i++) @@ -38,8 +38,8 @@ void InsertVirtualBoundaryLayer (Mesh & mesh) int snr = mesh.LineSegment(i).edgenr; if (snr != surfid) { - bndnodes.Clear (mesh.LineSegment(i).p1); - bndnodes.Clear (mesh.LineSegment(i).p2); + bndnodes.Clear (mesh.LineSegment(i)[0]); + bndnodes.Clear (mesh.LineSegment(i)[1]); } } @@ -66,8 +66,8 @@ void InsertVirtualBoundaryLayer (Mesh & mesh) int snr = mesh.LineSegment(i).edgenr; if (snr == surfid) { - int p1 = mesh.LineSegment(i).p1; - int p2 = mesh.LineSegment(i).p2; + int p1 = mesh.LineSegment(i)[0]; + int p2 = mesh.LineSegment(i)[1]; int p3 = mapto.Get (p1); if (!p3) p3 = p1; int p4 = mapto.Get (p2); diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 68f3202a..4e826dd6 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -686,8 +686,8 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, { SetThreadPercent(double(i)/mesh.GetNSeg()*100.); const Segment & seg = mesh[i]; - PointIndex pi1 = mesh[i].p1; - PointIndex pi2 = mesh[i].p2; + PointIndex pi1 = mesh[i][0]; + PointIndex pi2 = mesh[i][1]; bool swap = (pi1 > pi2); @@ -1087,7 +1087,7 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, if (info.order >= 2) { - if (mesh[info.elnr].p1 > mesh[info.elnr].p2) + if (mesh[info.elnr][0] > mesh[info.elnr][1]) xi = 1-xi; CalcEdgeShape (edgeorder[info.edgenr], 2*xi-1, &shapes(2)); } @@ -1131,7 +1131,7 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, if (info.order >= 2) { double fac = 2; - if (mesh[info.elnr].p1 > mesh[info.elnr].p2) + if (mesh[info.elnr][0] > mesh[info.elnr][1]) { xi = 1-xi; fac *= -1; @@ -1151,8 +1151,8 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, coefs.SetSize(info.ndof); - coefs[0] = Vec<3> (mesh[el.p1]); - coefs[1] = Vec<3> (mesh[el.p2]); + coefs[0] = Vec<3> (mesh[el[0]]); + coefs[1] = Vec<3> (mesh[el[1]]); if (info.order >= 2) { diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 09314b01..5ea6b239 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -1353,8 +1353,8 @@ namespace netgen case HP_SEGM: { Segment seg; - seg.p1 = hpel.pnums[0]; - seg.p2 = hpel.pnums[1]; + seg[0] = hpel.pnums[0]; + seg[1] = hpel.pnums[1]; // NOTE: only for less than 10000 elements (HACK) !!! seg.edgenr = hpel.index % 10000; seg.si = hpel.index / 10000; @@ -1593,8 +1593,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS for (int i = 1; i <= mesh.GetNSeg(); i++) if (mesh.LineSegment(i).singedge_left * levels >= act_ref) { - INDEX_2 i2 (mesh.LineSegment(i).p1, - mesh.LineSegment(i).p2); + INDEX_2 i2 (mesh.LineSegment(i)[0], + mesh.LineSegment(i)[1]); /* // before @@ -1712,8 +1712,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (seg.singedge_left * levels >= act_ref) { - INDEX_2 i2 (mesh.LineSegment(i).p1, - mesh.LineSegment(i).p2); + INDEX_2 i2 (mesh.LineSegment(i)[0], + mesh.LineSegment(i)[1]); edges.Set(i2,1); edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); @@ -1728,8 +1728,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (seg.singedge_right * levels >= act_ref) { - INDEX_2 i2 (mesh.LineSegment(i).p2, - mesh.LineSegment(i).p1); + INDEX_2 i2 (mesh.LineSegment(i)[1], + mesh.LineSegment(i)[0]); edges.Set (i2, 1); edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); @@ -1743,7 +1743,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS sing = 1; } - // (*testout) << "seg = " << ind << ", " << seg.p1 << "-" << seg.p2 << endl; + // (*testout) << "seg = " << ind << ", " << seg[0] << "-" << seg[1] << endl; if (seg.singedge_left * levels >= act_ref @@ -1751,7 +1751,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS { for (int j = 0; j < 2; j++) { - int pi = (j == 0) ? seg.p1 : seg.p2; + int pi = (j == 0) ? seg[0] : seg[1]; INDEX_3 & i3 = surfonpoint.Elem(pi); if (ind != i3.I1() && ind != i3.I2()) diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp index 8eeb7a22..9626c05d 100644 --- a/libsrc/meshing/improve2.cpp +++ b/libsrc/meshing/improve2.cpp @@ -490,7 +490,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh) SegmentIndex si; for (si = 0; si < mesh.GetNSeg(); si++) { - INDEX_2 i2(mesh[si].p1, mesh[si].p2); + INDEX_2 i2(mesh[si][0], mesh[si][1]); fixed[i2.I1()] = true; fixed[i2.I2()] = true; } diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 45f817f5..7874a08d 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -222,7 +222,7 @@ namespace netgen lock.Lock(); timestamp = NextTimeStamp(); - int maxn = max2 (s.p1, s.p2); + int maxn = max2 (s[0], s[1]); maxn += 1-PointIndex::BASE; /* @@ -234,16 +234,16 @@ namespace netgen ptyps[i] = INNERPOINT; } - if (ptyps[s.p1] > EDGEPOINT) ptyps[s.p1] = EDGEPOINT; - if (ptyps[s.p2] > EDGEPOINT) ptyps[s.p2] = EDGEPOINT; + if (ptyps[s[0]] > EDGEPOINT) ptyps[s[0]] = EDGEPOINT; + if (ptyps[s[1]] > EDGEPOINT) ptyps[s[1]] = EDGEPOINT; */ if (maxn <= points.Size()) { - if (points[s.p1].Type() > EDGEPOINT) - points[s.p1].SetType (EDGEPOINT); - if (points[s.p2].Type() > EDGEPOINT) - points[s.p2].SetType (EDGEPOINT); + if (points[s[0]].Type() > EDGEPOINT) + points[s[0]].SetType (EDGEPOINT); + if (points[s[1]].Type() > EDGEPOINT) + points[s[1]].SetType (EDGEPOINT); } /* else @@ -511,9 +511,9 @@ namespace netgen outfile.width(8); outfile << 0; outfile.width(8); - outfile << seg.p1; + outfile << seg[0]; outfile.width(8); - outfile << seg.p2; + outfile << seg[1]; outfile << " "; outfile.width(8); outfile << seg.geominfo[0].trignum; // stl dreiecke @@ -902,7 +902,7 @@ namespace netgen { Segment seg; int hi; - infile >> seg.si >> hi >> seg.p1 >> seg.p2; + infile >> seg.si >> hi >> seg[0] >> seg[1]; AddSegment (seg); } } @@ -916,7 +916,7 @@ namespace netgen { Segment seg; int hi; - infile >> seg.si >> hi >> seg.p1 >> seg.p2 + infile >> seg.si >> hi >> seg[0] >> seg[1] >> seg.geominfo[0].trignum >> seg.geominfo[1].trignum; AddSegment (seg); @@ -935,7 +935,7 @@ namespace netgen { Segment seg; int hi; - infile >> seg.si >> hi >> seg.p1 >> seg.p2 + infile >> seg.si >> hi >> seg[0] >> seg[1] >> seg.geominfo[0].trignum >> seg.geominfo[1].trignum >> seg.surfnr1 >> seg.surfnr2 @@ -1250,9 +1250,9 @@ namespace netgen { Segment seg; int hi; - infile >> seg.si >> hi >> seg.p1 >> seg.p2; - seg.p1 = seg.p1 + oldnp; - seg.p2 = seg.p2 + oldnp; + infile >> seg.si >> hi >> seg[0] >> seg[1]; + seg[0] = seg[0] + oldnp; + seg[1] = seg[1] + oldnp; AddSegment (seg); } } @@ -1266,11 +1266,11 @@ namespace netgen { Segment seg; int hi; - infile >> seg.si >> hi >> seg.p1 >> seg.p2 + infile >> seg.si >> hi >> seg[0] >> seg[1] >> seg.geominfo[0].trignum >> seg.geominfo[1].trignum; - seg.p1 = seg.p1 + oldnp; - seg.p2 = seg.p2 + oldnp; + seg[0] = seg[0] + oldnp; + seg[1] = seg[1] + oldnp; AddSegment (seg); } } @@ -1283,7 +1283,7 @@ namespace netgen { Segment seg; int hi; - infile >> seg.si >> hi >> seg.p1 >> seg.p2 + infile >> seg.si >> hi >> seg[0] >> seg[1] >> seg.geominfo[0].trignum >> seg.geominfo[1].trignum >> seg.surfnr1 >> seg.surfnr2 @@ -1298,8 +1298,8 @@ namespace netgen if(seg.surfnr1 >= 0) seg.surfnr1 = seg.surfnr1 + max_surfnr; if(seg.surfnr2 >= 0) seg.surfnr2 = seg.surfnr2 + max_surfnr; - seg.p1 = seg.p1 +oldnp; - seg.p2 = seg.p2 +oldnp; + seg[0] = seg[0] +oldnp; + seg[1] = seg[1] +oldnp; seg.edgenr = seg.edgenr + oldne; seg.epgeominfo[1].edgenr = seg.epgeominfo[1].edgenr + oldne; @@ -1458,7 +1458,7 @@ namespace netgen for (int i = 0; i < GetNSeg(); i++) { const Segment & seg = segments[i]; - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort(); boundaryedges -> Set (i2, 2); @@ -1587,7 +1587,7 @@ namespace netgen const Segment & seg = segments[i]; for (j = 1; j <= 2; j++) { - PointIndex hi = (j == 1) ? seg.p1 : seg.p2; + PointIndex hi = (j == 1) ? seg[0] : seg[1]; if (points[hi].Type() == INNERPOINT || points[hi].Type() == SURFACEPOINT) @@ -1624,7 +1624,7 @@ namespace netgen for (i = 0; i < GetNSeg(); i++) { const Segment & seg = segments[i]; - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort(); //boundaryedges -> Set (i2, 2); @@ -1968,7 +1968,7 @@ namespace netgen for (i = 1; i <= GetNSeg(); i++) { const Segment & seg = LineSegment(i); - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort(); if (!boundaryedges->Used (i2)) @@ -2009,7 +2009,7 @@ namespace netgen if (surfnr == 0 || seg.si == surfnr) { - INDEX_2 key(seg.p1, seg.p2); + INDEX_2 key(seg[0], seg[1]); INDEX_2 data(seg.si, -i); if (faceht.Used (key)) @@ -2029,7 +2029,7 @@ namespace netgen if (surfnr == 0 || seg.si == surfnr) { - INDEX_2 key(seg.p2, seg.p1); + INDEX_2 key(seg[1], seg[0]); if (!faceht.Used(key)) { cerr << "ERROR: Segment " << seg << " brother not used" << endl; @@ -2092,8 +2092,8 @@ namespace netgen if (data.I1()) // surfnr { Segment seg; - seg.p1 = i2.I1(); - seg.p2 = i2.I2(); + seg[0] = i2.I1(); + seg[1] = i2.I2(); seg.si = data.I1(); // find geomdata: @@ -2103,9 +2103,9 @@ namespace netgen const Element2d & el = SurfaceElement (data.I2()); for (k = 1; k <= el.GetNP(); k++) { - if (seg.p1 == el.PNum(k)) + if (seg[0] == el.PNum(k)) seg.geominfo[0] = el.GeomInfoPi(k); - if (seg.p2 == el.PNum(k)) + if (seg[1] == el.PNum(k)) seg.geominfo[1] = el.GeomInfoPi(k); } @@ -2121,8 +2121,8 @@ namespace netgen (*testout) << "line seg: "; } - (*testout) << seg.p1 << " - " << seg.p2 - << " len = " << Dist (Point(seg.p1), Point(seg.p2)) + (*testout) << seg[0] << " - " << seg[1] + << " len = " << Dist (Point(seg[0]), Point(seg[1])) << endl; opensegments.Append (seg); @@ -2145,14 +2145,14 @@ namespace netgen for (i = 1; i <= GetNSeg(); i++) { const Segment & seg = LineSegment (i); - ptyps.Elem(seg.p1) = EDGEPOINT; - ptyps.Elem(seg.p2) = EDGEPOINT; + ptyps.Elem(seg[0]) = EDGEPOINT; + ptyps.Elem(seg[1]) = EDGEPOINT; } for (i = 1; i <= GetNOpenSegments(); i++) { const Segment & seg = GetOpenSegment (i); - ptyps.Elem(seg.p1) = EDGEPOINT; - ptyps.Elem(seg.p2) = EDGEPOINT; + ptyps.Elem(seg[0]) = EDGEPOINT; + ptyps.Elem(seg[1]) = EDGEPOINT; } */ for (i = 1; i <= points.Size(); i++) @@ -2161,14 +2161,14 @@ namespace netgen for (i = 1; i <= GetNSeg(); i++) { const Segment & seg = LineSegment (i); - points[seg.p1].SetType(EDGEPOINT); - points[seg.p2].SetType(EDGEPOINT); + points[seg[0]].SetType(EDGEPOINT); + points[seg[1]].SetType(EDGEPOINT); } for (i = 1; i <= GetNOpenSegments(); i++) { const Segment & seg = GetOpenSegment (i); - points[seg.p1].SetType (EDGEPOINT); - points[seg.p2].SetType (EDGEPOINT); + points[seg[0]].SetType (EDGEPOINT); + points[seg[1]].SetType (EDGEPOINT); } @@ -2212,8 +2212,8 @@ namespace netgen for (i = 1; i <= GetNOpenSegments(); i++) { const Segment & seg = GetOpenSegment(i); - frontpoints.Set (seg.p1); - frontpoints.Set (seg.p2); + frontpoints.Set (seg[0]); + frontpoints.Set (seg[1]); } for (i = 1; i <= GetNSE(); i++) @@ -2534,15 +2534,15 @@ namespace netgen for (int i = 0; i < GetNSeg(); i++) { const Segment & seg = segments[i]; - const Point3d & p1 = points[seg.p1]; - const Point3d & p2 = points[seg.p2]; + const Point3d & p1 = points[seg[0]]; + const Point3d & p2 = points[seg[1]]; /* - INDEX_2 i21(seg.p1, seg.p2); - INDEX_2 i22(seg.p2, seg.p1); + INDEX_2 i21(seg[0], seg[1]); + INDEX_2 i22(seg[1], seg[0]); if (identifiedpoints) if (!identifiedpoints->Used (i21) && !identifiedpoints->Used (i22)) */ - if (!ident -> UsedSymmetric (seg.p1, seg.p2)) + if (!ident -> UsedSymmetric (seg[0], seg[1])) { lochfunc->SetH (Center (p1, p2), Dist (p1, p2)); } @@ -2648,7 +2648,7 @@ namespace netgen for (i = 1; i <= GetNSeg(); i++) { const Segment & seg = LineSegment(i); - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort(); bedges.Set (i2, 1); } @@ -2709,8 +2709,8 @@ namespace netgen for (i = 1; i <= GetNSeg(); i++) { const Segment & seg = LineSegment(i); - const Point3d & p1 = Point(seg.p1); - const Point3d & p2 = Point(seg.p2); + const Point3d & p1 = Point(seg[0]); + const Point3d & p2 = Point(seg[1]); RestrictLocalH (Center (p1, p2), Dist (p1, p2)); } @@ -2730,8 +2730,8 @@ namespace netgen linepoint.Clear(); for (i = 1; i <= nseg; i++) { - linepoint.Set (LineSegment(i).p1); - linepoint.Set (LineSegment(i).p2); + linepoint.Set (LineSegment(i)[0]); + linepoint.Set (LineSegment(i)[1]); } for (i = 1; i <= np; i++) @@ -2817,7 +2817,7 @@ namespace netgen case RESTRICTH_SEGMENT: { const Segment & seg = LineSegment(nr); - RestrictLocalHLine (Point (seg.p1), Point(seg.p2), loch); + RestrictLocalHLine (Point (seg[0]), Point(seg[1]), loch); break; } } @@ -3030,7 +3030,7 @@ namespace netgen } for (i = 0; i < segments.Size(); i++) - if (segments[i].p1 <= PointIndex::BASE-1) + if (segments[i][0] <= PointIndex::BASE-1) { segments.Delete(i); i--; @@ -3054,8 +3054,8 @@ namespace netgen for (i = 0; i < segments.Size(); i++) { const Segment & seg = segments[i]; - pused.Set (seg.p1); - pused.Set (seg.p2); + pused.Set (seg[0]); + pused.Set (seg[1]); } for (i = 0; i < openelements.Size(); i++) @@ -3121,8 +3121,8 @@ namespace netgen for (i = 0; i < segments.Size(); i++) { Segment & seg = segments[i]; - seg.p1 = op2np[seg.p1]; - seg.p2 = op2np[seg.p2]; + seg[0] = op2np[seg[0]]; + seg[1] = op2np[seg[1]]; } for (i = 1; i <= openelements.Size(); i++) diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 00474389..67d3e2a4 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -233,8 +233,8 @@ public: SegmentIndex AddSegment (const Segment & s); void DeleteSegment (int segnr) { - segments.Elem(segnr).p1 = PointIndex::BASE-1; - segments.Elem(segnr).p2 = PointIndex::BASE-1; + segments.Elem(segnr)[0] = PointIndex::BASE-1; + segments.Elem(segnr)[1] = PointIndex::BASE-1; } void FullDeleteSegment (int segnr) // von wem ist das ??? { diff --git a/libsrc/meshing/meshtool.cpp b/libsrc/meshing/meshtool.cpp index e85104ac..9ebc560e 100644 --- a/libsrc/meshing/meshtool.cpp +++ b/libsrc/meshing/meshtool.cpp @@ -640,7 +640,7 @@ namespace netgen { seg = &mesh.LineSegment(i); - of << seg->p2 << " " << seg->p1 << " " << seg->si << "\n"; + of << (*seg)[1] << " " << (*seg)[0] << " " << seg->si << "\n"; } } @@ -700,8 +700,8 @@ namespace netgen outfile << mesh2d.GetNSeg() << endl; for (i = 1; i <= mesh2d.GetNSeg(); i++) outfile << mesh2d.LineSegment(i).si << " " - << mesh2d.LineSegment(i).p1 << " " - << mesh2d.LineSegment(i).p2 << " " << endl; + << mesh2d.LineSegment(i)[0] << " " + << mesh2d.LineSegment(i)[1] << " " << endl; outfile << mesh2d.GetNSE() << endl; diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 1379be8c..4a139ad3 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -55,8 +55,8 @@ namespace netgen Segment :: Segment() { - p1 = -1; - p2 = -1; + pnums[0] = -1; + pnums[1] = -1; edgenr = -1; singedge_left = 0.; @@ -71,7 +71,7 @@ namespace netgen surfnr1 = -1; surfnr2 = -1; - pmid = -1; + pnums[2] = -1; meshdocval = 0; /* geominfo[0].trignum=-1; @@ -87,8 +87,7 @@ namespace netgen } Segment::Segment (const Segment & other) - : p1(other.p1), - p2(other.p2), + : edgenr(other.edgenr), singedge_left(other.singedge_left), singedge_right(other.singedge_right), @@ -101,10 +100,12 @@ namespace netgen surfnr1(other.surfnr1), surfnr2(other.surfnr2), epgeominfo(), - pmid(other.pmid), meshdocval(other.meshdocval), hp_elnr(other.hp_elnr) { + for (int j = 0; j < 3; j++) + pnums[j] = other.pnums[j]; + geominfo[0] = other.geominfo[0]; geominfo[1] = other.geominfo[1]; epgeominfo[0] = other.epgeominfo[0]; @@ -116,8 +117,8 @@ namespace netgen { if (&other != this) { - p1 = other.p1; - p2 = other.p2; + pnums[0] = other[0]; + pnums[1] = other[1]; edgenr = other.edgenr; singedge_left = other.singedge_left; singedge_right = other.singedge_right; @@ -132,7 +133,7 @@ namespace netgen surfnr2 = other.surfnr2; epgeominfo[0] = other.epgeominfo[0]; epgeominfo[1] = other.epgeominfo[1]; - pmid = other.pmid; + pnums[2] = other.pnums[2]; meshdocval = other.meshdocval; hp_elnr = other.hp_elnr; bcname = other.bcname; @@ -144,8 +145,8 @@ namespace netgen ostream & operator<<(ostream & s, const Segment & seg) { - s << seg.p1 << "(gi=" << seg.geominfo[0].trignum << ") - " - << seg.p2 << "(gi=" << seg.geominfo[1].trignum << ")" + s << seg[0] << "(gi=" << seg.geominfo[0].trignum << ") - " + << seg[1] << "(gi=" << seg.geominfo[1].trignum << ")" << " domin = " << seg.domin << ", domout = " << seg.domout << " si = " << seg.si << ", edgenr = " << seg.edgenr; return s; @@ -334,7 +335,7 @@ void Element2d :: Invert2() int Element2d::HasFace(const Element2d& el) const { - //nur für tets!!! hannes + //nur für tets!!! hannes for (int i = 1; i <= 3; i++) { if (PNumMod(i) == el[0] && diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 43e959be..e950a1df 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -233,7 +233,7 @@ class MeshPoint : public Point<3> public: MeshPoint () : layer(1), singular(0.), type(INNERPOINT) -{ + { #ifdef PARALLEL isghost = 0; #endif @@ -781,10 +781,12 @@ public: friend ostream & operator<<(ostream & s, const Segment & seg); + PointIndex pnums[3]; // p1, p2, pmid + /// point index 1 - PointIndex p1; + // PointIndex p1; /// point index 2 - PointIndex p2; + // PointIndex p2; /// edge nr int edgenr; /// @@ -810,7 +812,7 @@ public: /// EdgePointGeomInfo epgeominfo[2]; /// - int pmid; // for second order + // int pmid; // for second order /// int meshdocval; @@ -818,11 +820,13 @@ private: string* bcname; public: + /* PointIndex operator[] (int i) const { return (i == 0) ? p1 : p2; } PointIndex & operator[] (int i) { return (i == 0) ? p1 : p2; } + */ Segment& operator=(const Segment & other); @@ -849,7 +853,18 @@ public: return *bcname; } + int GetNP() const + { + return (pnums[2] < 0) ? 2 : 3; + } + ELEMENT_TYPE GetType() const + { + return (pnums[2] < 0) ? SEGMENT : SEGMENT3; + } + + PointIndex & operator[] (int i) { return pnums[i]; } + const PointIndex & operator[] (int i) const { return pnums[i]; } }; diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp index a04519de..b22e9e51 100644 --- a/libsrc/meshing/refine.cpp +++ b/libsrc/meshing/refine.cpp @@ -24,7 +24,7 @@ namespace netgen { const Segment & el = mesh.LineSegment(si); - INDEX_2 i2 = INDEX_2::Sort(el.p1, el.p2); + INDEX_2 i2 = INDEX_2::Sort(el[0], el[1]); PointIndex pinew; EdgePointGeomInfo ngi; @@ -37,8 +37,8 @@ namespace netgen { Point<3> pnew; - PointBetween (mesh.Point (el.p1), - mesh.Point (el.p2), 0.5, + PointBetween (mesh.Point (el[0]), + mesh.Point (el[1]), 0.5, el.surfnr1, el.surfnr2, el.epgeominfo[0], el.epgeominfo[1], pnew, ngi); @@ -54,9 +54,9 @@ namespace netgen Segment ns1 = el; Segment ns2 = el; - ns1.p2 = pinew; + ns1[1] = pinew; ns1.epgeominfo[1] = ngi; - ns2.p1 = pinew; + ns2[0] = pinew; ns2.epgeominfo[0] = ngi; mesh.LineSegment(si) = ns1; diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 6b3cd25a..4ec85074 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -30,22 +30,22 @@ namespace netgen { Segment & el = mesh.LineSegment(si); - INDEX_2 i2 = INDEX_2::Sort (el.p1, el.p2); + INDEX_2 i2 = INDEX_2::Sort (el[0], el[1]); if (between.Used(i2)) - el.pmid = between.Get(i2); + el[2] = between.Get(i2); else { Point<3> pb; EdgePointGeomInfo ngi; - PointBetween (mesh.Point (el.p1), - mesh.Point (el.p2), 0.5, + PointBetween (mesh.Point (el[0]), + mesh.Point (el[1]), 0.5, el.surfnr1, el.surfnr2, el.epgeominfo[0], el.epgeominfo[1], pb, ngi); - el.pmid = mesh.AddPoint (pb); - between.Set (i2, el.pmid); + el[2] = mesh.AddPoint (pb); + between.Set (i2, el[2]); } } diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 740cdd62..b662200a 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -120,16 +120,16 @@ void MeshTopology :: Update() for (int i = 1; i <= nseg; i++) { const Segment & seg = mesh.LineSegment(i); - cnt[seg.p1]++; - cnt[seg.p2]++; + cnt[seg[0]]++; + cnt[seg[1]]++; } vert2segment = new TABLE (cnt); for (int i = 1; i <= nseg; i++) { const Segment & seg = mesh.LineSegment(i); - vert2segment->AddSave (seg.p1, i); - vert2segment->AddSave (seg.p2, i); + vert2segment->AddSave (seg[0], i); + vert2segment->AddSave (seg[1], i); } if (buildedges) @@ -278,7 +278,7 @@ void MeshTopology :: Update() int elnr = vert2segment->Get(i,j); const Segment & el = mesh.LineSegment (elnr); - INDEX_2 edge(el.p1, el.p2); + INDEX_2 edge(el[0], el[1]); int edgedir = (edge.I1() > edge.I2()); if (edgedir) swap (edge.I1(), edge.I2()); @@ -348,7 +348,7 @@ void MeshTopology :: Update() { const Segment & el = mesh.LineSegment (i); - INDEX_2 edge(el.p1, el.p2); + INDEX_2 edge(el[0], el[1]); int edgedir = (edge.I1() > edge.I2()); if (edgedir) swap (edge.I1(), edge.I2()); diff --git a/libsrc/meshing/validate.cpp b/libsrc/meshing/validate.cpp index 01ace40f..ff11efc3 100644 --- a/libsrc/meshing/validate.cpp +++ b/libsrc/meshing/validate.cpp @@ -192,8 +192,8 @@ namespace netgen for(int i = 1; i <= mesh.GetNSeg(); i++) { const Segment & seg = mesh.LineSegment(i); - isedgepoint.Set(seg.p1); - isedgepoint.Set(seg.p2); + isedgepoint.Set(seg[0]); + isedgepoint.Set(seg[1]); } Array surfaceindex(np); diff --git a/libsrc/meshing/zrefine.cpp b/libsrc/meshing/zrefine.cpp index 00da8153..bb5a162e 100644 --- a/libsrc/meshing/zrefine.cpp +++ b/libsrc/meshing/zrefine.cpp @@ -33,7 +33,7 @@ namespace netgen const Segment & seg = mesh.LineSegment(i); if (seg.singedge_left || seg.singedge_right) { - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort(); singedges.Set (i2, 1); } @@ -224,8 +224,8 @@ namespace netgen for (j = 1; j <= 2; j++) { int pi = (j == 1) ? - mesh.LineSegment(i).p1 : - mesh.LineSegment(i).p2; + mesh.LineSegment(i)[0] : + mesh.LineSegment(i)[1]; edgesonpoint.Elem(pi)++; } } @@ -526,7 +526,7 @@ namespace netgen { const Segment & el = mesh.LineSegment(i); - INDEX_2 i2(el.p1, el.p2); + INDEX_2 i2(el[0], el[1]); i2.Sort(); int pnew; @@ -544,13 +544,13 @@ namespace netgen // Point3d pb; // /* - // geom->PointBetween (mesh.Point (el.p1), - // mesh.Point (el.p2), + // geom->PointBetween (mesh.Point (el[0]), + // mesh.Point (el[1]), // el.surfnr1, el.surfnr2, // el.epgeominfo[0], el.epgeominfo[1], // pb, ngi); // */ - // pb = Center (mesh.Point (el.p1), mesh.Point (el.p2)); + // pb = Center (mesh.Point (el[0]), mesh.Point (el[1])); // pnew = mesh.AddPoint (pb); @@ -563,9 +563,9 @@ namespace netgen Segment ns1 = el; Segment ns2 = el; - ns1.p2 = pnew; + ns1[1] = pnew; ns1.epgeominfo[1] = ngi; - ns2.p1 = pnew; + ns2[0] = pnew; ns2.epgeominfo[0] = ngi; mesh.LineSegment(i) = ns1; diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index 4baa4618..b7706348 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -311,8 +311,8 @@ namespace netgen edgenr++; Segment seg; - seg.p1 = pnums[i-1]; - seg.p2 = pnums[i]; + seg[0] = pnums[i-1]; + seg[1] = pnums[i]; seg.edgenr = edgenr; seg.si = facenr; seg.epgeominfo[0].dist = params[i-1]; @@ -362,7 +362,7 @@ namespace netgen if (edge.Orientation() == TopAbs_REVERSED) { - swap (seg.p1, seg.p2); + swap (seg[0], seg[1]); swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist); swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); @@ -379,7 +379,7 @@ namespace netgen // for(i=1; i<=mesh.GetNSeg(); i++) // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si - // << " p1 " << mesh.LineSegment(i).p1 << " p2 " << mesh.LineSegment(i).p2 << endl; + // << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl; // exit(10); mesh.CalcSurfacesOfNode(); @@ -471,7 +471,7 @@ namespace netgen { for (j = 1; j <= 2; j++) { - int pi = (j == 1) ? seg.p1 : seg.p2; + int pi = (j == 1) ? seg[0] : seg[1]; if (!glob2loc.Get(pi)) { meshing.AddPoint (mesh.Point(pi), pi); @@ -494,7 +494,7 @@ namespace netgen gi1.u = seg.epgeominfo[1].u; gi1.v = seg.epgeominfo[1].v; - meshing.AddBoundaryElement (glob2loc.Get(seg.p1), glob2loc.Get(seg.p2), gi0, gi1); + meshing.AddBoundaryElement (glob2loc.Get(seg[0]), glob2loc.Get(seg[1]), gi0, gi1); //(*testout) << gi0.u << " " << gi0.v << endl; //(*testout) << gi1.u << " " << gi1.v << endl; } @@ -543,7 +543,7 @@ namespace netgen if (locpnum[j] == 0) { - int pi = (j == 0) ? seg.p1 : seg.p2; + int pi = (j == 0) ? seg[0] : seg[1]; meshing.AddPoint (mesh.Point(pi), pi); gis.SetSize (gis.Size()+1); @@ -1417,8 +1417,8 @@ namespace netgen for (j = 1; j <= mesh->GetNSeg(); j++) { Segment & seg = mesh->LineSegment(j); - if (seg.p1 == i) seg.p1 = equalto[i-1]; - if (seg.p2 == i) seg.p2 = equalto[i-1]; + if (seg[0] == i) seg[0] = equalto[i-1]; + if (seg[1] == i) seg[1] = equalto[i-1]; } } @@ -1426,7 +1426,7 @@ namespace netgen for (j = 1; j <= mesh->GetNSeg(); j++) { Segment & seg = mesh->LineSegment(j); - if (seg.p1 == seg.p2) + if (seg[0] == seg[1]) { mesh->DeleteSegment(j); cout << "Deleting Segment " << j << endl; diff --git a/libsrc/stlgeom/meshstlsurface.cpp b/libsrc/stlgeom/meshstlsurface.cpp index 89264c98..c385a964 100644 --- a/libsrc/stlgeom/meshstlsurface.cpp +++ b/libsrc/stlgeom/meshstlsurface.cpp @@ -118,8 +118,8 @@ static void STLFindEdges (STLGeometry & geom, */ Point3d hp, hp2; Segment seg; - seg.p1 = p1; - seg.p2 = p2; + seg[0] = p1; + seg[1] = p2; seg.si = geom.GetTriangle(trig1).GetFaceNum(); seg.edgenr = i; @@ -140,7 +140,7 @@ static void STLFindEdges (STLGeometry & geom, /* geom.SelectChartOfTriangle (trig1); - hp = hp2 = mesh.Point (seg.p1); + hp = hp2 = mesh.Point (seg[0]); seg.geominfo[0].trignum = geom.Project (hp); (*testout) << "hp = " << hp2 << ", hp proj = " << hp << ", trignum = " << seg.geominfo[0].trignum << endl; @@ -150,7 +150,7 @@ static void STLFindEdges (STLGeometry & geom, } geom.SelectChartOfTriangle (trig1b); - hp = hp2 = mesh.Point (seg.p2); + hp = hp2 = mesh.Point (seg[1]); seg.geominfo[1].trignum = geom.Project (hp); (*testout) << "hp = " << hp2 << ", hp proj = " << hp << ", trignum = " << seg.geominfo[1].trignum << endl; @@ -161,12 +161,12 @@ static void STLFindEdges (STLGeometry & geom, */ - if (Dist (mesh.Point(seg.p1), mesh.Point(seg.p2)) < 1e-10) + if (Dist (mesh.Point(seg[0]), mesh.Point(seg[1])) < 1e-10) { (*testout) << "ERROR: Line segment of length 0" << endl; - (*testout) << "pi1, 2 = " << seg.p1 << ", " << seg.p2 << endl; - (*testout) << "p1, 2 = " << mesh.Point(seg.p1) - << ", " << mesh.Point(seg.p2) << endl; + (*testout) << "pi1, 2 = " << seg[0] << ", " << seg[1] << endl; + (*testout) << "p1, 2 = " << mesh.Point(seg[0]) + << ", " << mesh.Point(seg[1]) << endl; throw NgException ("Line segment of length 0"); } @@ -174,8 +174,8 @@ static void STLFindEdges (STLGeometry & geom, Segment seg2; - seg2.p1 = p2; - seg2.p2 = p1; + seg2[0] = p2; + seg2[1] = p1; seg2.si = geom.GetTriangle(trig2).GetFaceNum(); seg2.edgenr = i; @@ -196,7 +196,7 @@ static void STLFindEdges (STLGeometry & geom, /* geom.SelectChartOfTriangle (trig2); - hp = hp2 = mesh.Point (seg.p1); + hp = hp2 = mesh.Point (seg[0]); seg2.geominfo[0].trignum = geom.Project (hp); (*testout) << "hp = " << hp2 << ", hp proj = " << hp << ", trignum = " << seg.geominfo[0].trignum << endl; @@ -207,7 +207,7 @@ static void STLFindEdges (STLGeometry & geom, geom.SelectChartOfTriangle (trig2b); - hp = hp2 = mesh.Point (seg.p2); + hp = hp2 = mesh.Point (seg[1]); seg2.geominfo[1].trignum = geom.Project (hp); (*testout) << "hp = " << hp2 << ", hp proj = " << hp << ", trignum = " << seg.geominfo[1].trignum << endl; if (Dist (hp, hp2) > 1e-5 || seg2.geominfo[1].trignum == 0) @@ -301,7 +301,7 @@ int STLSurfaceMeshing (STLGeometry & geom, for (i = 1; i <= nopen; i++) { const Segment & seg = mesh.GetOpenSegment (i); - geom.AddMarkedSeg(mesh.Point(seg.p1),mesh.Point(seg.p2)); + geom.AddMarkedSeg(mesh.Point(seg[0]),mesh.Point(seg[1])); } geom.InitMarkedTrigs(); @@ -364,7 +364,7 @@ int STLSurfaceMeshing (STLGeometry & geom, for (i = 1; i <= mesh.GetNOpenSegments(); i++) { const Segment & seg = mesh.GetOpenSegment (i); - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort(); openseght.Set (i2, 1); } @@ -382,12 +382,12 @@ int STLSurfaceMeshing (STLGeometry & geom, for (i = 1; i <= nsegold; i++) { Segment seg = mesh.LineSegment(i); - INDEX_2 i2(seg.p1, seg.p2); + INDEX_2 i2(seg[0], seg[1]); i2.Sort(); if (openseght.Used (i2)) { // segment will be split - PrintMessage(7,"Split segment ", int(seg.p1), "-", int(seg.p2)); + PrintMessage(7,"Split segment ", int(seg[0]), "-", int(seg[1])); Segment nseg1, nseg2; EdgePointGeomInfo newgi; @@ -418,23 +418,23 @@ int STLSurfaceMeshing (STLGeometry & geom, nseg1 = seg; nseg2 = seg; - nseg1.p2 = newpi; + nseg1[1] = newpi; nseg1.epgeominfo[1] = newgi; - nseg2.p1 = newpi; + nseg2[0] = newpi; nseg2.epgeominfo[0] = newgi; mesh.LineSegment(i) = nseg1; mesh.AddSegment (nseg2); - mesh.RestrictLocalH (Center (mesh.Point(nseg1.p1), - mesh.Point(nseg1.p2)), - Dist (mesh.Point(nseg1.p1), - mesh.Point(nseg1.p2))); - mesh.RestrictLocalH (Center (mesh.Point(nseg2.p1), - mesh.Point(nseg2.p2)), - Dist (mesh.Point(nseg2.p1), - mesh.Point(nseg2.p2))); + mesh.RestrictLocalH (Center (mesh.Point(nseg1[0]), + mesh.Point(nseg1[1])), + Dist (mesh.Point(nseg1[0]), + mesh.Point(nseg1[1]))); + mesh.RestrictLocalH (Center (mesh.Point(nseg2[0]), + mesh.Point(nseg2[1])), + Dist (mesh.Point(nseg2[0]), + mesh.Point(nseg2[1]))); } } @@ -657,7 +657,7 @@ void STLSurfaceMeshing1 (STLGeometry & geom, { const Segment & seg = mesh.GetOpenSegment (i); if (seg.si == fnr) - meshing.AddBoundaryElement (seg.p1, seg.p2, seg.geominfo[0], seg.geominfo[1]); + meshing.AddBoundaryElement (seg[0], seg[1], seg.geominfo[0], seg.geominfo[1]); } diff --git a/libsrc/visualization/meshdoc.cpp b/libsrc/visualization/meshdoc.cpp index defc39b3..82c2af26 100644 --- a/libsrc/visualization/meshdoc.cpp +++ b/libsrc/visualization/meshdoc.cpp @@ -423,11 +423,11 @@ void VisualSceneMeshDoctor :: BuildScene (int zoomall) for (i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); - const Point3d & p1 = mesh->Point(seg.p1); - const Point3d & p2 = mesh->Point(seg.p2); + const Point3d & p1 = mesh->Point(seg[0]); + const Point3d & p2 = mesh->Point(seg[1]); - if (edgedist.Get(seg.p1) <= markedgedist && - edgedist.Get(seg.p2) <= markedgedist) + if (edgedist.Get(seg[0]) <= markedgedist && + edgedist.Get(seg[1]) <= markedgedist) { glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcolseledge); @@ -574,8 +574,8 @@ void VisualSceneMeshDoctor :: UpdateTables () for (i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); - if (seg.p1 == selpoint && seg.p2 == selpoint2 || - seg.p2 == selpoint && seg.p1 == selpoint2) + if (seg[0] == selpoint && seg[1] == selpoint2 || + seg[1] == selpoint && seg[0] == selpoint2) { edgedist.Elem(selpoint) = 1; edgedist.Elem(selpoint2) = 1; @@ -590,17 +590,17 @@ void VisualSceneMeshDoctor :: UpdateTables () { const Segment & seg = mesh->LineSegment(i); - int edist = min2 (edgedist.Get(seg.p1), edgedist.Get(seg.p2)); + int edist = min2 (edgedist.Get(seg[0]), edgedist.Get(seg[1])); edist++; - if (edgedist.Get(seg.p1) > edist) + if (edgedist.Get(seg[0]) > edist) { - edgedist.Elem(seg.p1) = edist; + edgedist.Elem(seg[0]) = edist; changed = 1; } - if (edgedist.Get(seg.p2) > edist) + if (edgedist.Get(seg[1]) > edist) { - edgedist.Elem(seg.p2) = edist; + edgedist.Elem(seg[1]) = edist; changed = 1; } } @@ -611,8 +611,8 @@ void VisualSceneMeshDoctor :: UpdateTables () int VisualSceneMeshDoctor :: IsSegmentMarked (int segnr) const { const Segment & seg = mesh->LineSegment(segnr); - return (edgedist.Get(seg.p1) <= markedgedist && - edgedist.Get(seg.p2) <= markedgedist); + return (edgedist.Get(seg[0]) <= markedgedist && + edgedist.Get(seg[1]) <= markedgedist); } } diff --git a/libsrc/visualization/mvdraw.cpp b/libsrc/visualization/mvdraw.cpp index 62b3cfaf..58b4173e 100644 --- a/libsrc/visualization/mvdraw.cpp +++ b/libsrc/visualization/mvdraw.cpp @@ -1282,10 +1282,10 @@ namespace netgen for (int i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh -> LineSegment (i); - glVertex3dv ( (*mesh)[seg.p1] ); - glVertex3dv ( (*mesh)[seg.p2] ); - // glVertex3dv ( &(*mesh)[seg.p1].X() ); - // glVertex3dv ( &(*mesh)[seg.p2].X() ); + glVertex3dv ( (*mesh)[seg[0]] ); + glVertex3dv ( (*mesh)[seg[1]] ); + // glVertex3dv ( &(*mesh)[seg[0]].X() ); + // glVertex3dv ( &(*mesh)[seg[1]].X() ); } glEnd(); } @@ -1316,8 +1316,8 @@ namespace netgen for (int i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh -> LineSegment (i); - const Point3d p1 = mesh -> Point (seg.p1); - const Point3d p2 = mesh -> Point (seg.p2); + const Point3d p1 = mesh -> Point (seg[0]); + const Point3d p2 = mesh -> Point (seg[1]); const Point3d p = Center (p1, p2); glRasterPos3d (p.X(), p.Y(), p.Z()); @@ -1421,11 +1421,11 @@ namespace netgen if (mesh->GetNSeg()) { - box.SetPoint (mesh->Point (mesh->LineSegment(1).p1)); + box.SetPoint (mesh->Point (mesh->LineSegment(1)[0])); for (int i = 1; i <= mesh->GetNSeg(); i++) { - box.AddPoint (mesh->Point (mesh->LineSegment(i).p1)); - box.AddPoint (mesh->Point (mesh->LineSegment(i).p2)); + box.AddPoint (mesh->Point (mesh->LineSegment(i)[0])); + box.AddPoint (mesh->Point (mesh->LineSegment(i)[1])); } } else if (specpoints.Size() >= 2) diff --git a/libsrc/visualization/vscsg.cpp b/libsrc/visualization/vscsg.cpp index 4d84d202..af379043 100644 --- a/libsrc/visualization/vscsg.cpp +++ b/libsrc/visualization/vscsg.cpp @@ -15,7 +15,7 @@ namespace netgen /* *********************** Draw Geometry **************** */ - + extern Array > project1, project2; extern AutoPtr geometry; @@ -117,6 +117,23 @@ void VisualSceneGeometry :: DrawScene () glDisable (GL_POLYGON_OFFSET_FILL); + /* + cout << "draw " << project1.Size() << " lines " << endl; + glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); + glLineWidth (1.0f); + glEnable (GL_COLOR_MATERIAL); + + glColor3f (1.0f, 0.0f, 0.0f); + + glBegin (GL_LINES); + for (int i = 0; i < project1.Size(); i++) + { + glVertex3dv (project1[i]); + glVertex3dv (project2[i]); + } + glEnd(); + */ + glPopMatrix(); glDisable(GL_CLIP_PLANE0); @@ -236,16 +253,11 @@ void VisualSceneGeometry :: BuildScene (int zoomall) glBegin (GL_TRIANGLES); for (int j = 0; j < ta.GetNT(); j++) { - for (int k = 0; k < 3; k++) { int pi = ta.GetTriangle(j)[k]; - glNormal3f (ta.GetNormal (pi)(0), - ta.GetNormal (pi)(1), - ta.GetNormal (pi)(2)); - glVertex3f (ta.GetPoint(pi)(0), - ta.GetPoint(pi)(1), - ta.GetPoint(pi)(2)); + glNormal3dv (ta.GetNormal (pi)); + glVertex3dv (ta.GetPoint(pi)); } } glEnd (); diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index 3c06ecb6..1d4b26f7 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -459,8 +459,8 @@ namespace netgen { const Segment & seg = (*mesh)[i]; - const Point3d & p1 = mesh->Point(seg.p1); - const Point3d & p2 = mesh->Point(seg.p2); + const Point3d & p1 = mesh->Point(seg[0]); + const Point3d & p2 = mesh->Point(seg[1]); const Point3d p = Center (p1, p2); glRasterPos3d (p.X(), p.Y(), p.Z()); @@ -1650,8 +1650,8 @@ namespace netgen for (int i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); - const Point3d & p1 = (*mesh)[seg.p1]; - const Point3d & p2 = (*mesh)[seg.p2]; + const Point3d & p1 = (*mesh)[seg[0]]; + const Point3d & p2 = (*mesh)[seg[1]]; if (seg.singedge_left || seg.singedge_right) glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, @@ -3239,8 +3239,8 @@ namespace netgen for (i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); - if (seg.p1 == selpoint && seg.p2 == selpoint2 || - seg.p2 == selpoint && seg.p1 == selpoint2) + if (seg[0] == selpoint && seg[1] == selpoint2 || + seg[1] == selpoint && seg[0] == selpoint2) { seledge = seg.edgenr; cout << "seledge = " << seledge << endl; @@ -3418,8 +3418,8 @@ namespace netgen for (i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); - if (seg.p1 == selpoint && seg.p2 == selpoint2 || - seg.p2 == selpoint && seg.p1 == selpoint2) + if (seg[0] == selpoint && seg[1] == selpoint2 || + seg[1] == selpoint && seg[0] == selpoint2) { seledge = seg.edgenr; cout << "seledge = " << seledge << endl; diff --git a/ng/Makefile.am b/ng/Makefile.am index 4d884a29..fec327c9 100644 --- a/ng/Makefile.am +++ b/ng/Makefile.am @@ -3,7 +3,7 @@ include_HEADERS = AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include -I$(top_srcdir)/libsrc/interface -DOPENGL -D$(TOGL_WINDOWINGSYSTEM) $(OCCFLAGS) $(TCL_INCLUDES) $(MPI_INCLUDES) bin_PROGRAMS = netgen -netgen_SOURCES = demoview.cpp ngappinit.cpp ngpkg.cpp onetcl.cpp nginterface.cpp parallelfunc.cpp parallelinterface.cpp demoview.hpp parallelfunc.hpp togl_1_7.h +netgen_SOURCES = demoview.cpp ngappinit.cpp ngpkg.cpp onetcl.cpp nginterface.cpp nginterface_v2.cpp parallelfunc.cpp parallelinterface.cpp demoview.hpp parallelfunc.hpp togl_1_7.h netgen_LDADD = $(top_builddir)/libsrc/visualization/libvisual.a \ @@ -17,8 +17,6 @@ netgen_LDADD = $(top_builddir)/libsrc/visualization/libvisual.a \ $(top_builddir)/libsrc/linalg/libla.la \ $(top_builddir)/libsrc/general/libgeneral.la \ $(OCCLIBS) -L$(TK_BIN_DIR)/Togl1.7 $(TOGLLIBDIR) -lTogl1.7 $(LIBGLU) $(TK_LIB_SPEC) $(TCL_LIB_SPEC) $(MPI_LIBS) $(PKG_LIBS) -#-lGLU -lGL - dist_bin_SCRIPTS = dialog.tcl menustat.tcl ngicon.tcl ng.tcl \ ngvisual.tcl sockets.tcl drawing.tcl nghelp.tcl ngshell.tcl \ diff --git a/ng/nginterface.cpp b/ng/nginterface.cpp index 1c0ee70c..a84ecf98 100644 --- a/ng/nginterface.cpp +++ b/ng/nginterface.cpp @@ -526,19 +526,19 @@ NG_ELEMENT_TYPE Ng_GetSurfaceElement (int ei, int * epi, int * np) { const Segment & seg = mesh->LineSegment (ei); - if (seg.pmid < 0) + if (seg[2] < 0) { - epi[0] = seg.p1; - epi[1] = seg.p2; + epi[0] = seg[0]; + epi[1] = seg[1]; if (np) *np = 2; return NG_SEGM; } else { - epi[0] = seg.p1; - epi[1] = seg.p2; - epi[2] = seg.pmid; + epi[0] = seg[0]; + epi[1] = seg[1]; + epi[2] = seg[2]; if (np) *np = 3; return NG_SEGM3; @@ -953,17 +953,17 @@ NG_ELEMENT_TYPE Ng_GetSegment (int ei, int * epi, int * np) { const Segment & seg = mesh->LineSegment (ei); - epi[0] = seg.p1; - epi[1] = seg.p2; + epi[0] = seg[0]; + epi[1] = seg[1]; - if (seg.pmid < 0) + if (seg[2] < 0) { if (np) *np = 2; return NG_SEGM; } else { - epi[2] = seg.pmid; + epi[2] = seg[2]; if (np) *np = 3; return NG_SEGM3; } @@ -2024,8 +2024,8 @@ int Ng_GetNPeriodicEdges (int idnr) for (SegmentIndex si = 0; si < nse; si++) { - PointIndex other1 = map[(*mesh)[si].p1]; - PointIndex other2 = map[(*mesh)[si].p2]; + PointIndex other1 = map[(*mesh)[si][0]]; + PointIndex other2 = map[(*mesh)[si][1]]; // (*testout) << "seg = " << (*mesh)[si] << "; other = " // << other1 << "-" << other2 << endl; if (other1 && other2 && mesh->IsSegment (other1, other2)) @@ -2052,8 +2052,8 @@ void Ng_GetPeriodicEdges (int idnr, int * pairs) for (SegmentIndex si = 0; si < nse; si++) { - PointIndex other1 = map[(*mesh)[si].p1]; - PointIndex other2 = map[(*mesh)[si].p2]; + PointIndex other1 = map[(*mesh)[si][0]]; + PointIndex other2 = map[(*mesh)[si][1]]; if (other1 && other2 && mesh->IsSegment (other1, other2)) { SegmentIndex otherseg = mesh->SegmentNr (other1, other2); diff --git a/ng/nginterface_v2.cpp b/ng/nginterface_v2.cpp new file mode 100644 index 00000000..579a319b --- /dev/null +++ b/ng/nginterface_v2.cpp @@ -0,0 +1,128 @@ +#include + + +#include +#include +#include +#include + + +#ifdef OCCGEOMETRY +#include +#endif + +#ifdef ACIS +#include +#endif + +#ifdef SOCKETS +#include "../sockets/sockets.hpp" +#endif + +#ifndef NOTCL +#include +#endif + + +#include "nginterface.h" +#include "nginterface_v2.hpp" + +// #include + + +// #include + + +namespace netgen +{ +#include "writeuser.hpp" + + extern AutoPtr mesh; +#ifndef NOTCL + extern VisualSceneMesh vsmesh; + extern Tcl_Interp * tcl_interp; +#endif + + extern AutoPtr geometry2d; + extern AutoPtr geometry; + extern STLGeometry * stlgeometry; + +#ifdef OCCGEOMETRY + extern OCCGeometry * occgeometry; +#endif +#ifdef ACIS + extern ACISGeometry * acisgeometry; +#endif + +#ifdef OPENGL + extern VisualSceneSolution vssolution; +#endif + extern CSGeometry * ParseCSG (istream & istr); + +#ifdef SOCKETS + extern AutoPtr clientsocket; + //extern Array< AutoPtr < ServerInfo > > servers; + extern Array< ServerInfo* > servers; +#endif +} + + +using namespace netgen; + + + +template <> int Ng_GetNElements<1> () +{ + return mesh->GetNSeg(); +} + +template <> int Ng_GetNElements<2> () +{ + return mesh->GetNSE(); +} + +template <> int Ng_GetNElements<3> () +{ + return mesh->GetNE(); +} + + + + +template <> Ng_Element Ng_GetElement<1> (int nr) +{ + const Segment & el = mesh->LineSegment (SegmentIndex(nr)); + + Ng_Element ret; + ret.type = NG_ELEMENT_TYPE(el.GetType()); + ret.npoints = el.GetNP(); + ret.points = (int*)&(el[0]); + + return ret; +} + +template <> +Ng_Element Ng_GetElement<2> (int nr) +{ + const Element2d & el = mesh->SurfaceElement (SurfaceElementIndex (nr)); + + Ng_Element ret; + ret.type = NG_ELEMENT_TYPE(el.GetType()); + ret.npoints = el.GetNP(); + ret.points = (int*)&el[0]; + return ret; +} + +template <> +Ng_Element Ng_GetElement<3> (int nr) +{ + const Element & el = mesh->VolumeElement (ElementIndex (nr)); + + Ng_Element ret; + ret.type = NG_ELEMENT_TYPE(el.GetType()); + ret.npoints = el.GetNP(); + ret.points = (int*)&el[0]; + return ret; +} + + diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index fa7e1b2d..b6cd784f 100644 --- a/nglib/nglib.cpp +++ b/nglib/nglib.cpp @@ -239,8 +239,8 @@ DLL_HEADER void Ng_AddBoundarySeg_2D (Ng_Mesh * mesh, int pi1, int pi2) Mesh * m = (Mesh*)mesh; Segment seg; - seg.p1 = pi1; - seg.p2 = pi2; + seg[0] = pi1; + seg[1] = pi2; m->AddSegment (seg); } @@ -285,8 +285,8 @@ DLL_HEADER void Ng_GetElement_2D (Ng_Mesh * mesh, int num, int * pi, int * matnu DLL_HEADER void Ng_GetSegment_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum) { const Segment & seg = ((Mesh*)mesh)->LineSegment(num); - pi[0] = seg.p1; - pi[1] = seg.p2; + pi[0] = seg[0]; + pi[1] = seg[1]; if (matnum) *matnum = seg.edgenr;