- libsrc/core/version.hpp | 7 ++++++-
- libsrc/meshing/meshclass.cpp | 6 +++---
- libsrc/occ/occgenmesh.cpp | 13 +++++++------
- libsrc/stlgeom/stlgeom.hpp | 2 +-
- libsrc/stlgeom/stltopology.hpp | 2 +-
- nglib/nglib.cpp | 2 +-
- 6 files changed, 19 insertions(+), 13 deletions(-)
+ libsrc/core/version.hpp | 7 ++-
+ libsrc/meshing/meshclass.cpp | 6 +-
+ libsrc/occ/occgenmesh.cpp | 17 +++---
+ libsrc/occ/occmeshsurf.cpp | 128 +++++++++++++----------------------------
+ libsrc/stlgeom/stlgeom.hpp | 2 +-
+ libsrc/stlgeom/stltopology.hpp | 2 +-
+ nglib/nglib.cpp | 2 +-
+ 7 files changed, 62 insertions(+), 102 deletions(-)
diff --git a/libsrc/core/version.hpp b/libsrc/core/version.hpp
index 3048ce5b..81524b58 100644
for (auto & fd : facedecoding)
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
-index e0164177..884d6d1d 100644
+index e0164177..431050e0 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/libsrc/occ/occgenmesh.cpp
@@ -7,6 +7,7 @@
#include <BRepAdaptor_Curve.hxx>
#include <BRepGProp.hxx>
#include <BRepLProp_CLProps.hxx>
+@@ -59,7 +60,7 @@ namespace netgen
+ double ComputeH (double kappa, const MeshingParameters & mparam)
+ {
+ kappa *= mparam.curvaturesafety;
+- /*
++ /**/
+ double hret;
+
+ if (mparam.maxh * kappa < 1)
+@@ -71,7 +72,7 @@ namespace netgen
+ hret = mparam.maxh;
+
+ return hret;
+- */
++ /**/
+ // return min(mparam.maxh, 1/kappa);
+ return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;
+ }
@@ -653,13 +654,13 @@ namespace netgen
for (int k = 1; k <=3; k++)
}
//double maxside = 0;
+diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp
+index 015d4abe..016d759d 100644
+--- a/libsrc/occ/occmeshsurf.cpp
++++ b/libsrc/occ/occmeshsurf.cpp
+@@ -163,50 +163,6 @@ namespace netgen
+ gp_Vec du, dv;
+ occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv);
+
+- // static Timer t("occ-defintangplane calculations");
+- // RegionTimer reg(t);
+-
+- Mat<3,2> D1_;
+- D1_(0,0) = du.X(); D1_(1,0) = du.Y(); D1_(2,0) = du.Z();
+- D1_(0,1) = dv.X(); D1_(1,1) = dv.Y(); D1_(2,1) = dv.Z();
+- auto D1T_ = Trans(D1_);
+- auto D1TD1_ = D1T_*D1_;
+- if (Det (D1TD1_) == 0) throw SingularMatrixException();
+- Mat<2,2> DDTinv_;
+- CalcInverse (D1TD1_, DDTinv_);
+-
+- Mat<3,2> Y_;
+- Vec<3> y1_ = (ap2-ap1).Normalize();
+- Vec<3> y2_ = Cross(n, y1_).Normalize();
+- for (int i = 0; i < 3; i++)
+- {
+- Y_(i,0) = y1_(i);
+- Y_(i,1) = y2_(i);
+- }
+-
+- auto A_ = DDTinv_ * D1T_ * Y_;
+- Mat<2,2> Ainv_;
+- if (Det(A_) == 0) throw SingularMatrixException();
+- CalcInverse (A_, Ainv_);
+-
+- Vec<2> temp_ = Ainv_ * (psp2-psp1);
+- double r_ = temp_.Length();
+- Mat<2,2> R_;
+- R_(0,0) = temp_(0)/r_;
+- R_(1,0) = temp_(1)/r_;
+- R_(0,1) = -R_(1,0);
+- R_(1,1) = R_(0,0);
+-
+- A_ = A_ * R_;
+- Ainv_ = Trans(R_) * Ainv_;
+-
+- Amat = A_;
+- Amatinv = Ainv_;
+-
+- // temp = Amatinv * (psp2-psp1);
+-
+-
+-#ifdef OLD
+ DenseMatrix D1(3,2), D1T(2,3), DDTinv(2,2);
+ D1(0,0) = du.X(); D1(1,0) = du.Y(); D1(2,0) = du.Z();
+ D1(0,1) = dv.X(); D1(1,1) = dv.Y(); D1(2,1) = dv.Z();
+@@ -279,8 +235,7 @@ namespace netgen
+ }
+ // cout << "=?= Ainv = " << endl << Ainv << endl;
+ temp = Amatinv * (psp2-psp1);
+- cout << " =?= Amatinv = " << Amatinv << endl;
+-#endif
++ //cout << " =?= Amatinv = " << Amatinv << endl;
+ };
+
+ }
+@@ -369,59 +324,58 @@ namespace netgen
+
+ double u = gi.u;
+ double v = gi.v;
+-#ifdef OLD
++ if (0) { // exclude Newton's method
+ // was a problem for pheres: got u-v parameters outside range of definition
+-
+- gp_Pnt x = occface->Value (u,v);
+
+- if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return;
++ gp_Pnt x = occface->Value(u, v);
+
+- gp_Vec du, dv;
+- occface->D1(u,v,x,du,dv);
++ if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return;
+
+- int count = 0;
+-
+- gp_Pnt xold;
+- gp_Vec n;
+- double det, lambda, mu;
+-
+- do
++ gp_Vec du, dv;
++ occface->D1(u, v, x, du, dv);
++
++ int count = 0;
++
++ gp_Pnt xold;
++ gp_Vec n;
++ double det, lambda, mu;
++
++ do
+ {
+ count++;
+-
+- n = du^dv;
+-
+- det = Det3 (n.X(), du.X(), dv.X(),
+- n.Y(), du.Y(), dv.Y(),
+- n.Z(), du.Z(), dv.Z());
+-
++
++ n = du ^ dv;
++
++ det = Det3(n.X(), du.X(), dv.X(),
++ n.Y(), du.Y(), dv.Y(),
++ n.Z(), du.Z(), dv.Z());
++
+ if (det < 1e-15)
+ break;
+-
+- lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
+- n.Y(), p.Y()-x.Y(), dv.Y(),
+- n.Z(), p.Z()-x.Z(), dv.Z())/det;
+-
+- mu = Det3 (n.X(), du.X(), p.X()-x.X(),
+- n.Y(), du.Y(), p.Y()-x.Y(),
+- n.Z(), du.Z(), p.Z()-x.Z())/det;
+-
++
++ lambda = Det3(n.X(), p.X() - x.X(), dv.X(),
++ n.Y(), p.Y() - x.Y(), dv.Y(),
++ n.Z(), p.Z() - x.Z(), dv.Z()) / det;
++
++ mu = Det3(n.X(), du.X(), p.X() - x.X(),
++ n.Y(), du.Y(), p.Y() - x.Y(),
++ n.Z(), du.Z(), p.Z() - x.Z()) / det;
++
+ u += lambda;
+ v += mu;
+-
++
+ xold = x;
+- occface->D1(u,v,x,du,dv);
+-
++ occface->D1(u, v, x, du, dv);
++
+ if (xold.SquareDistance(x) < sqr(PROJECTION_TOLERANCE))
+- {
+- ap = Point<3> (x.X(), x.Y(), x.Z());
+- gi.u = u;
+- gi.v = v;
+- return;
+- }
+- }
+- while (count < 20);
+-#endif
++ {
++ ap = Point<3>(x.X(), x.Y(), x.Z());
++ gi.u = u;
++ gi.v = v;
++ return;
++ }
++ } while (count < 20);
++ }
+
+ // Newton did not converge, use OCC projection
+
diff --git a/libsrc/stlgeom/stlgeom.hpp b/libsrc/stlgeom/stlgeom.hpp
index 8cc5739e..46bad1dc 100644
--- a/libsrc/stlgeom/stlgeom.hpp
