From 2584903baa6c5bb1acc778f42633587eba65b482 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 19 Apr 2009 21:15:26 +0000 Subject: [PATCH] extrusion fixes, reduce warnings --- libsrc/csg/csgeom.cpp | 4 +- libsrc/csg/csgeom.hpp | 2 - libsrc/csg/extrusion.cpp | 113 +- libsrc/csg/extrusion.hpp | 8 + libsrc/csg/identify.cpp | 34 +- libsrc/csg/revolution.cpp | 2 +- libsrc/csg/singularref.cpp | 14 +- libsrc/csg/specpoin.cpp | 2 +- libsrc/geom2d/genmesh2d.cpp | 9 +- libsrc/geom2d/splinegeometry.cpp | 8 +- libsrc/gprim/geomobjects.hpp | 2 +- libsrc/gprim/geomtest3d.cpp | 2 +- libsrc/interface/readtetmesh.cpp | 2 +- libsrc/interface/readuser.cpp | 33 +- libsrc/interface/writedolfin.cpp | 12 +- libsrc/interface/writeelmer.cpp | 3 +- libsrc/interface/wuchemnitz.cpp | 494 +- libsrc/meshing/bisect.cpp | 16 +- libsrc/meshing/classifyhpel.hpp | 4 + libsrc/meshing/curvedelems.cpp | 15 +- libsrc/meshing/findip.hpp | 2 +- libsrc/meshing/meshclass.cpp | 7357 +++++++++++++-------------- libsrc/meshing/meshing2.cpp | 13 +- libsrc/stlgeom/stlgeommesh.cpp | 2 +- libsrc/stlgeom/stltool.cpp | 6 +- libsrc/visualization/vssolution.cpp | 199 +- 26 files changed, 4192 insertions(+), 4166 deletions(-) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index 19b2ffb2..7e74bb51 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -53,14 +53,14 @@ namespace netgen CSGeometry :: CSGeometry () : boundingbox (default_boundingbox), - identicsurfaces (100), filename(""), ideps(1e-9) + identicsurfaces (100), ideps(1e-9), filename("") { ; } CSGeometry :: CSGeometry (const string & afilename) : boundingbox (default_boundingbox), - identicsurfaces (100), filename(afilename), ideps(1e-9) + identicsurfaces (100), ideps(1e-9), filename(afilename) { changeval++; } diff --git a/libsrc/csg/csgeom.hpp b/libsrc/csg/csgeom.hpp index a96672b0..a1927cb3 100644 --- a/libsrc/csg/csgeom.hpp +++ b/libsrc/csg/csgeom.hpp @@ -134,11 +134,9 @@ private: double ideps; - /// filename of inputfile string filename; - public: CSGeometry (); CSGeometry (const string & afilename); diff --git a/libsrc/csg/extrusion.cpp b/libsrc/csg/extrusion.cpp index abb3a6a5..9739fcf6 100644 --- a/libsrc/csg/extrusion.cpp +++ b/libsrc/csg/extrusion.cpp @@ -356,25 +356,27 @@ namespace netgen int seg; CalcProj (point, p2d, seg, t_path); + Point<3> phi; Vec<3> phip, phipp, phi_minus_point; path->GetSpline(seg).GetDerivatives(t_path, phi, phip, phipp); phi_minus_point = phi-point; - Vec<3> grad_t = (1.0/(phipp*phi_minus_point + phip*phip)) * phip; + Vec<3> grad_xbar, grad_ybar; + /* Vec<3> dphi_dX[3]; - Vec<3> dy_dir_dX[3]; - Vec<3> dx_dir_dX[3]; - for(int i = 0; i < 3; i++) dphi_dX[i] = grad_t(i)*phip; + Vec<3> dx_dir_dX[3]; + Vec<3> dy_dir_dX[3]; + Vec<3> dz_dir_dX[3]; + double lphip = phip.Length(); Vec<3> dy_dir_dt = (1./lphip) * phipp - ((phip*phipp)/pow(lphip,3)) * phip; - Vec<3> dx_dir_dt = Cross(dy_dir_dt,z_dir[seg]); for(int i = 0; i < 3; i++) @@ -383,22 +385,26 @@ namespace netgen dx_dir_dX[i] = grad_t(i) * dx_dir_dt; - Vec<3> grad_xbar; - for(int i=0; i<3; i++) grad_xbar(i) = -1.*(phi_minus_point * dx_dir_dX[i]) + x_dir[seg](i) - x_dir[seg] * dphi_dX[i]; double zy = z_dir[seg]*y_dir[seg]; - Vec<3> grad_ybar; Vec<3> aux = z_dir[seg] - zy*y_dir[seg]; for(int i=0; i<3; i++) grad_ybar(i) = ( (z_dir[seg]*dy_dir_dX[i])*y_dir[seg] + zy*dy_dir_dX[i] ) * phi_minus_point + aux[i] - aux * dphi_dX[i]; + */ + + // new version by JS + Vec<3> hex, hey, hez, dex, dey, dez; + CalcLocalCoordinatesDeriv (seg, t_path, hex, hey, hez, dex, dey, dez); + + grad_xbar = hex - (phi_minus_point*dex + hex*phip) * grad_t; + grad_ybar = hez - (phi_minus_point*dez + hez*phip) * grad_t; - double dFdxbar = 2.*profile_spline_coeff(0)*p2d(0) + profile_spline_coeff(2)*p2d(1) + profile_spline_coeff(3); @@ -409,31 +415,8 @@ namespace netgen grad = dFdxbar * grad_xbar + dFdybar * grad_ybar; + /* - { - cout << "grad_xbar = " << grad_xbar << endl; - cout << "grad_ybar = " << grad_ybar << endl; - Vec<3> numgradx, numgrady; - - for (int i = 0; i < 3; i++) - { - Point<3> hpl = point, hpr = point; - hpl(i) -= 1e-6; hpr(i) += 1e-6; - - Point<2> p2dl, p2dr; - CalcProj (hpl, p2dl, seg, t_path); - CalcProj (hpr, p2dr, seg, t_path); - - - numgradx(i) = (p2dr(0)-p2dl(0)) / (2e-6); - numgrady(i) = (p2dr(1)-p2dl(1)) / (2e-6); - } - cout << "num grad_xbar" << numgradx << endl; - cout << "num grad_ybar" << numgrady << endl; - } - - - { cout << "grad = " << grad << " =?= "; Vec<3> numgrad; @@ -449,24 +432,6 @@ namespace netgen } cout << " numgrad = " << numgrad << endl; - - for (int i = 0; i < 3; i++) - { - Point<3> hpl = point; - Point<3> hpr = point; - hpl(i) -= 1e-4; - hpr(i) += 1e-4; - double vall = CalcFunctionValue (hpl); - double valr = CalcFunctionValue (hpr); - numgrad(i) = (valr - vall) / (2e-4); - } - - cout << " numgrad2 = " << numgrad << endl; - } - static int cnt = 0; - cnt++; - if (cnt == 10000) - exit (0); */ } @@ -778,6 +743,52 @@ namespace netgen } + void ExtrusionFace :: + CalcLocalCoordinates (int seg, double t, + Vec<3> & ex, Vec<3> & ey, Vec<3> & ez) const + { + ey = path->GetSpline(seg).GetTangent(t); + ey /= ey.Length(); + ex = Cross (ey, glob_z_direction); + ex /= ex.Length(); + ez = Cross (ex, ey); + } + + void ExtrusionFace :: + CalcLocalCoordinatesDeriv (int seg, double t, + Vec<3> & ex, Vec<3> & ey, Vec<3> & ez, + Vec<3> & dex, Vec<3> & dey, Vec<3> & dez) const + { + Point<3> point; + Vec<3> first, second; + path->GetSpline(seg).GetDerivatives (t, point, first, second); + + ey = first; + ex = Cross (ey, glob_z_direction); + ez = Cross (ex, ey); + + dey = second; + dex = Cross (dey, glob_z_direction); + dez = Cross (dex, ey) + Cross (ex, dey); + + double lenx = ex.Length(); + double leny = ey.Length(); + double lenz = ez.Length(); + + ex /= lenx; + ey /= leny; + ez /= lenz; + + dex /= lenx; + dex -= (dex * ex) * ex; + + dey /= leny; + dey -= (dey * ey) * ey; + + dez /= lenz; + dez -= (dez * ez) * ez; + } + Extrusion :: Extrusion(const SplineGeometry<3> & path_in, const SplineGeometry<2> & profile_in, diff --git a/libsrc/csg/extrusion.hpp b/libsrc/csg/extrusion.hpp index 7a2d9355..eb2d0a6a 100644 --- a/libsrc/csg/extrusion.hpp +++ b/libsrc/csg/extrusion.hpp @@ -92,6 +92,14 @@ public: double GetProfilePar(void) const {return profile_par;} void GetRawData(Array & data) const; + + void CalcLocalCoordinates (int seg, double t, + Vec<3> & ex, Vec<3> & ey, Vec<3> & ez) const; + + void CalcLocalCoordinatesDeriv (int seg, double t, + Vec<3> & ex, Vec<3> & ey, Vec<3> & ez, + Vec<3> & dex, Vec<3> & dey, Vec<3> & dez) const; + }; diff --git a/libsrc/csg/identify.cpp b/libsrc/csg/identify.cpp index 27478160..d43c423c 100644 --- a/libsrc/csg/identify.cpp +++ b/libsrc/csg/identify.cpp @@ -256,7 +256,7 @@ Identifyable (const Point<3> & p1, const Point<3> & p2) const int PeriodicIdentification :: GetIdentifiedPoint (class Mesh & mesh, int pi) { - const Surface * sold, *snew; + const Surface *snew; const Point<3> & p = mesh.Point (pi); if (s1->PointOnSurface (p)) @@ -280,9 +280,8 @@ GetIdentifiedPoint (class Mesh & mesh, int pi) Point<3> hp = p; snew->Project (hp); - int i; int newpi = 0; - for (i = 1; i <= mesh.GetNP(); i++) + for (int i = 1; i <= mesh.GetNP(); i++) if (Dist2 (mesh.Point(i), hp) < 1e-12) { newpi = i; @@ -442,7 +441,7 @@ BuildSurfaceElements (Array & segs, Mesh & mesh, const Surface * surf) { int found = 0; - int fother; + int fother = -1; int facei = segs.Get(1).si; int surfnr = mesh.GetFaceDescriptor(facei).SurfNr(); @@ -876,7 +875,7 @@ ShortEdge (const SpecialPoint & sp1, const SpecialPoint & sp2) const int CloseSurfaceIdentification :: GetIdentifiedPoint (class Mesh & mesh, int pi) { - const Surface * sold, *snew; + const Surface *snew; const Point<3> & p = mesh.Point (pi); Array identmap(mesh.GetNP()); @@ -1062,7 +1061,7 @@ void CloseSurfaceIdentification :: IdentifyPoints (Mesh & mesh) void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh) { int fi1, fi2, side; - int s1rep, s2rep; + int s1rep = -1, s2rep = -1; for (int i = 0; i < geom.GetNSurf(); i++) { @@ -1305,9 +1304,9 @@ BuildSurfaceElements2 (Array & segs, if (!segs.Size()) return; bool found = 0; + int fother = -1; - int fother; - int facei = segs.Get(1).si; + int facei = segs[0].si; int surfnr = mesh.GetFaceDescriptor(facei).SurfNr(); @@ -1356,11 +1355,7 @@ BuildSurfaceElements2 (Array & segs, Element2d newel(sel.GetType()); newel.SetIndex (facei); for (int k = 1; k <= sel.GetNP(); k++) - { - newel.PNum(k) = - GetIdentifiedPoint (mesh, sel.PNum(k)); - // cout << "id-point = " << sel.PNum(k) << ", np = " << newel.PNum(k) << endl; - } + newel.PNum(k) = GetIdentifiedPoint (mesh, sel.PNum(k)); Vec<3> nt = Cross (Point<3> (mesh.Point (newel.PNum(2)))- Point<3> (mesh.Point (newel.PNum(1))), @@ -1572,12 +1567,9 @@ Identifyable (const SpecialPoint & sp1, const SpecialPoint & sp2, void CloseEdgesIdentification :: IdentifyPoints (Mesh & mesh) { - int i, j; - int i1, i2; - int np = mesh.GetNP(); - for (i1 = 1; i1 <= np; i1++) - for (i2 = 1; i2 <= np; i2++) + for (int i1 = 1; i1 <= np; i1++) + for (int i2 = 1; i2 <= np; i2++) { if (i2 == i1) continue; @@ -1617,15 +1609,13 @@ void CloseEdgesIdentification :: BuildSurfaceElements (Array & segs, Mesh & mesh, const Surface * surf) { - int i1, i2; int found = 0; - int i, j, k; if (surf != facet) return; - for (i1 = 1; i1 <= segs.Size(); i1++) - for (i2 = 1; i2 < i1; i2++) + for (int i1 = 1; i1 <= segs.Size(); i1++) + for (int i2 = 1; i2 < i1; i2++) { const Segment & s1 = segs.Get(i1); const Segment & s2 = segs.Get(i2); diff --git a/libsrc/csg/revolution.cpp b/libsrc/csg/revolution.cpp index b81198d5..7b1cbb36 100644 --- a/libsrc/csg/revolution.cpp +++ b/libsrc/csg/revolution.cpp @@ -45,7 +45,7 @@ namespace netgen bool first, bool last, const int id_in) : - spline(&spline_in), p0(p), v_axis(vec), isfirst(first), islast(last), id(id_in) + isfirst(first), islast(last), spline(&spline_in), p0(p), v_axis(vec), id(id_in) { deletable = false; Init(); diff --git a/libsrc/csg/singularref.cpp b/libsrc/csg/singularref.cpp index 9a37a146..0ad13db9 100644 --- a/libsrc/csg/singularref.cpp +++ b/libsrc/csg/singularref.cpp @@ -9,11 +9,11 @@ namespace netgen { SingularEdge :: SingularEdge (double abeta, int adomnr, - const CSGeometry & ageom, - const Solid * asol1, - const Solid * asol2, double sf, - const double maxh_at_initialization) - : geom(ageom), domnr(adomnr) + const CSGeometry & ageom, + const Solid * asol1, + const Solid * asol2, double sf, + const double maxh_at_initialization) + : domnr(adomnr), geom(ageom) { beta = abeta; maxhinit = maxh_at_initialization; @@ -86,8 +86,8 @@ void SingularEdge :: FindPointsOnEdge (class Mesh & mesh) int surfi1 = geom.GetSurfaceClassRepresentant(mesh[si].surfnr1); int surfi2 = geom.GetSurfaceClassRepresentant(mesh[si].surfnr2); - if (si1.Contains(surfi1) && si2.Contains(surfi2) || - si1.Contains(surfi2) && si2.Contains(surfi1)) + if ( (si1.Contains(surfi1) && si2.Contains(surfi2)) || + (si1.Contains(surfi2) && si2.Contains(surfi1)) ) // if (onedge) { diff --git a/libsrc/csg/specpoin.cpp b/libsrc/csg/specpoin.cpp index afd6176b..7f80c724 100644 --- a/libsrc/csg/specpoin.cpp +++ b/libsrc/csg/specpoin.cpp @@ -1368,7 +1368,7 @@ namespace netgen Array surfind, rep_surfind, surfind2, rep_surfind2, surfind3; Array > normalvecs; - Vec<3> nsurf; + Vec<3> nsurf = 0.0; Array specpoint2point; specpoints.SetSize (0); diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 63679c65..8fed2a10 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -80,14 +80,11 @@ namespace netgen mesh->SetNBCNames(maxsegmentindex); for ( int sindex = 0; sindex < maxsegmentindex; sindex++ ) - { - mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) ); - } + mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) ); for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++) - { - (*mesh)[si].SetBCName ( (*mesh).GetBCNamePtr( (*mesh)[si].si-1 ) ); - } + (*mesh)[si].SetBCName ( (*mesh).GetBCNamePtr( (*mesh)[si].si-1 ) ); + Point3d pmin(bbox.PMin()(0), bbox.PMin()(1), -bbox.Diam()); Point3d pmax(bbox.PMax()(0), bbox.PMax()(1), bbox.Diam()); diff --git a/libsrc/geom2d/splinegeometry.cpp b/libsrc/geom2d/splinegeometry.cpp index a9052956..7d3ab5e6 100644 --- a/libsrc/geom2d/splinegeometry.cpp +++ b/libsrc/geom2d/splinegeometry.cpp @@ -1226,12 +1226,10 @@ void SplineGeometry :: AppendDiscretePointsSegment (const Array< Point > & template string SplineGeometry :: GetBCName( const int bcnr ) const { - string bcname; if ( bcnames.Size() >= bcnr) - if ( bcnames[bcnr-1] ) - bcname = *bcnames[bcnr-1]; - else bcname = "default"; - return bcname; + if (bcnames[bcnr-1] ) + return *bcnames[bcnr-1]; + return "default"; } template diff --git a/libsrc/gprim/geomobjects.hpp b/libsrc/gprim/geomobjects.hpp index 87c38d9b..85059c13 100644 --- a/libsrc/gprim/geomobjects.hpp +++ b/libsrc/gprim/geomobjects.hpp @@ -22,7 +22,7 @@ protected: public: Point () { ; } - Point (double ax) { x[0] = ax; } + Point (double ax) { for (int i = 0; i < D; i++) x[i] = ax; } Point (double ax, double ay) { x[0] = ax; x[1] = ay; } Point (double ax, double ay, double az) { x[0] = ax; x[1] = ay; x[2] = az; } diff --git a/libsrc/gprim/geomtest3d.cpp b/libsrc/gprim/geomtest3d.cpp index 20a3be0b..97501d7e 100644 --- a/libsrc/gprim/geomtest3d.cpp +++ b/libsrc/gprim/geomtest3d.cpp @@ -429,7 +429,7 @@ int IntersectTetTriangleRef (const Point<3> ** tri, const int * tripi) { int i, j; double eps = 1e-8; - double eps2 = eps * eps; + // double eps2 = eps * eps; static Point<3> rtetp1(0, 0, 0); static Point<3> rtetp2(1, 0, 0); diff --git a/libsrc/interface/readtetmesh.cpp b/libsrc/interface/readtetmesh.cpp index 524b0bcc..1258068b 100644 --- a/libsrc/interface/readtetmesh.cpp +++ b/libsrc/interface/readtetmesh.cpp @@ -50,7 +50,7 @@ namespace netgen int numObj3D,numObj2D,numObj1D,numObj0D; bool nullstarted; Array eldom; - int minId3D,minId2D; + int minId3D = -1, minId2D = -1; int maxId3D(-1), maxId2D(-1), maxId1D(-1), maxId0D(-1); Array *> segmentdata; Array tris; diff --git a/libsrc/interface/readuser.cpp b/libsrc/interface/readuser.cpp index 847ffb26..98f325d9 100644 --- a/libsrc/interface/readuser.cpp +++ b/libsrc/interface/readuser.cpp @@ -22,8 +22,6 @@ namespace netgen const char * filename = hfilename.c_str(); - int i, j; - char reco[100]; int np, nbe; @@ -41,7 +39,7 @@ namespace netgen in >> reco; in >> np; - for (i = 1; i <= np; i++) + for (int i = 1; i <= np; i++) { Point3d p; in >> p.X() >> p.Y() >> p.Z(); @@ -53,15 +51,12 @@ namespace netgen in >> nbe; // int invert = globflags.GetDefineFlag ("invertsurfacemesh"); - for (i = 1; i <= nbe; i++) + for (int i = 1; i <= nbe; i++) { Element2d el; - int hi; - el.SetIndex(1); - - // in >> hi; - for (j = 1; j <= 3; j++) + + for (int j = 1; j <= 3; j++) { in >> el.PNum(j); // el.PNum(j)++; @@ -93,11 +88,11 @@ namespace netgen if ( (strlen (filename) > 4) && strcmp (&filename[strlen (filename)-4], ".unv") == 0 ) { - int i, j, k; + // int i, j, k; - double h; + // double h; char reco[100]; - int np, nbe; + // int np, nbe; int invert; @@ -116,13 +111,13 @@ namespace netgen if (strcmp (reco, "NODES") == 0) { cout << "nodes found" << endl; - for (j = 1; j <= 4; j++) + for (int j = 1; j <= 4; j++) in >> reco; // read dummy while (1) { int pi, hi; - double x, y, z; + // double x, y, z; Point3d p; in >> pi; @@ -144,7 +139,7 @@ namespace netgen if (strcmp (reco, "ELEMENTS") == 0) { cout << "elements found" << endl; - for (j = 1; j <= 4; j++) + for (int j = 1; j <= 4; j++) in >> reco; // read dummy while (1) @@ -152,7 +147,7 @@ namespace netgen int hi; in >> hi; if (hi == -1) break; - for (j = 1; j <= 7; j++) + for (int j = 1; j <= 7; j++) in >> hi; Element2d el; @@ -163,7 +158,7 @@ namespace netgen swap (el.PNum(2), el.PNum(3)); mesh.AddSurfaceElement (el); - for (j = 1; j <= 5; j++) + for (int j = 1; j <= 5; j++) in >> hi; } } @@ -260,7 +255,7 @@ namespace netgen in >> nse; for (i = 1; i <= nse; i++) { - int mat, nelp; + int mat; // , nelp; in >> mat; Element2d el (TRIG); el.SetIndex (mat); @@ -301,7 +296,7 @@ namespace netgen cout << "pktfile = " << pktfile << endl; int np, nse, i; - int num, bcprop; + int bcprop; ifstream inpkt (pktfile.c_str()); inpkt >> np; Array values(np); diff --git a/libsrc/interface/writedolfin.cpp b/libsrc/interface/writedolfin.cpp index 0fe36e94..8fd9e749 100644 --- a/libsrc/interface/writedolfin.cpp +++ b/libsrc/interface/writedolfin.cpp @@ -25,14 +25,14 @@ namespace netgen int np = mesh.GetNP(); int ne = mesh.GetNE(); - int nse = mesh.GetNSE(); + // int nse = mesh.GetNSE(); int nsd = mesh.GetDimension(); - int invertsurf = mparam.inverttrigs; - int i, j; + // int invertsurf = mparam.inverttrigs; + // int i, j; ofstream outfile (filename.c_str()); - char str[100]; + // char str[100]; outfile.precision(8); outfile.setf (ios::fixed, ios::floatfield); outfile.setf (ios::showpoint); @@ -45,7 +45,7 @@ namespace netgen outfile << ""<" <"<"<"<"< points; -static Array volelements; -static Array surfelements; + static Array points; + static Array volelements; + static Array surfelements; -static Array faces; -static Array edges; + static Array faces; + static Array edges; -void ReadFile (char * filename) + void ReadFile (char * filename) { - int i, n; - ifstream infile(filename); - char reco[100]; + int i, n; + ifstream infile(filename); + char reco[100]; - infile >> reco; // file format recognition + infile >> reco; // file format recognition - infile >> n; // number of surface elements - cout << n << " Surface elements" << endl; + infile >> n; // number of surface elements + cout << n << " Surface elements" << endl; - for (i = 1; i <= n; i++) - { - SURFELEMENT sel; - infile >> sel.snr >> sel.p1 >> sel.p2 >> sel.p3; - surfelements.Append (sel); - } - - infile >> n; // number of volume elements - cout << n << " Volume elements" << endl; - - for (i = 1; i <= n; i++) - { - VOLELEMENT el; - infile >> el.p1 >> el.p2 >> el.p3 >> el.p4; - volelements.Append (el); - } - - infile >> n; // number of points - cout << n << " Points" << endl; - - for (i = 1; i <= n; i++) - { - POINT3D p; - infile >> p.x >> p.y >> p.z; - points.Append (p); - } - } - - - -void ReadFileMesh (const Mesh & mesh) -{ - int i, n; - - n = mesh.GetNSE(); // number of surface elements - cout << n << " Surface elements" << endl; - - for (i = 1; i <= n; i++) - { - SURFELEMENT sel; - const Element2d & el = mesh.SurfaceElement(i); - sel.snr = el.GetIndex(); - sel.p1 = el.PNum(1); - sel.p2 = el.PNum(2); - sel.p3 = el.PNum(3); - surfelements.Append (sel); - } - - n = mesh.GetNE(); // number of volume elements - cout << n << " Volume elements" << endl; - - for (i = 1; i <= n; i++) - { - VOLELEMENT el; - const Element & nel = mesh.VolumeElement(i); - el.p1 = nel.PNum(1); - el.p2 = nel.PNum(2); - el.p3 = nel.PNum(3); - el.p4 = nel.PNum(4); - // infile >> el.p1 >> el.p2 >> el.p3 >> el.p4; - volelements.Append (el); - } - - n = mesh.GetNP(); // number of points - cout << n << " Points" << endl; - - for (i = 1; i <= n; i++) - { - POINT3D p; - Point3d mp = mesh.Point(i); - p.x = mp.X(); - p.y = mp.Y(); - p.z = mp.Z(); - // infile >> p.x >> p.y >> p.z; - points.Append (p); - } - } - - - - -void Convert () - { - int i, j, facei, edgei; - INDEX_3 i3; - INDEX_2 i2; - - INDEX_3_HASHTABLE faceindex(volelements.Size()/5 + 1); - INDEX_2_HASHTABLE edgeindex(volelements.Size()/5 + 1); - - for (i = 1; i <= volelements.Size(); i++) - { - for (j = 1; j <= 4; j++) + for (i = 1; i <= n; i++) { - switch (j) - { - case 1: - i3.I1() = volelements.Get(i).p2; - i3.I2() = volelements.Get(i).p3; - i3.I3() = volelements.Get(i).p4; - break; - case 2: - i3.I1() = volelements.Get(i).p1; - i3.I2() = volelements.Get(i).p3; - i3.I3() = volelements.Get(i).p4; - break; - case 3: - i3.I1() = volelements.Get(i).p1; - i3.I2() = volelements.Get(i).p2; - i3.I3() = volelements.Get(i).p4; - break; - case 4: - i3.I1() = volelements.Get(i).p1; - i3.I2() = volelements.Get(i).p2; - i3.I3() = volelements.Get(i).p3; - break; - default: - i3.I1()=i3.I2()=i3.I3()=0; - } - i3.Sort(); - if (faceindex.Used (i3)) - facei = faceindex.Get(i3); - else - { - FACE fa; - fa.p1 = i3.I1(); - fa.p2 = i3.I2(); - fa.p3 = i3.I3(); - facei = faces.Append (fa); - faceindex.Set (i3, facei); - } - - volelements.Elem(i).faces[j-1] = facei; - } + SURFELEMENT sel; + infile >> sel.snr >> sel.p1 >> sel.p2 >> sel.p3; + surfelements.Append (sel); + } - } + infile >> n; // number of volume elements + cout << n << " Volume elements" << endl; + + for (i = 1; i <= n; i++) + { + VOLELEMENT el; + infile >> el.p1 >> el.p2 >> el.p3 >> el.p4; + volelements.Append (el); + } + + infile >> n; // number of points + cout << n << " Points" << endl; + + for (i = 1; i <= n; i++) + { + POINT3D p; + infile >> p.x >> p.y >> p.z; + points.Append (p); + } + } + + + + void ReadFileMesh (const Mesh & mesh) + { + int i, n; + + n = mesh.GetNSE(); // number of surface elements + cout << n << " Surface elements" << endl; + + for (i = 1; i <= n; i++) + { + SURFELEMENT sel; + const Element2d & el = mesh.SurfaceElement(i); + sel.snr = el.GetIndex(); + sel.p1 = el.PNum(1); + sel.p2 = el.PNum(2); + sel.p3 = el.PNum(3); + surfelements.Append (sel); + } + + n = mesh.GetNE(); // number of volume elements + cout << n << " Volume elements" << endl; + + for (i = 1; i <= n; i++) + { + VOLELEMENT el; + const Element & nel = mesh.VolumeElement(i); + el.p1 = nel.PNum(1); + el.p2 = nel.PNum(2); + el.p3 = nel.PNum(3); + el.p4 = nel.PNum(4); + // infile >> el.p1 >> el.p2 >> el.p3 >> el.p4; + volelements.Append (el); + } + + n = mesh.GetNP(); // number of points + cout << n << " Points" << endl; + + for (i = 1; i <= n; i++) + { + POINT3D p; + Point3d mp = mesh.Point(i); + p.x = mp.X(); + p.y = mp.Y(); + p.z = mp.Z(); + // infile >> p.x >> p.y >> p.z; + points.Append (p); + } + } + + + + + void Convert () + { + int i, j, facei, edgei; + INDEX_3 i3; + INDEX_2 i2; + + INDEX_3_HASHTABLE faceindex(volelements.Size()/5 + 1); + INDEX_2_HASHTABLE edgeindex(volelements.Size()/5 + 1); + + for (i = 1; i <= volelements.Size(); i++) + { + for (j = 1; j <= 4; j++) + { + switch (j) + { + case 1: + i3.I1() = volelements.Get(i).p2; + i3.I2() = volelements.Get(i).p3; + i3.I3() = volelements.Get(i).p4; + break; + case 2: + i3.I1() = volelements.Get(i).p1; + i3.I2() = volelements.Get(i).p3; + i3.I3() = volelements.Get(i).p4; + break; + case 3: + i3.I1() = volelements.Get(i).p1; + i3.I2() = volelements.Get(i).p2; + i3.I3() = volelements.Get(i).p4; + break; + case 4: + i3.I1() = volelements.Get(i).p1; + i3.I2() = volelements.Get(i).p2; + i3.I3() = volelements.Get(i).p3; + break; + default: + i3.I1()=i3.I2()=i3.I3()=0; + } + i3.Sort(); + if (faceindex.Used (i3)) + facei = faceindex.Get(i3); + else + { + FACE fa; + fa.p1 = i3.I1(); + fa.p2 = i3.I2(); + fa.p3 = i3.I3(); + facei = faces.Append (fa); + faceindex.Set (i3, facei); + } + + volelements.Elem(i).faces[j-1] = facei; + } + + } - for (i = 1; i <= faces.Size(); i++) - { - for (j = 1; j <= 3; j++) + for (i = 1; i <= faces.Size(); i++) { - switch (j) - { - case 1: - i2.I1() = faces.Get(i).p2; - i2.I2() = faces.Get(i).p3; - break; - case 2: - i2.I1() = faces.Get(i).p1; - i2.I2() = faces.Get(i).p3; - break; - case 3: - i2.I1() = faces.Get(i).p1; - i2.I2() = faces.Get(i).p2; - break; - default: - i2.I1()=i2.I2()=0; - } - if (i2.I1() > i2.I2()) swap (i2.I1(), i2.I2()); - if (edgeindex.Used (i2)) - edgei = edgeindex.Get(i2); - else - { - EDGE ed; - ed.p1 = i2.I1(); - ed.p2 = i2.I2(); - edgei = edges.Append (ed); - edgeindex.Set (i2, edgei); - } + for (j = 1; j <= 3; j++) + { + switch (j) + { + case 1: + i2.I1() = faces.Get(i).p2; + i2.I2() = faces.Get(i).p3; + break; + case 2: + i2.I1() = faces.Get(i).p1; + i2.I2() = faces.Get(i).p3; + break; + case 3: + i2.I1() = faces.Get(i).p1; + i2.I2() = faces.Get(i).p2; + break; + default: + i2.I1()=i2.I2()=0; + } + if (i2.I1() > i2.I2()) swap (i2.I1(), i2.I2()); + if (edgeindex.Used (i2)) + edgei = edgeindex.Get(i2); + else + { + EDGE ed; + ed.p1 = i2.I1(); + ed.p2 = i2.I2(); + edgei = edges.Append (ed); + edgeindex.Set (i2, edgei); + } - faces.Elem(i).edges[j-1] = edgei; - } + faces.Elem(i).edges[j-1] = edgei; + } - } + } } -void WriteFile (ostream & outfile) + void WriteFile (ostream & outfile) { - int i; + int i; - outfile - << "#VERSION: 1.0" << endl - << "#PROGRAM: NETGEN" << endl - << "#EQN_TYPE: POISSON" << endl - << "#DIMENSION: 3D" << endl - << "#DEG_OF_FREE: 1" << endl - << "#DESCRIPTION: I don't know" << endl - << "##RENUM: not done" << endl - << "#USER: Kleinzen" << endl - << "DATE: 10.06.1996" << endl; + outfile + << "#VERSION: 1.0" << endl + << "#PROGRAM: NETGEN" << endl + << "#EQN_TYPE: POISSON" << endl + << "#DIMENSION: 3D" << endl + << "#DEG_OF_FREE: 1" << endl + << "#DESCRIPTION: I don't know" << endl + << "##RENUM: not done" << endl + << "#USER: Kleinzen" << endl + << "DATE: 10.06.1996" << endl; - outfile << "#HEADER: 8" << endl - << points.Size() << " " << edges.Size() << " " - << faces.Size() << " " << volelements.Size() << " 0 0 0 0" << endl; + outfile << "#HEADER: 8" << endl + << points.Size() << " " << edges.Size() << " " + << faces.Size() << " " << volelements.Size() << " 0 0 0 0" << endl; - outfile << "#VERTEX: " << points.Size() << endl; - for (i = 1; i <= points.Size(); i++) - outfile << " " << i << " " << points.Get(i).x << " " << points.Get(i).y - << " " << points.Get(i).z << endl; + outfile << "#VERTEX: " << points.Size() << endl; + for (i = 1; i <= points.Size(); i++) + outfile << " " << i << " " << points.Get(i).x << " " << points.Get(i).y + << " " << points.Get(i).z << endl; - outfile << "#EDGE: " << edges.Size() << endl; - for (i = 1; i <= edges.Size(); i++) - outfile << " " << i << " 1 " - << edges.Get(i).p1 << " " - << edges.Get(i).p2 - << " 0" << endl; + outfile << "#EDGE: " << edges.Size() << endl; + for (i = 1; i <= edges.Size(); i++) + outfile << " " << i << " 1 " + << edges.Get(i).p1 << " " + << edges.Get(i).p2 + << " 0" << endl; - outfile << "#FACE: " << faces.Size() << endl; - for (i = 1; i <= faces.Size(); i++) - outfile << " " << i << " 1 3 " - << faces.Get(i).edges[0] << " " - << faces.Get(i).edges[1] << " " - << faces.Get(i).edges[2] << endl; + outfile << "#FACE: " << faces.Size() << endl; + for (i = 1; i <= faces.Size(); i++) + outfile << " " << i << " 1 3 " + << faces.Get(i).edges[0] << " " + << faces.Get(i).edges[1] << " " + << faces.Get(i).edges[2] << endl; - outfile << "#SOLID: " << volelements.Size() << endl; - for (i = 1; i <= volelements.Size(); i++) - outfile << " " << i << " 1 4 " - << volelements.Get(i).faces[0] << " " - << volelements.Get(i).faces[1] << " " - << volelements.Get(i).faces[2] << " " - << volelements.Get(i).faces[3] << endl; + outfile << "#SOLID: " << volelements.Size() << endl; + for (i = 1; i <= volelements.Size(); i++) + outfile << " " << i << " 1 4 " + << volelements.Get(i).faces[0] << " " + << volelements.Get(i).faces[1] << " " + << volelements.Get(i).faces[2] << " " + << volelements.Get(i).faces[3] << endl; - outfile << "#END_OF_DATA" << endl; + outfile << "#END_OF_DATA" << endl; } -void WriteUserChemnitz (const Mesh & mesh, - const string & filename) -{ - ofstream outfile (filename.c_str()); + void WriteUserChemnitz (const Mesh & mesh, + const string & filename) + { + ofstream outfile (filename.c_str()); - ReadFileMesh (mesh); - Convert (); + ReadFileMesh (mesh); + Convert (); - WriteFile (outfile); - cout << "Wrote Chemnitz standard file" << endl; -} + WriteFile (outfile); + cout << "Wrote Chemnitz standard file" << endl; + } } diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 712a970d..163996eb 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -362,6 +362,8 @@ namespace netgen ned = 6; break; } + default: + throw NgException("Bisect, element type not handled in switch"); } for (j = 0; j < ned; j++) @@ -462,7 +464,7 @@ namespace netgen { { 1, 2, 4, 3 }, { 1, 5, 4, 5 }, { 2, 5, 3, 5 } }; - + int (*pairs)[4] = NULL; switch (el.GetType()) { @@ -477,6 +479,8 @@ namespace netgen pairs = pyramidpairs; break; } + default: + throw NgException("Bisect, element type not handled in switch, 2"); } for (j = 0; j < 3; j++) @@ -744,6 +748,8 @@ namespace netgen ned = 6; break; } + default: + throw NgException("Bisect, element type not handled in switch, 3"); } for (j = 0; j < ned; j++) @@ -802,6 +808,8 @@ namespace netgen pairs = pyramidpairs; break; } + default: + throw NgException("Bisect, element type not handled in switch, 3a"); } for (j = 0; j < 3; j++) @@ -2100,6 +2108,8 @@ namespace netgen mprisms.Append (mp); break; } + default: + throw NgException("Bisect, element type not handled in switch, 4"); } } @@ -2523,6 +2533,8 @@ namespace netgen << endl; break; } + default: + throw NgException("Bisect, element type not handled in switch, 6"); } } @@ -2563,6 +2575,8 @@ namespace netgen mquads.Append (mt); break; } + default: + throw NgException("Bisect, element type not handled in switch, 5"); } diff --git a/libsrc/meshing/classifyhpel.hpp b/libsrc/meshing/classifyhpel.hpp index 7808d199..254ef19c 100644 --- a/libsrc/meshing/classifyhpel.hpp +++ b/libsrc/meshing/classifyhpel.hpp @@ -792,6 +792,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge if(cornerpoint.Test(el.PNum(p[k]))) point_sing[p[k]-1] = 3; + *testout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl; if(edge_sing[0] + edge_sing[1] + edge_sing[2] == 0) { @@ -851,8 +852,11 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge // cout << " run for " << j << " gives type " << type << endl; //*testout << " run for " << j << " gives type " << type << endl; + if(type!=HP_NONE) break; } + + *testout << "type = " << type << endl; for(int k=0;k<3;k++) el[k] = pnums[k]; /*if(type != HP_NONE) diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 4e826dd6..c7aab2fb 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -367,7 +367,6 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, AutoDiff<3> adt(t, 2); AutoDiff<3> res[2000]; CalcScaledTrigShape (n, adx, ady, adt, &res[0]); - double dshape1[6000]; int ndof = (n-1)*(n-2)/2; for (int i = 0; i < ndof; i++) { @@ -377,6 +376,7 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, } /* + double dshape1[6000]; if (n < 3) return; double hvl[1000], hvr[1000]; int nd = (n-1)*(n-2)/2; @@ -1486,6 +1486,9 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, break; } + + default: + throw NgException("CurvedElements::CalcShape 2d, element type not handled"); }; } @@ -1746,6 +1749,9 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, */ break; } + default: + throw NgException("CurvedElements::CalcDShape 2d, element type not handled"); + }; } @@ -2189,6 +2195,10 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, shapes[7] = (1-x)* y *(z); break; } + + default: + throw NgException("CurvedElements::CalcShape 3d, element type not handled"); + }; } @@ -2616,6 +2626,9 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, break; } + + default: + throw NgException("CurvedElements::CalcDShape 3d, element type not handled"); } /* diff --git a/libsrc/meshing/findip.hpp b/libsrc/meshing/findip.hpp index 0ce90597..f5d8e424 100644 --- a/libsrc/meshing/findip.hpp +++ b/libsrc/meshing/findip.hpp @@ -78,7 +78,7 @@ inline int FindInnerPoint (POINTArray & points, Array a; Array c; Mat<3> m, inv; - Vec<3> rs, x, center; + Vec<3> rs, x = 0.0, center; double f; int nf = faces.Size(); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 7874a08d..e49a554c 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -5,227 +5,227 @@ namespace netgen { - Mesh :: Mesh () - { - volelements.SetName ("vol elements"); - surfelements.SetName ("surf elements"); - points.SetName ("meshpoints"); + Mesh :: Mesh () + { + volelements.SetName ("vol elements"); + surfelements.SetName ("surf elements"); + points.SetName ("meshpoints"); - boundaryedges = NULL; - surfelementht = NULL; - segmentht = NULL; + boundaryedges = NULL; + surfelementht = NULL; + segmentht = NULL; - lochfunc = NULL; - mglevels = 1; - elementsearchtree = NULL; - elementsearchtreets = NextTimeStamp(); - majortimestamp = timestamp = NextTimeStamp(); - hglob = 1e10; - hmin = 0; - numvertices = -1; - dimension = 3; + lochfunc = NULL; + mglevels = 1; + elementsearchtree = NULL; + elementsearchtreets = NextTimeStamp(); + majortimestamp = timestamp = NextTimeStamp(); + hglob = 1e10; + hmin = 0; + numvertices = -1; + dimension = 3; - topology = new MeshTopology (*this); - curvedelems = new CurvedElements (*this); - clusters = new AnisotropicClusters (*this); - ident = new Identifications (*this); + topology = new MeshTopology (*this); + curvedelems = new CurvedElements (*this); + clusters = new AnisotropicClusters (*this); + ident = new Identifications (*this); - hpelements = NULL; - coarsemesh = NULL; + hpelements = NULL; + coarsemesh = NULL; - ps_startelement = 0; + ps_startelement = 0; - geomtype = NO_GEOM; + geomtype = NO_GEOM; - bcnames.SetSize(0); + bcnames.SetSize(0); #ifdef PARALLEL - paralleltop = new ParallelMeshTopology (*this); + paralleltop = new ParallelMeshTopology (*this); #endif - } + } - Mesh :: ~Mesh() - { - delete lochfunc; - delete boundaryedges; - delete surfelementht; - delete segmentht; - delete curvedelems; - delete clusters; - delete topology; - delete ident; - delete elementsearchtree; + Mesh :: ~Mesh() + { + delete lochfunc; + delete boundaryedges; + delete surfelementht; + delete segmentht; + delete curvedelems; + delete clusters; + delete topology; + delete ident; + delete elementsearchtree; - delete coarsemesh; - delete hpelements; + delete coarsemesh; + delete hpelements; - for (int i = 0; i < materials.Size(); i++) - delete [] materials[i]; + for (int i = 0; i < materials.Size(); i++) + delete [] materials[i]; - for(int i = 0; i < userdata_int.Size(); i++) - delete userdata_int[i]; - for(int i = 0; i < userdata_double.Size(); i++) - delete userdata_double[i]; + for(int i = 0; i < userdata_int.Size(); i++) + delete userdata_int[i]; + for(int i = 0; i < userdata_double.Size(); i++) + delete userdata_double[i]; - for (int i = 0; i < bcnames.Size(); i++ ) - if ( bcnames[i] ) delete bcnames[i]; + for (int i = 0; i < bcnames.Size(); i++ ) + if ( bcnames[i] ) delete bcnames[i]; #ifdef PARALLEL - delete paralleltop; + delete paralleltop; #endif - } + } - Mesh & Mesh :: operator= (const Mesh & mesh2) - { - points = mesh2.points; - // eltyps = mesh2.eltyps; - segments = mesh2.segments; - surfelements = mesh2.surfelements; - volelements = mesh2.volelements; - lockedpoints = mesh2.lockedpoints; - facedecoding = mesh2.facedecoding; - dimension = mesh2.dimension; + Mesh & Mesh :: operator= (const Mesh & mesh2) + { + points = mesh2.points; + // eltyps = mesh2.eltyps; + segments = mesh2.segments; + surfelements = mesh2.surfelements; + volelements = mesh2.volelements; + lockedpoints = mesh2.lockedpoints; + facedecoding = mesh2.facedecoding; + dimension = mesh2.dimension; - bcnames.SetSize( mesh2.bcnames.Size() ); - for ( int i = 0; i < mesh2.bcnames.Size(); i++ ) - if ( mesh2.bcnames[i] ) bcnames[i] = new string ( *mesh2.bcnames[i] ); - else bcnames[i] = 0; + bcnames.SetSize( mesh2.bcnames.Size() ); + for ( int i = 0; i < mesh2.bcnames.Size(); i++ ) + if ( mesh2.bcnames[i] ) bcnames[i] = new string ( *mesh2.bcnames[i] ); + else bcnames[i] = 0; - return *this; - } + return *this; + } - void Mesh :: DeleteMesh() - { - NgLock lock(mutex); - lock.Lock(); + void Mesh :: DeleteMesh() + { + NgLock lock(mutex); + lock.Lock(); - points.SetSize(0); - segments.SetSize(0); - surfelements.SetSize(0); - volelements.SetSize(0); - lockedpoints.SetSize(0); - surfacesonnode.SetSize(0); + points.SetSize(0); + segments.SetSize(0); + surfelements.SetSize(0); + volelements.SetSize(0); + lockedpoints.SetSize(0); + surfacesonnode.SetSize(0); - delete boundaryedges; - boundaryedges = NULL; + delete boundaryedges; + boundaryedges = NULL; - openelements.SetSize(0); - facedecoding.SetSize(0); + openelements.SetSize(0); + facedecoding.SetSize(0); - delete ident; - ident = new Identifications (*this); - delete topology; - topology = new MeshTopology (*this); - delete curvedelems; - curvedelems = new CurvedElements (*this); - delete clusters; - clusters = new AnisotropicClusters (*this); + delete ident; + ident = new Identifications (*this); + delete topology; + topology = new MeshTopology (*this); + delete curvedelems; + curvedelems = new CurvedElements (*this); + delete clusters; + clusters = new AnisotropicClusters (*this); - for ( int i = 0; i < bcnames.Size(); i++ ) - if ( bcnames[i] ) delete bcnames[i]; + for ( int i = 0; i < bcnames.Size(); i++ ) + if ( bcnames[i] ) delete bcnames[i]; #ifdef PARALLEL - delete paralleltop; - paralleltop = new ParallelMeshTopology (*this); + delete paralleltop; + paralleltop = new ParallelMeshTopology (*this); #endif - lock.UnLock(); + lock.UnLock(); - timestamp = NextTimeStamp(); - } + timestamp = NextTimeStamp(); + } - PointIndex Mesh :: AddPoint (const Point3d & p, int layer) - { - NgLock lock(mutex); - lock.Lock(); + PointIndex Mesh :: AddPoint (const Point3d & p, int layer) + { + NgLock lock(mutex); + lock.Lock(); - timestamp = NextTimeStamp(); + timestamp = NextTimeStamp(); - PointIndex pi = points.Size() + PointIndex::BASE; - points.Append ( MeshPoint (p, layer, INNERPOINT) ); + PointIndex pi = points.Size() + PointIndex::BASE; + points.Append ( MeshPoint (p, layer, INNERPOINT) ); #ifdef PARALLEL - points.Last().SetGhost(0); + points.Last().SetGhost(0); #endif - lock.UnLock(); + lock.UnLock(); - return pi; - } + return pi; + } - PointIndex Mesh :: AddPoint (const Point3d & p, int layer, POINTTYPE type) - { - NgLock lock(mutex); - lock.Lock(); + PointIndex Mesh :: AddPoint (const Point3d & p, int layer, POINTTYPE type) + { + NgLock lock(mutex); + lock.Lock(); - timestamp = NextTimeStamp(); + timestamp = NextTimeStamp(); - PointIndex pi = points.Size() + PointIndex::BASE; - points.Append ( MeshPoint (p, layer, type) ); + PointIndex pi = points.Size() + PointIndex::BASE; + points.Append ( MeshPoint (p, layer, type) ); #ifdef PARALLEL - points.Last().SetGhost(0); + points.Last().SetGhost(0); #endif - lock.UnLock(); + lock.UnLock(); - return pi; - } + return pi; + } #ifdef PARALLEL - PointIndex Mesh :: AddPoint (const Point3d & p, bool isghost, int layer) - { - NgLock lock(mutex); - lock.Lock(); + PointIndex Mesh :: AddPoint (const Point3d & p, bool isghost, int layer) + { + NgLock lock(mutex); + lock.Lock(); - timestamp = NextTimeStamp(); + timestamp = NextTimeStamp(); - PointIndex pi = points.Size() + PointIndex::BASE; - points.Append ( MeshPoint (p, layer, INNERPOINT) ); + PointIndex pi = points.Size() + PointIndex::BASE; + points.Append ( MeshPoint (p, layer, INNERPOINT) ); - points.Last().SetGhost(isghost); + points.Last().SetGhost(isghost); - lock.UnLock(); + lock.UnLock(); - return pi; - } + return pi; + } - PointIndex Mesh :: AddPoint (const Point3d & p, bool isghost, int layer, POINTTYPE type) - { - NgLock lock(mutex); - lock.Lock(); + PointIndex Mesh :: AddPoint (const Point3d & p, bool isghost, int layer, POINTTYPE type) + { + NgLock lock(mutex); + lock.Lock(); - timestamp = NextTimeStamp(); + timestamp = NextTimeStamp(); - PointIndex pi = points.Size() + PointIndex::BASE; - points.Append ( MeshPoint (p, layer, type) ); + PointIndex pi = points.Size() + PointIndex::BASE; + points.Append ( MeshPoint (p, layer, type) ); - points.Last().SetGhost(isghost); + points.Last().SetGhost(isghost); - lock.UnLock(); + lock.UnLock(); - return pi; - } + return pi; + } #endif - SegmentIndex Mesh :: AddSegment (const Segment & s) - { - NgLock lock(mutex); - lock.Lock(); - timestamp = NextTimeStamp(); + SegmentIndex Mesh :: AddSegment (const Segment & s) + { + NgLock lock(mutex); + lock.Lock(); + timestamp = NextTimeStamp(); - int maxn = max2 (s[0], s[1]); - maxn += 1-PointIndex::BASE; + int maxn = max2 (s[0], s[1]); + maxn += 1-PointIndex::BASE; - /* + /* if (maxn > ptyps.Size()) { int maxo = ptyps.Size(); @@ -236,42 +236,42 @@ namespace netgen if (ptyps[s[0]] > EDGEPOINT) ptyps[s[0]] = EDGEPOINT; if (ptyps[s[1]] > EDGEPOINT) ptyps[s[1]] = EDGEPOINT; - */ + */ - if (maxn <= points.Size()) + if (maxn <= points.Size()) { - if (points[s[0]].Type() > EDGEPOINT) - points[s[0]].SetType (EDGEPOINT); - if (points[s[1]].Type() > EDGEPOINT) - points[s[1]].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 { cerr << "edge points nrs > points.Size" << endl; } - */ + */ - SegmentIndex si = segments.Size(); - segments.Append (s); + SegmentIndex si = segments.Size(); + segments.Append (s); - lock.UnLock(); - return si; - } + lock.UnLock(); + return si; + } - SurfaceElementIndex Mesh :: AddSurfaceElement (const Element2d & el) - { - NgLock lock(mutex); - lock.Lock(); - timestamp = NextTimeStamp(); + SurfaceElementIndex Mesh :: AddSurfaceElement (const Element2d & el) + { + NgLock lock(mutex); + lock.Lock(); + timestamp = NextTimeStamp(); - int maxn = el[0]; - for (int i = 1; i < el.GetNP(); i++) - if (el[i] > maxn) maxn = el[i]; + int maxn = el[0]; + for (int i = 1; i < el.GetNP(); i++) + if (el[i] > maxn) maxn = el[i]; - maxn += 1-PointIndex::BASE; + maxn += 1-PointIndex::BASE; - /* + /* if (maxn > ptyps.Size()) { int maxo = ptyps.Size(); @@ -281,50 +281,50 @@ namespace netgen ptyps[i] = INNERPOINT; } - */ - if (maxn <= points.Size()) + */ + if (maxn <= points.Size()) { - for (int i = 0; i < el.GetNP(); i++) - if (points[el[i]].Type() > SURFACEPOINT) - points[el[i]].SetType(SURFACEPOINT); + for (int i = 0; i < el.GetNP(); i++) + if (points[el[i]].Type() > SURFACEPOINT) + points[el[i]].SetType(SURFACEPOINT); } - /* + /* else { cerr << "surf points nrs > points.Size" << endl; } - */ + */ - SurfaceElementIndex si = surfelements.Size(); - surfelements.Append (el); + SurfaceElementIndex si = surfelements.Size(); + surfelements.Append (el); - if (el.index > facedecoding.Size()) - cerr << "has no facedecoding: fd.size = " << facedecoding.Size() << ", ind = " << el.index << endl; + if (el.index > facedecoding.Size()) + cerr << "has no facedecoding: fd.size = " << facedecoding.Size() << ", ind = " << el.index << endl; - surfelements.Last().next = facedecoding[el.index-1].firstelement; - facedecoding[el.index-1].firstelement = si; + surfelements.Last().next = facedecoding[el.index-1].firstelement; + facedecoding[el.index-1].firstelement = si; #ifdef PARALLEL - surfelements.Last().SetGhost ( el.IsGhost() ); + surfelements.Last().SetGhost ( el.IsGhost() ); #endif - lock.UnLock(); - return si; - } + lock.UnLock(); + return si; + } - ElementIndex Mesh :: AddVolumeElement (const Element & el) - { - NgLock lock(mutex); - lock.Lock(); + ElementIndex Mesh :: AddVolumeElement (const Element & el) + { + NgLock lock(mutex); + lock.Lock(); - int maxn = el[0]; - for (int i = 1; i < el.GetNP(); i++) - if (el[i] > maxn) maxn = el[i]; + int maxn = el[0]; + for (int i = 1; i < el.GetNP(); i++) + if (el[i] > maxn) maxn = el[i]; - maxn += 1-PointIndex::BASE; + maxn += 1-PointIndex::BASE; - /* + /* if (maxn > ptyps.Size()) { int maxo = ptyps.Size(); @@ -333,88 +333,88 @@ namespace netgen i < maxn+PointIndex::BASE; i++) ptyps[i] = INNERPOINT; } - */ - /* + */ + /* if (maxn > points.Size()) { cerr << "add vol element before point" << endl; } - */ + */ - int ve = volelements.Size(); + int ve = volelements.Size(); - volelements.Append (el); - volelements.Last().flags.illegal_valid = 0; + volelements.Append (el); + volelements.Last().flags.illegal_valid = 0; #ifdef PARALLEL - volelements.Last().SetGhost ( el.IsGhost() ); + volelements.Last().SetGhost ( el.IsGhost() ); #endif - // while (volelements.Size() > eltyps.Size()) - // eltyps.Append (FREEELEMENT); + // while (volelements.Size() > eltyps.Size()) + // eltyps.Append (FREEELEMENT); - timestamp = NextTimeStamp(); + timestamp = NextTimeStamp(); - lock.UnLock(); - return ve; - } + lock.UnLock(); + return ve; + } - void Mesh :: Save (const string & filename) const - { + void Mesh :: Save (const string & filename) const + { - ofstream outfile(filename.c_str()); + ofstream outfile(filename.c_str()); - Save(outfile); - } + Save(outfile); + } - void Mesh :: Save (ostream & outfile) const - { - int i, j; + void Mesh :: Save (ostream & outfile) const + { + int i, j; - double scale = 1; // globflags.GetNumFlag ("scale", 1); - int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); - int invertsurf = 0; // globflags.GetDefineFlag ("invertsurfacemesh"); + double scale = 1; // globflags.GetNumFlag ("scale", 1); + int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); + int invertsurf = 0; // globflags.GetDefineFlag ("invertsurfacemesh"); - outfile << "mesh3d" << "\n"; + outfile << "mesh3d" << "\n"; - outfile << "dimension\n" << GetDimension() << "\n"; + outfile << "dimension\n" << GetDimension() << "\n"; - outfile << "geomtype\n" << int(geomtype) << "\n"; + outfile << "geomtype\n" << int(geomtype) << "\n"; - outfile << "\n"; - outfile << "# surfnr bcnr domin domout np p1 p2 p3" - << "\n"; + outfile << "\n"; + outfile << "# surfnr bcnr domin domout np p1 p2 p3" + << "\n"; - switch (geomtype) + switch (geomtype) { case GEOM_STL: - outfile << "surfaceelementsgi" << "\n"; - break; + outfile << "surfaceelementsgi" << "\n"; + break; case GEOM_OCC: case GEOM_ACIS: - outfile << "surfaceelementsuv" << "\n"; - break; + outfile << "surfaceelementsuv" << "\n"; + break; default: - outfile << "surfaceelements" << "\n"; + outfile << "surfaceelements" << "\n"; } - outfile << GetNSE() << "\n"; + outfile << GetNSE() << "\n"; - SurfaceElementIndex sei; - for (sei = 0; sei < GetNSE(); sei++) + SurfaceElementIndex sei; + for (sei = 0; sei < GetNSE(); sei++) { - if ((*this)[sei].GetIndex()) - { + if ((*this)[sei].GetIndex()) + { outfile.width(8); outfile << GetFaceDescriptor((*this)[sei].GetIndex ()).SurfNr()+1; outfile.width(8); @@ -423,225 +423,225 @@ namespace netgen outfile << GetFaceDescriptor((*this)[sei].GetIndex ()).DomainIn(); outfile.width(8); outfile << GetFaceDescriptor((*this)[sei].GetIndex ()).DomainOut(); - } - else - outfile << " 0 0 0"; + } + else + outfile << " 0 0 0"; - Element2d sel = (*this)[sei]; - if (invertsurf) - sel.Invert(); + Element2d sel = (*this)[sei]; + if (invertsurf) + sel.Invert(); - outfile.width(8); - outfile << sel.GetNP(); + outfile.width(8); + outfile << sel.GetNP(); - for (j = 0; j < sel.GetNP(); j++) - { + for (j = 0; j < sel.GetNP(); j++) + { outfile.width(8); outfile << sel[j]; - } + } - switch (geomtype) - { - case GEOM_STL: + switch (geomtype) + { + case GEOM_STL: for (j = 1; j <= sel.GetNP(); j++) - { - outfile.width(7); - outfile << " " << sel.GeomInfoPi(j).trignum; - } + { + outfile.width(7); + outfile << " " << sel.GeomInfoPi(j).trignum; + } break; - case GEOM_OCC: case GEOM_ACIS: + case GEOM_OCC: case GEOM_ACIS: for (j = 1; j <= sel.GetNP(); j++) - { - outfile.width(7); - outfile << " " << sel.GeomInfoPi(j).u; - outfile << " " << sel.GeomInfoPi(j).v; - } + { + outfile.width(7); + outfile << " " << sel.GeomInfoPi(j).u; + outfile << " " << sel.GeomInfoPi(j).v; + } break; - default: + default: outfile << "\n"; - } + } - outfile << endl; + outfile << endl; } - outfile << "\n" << "\n"; - outfile << "# matnr np p1 p2 p3 p4" << "\n"; - outfile << "volumeelements" << "\n"; - outfile << GetNE() << "\n"; + outfile << "\n" << "\n"; + outfile << "# matnr np p1 p2 p3 p4" << "\n"; + outfile << "volumeelements" << "\n"; + outfile << GetNE() << "\n"; - for (ElementIndex ei = 0; ei < GetNE(); ei++) + for (ElementIndex ei = 0; ei < GetNE(); ei++) { - outfile.width(8); - outfile << (*this)[ei].GetIndex(); - outfile.width(8); - outfile << (*this)[ei].GetNP(); + outfile.width(8); + outfile << (*this)[ei].GetIndex(); + outfile.width(8); + outfile << (*this)[ei].GetNP(); - Element el = (*this)[ei]; - if (inverttets) - el.Invert(); + Element el = (*this)[ei]; + if (inverttets) + el.Invert(); - /* - for (j = 0; j < el.GetNP(); j++) - for (int k = 0; k < el.GetNP()-1; k++) - if (el[k] > el[k+1]) swap (el[k], el[k+1]); - */ + /* + for (j = 0; j < el.GetNP(); j++) + for (int k = 0; k < el.GetNP()-1; k++) + if (el[k] > el[k+1]) swap (el[k], el[k+1]); + */ - for (j = 0; j < el.GetNP(); j++) - { + for (j = 0; j < el.GetNP(); j++) + { outfile.width(8); outfile << el[j]; - } - outfile << "\n"; + } + outfile << "\n"; } - outfile << "\n" << "\n"; - // outfile << " surf1 surf2 p1 p2" << "\n"; - outfile << "# surfid 0 p1 p2 trignum1 trignum2 domin/surfnr1 domout/surfnr2 ednr1 dist1 ednr2 dist2 \n"; - outfile << "edgesegmentsgi2" << "\n"; - outfile << GetNSeg() << "\n"; + outfile << "\n" << "\n"; + // outfile << " surf1 surf2 p1 p2" << "\n"; + outfile << "# surfid 0 p1 p2 trignum1 trignum2 domin/surfnr1 domout/surfnr2 ednr1 dist1 ednr2 dist2 \n"; + outfile << "edgesegmentsgi2" << "\n"; + outfile << GetNSeg() << "\n"; - for (i = 1; i <= GetNSeg(); i++) + for (i = 1; i <= GetNSeg(); i++) { - const Segment & seg = LineSegment (i); - outfile.width(8); - outfile << seg.si; // 2D: bc number, 3D: wievielte Kante - outfile.width(8); - outfile << 0; - outfile.width(8); - outfile << seg[0]; - outfile.width(8); - outfile << seg[1]; - outfile << " "; - outfile.width(8); - outfile << seg.geominfo[0].trignum; // stl dreiecke - outfile << " "; - outfile.width(8); - outfile << seg.geominfo[1].trignum; // << endl; // stl dreieck + const Segment & seg = LineSegment (i); + outfile.width(8); + outfile << seg.si; // 2D: bc number, 3D: wievielte Kante + outfile.width(8); + outfile << 0; + outfile.width(8); + outfile << seg[0]; + outfile.width(8); + outfile << seg[1]; + outfile << " "; + outfile.width(8); + outfile << seg.geominfo[0].trignum; // stl dreiecke + outfile << " "; + outfile.width(8); + outfile << seg.geominfo[1].trignum; // << endl; // stl dreieck - if (dimension == 3) - { + if (dimension == 3) + { outfile << " "; outfile.width(8); outfile << seg.surfnr1+1; outfile << " "; outfile.width(8); outfile << seg.surfnr2+1; - } - else - { + } + else + { outfile << " "; outfile.width(8); outfile << seg.domin; outfile << " "; outfile.width(8); outfile << seg.domout; - } + } - outfile << " "; - outfile.width(8); - outfile << seg.edgenr; - outfile << " "; - outfile.width(12); - outfile.precision(16); - outfile << seg.epgeominfo[0].dist; // splineparameter (2D) - outfile << " "; - outfile.width(8); - outfile.precision(16); - outfile << seg.epgeominfo[1].edgenr; // geometry dependent - outfile << " "; - outfile.width(12); - outfile << seg.epgeominfo[1].dist; + outfile << " "; + outfile.width(8); + outfile << seg.edgenr; + outfile << " "; + outfile.width(12); + outfile.precision(16); + outfile << seg.epgeominfo[0].dist; // splineparameter (2D) + outfile << " "; + outfile.width(8); + outfile.precision(16); + outfile << seg.epgeominfo[1].edgenr; // geometry dependent + outfile << " "; + outfile.width(12); + outfile << seg.epgeominfo[1].dist; - outfile << "\n"; + outfile << "\n"; } - outfile << "\n" << "\n"; - outfile << "# X Y Z" << "\n"; - outfile << "points" << "\n"; - outfile << GetNP() << "\n"; - outfile.precision(16); - outfile.setf (ios::fixed, ios::floatfield); - outfile.setf (ios::showpoint); + outfile << "\n" << "\n"; + outfile << "# X Y Z" << "\n"; + outfile << "points" << "\n"; + outfile << GetNP() << "\n"; + outfile.precision(16); + outfile.setf (ios::fixed, ios::floatfield); + outfile.setf (ios::showpoint); - PointIndex pi; - for (pi = PointIndex::BASE; + PointIndex pi; + for (pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) { - outfile.width(22); - outfile << (*this)[pi](0)/scale << " "; - outfile.width(22); - outfile << (*this)[pi](1)/scale << " "; - outfile.width(22); - outfile << (*this)[pi](2)/scale << "\n"; + outfile.width(22); + outfile << (*this)[pi](0)/scale << " "; + outfile.width(22); + outfile << (*this)[pi](1)/scale << " "; + outfile.width(22); + outfile << (*this)[pi](2)/scale << "\n"; } - if (ident -> GetMaxNr() > 0) + if (ident -> GetMaxNr() > 0) { - outfile << "identifications\n"; - Array identpairs; - int cnt = 0; - for (i = 1; i <= ident -> GetMaxNr(); i++) - { + outfile << "identifications\n"; + Array identpairs; + int cnt = 0; + for (i = 1; i <= ident -> GetMaxNr(); i++) + { ident -> GetPairs (i, identpairs); cnt += identpairs.Size(); - } - outfile << cnt << "\n"; - for (i = 1; i <= ident -> GetMaxNr(); i++) - { + } + outfile << cnt << "\n"; + for (i = 1; i <= ident -> GetMaxNr(); i++) + { ident -> GetPairs (i, identpairs); for (j = 1; j <= identpairs.Size(); j++) - { - outfile.width (8); - outfile << identpairs.Get(j).I1(); - outfile.width (8); - outfile << identpairs.Get(j).I2(); - outfile.width (8); - outfile << i << "\n"; - } - } + { + outfile.width (8); + outfile << identpairs.Get(j).I1(); + outfile.width (8); + outfile << identpairs.Get(j).I2(); + outfile.width (8); + outfile << i << "\n"; + } + } - outfile << "identificationtypes\n"; - outfile << ident -> GetMaxNr() << "\n"; - for (i = 1; i <= ident -> GetMaxNr(); i++) - { + outfile << "identificationtypes\n"; + outfile << ident -> GetMaxNr() << "\n"; + for (i = 1; i <= ident -> GetMaxNr(); i++) + { int type = ident -> GetType(i); outfile << " " << type; - } - outfile << "\n"; + } + outfile << "\n"; } - int cntmat = 0; - for (i = 1; i <= materials.Size(); i++) - if (materials.Get(i) && strlen (materials.Get(i))) - cntmat++; + int cntmat = 0; + for (i = 1; i <= materials.Size(); i++) + if (materials.Get(i) && strlen (materials.Get(i))) + cntmat++; - if (cntmat) + if (cntmat) { - outfile << "materials" << endl; - outfile << cntmat << endl; - for (i = 1; i <= materials.Size(); i++) - if (materials.Get(i) && strlen (materials.Get(i))) - outfile << i << " " << materials.Get(i) << endl; + outfile << "materials" << endl; + outfile << cntmat << endl; + for (i = 1; i <= materials.Size(); i++) + if (materials.Get(i) && strlen (materials.Get(i))) + outfile << i << " " << materials.Get(i) << endl; } - int cntbcnames = 0; - for ( int ii = 0; ii < bcnames.Size(); ii++ ) - if ( bcnames[ii] ) cntbcnames++; + int cntbcnames = 0; + for ( int ii = 0; ii < bcnames.Size(); ii++ ) + if ( bcnames[ii] ) cntbcnames++; - if ( cntbcnames ) + if ( cntbcnames ) { - outfile << "\n\nbcnames" << endl << bcnames.Size() << endl; - for ( i = 0; i < bcnames.Size(); i++ ) - outfile << i+1 << "\t" << GetBCName(i) << endl; - outfile << endl << endl; + outfile << "\n\nbcnames" << endl << bcnames.Size() << endl; + for ( i = 0; i < bcnames.Size(); i++ ) + outfile << i+1 << "\t" << GetBCName(i) << endl; + outfile << endl << endl; } - /* + /* if ( GetDimension() == 2 ) { for (i = 1; i <= GetNSeg(); i++) @@ -707,224 +707,224 @@ namespace netgen } outfile << endl << endl; } - */ + */ - int cnt_sing = 0; - for (PointIndex pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) - if ((*this)[pi].Singularity()>=1.) cnt_sing++; + int cnt_sing = 0; + for (PointIndex pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) + if ((*this)[pi].Singularity()>=1.) cnt_sing++; - if (cnt_sing) + if (cnt_sing) { - outfile << "singular_points" << endl << cnt_sing << endl; - for (PointIndex pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) - if ((*this)[pi].Singularity()>=1.) - outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl; + outfile << "singular_points" << endl << cnt_sing << endl; + for (PointIndex pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) + if ((*this)[pi].Singularity()>=1.) + outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl; } - cnt_sing = 0; - for (SegmentIndex si = 0; si < GetNSeg(); si++) - if ( segments[si].singedge_left ) cnt_sing++; - if (cnt_sing) + cnt_sing = 0; + for (SegmentIndex si = 0; si < GetNSeg(); si++) + if ( segments[si].singedge_left ) cnt_sing++; + if (cnt_sing) { - outfile << "singular_edge_left" << endl << cnt_sing << endl; - for (SegmentIndex si = 0; si < GetNSeg(); si++) - if ( segments[si].singedge_left ) - outfile << int(si) << "\t" << segments[si].singedge_left << endl; + outfile << "singular_edge_left" << endl << cnt_sing << endl; + for (SegmentIndex si = 0; si < GetNSeg(); si++) + if ( segments[si].singedge_left ) + outfile << int(si) << "\t" << segments[si].singedge_left << endl; } - cnt_sing = 0; - for (SegmentIndex si = 0; si < GetNSeg(); si++) - if ( segments[si].singedge_right ) cnt_sing++; - if (cnt_sing) + cnt_sing = 0; + for (SegmentIndex si = 0; si < GetNSeg(); si++) + if ( segments[si].singedge_right ) cnt_sing++; + if (cnt_sing) { - outfile << "singular_edge_right" << endl << cnt_sing << endl; - for (SegmentIndex si = 0; si < GetNSeg(); si++) - if ( segments[si].singedge_right ) - outfile << int(si) << "\t" << segments[si].singedge_right << endl; + outfile << "singular_edge_right" << endl << cnt_sing << endl; + for (SegmentIndex si = 0; si < GetNSeg(); si++) + if ( segments[si].singedge_right ) + outfile << int(si) << "\t" << segments[si].singedge_right << endl; } - cnt_sing = 0; - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) - if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular) - cnt_sing++; + cnt_sing = 0; + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular) + cnt_sing++; - if (cnt_sing) + if (cnt_sing) { - outfile << "singular_face_inside" << endl << cnt_sing << endl; - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) - if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular) - outfile << int(sei) << "\t" << - GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular << endl; + outfile << "singular_face_inside" << endl << cnt_sing << endl; + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular) + outfile << int(sei) << "\t" << + GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular << endl; } - cnt_sing = 0; - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) - if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular) cnt_sing++; - if (cnt_sing) + cnt_sing = 0; + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular) cnt_sing++; + if (cnt_sing) { - outfile << "singular_face_outside" << endl << cnt_sing << endl; - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) - if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular) - outfile << int(sei) << "\t" - << GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular << endl; + outfile << "singular_face_outside" << endl << cnt_sing << endl; + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular) + outfile << int(sei) << "\t" + << GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular << endl; } - } + } - void Mesh :: Load (const string & filename) - { + void Mesh :: Load (const string & filename) + { - ifstream infile(filename.c_str()); - if (!infile.good()) - throw NgException ("mesh file not found"); + ifstream infile(filename.c_str()); + if (!infile.good()) + throw NgException ("mesh file not found"); - Load(infile); - } + Load(infile); + } - void Mesh :: Load (istream & infile) - { + void Mesh :: Load (istream & infile) + { - char str[100]; - int i, n; + char str[100]; + int i, n; - double scale = 1; // globflags.GetNumFlag ("scale", 1); - int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); - int invertsurf = 0; // globflags.GetDefineFlag ("invertsurfacemesh"); + double scale = 1; // globflags.GetNumFlag ("scale", 1); + int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); + int invertsurf = 0; // globflags.GetDefineFlag ("invertsurfacemesh"); - facedecoding.SetSize(0); + facedecoding.SetSize(0); - bool endmesh = false; + bool endmesh = false; - while (infile.good() && !endmesh) + while (infile.good() && !endmesh) { - infile >> str; + infile >> str; - if (strcmp (str, "dimension") == 0) - { + if (strcmp (str, "dimension") == 0) + { infile >> dimension; - } + } - if (strcmp (str, "geomtype") == 0) - { + if (strcmp (str, "geomtype") == 0) + { int hi; infile >> hi; geomtype = GEOM_TYPE(hi); - } + } - if (strcmp (str, "surfaceelements") == 0 || strcmp (str, "surfaceelementsgi")==0 || strcmp (str, "surfaceelementsuv") == 0) - { + if (strcmp (str, "surfaceelements") == 0 || strcmp (str, "surfaceelementsgi")==0 || strcmp (str, "surfaceelementsuv") == 0) + { infile >> n; PrintMessage (3, n, " surface elements"); for (i = 1; i <= n; i++) - { - int j; - int surfnr, bcp, domin, domout, nep, faceind = 0; + { + int j; + int surfnr, bcp, domin, domout, nep, faceind = 0; - infile >> surfnr >> bcp >> domin >> domout; - surfnr--; + infile >> surfnr >> bcp >> domin >> domout; + surfnr--; - for (j = 1; j <= facedecoding.Size(); j++) + for (j = 1; j <= facedecoding.Size(); j++) if (GetFaceDescriptor(j).SurfNr() == surfnr && - GetFaceDescriptor(j).BCProperty() == bcp && - GetFaceDescriptor(j).DomainIn() == domin && - GetFaceDescriptor(j).DomainOut() == domout) - faceind = j; + GetFaceDescriptor(j).BCProperty() == bcp && + GetFaceDescriptor(j).DomainIn() == domin && + GetFaceDescriptor(j).DomainOut() == domout) + faceind = j; - if (!faceind) - { - faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, 0)); - GetFaceDescriptor(faceind).SetBCProperty (bcp); - } + if (!faceind) + { + faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, 0)); + GetFaceDescriptor(faceind).SetBCProperty (bcp); + } - infile >> nep; - if (!nep) nep = 3; + infile >> nep; + if (!nep) nep = 3; - Element2d tri(nep); - tri.SetIndex(faceind); + Element2d tri(nep); + tri.SetIndex(faceind); - for (j = 1; j <= nep; j++) + for (j = 1; j <= nep; j++) infile >> tri.PNum(j); - if (strcmp (str, "surfaceelementsgi") == 0) + if (strcmp (str, "surfaceelementsgi") == 0) for (j = 1; j <= nep; j++) - infile >> tri.GeomInfoPi(j).trignum; + infile >> tri.GeomInfoPi(j).trignum; - if (strcmp (str, "surfaceelementsuv") == 0) + if (strcmp (str, "surfaceelementsuv") == 0) for (j = 1; j <= nep; j++) - infile >> tri.GeomInfoPi(j).u >> tri.GeomInfoPi(j).v; + infile >> tri.GeomInfoPi(j).u >> tri.GeomInfoPi(j).v; - if (invertsurf) + if (invertsurf) tri.Invert(); - AddSurfaceElement (tri); - } - } + AddSurfaceElement (tri); + } + } - if (strcmp (str, "volumeelements") == 0) - { + if (strcmp (str, "volumeelements") == 0) + { infile >> n; PrintMessage (3, n, " volume elements"); for (i = 1; i <= n; i++) - { - Element el; - int hi, nep; - infile >> hi; - if (hi == 0) hi = 1; - el.SetIndex(hi); - infile >> nep; - el.SetNP(nep); + { + Element el; + int hi, nep; + infile >> hi; + if (hi == 0) hi = 1; + el.SetIndex(hi); + infile >> nep; + el.SetNP(nep); - for (int j = 0; j < nep; j++) + for (int j = 0; j < nep; j++) infile >> (int&)(el[j]); - if (inverttets) + if (inverttets) el.Invert(); - AddVolumeElement (el); - } - } + AddVolumeElement (el); + } + } - if (strcmp (str, "edgesegments") == 0) - { + if (strcmp (str, "edgesegments") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - Segment seg; - int hi; - infile >> seg.si >> hi >> seg[0] >> seg[1]; - AddSegment (seg); - } - } + { + Segment seg; + int hi; + infile >> seg.si >> hi >> seg[0] >> seg[1]; + AddSegment (seg); + } + } - if (strcmp (str, "edgesegmentsgi") == 0) - { + if (strcmp (str, "edgesegmentsgi") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - Segment seg; - int hi; - infile >> seg.si >> hi >> seg[0] >> seg[1] - >> seg.geominfo[0].trignum - >> seg.geominfo[1].trignum; - AddSegment (seg); - } - } + { + Segment seg; + int hi; + infile >> seg.si >> hi >> seg[0] >> seg[1] + >> seg.geominfo[0].trignum + >> seg.geominfo[1].trignum; + AddSegment (seg); + } + } - if (strcmp (str, "edgesegmentsgi2") == 0) - { + if (strcmp (str, "edgesegmentsgi2") == 0) + { int a; infile >> a; n=a; @@ -932,451 +932,450 @@ namespace netgen PrintMessage (3, n, " curve elements"); for (i = 1; i <= n; i++) - { - Segment seg; - int hi; - infile >> seg.si >> hi >> seg[0] >> seg[1] - >> seg.geominfo[0].trignum - >> seg.geominfo[1].trignum - >> seg.surfnr1 >> seg.surfnr2 - >> seg.edgenr - >> seg.epgeominfo[0].dist - >> seg.epgeominfo[1].edgenr - >> seg.epgeominfo[1].dist; + { + Segment seg; + int hi; + infile >> seg.si >> hi >> seg[0] >> seg[1] + >> seg.geominfo[0].trignum + >> seg.geominfo[1].trignum + >> seg.surfnr1 >> seg.surfnr2 + >> seg.edgenr + >> seg.epgeominfo[0].dist + >> seg.epgeominfo[1].edgenr + >> seg.epgeominfo[1].dist; - seg.epgeominfo[0].edgenr = seg.epgeominfo[1].edgenr; + seg.epgeominfo[0].edgenr = seg.epgeominfo[1].edgenr; - seg.domin = seg.surfnr1; - seg.domout = seg.surfnr2; + seg.domin = seg.surfnr1; + seg.domout = seg.surfnr2; - seg.surfnr1--; - seg.surfnr2--; + seg.surfnr1--; + seg.surfnr2--; - AddSegment (seg); - } - } + AddSegment (seg); + } + } - if (strcmp (str, "points") == 0) - { + if (strcmp (str, "points") == 0) + { infile >> n; PrintMessage (3, n, " points"); for (i = 1; i <= n; i++) - { - Point3d p; - infile >> p.X() >> p.Y() >> p.Z(); - p.X() *= scale; - p.Y() *= scale; - p.Z() *= scale; - AddPoint (p); - } - } + { + Point3d p; + infile >> p.X() >> p.Y() >> p.Z(); + p.X() *= scale; + p.Y() *= scale; + p.Z() *= scale; + AddPoint (p); + } + } - if (strcmp (str, "identifications") == 0) - { + if (strcmp (str, "identifications") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - PointIndex pi1, pi2; - int ind; - infile >> pi1 >> pi2 >> ind; - ident -> Add (pi1, pi2, ind); - } - } + { + PointIndex pi1, pi2; + int ind; + infile >> pi1 >> pi2 >> ind; + ident -> Add (pi1, pi2, ind); + } + } - if (strcmp (str, "identificationtypes") == 0) - { + if (strcmp (str, "identificationtypes") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - int type; - infile >> type; - ident -> SetType(i,Identifications::ID_TYPE(type)); - } - } + { + int type; + infile >> type; + ident -> SetType(i,Identifications::ID_TYPE(type)); + } + } - if (strcmp (str, "materials") == 0) - { + if (strcmp (str, "materials") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - int nr; - string mat; - infile >> nr >> mat; - SetMaterial (nr, mat.c_str()); - } - } + { + int nr; + string mat; + infile >> nr >> mat; + SetMaterial (nr, mat.c_str()); + } + } - if ( strcmp (str, "bcnames" ) == 0 ) - { + if ( strcmp (str, "bcnames" ) == 0 ) + { infile >> n; Array bcnrs(n); - SetNBCNames(n); for ( i = 1; i <= n; i++ ) - { - string nextbcname; - infile >> bcnrs[i-1] >> nextbcname; - bcnames[bcnrs[i-1]-1] = new string(nextbcname); - } + { + string nextbcname; + infile >> bcnrs[i-1] >> nextbcname; + bcnames[bcnrs[i-1]-1] = new string(nextbcname); + } if ( GetDimension() == 2 ) - { - for (i = 1; i <= GetNSeg(); i++) - { - Segment & seg = LineSegment (i); - if ( seg.si <= n ) - seg.SetBCName (bcnames[seg.si-1]); - else - seg.SetBCName(0); - } - } - else - { - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) - { - if ((*this)[sei].GetIndex()) + { + for (i = 1; i <= GetNSeg(); i++) { - int bcp = GetFaceDescriptor((*this)[sei].GetIndex ()).BCProperty(); - if ( bcp <= n ) - GetFaceDescriptor((*this)[sei].GetIndex ()).SetBCName(bcnames[bcp-1]); - else - GetFaceDescriptor((*this)[sei].GetIndex ()).SetBCName(0); - + Segment & seg = LineSegment (i); + if ( seg.si <= n ) + seg.SetBCName (bcnames[seg.si-1]); + else + seg.SetBCName(0); } - } + } + else + { + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + { + if ((*this)[sei].GetIndex()) + { + int bcp = GetFaceDescriptor((*this)[sei].GetIndex ()).BCProperty(); + if ( bcp <= n ) + GetFaceDescriptor((*this)[sei].GetIndex ()).SetBCName(bcnames[bcp-1]); + else + GetFaceDescriptor((*this)[sei].GetIndex ()).SetBCName(0); - } + } + } + + } - } + } - if (strcmp (str, "singular_points") == 0) - { + if (strcmp (str, "singular_points") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - PointIndex pi; - double s; - infile >> pi; - infile >> s; - (*this)[pi].Singularity (s); - } - } + { + PointIndex pi; + double s; + infile >> pi; + infile >> s; + (*this)[pi].Singularity (s); + } + } - if (strcmp (str, "singular_edge_left") == 0) - { + if (strcmp (str, "singular_edge_left") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - SegmentIndex si; - double s; - infile >> si; - infile >> s; - (*this)[si].singedge_left = s; - } - } - if (strcmp (str, "singular_edge_right") == 0) - { + { + SegmentIndex si; + double s; + infile >> si; + infile >> s; + (*this)[si].singedge_left = s; + } + } + if (strcmp (str, "singular_edge_right") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - SegmentIndex si; - double s; - infile >> si; - infile >> s; - (*this)[si].singedge_right = s; - } - } + { + SegmentIndex si; + double s; + infile >> si; + infile >> s; + (*this)[si].singedge_right = s; + } + } - if (strcmp (str, "singular_face_inside") == 0) - { + if (strcmp (str, "singular_face_inside") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - SurfaceElementIndex sei; - double s; - infile >> sei; - infile >> s; - GetFaceDescriptor((*this)[sei].GetIndex()).domin_singular = s; - } - } + { + SurfaceElementIndex sei; + double s; + infile >> sei; + infile >> s; + GetFaceDescriptor((*this)[sei].GetIndex()).domin_singular = s; + } + } - if (strcmp (str, "singular_face_outside") == 0) - { + if (strcmp (str, "singular_face_outside") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - SurfaceElementIndex sei; - double s; - infile >> sei; - infile >> s; - GetFaceDescriptor((*this)[sei].GetIndex()).domout_singular = s; - } - } + { + SurfaceElementIndex sei; + double s; + infile >> sei; + infile >> s; + GetFaceDescriptor((*this)[sei].GetIndex()).domout_singular = s; + } + } - if (strcmp (str, "endmesh") == 0) - endmesh = true; + if (strcmp (str, "endmesh") == 0) + endmesh = true; - strcpy (str, ""); + strcpy (str, ""); } - CalcSurfacesOfNode (); - // BuildConnectedNodes (); - topology -> Update(); - clusters -> Update(); + CalcSurfacesOfNode (); + // BuildConnectedNodes (); + topology -> Update(); + clusters -> Update(); - SetNextMajorTimeStamp(); - // PrintMemInfo (cout); + SetNextMajorTimeStamp(); + // PrintMemInfo (cout); #ifdef PARALLEL - if ( ntasks > 1 ) + if ( ntasks > 1 ) { - // for parallel processing - Distribute (); - return; + // for parallel processing + Distribute (); + return; } #endif - } + } - void Mesh :: Merge (const string & filename, const int surfindex_offset) - { - ifstream infile(filename.c_str()); - if (!infile.good()) - throw NgException ("mesh file not found"); + void Mesh :: Merge (const string & filename, const int surfindex_offset) + { + ifstream infile(filename.c_str()); + if (!infile.good()) + throw NgException ("mesh file not found"); - Merge(infile,surfindex_offset); + Merge(infile,surfindex_offset); - } + } - void Mesh :: Merge (istream & infile, const int surfindex_offset) - { - char str[100]; - int i, n; + void Mesh :: Merge (istream & infile, const int surfindex_offset) + { + char str[100]; + int i, n; - int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); + int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); - int oldnp = GetNP(); - int oldne = GetNSeg(); - int oldnd = GetNDomains(); + int oldnp = GetNP(); + int oldne = GetNSeg(); + int oldnd = GetNDomains(); - for(SurfaceElementIndex si = 0; si < GetNSE(); si++) - for(int j=1; j<=(*this)[si].GetNP(); j++) (*this)[si].GeomInfoPi(j).trignum = -1; + for(SurfaceElementIndex si = 0; si < GetNSE(); si++) + for(int j=1; j<=(*this)[si].GetNP(); j++) (*this)[si].GeomInfoPi(j).trignum = -1; - int max_surfnr = 0; - for (i = 1; i <= GetNFD(); i++) - max_surfnr = max2 (max_surfnr, GetFaceDescriptor(i).SurfNr()); - max_surfnr++; + int max_surfnr = 0; + for (i = 1; i <= GetNFD(); i++) + max_surfnr = max2 (max_surfnr, GetFaceDescriptor(i).SurfNr()); + max_surfnr++; - if(max_surfnr < surfindex_offset) max_surfnr = surfindex_offset; + if(max_surfnr < surfindex_offset) max_surfnr = surfindex_offset; - bool endmesh = false; + bool endmesh = false; - while (infile.good() && !endmesh) + while (infile.good() && !endmesh) { - infile >> str; + infile >> str; - if (strcmp (str, "surfaceelementsgi") == 0 || strcmp (str, "surfaceelements") == 0) - { + if (strcmp (str, "surfaceelementsgi") == 0 || strcmp (str, "surfaceelements") == 0) + { infile >> n; PrintMessage (3, n, " surface elements"); for (i = 1; i <= n; i++) - { - int j; - int surfnr, bcp, domin, domout, nep, faceind = 0; - infile >> surfnr >> bcp >> domin >> domout; + { + int j; + int surfnr, bcp, domin, domout, nep, faceind = 0; + infile >> surfnr >> bcp >> domin >> domout; - surfnr--; + surfnr--; - if(domin > 0) domin += oldnd; - if(domout > 0) domout += oldnd; - surfnr += max_surfnr; + if(domin > 0) domin += oldnd; + if(domout > 0) domout += oldnd; + surfnr += max_surfnr; - for (j = 1; j <= facedecoding.Size(); j++) + for (j = 1; j <= facedecoding.Size(); j++) if (GetFaceDescriptor(j).SurfNr() == surfnr && - GetFaceDescriptor(j).BCProperty() == bcp && - GetFaceDescriptor(j).DomainIn() == domin && - GetFaceDescriptor(j).DomainOut() == domout) - faceind = j; + GetFaceDescriptor(j).BCProperty() == bcp && + GetFaceDescriptor(j).DomainIn() == domin && + GetFaceDescriptor(j).DomainOut() == domout) + faceind = j; - if (!faceind) - { - faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, 0)); - if(GetDimension() == 2) bcp++; - GetFaceDescriptor(faceind).SetBCProperty (bcp); - } - - infile >> nep; - if (!nep) nep = 3; - - Element2d tri(nep); - tri.SetIndex(faceind); - - for (j = 1; j <= nep; j++) - { - infile >> tri.PNum(j); - tri.PNum(j) = tri.PNum(j) + oldnp; - } - - - if (strcmp (str, "surfaceelementsgi") == 0) - for (j = 1; j <= nep; j++) + if (!faceind) { - infile >> tri.GeomInfoPi(j).trignum; - tri.GeomInfoPi(j).trignum = -1; + faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, 0)); + if(GetDimension() == 2) bcp++; + GetFaceDescriptor(faceind).SetBCProperty (bcp); } - AddSurfaceElement (tri); - } - } + infile >> nep; + if (!nep) nep = 3; + + Element2d tri(nep); + tri.SetIndex(faceind); + + for (j = 1; j <= nep; j++) + { + infile >> tri.PNum(j); + tri.PNum(j) = tri.PNum(j) + oldnp; + } - if (strcmp (str, "edgesegments") == 0) - { + if (strcmp (str, "surfaceelementsgi") == 0) + for (j = 1; j <= nep; j++) + { + infile >> tri.GeomInfoPi(j).trignum; + tri.GeomInfoPi(j).trignum = -1; + } + + AddSurfaceElement (tri); + } + } + + + if (strcmp (str, "edgesegments") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - Segment seg; - int hi; - infile >> seg.si >> hi >> seg[0] >> seg[1]; - seg[0] = seg[0] + oldnp; - seg[1] = seg[1] + oldnp; - AddSegment (seg); - } - } + { + Segment seg; + int hi; + infile >> seg.si >> hi >> seg[0] >> seg[1]; + seg[0] = seg[0] + oldnp; + seg[1] = seg[1] + oldnp; + AddSegment (seg); + } + } - if (strcmp (str, "edgesegmentsgi") == 0) - { + if (strcmp (str, "edgesegmentsgi") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - Segment seg; - int hi; - infile >> seg.si >> hi >> seg[0] >> seg[1] - >> seg.geominfo[0].trignum - >> seg.geominfo[1].trignum; - seg[0] = seg[0] + oldnp; - seg[1] = seg[1] + oldnp; - AddSegment (seg); - } - } - if (strcmp (str, "edgesegmentsgi2") == 0) - { + { + Segment seg; + int hi; + infile >> seg.si >> hi >> seg[0] >> seg[1] + >> seg.geominfo[0].trignum + >> seg.geominfo[1].trignum; + seg[0] = seg[0] + oldnp; + seg[1] = seg[1] + oldnp; + AddSegment (seg); + } + } + if (strcmp (str, "edgesegmentsgi2") == 0) + { infile >> n; PrintMessage (3, n, " curve elements"); for (i = 1; i <= n; i++) - { - Segment seg; - int hi; - infile >> seg.si >> hi >> seg[0] >> seg[1] - >> seg.geominfo[0].trignum - >> seg.geominfo[1].trignum - >> seg.surfnr1 >> seg.surfnr2 - >> seg.edgenr - >> seg.epgeominfo[0].dist - >> seg.epgeominfo[1].edgenr - >> seg.epgeominfo[1].dist; - seg.epgeominfo[0].edgenr = seg.epgeominfo[1].edgenr; + { + Segment seg; + int hi; + infile >> seg.si >> hi >> seg[0] >> seg[1] + >> seg.geominfo[0].trignum + >> seg.geominfo[1].trignum + >> seg.surfnr1 >> seg.surfnr2 + >> seg.edgenr + >> seg.epgeominfo[0].dist + >> seg.epgeominfo[1].edgenr + >> seg.epgeominfo[1].dist; + seg.epgeominfo[0].edgenr = seg.epgeominfo[1].edgenr; - seg.surfnr1--; - seg.surfnr2--; + seg.surfnr1--; + seg.surfnr2--; - if(seg.surfnr1 >= 0) seg.surfnr1 = seg.surfnr1 + max_surfnr; - if(seg.surfnr2 >= 0) seg.surfnr2 = seg.surfnr2 + max_surfnr; - seg[0] = seg[0] +oldnp; - seg[1] = seg[1] +oldnp; - seg.edgenr = seg.edgenr + oldne; - seg.epgeominfo[1].edgenr = seg.epgeominfo[1].edgenr + oldne; + if(seg.surfnr1 >= 0) seg.surfnr1 = seg.surfnr1 + max_surfnr; + if(seg.surfnr2 >= 0) seg.surfnr2 = seg.surfnr2 + max_surfnr; + seg[0] = seg[0] +oldnp; + seg[1] = seg[1] +oldnp; + seg.edgenr = seg.edgenr + oldne; + seg.epgeominfo[1].edgenr = seg.epgeominfo[1].edgenr + oldne; - AddSegment (seg); - } - } + AddSegment (seg); + } + } - if (strcmp (str, "volumeelements") == 0) - { + if (strcmp (str, "volumeelements") == 0) + { infile >> n; PrintMessage (3, n, " volume elements"); for (i = 1; i <= n; i++) - { - Element el; - int hi, nep; - infile >> hi; - if (hi == 0) hi = 1; - el.SetIndex(hi+oldnd); - infile >> nep; - el.SetNP(nep); + { + Element el; + int hi, nep; + infile >> hi; + if (hi == 0) hi = 1; + el.SetIndex(hi+oldnd); + infile >> nep; + el.SetNP(nep); - for (int j = 0; j < nep; j++) - { - infile >> (int&)(el[j]); - el[j] = el[j]+oldnp; - } + for (int j = 0; j < nep; j++) + { + infile >> (int&)(el[j]); + el[j] = el[j]+oldnp; + } - if (inverttets) + if (inverttets) el.Invert(); - AddVolumeElement (el); - } - } + AddVolumeElement (el); + } + } - if (strcmp (str, "points") == 0) - { + if (strcmp (str, "points") == 0) + { infile >> n; PrintMessage (3, n, " points"); for (i = 1; i <= n; i++) - { - Point3d p; - infile >> p.X() >> p.Y() >> p.Z(); - AddPoint (p); - } - } + { + Point3d p; + infile >> p.X() >> p.Y() >> p.Z(); + AddPoint (p); + } + } - if (strcmp (str, "endmesh") == 0) - { + if (strcmp (str, "endmesh") == 0) + { endmesh = true; - } + } - if (strcmp (str, "materials") == 0) - { + if (strcmp (str, "materials") == 0) + { infile >> n; for (i = 1; i <= n; i++) - { - int nr; - string mat; - infile >> nr >> mat; - SetMaterial (nr+oldnd, mat.c_str()); - } - } + { + int nr; + string mat; + infile >> nr >> mat; + SetMaterial (nr+oldnd, mat.c_str()); + } + } - strcpy (str, ""); + strcpy (str, ""); } - CalcSurfacesOfNode (); + CalcSurfacesOfNode (); - topology -> Update(); - clusters -> Update(); + topology -> Update(); + clusters -> Update(); - SetNextMajorTimeStamp(); - } + SetNextMajorTimeStamp(); + } @@ -1387,64 +1386,64 @@ namespace netgen - bool Mesh :: TestOk () const - { - for (ElementIndex ei = 0; ei < volelements.Size(); ei++) + bool Mesh :: TestOk () const + { + for (ElementIndex ei = 0; ei < volelements.Size(); ei++) { - for (int j = 0; j < 4; j++) - if ( (*this)[ei][j] <= PointIndex::BASE-1) + for (int j = 0; j < 4; j++) + if ( (*this)[ei][j] <= PointIndex::BASE-1) { - (*testout) << "El " << ei << " has 0 nodes: "; - for (int k = 0; k < 4; k++) - (*testout) << (*this)[ei][k]; - break; + (*testout) << "El " << ei << " has 0 nodes: "; + for (int k = 0; k < 4; k++) + (*testout) << (*this)[ei][k]; + break; } } - CheckMesh3D (*this); - return 1; - } + CheckMesh3D (*this); + return 1; + } - void Mesh :: SetAllocSize(int nnodes, int nsegs, int nsel, int nel) - { - points.SetAllocSize(nnodes); - segments.SetAllocSize(nsegs); - surfelements.SetAllocSize(nsel); - volelements.SetAllocSize(nel); - } + void Mesh :: SetAllocSize(int nnodes, int nsegs, int nsel, int nel) + { + points.SetAllocSize(nnodes); + segments.SetAllocSize(nsegs); + surfelements.SetAllocSize(nsel); + volelements.SetAllocSize(nel); + } - void Mesh :: BuildBoundaryEdges(void) - { - delete boundaryedges; + void Mesh :: BuildBoundaryEdges(void) + { + delete boundaryedges; - boundaryedges = new INDEX_2_CLOSED_HASHTABLE - (3 * (GetNSE() + GetNOpenElements()) + GetNSeg() + 1); + boundaryedges = new INDEX_2_CLOSED_HASHTABLE + (3 * (GetNSE() + GetNOpenElements()) + GetNSeg() + 1); - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { - const Element2d & sel = surfelements[sei]; - if (sel.IsDeleted()) continue; + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; - int si = sel.GetIndex(); + // int si = sel.GetIndex(); - for (int j = 0; j < sel.GetNP(); j++) - { + for (int j = 0; j < sel.GetNP(); j++) + { INDEX_2 i2; i2.I1() = sel.PNumMod(j+1); i2.I2() = sel.PNumMod(j+2); i2.Sort(); if (sel.GetNP() <= 4) - boundaryedges->Set (i2, 1); - } + boundaryedges->Set (i2, 1); + } } - for (int i = 0; i < openelements.Size(); i++) + for (int i = 0; i < openelements.Size(); i++) { - const Element2d & sel = openelements[i]; - for (int j = 0; j < sel.GetNP(); j++) - { + const Element2d & sel = openelements[i]; + for (int j = 0; j < sel.GetNP(); j++) + { INDEX_2 i2; i2.I1() = sel.PNumMod(j+1); i2.I2() = sel.PNumMod(j+2); @@ -1452,67 +1451,67 @@ namespace netgen boundaryedges->Set (i2, 1); points[sel[j]].SetType(FIXEDPOINT); - } + } } - for (int i = 0; i < GetNSeg(); i++) + for (int i = 0; i < GetNSeg(); i++) { - const Segment & seg = segments[i]; - INDEX_2 i2(seg[0], seg[1]); - i2.Sort(); + const Segment & seg = segments[i]; + INDEX_2 i2(seg[0], seg[1]); + i2.Sort(); - boundaryedges -> Set (i2, 2); - //segmentht -> Set (i2, i); + boundaryedges -> Set (i2, 2); + //segmentht -> Set (i2, i); } - } + } - void Mesh :: CalcSurfacesOfNode () - { - int i, j, k; - SurfaceElementIndex sei; + void Mesh :: CalcSurfacesOfNode () + { + int i, j, k; + SurfaceElementIndex sei; - surfacesonnode.SetSize (GetNP()); + surfacesonnode.SetSize (GetNP()); - delete boundaryedges; - boundaryedges = NULL; + delete boundaryedges; + boundaryedges = NULL; - delete surfelementht; - delete segmentht; + delete surfelementht; + delete segmentht; - /* + /* surfelementht = new INDEX_3_HASHTABLE (GetNSE()/4 + 1); segmentht = new INDEX_2_HASHTABLE (GetNSeg() + 1); - */ + */ - surfelementht = new INDEX_3_CLOSED_HASHTABLE (3*GetNSE() + 1); - segmentht = new INDEX_2_CLOSED_HASHTABLE (3*GetNSeg() + 1); + surfelementht = new INDEX_3_CLOSED_HASHTABLE (3*GetNSE() + 1); + segmentht = new INDEX_2_CLOSED_HASHTABLE (3*GetNSeg() + 1); - for (sei = 0; sei < GetNSE(); sei++) + for (sei = 0; sei < GetNSE(); sei++) { - const Element2d & sel = surfelements[sei]; - if (sel.IsDeleted()) continue; + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; - int si = sel.GetIndex(); + int si = sel.GetIndex(); - for (j = 0; j < sel.GetNP(); j++) - { + for (j = 0; j < sel.GetNP(); j++) + { PointIndex pi = sel[j]; bool found = 0; for (k = 0; k < surfacesonnode[pi].Size(); k++) - if (surfacesonnode[pi][k] == si) - { + if (surfacesonnode[pi][k] == si) + { found = 1; break; - } + } - if (!found) - surfacesonnode.Add (pi, si); + if (!found) + surfacesonnode.Add (pi, si); - } + } } - /* + /* for (sei = 0; sei < GetNSE(); sei++) { const Element2d & sel = surfelements[sei]; @@ -1527,81 +1526,81 @@ namespace netgen } surfelementht -> AllocateElements(); - */ - for (sei = 0; sei < GetNSE(); sei++) + */ + for (sei = 0; sei < GetNSE(); sei++) { - const Element2d & sel = surfelements[sei]; - if (sel.IsDeleted()) continue; + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; - INDEX_3 i3; - i3.I1() = sel.PNum(1); - i3.I2() = sel.PNum(2); - i3.I3() = sel.PNum(3); - i3.Sort(); - surfelementht -> Set (i3, sei); // war das wichtig ??? sel.GetIndex()); + INDEX_3 i3; + i3.I1() = sel.PNum(1); + i3.I2() = sel.PNum(2); + i3.I3() = sel.PNum(3); + i3.Sort(); + surfelementht -> Set (i3, sei); // war das wichtig ??? sel.GetIndex()); } - int np = GetNP(); + int np = GetNP(); - if (dimension == 3) + if (dimension == 3) { - for (PointIndex pi = PointIndex::BASE; - pi < np+PointIndex::BASE; pi++) - points[pi].SetType (INNERPOINT); + for (PointIndex pi = PointIndex::BASE; + pi < np+PointIndex::BASE; pi++) + points[pi].SetType (INNERPOINT); - if (GetNFD() == 0) - { + if (GetNFD() == 0) + { for (sei = 0; sei < GetNSE(); sei++) - { - const Element2d & sel = surfelements[sei]; - if (sel.IsDeleted()) continue; - for (j = 0; j < sel.GetNP(); j++) - { - PointIndex pi = SurfaceElement(sei)[j]; - points[pi].SetType(FIXEDPOINT); - } - } - } - else - { + { + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; + for (j = 0; j < sel.GetNP(); j++) + { + PointIndex pi = SurfaceElement(sei)[j]; + points[pi].SetType(FIXEDPOINT); + } + } + } + else + { for (sei = 0; sei < GetNSE(); sei++) - { - const Element2d & sel = surfelements[sei]; - if (sel.IsDeleted()) continue; - for (j = 0; j < sel.GetNP(); j++) - { - PointIndex pi = sel[j]; - int ns = surfacesonnode[pi].Size(); - if (ns == 1) - points[pi].SetType(SURFACEPOINT); - if (ns == 2) - points[pi].SetType(EDGEPOINT); - if (ns >= 3) - points[pi].SetType(FIXEDPOINT); - } - } - } + { + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; + for (j = 0; j < sel.GetNP(); j++) + { + PointIndex pi = sel[j]; + int ns = surfacesonnode[pi].Size(); + if (ns == 1) + points[pi].SetType(SURFACEPOINT); + if (ns == 2) + points[pi].SetType(EDGEPOINT); + if (ns >= 3) + points[pi].SetType(FIXEDPOINT); + } + } + } - for (i = 0; i < segments.Size(); i++) - { + for (i = 0; i < segments.Size(); i++) + { const Segment & seg = segments[i]; for (j = 1; j <= 2; j++) - { - PointIndex hi = (j == 1) ? seg[0] : seg[1]; + { + PointIndex hi = (j == 1) ? seg[0] : seg[1]; - if (points[hi].Type() == INNERPOINT || - points[hi].Type() == SURFACEPOINT) + if (points[hi].Type() == INNERPOINT || + points[hi].Type() == SURFACEPOINT) points[hi].SetType(EDGEPOINT); - } - } + } + } - for (i = 0; i < lockedpoints.Size(); i++) - points[lockedpoints[i]].SetType(FIXEDPOINT); + for (i = 0; i < lockedpoints.Size(); i++) + points[lockedpoints[i]].SetType(FIXEDPOINT); } - /* + /* for (i = 0; i < openelements.Size(); i++) { const Element2d & sel = openelements[i]; @@ -1616,938 +1615,938 @@ namespace netgen points[sel[j]].SetType(FIXEDPOINT); } } - */ + */ - // eltyps.SetSize (GetNE()); - // eltyps = FREEELEMENT; + // eltyps.SetSize (GetNE()); + // eltyps = FREEELEMENT; - for (i = 0; i < GetNSeg(); i++) + for (i = 0; i < GetNSeg(); i++) { - const Segment & seg = segments[i]; - INDEX_2 i2(seg[0], seg[1]); - i2.Sort(); + const Segment & seg = segments[i]; + INDEX_2 i2(seg[0], seg[1]); + i2.Sort(); - //boundaryedges -> Set (i2, 2); - segmentht -> Set (i2, i); + //boundaryedges -> Set (i2, 2); + segmentht -> Set (i2, i); } - } + } - void Mesh :: FixPoints (const BitArray & fixpoints) - { - if (fixpoints.Size() != GetNP()) + void Mesh :: FixPoints (const BitArray & fixpoints) + { + if (fixpoints.Size() != GetNP()) { - cerr << "Mesh::FixPoints: sizes don't fit" << endl; - return; + cerr << "Mesh::FixPoints: sizes don't fit" << endl; + return; } - int np = GetNP(); - for (int i = 1; i <= np; i++) - if (fixpoints.Test(i)) - { - points.Elem(i).SetType (FIXEDPOINT); - } - } + int np = GetNP(); + for (int i = 1; i <= np; i++) + if (fixpoints.Test(i)) + { + points.Elem(i).SetType (FIXEDPOINT); + } + } - void Mesh :: FindOpenElements (int dom) - { - static int timer = NgProfiler::CreateTimer ("Mesh::FindOpenElements"); - static int timera = NgProfiler::CreateTimer ("Mesh::FindOpenElements A"); - static int timerb = NgProfiler::CreateTimer ("Mesh::FindOpenElements B"); - static int timerc = NgProfiler::CreateTimer ("Mesh::FindOpenElements C"); - static int timerd = NgProfiler::CreateTimer ("Mesh::FindOpenElements D"); - static int timere = NgProfiler::CreateTimer ("Mesh::FindOpenElements E"); + void Mesh :: FindOpenElements (int dom) + { + static int timer = NgProfiler::CreateTimer ("Mesh::FindOpenElements"); + static int timera = NgProfiler::CreateTimer ("Mesh::FindOpenElements A"); + static int timerb = NgProfiler::CreateTimer ("Mesh::FindOpenElements B"); + static int timerc = NgProfiler::CreateTimer ("Mesh::FindOpenElements C"); + static int timerd = NgProfiler::CreateTimer ("Mesh::FindOpenElements D"); + static int timere = NgProfiler::CreateTimer ("Mesh::FindOpenElements E"); - NgProfiler::RegionTimer reg (timer); + NgProfiler::RegionTimer reg (timer); - int np = GetNP(); - int ne = GetNE(); - int nse = GetNSE(); + int np = GetNP(); + int ne = GetNE(); + int nse = GetNSE(); - Array numonpoint(np); + Array numonpoint(np); - numonpoint = 0; + numonpoint = 0; - NgProfiler::StartTimer (timera); - for (ElementIndex ei = 0; ei < ne; ei++) + NgProfiler::StartTimer (timera); + for (ElementIndex ei = 0; ei < ne; ei++) { - const Element & el = (*this)[ei]; - if (dom == 0 || dom == el.GetIndex()) - { + const Element & el = (*this)[ei]; + if (dom == 0 || dom == el.GetIndex()) + { if (el.GetNP() == 4) - { - INDEX_4 i4(el[0], el[1], el[2], el[3]); - i4.Sort(); - numonpoint[i4.I1()]++; - numonpoint[i4.I2()]++; - } + { + INDEX_4 i4(el[0], el[1], el[2], el[3]); + i4.Sort(); + numonpoint[i4.I1()]++; + numonpoint[i4.I2()]++; + } else - for (int j = 0; j < el.GetNP(); j++) - numonpoint[el[j]]++; - } + for (int j = 0; j < el.GetNP(); j++) + numonpoint[el[j]]++; + } } - TABLE elsonpoint(numonpoint); - for (ElementIndex ei = 0; ei < ne; ei++) + TABLE elsonpoint(numonpoint); + for (ElementIndex ei = 0; ei < ne; ei++) { - const Element & el = (*this)[ei]; - if (dom == 0 || dom == el.GetIndex()) - { + const Element & el = (*this)[ei]; + if (dom == 0 || dom == el.GetIndex()) + { if (el.GetNP() == 4) - { - INDEX_4 i4(el[0], el[1], el[2], el[3]); - i4.Sort(); - elsonpoint.Add (i4.I1(), ei); - elsonpoint.Add (i4.I2(), ei); - } + { + INDEX_4 i4(el[0], el[1], el[2], el[3]); + i4.Sort(); + elsonpoint.Add (i4.I1(), ei); + elsonpoint.Add (i4.I2(), ei); + } else - for (int j = 0; j < el.GetNP(); j++) - elsonpoint.Add (el[j], ei); - } + for (int j = 0; j < el.GetNP(); j++) + elsonpoint.Add (el[j], ei); + } } - NgProfiler::StopTimer (timera); + NgProfiler::StopTimer (timera); - NgProfiler::StartTimer (timerb); + NgProfiler::StartTimer (timerb); - Array hasface(GetNFD()); + Array hasface(GetNFD()); - int i; - for (i = 1; i <= GetNFD(); i++) + int i; + for (i = 1; i <= GetNFD(); i++) { - int domin = GetFaceDescriptor(i).DomainIn(); - int domout = GetFaceDescriptor(i).DomainOut(); - hasface[i] = - dom == 0 && (domin != 0 || domout != 0) || - dom != 0 && (domin == dom || domout == dom); + int domin = GetFaceDescriptor(i).DomainIn(); + int domout = GetFaceDescriptor(i).DomainOut(); + hasface[i] = + ( dom == 0 && (domin != 0 || domout != 0) ) || + ( dom != 0 && (domin == dom || domout == dom) ); } - numonpoint = 0; - for (SurfaceElementIndex sii = 0; sii < nse; sii++) + numonpoint = 0; + for (SurfaceElementIndex sii = 0; sii < nse; sii++) { - int ind = surfelements[sii].GetIndex(); - /* - if ( - GetFaceDescriptor(ind).DomainIn() && - (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) - || - GetFaceDescriptor(ind).DomainOut() && - (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) - ) - */ - if (hasface[ind]) - { + int ind = surfelements[sii].GetIndex(); + /* + if ( + GetFaceDescriptor(ind).DomainIn() && + (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) + || + GetFaceDescriptor(ind).DomainOut() && + (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) + ) + */ + if (hasface[ind]) + { /* - Element2d hel = surfelements[i]; - hel.NormalizeNumbering(); - numonpoint[hel[0]]++; + Element2d hel = surfelements[i]; + hel.NormalizeNumbering(); + numonpoint[hel[0]]++; */ const Element2d & hel = surfelements[sii]; int mini = 0; for (int j = 1; j < hel.GetNP(); j++) - if (hel[j] < hel[mini]) - mini = j; + if (hel[j] < hel[mini]) + mini = j; numonpoint[hel[mini]]++; - } + } } - TABLE selsonpoint(numonpoint); - for (SurfaceElementIndex sii = 0; sii < nse; sii++) + TABLE selsonpoint(numonpoint); + for (SurfaceElementIndex sii = 0; sii < nse; sii++) { - int ind = surfelements[sii].GetIndex(); + int ind = surfelements[sii].GetIndex(); - /* - if ( - GetFaceDescriptor(ind).DomainIn() && - (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) - || - GetFaceDescriptor(ind).DomainOut() && - (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) - ) - */ - if (hasface[ind]) - { + /* + if ( + GetFaceDescriptor(ind).DomainIn() && + (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) + || + GetFaceDescriptor(ind).DomainOut() && + (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) + ) + */ + if (hasface[ind]) + { /* - Element2d hel = surfelements[i]; - hel.NormalizeNumbering(); - selsonpoint.Add (hel[0], i); + Element2d hel = surfelements[i]; + hel.NormalizeNumbering(); + selsonpoint.Add (hel[0], i); */ const Element2d & hel = surfelements[sii]; int mini = 0; for (int j = 1; j < hel.GetNP(); j++) - if (hel[j] < hel[mini]) - mini = j; + if (hel[j] < hel[mini]) + mini = j; selsonpoint.Add (hel[mini], sii); - } + } } - NgProfiler::StopTimer (timerb); + NgProfiler::StopTimer (timerb); - int ii, j, k, l; - PointIndex pi; - SurfaceElementIndex sei; - Element2d hel; + int ii; + PointIndex pi; + SurfaceElementIndex sei; + Element2d hel; - NgProfiler::RegionTimer regc (timerc); + NgProfiler::RegionTimer regc (timerc); - INDEX_3_CLOSED_HASHTABLE faceht(100); - openelements.SetSize(0); + INDEX_3_CLOSED_HASHTABLE faceht(100); + openelements.SetSize(0); - for (pi = PointIndex::BASE; pi < np+PointIndex::BASE; pi++) - if (selsonpoint[pi].Size()+elsonpoint[pi].Size()) - { - faceht.SetSize (2 * selsonpoint[pi].Size() + 4 * elsonpoint[pi].Size()); + for (PointIndex pi = PointIndex::BASE; pi < np+PointIndex::BASE; pi++) + if (selsonpoint[pi].Size()+elsonpoint[pi].Size()) + { + faceht.SetSize (2 * selsonpoint[pi].Size() + 4 * elsonpoint[pi].Size()); - FlatArray row = selsonpoint[pi]; - for (ii = 0; ii < row.Size(); ii++) + FlatArray row = selsonpoint[pi]; + for (ii = 0; ii < row.Size(); ii++) { - hel = SurfaceElement(row[ii]); - int ind = hel.GetIndex(); + hel = SurfaceElement(row[ii]); + int ind = hel.GetIndex(); - if (GetFaceDescriptor(ind).DomainIn() && + if (GetFaceDescriptor(ind).DomainIn() && (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) ) - { + { hel.NormalizeNumbering(); if (hel.PNum(1) == pi) - { - INDEX_3 i3(hel[0], hel[1], hel[2]); - INDEX_2 i2 (GetFaceDescriptor(ind).DomainIn(), - (hel.GetNP() == 3) - ? PointIndex (PointIndex::BASE-1) - : hel.PNum(4)); - faceht.Set (i3, i2); - } - } - if (GetFaceDescriptor(ind).DomainOut() && + { + INDEX_3 i3(hel[0], hel[1], hel[2]); + INDEX_2 i2 (GetFaceDescriptor(ind).DomainIn(), + (hel.GetNP() == 3) + ? PointIndex (PointIndex::BASE-1) + : hel.PNum(4)); + faceht.Set (i3, i2); + } + } + if (GetFaceDescriptor(ind).DomainOut() && (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) ) - { + { hel.Invert(); hel.NormalizeNumbering(); if (hel.PNum(1) == pi) - { - INDEX_3 i3(hel[0], hel[1], hel[2]); - INDEX_2 i2 (GetFaceDescriptor(ind).DomainOut(), - (hel.GetNP() == 3) - ? PointIndex (PointIndex::BASE-1) - : hel.PNum(4)); - faceht.Set (i3, i2); - } - } + { + INDEX_3 i3(hel[0], hel[1], hel[2]); + INDEX_2 i2 (GetFaceDescriptor(ind).DomainOut(), + (hel.GetNP() == 3) + ? PointIndex (PointIndex::BASE-1) + : hel.PNum(4)); + faceht.Set (i3, i2); + } + } } - FlatArray rowel = elsonpoint[pi]; - for (ii = 0; ii < rowel.Size(); ii++) + FlatArray rowel = elsonpoint[pi]; + for (ii = 0; ii < rowel.Size(); ii++) { - const Element & el = VolumeElement(rowel[ii]); + const Element & el = VolumeElement(rowel[ii]); - if (dom == 0 || el.GetIndex() == dom) - { - for (j = 1; j <= el.GetNFaces(); j++) - { - el.GetFace (j, hel); - hel.Invert(); - hel.NormalizeNumbering(); + if (dom == 0 || el.GetIndex() == dom) + { + for (int j = 1; j <= el.GetNFaces(); j++) + { + el.GetFace (j, hel); + hel.Invert(); + hel.NormalizeNumbering(); - if (hel[0] == pi) - { - INDEX_3 i3(hel[0], hel[1], hel[2]); - - if (faceht.Used (i3)) + if (hel[0] == pi) { - INDEX_2 i2 = faceht.Get(i3); - if (i2.I1() == el.GetIndex()) - { - i2.I1() = PointIndex::BASE-1; + INDEX_3 i3(hel[0], hel[1], hel[2]); + + if (faceht.Used (i3)) + { + INDEX_2 i2 = faceht.Get(i3); + if (i2.I1() == el.GetIndex()) + { + i2.I1() = PointIndex::BASE-1; + faceht.Set (i3, i2); + } + else + { + if (i2.I1() == 0) + { + PrintSysError ("more elements on face"); + (*testout) << "more elements on face!!!" << endl; + (*testout) << "el = " << el << endl; + (*testout) << "hel = " << hel << endl; + (*testout) << "face = " << i3 << endl; + (*testout) << "points = " << endl; + for (int jj = 1; jj <= 3; jj++) + (*testout) << "p = " << Point(i3.I(jj)) << endl; + } + } + } + else + { + hel.Invert(); + hel.NormalizeNumbering(); + INDEX_3 i3(hel[0], hel[1], hel[2]); + INDEX_2 i2(el.GetIndex(), + (hel.GetNP() == 3) + ? PointIndex (PointIndex::BASE-1) + : hel[3]); faceht.Set (i3, i2); - } - else - { - if (i2.I1() == 0) - { - PrintSysError ("more elements on face"); - (*testout) << "more elements on face!!!" << endl; - (*testout) << "el = " << el << endl; - (*testout) << "hel = " << hel << endl; - (*testout) << "face = " << i3 << endl; - (*testout) << "points = " << endl; - for (int jj = 1; jj <= 3; jj++) - (*testout) << "p = " << Point(i3.I(jj)) << endl; - } - } + } } - else - { - hel.Invert(); - hel.NormalizeNumbering(); - INDEX_3 i3(hel[0], hel[1], hel[2]); - INDEX_2 i2(el.GetIndex(), - (hel.GetNP() == 3) - ? PointIndex (PointIndex::BASE-1) - : hel[3]); - faceht.Set (i3, i2); - } - } - } - } + } + } } - for (i = 0; i < faceht.Size(); i++) - if (faceht.UsedPos (i)) - { - INDEX_3 i3; - INDEX_2 i2; - faceht.GetData (i, i3, i2); - if (i2.I1() != PointIndex::BASE-1) + for (int i = 0; i < faceht.Size(); i++) + if (faceht.UsedPos (i)) + { + INDEX_3 i3; + INDEX_2 i2; + faceht.GetData (i, i3, i2); + if (i2.I1() != PointIndex::BASE-1) { - Element2d tri; - tri.SetType ( (i2.I2() == PointIndex::BASE-1) ? TRIG : QUAD); - for (l = 0; l < 3; l++) - tri[l] = i3.I(l+1); - tri.PNum(4) = i2.I2(); - tri.SetIndex (i2.I1()); + Element2d tri; + tri.SetType ( (i2.I2() == PointIndex::BASE-1) ? TRIG : QUAD); + for (int l = 0; l < 3; l++) + tri[l] = i3.I(l+1); + tri.PNum(4) = i2.I2(); + tri.SetIndex (i2.I1()); - // tri.Invert(); + // tri.Invert(); - openelements.Append (tri); + openelements.Append (tri); } - } - } + } + } - int cnt3 = 0; - for (i = 0; i < openelements.Size(); i++) - if (openelements[i].GetNP() == 3) - cnt3++; + int cnt3 = 0; + for (i = 0; i < openelements.Size(); i++) + if (openelements[i].GetNP() == 3) + cnt3++; - int cnt4 = openelements.Size() - cnt3; + int cnt4 = openelements.Size() - cnt3; - MyStr treequad; - if (cnt4) - treequad = MyStr(" (") + MyStr(cnt3) + MyStr (" + ") + - MyStr(cnt4) + MyStr(")"); + MyStr treequad; + if (cnt4) + treequad = MyStr(" (") + MyStr(cnt3) + MyStr (" + ") + + MyStr(cnt4) + MyStr(")"); - PrintMessage (5, openelements.Size(), treequad, " open elements"); + PrintMessage (5, openelements.Size(), treequad, " open elements"); - BuildBoundaryEdges(); + BuildBoundaryEdges(); - NgProfiler::RegionTimer regd (timerd); + NgProfiler::RegionTimer regd (timerd); - for (i = 1; i <= openelements.Size(); i++) - { - const Element2d & sel = openelements.Get(i); + for (int i = 1; i <= openelements.Size(); i++) + { + const Element2d & sel = openelements.Get(i); - if (boundaryedges) - for (j = 1; j <= sel.GetNP(); j++) - { - INDEX_2 i2; - i2.I1() = sel.PNumMod(j); - i2.I2() = sel.PNumMod(j+1); - i2.Sort(); - boundaryedges->Set (i2, 1); - } + if (boundaryedges) + for (int j = 1; j <= sel.GetNP(); j++) + { + INDEX_2 i2; + i2.I1() = sel.PNumMod(j); + i2.I2() = sel.PNumMod(j+1); + i2.Sort(); + boundaryedges->Set (i2, 1); + } - for (j = 1; j <= 3; j++) - { - int pi = sel.PNum(j); - if (pi < points.Size()+PointIndex::BASE) - points[pi].SetType (FIXEDPOINT); - } - } + for (int j = 1; j <= 3; j++) + { + int pi = sel.PNum(j); + if (pi < points.Size()+PointIndex::BASE) + points[pi].SetType (FIXEDPOINT); + } + } - NgProfiler::RegionTimer rege (timere); + NgProfiler::RegionTimer rege (timere); - /* - for (i = 1; i <= GetNSeg(); i++) - { - const Segment & seg = LineSegment(i); - INDEX_2 i2(seg[0], seg[1]); - i2.Sort(); - - if (!boundaryedges->Used (i2)) - cerr << "WARNING: no boundedge, but seg edge: " << i2 << endl; - - boundaryedges -> Set (i2, 2); - segmentht -> Set (i2, i-1); - } - */ - } - - bool Mesh :: HasOpenQuads () const - { - int no = GetNOpenElements(); - for (int i = 0; i < no; i++) - if (openelements[i].GetNP() == 4) - return true; - return false; - } - - - - - - void Mesh :: FindOpenSegments (int surfnr) - { - int i, j, k; - - // new version, general elemetns - // hash index: pnum1-2 - // hash data : surfnr, surfel-nr (pos) or segment nr(neg) - INDEX_2_HASHTABLE faceht(4 * GetNSE()+GetNSeg()+1); - - PrintMessage (5, "Test Opensegments"); + /* for (i = 1; i <= GetNSeg(); i++) { - const Segment & seg = LineSegment (i); + const Segment & seg = LineSegment(i); + INDEX_2 i2(seg[0], seg[1]); + i2.Sort(); - if (surfnr == 0 || seg.si == surfnr) - { + if (!boundaryedges->Used (i2)) + cerr << "WARNING: no boundedge, but seg edge: " << i2 << endl; + + boundaryedges -> Set (i2, 2); + segmentht -> Set (i2, i-1); + } + */ + } + + bool Mesh :: HasOpenQuads () const + { + int no = GetNOpenElements(); + for (int i = 0; i < no; i++) + if (openelements[i].GetNP() == 4) + return true; + return false; + } + + + + + + void Mesh :: FindOpenSegments (int surfnr) + { + int i, j, k; + + // new version, general elemetns + // hash index: pnum1-2 + // hash data : surfnr, surfel-nr (pos) or segment nr(neg) + INDEX_2_HASHTABLE faceht(4 * GetNSE()+GetNSeg()+1); + + PrintMessage (5, "Test Opensegments"); + for (i = 1; i <= GetNSeg(); i++) + { + const Segment & seg = LineSegment (i); + + if (surfnr == 0 || seg.si == surfnr) + { INDEX_2 key(seg[0], seg[1]); INDEX_2 data(seg.si, -i); if (faceht.Used (key)) - { - cerr << "ERROR: Segment " << seg << " already used" << endl; - (*testout) << "ERROR: Segment " << seg << " already used" << endl; - } + { + cerr << "ERROR: Segment " << seg << " already used" << endl; + (*testout) << "ERROR: Segment " << seg << " already used" << endl; + } faceht.Set (key, data); - } + } } - for (i = 1; i <= GetNSeg(); i++) + for (i = 1; i <= GetNSeg(); i++) { - const Segment & seg = LineSegment (i); + const Segment & seg = LineSegment (i); - if (surfnr == 0 || seg.si == surfnr) - { + if (surfnr == 0 || seg.si == surfnr) + { INDEX_2 key(seg[1], seg[0]); if (!faceht.Used(key)) - { - cerr << "ERROR: Segment " << seg << " brother not used" << endl; - (*testout) << "ERROR: Segment " << seg << " brother not used" << endl; - } - } + { + cerr << "ERROR: Segment " << seg << " brother not used" << endl; + (*testout) << "ERROR: Segment " << seg << " brother not used" << endl; + } + } } - for (i = 1; i <= GetNSE(); i++) + for (i = 1; i <= GetNSE(); i++) { - const Element2d & el = SurfaceElement(i); - if (el.IsDeleted()) continue; + const Element2d & el = SurfaceElement(i); + if (el.IsDeleted()) continue; - if (surfnr == 0 || el.GetIndex() == surfnr) - { + if (surfnr == 0 || el.GetIndex() == surfnr) + { for (j = 1; j <= el.GetNP(); j++) - { - INDEX_2 seg (el.PNumMod(j), el.PNumMod(j+1)); - INDEX_2 data; + { + INDEX_2 seg (el.PNumMod(j), el.PNumMod(j+1)); + INDEX_2 data; - if (seg.I1() <= 0 || seg.I2() <= 0) + if (seg.I1() <= 0 || seg.I2() <= 0) cerr << "seg = " << seg << endl; - if (faceht.Used(seg)) - { - data = faceht.Get(seg); - if (data.I1() == el.GetIndex()) + if (faceht.Used(seg)) { - data.I1() = 0; - faceht.Set (seg, data); + data = faceht.Get(seg); + if (data.I1() == el.GetIndex()) + { + data.I1() = 0; + faceht.Set (seg, data); + } + else + { + PrintSysError ("hash table si not fitting for segment: ", + seg.I1(), "-", seg.I2(), " other = ", + data.I2()); + } } - else + else { - PrintSysError ("hash table si not fitting for segment: ", - seg.I1(), "-", seg.I2(), " other = ", - data.I2()); - } - } - else - { - Swap (seg.I1(), seg.I2()); - data.I1() = el.GetIndex(); - data.I2() = i; + Swap (seg.I1(), seg.I2()); + data.I1() = el.GetIndex(); + data.I2() = i; - faceht.Set (seg, data); - } - } - } + faceht.Set (seg, data); + } + } + } } - (*testout) << "open segments: " << endl; - opensegments.SetSize(0); - for (i = 1; i <= faceht.GetNBags(); i++) - for (j = 1; j <= faceht.GetBagSize(i); j++) - { - INDEX_2 i2; - INDEX_2 data; - faceht.GetData (i, j, i2, data); - if (data.I1()) // surfnr + (*testout) << "open segments: " << endl; + opensegments.SetSize(0); + for (i = 1; i <= faceht.GetNBags(); i++) + for (j = 1; j <= faceht.GetBagSize(i); j++) + { + INDEX_2 i2; + INDEX_2 data; + faceht.GetData (i, j, i2, data); + if (data.I1()) // surfnr { - Segment seg; - seg[0] = i2.I1(); - seg[1] = i2.I2(); - seg.si = data.I1(); + Segment seg; + seg[0] = i2.I1(); + seg[1] = i2.I2(); + seg.si = data.I1(); - // find geomdata: - if (data.I2() > 0) - { + // find geomdata: + if (data.I2() > 0) + { // segment due to triangle const Element2d & el = SurfaceElement (data.I2()); for (k = 1; k <= el.GetNP(); k++) - { - if (seg[0] == el.PNum(k)) + { + if (seg[0] == el.PNum(k)) seg.geominfo[0] = el.GeomInfoPi(k); - if (seg[1] == el.PNum(k)) + if (seg[1] == el.PNum(k)) seg.geominfo[1] = el.GeomInfoPi(k); - } + } (*testout) << "trig seg: "; - } - else - { + } + else + { // segment due to line const Segment & lseg = LineSegment (-data.I2()); seg.geominfo[0] = lseg.geominfo[0]; seg.geominfo[1] = lseg.geominfo[1]; (*testout) << "line seg: "; - } + } - (*testout) << seg[0] << " - " << seg[1] - << " len = " << Dist (Point(seg[0]), Point(seg[1])) - << endl; + (*testout) << seg[0] << " - " << seg[1] + << " len = " << Dist (Point(seg[0]), Point(seg[1])) + << endl; - opensegments.Append (seg); - if (seg.geominfo[0].trignum <= 0 || seg.geominfo[1].trignum <= 0) - { + opensegments.Append (seg); + if (seg.geominfo[0].trignum <= 0 || seg.geominfo[1].trignum <= 0) + { (*testout) << "Problem with open segment: " << seg << endl; - } + } } - } + } - PrintMessage (3, opensegments.Size(), " open segments found"); - (*testout) << opensegments.Size() << " open segments found" << endl; + PrintMessage (3, opensegments.Size(), " open segments found"); + (*testout) << opensegments.Size() << " open segments found" << endl; - /* - ptyps.SetSize (GetNP()); - for (i = 1; i <= ptyps.Size(); i++) - ptyps.Elem(i) = SURFACEPOINT; + /* + ptyps.SetSize (GetNP()); + for (i = 1; i <= ptyps.Size(); i++) + ptyps.Elem(i) = SURFACEPOINT; - for (i = 1; i <= GetNSeg(); i++) - { - const Segment & seg = LineSegment (i); - ptyps.Elem(seg[0]) = EDGEPOINT; - ptyps.Elem(seg[1]) = EDGEPOINT; - } - for (i = 1; i <= GetNOpenSegments(); i++) - { - const Segment & seg = GetOpenSegment (i); - ptyps.Elem(seg[0]) = EDGEPOINT; - ptyps.Elem(seg[1]) = EDGEPOINT; - } - */ - for (i = 1; i <= points.Size(); i++) - points.Elem(i).SetType(SURFACEPOINT); - - for (i = 1; i <= GetNSeg(); i++) - { - const Segment & seg = LineSegment (i); - points[seg[0]].SetType(EDGEPOINT); - points[seg[1]].SetType(EDGEPOINT); - } - for (i = 1; i <= GetNOpenSegments(); i++) - { - const Segment & seg = GetOpenSegment (i); - points[seg[0]].SetType (EDGEPOINT); - points[seg[1]].SetType (EDGEPOINT); - } - - - - /* - - for (i = 1; i <= openelements.Size(); i++) - { - const Element2d & sel = openelements.Get(i); - - if (boundaryedges) - for (j = 1; j <= sel.GetNP(); j++) - { - INDEX_2 i2; - i2.I1() = sel.PNumMod(j); - i2.I2() = sel.PNumMod(j+1); - i2.Sort(); - boundaryedges->Set (i2, 1); - } - - for (j = 1; j <= 3; j++) - { - int pi = sel.PNum(j); - if (pi <= ptyps.Size()) - ptyps.Elem(pi) = FIXEDPOINT; - } - } - */ - } - - - void Mesh :: RemoveOneLayerSurfaceElements () - { - int i, j; - int np = GetNP(); - - FindOpenSegments(); - BitArray frontpoints(np); - - frontpoints.Clear(); + for (i = 1; i <= GetNSeg(); i++) + { + const Segment & seg = LineSegment (i); + ptyps.Elem(seg[0]) = EDGEPOINT; + ptyps.Elem(seg[1]) = EDGEPOINT; + } for (i = 1; i <= GetNOpenSegments(); i++) { - const Segment & seg = GetOpenSegment(i); - frontpoints.Set (seg[0]); - frontpoints.Set (seg[1]); + const Segment & seg = GetOpenSegment (i); + ptyps.Elem(seg[0]) = EDGEPOINT; + ptyps.Elem(seg[1]) = EDGEPOINT; + } + */ + for (i = 1; i <= points.Size(); i++) + points.Elem(i).SetType(SURFACEPOINT); + + for (i = 1; i <= GetNSeg(); i++) + { + const Segment & seg = LineSegment (i); + points[seg[0]].SetType(EDGEPOINT); + points[seg[1]].SetType(EDGEPOINT); + } + for (i = 1; i <= GetNOpenSegments(); i++) + { + const Segment & seg = GetOpenSegment (i); + points[seg[0]].SetType (EDGEPOINT); + points[seg[1]].SetType (EDGEPOINT); } - for (i = 1; i <= GetNSE(); i++) + + + /* + + for (i = 1; i <= openelements.Size(); i++) + { + const Element2d & sel = openelements.Get(i); + + if (boundaryedges) + for (j = 1; j <= sel.GetNP(); j++) + { + INDEX_2 i2; + i2.I1() = sel.PNumMod(j); + i2.I2() = sel.PNumMod(j+1); + i2.Sort(); + boundaryedges->Set (i2, 1); + } + + for (j = 1; j <= 3; j++) + { + int pi = sel.PNum(j); + if (pi <= ptyps.Size()) + ptyps.Elem(pi) = FIXEDPOINT; + } + } + */ + } + + + void Mesh :: RemoveOneLayerSurfaceElements () + { + int i, j; + int np = GetNP(); + + FindOpenSegments(); + BitArray frontpoints(np); + + frontpoints.Clear(); + for (i = 1; i <= GetNOpenSegments(); i++) { - Element2d & sel = surfelements.Elem(i); - int remove = 0; - for (j = 1; j <= sel.GetNP(); j++) - if (frontpoints.Test(sel.PNum(j))) - remove = 1; - if (remove) - sel.PNum(1) = 0; + const Segment & seg = GetOpenSegment(i); + frontpoints.Set (seg[0]); + frontpoints.Set (seg[1]); } - for (i = surfelements.Size(); i >= 1; i--) + for (i = 1; i <= GetNSE(); i++) { - if (surfelements.Elem(i).PNum(1) == 0) - { + Element2d & sel = surfelements.Elem(i); + int remove = 0; + for (j = 1; j <= sel.GetNP(); j++) + if (frontpoints.Test(sel.PNum(j))) + remove = 1; + if (remove) + sel.PNum(1) = 0; + } + + for (i = surfelements.Size(); i >= 1; i--) + { + if (surfelements.Elem(i).PNum(1) == 0) + { surfelements.Elem(i) = surfelements.Last(); surfelements.DeleteLast(); - } + } } - for (int i = 0; i < facedecoding.Size(); i++) - facedecoding[i].firstelement = -1; - for (int i = surfelements.Size()-1; i >= 0; i--) + for (int i = 0; i < facedecoding.Size(); i++) + facedecoding[i].firstelement = -1; + for (int i = surfelements.Size()-1; i >= 0; i--) { - int ind = surfelements[i].GetIndex(); - surfelements[i].next = facedecoding[ind-1].firstelement; - facedecoding[ind-1].firstelement = i; + int ind = surfelements[i].GetIndex(); + surfelements[i].next = facedecoding[ind-1].firstelement; + facedecoding[ind-1].firstelement = i; } - timestamp = NextTimeStamp(); - // Compress(); - } + timestamp = NextTimeStamp(); + // Compress(); + } - void Mesh :: FreeOpenElementsEnvironment (int layers) - { - int i, j, k; - PointIndex pi; - const int large = 9999; - Array dist(GetNP()); + void Mesh :: FreeOpenElementsEnvironment (int layers) + { + int i, j, k; + PointIndex pi; + const int large = 9999; + Array dist(GetNP()); - dist = large; + dist = large; - for (int i = 1; i <= GetNOpenElements(); i++) + for (int i = 1; i <= GetNOpenElements(); i++) { - const Element2d & face = OpenElement(i); - for (j = 0; j < face.GetNP(); j++) - dist[face[j]] = 1; + const Element2d & face = OpenElement(i); + for (j = 0; j < face.GetNP(); j++) + dist[face[j]] = 1; } - for (k = 1; k <= layers; k++) - for (i = 1; i <= GetNE(); i++) - { - const Element & el = VolumeElement(i); - if (el[0] == -1 || el.IsDeleted()) continue; + for (k = 1; k <= layers; k++) + for (i = 1; i <= GetNE(); i++) + { + const Element & el = VolumeElement(i); + if (el[0] == -1 || el.IsDeleted()) continue; - int elmin = large; - for (j = 0; j < el.GetNP(); j++) - if (dist[el[j]] < elmin) - elmin = dist[el[j]]; + int elmin = large; + for (j = 0; j < el.GetNP(); j++) + if (dist[el[j]] < elmin) + elmin = dist[el[j]]; - if (elmin < large) + if (elmin < large) { - for (j = 0; j < el.GetNP(); j++) - if (dist[el[j]] > elmin+1) - dist[el[j]] = elmin+1; + for (j = 0; j < el.GetNP(); j++) + if (dist[el[j]] > elmin+1) + dist[el[j]] = elmin+1; } - } + } - int cntfree = 0; - for (i = 1; i <= GetNE(); i++) - { - Element & el = VolumeElement(i); - if (el[0] == -1 || el.IsDeleted()) continue; - - int elmin = large; - for (j = 0; j < el.GetNP(); j++) - if (dist[el[j]] < elmin) - elmin = dist[el[j]]; - - el.flags.fixed = elmin > layers; - // eltyps.Elem(i) = (elmin <= layers) ? - // FREEELEMENT : FIXEDELEMENT; - if (elmin <= layers) - cntfree++; - } - - PrintMessage (5, "free: ", cntfree, ", fixed: ", GetNE()-cntfree); - (*testout) << "free: " << cntfree << ", fixed: " << GetNE()-cntfree << endl; - - for (pi = PointIndex::BASE; - pi < GetNP()+PointIndex::BASE; pi++) - { - if (dist[pi] > layers+1) - points[pi].SetType(FIXEDPOINT); - } - } - - - - void Mesh :: SetLocalH (const Point3d & pmin, const Point3d & pmax, double grading) - { - Point3d c = Center (pmin, pmax); - double d = max3 (pmax.X()-pmin.X(), - pmax.Y()-pmin.Y(), - pmax.Z()-pmin.Z()); - d /= 2; - Point3d pmin2 = c - Vec3d (d, d, d); - Point3d pmax2 = c + Vec3d (d, d, d); - - - delete lochfunc; - lochfunc = new LocalH (pmin2, pmax2, grading); - } - - void Mesh :: RestrictLocalH (const Point3d & p, double hloc) - { - if(hloc < hmin) - hloc = hmin; - - //cout << "restrict h in " << p << " to " << hloc << endl; - if (!lochfunc) + int cntfree = 0; + for (i = 1; i <= GetNE(); i++) { - PrintWarning("RestrictLocalH called, creating mesh-size tree"); + Element & el = VolumeElement(i); + if (el[0] == -1 || el.IsDeleted()) continue; - Point3d boxmin, boxmax; - GetBox (boxmin, boxmax); - SetLocalH (boxmin, boxmax, 0.8); + int elmin = large; + for (j = 0; j < el.GetNP(); j++) + if (dist[el[j]] < elmin) + elmin = dist[el[j]]; + + el.flags.fixed = elmin > layers; + // eltyps.Elem(i) = (elmin <= layers) ? + // FREEELEMENT : FIXEDELEMENT; + if (elmin <= layers) + cntfree++; } - lochfunc -> SetH (p, hloc); - } + PrintMessage (5, "free: ", cntfree, ", fixed: ", GetNE()-cntfree); + (*testout) << "free: " << cntfree << ", fixed: " << GetNE()-cntfree << endl; - void Mesh :: RestrictLocalHLine (const Point3d & p1, - const Point3d & p2, - double hloc) - { - if(hloc < hmin) - hloc = hmin; - - // cout << "restrict h along " << p1 << " - " << p2 << " to " << hloc << endl; - int i; - int steps = int (Dist (p1, p2) / hloc) + 2; - Vec3d v(p1, p2); - - for (i = 0; i <= steps; i++) + for (pi = PointIndex::BASE; + pi < GetNP()+PointIndex::BASE; pi++) { - Point3d p = p1 + (double(i)/double(steps) * v); - RestrictLocalH (p, hloc); + if (dist[pi] > layers+1) + points[pi].SetType(FIXEDPOINT); } - } + } - void Mesh :: SetMinimalH (double h) - { - hmin = h; - } + + void Mesh :: SetLocalH (const Point3d & pmin, const Point3d & pmax, double grading) + { + Point3d c = Center (pmin, pmax); + double d = max3 (pmax.X()-pmin.X(), + pmax.Y()-pmin.Y(), + pmax.Z()-pmin.Z()); + d /= 2; + Point3d pmin2 = c - Vec3d (d, d, d); + Point3d pmax2 = c + Vec3d (d, d, d); - void Mesh :: SetGlobalH (double h) - { - hglob = h; - } + delete lochfunc; + lochfunc = new LocalH (pmin2, pmax2, grading); + } - double Mesh :: MaxHDomain (int dom) const - { - if (maxhdomain.Size()) - return maxhdomain.Get(dom); - else - return 1e10; - } + void Mesh :: RestrictLocalH (const Point3d & p, double hloc) + { + if(hloc < hmin) + hloc = hmin; - void Mesh :: SetMaxHDomain (const Array & mhd) - { - maxhdomain.SetSize(mhd.Size()); - for (int i = 1; i <= mhd.Size(); i++) - maxhdomain.Elem(i) = mhd.Get(i); - } - - - double Mesh :: GetH (const Point3d & p) const - { - double hmin = hglob; - if (lochfunc) + //cout << "restrict h in " << p << " to " << hloc << endl; + if (!lochfunc) { - double hl = lochfunc->GetH (p); - if (hl < hglob) - hmin = hl; + PrintWarning("RestrictLocalH called, creating mesh-size tree"); + + Point3d boxmin, boxmax; + GetBox (boxmin, boxmax); + SetLocalH (boxmin, boxmax, 0.8); } - return hmin; - } - double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax) - { - double hmin = hglob; - if (lochfunc) + lochfunc -> SetH (p, hloc); + } + + void Mesh :: RestrictLocalHLine (const Point3d & p1, + const Point3d & p2, + double hloc) + { + if(hloc < hmin) + hloc = hmin; + + // cout << "restrict h along " << p1 << " - " << p2 << " to " << hloc << endl; + int i; + int steps = int (Dist (p1, p2) / hloc) + 2; + Vec3d v(p1, p2); + + for (i = 0; i <= steps; i++) { - double hl = lochfunc->GetMinH (pmin, pmax); - if (hl < hmin) - hmin = hl; + Point3d p = p1 + (double(i)/double(steps) * v); + RestrictLocalH (p, hloc); } - return hmin; - } + } + void Mesh :: SetMinimalH (double h) + { + hmin = h; + } + void Mesh :: SetGlobalH (double h) + { + hglob = h; + } - double Mesh :: AverageH (int surfnr) const - { - int i, j, n; - double hi, hsum; - double maxh = 0, minh = 1e10; + double Mesh :: MaxHDomain (int dom) const + { + if (maxhdomain.Size()) + return maxhdomain.Get(dom); + else + return 1e10; + } - hsum = 0; - n = 0; - for (i = 1; i <= GetNSE(); i++) + void Mesh :: SetMaxHDomain (const Array & mhd) + { + maxhdomain.SetSize(mhd.Size()); + for (int i = 1; i <= mhd.Size(); i++) + maxhdomain.Elem(i) = mhd.Get(i); + } + + + double Mesh :: GetH (const Point3d & p) const + { + double hmin = hglob; + if (lochfunc) { - const Element2d & el = SurfaceElement(i); - if (surfnr == 0 || el.GetIndex() == surfnr) - { + double hl = lochfunc->GetH (p); + if (hl < hglob) + hmin = hl; + } + return hmin; + } + + double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax) + { + double hmin = hglob; + if (lochfunc) + { + double hl = lochfunc->GetMinH (pmin, pmax); + if (hl < hmin) + hmin = hl; + } + return hmin; + } + + + + + + double Mesh :: AverageH (int surfnr) const + { + int i, j, n; + double hi, hsum; + double maxh = 0, minh = 1e10; + + hsum = 0; + n = 0; + for (i = 1; i <= GetNSE(); i++) + { + const Element2d & el = SurfaceElement(i); + if (surfnr == 0 || el.GetIndex() == surfnr) + { for (j = 1; j <= 3; j++) - { - hi = Dist (Point (el.PNumMod(j)), - Point (el.PNumMod(j+1))); + { + hi = Dist (Point (el.PNumMod(j)), + Point (el.PNumMod(j+1))); - hsum += hi; + hsum += hi; - if (hi > maxh) maxh = hi; - if (hi < minh) minh = hi; - n++; - } - } + if (hi > maxh) maxh = hi; + if (hi < minh) minh = hi; + n++; + } + } } - PrintMessage (5, "minh = ", minh, " avh = ", (hsum/n), " maxh = ", maxh); - return (hsum / n); - } + PrintMessage (5, "minh = ", minh, " avh = ", (hsum/n), " maxh = ", maxh); + return (hsum / n); + } - void Mesh :: CalcLocalH () - { - if (!lochfunc) + void Mesh :: CalcLocalH () + { + if (!lochfunc) { - Point3d pmin, pmax; - GetBox (pmin, pmax); - SetLocalH (pmin, pmax, mparam.grading); + Point3d pmin, pmax; + GetBox (pmin, pmax); + SetLocalH (pmin, pmax, mparam.grading); } - PrintMessage (3, - "CalcLocalH: ", - GetNP(), " Points ", - GetNE(), " Elements ", - GetNSE(), " Surface Elements"); + PrintMessage (3, + "CalcLocalH: ", + GetNP(), " Points ", + GetNE(), " Elements ", + GetNSE(), " Surface Elements"); - for (int i = 0; i < GetNSE(); i++) + for (int i = 0; i < GetNSE(); i++) { - const Element2d & el = surfelements[i]; - int j; + const Element2d & el = surfelements[i]; + int j; - if (el.GetNP() == 3) - { + if (el.GetNP() == 3) + { double hel = -1; for (j = 1; j <= 3; j++) - { - const Point3d & p1 = points[el.PNumMod(j)]; - const Point3d & p2 = points[el.PNumMod(j+1)]; + { + const Point3d & p1 = points[el.PNumMod(j)]; + const Point3d & p2 = points[el.PNumMod(j+1)]; - /* - INDEX_2 i21(el.PNumMod(j), el.PNumMod(j+1)); - INDEX_2 i22(el.PNumMod(j+1), el.PNumMod(j)); - if (! identifiedpoints->Used (i21) && - ! identifiedpoints->Used (i22) ) - */ - if (!ident -> UsedSymmetric (el.PNumMod(j), - el.PNumMod(j+1))) - { - double hedge = Dist (p1, p2); - if (hedge > hel) - hel = hedge; - // lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); - // (*testout) << "trigseth, p1,2 = " << el.PNumMod(j) << ", " << el.PNumMod(j+1) - // << " h = " << (2 * Dist(p1, p2)) << endl; - } - } + /* + INDEX_2 i21(el.PNumMod(j), el.PNumMod(j+1)); + INDEX_2 i22(el.PNumMod(j+1), el.PNumMod(j)); + if (! identifiedpoints->Used (i21) && + ! identifiedpoints->Used (i22) ) + */ + if (!ident -> UsedSymmetric (el.PNumMod(j), + el.PNumMod(j+1))) + { + double hedge = Dist (p1, p2); + if (hedge > hel) + hel = hedge; + // lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); + // (*testout) << "trigseth, p1,2 = " << el.PNumMod(j) << ", " << el.PNumMod(j+1) + // << " h = " << (2 * Dist(p1, p2)) << endl; + } + } if (hel > 0) + { + const Point3d & p1 = points[el.PNum(1)]; + const Point3d & p2 = points[el.PNum(2)]; + const Point3d & p3 = points[el.PNum(3)]; + lochfunc->SetH (Center (p1, p2, p3), hel); + } + } + else + { { - const Point3d & p1 = points[el.PNum(1)]; - const Point3d & p2 = points[el.PNum(2)]; - const Point3d & p3 = points[el.PNum(3)]; - lochfunc->SetH (Center (p1, p2, p3), hel); - } - } - else - { - { - const Point3d & p1 = points[el.PNum(1)]; - const Point3d & p2 = points[el.PNum(2)]; - lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); + const Point3d & p1 = points[el.PNum(1)]; + const Point3d & p2 = points[el.PNum(2)]; + lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); } { - const Point3d & p1 = points[el.PNum(3)]; - const Point3d & p2 = points[el.PNum(4)]; - lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); + const Point3d & p1 = points[el.PNum(3)]; + const Point3d & p2 = points[el.PNum(4)]; + lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); } - } + } } - for (int i = 0; i < GetNSeg(); i++) + for (int i = 0; i < GetNSeg(); i++) { - const Segment & seg = segments[i]; - const Point3d & p1 = points[seg[0]]; - const Point3d & p2 = points[seg[1]]; - /* - 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[0], seg[1])) - { + const Segment & seg = segments[i]; + const Point3d & p1 = points[seg[0]]; + const Point3d & p2 = points[seg[1]]; + /* + 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[0], seg[1])) + { lochfunc->SetH (Center (p1, p2), Dist (p1, p2)); - } + } } - /* + /* cerr << "do vol" << endl; for (i = 1; i <= GetNE(); i++) { @@ -2566,9 +2565,9 @@ namespace netgen } } } - */ + */ - /* + /* const char * meshsizefilename = globflags.GetStringFlag ("meshsize", NULL); if (meshsizefilename) @@ -2588,421 +2587,421 @@ namespace netgen } } } - */ - // lochfunc -> Convexify(); - // lochfunc -> PrintMemInfo (cout); - } + */ + // lochfunc -> Convexify(); + // lochfunc -> PrintMemInfo (cout); + } - void Mesh :: CalcLocalHFromPointDistances(void) - { - PrintMessage (3, "Calculating local h from point distances"); + void Mesh :: CalcLocalHFromPointDistances(void) + { + PrintMessage (3, "Calculating local h from point distances"); - if (!lochfunc) + if (!lochfunc) { - Point3d pmin, pmax; - GetBox (pmin, pmax); + Point3d pmin, pmax; + GetBox (pmin, pmax); - SetLocalH (pmin, pmax, mparam.grading); + SetLocalH (pmin, pmax, mparam.grading); } - PointIndex i,j; - double hl; + PointIndex i,j; + double hl; - for (i = PointIndex::BASE; + for (i = PointIndex::BASE; i < GetNP()+PointIndex::BASE; i++) { - for(j=i+1; j edges(3 * GetNP() + 2); - INDEX_2_HASHTABLE bedges(GetNSeg() + 2); - int i, j; + INDEX_2_HASHTABLE edges(3 * GetNP() + 2); + INDEX_2_HASHTABLE bedges(GetNSeg() + 2); + int i, j; - for (i = 1; i <= GetNSeg(); i++) + for (i = 1; i <= GetNSeg(); i++) { - const Segment & seg = LineSegment(i); - INDEX_2 i2(seg[0], seg[1]); - i2.Sort(); - bedges.Set (i2, 1); + const Segment & seg = LineSegment(i); + INDEX_2 i2(seg[0], seg[1]); + i2.Sort(); + bedges.Set (i2, 1); } - for (i = 1; i <= GetNSE(); i++) + for (i = 1; i <= GetNSE(); i++) { - const Element2d & sel = SurfaceElement(i); - if (!sel.PNum(1)) - continue; - for (j = 1; j <= 3; j++) - { + const Element2d & sel = SurfaceElement(i); + if (!sel.PNum(1)) + continue; + for (j = 1; j <= 3; j++) + { INDEX_2 i2(sel.PNumMod(j), sel.PNumMod(j+1)); i2.Sort(); if (bedges.Used(i2)) continue; if (edges.Used(i2)) - { - int other = edges.Get(i2); + { + int other = edges.Get(i2); - const Element2d & elother = SurfaceElement(other); + const Element2d & elother = SurfaceElement(other); - int pi3 = 1; - while ( (sel.PNum(pi3) == i2.I1()) || - (sel.PNum(pi3) == i2.I2())) + int pi3 = 1; + while ( (sel.PNum(pi3) == i2.I1()) || + (sel.PNum(pi3) == i2.I2())) pi3++; - pi3 = sel.PNum(pi3); + pi3 = sel.PNum(pi3); - int pi4 = 1; - while ( (elother.PNum(pi4) == i2.I1()) || - (elother.PNum(pi4) == i2.I2())) + int pi4 = 1; + while ( (elother.PNum(pi4) == i2.I1()) || + (elother.PNum(pi4) == i2.I2())) pi4++; - pi4 = elother.PNum(pi4); + pi4 = elother.PNum(pi4); - double rad = ComputeCylinderRadius (Point (i2.I1()), - Point (i2.I2()), - Point (pi3), - Point (pi4)); + double rad = ComputeCylinderRadius (Point (i2.I1()), + Point (i2.I2()), + Point (pi3), + Point (pi4)); - RestrictLocalHLine (Point(i2.I1()), Point(i2.I2()), rad/elperr); + RestrictLocalHLine (Point(i2.I1()), Point(i2.I2()), rad/elperr); - /* - (*testout) << "pi1,2, 3, 4 = " << i2.I1() << ", " << i2.I2() << ", " << pi3 << ", " << pi4 - << " p1 = " << Point(i2.I1()) - << ", p2 = " << Point(i2.I2()) - // << ", p3 = " << Point(pi3) - // << ", p4 = " << Point(pi4) - << ", rad = " << rad << endl; - */ - } + /* + (*testout) << "pi1,2, 3, 4 = " << i2.I1() << ", " << i2.I2() << ", " << pi3 << ", " << pi4 + << " p1 = " << Point(i2.I1()) + << ", p2 = " << Point(i2.I2()) + // << ", p3 = " << Point(pi3) + // << ", p4 = " << Point(pi4) + << ", rad = " << rad << endl; + */ + } else - edges.Set (i2, i); - } + edges.Set (i2, i); + } } - // Restrict h due to line segments + // Restrict h due to line segments - for (i = 1; i <= GetNSeg(); i++) + for (i = 1; i <= GetNSeg(); i++) { - const Segment & seg = LineSegment(i); - const Point3d & p1 = Point(seg[0]); - const Point3d & p2 = Point(seg[1]); - RestrictLocalH (Center (p1, p2), Dist (p1, p2)); + const Segment & seg = LineSegment(i); + const Point3d & p1 = Point(seg[0]); + const Point3d & p2 = Point(seg[1]); + RestrictLocalH (Center (p1, p2), Dist (p1, p2)); } - /* + /* - int i, j; - int np = GetNP(); - int nseg = GetNSeg(); - int nse = GetNSE(); + int i, j; + int np = GetNP(); + int nseg = GetNSeg(); + int nse = GetNSE(); - Array normals(np); - BitArray linepoint(np); + Array normals(np); + BitArray linepoint(np); - linepoint.Clear(); - for (i = 1; i <= nseg; i++) - { - linepoint.Set (LineSegment(i)[0]); - linepoint.Set (LineSegment(i)[1]); - } + linepoint.Clear(); + for (i = 1; i <= nseg; i++) + { + linepoint.Set (LineSegment(i)[0]); + linepoint.Set (LineSegment(i)[1]); + } - for (i = 1; i <= np; i++) - normals.Elem(i) = Vec3d(0,0,0); + for (i = 1; i <= np; i++) + normals.Elem(i) = Vec3d(0,0,0); - for (i = 1; i <= nse; i++) - { - Element2d & el = SurfaceElement(i); - Vec3d nf = Cross (Vec3d (Point (el.PNum(1)), Point(el.PNum(2))), - Vec3d (Point (el.PNum(1)), Point(el.PNum(3)))); - for (j = 1; j <= 3; j++) - normals.Elem(el.PNum(j)) += nf; - } + for (i = 1; i <= nse; i++) + { + Element2d & el = SurfaceElement(i); + Vec3d nf = Cross (Vec3d (Point (el.PNum(1)), Point(el.PNum(2))), + Vec3d (Point (el.PNum(1)), Point(el.PNum(3)))); + for (j = 1; j <= 3; j++) + normals.Elem(el.PNum(j)) += nf; + } - for (i = 1; i <= np; i++) - normals.Elem(i) /= (1e-12 + normals.Elem(i).Length()); + for (i = 1; i <= np; i++) + normals.Elem(i) /= (1e-12 + normals.Elem(i).Length()); - for (i = 1; i <= nse; i++) - { - Element2d & el = SurfaceElement(i); - Vec3d nf = Cross (Vec3d (Point (el.PNum(1)), Point(el.PNum(2))), - Vec3d (Point (el.PNum(1)), Point(el.PNum(3)))); - nf /= nf.Length(); - Point3d c = Center (Point(el.PNum(1)), - Point(el.PNum(2)), - Point(el.PNum(3))); + for (i = 1; i <= nse; i++) + { + Element2d & el = SurfaceElement(i); + Vec3d nf = Cross (Vec3d (Point (el.PNum(1)), Point(el.PNum(2))), + Vec3d (Point (el.PNum(1)), Point(el.PNum(3)))); + nf /= nf.Length(); + Point3d c = Center (Point(el.PNum(1)), + Point(el.PNum(2)), + Point(el.PNum(3))); - for (j = 1; j <= 3; j++) - { - if (!linepoint.Test (el.PNum(j))) - { - double dist = Dist (c, Point(el.PNum(j))); - double dn = (nf - normals.Get(el.PNum(j))).Length(); + for (j = 1; j <= 3; j++) + { + if (!linepoint.Test (el.PNum(j))) + { + double dist = Dist (c, Point(el.PNum(j))); + double dn = (nf - normals.Get(el.PNum(j))).Length(); - RestrictLocalH (Point(el.PNum(j)), dist / (dn+1e-12) /elperr); - } - } - } - */ - } + RestrictLocalH (Point(el.PNum(j)), dist / (dn+1e-12) /elperr); + } + } + } + */ + } - void Mesh :: RestrictLocalH (resthtype rht, int nr, double loch) - { - int i; - switch (rht) + void Mesh :: RestrictLocalH (resthtype rht, int nr, double loch) + { + int i; + switch (rht) { case RESTRICTH_FACE: - { - for (i = 1; i <= GetNSE(); i++) + { + for (i = 1; i <= GetNSE(); i++) { - const Element2d & sel = SurfaceElement(i); - if (sel.GetIndex() == nr) - RestrictLocalH (RESTRICTH_SURFACEELEMENT, i, loch); + const Element2d & sel = SurfaceElement(i); + if (sel.GetIndex() == nr) + RestrictLocalH (RESTRICTH_SURFACEELEMENT, i, loch); } - break; - } + break; + } case RESTRICTH_EDGE: - { - for (i = 1; i <= GetNSeg(); i++) + { + for (i = 1; i <= GetNSeg(); i++) { - const Segment & seg = LineSegment(i); - if (seg.edgenr == nr) - RestrictLocalH (RESTRICTH_SEGMENT, i, loch); + const Segment & seg = LineSegment(i); + if (seg.edgenr == nr) + RestrictLocalH (RESTRICTH_SEGMENT, i, loch); } - break; - } + break; + } case RESTRICTH_POINT: - { - RestrictLocalH (Point (nr), loch); - break; - } + { + RestrictLocalH (Point (nr), loch); + break; + } case RESTRICTH_SURFACEELEMENT: - { - const Element2d & sel = SurfaceElement(nr); - Point3d p = Center (Point(sel.PNum(1)), - Point(sel.PNum(2)), - Point(sel.PNum(3))); - RestrictLocalH (p, loch); - break; - } + { + const Element2d & sel = SurfaceElement(nr); + Point3d p = Center (Point(sel.PNum(1)), + Point(sel.PNum(2)), + Point(sel.PNum(3))); + RestrictLocalH (p, loch); + break; + } case RESTRICTH_SEGMENT: - { - const Segment & seg = LineSegment(nr); - RestrictLocalHLine (Point (seg[0]), Point(seg[1]), loch); - break; - } + { + const Segment & seg = LineSegment(nr); + RestrictLocalHLine (Point (seg[0]), Point(seg[1]), loch); + break; + } } - } + } - void Mesh :: LoadLocalMeshSize (const char * meshsizefilename) - { - // Philippose - 10/03/2009 - // Improve error checking when loading and reading - // the local mesh size file + void Mesh :: LoadLocalMeshSize (const char * meshsizefilename) + { + // Philippose - 10/03/2009 + // Improve error checking when loading and reading + // the local mesh size file - if (!meshsizefilename) return; + if (!meshsizefilename) return; - ifstream msf(meshsizefilename); + ifstream msf(meshsizefilename); - // Philippose - 09/03/2009 - // Adding print message information in case the specified - // does not exist, or does not load successfully due to - // other reasons such as access rights, etc... - if (!msf) + // Philippose - 09/03/2009 + // Adding print message information in case the specified + // does not exist, or does not load successfully due to + // other reasons such as access rights, etc... + if (!msf) { - PrintMessage(3, "Error loading mesh size file: ", meshsizefilename, "....","Skipping!"); - return; + PrintMessage(3, "Error loading mesh size file: ", meshsizefilename, "....","Skipping!"); + return; } - PrintMessage (3, "Load local mesh-size file: ", meshsizefilename); + PrintMessage (3, "Load local mesh-size file: ", meshsizefilename); - int nmsp = 0; - int nmsl = 0; + int nmsp = 0; + int nmsl = 0; - msf >> nmsp; - if(nmsp > 0) + msf >> nmsp; + if(nmsp > 0) { - if(!msf.good()) - throw NgException ("Mesh-size file error: No points found\n"); - PrintMessage (4, "Number of mesh-size restriction points: ", nmsp); + if(!msf.good()) + throw NgException ("Mesh-size file error: No points found\n"); + PrintMessage (4, "Number of mesh-size restriction points: ", nmsp); } - else + else { - msf.close(); - return; + msf.close(); + return; } - for (int i = 0; i < nmsp; i++) + for (int i = 0; i < nmsp; i++) { - Point3d pi; - double hi; - msf >> pi.X() >> pi.Y() >> pi.Z(); - msf >> hi; - if (!msf.good()) - throw NgException ("Mesh-size file error: Number of points don't match specified list size\n"); - RestrictLocalH (pi, hi); + Point3d pi; + double hi; + msf >> pi.X() >> pi.Y() >> pi.Z(); + msf >> hi; + if (!msf.good()) + throw NgException ("Mesh-size file error: Number of points don't match specified list size\n"); + RestrictLocalH (pi, hi); } - msf >> nmsl; - if(nmsl > 0) + msf >> nmsl; + if(nmsl > 0) { - cout << "Number of line definitions expected = " << nmsl << endl; - if(!msf.good()) - throw NgException ("Mesh-size file error: No line definitions found\n"); - PrintMessage (4, "Number of mesh-size restriction lines: ", nmsl); + cout << "Number of line definitions expected = " << nmsl << endl; + if(!msf.good()) + throw NgException ("Mesh-size file error: No line definitions found\n"); + PrintMessage (4, "Number of mesh-size restriction lines: ", nmsl); } - else + else { - msf.close(); - return; + msf.close(); + return; } - for (int i = 0; i < nmsl; i++) + for (int i = 0; i < nmsl; i++) { - Point3d p1, p2; - double hi; - msf >> p1.X() >> p1.Y() >> p1.Z(); - msf >> p2.X() >> p2.Y() >> p2.Z(); - msf >> hi; - if (!msf.good()) - throw NgException ("Mesh-size file error: Number of line definitions don't match specified list size\n"); - RestrictLocalHLine (p1, p2, hi); + Point3d p1, p2; + double hi; + msf >> p1.X() >> p1.Y() >> p1.Z(); + msf >> p2.X() >> p2.Y() >> p2.Z(); + msf >> hi; + if (!msf.good()) + throw NgException ("Mesh-size file error: Number of line definitions don't match specified list size\n"); + RestrictLocalHLine (p1, p2, hi); } - msf.close(); - } + msf.close(); + } - void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, int dom) const - { - if (points.Size() == 0) + void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, int dom) const + { + if (points.Size() == 0) { - pmin = pmax = Point3d(0,0,0); - return; + pmin = pmax = Point3d(0,0,0); + return; } - if (dom <= 0) + if (dom <= 0) { - pmin = Point3d (1e10, 1e10, 1e10); - pmax = Point3d (-1e10, -1e10, -1e10); + pmin = Point3d (1e10, 1e10, 1e10); + pmax = Point3d (-1e10, -1e10, -1e10); - for (PointIndex pi = PointIndex::BASE; - pi < GetNP()+PointIndex::BASE; pi++) - { + for (PointIndex pi = PointIndex::BASE; + pi < GetNP()+PointIndex::BASE; pi++) + { pmin.SetToMin ( (*this) [pi] ); pmax.SetToMax ( (*this) [pi] ); - } + } } - else + else { - int j, nse = GetNSE(); - SurfaceElementIndex sei; + int j, nse = GetNSE(); + SurfaceElementIndex sei; - pmin = Point3d (1e10, 1e10, 1e10); - pmax = Point3d (-1e10, -1e10, -1e10); - for (sei = 0; sei < nse; sei++) - { + pmin = Point3d (1e10, 1e10, 1e10); + pmax = Point3d (-1e10, -1e10, -1e10); + for (sei = 0; sei < nse; sei++) + { const Element2d & el = (*this)[sei]; if (el.IsDeleted() ) continue; if (dom == -1 || el.GetIndex() == dom) - { - for (j = 0; j < 3; j++) - { - pmin.SetToMin ( (*this) [el[j]] ); - pmax.SetToMax ( (*this) [el[j]] ); - } - } - } + { + for (j = 0; j < 3; j++) + { + pmin.SetToMin ( (*this) [el[j]] ); + pmax.SetToMax ( (*this) [el[j]] ); + } + } + } } - if (pmin.X() > 0.5e10) + if (pmin.X() > 0.5e10) { - pmin = pmax = Point3d(0,0,0); + pmin = pmax = Point3d(0,0,0); } - } + } - void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, POINTTYPE ptyp) const - { - if (points.Size() == 0) + void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, POINTTYPE ptyp) const + { + if (points.Size() == 0) { - pmin = pmax = Point3d(0,0,0); - return; + pmin = pmax = Point3d(0,0,0); + return; } - pmin = Point3d (1e10, 1e10, 1e10); - pmax = Point3d (-1e10, -1e10, -1e10); + pmin = Point3d (1e10, 1e10, 1e10); + pmax = Point3d (-1e10, -1e10, -1e10); - for (PointIndex pi = PointIndex::BASE; + for (PointIndex pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) - if (points[pi].Type() <= ptyp) - { - pmin.SetToMin ( (*this) [pi] ); - pmax.SetToMax ( (*this) [pi] ); - } - } + if (points[pi].Type() <= ptyp) + { + pmin.SetToMin ( (*this) [pi] ); + pmax.SetToMax ( (*this) [pi] ); + } + } - double Mesh :: ElementError (int eli) const - { - const Element & el = volelements.Get(eli); - return CalcTetBadness (points.Get(el[0]), points.Get(el[1]), - points.Get(el[2]), points.Get(el[3]), -1); - } + double Mesh :: ElementError (int eli) const + { + const Element & el = volelements.Get(eli); + return CalcTetBadness (points.Get(el[0]), points.Get(el[1]), + points.Get(el[2]), points.Get(el[3]), -1); + } - void Mesh :: AddLockedPoint (PointIndex pi) - { - lockedpoints.Append (pi); - } + void Mesh :: AddLockedPoint (PointIndex pi) + { + lockedpoints.Append (pi); + } - void Mesh :: ClearLockedPoints () - { - lockedpoints.SetSize (0); - } + void Mesh :: ClearLockedPoints () + { + lockedpoints.SetSize (0); + } - void Mesh :: Compress () - { - int i, j; - Array op2np(GetNP()); - Array hpoints; - BitArrayChar pused(GetNP()); + void Mesh :: Compress () + { + int i, j; + Array op2np(GetNP()); + Array hpoints; + BitArrayChar pused(GetNP()); - /* + /* (*testout) << "volels: " << endl; for (i = 1; i <= volelements.Size(); i++) { @@ -3011,680 +3010,679 @@ namespace netgen (*testout) << endl; } (*testout) << "np: " << GetNP() << endl; - */ + */ - for (i = 0; i < volelements.Size(); i++) - if (volelements[i][0] <= PointIndex::BASE-1 || - volelements[i].IsDeleted()) - { - volelements.Delete(i); - i--; - } + for (i = 0; i < volelements.Size(); i++) + if (volelements[i][0] <= PointIndex::BASE-1 || + volelements[i].IsDeleted()) + { + volelements.Delete(i); + i--; + } - for (i = 0; i < surfelements.Size(); i++) - if (surfelements[i].IsDeleted()) - { - surfelements.Delete(i); - i--; - } + for (i = 0; i < surfelements.Size(); i++) + if (surfelements[i].IsDeleted()) + { + surfelements.Delete(i); + i--; + } - for (i = 0; i < segments.Size(); i++) - if (segments[i][0] <= PointIndex::BASE-1) - { - segments.Delete(i); - i--; - } + for (i = 0; i < segments.Size(); i++) + if (segments[i][0] <= PointIndex::BASE-1) + { + segments.Delete(i); + i--; + } - pused.Clear(); - for (i = 0; i < volelements.Size(); i++) - { - const Element & el = volelements[i]; - for (j = 0; j < el.GetNP(); j++) - pused.Set (el[j]); - } - - for (i = 0; i < surfelements.Size(); i++) - { - const Element2d & el = surfelements[i]; - for (j = 0; j < el.GetNP(); j++) - pused.Set (el[j]); - } - - for (i = 0; i < segments.Size(); i++) - { - const Segment & seg = segments[i]; - pused.Set (seg[0]); - pused.Set (seg[1]); - } - - for (i = 0; i < openelements.Size(); i++) - { - const Element2d & el = openelements[i]; - for (j = 0; j < el.GetNP(); j++) - pused.Set(el[j]); - } - - for (i = 0; i < lockedpoints.Size(); i++) - pused.Set (lockedpoints[i]); - - - /* - // compress points doesnt work for identified points ! - if (identifiedpoints) - { - for (i = 1; i <= identifiedpoints->GetNBags(); i++) - if (identifiedpoints->GetBagSize(i)) - { - pused.Set (); - break; - } - } - */ - // pused.Set(); - - - int npi = PointIndex::BASE-1; - - for (i = PointIndex::BASE; - i < points.Size()+PointIndex::BASE; i++) - if (pused.Test(i)) - { - npi++; - op2np[i] = npi; - hpoints.Append (points[i]); - } - else - op2np[i] = -1; - - - - points.SetSize(0); - for (i = 0; i < hpoints.Size(); i++) - points.Append (hpoints[i]); - - - for (i = 1; i <= volelements.Size(); i++) - { - Element & el = VolumeElement(i); - for (j = 0; j < el.GetNP(); j++) - el[j] = op2np[el[j]]; - } - - for (i = 1; i <= surfelements.Size(); i++) - { - Element2d & el = SurfaceElement(i); - for (j = 0; j < el.GetNP(); j++) - el[j] = op2np[el[j]]; - } - - for (i = 0; i < segments.Size(); i++) - { - Segment & seg = segments[i]; - seg[0] = op2np[seg[0]]; - seg[1] = op2np[seg[1]]; - } - - for (i = 1; i <= openelements.Size(); i++) - { - Element2d & el = openelements.Elem(i); - for (j = 0; j < el.GetNP(); j++) - el[j] = op2np[el[j]]; - } - - - for (i = 0; i < lockedpoints.Size(); i++) - lockedpoints[i] = op2np[lockedpoints[i]]; - - for (int i = 0; i < facedecoding.Size(); i++) - facedecoding[i].firstelement = -1; - for (int i = surfelements.Size()-1; i >= 0; i--) - { - int ind = surfelements[i].GetIndex(); - surfelements[i].next = facedecoding[ind-1].firstelement; - facedecoding[ind-1].firstelement = i; - } - - - CalcSurfacesOfNode(); - - - // FindOpenElements(); - timestamp = NextTimeStamp(); - - /* - (*testout) << "compress, done" << endl - << "np = " << points.Size() - << "ne = " << volelements.Size() << ", type.size = " << eltyps.Size() - << "volelements = " << volelements << endl; - */ - } - - - int Mesh :: CheckConsistentBoundary () const - { - int nf = GetNOpenElements(); - INDEX_2_HASHTABLE edges(nf+2); - INDEX_2 i2, i2s, edge; - int err = 0; - - for (int i = 1; i <= nf; i++) + pused.Clear(); + for (i = 0; i < volelements.Size(); i++) { - const Element2d & sel = OpenElement(i); + const Element & el = volelements[i]; + for (j = 0; j < el.GetNP(); j++) + pused.Set (el[j]); + } - for (int j = 1; j <= sel.GetNP(); j++) - { + for (i = 0; i < surfelements.Size(); i++) + { + const Element2d & el = surfelements[i]; + for (j = 0; j < el.GetNP(); j++) + pused.Set (el[j]); + } + + for (i = 0; i < segments.Size(); i++) + { + const Segment & seg = segments[i]; + pused.Set (seg[0]); + pused.Set (seg[1]); + } + + for (i = 0; i < openelements.Size(); i++) + { + const Element2d & el = openelements[i]; + for (j = 0; j < el.GetNP(); j++) + pused.Set(el[j]); + } + + for (i = 0; i < lockedpoints.Size(); i++) + pused.Set (lockedpoints[i]); + + + /* + // compress points doesnt work for identified points ! + if (identifiedpoints) + { + for (i = 1; i <= identifiedpoints->GetNBags(); i++) + if (identifiedpoints->GetBagSize(i)) + { + pused.Set (); + break; + } + } + */ + // pused.Set(); + + + int npi = PointIndex::BASE-1; + + for (i = PointIndex::BASE; + i < points.Size()+PointIndex::BASE; i++) + if (pused.Test(i)) + { + npi++; + op2np[i] = npi; + hpoints.Append (points[i]); + } + else + op2np[i] = -1; + + + + points.SetSize(0); + for (i = 0; i < hpoints.Size(); i++) + points.Append (hpoints[i]); + + + for (i = 1; i <= volelements.Size(); i++) + { + Element & el = VolumeElement(i); + for (j = 0; j < el.GetNP(); j++) + el[j] = op2np[el[j]]; + } + + for (i = 1; i <= surfelements.Size(); i++) + { + Element2d & el = SurfaceElement(i); + for (j = 0; j < el.GetNP(); j++) + el[j] = op2np[el[j]]; + } + + for (i = 0; i < segments.Size(); i++) + { + Segment & seg = segments[i]; + seg[0] = op2np[seg[0]]; + seg[1] = op2np[seg[1]]; + } + + for (i = 1; i <= openelements.Size(); i++) + { + Element2d & el = openelements.Elem(i); + for (j = 0; j < el.GetNP(); j++) + el[j] = op2np[el[j]]; + } + + + for (i = 0; i < lockedpoints.Size(); i++) + lockedpoints[i] = op2np[lockedpoints[i]]; + + for (int i = 0; i < facedecoding.Size(); i++) + facedecoding[i].firstelement = -1; + for (int i = surfelements.Size()-1; i >= 0; i--) + { + int ind = surfelements[i].GetIndex(); + surfelements[i].next = facedecoding[ind-1].firstelement; + facedecoding[ind-1].firstelement = i; + } + + + CalcSurfacesOfNode(); + + + // FindOpenElements(); + timestamp = NextTimeStamp(); + + /* + (*testout) << "compress, done" << endl + << "np = " << points.Size() + << "ne = " << volelements.Size() << ", type.size = " << eltyps.Size() + << "volelements = " << volelements << endl; + */ + } + + + int Mesh :: CheckConsistentBoundary () const + { + int nf = GetNOpenElements(); + INDEX_2_HASHTABLE edges(nf+2); + INDEX_2 i2, i2s, edge; + int err = 0; + + for (int i = 1; i <= nf; i++) + { + const Element2d & sel = OpenElement(i); + + for (int j = 1; j <= sel.GetNP(); j++) + { i2.I1() = sel.PNumMod(j); i2.I2() = sel.PNumMod(j+1); int sign = (i2.I2() > i2.I1()) ? 1 : -1; i2.Sort(); if (!edges.Used (i2)) - edges.Set (i2, 0); + edges.Set (i2, 0); edges.Set (i2, edges.Get(i2) + sign); - } + } } - for (int i = 1; i <= edges.GetNBags(); i++) - for (int j = 1; j <= edges.GetBagSize(i); j++) - { - int cnt = 0; - edges.GetData (i, j, i2, cnt); - if (cnt) + for (int i = 1; i <= edges.GetNBags(); i++) + for (int j = 1; j <= edges.GetBagSize(i); j++) + { + int cnt = 0; + edges.GetData (i, j, i2, cnt); + if (cnt) { - PrintError ("Edge ", i2.I1() , " - ", i2.I2(), " multiple times in surface mesh"); + PrintError ("Edge ", i2.I1() , " - ", i2.I2(), " multiple times in surface mesh"); - (*testout) << "Edge " << i2 << " multiple times in surface mesh" << endl; - i2s = i2; - i2s.Sort(); - for (int k = 1; k <= nf; k++) - { + (*testout) << "Edge " << i2 << " multiple times in surface mesh" << endl; + i2s = i2; + i2s.Sort(); + for (int k = 1; k <= nf; k++) + { const Element2d & sel = OpenElement(k); for (int l = 1; l <= sel.GetNP(); l++) - { - edge.I1() = sel.PNumMod(l); - edge.I2() = sel.PNumMod(l+1); - edge.Sort(); + { + edge.I1() = sel.PNumMod(l); + edge.I2() = sel.PNumMod(l+1); + edge.Sort(); - if (edge == i2s) + if (edge == i2s) (*testout) << "edge of element " << sel << endl; - } - } + } + } - err = 2; + err = 2; } - } + } - return err; - } + return err; + } - int Mesh :: CheckOverlappingBoundary () - { - int i, j, k; + int Mesh :: CheckOverlappingBoundary () + { + int i, j, k; - Point3d pmin, pmax; - GetBox (pmin, pmax); - Box3dTree setree(pmin, pmax); - Array inters; + Point3d pmin, pmax; + GetBox (pmin, pmax); + Box3dTree setree(pmin, pmax); + Array inters; - bool overlap = 0; - bool incons_layers = 0; + bool overlap = 0; + bool incons_layers = 0; - for (i = 1; i <= GetNSE(); i++) - SurfaceElement(i).badel = 0; + for (i = 1; i <= GetNSE(); i++) + SurfaceElement(i).badel = 0; - for (i = 1; i <= GetNSE(); i++) + for (i = 1; i <= GetNSE(); i++) { - const Element2d & tri = SurfaceElement(i); + const Element2d & tri = SurfaceElement(i); - Point3d tpmin (Point(tri[0])); - Point3d tpmax (tpmin); + Point3d tpmin (Point(tri[0])); + Point3d tpmax (tpmin); - for (k = 1; k < tri.GetNP(); k++) - { + for (k = 1; k < tri.GetNP(); k++) + { tpmin.SetToMin (Point (tri[k])); tpmax.SetToMax (Point (tri[k])); - } - Vec3d diag(tpmin, tpmax); + } + Vec3d diag(tpmin, tpmax); - tpmax = tpmax + 0.1 * diag; - tpmin = tpmin - 0.1 * diag; + tpmax = tpmax + 0.1 * diag; + tpmin = tpmin - 0.1 * diag; - setree.Insert (tpmin, tpmax, i); + setree.Insert (tpmin, tpmax, i); } - for (i = 1; i <= GetNSE(); i++) + for (i = 1; i <= GetNSE(); i++) { - const Element2d & tri = SurfaceElement(i); + const Element2d & tri = SurfaceElement(i); - Point3d tpmin (Point(tri[0])); - Point3d tpmax (tpmin); + Point3d tpmin (Point(tri[0])); + Point3d tpmax (tpmin); - for (k = 1; k < tri.GetNP(); k++) - { + for (k = 1; k < tri.GetNP(); k++) + { tpmin.SetToMin (Point (tri[k])); tpmax.SetToMax (Point (tri[k])); - } + } - setree.GetIntersecting (tpmin, tpmax, inters); + setree.GetIntersecting (tpmin, tpmax, inters); - for (j = 1; j <= inters.Size(); j++) - { + for (j = 1; j <= inters.Size(); j++) + { const Element2d & tri2 = SurfaceElement(inters.Get(j)); if ( (*this)[tri[0]].GetLayer() != (*this)[tri2[0]].GetLayer()) - continue; + continue; if ( (*this)[tri[0]].GetLayer() != (*this)[tri[1]].GetLayer() || - (*this)[tri[0]].GetLayer() != (*this)[tri[2]].GetLayer()) - { - incons_layers = 1; - cout << "inconsistent layers in triangle" << endl; - } + (*this)[tri[0]].GetLayer() != (*this)[tri[2]].GetLayer()) + { + incons_layers = 1; + cout << "inconsistent layers in triangle" << endl; + } const netgen::Point<3> *trip1[3], *trip2[3]; for (k = 1; k <= 3; k++) - { - trip1[k-1] = &Point (tri.PNum(k)); - trip2[k-1] = &Point (tri2.PNum(k)); - } + { + trip1[k-1] = &Point (tri.PNum(k)); + trip2[k-1] = &Point (tri2.PNum(k)); + } if (IntersectTriangleTriangle (&trip1[0], &trip2[0])) - { - overlap = 1; - PrintWarning ("Intersecting elements " - ,i, " and ", inters.Get(j)); + { + overlap = 1; + PrintWarning ("Intersecting elements " + ,i, " and ", inters.Get(j)); - (*testout) << "Intersecting: " << endl; - (*testout) << "openelement " << i << " with open element " << inters.Get(j) << endl; + (*testout) << "Intersecting: " << endl; + (*testout) << "openelement " << i << " with open element " << inters.Get(j) << endl; - cout << "el1 = " << tri << endl; - cout << "el2 = " << tri2 << endl; - cout << "layer1 = " << (*this)[tri[0]].GetLayer() << endl; - cout << "layer2 = " << (*this)[tri2[0]].GetLayer() << endl; + cout << "el1 = " << tri << endl; + cout << "el2 = " << tri2 << endl; + cout << "layer1 = " << (*this)[tri[0]].GetLayer() << endl; + cout << "layer2 = " << (*this)[tri2[0]].GetLayer() << endl; - for (k = 1; k <= 3; k++) + for (k = 1; k <= 3; k++) (*testout) << tri.PNum(k) << " "; - (*testout) << endl; - for (k = 1; k <= 3; k++) + (*testout) << endl; + for (k = 1; k <= 3; k++) (*testout) << tri2.PNum(k) << " "; - (*testout) << endl; + (*testout) << endl; - for (k = 0; k <= 2; k++) + for (k = 0; k <= 2; k++) (*testout) << *trip1[k] << " "; - (*testout) << endl; - for (k = 0; k <= 2; k++) + (*testout) << endl; + for (k = 0; k <= 2; k++) (*testout) << *trip2[k] << " "; - (*testout) << endl; + (*testout) << endl; - (*testout) << "Face1 = " << GetFaceDescriptor(tri.GetIndex()) << endl; - (*testout) << "Face1 = " << GetFaceDescriptor(tri2.GetIndex()) << endl; + (*testout) << "Face1 = " << GetFaceDescriptor(tri.GetIndex()) << endl; + (*testout) << "Face1 = " << GetFaceDescriptor(tri2.GetIndex()) << endl; - /* - INDEX_3 i3(tri.PNum(1), tri.PNum(2), tri.PNum(3)); - i3.Sort(); - for (k = 1; k <= GetNSE(); k++) - { - const Element2d & el2 = SurfaceElement(k); - INDEX_3 i3b(el2.PNum(1), el2.PNum(2), el2.PNum(3)); - i3b.Sort(); - if (i3 == i3b) - { - SurfaceElement(k).badel = 1; - } - } - */ - SurfaceElement(i).badel = 1; - SurfaceElement(inters.Get(j)).badel = 1; - } - } + /* + INDEX_3 i3(tri.PNum(1), tri.PNum(2), tri.PNum(3)); + i3.Sort(); + for (k = 1; k <= GetNSE(); k++) + { + const Element2d & el2 = SurfaceElement(k); + INDEX_3 i3b(el2.PNum(1), el2.PNum(2), el2.PNum(3)); + i3b.Sort(); + if (i3 == i3b) + { + SurfaceElement(k).badel = 1; + } + } + */ + SurfaceElement(i).badel = 1; + SurfaceElement(inters.Get(j)).badel = 1; + } + } } - // bug 'fix' - if (incons_layers) overlap = 0; + // bug 'fix' + if (incons_layers) overlap = 0; - return overlap; - } + return overlap; + } - int Mesh :: CheckVolumeMesh () const - { - PrintMessage (3, "Checking volume mesh"); + int Mesh :: CheckVolumeMesh () const + { + PrintMessage (3, "Checking volume mesh"); - int ne = GetNE(); - DenseMatrix dtrans(3,3); - int i, j; + int ne = GetNE(); + DenseMatrix dtrans(3,3); + int i, j; - PrintMessage (5, "elements: ", ne); - for (i = 1; i <= ne; i++) + PrintMessage (5, "elements: ", ne); + for (i = 1; i <= ne; i++) { - Element & el = (Element&) VolumeElement(i); - el.flags.badel = 0; - int nip = el.GetNIP(); - for (j = 1; j <= nip; j++) - { + Element & el = (Element&) VolumeElement(i); + el.flags.badel = 0; + int nip = el.GetNIP(); + for (j = 1; j <= nip; j++) + { el.GetTransformation (j, Points(), dtrans); double det = dtrans.Det(); if (det > 0) - { - PrintError ("Element ", i , " has wrong orientation"); - el.flags.badel = 1; - } - } + { + PrintError ("Element ", i , " has wrong orientation"); + el.flags.badel = 1; + } + } } - return 0; - } + return 0; + } - bool Mesh :: LegalTrig (const Element2d & el) const - { - return 1; - if ( /* hp */ 1) // needed for old, simple hp-refinement + bool Mesh :: LegalTrig (const Element2d & el) const + { + return 1; + if ( /* hp */ 1) // needed for old, simple hp-refinement { - // trigs with 2 or more segments are illegal - int i; - int nseg = 0; + // trigs with 2 or more segments are illegal + int i; + int nseg = 0; - if (!segmentht) - { + if (!segmentht) + { cerr << "no segmentht allocated" << endl; return 0; - } + } - // Point3d cp(0.5, 0.5, 0.5); - for (i = 1; i <= 3; i++) - { + // Point3d cp(0.5, 0.5, 0.5); + for (i = 1; i <= 3; i++) + { INDEX_2 i2(el.PNumMod (i), el.PNumMod (i+1)); i2.Sort(); if (segmentht -> Used (i2)) - nseg++; - } - if (nseg >= 2) - return 0; + nseg++; + } + if (nseg >= 2) + return 0; } - return 1; - } + return 1; + } - /// - bool Mesh :: LegalTet2 (Element & el) const - { - // static int timer1 = NgProfiler::CreateTimer ("Legaltet2"); + /// + bool Mesh :: LegalTet2 (Element & el) const + { + // static int timer1 = NgProfiler::CreateTimer ("Legaltet2"); - // Test, whether 4 points have a common surface plus - // at least 4 edges at the boundary + // Test, whether 4 points have a common surface plus + // at least 4 edges at the boundary - if(!boundaryedges) - const_cast(this)->BuildBoundaryEdges(); + if(!boundaryedges) + const_cast(this)->BuildBoundaryEdges(); - // non-tets are always legal - if (el.GetType() != TET) + // non-tets are always legal + if (el.GetType() != TET) { - el.SetLegal (1); - return 1; + el.SetLegal (1); + return 1; } - POINTTYPE pointtype[4]; - for(int i = 0; i < 4; i++) - pointtype[i] = (*this)[el[i]].Type(); + POINTTYPE pointtype[4]; + for(int i = 0; i < 4; i++) + pointtype[i] = (*this)[el[i]].Type(); - // element has at least 2 inner points ---> legal - int cnti = 0; - for (int j = 0; j < 4; j++) - if ( pointtype[j] == INNERPOINT) - { - cnti++; - if (cnti >= 2) + // element has at least 2 inner points ---> legal + int cnti = 0; + for (int j = 0; j < 4; j++) + if ( pointtype[j] == INNERPOINT) + { + cnti++; + if (cnti >= 2) { - el.SetLegal (1); - return 1; + el.SetLegal (1); + return 1; } - } + } - // which faces are boundary faces ? - int bface[4]; - for (int i = 0; i < 4; i++) - { - bface[i] = surfelementht->Used (INDEX_3::Sort(el[gftetfacesa[i][0]], - el[gftetfacesa[i][1]], - el[gftetfacesa[i][2]])); - } + // which faces are boundary faces ? + int bface[4]; + for (int i = 0; i < 4; i++) + { + bface[i] = surfelementht->Used (INDEX_3::Sort(el[gftetfacesa[i][0]], + el[gftetfacesa[i][1]], + el[gftetfacesa[i][2]])); + } - int bedge[4][4]; - int segedge[4][4]; - static const int pi3map[4][4] = { { -1, 2, 1, 1 }, - { 2, -1, 0, 0 }, - { 1, 0, -1, 0 }, - { 1, 0, 0, -1 } }; + int bedge[4][4]; + int segedge[4][4]; + static const int pi3map[4][4] = { { -1, 2, 1, 1 }, + { 2, -1, 0, 0 }, + { 1, 0, -1, 0 }, + { 1, 0, 0, -1 } }; - static const int pi4map[4][4] = { { -1, 3, 3, 2 }, - { 3, -1, 3, 2 }, - { 3, 3, -1, 1 }, - { 2, 2, 1, -1 } }; + static const int pi4map[4][4] = { { -1, 3, 3, 2 }, + { 3, -1, 3, 2 }, + { 3, 3, -1, 1 }, + { 2, 2, 1, -1 } }; - for (int i = 0; i < 4; i++) - for (int j = 0; j < i; j++) + for (int i = 0; i < 4; i++) + for (int j = 0; j < i; j++) + { + bool sege = false, be = false; + + int pos = boundaryedges -> Position(INDEX_2::Sort(el[i], el[j])); + if (pos) { - bool sege = false, be = false; - - int pos = boundaryedges -> Position(INDEX_2::Sort(el[i], el[j])); - if (pos) - { - be = true; - if (boundaryedges -> GetData(pos) == 2) - sege = true; - } - - segedge[j][i] = segedge[i][j] = sege; - bedge[j][i] = bedge[i][j] = be; + be = true; + if (boundaryedges -> GetData(pos) == 2) + sege = true; } - // two boundary faces and no edge is illegal - for (int i = 0; i < 3; i++) - for (int j = i+1; j < 4; j++) - { - if (bface[i] && bface[j]) - if (!segedge[pi3map[i][j]][pi4map[i][j]]) - { - // 2 boundary faces withoud edge in between - el.SetLegal (0); - return 0; - } - } + segedge[j][i] = segedge[i][j] = sege; + bedge[j][i] = bedge[i][j] = be; + } - // three boundary edges meeting in a Surface point - for (int i = 0; i < 4; i++) - { - bool alledges = 1; - if ( pointtype[i] == SURFACEPOINT) - { - bool alledges = 1; - for (int j = 0; j < 4; j++) - if (j != i && !bedge[i][j]) - { - alledges = 0; - break; - } - if (alledges) - { - // cout << "tet illegal due to unmarked node" << endl; - el.SetLegal (0); - return 0; - } - } - } + // two boundary faces and no edge is illegal + for (int i = 0; i < 3; i++) + for (int j = i+1; j < 4; j++) + { + if (bface[i] && bface[j]) + if (!segedge[pi3map[i][j]][pi4map[i][j]]) + { + // 2 boundary faces withoud edge in between + el.SetLegal (0); + return 0; + } + } - - - for (int fnr = 0; fnr < 4; fnr++) - if (!bface[fnr]) - for (int i = 0; i < 4; i++) - if (i != fnr) - { - int pi1 = pi3map[i][fnr]; - int pi2 = pi4map[i][fnr]; - - if ( pointtype[i] == SURFACEPOINT) - { - // two connected edges on surface, but no face - if (bedge[i][pi1] && bedge[i][pi2]) - { - el.SetLegal (0); - return 0; - } - } - - if ( pointtype[i] == EDGEPOINT) - { - // connected surface edge and edge edge, but no face - if (bedge[i][pi1] && segedge[i][pi2] || - bedge[i][pi2] && segedge[i][pi1]) - { - el.SetLegal (0); - return 0; - } - } - - } - - - el.SetLegal (1); - return 1; - - } - - - - int Mesh :: GetNDomains() const - { - int ndom = 0; - - for (int k = 0; k < facedecoding.Size(); k++) + // three boundary edges meeting in a Surface point + for (int i = 0; i < 4; i++) { - if (facedecoding[k].DomainIn() > ndom) - ndom = facedecoding[k].DomainIn(); - if (facedecoding[k].DomainOut() > ndom) - ndom = facedecoding[k].DomainOut(); + if ( pointtype[i] == SURFACEPOINT) + { + bool alledges = 1; + for (int j = 0; j < 4; j++) + if (j != i && !bedge[i][j]) + { + alledges = 0; + break; + } + if (alledges) + { + // cout << "tet illegal due to unmarked node" << endl; + el.SetLegal (0); + return 0; + } + } } - return ndom; - } + + + for (int fnr = 0; fnr < 4; fnr++) + if (!bface[fnr]) + for (int i = 0; i < 4; i++) + if (i != fnr) + { + int pi1 = pi3map[i][fnr]; + int pi2 = pi4map[i][fnr]; + + if ( pointtype[i] == SURFACEPOINT) + { + // two connected edges on surface, but no face + if (bedge[i][pi1] && bedge[i][pi2]) + { + el.SetLegal (0); + return 0; + } + } + + if ( pointtype[i] == EDGEPOINT) + { + // connected surface edge and edge edge, but no face + if ( (bedge[i][pi1] && segedge[i][pi2]) || + (bedge[i][pi2] && segedge[i][pi1]) ) + { + el.SetLegal (0); + return 0; + } + } + + } + + + el.SetLegal (1); + return 1; + + } - void Mesh :: SurfaceMeshOrientation () - { - int i, j; - int nse = GetNSE(); + int Mesh :: GetNDomains() const + { + int ndom = 0; - BitArray used(nse); - used.Clear(); - INDEX_2_HASHTABLE edges(nse+1); - - bool haschanged = 0; - - - const Element2d & tri = SurfaceElement(1); - for (j = 1; j <= 3; j++) + for (int k = 0; k < facedecoding.Size(); k++) { - INDEX_2 i2(tri.PNumMod(j), tri.PNumMod(j+1)); - edges.Set (i2, 1); + if (facedecoding[k].DomainIn() > ndom) + ndom = facedecoding[k].DomainIn(); + if (facedecoding[k].DomainOut() > ndom) + ndom = facedecoding[k].DomainOut(); } - used.Set(1); - bool unused; - do + return ndom; + } + + + + void Mesh :: SurfaceMeshOrientation () + { + int i, j; + int nse = GetNSE(); + + BitArray used(nse); + used.Clear(); + INDEX_2_HASHTABLE edges(nse+1); + + bool haschanged = 0; + + + const Element2d & tri = SurfaceElement(1); + for (j = 1; j <= 3; j++) { - bool changed; - do - { + INDEX_2 i2(tri.PNumMod(j), tri.PNumMod(j+1)); + edges.Set (i2, 1); + } + used.Set(1); + + bool unused; + do + { + bool changed; + do + { changed = 0; for (i = 1; i <= nse; i++) - if (!used.Test(i)) - { + if (!used.Test(i)) + { Element2d & el = surfelements.Elem(i); int found = 0, foundrev = 0; for (j = 1; j <= 3; j++) - { - INDEX_2 i2(el.PNumMod(j), el.PNumMod(j+1)); - if (edges.Used(i2)) + { + INDEX_2 i2(el.PNumMod(j), el.PNumMod(j+1)); + if (edges.Used(i2)) foundrev = 1; - swap (i2.I1(), i2.I2()); - if (edges.Used(i2)) + swap (i2.I1(), i2.I2()); + if (edges.Used(i2)) found = 1; - } + } if (found || foundrev) - { - if (foundrev) + { + if (foundrev) swap (el.PNum(2), el.PNum(3)); - changed = 1; - for (j = 1; j <= 3; j++) - { - INDEX_2 i2(el.PNumMod(j), el.PNumMod(j+1)); - edges.Set (i2, 1); - } - used.Set (i); - } - } - if (changed) - haschanged = 1; - } - while (changed); + changed = 1; + for (j = 1; j <= 3; j++) + { + INDEX_2 i2(el.PNumMod(j), el.PNumMod(j+1)); + edges.Set (i2, 1); + } + used.Set (i); + } + } + if (changed) + haschanged = 1; + } + while (changed); - unused = 0; - for (i = 1; i <= nse; i++) - if (!used.Test(i)) + unused = 0; + for (i = 1; i <= nse; i++) + if (!used.Test(i)) { - unused = 1; - const Element2d & tri = SurfaceElement(i); - for (j = 1; j <= 3; j++) - { + unused = 1; + const Element2d & tri = SurfaceElement(i); + for (j = 1; j <= 3; j++) + { INDEX_2 i2(tri.PNumMod(j), tri.PNumMod(j+1)); edges.Set (i2, 1); - } - used.Set(i); - break; + } + used.Set(i); + break; } } - while (unused); + while (unused); - if (haschanged) - timestamp = NextTimeStamp(); - } + if (haschanged) + timestamp = NextTimeStamp(); + } - void Mesh :: Split2Tets() - { - PrintMessage (1, "Split To Tets"); - bool has_prisms = 0; + void Mesh :: Split2Tets() + { + PrintMessage (1, "Split To Tets"); + bool has_prisms = 0; - int oldne = GetNE(); - for (int i = 1; i <= oldne; i++) + int oldne = GetNE(); + for (int i = 1; i <= oldne; i++) { - Element el = VolumeElement(i); + Element el = VolumeElement(i); - if (el.GetType() == PRISM) - { + if (el.GetType() == PRISM) + { // prism, to 3 tets // make minimal node to node 1 @@ -3693,94 +3691,94 @@ namespace netgen minpnum = GetNP() + 1; for (int j = 1; j <= 6; j++) - { - if (el.PNum(j) < minpnum) - { - minpnum = el.PNum(j); - minpi = j; - } - } + { + if (el.PNum(j) < minpnum) + { + minpnum = el.PNum(j); + minpi = j; + } + } if (minpi >= 4) - { - for (int j = 1; j <= 3; j++) + { + for (int j = 1; j <= 3; j++) swap (el.PNum(j), el.PNum(j+3)); - minpi -= 3; - } + minpi -= 3; + } while (minpi > 1) - { - int hi = 0; - for (int j = 0; j <= 3; j+= 3) - { - hi = el.PNum(1+j); - el.PNum(1+j) = el.PNum(2+j); - el.PNum(2+j) = el.PNum(3+j); - el.PNum(3+j) = hi; - } - minpi--; - } + { + int hi = 0; + for (int j = 0; j <= 3; j+= 3) + { + hi = el.PNum(1+j); + el.PNum(1+j) = el.PNum(2+j); + el.PNum(2+j) = el.PNum(3+j); + el.PNum(3+j) = hi; + } + minpi--; + } /* - version 1: edge from pi2 to pi6, - version 2: edge from pi3 to pi5, + version 1: edge from pi2 to pi6, + version 2: edge from pi3 to pi5, */ static const int ntets[2][12] = - { { 1, 4, 5, 6, 1, 2, 3, 6, 1, 2, 5, 6 }, - { 1, 4, 5, 6, 1, 2, 3, 5, 3, 1, 5, 6 } }; + { { 1, 4, 5, 6, 1, 2, 3, 6, 1, 2, 5, 6 }, + { 1, 4, 5, 6, 1, 2, 3, 5, 3, 1, 5, 6 } }; const int * min2pi; if (min2 (el.PNum(2), el.PNum(6)) < - min2 (el.PNum(3), el.PNum(5))) - { - min2pi = &ntets[0][0]; - // (*testout) << "version 1 "; - } + min2 (el.PNum(3), el.PNum(5))) + { + min2pi = &ntets[0][0]; + // (*testout) << "version 1 "; + } else - { - min2pi = &ntets[1][0]; - // (*testout) << "version 2 "; - } + { + min2pi = &ntets[1][0]; + // (*testout) << "version 2 "; + } int firsttet = 1; for (int j = 1; j <= 3; j++) - { - Element nel(TET); - for (int k = 1; k <= 4; k++) + { + Element nel(TET); + for (int k = 1; k <= 4; k++) nel.PNum(k) = el.PNum(min2pi[4 * j + k - 5]); - nel.SetIndex (el.GetIndex()); + nel.SetIndex (el.GetIndex()); - int legal = 1; - for (int k = 1; k <= 3; k++) + int legal = 1; + for (int k = 1; k <= 3; k++) for (int l = k+1; l <= 4; l++) - if (nel.PNum(k) == nel.PNum(l)) - legal = 0; + if (nel.PNum(k) == nel.PNum(l)) + legal = 0; - // (*testout) << nel << " "; - if (legal) - { - if (firsttet) + // (*testout) << nel << " "; + if (legal) { - VolumeElement(i) = nel; - firsttet = 0; + if (firsttet) + { + VolumeElement(i) = nel; + firsttet = 0; + } + else + { + AddVolumeElement(nel); + } } - else - { - AddVolumeElement(nel); - } - } - } + } if (firsttet) cout << "no legal"; (*testout) << endl; - } + } - else if (el.GetType() == HEX) - { + else if (el.GetType() == HEX) + { // hex to A) 2 prisms or B) to 5 tets // make minimal node to node 1 @@ -3789,216 +3787,216 @@ namespace netgen minpnum = GetNP() + 1; for (int j = 1; j <= 8; j++) - { - if (el.PNum(j) < minpnum) - { - minpnum = el.PNum(j); - minpi = j; - } - } + { + if (el.PNum(j) < minpnum) + { + minpnum = el.PNum(j); + minpi = j; + } + } if (minpi >= 5) - { - for (int j = 1; j <= 4; j++) + { + for (int j = 1; j <= 4; j++) swap (el.PNum(j), el.PNum(j+4)); - minpi -= 4; - } + minpi -= 4; + } while (minpi > 1) - { - int hi = 0; - for (int j = 0; j <= 4; j+= 4) - { - hi = el.PNum(1+j); - el.PNum(1+j) = el.PNum(2+j); - el.PNum(2+j) = el.PNum(3+j); - el.PNum(3+j) = el.PNum(4+j); - el.PNum(4+j) = hi; - } - minpi--; - } + { + int hi = 0; + for (int j = 0; j <= 4; j+= 4) + { + hi = el.PNum(1+j); + el.PNum(1+j) = el.PNum(2+j); + el.PNum(2+j) = el.PNum(3+j); + el.PNum(3+j) = el.PNum(4+j); + el.PNum(4+j) = hi; + } + minpi--; + } static const int to_prisms[3][12] = - { { 0, 1, 2, 4, 5, 6, 0, 2, 3, 4, 6, 7 }, - { 0, 1, 5, 3, 2, 6, 0, 5, 4, 3, 6, 7 }, - { 0, 7, 4, 1, 6, 5, 0, 3, 7, 1, 2, 6 }, - }; + { { 0, 1, 2, 4, 5, 6, 0, 2, 3, 4, 6, 7 }, + { 0, 1, 5, 3, 2, 6, 0, 5, 4, 3, 6, 7 }, + { 0, 7, 4, 1, 6, 5, 0, 3, 7, 1, 2, 6 }, + }; const int * min2pi = 0; if (min2 (el[4], el[6]) < min2 (el[5], el[7])) - min2pi = &to_prisms[0][0]; + min2pi = &to_prisms[0][0]; else if (min2 (el[3], el[6]) < min2 (el[2], el[7])) - min2pi = &to_prisms[1][0]; + min2pi = &to_prisms[1][0]; else if (min2 (el[1], el[6]) < min2 (el[2], el[5])) - min2pi = &to_prisms[2][0]; + min2pi = &to_prisms[2][0]; if (min2pi) - { - has_prisms = 1; - for (int j = 0; j < 2; j++) - { - Element nel(PRISM); - for (int k = 0; k < 6; k++) - nel[k] = el[min2pi[6*j + k]]; - nel.SetIndex (el.GetIndex()); + { + has_prisms = 1; + for (int j = 0; j < 2; j++) + { + Element nel(PRISM); + for (int k = 0; k < 6; k++) + nel[k] = el[min2pi[6*j + k]]; + nel.SetIndex (el.GetIndex()); - if (j == 0) - VolumeElement(i) = nel; - else - AddVolumeElement(nel); - } - } + if (j == 0) + VolumeElement(i) = nel; + else + AddVolumeElement(nel); + } + } else - { - // split to 5 tets + { + // split to 5 tets - static const int to_tets[20] = - { - 1, 2, 0, 5, - 3, 0, 2, 7, - 4, 5, 7, 0, - 6, 7, 5, 2, - 0, 2, 7, 5 - }; + static const int to_tets[20] = + { + 1, 2, 0, 5, + 3, 0, 2, 7, + 4, 5, 7, 0, + 6, 7, 5, 2, + 0, 2, 7, 5 + }; - for (int j = 0; j < 5; j++) - { - Element nel(TET); - for (int k = 0; k < 4; k++) - nel[k] = el[to_tets[4*j + k]]; - nel.SetIndex (el.GetIndex()); + for (int j = 0; j < 5; j++) + { + Element nel(TET); + for (int k = 0; k < 4; k++) + nel[k] = el[to_tets[4*j + k]]; + nel.SetIndex (el.GetIndex()); - if (j == 0) - VolumeElement(i) = nel; - else - AddVolumeElement(nel); - } + if (j == 0) + VolumeElement(i) = nel; + else + AddVolumeElement(nel); + } - } - } + } + } - else if (el.GetType() == PYRAMID) - { + else if (el.GetType() == PYRAMID) + { // pyramid, to 2 tets // cout << "pyramid: " << el << endl; static const int ntets[2][8] = - { { 1, 2, 3, 5, 1, 3, 4, 5 }, - { 1, 2, 4, 5, 4, 2, 3, 5 }}; + { { 1, 2, 3, 5, 1, 3, 4, 5 }, + { 1, 2, 4, 5, 4, 2, 3, 5 }}; const int * min2pi; if (min2 (el[0], el[2]) < min2 (el[1], el[3])) - min2pi = &ntets[0][0]; + min2pi = &ntets[0][0]; else - min2pi = &ntets[1][0]; + min2pi = &ntets[1][0]; bool firsttet = 1; for (int j = 0; j < 2; j++) - { - Element nel(TET); - for (int k = 0; k < 4; k++) + { + Element nel(TET); + for (int k = 0; k < 4; k++) nel[k] = el[min2pi[4*j + k]-1]; - nel.SetIndex (el.GetIndex()); + nel.SetIndex (el.GetIndex()); - // cout << "pyramid-tet: " << nel << endl; + // cout << "pyramid-tet: " << nel << endl; - bool legal = 1; - for (int k = 0; k < 3; k++) + bool legal = 1; + for (int k = 0; k < 3; k++) for (int l = k+1; l < 4; l++) - if (nel[k] == nel[l]) - legal = 0; + if (nel[k] == nel[l]) + legal = 0; - if (legal) - { - (*testout) << nel << " "; - if (firsttet) - VolumeElement(i) = nel; - else - AddVolumeElement(nel); + if (legal) + { + (*testout) << nel << " "; + if (firsttet) + VolumeElement(i) = nel; + else + AddVolumeElement(nel); - firsttet = 0; - } - } + firsttet = 0; + } + } if (firsttet) cout << "no legal"; (*testout) << endl; - } + } } - int oldnse = GetNSE(); - for (int i = 1; i <= oldnse; i++) + int oldnse = GetNSE(); + for (int i = 1; i <= oldnse; i++) { - Element2d el = SurfaceElement(i); - if (el.GetNP() == 4) - { + Element2d el = SurfaceElement(i); + if (el.GetNP() == 4) + { (*testout) << "split el: " << el << " to "; static const int ntris[2][6] = - { { 1, 2, 3, 1, 3, 4 }, - { 1, 2, 4, 4, 2, 3 }}; + { { 1, 2, 3, 1, 3, 4 }, + { 1, 2, 4, 4, 2, 3 }}; const int * min2pi; if (min2 (el.PNum(1), el.PNum(3)) < - min2 (el.PNum(2), el.PNum(4))) - min2pi = &ntris[0][0]; + min2 (el.PNum(2), el.PNum(4))) + min2pi = &ntris[0][0]; else - min2pi = &ntris[1][0]; + min2pi = &ntris[1][0]; for (int j = 0; j <6; j++) - (*testout) << min2pi[j] << " "; + (*testout) << min2pi[j] << " "; int firsttri = 1; for (int j = 1; j <= 2; j++) - { - Element2d nel(3); - for (int k = 1; k <= 3; k++) + { + Element2d nel(3); + for (int k = 1; k <= 3; k++) nel.PNum(k) = el.PNum(min2pi[3 * j + k - 4]); - nel.SetIndex (el.GetIndex()); + nel.SetIndex (el.GetIndex()); - int legal = 1; - for (int k = 1; k <= 2; k++) + int legal = 1; + for (int k = 1; k <= 2; k++) for (int l = k+1; l <= 3; l++) - if (nel.PNum(k) == nel.PNum(l)) - legal = 0; + if (nel.PNum(k) == nel.PNum(l)) + legal = 0; - if (legal) - { - (*testout) << nel << " "; - if (firsttri) + if (legal) { - SurfaceElement(i) = nel; - firsttri = 0; + (*testout) << nel << " "; + if (firsttri) + { + SurfaceElement(i) = nel; + firsttri = 0; + } + else + { + AddSurfaceElement(nel); + } } - else - { - AddSurfaceElement(nel); - } - } - } + } (*testout) << endl; - } + } } - if (has_prisms) + if (has_prisms) - Split2Tets(); + Split2Tets(); - else + else { - for (int i = 1; i <= GetNE(); i++) - { + for (int i = 1; i <= GetNE(); i++) + { Element & el = VolumeElement(i); const Point3d & p1 = Point (el.PNum(1)); const Point3d & p2 = Point (el.PNum(2)); @@ -4006,192 +4004,192 @@ namespace netgen const Point3d & p4 = Point (el.PNum(4)); double vol = (Vec3d (p1, p2) * - Cross (Vec3d (p1, p3), Vec3d(p1, p4))); + Cross (Vec3d (p1, p3), Vec3d(p1, p4))); if (vol > 0) - swap (el.PNum(3), el.PNum(4)); - } + swap (el.PNum(3), el.PNum(4)); + } - UpdateTopology(); - timestamp = NextTimeStamp(); + UpdateTopology(); + timestamp = NextTimeStamp(); } - } + } - void Mesh :: BuildElementSearchTree () - { - if (elementsearchtreets == GetTimeStamp()) - return; + void Mesh :: BuildElementSearchTree () + { + if (elementsearchtreets == GetTimeStamp()) + return; - NgLock lock(mutex); - lock.Lock(); + NgLock lock(mutex); + lock.Lock(); - PrintMessage (4, "Rebuild element searchtree"); + PrintMessage (4, "Rebuild element searchtree"); - if (elementsearchtree) - delete elementsearchtree; - elementsearchtree = NULL; + if (elementsearchtree) + delete elementsearchtree; + elementsearchtree = NULL; - Box3d box; - int i, j; - int ne = GetNE(); - if (!ne) + Box3d box; + int i, j; + int ne = GetNE(); + if (!ne) { - lock.UnLock(); - return; + lock.UnLock(); + return; } - box.SetPoint (Point (VolumeElement(1).PNum(1))); - for (i = 1; i <= ne; i++) + box.SetPoint (Point (VolumeElement(1).PNum(1))); + for (i = 1; i <= ne; i++) { - const Element & el = VolumeElement(i); - for (j = 1; j <= el.GetNP(); j++) - box.AddPoint (Point (el.PNum(j))); + const Element & el = VolumeElement(i); + for (j = 1; j <= el.GetNP(); j++) + box.AddPoint (Point (el.PNum(j))); } - box.Increase (1.01 * box.CalcDiam()); - elementsearchtree = new Box3dTree (box.PMin(), box.PMax()); + box.Increase (1.01 * box.CalcDiam()); + elementsearchtree = new Box3dTree (box.PMin(), box.PMax()); - for (i = 1; i <= ne; i++) + for (i = 1; i <= ne; i++) { - const Element & el = VolumeElement(i); - box.SetPoint (Point (el.PNum(1))); - for (j = 1; j <= el.GetNP(); j++) - box.AddPoint (Point (el.PNum(j))); + const Element & el = VolumeElement(i); + box.SetPoint (Point (el.PNum(1))); + for (j = 1; j <= el.GetNP(); j++) + box.AddPoint (Point (el.PNum(j))); - elementsearchtree -> Insert (box.PMin(), box.PMax(), i); + elementsearchtree -> Insert (box.PMin(), box.PMax(), i); } - elementsearchtreets = GetTimeStamp(); + elementsearchtreets = GetTimeStamp(); - lock.UnLock(); - } + lock.UnLock(); + } - bool Mesh :: PointContainedIn2DElement(const Point3d & p, - double lami[3], - const int element, - bool consider3D) const - { - static Vec3d col1, col2, col3; - static Vec3d rhs, sol; - const double eps = 1e-6; + bool Mesh :: PointContainedIn2DElement(const Point3d & p, + double lami[3], + const int element, + bool consider3D) const + { + static Vec3d col1, col2, col3; + static Vec3d rhs, sol; + const double eps = 1e-6; - static Array loctrigs; + static Array loctrigs; - //SZ - if(SurfaceElement(element).GetType()==QUAD) + //SZ + if(SurfaceElement(element).GetType()==QUAD) { - const Element2d & el = SurfaceElement(element); + const Element2d & el = SurfaceElement(element); - const Point3d & p1 = Point(el.PNum(1)); - const Point3d & p2 = Point(el.PNum(2)); - const Point3d & p3 = Point(el.PNum(3)); - const Point3d & p4 = Point(el.PNum(4)); + const Point3d & p1 = Point(el.PNum(1)); + const Point3d & p2 = Point(el.PNum(2)); + const Point3d & p3 = Point(el.PNum(3)); + const Point3d & p4 = Point(el.PNum(4)); - // Coefficients of Bilinear Mapping from Ref-Elem to global Elem - // X = a + b x + c y + d x y - Vec3d a = p1; - Vec3d b = p2 - a; - Vec3d c = p4 - a; - Vec3d d = p3 - a - b - c; + // Coefficients of Bilinear Mapping from Ref-Elem to global Elem + // X = a + b x + c y + d x y + Vec3d a = p1; + Vec3d b = p2 - a; + Vec3d c = p4 - a; + Vec3d d = p3 - a - b - c; - double dxb = d.X()*b.Y()-d.Y()*b.X(); - double dxc = d.X()*c.Y()-d.Y()*c.X(); - double dxa = d.X()*a.Y()-d.Y()*a.X(); - double dxp = d.X()*p.Y()-d.Y()*p.X(); + double dxb = d.X()*b.Y()-d.Y()*b.X(); + double dxc = d.X()*c.Y()-d.Y()*c.X(); + double dxa = d.X()*a.Y()-d.Y()*a.X(); + double dxp = d.X()*p.Y()-d.Y()*p.X(); - double c0,c1,c2,rt; - lami[2]=0.; - double eps = 1.E-12; + double c0,c1,c2; // ,rt; + lami[2]=0.; + double eps = 1.E-12; - if(fabs(d.X()) <= eps && fabs(d.Y())<= eps) - { + if(fabs(d.X()) <= eps && fabs(d.Y())<= eps) + { //Solve Linear System lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/ - (b.X()*c.Y() -b.Y()*c.X()); + (b.X()*c.Y() -b.Y()*c.X()); lami[1]=(-b.Y()*(p.X()-a.X())+b.X()*(p.Y()-a.Y()))/ - (b.X()*c.Y() -b.Y()*c.X()); - } - else - if(fabs(dxb) <= eps) + (b.X()*c.Y() -b.Y()*c.X()); + } + else + if(fabs(dxb) <= eps) { - lami[1] = (dxp-dxa)/dxc; - if(fabs(b.X()-d.X()*lami[1])>=eps) - lami[0] = (p.X()-a.X() - c.X()*lami[1])/(b.X()+d.X()*lami[1]); - else - lami[0] = (p.Y()-a.Y() - c.Y()*lami[1])/(b.Y()+d.Y()*lami[1]); + lami[1] = (dxp-dxa)/dxc; + if(fabs(b.X()-d.X()*lami[1])>=eps) + lami[0] = (p.X()-a.X() - c.X()*lami[1])/(b.X()+d.X()*lami[1]); + else + lami[0] = (p.Y()-a.Y() - c.Y()*lami[1])/(b.Y()+d.Y()*lami[1]); } + else + if(fabs(dxc) <= eps) + { + lami[0] = (dxp-dxa)/dxb; + if(fabs(c.X()-d.X()*lami[0])>=eps) + lami[1] = (p.X()-a.X() - b.X()*lami[0])/(c.X()+d.X()*lami[0]); + else + lami[1] = (p.Y()-a.Y() - b.Y()*lami[0])/(c.Y()+d.Y()*lami[0]); + } + else //Solve quadratic equation + { + if(fabs(d.X()) >= eps) + { + c2 = d.X()*dxc; + c1 = d.X()*dxc - c.X()*dxb - d.X()*(dxp-dxa); + c0 = -b.X()*(dxp -dxa) - (a.X()-p.X())*dxb; + } + else + { + c2 = d.Y()*dxc; + c1 = d.Y()*dxc - c.Y()*dxb - d.Y()*(dxp-dxa); + c0 = -b.Y()*(dxp -dxa) - (a.Y()-p.Y())*dxb; + } + + double rt = c1*c1 - 4*c2*c0; + if (rt < 0.) return false; + lami[1] = (-c1 + sqrt(rt))/2/c2; + if(lami[1]<=1. && lami[1]>=0.) + { + lami[0] = (dxp - dxa -dxc*lami[1])/dxb; + if(lami[0]<=1. && lami[0]>=0.) + return true; + } + + lami[1] = (-c1 - sqrt(rt))/2/c2; + lami[0] = (dxp - dxa -dxc*lami[1])/dxb; + } + + if( lami[0] <= 1.+eps && lami[0] >= -eps && lami[1]<=1.+eps && lami[1]>=-eps) + { + if(consider3D) + { + Vec3d n = Cross(b,c); + lami[2] = 0; + for(int i=1; i<=3; i++) + lami[2] +=(p.X(i)-a.X(i)-lami[0]*b.X(i)-lami[1]*c.X(i)) * n.X(i); + if(lami[2] >= -eps && lami[2] <= eps) + return true; + } else - if(fabs(dxc) <= eps) - { - lami[0] = (dxp-dxa)/dxb; - if(fabs(c.X()-d.X()*lami[0])>=eps) - lami[1] = (p.X()-a.X() - b.X()*lami[0])/(c.X()+d.X()*lami[0]); - else - lami[1] = (p.Y()-a.Y() - b.Y()*lami[0])/(c.Y()+d.Y()*lami[0]); - } - else //Solve quadratic equation - { - if(fabs(d.X()) >= eps) - { - c2 = d.X()*dxc; - c1 = d.X()*dxc - c.X()*dxb - d.X()*(dxp-dxa); - c0 = -b.X()*(dxp -dxa) - (a.X()-p.X())*dxb; - } - else - { - c2 = d.Y()*dxc; - c1 = d.Y()*dxc - c.Y()*dxb - d.Y()*(dxp-dxa); - c0 = -b.Y()*(dxp -dxa) - (a.Y()-p.Y())*dxb; - } + return true; + } - double rt = c1*c1 - 4*c2*c0; - if (rt < 0.) return false; - lami[1] = (-c1 + sqrt(rt))/2/c2; - if(lami[1]<=1. && lami[1]>=0.) - { - lami[0] = (dxp - dxa -dxc*lami[1])/dxb; - if(lami[0]<=1. && lami[0]>=0.) - return true; - } - - lami[1] = (-c1 - sqrt(rt))/2/c2; - lami[0] = (dxp - dxa -dxc*lami[1])/dxb; - } - - if( lami[0] <= 1.+eps && lami[0] >= -eps && lami[1]<=1.+eps && lami[1]>=-eps) - { - if(consider3D) - { - Vec3d n = Cross(b,c); - lami[2] = 0; - for(int i=1; i<=3; i++) - lami[2] +=(p.X(i)-a.X(i)-lami[0]*b.X(i)-lami[1]*c.X(i)) * n.X(i); - if(lami[2] >= -eps && lami[2] <= eps) - return true; - } - else - return true; - } - - return false; + return false; } - else + else { - // SurfaceElement(element).GetTets (loctets); - loctrigs.SetSize(1); - loctrigs.Elem(1) = SurfaceElement(element); + // SurfaceElement(element).GetTets (loctets); + loctrigs.SetSize(1); + loctrigs.Elem(1) = SurfaceElement(element); - for (int j = 1; j <= loctrigs.Size(); j++) - { + for (int j = 1; j <= loctrigs.Size(); j++) + { const Element2d & el = loctrigs.Get(j); @@ -4199,13 +4197,13 @@ namespace netgen const Point3d & p2 = Point(el.PNum(2)); const Point3d & p3 = Point(el.PNum(3)); /* - Box3d box; - box.SetPoint (p1); - box.AddPoint (p2); - box.AddPoint (p3); - box.AddPoint (p4); - if (!box.IsIn (p)) - continue; + Box3d box; + box.SetPoint (p1); + box.AddPoint (p2); + box.AddPoint (p3); + box.AddPoint (p4); + if (!box.IsIn (p)) + continue; */ col1 = p2-p1; col2 = p3-p1; @@ -4213,7 +4211,8 @@ namespace netgen //col3 = Vec3d(0, 0, 1); rhs = p - p1; - int retval = SolveLinearSystem (col1, col2, col3, rhs, sol); + // int retval = + SolveLinearSystem (col1, col2, col3, rhs, sol); //(*testout) << "retval " << retval << endl; @@ -4221,176 +4220,176 @@ namespace netgen //(*testout) << "sol " << sol << endl; if (sol.X() >= -eps && sol.Y() >= -eps && - sol.X() + sol.Y() <= 1+eps) - { - if(!consider3D || (sol.Z() >= -eps && sol.Z() <= eps)) - { - lami[0] = sol.X(); - lami[1] = sol.Y(); - lami[2] = sol.Z(); + sol.X() + sol.Y() <= 1+eps) + { + if(!consider3D || (sol.Z() >= -eps && sol.Z() <= eps)) + { + lami[0] = sol.X(); + lami[1] = sol.Y(); + lami[2] = sol.Z(); - return true; - } - } - } + return true; + } + } + } } + return false; + + } + + + + + bool Mesh :: PointContainedIn3DElement(const Point3d & p, + double lami[3], + const int element) const + { + //bool oldresult = PointContainedIn3DElementOld(p,lami,element); + //(*testout) << "old result: " << oldresult + // << " lam " << lami[0] << " " << lami[1] << " " << lami[2] << endl; + + //if(!curvedelems->IsElementCurved(element-1)) + // return PointContainedIn3DElementOld(p,lami,element); + + + const double eps = 1.e-4; + const Element & el = VolumeElement(element); + + netgen::Point<3> lam = 0.0; + + if (el.GetType() == TET) + { + lam = 0.25; + } + else if (el.GetType() == PRISM) + { + lam(0) = 0.33; lam(1) = 0.33; lam(2) = 0.5; + } + else if (el.GetType() == PYRAMID) + { + lam(0) = 0.4; lam(1) = 0.4; lam(2) = 0.2; + } + else if (el.GetType() == HEX) + { + lam = 0.5; + } + + + Vec<3> deltalam,rhs; + netgen::Point<3> x; + Mat<3,3> Jac,Jact; + + double delta=1; + + bool retval; + + int i = 0; + + const int maxits = 30; + + while(delta > 1e-16 && iCalcElementTransformation(lam,element-1,x,Jac); + + rhs = p-x; + Jac.Solve(rhs,deltalam); + + lam += deltalam; + + delta = deltalam.Length2(); + + i++; + //(*testout) << "pcie i " << i << " delta " << delta << " p " << p << " x " << x << " lam " << lam << endl; + //<< "Jac " << Jac << endl; + } + + if(i==maxits) return false; - } + + for(i=0; i<3; i++) + lami[i] = lam(i); - - bool Mesh :: PointContainedIn3DElement(const Point3d & p, - double lami[3], - const int element) const - { - //bool oldresult = PointContainedIn3DElementOld(p,lami,element); - //(*testout) << "old result: " << oldresult - // << " lam " << lami[0] << " " << lami[1] << " " << lami[2] << endl; - - //if(!curvedelems->IsElementCurved(element-1)) - // return PointContainedIn3DElementOld(p,lami,element); - - - const double eps = 1.e-4; - const Element & el = VolumeElement(element); - - netgen::Point<3> lam; - - if (el.GetType() == TET) + if (el.GetType() == TET) { - lam = 0.25; + retval = (lam(0) > -eps && + lam(1) > -eps && + lam(2) > -eps && + lam(0) + lam(1) + lam(2) < 1+eps); } - else if (el.GetType() == PRISM) + else if (el.GetType() == PRISM) { - lam(0) = 0.33; lam(1) = 0.33; lam(2) = 0.5; + retval = (lam(0) > -eps && + lam(1) > -eps && + lam(2) > -eps && + lam(2) < 1+eps && + lam(0) + lam(1) < 1+eps); } - else if (el.GetType() == PYRAMID) + else if (el.GetType() == PYRAMID) { - lam(0) = 0.4; lam(1) = 0.4; lam(2) = 0.2; + retval = (lam(0) > -eps && + lam(1) > -eps && + lam(2) > -eps && + lam(0) + lam(2) < 1+eps && + lam(1) + lam(2) < 1+eps); } - else if (el.GetType() == HEX) + else if (el.GetType() == HEX) { - lam = 0.5; + retval = (lam(0) > -eps && lam(0) < 1+eps && + lam(1) > -eps && lam(1) < 1+eps && + lam(2) > -eps && lam(2) < 1+eps); } + else + throw NgException("Da haun i wos vagessn"); + + return retval; + } - Vec<3> deltalam,rhs; - netgen::Point<3> x; - Mat<3,3> Jac,Jact; - double delta=1; + bool Mesh :: PointContainedIn3DElementOld(const Point3d & p, + double lami[3], + const int element) const + { - bool retval; + static Vec3d col1, col2, col3; + static Vec3d rhs, sol; + const double eps = 1.e-4; - int i = 0; + static Array loctets; - const int maxits = 30; + VolumeElement(element).GetTets (loctets); - while(delta > 1e-16 && iCalcElementTransformation(lam,element-1,x,Jac); + const Element & el = loctets.Get(j); - rhs = p-x; - Jac.Solve(rhs,deltalam); + const Point3d & p1 = Point(el.PNum(1)); + const Point3d & p2 = Point(el.PNum(2)); + const Point3d & p3 = Point(el.PNum(3)); + const Point3d & p4 = Point(el.PNum(4)); - lam += deltalam; + Box3d box; + box.SetPoint (p1); + box.AddPoint (p2); + box.AddPoint (p3); + box.AddPoint (p4); + if (!box.IsIn (p)) + continue; - delta = deltalam.Length2(); + col1 = p2-p1; + col2 = p3-p1; + col3 = p4-p1; + rhs = p - p1; - i++; - //(*testout) << "pcie i " << i << " delta " << delta << " p " << p << " x " << x << " lam " << lam << endl; - //<< "Jac " << Jac << endl; - } + SolveLinearSystem (col1, col2, col3, rhs, sol); - if(i==maxits) - return false; - - - for(i=0; i<3; i++) - lami[i] = lam(i); - - - - if (el.GetType() == TET) - { - retval = (lam(0) > -eps && - lam(1) > -eps && - lam(2) > -eps && - lam(0) + lam(1) + lam(2) < 1+eps); - } - else if (el.GetType() == PRISM) - { - retval = (lam(0) > -eps && - lam(1) > -eps && - lam(2) > -eps && - lam(2) < 1+eps && - lam(0) + lam(1) < 1+eps); - } - else if (el.GetType() == PYRAMID) - { - retval = (lam(0) > -eps && - lam(1) > -eps && - lam(2) > -eps && - lam(0) + lam(2) < 1+eps && - lam(1) + lam(2) < 1+eps); - } - else if (el.GetType() == HEX) - { - retval = (lam(0) > -eps && lam(0) < 1+eps && - lam(1) > -eps && lam(1) < 1+eps && - lam(2) > -eps && lam(2) < 1+eps); - } - else - throw NgException("Da haun i wos vagessn"); - - return retval; - } - - - - bool Mesh :: PointContainedIn3DElementOld(const Point3d & p, - double lami[3], - const int element) const - { - - static Vec3d col1, col2, col3; - static Vec3d rhs, sol; - const double eps = 1.e-4; - - static Array loctets; - - VolumeElement(element).GetTets (loctets); - - for (int j = 1; j <= loctets.Size(); j++) - { - const Element & el = loctets.Get(j); - - const Point3d & p1 = Point(el.PNum(1)); - const Point3d & p2 = Point(el.PNum(2)); - const Point3d & p3 = Point(el.PNum(3)); - const Point3d & p4 = Point(el.PNum(4)); - - Box3d box; - box.SetPoint (p1); - box.AddPoint (p2); - box.AddPoint (p3); - box.AddPoint (p4); - if (!box.IsIn (p)) - continue; - - col1 = p2-p1; - col2 = p3-p1; - col3 = p4-p1; - rhs = p - p1; - - SolveLinearSystem (col1, col2, col3, rhs, sol); - - if (sol.X() >= -eps && sol.Y() >= -eps && sol.Z() >= -eps && + if (sol.X() >= -eps && sol.Y() >= -eps && sol.Z() >= -eps && sol.X() + sol.Y() + sol.Z() <= 1+eps) - { + { Array loctetsloc; Array > pointsloc; @@ -4401,342 +4400,340 @@ namespace netgen Point3d pp = - pointsloc.Get(le.PNum(1)) - + sol.X() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(2))) - + sol.Y() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(3))) - + sol.Z() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(4))) ; + pointsloc.Get(le.PNum(1)) + + sol.X() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(2))) + + sol.Y() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(3))) + + sol.Z() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(4))) ; lami[0] = pp.X(); lami[1] = pp.Y(); lami[2] = pp.Z(); return true; - } + } } - return false; - } + return false; + } - int Mesh :: GetElementOfPoint (const Point3d & p, - double lami[3], - bool build_searchtree, - const int index, - const bool allowindex) const - { - if(index != -1) + int Mesh :: GetElementOfPoint (const Point3d & p, + double lami[3], + bool build_searchtree, + const int index, + const bool allowindex) const + { + if(index != -1) { - Array dummy(1); - dummy[0] = index; - return GetElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); + Array dummy(1); + dummy[0] = index; + return GetElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); } - else - return GetElementOfPoint(p,lami,NULL,build_searchtree,allowindex); - } + else + return GetElementOfPoint(p,lami,NULL,build_searchtree,allowindex); + } - int Mesh :: GetElementOfPoint (const Point3d & p, - double lami[3], - const Array * const indices, - bool build_searchtree, - const bool allowindex) const - { - if (dimension == 2) + int Mesh :: GetElementOfPoint (const Point3d & p, + double lami[3], + const Array * const indices, + bool build_searchtree, + const bool allowindex) const + { + if (dimension == 2) { - int i, j; - int ne; + int ne; - if(ps_startelement != 0 && ps_startelement <= GetNSE() && PointContainedIn2DElement(p,lami,ps_startelement)) - return ps_startelement; + if(ps_startelement != 0 && ps_startelement <= GetNSE() && PointContainedIn2DElement(p,lami,ps_startelement)) + return ps_startelement; - Array locels; - if (0) - { + Array locels; + if (0) + { elementsearchtree->GetIntersecting (p, p, locels); ne = locels.Size(); - } - else - ne = GetNSE(); + } + else + ne = GetNSE(); - for (i = 1; i <= ne; i++) - { + for (int i = 1; i <= ne; i++) + { int ii; if (0) - ii = locels.Get(i); + ii = locels.Get(i); else - ii = i; + ii = i; if(ii == ps_startelement) continue; if(indices != NULL && indices->Size() > 0) - { - bool contained = indices->Contains(SurfaceElement(ii).GetIndex()); - if((allowindex && !contained) || (!allowindex && contained)) continue; - } + { + bool contained = indices->Contains(SurfaceElement(ii).GetIndex()); + if((allowindex && !contained) || (!allowindex && contained)) continue; + } if(PointContainedIn2DElement(p,lami,ii)) return ii; - } - return 0; + } + return 0; } - else + else { - int i, j; - int ne; + // int i, j; + int ne; - if(ps_startelement != 0 && PointContainedIn3DElement(p,lami,ps_startelement)) - return ps_startelement; + if(ps_startelement != 0 && PointContainedIn3DElement(p,lami,ps_startelement)) + return ps_startelement; - Array locels; - if (elementsearchtree || build_searchtree) - { + Array locels; + if (elementsearchtree || build_searchtree) + { // update if necessary: const_cast(*this).BuildElementSearchTree (); elementsearchtree->GetIntersecting (p, p, locels); ne = locels.Size(); - } - else - ne = GetNE(); + } + else + ne = GetNE(); - for (i = 1; i <= ne; i++) - { + for (int i = 1; i <= ne; i++) + { int ii; if (elementsearchtree) - ii = locels.Get(i); + ii = locels.Get(i); else - ii = i; + ii = i; if(ii == ps_startelement) continue; if(indices != NULL && indices->Size() > 0) - { - bool contained = indices->Contains(VolumeElement(ii).GetIndex()); - if((allowindex && !contained) || (!allowindex && contained)) continue; - } - + { + bool contained = indices->Contains(VolumeElement(ii).GetIndex()); + if((allowindex && !contained) || (!allowindex && contained)) continue; + } if(PointContainedIn3DElement(p,lami,ii)) - { - ps_startelement = ii; - return ii; - } - } + { + ps_startelement = ii; + return ii; + } + } - // Not found, try uncurved variant: - for (i = 1; i <= ne; i++) - { + // Not found, try uncurved variant: + for (int i = 1; i <= ne; i++) + { int ii; if (elementsearchtree) - ii = locels.Get(i); + ii = locels.Get(i); else - ii = i; + ii = i; if(indices != NULL && indices->Size() > 0) - { - bool contained = indices->Contains(VolumeElement(ii).GetIndex()); - if((allowindex && !contained) || (!allowindex && contained)) continue; - } + { + bool contained = indices->Contains(VolumeElement(ii).GetIndex()); + if((allowindex && !contained) || (!allowindex && contained)) continue; + } if(PointContainedIn3DElementOld(p,lami,ii)) - { - ps_startelement = ii; - (*testout) << "WARNING: found element of point " << p <<" only for uncurved mesh" << endl; - return ii; - } - } + { + ps_startelement = ii; + (*testout) << "WARNING: found element of point " << p <<" only for uncurved mesh" << endl; + return ii; + } + } - return 0; + return 0; } - } + } - int Mesh :: GetSurfaceElementOfPoint (const Point3d & p, - double lami[3], - bool build_searchtree, - const int index, - const bool allowindex) const - { - if(index != -1) + int Mesh :: GetSurfaceElementOfPoint (const Point3d & p, + double lami[3], + bool build_searchtree, + const int index, + const bool allowindex) const + { + if(index != -1) { - Array dummy(1); - dummy[0] = index; - return GetSurfaceElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); + Array dummy(1); + dummy[0] = index; + return GetSurfaceElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); } - else - return GetSurfaceElementOfPoint(p,lami,NULL,build_searchtree,allowindex); - } + else + return GetSurfaceElementOfPoint(p,lami,NULL,build_searchtree,allowindex); + } - int Mesh :: GetSurfaceElementOfPoint (const Point3d & p, - double lami[3], - const Array * const indices, - bool build_searchtree, - const bool allowindex) const - { - if (dimension == 2) + int Mesh :: GetSurfaceElementOfPoint (const Point3d & p, + double lami[3], + const Array * const indices, + bool build_searchtree, + const bool allowindex) const + { + if (dimension == 2) { - throw NgException("GetSurfaceElementOfPoint not yet implemented for 2D meshes"); + throw NgException("GetSurfaceElementOfPoint not yet implemented for 2D meshes"); } - else + else { - double vlam[3]; - int velement = GetElementOfPoint(p,vlam,NULL,build_searchtree,allowindex); + double vlam[3]; + int velement = GetElementOfPoint(p,vlam,NULL,build_searchtree,allowindex); - //(*testout) << "p " << p << endl; - //(*testout) << "velement " << velement << endl; + //(*testout) << "p " << p << endl; + //(*testout) << "velement " << velement << endl; - Array faces; - topology->GetElementFaces(velement,faces); + Array faces; + topology->GetElementFaces(velement,faces); - //(*testout) << "faces " << faces << endl; + //(*testout) << "faces " << faces << endl; - for(int i=0; iGetFace2SurfaceElement(faces[i]); + for(int i=0; iGetFace2SurfaceElement(faces[i]); - //(*testout) << "surfel " << faces << endl; + //(*testout) << "surfel " << faces << endl; - for(int i=0; iSize() != 0) - { - if(indices->Contains(SurfaceElement(faces[i]).GetIndex()) && - PointContainedIn2DElement(p,lami,faces[i],true)) + { + if(indices->Contains(SurfaceElement(faces[i]).GetIndex()) && + PointContainedIn2DElement(p,lami,faces[i],true)) return faces[i]; - } + } else - { - if(PointContainedIn2DElement(p,lami,faces[i],true)) - { - //(*testout) << "found point " << p << " in sel " << faces[i] - // << ", lam " << lami[0] << ", " << lami[1] << ", " << lami[2] << endl; - return faces[i]; - } - } - } + { + if(PointContainedIn2DElement(p,lami,faces[i],true)) + { + //(*testout) << "found point " << p << " in sel " << faces[i] + // << ", lam " << lami[0] << ", " << lami[1] << ", " << lami[2] << endl; + return faces[i]; + } + } + } } - return 0; - } + return 0; + } - void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, - Array & locels) const - { - elementsearchtree->GetIntersecting (p1, p2, locels); - } + void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, + Array & locels) const + { + elementsearchtree->GetIntersecting (p1, p2, locels); + } - void Mesh :: SplitIntoParts() - { - int i, j, dom; - int ne = GetNE(); - int np = GetNP(); - int nse = GetNSE(); + void Mesh :: SplitIntoParts() + { + int i, j, dom; + int ne = GetNE(); + int np = GetNP(); + int nse = GetNSE(); - BitArray surfused(nse); - BitArray pused (np); + BitArray surfused(nse); + BitArray pused (np); - surfused.Clear(); + surfused.Clear(); - dom = 0; + dom = 0; - while (1) + while (1) { - int cntd = 1; + int cntd = 1; - dom++; + dom++; - pused.Clear(); + pused.Clear(); - int found = 0; - for (i = 1; i <= nse; i++) - if (!surfused.Test(i)) + int found = 0; + for (i = 1; i <= nse; i++) + if (!surfused.Test(i)) { - SurfaceElement(i).SetIndex (dom); - for (j = 1; j <= 3; j++) - pused.Set (SurfaceElement(i).PNum(j)); - found = 1; - cntd = 1; - surfused.Set(i); - break; + SurfaceElement(i).SetIndex (dom); + for (j = 1; j <= 3; j++) + pused.Set (SurfaceElement(i).PNum(j)); + found = 1; + cntd = 1; + surfused.Set(i); + break; } - if (!found) - break; + if (!found) + break; - int change; - do - { - change = 0; - for (i = 1; i <= nse; i++) - { - int is = 0, isnot = 0; - for (j = 1; j <= 3; j++) - if (pused.Test(SurfaceElement(i).PNum(j))) - is = 1; - else - isnot = 1; + int change; + do + { + change = 0; + for (i = 1; i <= nse; i++) + { + int is = 0, isnot = 0; + for (j = 1; j <= 3; j++) + if (pused.Test(SurfaceElement(i).PNum(j))) + is = 1; + else + isnot = 1; - if (is && isnot) + if (is && isnot) { - change = 1; - for (j = 1; j <= 3; j++) - pused.Set (SurfaceElement(i).PNum(j)); + change = 1; + for (j = 1; j <= 3; j++) + pused.Set (SurfaceElement(i).PNum(j)); } - if (is) + if (is) { - if (!surfused.Test(i)) - { + if (!surfused.Test(i)) + { surfused.Set(i); SurfaceElement(i).SetIndex (dom); cntd++; - } + } } - } + } - for (i = 1; i <= ne; i++) - { - int is = 0, isnot = 0; - for (j = 1; j <= 4; j++) - if (pused.Test(VolumeElement(i).PNum(j))) - is = 1; - else - isnot = 1; + for (i = 1; i <= ne; i++) + { + int is = 0, isnot = 0; + for (j = 1; j <= 4; j++) + if (pused.Test(VolumeElement(i).PNum(j))) + is = 1; + else + isnot = 1; - if (is && isnot) + if (is && isnot) { - change = 1; - for (j = 1; j <= 4; j++) - pused.Set (VolumeElement(i).PNum(j)); + change = 1; + for (j = 1; j <= 4; j++) + pused.Set (VolumeElement(i).PNum(j)); } - if (is) + if (is) { - VolumeElement(i).SetIndex (dom); + VolumeElement(i).SetIndex (dom); } - } - } - while (change); + } + } + while (change); - PrintMessage (3, "domain ", dom, " has ", cntd, " surfaceelements"); + PrintMessage (3, "domain ", dom, " has ", cntd, " surfaceelements"); } - /* + /* facedecoding.SetSize (dom); for (i = 1; i <= dom; i++) { @@ -4744,106 +4741,106 @@ namespace netgen facedecoding.Elem(i).domin = i; facedecoding.Elem(i).domout = 0; } - */ - ClearFaceDescriptors(); - for (i = 1; i <= dom; i++) - AddFaceDescriptor (FaceDescriptor (0, i, 0, 0)); - CalcSurfacesOfNode(); - timestamp = NextTimeStamp(); - } + */ + ClearFaceDescriptors(); + for (i = 1; i <= dom; i++) + AddFaceDescriptor (FaceDescriptor (0, i, 0, 0)); + CalcSurfacesOfNode(); + timestamp = NextTimeStamp(); + } - void Mesh :: SplitSeparatedFaces () - { - PrintMessage (3, "SplitSeparateFaces"); - int fdi; - int np = GetNP(); + void Mesh :: SplitSeparatedFaces () + { + PrintMessage (3, "SplitSeparateFaces"); + int fdi; + int np = GetNP(); - BitArray usedp(np); - Array els_of_face; + BitArray usedp(np); + Array els_of_face; - fdi = 1; - while (fdi <= GetNFD()) + fdi = 1; + while (fdi <= GetNFD()) { - GetSurfaceElementsOfFace (fdi, els_of_face); + GetSurfaceElementsOfFace (fdi, els_of_face); - if (els_of_face.Size() == 0) continue; + if (els_of_face.Size() == 0) continue; - SurfaceElementIndex firstel = els_of_face[0]; + SurfaceElementIndex firstel = els_of_face[0]; - usedp.Clear(); - for (int j = 1; j <= SurfaceElement(firstel).GetNP(); j++) - usedp.Set (SurfaceElement(firstel).PNum(j)); + usedp.Clear(); + for (int j = 1; j <= SurfaceElement(firstel).GetNP(); j++) + usedp.Set (SurfaceElement(firstel).PNum(j)); - bool changed; - do - { + bool changed; + do + { changed = false; for (int i = 0; i < els_of_face.Size(); i++) - { - const Element2d & el = SurfaceElement(els_of_face[i]); + { + const Element2d & el = SurfaceElement(els_of_face[i]); - bool has = 0; - bool hasno = 0; - for (int j = 0; j < el.GetNP(); j++) - { - if (usedp.Test(el[j])) - has = true; - else - hasno = true; - } + bool has = 0; + bool hasno = 0; + for (int j = 0; j < el.GetNP(); j++) + { + if (usedp.Test(el[j])) + has = true; + else + hasno = true; + } - if (has && hasno) + if (has && hasno) changed = true; - if (has) + if (has) for (int j = 0; j < el.GetNP(); j++) - usedp.Set (el[j]); - } - } - while (changed); + usedp.Set (el[j]); + } + } + while (changed); - int nface = 0; - for (int i = 0; i < els_of_face.Size(); i++) - { + int nface = 0; + for (int i = 0; i < els_of_face.Size(); i++) + { Element2d & el = SurfaceElement(els_of_face[i]); int hasno = 0; for (int j = 1; j <= el.GetNP(); j++) - if (!usedp.Test(el.PNum(j))) - hasno = 1; + if (!usedp.Test(el.PNum(j))) + hasno = 1; if (hasno) - { - if (!nface) - { - FaceDescriptor nfd = GetFaceDescriptor(fdi); - nface = AddFaceDescriptor (nfd); - } + { + if (!nface) + { + FaceDescriptor nfd = GetFaceDescriptor(fdi); + nface = AddFaceDescriptor (nfd); + } - el.SetIndex (nface); - } - } + el.SetIndex (nface); + } + } - // reconnect list - if (nface) - { + // reconnect list + if (nface) + { facedecoding[nface-1].firstelement = -1; facedecoding[fdi-1].firstelement = -1; for (int i = 0; i < els_of_face.Size(); i++) - { - int ind = SurfaceElement(els_of_face[i]).GetIndex(); - SurfaceElement(els_of_face[i]).next = facedecoding[ind-1].firstelement; - facedecoding[ind-1].firstelement = els_of_face[i]; - } - } + { + int ind = SurfaceElement(els_of_face[i]).GetIndex(); + SurfaceElement(els_of_face[i]).next = facedecoding[ind-1].firstelement; + facedecoding[ind-1].firstelement = els_of_face[i]; + } + } - fdi++; + fdi++; } - /* + /* fdi = 1; while (fdi <= GetNFD()) { @@ -4916,16 +4913,16 @@ namespace netgen } fdi++; } - */ - } + */ + } - void Mesh :: GetSurfaceElementsOfFace (int facenr, Array & sei) const - { - static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace"); - NgProfiler::RegionTimer reg (timer); + void Mesh :: GetSurfaceElementsOfFace (int facenr, Array & sei) const + { + static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace"); + NgProfiler::RegionTimer reg (timer); - /* + /* sei.SetSize (0); for (SurfaceElementIndex i = 0; i < GetNSE(); i++) if ( (*this)[i].GetIndex () == facenr && (*this)[i][0] >= PointIndex::BASE && @@ -4933,614 +4930,614 @@ namespace netgen sei.Append (i); int size1 = sei.Size(); - */ + */ - sei.SetSize(0); + sei.SetSize(0); - SurfaceElementIndex si = facedecoding[facenr-1].firstelement; - while (si != -1) + SurfaceElementIndex si = facedecoding[facenr-1].firstelement; + while (si != -1) { - if ( (*this)[si].GetIndex () == facenr && (*this)[si][0] >= PointIndex::BASE && - !(*this)[si].IsDeleted() ) - { + if ( (*this)[si].GetIndex () == facenr && (*this)[si][0] >= PointIndex::BASE && + !(*this)[si].IsDeleted() ) + { sei.Append (si); - } + } - si = (*this)[si].next; + si = (*this)[si].next; } - /* - // *testout << "with list = " << endl << sei << endl; + /* + // *testout << "with list = " << endl << sei << endl; - if (size1 != sei.Size()) + if (size1 != sei.Size()) + { + cout << "size mismatch" << endl; + exit(1); + } + */ + } + + + + + void Mesh :: CalcMinMaxAngle (double badellimit, double * retvalues) + { + int i, j; + int lpi1, lpi2, lpi3, lpi4; + double phimax = 0, phimin = 10; + double facephimax = 0, facephimin = 10; + int illegaltets = 0, negativetets = 0, badtets = 0; + + for (i = 1; i <= GetNE(); i++) { - cout << "size mismatch" << endl; - exit(1); - } - */ - } + int badel = 0; + Element & el = VolumeElement(i); - - - void Mesh :: CalcMinMaxAngle (double badellimit, double * retvalues) - { - int i, j; - int lpi1, lpi2, lpi3, lpi4; - double phimax = 0, phimin = 10; - double facephimax = 0, facephimin = 10; - int illegaltets = 0, negativetets = 0, badtets = 0; - - for (i = 1; i <= GetNE(); i++) - { - int badel = 0; - - Element & el = VolumeElement(i); - - if (el.GetType() != TET) - { + if (el.GetType() != TET) + { VolumeElement(i).flags.badel = 0; continue; - } + } - if (el.Volume(Points()) < 0) - { + if (el.Volume(Points()) < 0) + { badel = 1; negativetets++; - } + } - if (!LegalTet (el)) - { + if (!LegalTet (el)) + { badel = 1; illegaltets++; (*testout) << "illegal tet: " << i << " "; for (j = 1; j <= el.GetNP(); j++) - (*testout) << el.PNum(j) << " "; + (*testout) << el.PNum(j) << " "; (*testout) << endl; - } + } - // angles between faces - for (lpi1 = 1; lpi1 <= 3; lpi1++) - for (lpi2 = lpi1+1; lpi2 <= 4; lpi2++) + // angles between faces + for (lpi1 = 1; lpi1 <= 3; lpi1++) + for (lpi2 = lpi1+1; lpi2 <= 4; lpi2++) { - lpi3 = 1; - while (lpi3 == lpi1 || lpi3 == lpi2) - lpi3++; - lpi4 = 10 - lpi1 - lpi2 - lpi3; + lpi3 = 1; + while (lpi3 == lpi1 || lpi3 == lpi2) + lpi3++; + lpi4 = 10 - lpi1 - lpi2 - lpi3; - const Point3d & p1 = Point (el.PNum(lpi1)); - const Point3d & p2 = Point (el.PNum(lpi2)); - const Point3d & p3 = Point (el.PNum(lpi3)); - const Point3d & p4 = Point (el.PNum(lpi4)); + const Point3d & p1 = Point (el.PNum(lpi1)); + const Point3d & p2 = Point (el.PNum(lpi2)); + const Point3d & p3 = Point (el.PNum(lpi3)); + const Point3d & p4 = Point (el.PNum(lpi4)); - Vec3d n(p1, p2); - n /= n.Length(); - Vec3d v1(p1, p3); - Vec3d v2(p1, p4); + Vec3d n(p1, p2); + n /= n.Length(); + Vec3d v1(p1, p3); + Vec3d v2(p1, p4); - v1 -= (n * v1) * n; - v2 -= (n * v2) * n; + v1 -= (n * v1) * n; + v2 -= (n * v2) * n; - double cosphi = (v1 * v2) / (v1.Length() * v2.Length()); - double phi = acos (cosphi); - if (phi > phimax) phimax = phi; - if (phi < phimin) phimin = phi; + double cosphi = (v1 * v2) / (v1.Length() * v2.Length()); + double phi = acos (cosphi); + if (phi > phimax) phimax = phi; + if (phi < phimin) phimin = phi; - if ((180/M_PI) * phi > badellimit) + if ((180/M_PI) * phi > badellimit) + badel = 1; + } + + + // angles in faces + for (j = 1; j <= 4; j++) + { + Element2d face; + el.GetFace (j, face); + for (lpi1 = 1; lpi1 <= 3; lpi1++) + { + lpi2 = lpi1 % 3 + 1; + lpi3 = lpi2 % 3 + 1; + + const Point3d & p1 = Point (el.PNum(lpi1)); + const Point3d & p2 = Point (el.PNum(lpi2)); + const Point3d & p3 = Point (el.PNum(lpi3)); + + Vec3d v1(p1, p2); + Vec3d v2(p1, p3); + double cosphi = (v1 * v2) / (v1.Length() * v2.Length()); + double phi = acos (cosphi); + if (phi > facephimax) facephimax = phi; + if (phi < facephimin) facephimin = phi; + + if ((180/M_PI) * phi > badellimit) badel = 1; - } + + } + } - // angles in faces - for (j = 1; j <= 4; j++) - { - Element2d face; - el.GetFace (j, face); - for (lpi1 = 1; lpi1 <= 3; lpi1++) - { - lpi2 = lpi1 % 3 + 1; - lpi3 = lpi2 % 3 + 1; - - const Point3d & p1 = Point (el.PNum(lpi1)); - const Point3d & p2 = Point (el.PNum(lpi2)); - const Point3d & p3 = Point (el.PNum(lpi3)); - - Vec3d v1(p1, p2); - Vec3d v2(p1, p3); - double cosphi = (v1 * v2) / (v1.Length() * v2.Length()); - double phi = acos (cosphi); - if (phi > facephimax) facephimax = phi; - if (phi < facephimin) facephimin = phi; - - if ((180/M_PI) * phi > badellimit) - badel = 1; - - } - } - - - VolumeElement(i).flags.badel = badel; - if (badel) badtets++; + VolumeElement(i).flags.badel = badel; + if (badel) badtets++; } - if (!GetNE()) + if (!GetNE()) { - phimin = phimax = facephimin = facephimax = 0; + phimin = phimax = facephimin = facephimax = 0; } - if (!retvalues) + if (!retvalues) { - PrintMessage (1, ""); - PrintMessage (1, "between planes: phimin = ", (180/M_PI) * phimin, - " phimax = ", (180/M_PI) *phimax); - PrintMessage (1, "inside planes: phimin = ", (180/M_PI) * facephimin, - " phimax = ", (180/M_PI) * facephimax); - PrintMessage (1, ""); + PrintMessage (1, ""); + PrintMessage (1, "between planes: phimin = ", (180/M_PI) * phimin, + " phimax = ", (180/M_PI) *phimax); + PrintMessage (1, "inside planes: phimin = ", (180/M_PI) * facephimin, + " phimax = ", (180/M_PI) * facephimax); + PrintMessage (1, ""); } - else + else { - retvalues[0] = (180/M_PI) * facephimin; - retvalues[1] = (180/M_PI) * facephimax; - retvalues[2] = (180/M_PI) * phimin; - retvalues[3] = (180/M_PI) * phimax; + retvalues[0] = (180/M_PI) * facephimin; + retvalues[1] = (180/M_PI) * facephimax; + retvalues[2] = (180/M_PI) * phimin; + retvalues[3] = (180/M_PI) * phimax; } - PrintMessage (3, "negative tets: ", negativetets); - PrintMessage (3, "illegal tets: ", illegaltets); - PrintMessage (3, "bad tets: ", badtets); - } + PrintMessage (3, "negative tets: ", negativetets); + PrintMessage (3, "illegal tets: ", illegaltets); + PrintMessage (3, "bad tets: ", badtets); + } - int Mesh :: MarkIllegalElements () - { - int cnt = 0; - int i; + int Mesh :: MarkIllegalElements () + { + int cnt = 0; + int i; - for (i = 1; i <= GetNE(); i++) + for (i = 1; i <= GetNE(); i++) { - LegalTet (VolumeElement(i)); + LegalTet (VolumeElement(i)); - /* - Element & el = VolumeElement(i); - int leg1 = LegalTet (el); - el.flags.illegal_valid = 0; - int leg2 = LegalTet (el); + /* + Element & el = VolumeElement(i); + int leg1 = LegalTet (el); + el.flags.illegal_valid = 0; + int leg2 = LegalTet (el); - if (leg1 != leg2) - { - cerr << "legal differs!!" << endl; - (*testout) << "legal differs" << endl; - (*testout) << "elnr = " << i << ", el = " << el - << " leg1 = " << leg1 << ", leg2 = " << leg2 << endl; - } + if (leg1 != leg2) + { + cerr << "legal differs!!" << endl; + (*testout) << "legal differs" << endl; + (*testout) << "elnr = " << i << ", el = " << el + << " leg1 = " << leg1 << ", leg2 = " << leg2 << endl; + } - // el.flags.illegal = !LegalTet (el); - */ - cnt += VolumeElement(i).Illegal(); + // el.flags.illegal = !LegalTet (el); + */ + cnt += VolumeElement(i).Illegal(); } - return cnt; - } + return cnt; + } - // #ifdef NONE - // void Mesh :: AddIdentification (int pi1, int pi2, int identnr) - // { - // INDEX_2 pair(pi1, pi2); - // // pair.Sort(); - // identifiedpoints->Set (pair, identnr); - // if (identnr > maxidentnr) - // maxidentnr = identnr; - // timestamp = NextTimeStamp(); - // } + // #ifdef NONE + // void Mesh :: AddIdentification (int pi1, int pi2, int identnr) + // { + // INDEX_2 pair(pi1, pi2); + // // pair.Sort(); + // identifiedpoints->Set (pair, identnr); + // if (identnr > maxidentnr) + // maxidentnr = identnr; + // timestamp = NextTimeStamp(); + // } - // int Mesh :: GetIdentification (int pi1, int pi2) const - // { - // INDEX_2 pair(pi1, pi2); - // if (identifiedpoints->Used (pair)) - // return identifiedpoints->Get(pair); - // else - // return 0; - // } + // int Mesh :: GetIdentification (int pi1, int pi2) const + // { + // INDEX_2 pair(pi1, pi2); + // if (identifiedpoints->Used (pair)) + // return identifiedpoints->Get(pair); + // else + // return 0; + // } - // int Mesh :: GetIdentificationSym (int pi1, int pi2) const - // { - // INDEX_2 pair(pi1, pi2); - // if (identifiedpoints->Used (pair)) - // return identifiedpoints->Get(pair); + // int Mesh :: GetIdentificationSym (int pi1, int pi2) const + // { + // INDEX_2 pair(pi1, pi2); + // if (identifiedpoints->Used (pair)) + // return identifiedpoints->Get(pair); - // pair = INDEX_2 (pi2, pi1); - // if (identifiedpoints->Used (pair)) - // return identifiedpoints->Get(pair); + // pair = INDEX_2 (pi2, pi1); + // if (identifiedpoints->Used (pair)) + // return identifiedpoints->Get(pair); - // return 0; - // } + // return 0; + // } - // void Mesh :: GetIdentificationMap (int identnr, Array & identmap) const - // { - // int i, j; + // void Mesh :: GetIdentificationMap (int identnr, Array & identmap) const + // { + // int i, j; - // identmap.SetSize (GetNP()); - // for (i = 1; i <= identmap.Size(); i++) - // identmap.Elem(i) = 0; + // identmap.SetSize (GetNP()); + // for (i = 1; i <= identmap.Size(); i++) + // identmap.Elem(i) = 0; - // for (i = 1; i <= identifiedpoints->GetNBags(); i++) - // for (j = 1; j <= identifiedpoints->GetBagSize(i); j++) - // { - // INDEX_2 i2; - // int nr; - // identifiedpoints->GetData (i, j, i2, nr); + // for (i = 1; i <= identifiedpoints->GetNBags(); i++) + // for (j = 1; j <= identifiedpoints->GetBagSize(i); j++) + // { + // INDEX_2 i2; + // int nr; + // identifiedpoints->GetData (i, j, i2, nr); - // if (nr == identnr) - // { - // identmap.Elem(i2.I1()) = i2.I2(); - // } - // } - // } + // if (nr == identnr) + // { + // identmap.Elem(i2.I1()) = i2.I2(); + // } + // } + // } - // void Mesh :: GetIdentificationPairs (int identnr, Array & identpairs) const - // { - // int i, j; + // void Mesh :: GetIdentificationPairs (int identnr, Array & identpairs) const + // { + // int i, j; - // identpairs.SetSize(0); + // identpairs.SetSize(0); - // for (i = 1; i <= identifiedpoints->GetNBags(); i++) - // for (j = 1; j <= identifiedpoints->GetBagSize(i); j++) - // { - // INDEX_2 i2; - // int nr; - // identifiedpoints->GetData (i, j, i2, nr); + // for (i = 1; i <= identifiedpoints->GetNBags(); i++) + // for (j = 1; j <= identifiedpoints->GetBagSize(i); j++) + // { + // INDEX_2 i2; + // int nr; + // identifiedpoints->GetData (i, j, i2, nr); - // if (identnr == 0 || nr == identnr) - // identpairs.Append (i2); - // } - // } - // #endif + // if (identnr == 0 || nr == identnr) + // identpairs.Append (i2); + // } + // } + // #endif - void Mesh :: InitPointCurve(double red, double green, double blue) const - { - pointcurves_startpoint.Append(pointcurves.Size()); - pointcurves_red.Append(red); - pointcurves_green.Append(green); - pointcurves_blue.Append(blue); - } - void Mesh :: AddPointCurvePoint(const Point3d & pt) const - { - pointcurves.Append(pt); - } - int Mesh :: GetNumPointCurves(void) const - { - return pointcurves_startpoint.Size(); - } - int Mesh :: GetNumPointsOfPointCurve(int curve) const - { - if(curve == pointcurves_startpoint.Size()-1) - return (pointcurves.Size() - pointcurves_startpoint.Last()); - else - return (pointcurves_startpoint[curve+1]-pointcurves_startpoint[curve]); - } + void Mesh :: InitPointCurve(double red, double green, double blue) const + { + pointcurves_startpoint.Append(pointcurves.Size()); + pointcurves_red.Append(red); + pointcurves_green.Append(green); + pointcurves_blue.Append(blue); + } + void Mesh :: AddPointCurvePoint(const Point3d & pt) const + { + pointcurves.Append(pt); + } + int Mesh :: GetNumPointCurves(void) const + { + return pointcurves_startpoint.Size(); + } + int Mesh :: GetNumPointsOfPointCurve(int curve) const + { + if(curve == pointcurves_startpoint.Size()-1) + return (pointcurves.Size() - pointcurves_startpoint.Last()); + else + return (pointcurves_startpoint[curve+1]-pointcurves_startpoint[curve]); + } - Point3d & Mesh :: GetPointCurvePoint(int curve, int n) const - { - return pointcurves[pointcurves_startpoint[curve]+n]; - } + Point3d & Mesh :: GetPointCurvePoint(int curve, int n) const + { + return pointcurves[pointcurves_startpoint[curve]+n]; + } - void Mesh :: GetPointCurveColor(int curve, double & red, double & green, double & blue) const - { - red = pointcurves_red[curve]; - green = pointcurves_green[curve]; - blue = pointcurves_blue[curve]; - } + void Mesh :: GetPointCurveColor(int curve, double & red, double & green, double & blue) const + { + red = pointcurves_red[curve]; + green = pointcurves_green[curve]; + blue = pointcurves_blue[curve]; + } - void Mesh :: ComputeNVertices () - { - int i, j, nv; - int ne = GetNE(); - int nse = GetNSE(); + void Mesh :: ComputeNVertices () + { + int i, j, nv; + int ne = GetNE(); + int nse = GetNSE(); - numvertices = 0; - for (i = 1; i <= ne; i++) + numvertices = 0; + for (i = 1; i <= ne; i++) { - const Element & el = VolumeElement(i); - nv = el.GetNV(); - for (j = 0; j < nv; j++) - if (el[j] > numvertices) - numvertices = el[j]; + const Element & el = VolumeElement(i); + nv = el.GetNV(); + for (j = 0; j < nv; j++) + if (el[j] > numvertices) + numvertices = el[j]; } - for (i = 1; i <= nse; i++) + for (i = 1; i <= nse; i++) { - const Element2d & el = SurfaceElement(i); - nv = el.GetNV(); - for (j = 1; j <= nv; j++) - if (el.PNum(j) > numvertices) - numvertices = el.PNum(j); + const Element2d & el = SurfaceElement(i); + nv = el.GetNV(); + for (j = 1; j <= nv; j++) + if (el.PNum(j) > numvertices) + numvertices = el.PNum(j); } - numvertices += 1- PointIndex::BASE; - } + numvertices += 1- PointIndex::BASE; + } - int Mesh :: GetNV () const - { - if (numvertices < 0) - return GetNP(); - else - return numvertices; - } + int Mesh :: GetNV () const + { + if (numvertices < 0) + return GetNP(); + else + return numvertices; + } - void Mesh :: SetNP (int np) - { - points.SetSize(np); - // ptyps.SetSize(np); + void Mesh :: SetNP (int np) + { + points.SetSize(np); + // ptyps.SetSize(np); - int mlold = mlbetweennodes.Size(); - mlbetweennodes.SetSize(np); - if (np > mlold) - for (int i = mlold+PointIndex::BASE; - i < np+PointIndex::BASE; i++) - { - mlbetweennodes[i].I1() = PointIndex::BASE-1; - mlbetweennodes[i].I2() = PointIndex::BASE-1; - } + int mlold = mlbetweennodes.Size(); + mlbetweennodes.SetSize(np); + if (np > mlold) + for (int i = mlold+PointIndex::BASE; + i < np+PointIndex::BASE; i++) + { + mlbetweennodes[i].I1() = PointIndex::BASE-1; + mlbetweennodes[i].I2() = PointIndex::BASE-1; + } - GetIdentifications().SetMaxPointNr (np + PointIndex::BASE-1); - } + GetIdentifications().SetMaxPointNr (np + PointIndex::BASE-1); + } - /* - void Mesh :: BuildConnectedNodes () - { - if (PureTetMesh()) - { - connectedtonode.SetSize(0); - return; - } + /* + void Mesh :: BuildConnectedNodes () + { + if (PureTetMesh()) + { + connectedtonode.SetSize(0); + return; + } - int i, j, k; - int np = GetNP(); - int ne = GetNE(); - TABLE conto(np); - for (i = 1; i <= ne; i++) - { - const Element & el = VolumeElement(i); + int i, j, k; + int np = GetNP(); + int ne = GetNE(); + TABLE conto(np); + for (i = 1; i <= ne; i++) + { + const Element & el = VolumeElement(i); - if (el.GetType() == PRISM) - { - for (j = 1; j <= 6; j++) - { - int n1 = el.PNum (j); - int n2 = el.PNum ((j+2)%6+1); - // if (n1 != n2) - { - int found = 0; - for (k = 1; k <= conto.EntrySize(n1); k++) - if (conto.Get(n1, k) == n2) - { - found = 1; - break; - } - if (!found) - conto.Add (n1, n2); - } - } - } - else if (el.GetType() == PYRAMID) - { - for (j = 1; j <= 4; j++) - { - int n1, n2; - switch (j) - { - case 1: n1 = 1; n2 = 4; break; - case 2: n1 = 4; n2 = 1; break; - case 3: n1 = 2; n2 = 3; break; - case 4: n1 = 3; n2 = 2; break; - } + if (el.GetType() == PRISM) + { + for (j = 1; j <= 6; j++) + { + int n1 = el.PNum (j); + int n2 = el.PNum ((j+2)%6+1); + // if (n1 != n2) + { + int found = 0; + for (k = 1; k <= conto.EntrySize(n1); k++) + if (conto.Get(n1, k) == n2) + { + found = 1; + break; + } + if (!found) + conto.Add (n1, n2); + } + } + } + else if (el.GetType() == PYRAMID) + { + for (j = 1; j <= 4; j++) + { + int n1, n2; + switch (j) + { + case 1: n1 = 1; n2 = 4; break; + case 2: n1 = 4; n2 = 1; break; + case 3: n1 = 2; n2 = 3; break; + case 4: n1 = 3; n2 = 2; break; + } - int found = 0; - for (k = 1; k <= conto.EntrySize(n1); k++) - if (conto.Get(n1, k) == n2) - { - found = 1; - break; - } - if (!found) - conto.Add (n1, n2); - } - } - } + int found = 0; + for (k = 1; k <= conto.EntrySize(n1); k++) + if (conto.Get(n1, k) == n2) + { + found = 1; + break; + } + if (!found) + conto.Add (n1, n2); + } + } + } - connectedtonode.SetSize(np); - for (i = 1; i <= np; i++) - connectedtonode.Elem(i) = 0; + connectedtonode.SetSize(np); + for (i = 1; i <= np; i++) + connectedtonode.Elem(i) = 0; - for (i = 1; i <= np; i++) - if (connectedtonode.Elem(i) == 0) - { - connectedtonode.Elem(i) = i; - ConnectToNodeRec (i, i, conto); - } + for (i = 1; i <= np; i++) + if (connectedtonode.Elem(i) == 0) + { + connectedtonode.Elem(i) = i; + ConnectToNodeRec (i, i, conto); + } - } + } - void Mesh :: ConnectToNodeRec (int node, int tonode, - const TABLE & conto) - { - int i, n2; - // (*testout) << "connect " << node << " to " << tonode << endl; - for (i = 1; i <= conto.EntrySize(node); i++) - { - n2 = conto.Get(node, i); - if (!connectedtonode.Get(n2)) - { - connectedtonode.Elem(n2) = tonode; - ConnectToNodeRec (n2, tonode, conto); - } - } - } - */ + void Mesh :: ConnectToNodeRec (int node, int tonode, + const TABLE & conto) + { + int i, n2; + // (*testout) << "connect " << node << " to " << tonode << endl; + for (i = 1; i <= conto.EntrySize(node); i++) + { + n2 = conto.Get(node, i); + if (!connectedtonode.Get(n2)) + { + connectedtonode.Elem(n2) = tonode; + ConnectToNodeRec (n2, tonode, conto); + } + } + } + */ - bool Mesh :: PureTrigMesh (int faceindex) const - { - if (!faceindex) - return !mparam.quad; + bool Mesh :: PureTrigMesh (int faceindex) const + { + if (!faceindex) + return !mparam.quad; - int i; - for (i = 1; i <= GetNSE(); i++) - if (SurfaceElement(i).GetIndex() == faceindex && - SurfaceElement(i).GetNP() != 3) - return 0; - return 1; - } + int i; + for (i = 1; i <= GetNSE(); i++) + if (SurfaceElement(i).GetIndex() == faceindex && + SurfaceElement(i).GetNP() != 3) + return 0; + return 1; + } - bool Mesh :: PureTetMesh () const - { - for (ElementIndex ei = 0; ei < GetNE(); ei++) - if (VolumeElement(ei).GetNP() != 4) - return 0; - return 1; - } + bool Mesh :: PureTetMesh () const + { + for (ElementIndex ei = 0; ei < GetNE(); ei++) + if (VolumeElement(ei).GetNP() != 4) + return 0; + return 1; + } - void Mesh :: UpdateTopology() - { - topology->Update(); - clusters->Update(); - } + void Mesh :: UpdateTopology() + { + topology->Update(); + clusters->Update(); + } - void Mesh :: SetMaterial (int domnr, const char * mat) - { - if (domnr > materials.Size()) + void Mesh :: SetMaterial (int domnr, const char * mat) + { + if (domnr > materials.Size()) { - int olds = materials.Size(); - materials.SetSize (domnr); - for (int i = olds; i < domnr; i++) - materials[i] = 0; + int olds = materials.Size(); + materials.SetSize (domnr); + for (int i = olds; i < domnr; i++) + materials[i] = 0; } - materials.Elem(domnr) = new char[strlen(mat)+1]; - strcpy (materials.Elem(domnr), mat); - } + materials.Elem(domnr) = new char[strlen(mat)+1]; + strcpy (materials.Elem(domnr), mat); + } - const char * Mesh :: GetMaterial (int domnr) const - { - if (domnr <= materials.Size()) - return materials.Get(domnr); - return 0; - } + const char * Mesh :: GetMaterial (int domnr) const + { + if (domnr <= materials.Size()) + return materials.Get(domnr); + return 0; + } - void Mesh ::SetNBCNames ( int nbcn ) - { - if ( bcnames.Size() ) - for ( int i = 0; i < bcnames.Size(); i++) - if ( bcnames[i] ) delete bcnames[i]; - bcnames.SetSize(nbcn); - bcnames = 0; - } + void Mesh ::SetNBCNames ( int nbcn ) + { + if ( bcnames.Size() ) + for ( int i = 0; i < bcnames.Size(); i++) + if ( bcnames[i] ) delete bcnames[i]; + bcnames.SetSize(nbcn); + bcnames = 0; + } - void Mesh ::SetBCName ( int bcnr, const string & abcname ) - { - if ( bcnames[bcnr] ) delete bcnames[bcnr]; - if ( abcname != "default" ) - bcnames[bcnr] = new string ( abcname ); - else - bcnames[bcnr] = 0; - } + void Mesh ::SetBCName ( int bcnr, const string & abcname ) + { + if ( bcnames[bcnr] ) delete bcnames[bcnr]; + if ( abcname != "default" ) + bcnames[bcnr] = new string ( abcname ); + else + bcnames[bcnr] = 0; + } - string Mesh ::GetBCName ( int bcnr ) const - { - if ( !bcnames.Size() ) - return "default"; - if ( bcnames[bcnr] ) - return *bcnames[bcnr]; - else - return "default"; - } + string Mesh ::GetBCName ( int bcnr ) const + { + if ( !bcnames.Size() ) + return "default"; + if ( bcnames[bcnr] ) + return *bcnames[bcnr]; + else + return "default"; + } - void Mesh :: SetUserData(const char * id, Array & data) - { - if(userdata_int.Used(id)) - delete userdata_int.Get(id); + void Mesh :: SetUserData(const char * id, Array & data) + { + if(userdata_int.Used(id)) + delete userdata_int.Get(id); - Array * newdata = new Array(data); + Array * newdata = new Array(data); - userdata_int.Set(id,newdata); - } - bool Mesh :: GetUserData(const char * id, Array & data, int shift) const - { - if(userdata_int.Used(id)) + userdata_int.Set(id,newdata); + } + bool Mesh :: GetUserData(const char * id, Array & data, int shift) const + { + if(userdata_int.Used(id)) { - if(data.Size() < (*userdata_int.Get(id)).Size()+shift) - data.SetSize((*userdata_int.Get(id)).Size()+shift); - for(int i=0; i<(*userdata_int.Get(id)).Size(); i++) - data[i+shift] = (*userdata_int.Get(id))[i]; - return true; + if(data.Size() < (*userdata_int.Get(id)).Size()+shift) + data.SetSize((*userdata_int.Get(id)).Size()+shift); + for(int i=0; i<(*userdata_int.Get(id)).Size(); i++) + data[i+shift] = (*userdata_int.Get(id))[i]; + return true; } - else + else { - data.SetSize(0); - return false; + data.SetSize(0); + return false; } - } - void Mesh :: SetUserData(const char * id, Array & data) - { - if(userdata_double.Used(id)) - delete userdata_double.Get(id); + } + void Mesh :: SetUserData(const char * id, Array & data) + { + if(userdata_double.Used(id)) + delete userdata_double.Get(id); - Array * newdata = new Array(data); + Array * newdata = new Array(data); - userdata_double.Set(id,newdata); - } - bool Mesh :: GetUserData(const char * id, Array & data, int shift) const - { - if(userdata_double.Used(id)) + userdata_double.Set(id,newdata); + } + bool Mesh :: GetUserData(const char * id, Array & data, int shift) const + { + if(userdata_double.Used(id)) { - if(data.Size() < (*userdata_double.Get(id)).Size()+shift) - data.SetSize((*userdata_double.Get(id)).Size()+shift); - for(int i=0; i<(*userdata_double.Get(id)).Size(); i++) - data[i+shift] = (*userdata_double.Get(id))[i]; - return true; + if(data.Size() < (*userdata_double.Get(id)).Size()+shift) + data.SetSize((*userdata_double.Get(id)).Size()+shift); + for(int i=0; i<(*userdata_double.Get(id)).Size(); i++) + data[i+shift] = (*userdata_double.Get(id))[i]; + return true; } - else + else { - data.SetSize(0); - return false; + data.SetSize(0); + return false; } - } + } - void Mesh :: PrintMemInfo (ostream & ost) const - { - ost << "Mesh Mem:" << endl; + void Mesh :: PrintMemInfo (ostream & ost) const + { + ost << "Mesh Mem:" << endl; - ost << GetNP() << " Points, of size " - << sizeof (Point3d) << " + " << sizeof(POINTTYPE) << " = " - << GetNP() * (sizeof (Point3d) + sizeof(POINTTYPE)) << endl; + ost << GetNP() << " Points, of size " + << sizeof (Point3d) << " + " << sizeof(POINTTYPE) << " = " + << GetNP() * (sizeof (Point3d) + sizeof(POINTTYPE)) << endl; - ost << GetNSE() << " Surface elements, of size " - << sizeof (Element2d) << " = " - << GetNSE() * sizeof(Element2d) << endl; + ost << GetNSE() << " Surface elements, of size " + << sizeof (Element2d) << " = " + << GetNSE() * sizeof(Element2d) << endl; - ost << GetNE() << " Volume elements, of size " - << sizeof (Element) << " = " - << GetNE() * sizeof(Element) << endl; + ost << GetNE() << " Volume elements, of size " + << sizeof (Element) << " = " + << GetNE() * sizeof(Element) << endl; - ost << "surfs on node:"; - surfacesonnode.PrintMemInfo (cout); + ost << "surfs on node:"; + surfacesonnode.PrintMemInfo (cout); - ost << "boundaryedges: "; - if (boundaryedges) - boundaryedges->PrintMemInfo (cout); + ost << "boundaryedges: "; + if (boundaryedges) + boundaryedges->PrintMemInfo (cout); - ost << "surfelementht: "; - if (surfelementht) - surfelementht->PrintMemInfo (cout); - } + ost << "surfelementht: "; + if (surfelementht) + surfelementht->PrintMemInfo (cout); + } } diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp index 60bebf35..1d159c86 100644 --- a/libsrc/meshing/meshing2.cpp +++ b/libsrc/meshing/meshing2.cpp @@ -13,7 +13,7 @@ namespace netgen static Array plainpoints; static Array plainzones; static Array loclines; - static int geomtrig; + // static int geomtrig; //static const char * rname; static int cntelem, trials, nfaces; static int oldnl; @@ -400,13 +400,10 @@ namespace netgen debugflag = debugparam.haltsegment && - ( (debugparam.haltsegmentp1 == gpi1) && - (debugparam.haltsegmentp2 == gpi2) || - (debugparam.haltsegmentp1 == gpi2) && - (debugparam.haltsegmentp2 == gpi1)) || + ( ((debugparam.haltsegmentp1 == gpi1) && (debugparam.haltsegmentp2 == gpi2)) || + ((debugparam.haltsegmentp1 == gpi2) && (debugparam.haltsegmentp2 == gpi1))) || debugparam.haltnode && - ( (debugparam.haltsegmentp1 == gpi1) || - (debugparam.haltsegmentp2 == gpi1)); + ( (debugparam.haltsegmentp1 == gpi1) || (debugparam.haltsegmentp2 == gpi1)); if (debugparam.haltface && debugparam.haltfacenr == facenr) @@ -566,7 +563,7 @@ namespace netgen } } - else if (z1 > 0 && z2 > 0 && (z1 != z2) || (z1 < 0) && (z2 < 0) ) + else if ( (z1 > 0 && z2 > 0 && (z1 != z2)) || ((z1 < 0) && (z2 < 0)) ) { loclines.DeleteElement(i); lindex.DeleteElement(i); diff --git a/libsrc/stlgeom/stlgeommesh.cpp b/libsrc/stlgeom/stlgeommesh.cpp index cda9b176..4eec9043 100644 --- a/libsrc/stlgeom/stlgeommesh.cpp +++ b/libsrc/stlgeom/stlgeommesh.cpp @@ -678,7 +678,7 @@ int STLGeometry :: ProjectOnWholeSurface(Point<3> & p3d) const int STLGeometry :: ProjectNearest(Point<3> & p3d) const { - Point<3> p, pf; + Point<3> p, pf = 0.0; //set new chart const STLChart& chart = GetChart(meshchart); diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp index 432a621e..655e2629 100644 --- a/libsrc/stlgeom/stltool.cpp +++ b/libsrc/stlgeom/stltool.cpp @@ -569,12 +569,10 @@ double STLTriangle :: GetNearestPoint(const Array >& ap, if (PointInside(ap, p)) {p3d = p; return dist;} else { - Point<3> pf; + Point<3> pf = 0.0; double nearest = 1E50; //int fi = 0; - int j; - - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { p = p3d; dist = GetDistFromLine(ap.Get(PNum(j)), ap.Get(PNumMod(j+1)), p); diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index d6e9ef75..1e900c09 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -2166,6 +2166,7 @@ namespace netgen } else if (el.GetType() == QUAD) { + /* Array < Point<3> > lp(3); lp[0] = mesh->Point(el[0]); @@ -2179,106 +2180,106 @@ namespace netgen lp[2] = mesh->Point(el[3]); DrawTrigSurfaceVectors(lp,pmin,pmax,sei,vsol); - - /* - Point<3> lp[4]; - Point<2> p2d[4]; - - for (k = 0; k < 4; k++) - lp[k] = mesh->Point (el[k]); - - - Vec<3> n = Cross (lp[1]-lp[0], lp[2]-lp[0]); - Vec<3> na (fabs (n(0)), fabs(n(1)), fabs(n(2))); - if (na(0) > na(1) && na(0) > na(2)) - dir = 1; - else if (na(1) > na(2)) - dir = 2; - else - dir = 3; - - dir1 = (dir % 3) + 1; - dir2 = (dir1 % 3) + 1; - - for (k = 0; k < 4; k++) - { - p2d[k] = Point<2> ((lp[k](dir1-1) - pmin(dir1-1)) / (2*rad), - (lp[k](dir2-1) - pmin(dir2-1)) / (2*rad)); - } - - double minx2d, maxx2d, miny2d, maxy2d; - minx2d = maxx2d = p2d[0](0); - miny2d = maxy2d = p2d[0](1); - for (k = 1; k < 4; k++) - { - minx2d = min2 (minx2d, p2d[k](0)); - maxx2d = max2 (maxx2d, p2d[k](0)); - miny2d = min2 (miny2d, p2d[k](1)); - maxy2d = max2 (maxy2d, p2d[k](1)); - } - - for (s = xoffset/gridsize; s <= 1+xoffset/gridsize; s += 1.0 / gridsize) - if (s >= minx2d && s <= maxx2d) - for (t = yoffset/gridsize; t <= 1+yoffset/gridsize; t += 1.0 / gridsize) - if (t >= miny2d && t <= maxy2d) - { - double lami[3]; - Point3d p3d(2*rad*s+pmin(0), 2*rad*t+pmin(1),0); - - if (mesh->PointContainedIn2DElement (p3d, lami, sei+1)) - { - Point<3> cp = p3d; - double lam1 = lami[0]; - double lam2 = lami[1]; - - //for (k = 0; k < 3; k++) - //cp(k) = lp[0](k) + - //lam1 * (lp[1](k)-lp[0](k)) + - //lam2 * (lp[2](k)-lp[0](k)); - - - Vec<3> v; - double values[6]; - drawelem = GetSurfValues (vsol, sei, lam1, lam2, values); - (*testout) << "sei " << sei << " lam1 " << lam1 << " lam2 " << lam2 << " drawelem " << drawelem << endl; - - if (!vsol->iscomplex) - for (k = 0; k < 3; k++) - v(k) = values[k]; - else - { - if (!imag_part) - for (k = 0; k < 3; k++) - v(k) = values[2*k]; - else - for (k = 0; k < 3; k++) - v(k) = values[2*k+1]; - } - - if (mesh->GetDimension() == 2) - if ( (!vsol->iscomplex && vsol->components != 3) || - (vsol->iscomplex && vsol->components != 6) ) - v(2) = 0; - - double val = v.Length(); - SetOpenGlColor (val, minval, maxval, logscale); - - (*testout) << "v " << v << endl; - - if (val > 1e-10 * maxval) - v *= (rad / val / gridsize * 0.5); - - (*testout) << "v " << v << endl; - - if ( drawelem ) - { - DrawCone (cp, cp+4*v, 0.8*rad / gridsize); - (*testout) << "cp " << cp << " rad " << rad << " gridsize " << gridsize << endl; - } - - } - } */ + + Point<3> lp[4]; + Point<2> p2d[4]; + + for (int k = 0; k < 4; k++) + lp[k] = mesh->Point (el[k]); + + + Vec<3> n = Cross (lp[1]-lp[0], lp[2]-lp[0]); + Vec<3> na (fabs (n(0)), fabs(n(1)), fabs(n(2))); + int dir, dir1, dir2; + if (na(0) > na(1) && na(0) > na(2)) + dir = 1; + else if (na(1) > na(2)) + dir = 2; + else + dir = 3; + + dir1 = (dir % 3) + 1; + dir2 = (dir1 % 3) + 1; + + for (int k = 0; k < 4; k++) + { + p2d[k] = Point<2> ((lp[k](dir1-1) - pmin(dir1-1)) / (2*rad), + (lp[k](dir2-1) - pmin(dir2-1)) / (2*rad)); + } + + double minx2d, maxx2d, miny2d, maxy2d; + minx2d = maxx2d = p2d[0](0); + miny2d = maxy2d = p2d[0](1); + for (int k = 1; k < 4; k++) + { + minx2d = min2 (minx2d, p2d[k](0)); + maxx2d = max2 (maxx2d, p2d[k](0)); + miny2d = min2 (miny2d, p2d[k](1)); + maxy2d = max2 (maxy2d, p2d[k](1)); + } + + for (double s = xoffset/gridsize; s <= 1+xoffset/gridsize; s += 1.0 / gridsize) + if (s >= minx2d && s <= maxx2d) + for (double t = yoffset/gridsize; t <= 1+yoffset/gridsize; t += 1.0 / gridsize) + if (t >= miny2d && t <= maxy2d) + { + double lami[3]; + Point3d p3d(2*rad*s+pmin(0), 2*rad*t+pmin(1),0); + + if (mesh->PointContainedIn2DElement (p3d, lami, sei+1)) + { + Point<3> cp = p3d; + double lam1 = lami[0]; + double lam2 = lami[1]; + + //for (k = 0; k < 3; k++) + //cp(k) = lp[0](k) + + //lam1 * (lp[1](k)-lp[0](k)) + + //lam2 * (lp[2](k)-lp[0](k)); + + + Vec<3> v; + double values[6]; + bool drawelem = GetSurfValues (vsol, sei, lam1, lam2, values); + (*testout) << "sei " << sei << " lam1 " << lam1 << " lam2 " << lam2 << " drawelem " << drawelem << endl; + + if (!vsol->iscomplex) + for (int k = 0; k < 3; k++) + v(k) = values[k]; + else + { + if (!imag_part) + for (int k = 0; k < 3; k++) + v(k) = values[2*k]; + else + for (int k = 0; k < 3; k++) + v(k) = values[2*k+1]; + } + + if (mesh->GetDimension() == 2) + if ( (!vsol->iscomplex && vsol->components != 3) || + (vsol->iscomplex && vsol->components != 6) ) + v(2) = 0; + + double val = v.Length(); + SetOpenGlColor (val, minval, maxval, logscale); + + (*testout) << "v " << v << endl; + + if (val > 1e-10 * maxval) + v *= (rad / val / gridsize * 0.5); + + (*testout) << "v " << v << endl; + + if ( drawelem ) + { + DrawCone (cp, cp+4*v, 0.8*rad / gridsize); + (*testout) << "cp " << cp << " rad " << rad << " gridsize " << gridsize << endl; + } + + } + } } } }