use VecInSolid2 in SpecialPoint analysis, VecInSolid2 return also does_intersect

This commit is contained in:
Joachim Schöberl 2020-10-17 13:01:07 +02:00
parent 01c1411d65
commit decb6c6e90
5 changed files with 89 additions and 10 deletions

View File

@ -540,7 +540,7 @@ namespace netgen
double lam3 = 1-lam1-lam2; double lam3 = 1-lam1-lam2;
if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam3 >= -eps_base1) if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam3 >= -eps_base1)
{ // point is close to trianlge, perturbe by alpha*v { // point is close to trianlge, perturbe by alpha*v
double dlamn = face.nn*v; double dlamn = face.nn*v;
if (fabs(dlamn) < 1e-8) // vec also in plane if (fabs(dlamn) < 1e-8) // vec also in plane

View File

@ -285,6 +285,7 @@ namespace netgen
} }
/*
bool Solid::VectorIn2 (const Point<3> & p, const Vec<3> & v1, bool Solid::VectorIn2 (const Point<3> & p, const Vec<3> & v1,
const Vec<3> & v2, double eps) const const Vec<3> & v2, double eps) const
{ {
@ -317,10 +318,51 @@ namespace netgen
} }
return 0; return 0;
} }
*/
bool Solid::VectorIn2 (const Point<3> & p, const Vec<3> & v1,
const Vec<3> & v2, double eps) const
{
switch (op)
{
case TERM: case TERM_REF:
{
auto res = prim->VecInSolid2 (p, v1, v2, eps);
return res != IS_OUTSIDE;
}
case SECTION:
return s1->VectorIn2 (p, v1, v2, eps) && s2->VectorIn2 (p, v1, v2, eps);
case UNION:
return s1->VectorIn2 (p, v1, v2, eps) || s2->VectorIn2 (p, v1, v2, eps);
case SUB:
return !s1->VectorStrictIn2 (p, v1, v2, eps);
case ROOT:
return s1->VectorIn2 (p, v1, v2, eps);
}
// return 0;
}
bool Solid :: VectorStrictIn2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
double eps) const
{
switch (op)
{
case TERM: case TERM_REF:
{
auto res = prim->VecInSolid2 (p, v1, v2, eps);
return (res == IS_INSIDE);
}
case SECTION:
return s1->VectorStrictIn2 (p, v1, v2, eps) && s2->VectorStrictIn2 (p, v1, v2, eps);
case UNION:
return s1->VectorStrictIn2 (p, v1, v2, eps) || s2->VectorStrictIn2 (p, v1, v2, eps);
case SUB:
return !s1->VectorIn2 (p, v1, v2, eps);
case ROOT:
return s1->VectorStrictIn2 (p, v1, v2, eps);
}
}
void Solid :: Print (ostream & str) const void Solid :: Print (ostream & str) const

View File

@ -109,9 +109,12 @@ namespace netgen
bool VectorIn2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, bool VectorIn2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
double eps) const; double eps) const;
/*
bool VectorIn2Rec (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, bool VectorIn2Rec (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
double eps) const; double eps) const;
*/
bool VectorStrictIn2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
double eps) const;
/// compute localization in point p /// compute localization in point p
unique_ptr<Solid> TangentialSolid (const Point<3> & p, NgArray<int> & surfids, double eps) const; unique_ptr<Solid> TangentialSolid (const Point<3> & p, NgArray<int> & surfids, double eps) const;

View File

@ -2076,16 +2076,43 @@ namespace netgen
continue; continue;
Vec<3> s = Cross (normalvecs[m], t); Vec<3> s = Cross (normalvecs[m], t);
Vec<3> t2a = t + 0.01 *s; Vec<3> t2a = t + 0.01 *s;
Vec<3> t2b = t - 0.01 *s; Vec<3> t2b = t - 0.01 *s;
bool isface = bool isfaceold =
(locsol->VectorIn (p, t2a, 1e-6*geomsize) && (locsol->VectorIn (p, t2a, 1e-6*geomsize) &&
!locsol->VectorStrictIn (p, t2a, 1e-6*geomsize)) !locsol->VectorStrictIn (p, t2a, 1e-6*geomsize))
|| ||
(locsol->VectorIn (p, t2b, 1e-6*geomsize) && (locsol->VectorIn (p, t2b, 1e-6*geomsize) &&
!locsol->VectorStrictIn (p, t2b, 1e-6*geomsize)); !locsol->VectorStrictIn (p, t2b, 1e-6*geomsize));
bool isfacenew =
(locsol->VectorIn2 (p, t, s, 1e-6*geomsize) && !locsol->VectorStrictIn2 (p, t, s, 1e-6*geomsize)) ||
(locsol->VectorIn2 (p, t, -s, 1e-6*geomsize) && !locsol->VectorStrictIn2 (p, t, -s, 1e-6*geomsize));
bool isface = isfacenew;
if (isfaceold != isfacenew)
{
*testout << "different, p = " << p << ", t = " << t << ", s = " << s << endl;
*testout << "tlo = " << si << endl;
*testout << "isface, old = " << isface << ", isfacenew = " << isfacenew << endl;
*testout << "t2a = " << t2a << endl;
*testout << "vecin(p,t2a) = " << locsol->VectorIn (p, t2a, 1e-6*geomsize) << endl;
*testout << "vecstrictin(p,t2a) = " << locsol->VectorStrictIn (p, t2a, 1e-6*geomsize) << endl;
*testout << "vectorin2 = " << locsol->VectorIn2 (p, t, s, 1e-6*geomsize) << endl;
*testout << "vectorstrictin2 = " << locsol->VectorStrictIn2 (p, t, s, 1e-6*geomsize) << endl;
*testout << "t2b = " << t2b << endl;
*testout << "vecin(p,t2b) = " << locsol->VectorIn (p, t2b, 1e-6*geomsize) << endl;
*testout << "vecstrictin(p,t2b) = " << locsol->VectorStrictIn (p, t2b, 1e-6*geomsize) << endl;
*testout << "vectorin2- = " << locsol->VectorIn2 (p, t, -s, 1e-6*geomsize) << endl;
*testout << "vectorstrictin2- = " << locsol->VectorStrictIn2 (p, t, -s, 1e-6*geomsize) << endl;
}
/* /*
bool isface = bool isface =
(locsol->VectorIn (p, t2a) && (locsol->VectorIn (p, t2a) &&
@ -2095,10 +2122,8 @@ namespace netgen
!locsol->VectorStrictIn (p, t2b)); !locsol->VectorStrictIn (p, t2b));
*/ */
if (isface) if (isface)
{ cnts++;
cnts++;
}
} }
if (cnts < 2) isedge = 0; if (cnts < 2) isedge = 0;
} }

View File

@ -427,11 +427,20 @@ VecInSolid2 (const Point<3> & p,
if (hv1 >= eps) if (hv1 >= eps)
return IS_OUTSIDE; return IS_OUTSIDE;
double hv2 = v2 * hv;
if (hv2 <= -eps)
return IS_INSIDE;
if (hv2 >= eps)
return IS_OUTSIDE;
return DOES_INTERSECT;
/*
double hv2 = v2 * hv; double hv2 = v2 * hv;
if (hv2 <= 0) if (hv2 <= 0)
return IS_INSIDE; return IS_INSIDE;
else else
return IS_OUTSIDE; return IS_OUTSIDE;
*/
} }