netgen/libsrc/occ/occmeshsurf.cpp
2019-07-27 19:05:43 +02:00

738 lines
17 KiB
C++

#ifdef OCCGEOMETRY
#include <mystdlib.h>
#include <occgeom.hpp>
#include <meshing.hpp>
#include <GeomLProp_SLProps.hxx>
#include <ShapeAnalysis_Surface.hxx>
namespace netgen
{
#include "occmeshsurf.hpp"
bool glob_testout(false);
void OCCSurface :: GetNormalVector (const Point<3> & p,
const PointGeomInfo & geominfo,
Vec<3> & n) const
{
gp_Pnt pnt;
gp_Vec du, dv;
/*
double gu = geominfo.u;
double gv = geominfo.v;
if (fabs (gu) < 1e-3) gu = 0;
if (fabs (gv) < 1e-3) gv = 0;
occface->D1(gu,gv,pnt,du,dv);
*/
/*
occface->D1(geominfo.u,geominfo.v,pnt,du,dv);
n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
Vec<3>(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
*/
GeomLProp_SLProps lprop(occface,geominfo.u,geominfo.v,1,1e-5);
double setu=geominfo.u,setv=geominfo.v;
if(lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5)
{
double ustep = 0.01*(umax-umin);
double vstep = 0.01*(vmax-vmin);
n=0;
while(setu < umax && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5))
setu += ustep;
if(setu < umax)
{
lprop.SetParameters(setu,setv);
n(0)+=lprop.Normal().X();
n(1)+=lprop.Normal().Y();
n(2)+=lprop.Normal().Z();
}
setu = geominfo.u;
while(setu > umin && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5))
setu -= ustep;
if(setu > umin)
{
lprop.SetParameters(setu,setv);
n(0)+=lprop.Normal().X();
n(1)+=lprop.Normal().Y();
n(2)+=lprop.Normal().Z();
}
setu = geominfo.u;
while(setv < vmax && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5))
setv += ustep;
if(setv < vmax)
{
lprop.SetParameters(setu,setv);
n(0)+=lprop.Normal().X();
n(1)+=lprop.Normal().Y();
n(2)+=lprop.Normal().Z();
}
setv = geominfo.v;
while(setv > vmin && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5))
setv -= ustep;
if(setv > vmin)
{
lprop.SetParameters(setu,setv);
n(0)+=lprop.Normal().X();
n(1)+=lprop.Normal().Y();
n(2)+=lprop.Normal().Z();
}
setv = geominfo.v;
n.Normalize();
}
else
{
n(0)=lprop.Normal().X();
n(1)=lprop.Normal().Y();
n(2)=lprop.Normal().Z();
}
if(glob_testout)
{
(*testout) << "u " << geominfo.u << " v " << geominfo.v
<< " du " << lprop.D1U().X() << " "<< lprop.D1U().Y() << " "<< lprop.D1U().Z()
<< " dv " << lprop.D1V().X() << " "<< lprop.D1V().Y() << " "<< lprop.D1V().Z() << endl;
}
if (orient == TopAbs_REVERSED) n = -1*n;
// (*testout) << "GetNormalVector" << endl;
}
void OCCSurface :: DefineTangentialPlane (const Point<3> & ap1,
const PointGeomInfo & geominfo1,
const Point<3> & ap2,
const PointGeomInfo & geominfo2)
{
if (projecttype == PLANESPACE)
{
p1 = ap1; p2 = ap2;
//cout << "p1 = " << p1 << endl;
//cout << "p2 = " << p2 << endl;
GetNormalVector (p1, geominfo1, ez);
ex = p2 - p1;
ex -= (ex * ez) * ez;
ex.Normalize();
ey = Cross (ez, ex);
GetNormalVector (p2, geominfo2, n2);
nmid = 0.5*(n2+ez);
ez = nmid;
ez.Normalize();
ex = (p2 - p1).Normalize();
ez -= (ez * ex) * ex;
ez.Normalize();
ey = Cross (ez, ex);
nmid = ez;
//cout << "ex " << ex << " ey " << ey << " ez " << ez << endl;
}
else
{
if ( (geominfo1.u < umin) ||
(geominfo1.u > umax) ||
(geominfo2.u < umin) ||
(geominfo2.u > umax) ||
(geominfo1.v < vmin) ||
(geominfo1.v > vmax) ||
(geominfo2.v < vmin) ||
(geominfo2.v > vmax) ) throw UVBoundsException();
p1 = ap1; p2 = ap2;
psp1 = Point<2>(geominfo1.u, geominfo1.v);
psp2 = Point<2>(geominfo2.u, geominfo2.v);
Vec<3> n;
GetNormalVector (p1, geominfo1, n);
gp_Pnt pnt;
gp_Vec du, dv;
occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv);
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();
/*
(*testout) << "DefineTangentialPlane" << endl
<< "---------------------" << endl;
(*testout) << "D1 = " << endl << D1 << endl;
*/
Transpose (D1, D1T);
DenseMatrix D1TD1(3,3);
D1TD1 = D1T*D1;
if (D1TD1.Det() == 0) throw SingularMatrixException();
CalcInverse (D1TD1, DDTinv);
DenseMatrix Y(3,2);
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);
}
DenseMatrix A(2,2);
A = DDTinv * D1T * Y;
DenseMatrix Ainv(2,2);
if (A.Det() == 0) throw SingularMatrixException();
CalcInverse (A, Ainv);
for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; j++)
{
Amat(i,j) = A(i,j);
Amatinv(i,j) = Ainv(i,j);
}
Vec<2> temp = Amatinv * (psp2-psp1);
double r = temp.Length();
// double alpha = -acos (temp(0)/r);
double alpha = -atan2 (temp(1),temp(0));
DenseMatrix R(2,2);
R(0,0) = cos (alpha);
R(1,0) = -sin (alpha);
R(0,1) = sin (alpha);
R(1,1) = cos (alpha);
A = A*R;
if (A.Det() == 0) throw SingularMatrixException();
CalcInverse (A, Ainv);
for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; j++)
{
Amat(i,j) = A(i,j);
Amatinv(i,j) = Ainv(i,j);
}
temp = Amatinv * (psp2-psp1);
};
}
void OCCSurface :: ToPlane (const Point<3> & p3d,
const PointGeomInfo & geominfo,
Point<2> & pplane,
double h, int & zone) const
{
if (projecttype == PLANESPACE)
{
Vec<3> p1p, n;
GetNormalVector (p3d, geominfo, n);
p1p = p3d - p1;
pplane(0) = (p1p * ex) / h;
pplane(1) = (p1p * ey) / h;
if (n * nmid < 0)
zone = -1;
else
zone = 0;
/*
if(zone == -1)
{
(*testout) << "zone = -1 for " << p3d << " 2D: " << pplane << " n " << n << " nmid " << nmid << endl;
glob_testout = true;
GetNormalVector (p3d, geominfo, n);
glob_testout = false;
}
*/
}
else
{
pplane = Point<2>(geominfo.u, geominfo.v);
// (*testout) << "(u,v) = " << geominfo.u << ", " << geominfo.v << endl;
pplane = Point<2> (1/h * (Amatinv * (pplane-psp1)));
// pplane = Point<2> (h * (Amatinv * (pplane-psp1)));
// pplane = Point<2> (1/h * ((pplane-psp1)));
zone = 0;
};
}
void OCCSurface :: FromPlane (const Point<2> & pplane,
Point<3> & p3d,
PointGeomInfo & gi,
double h)
{
if (projecttype == PLANESPACE)
{
// cout << "2d : " << pplane << endl;
p3d = p1 + (h * pplane(0)) * ex + (h * pplane(1)) * ey;
// cout << "3d : " << p3d << endl;
Project (p3d, gi);
// cout << "proj : " << p3d << endl;
}
else
{
// Point<2> pspnew = Point<2>(1/h * (Amat * Vec<2>(pplane)) + Vec<2>(psp1));
Point<2> pspnew = Point<2>(h * (Amat * Vec<2>(pplane)) + Vec<2>(psp1));
// Point<2> pspnew = Point<2>(h * (Vec<2>(pplane)) + Vec<2>(psp1));
gi.u = pspnew(0);
gi.v = pspnew(1);
gi.trignum = 1;
gp_Pnt val = occface->Value (gi.u, gi.v);
p3d = Point<3> (val.X(), val.Y(), val.Z());
};
}
void OCCSurface :: Project (Point<3> & p, PointGeomInfo & gi)
{
// static int cnt = 0;
// if (cnt++ % 1000 == 0) cout << "********************************************** OCCSurfce :: Project, cnt = " << cnt << endl;
gp_Pnt pnt(p(0), p(1), p(2));
//(*testout) << "pnt = " << pnt.X() << ", " << pnt.Y() << ", " << pnt.Z() << endl;
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
if (!proj.NbPoints())
{
cout << "Project Point on Surface FAIL" << endl;
throw UVBoundsException();
}
*/
/*
cout << "NP = " << proj.NbPoints() << endl;
for (int i = 1; i <= proj.NbPoints(); i++)
{
gp_Pnt pnt2 = proj.Point(i);
Point<3> p2 = Point<3> (pnt2.X(), pnt2.Y(), pnt2.Z());
cout << i << ". p = " << p2 << ", dist = " << (p2-p).Length() << endl;
}
*/
/*
pnt = proj.NearestPoint();
proj.LowerDistanceParameters (gi.u, gi.v);
*/
double u,v;
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( topods_face ) );
suval.Coord( u, v);
pnt = occface->Value( u, v );
//(*testout) << "pnt(proj) = " << pnt.X() << ", " << pnt.Y() << ", " << pnt.Z() << endl;
gi.u = u;
gi.v = v;
gi.trignum = 1;
p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
}
Meshing2OCCSurfaces :: Meshing2OCCSurfaces (const TopoDS_Shape & asurf,
const Box<3> & abb, int aprojecttype,
const MeshingParameters & mparam)
: Meshing2(mparam, Box<3>(abb.PMin(), abb.PMax())), surface(TopoDS::Face(asurf), aprojecttype)
{
;
}
void Meshing2OCCSurfaces :: DefineTransformation (const Point3d & p1, const Point3d & p2,
const PointGeomInfo * geominfo1,
const PointGeomInfo * geominfo2)
{
((OCCSurface&)surface).DefineTangentialPlane (p1, *geominfo1, p2, *geominfo2);
}
void Meshing2OCCSurfaces :: TransformToPlain (const Point3d & locpoint,
const MultiPointGeomInfo & geominfo,
Point2d & planepoint,
double h, int & zone)
{
Point<2> hp;
surface.ToPlane (locpoint, geominfo.GetPGI(1), hp, h, zone);
planepoint.X() = hp(0);
planepoint.Y() = hp(1);
}
int Meshing2OCCSurfaces :: TransformFromPlain (Point2d & planepoint,
Point3d & locpoint,
PointGeomInfo & gi,
double h)
{
Point<3> hp;
Point<2> hp2 (planepoint.X(), planepoint.Y());
surface.FromPlane (hp2, hp, gi, h);
locpoint = hp;
return 0;
}
double Meshing2OCCSurfaces :: CalcLocalH (const Point3d & p, double gh) const
{
return gh;
}
MeshOptimize2dOCCSurfaces :: MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry)
: MeshOptimize2d(), geometry(ageometry)
{
;
}
void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
geometry.Project (surfind, p);
}
int MeshOptimize2dOCCSurfaces :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
{
double u = gi.u;
double v = gi.v;
Point<3> hp = p;
if (geometry.FastProject (surfind, hp, u, v))
{
p = hp;
return 1;
}
ProjectPoint (surfind, p);
return CalcPointGeomInfo (surfind, gi, p);
}
void MeshOptimize2dOCCSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2,
Point<3> & p) const
{
TopExp_Explorer exp0, exp1;
bool done = false;
Handle(Geom_Curve) c;
for (exp0.Init(geometry.fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next())
for (exp1.Init(geometry.fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next())
{
if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current())))
{
done = true;
double s0, s1;
c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1);
}
}
gp_Pnt pnt(p(0), p(1), p(2));
GeomAPI_ProjectPointOnCurve proj(pnt, c);
pnt = proj.NearestPoint();
p(0) = pnt.X();
p(1) = pnt.Y();
p(2) = pnt.Z();
}
void MeshOptimize2dOCCSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & geominfo, Vec<3> & n) const
{
gp_Pnt pnt;
gp_Vec du, dv;
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
occface->D1(geominfo.u,geominfo.v,pnt,du,dv);
n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
Vec<3>(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n;
// GetNormalVector (surfind, p, n);
}
void MeshOptimize2dOCCSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const
{
// static int cnt = 0;
// if (cnt++ % 1000 == 0) cout << "GetNV cnt = " << cnt << endl;
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) );
suval.Coord( u, v);
pnt = occface->Value( u, v );
gp_Vec du, dv;
occface->D1(u,v,pnt,du,dv);
/*
if (!occface->IsCNu (1) || !occface->IsCNv (1))
(*testout) << "SurfOpt: Differentiation FAIL" << endl;
*/
n = Cross (Vec3d(du.X(), du.Y(), du.Z()),
Vec3d(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n;
}
int MeshOptimize2dOCCSurfaces ::
CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const
{
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return 0;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) );
suval.Coord( u, v);
//pnt = occface->Value( u, v );
gi.u = u;
gi.v = v;
return 1;
}
OCCRefinementSurfaces :: OCCRefinementSurfaces (const OCCGeometry & ageometry)
: Refinement(), geometry(ageometry)
{
;
}
OCCRefinementSurfaces :: ~OCCRefinementSurfaces ()
{
;
}
/*
inline double Det3 (double a00, double a01, double a02,
double a10, double a11, double a12,
double a20, double a21, double a22)
{
return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00;
}
bool ProjectToSurface (gp_Pnt & p, Handle(Geom_Surface) surface, double& u, double& v)
{
gp_Pnt x = surface->Value (u,v);
if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
gp_Vec du, dv;
surface->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());
if (det < 1e-15) return false;
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;
surface->D1(u,v,x,du,dv);
} while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) || count > 50);
if (count > 50) return false;
p = x;
return true;
}
*/
void OCCRefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi > 0)
{
double u = gi1.u+secpoint*(gi2.u-gi1.u);
double v = gi1.v+secpoint*(gi2.v-gi1.v);
auto savept = hnewp;
if (!geometry.FastProject (surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2))
{
// cout << "Fast projection to surface fails! Using OCC projection" << endl;
hnewp = savept;
geometry.Project (surfi, hnewp);
}
newgi.trignum = 1;
newgi.u = u;
newgi.v = v;
}
newp = hnewp;
}
void OCCRefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
double s0, s1;
Point<3> hnewp = p1+secpoint*(p2-p1);
gp_Pnt pnt(hnewp(0), hnewp(1), hnewp(2));
GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(geometry.emap(ap1.edgenr)), s0, s1));
pnt = proj.NearestPoint();
hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
newp = hnewp;
newgi = ap1;
};
void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi) const
{
if (surfi > 0)
geometry.Project (surfi, p);
};
void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const
{
if (surfi > 0)
if (!geometry.FastProject (surfi, p, gi.u, gi.v))
{
cout << "Fast projection to surface fails! Using OCC projection" << endl;
geometry.Project (surfi, p);
}
};
}
#endif