ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-26 05:50:32 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'pointoutofquad' into 'master'

Browse Source 
point out of quad fix

Now the error "point out of domain" does not appear any more, but the values are wrong. [pointoutofquad.py](/uploads/91e142ab369491ac5a4e4a8d9c93d671/pointoutofquad.py)

See merge request !28
This commit is contained in:
Joachim Schöberl 2017-01-06 22:16:53 +01:00
commit f44a32e4de
1 changed files with 87 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

134
libsrc/meshing/meshclass.cpp
View File

@ -4299,26 +4299,26 @@ namespace netgen
//SZ //SZ
if(SurfaceElement(element).GetType()==QUAD) if(SurfaceElement(element).GetType()==QUAD)
{ {
//cout << "check quad element" << endl;
const Element2d & el = SurfaceElement(element); const Element2d & el = SurfaceElement(element);
const Point3d & p1 = Point(el.PNum(1)); const Point3d & p1 = Point(el.PNum(1));
const Point3d & p2 = Point(el.PNum(2)); const Point3d & p2 = Point(el.PNum(2));
const Point3d & p3 = Point(el.PNum(3)); const Point3d & p3 = Point(el.PNum(3));
const Point3d & p4 = Point(el.PNum(4)); const Point3d & p4 = Point(el.PNum(4));
// Coefficients of Bilinear Mapping from Ref-Elem to global Elem // Coefficients of Bilinear Mapping from Ref-Elem to global Elem
// X = a + b x + c y + d x y // X = a + b x + c y + d x y
Vec3d a = p1; Vec3d a = p1;
Vec3d b = p2 - a; Vec3d b = p2 - a;
Vec3d c = p4 - a; Vec3d c = p4 - a;
Vec3d d = p3 - a - b - c; Vec3d d = p3 - a - b - c;
double dxb = d.X()*b.Y()-d.Y()*b.X(); double dxb = d.X()*b.Y()-d.Y()*b.X();
double dxc = d.X()*c.Y()-d.Y()*c.X(); double dxc = d.X()*c.Y()-d.Y()*c.X();
double dxa = d.X()*a.Y()-d.Y()*a.X(); double dxa = d.X()*a.Y()-d.Y()*a.X();
double dxp = d.X()*p.Y()-d.Y()*p.X(); double dxp = d.X()*p.Y()-d.Y()*p.X();
double c0,c1,c2; // ,rt; double c0,c1,c2; // ,rt;
lami[2]=0.; lami[2]=0.;
@ -4328,17 +4328,20 @@ namespace netgen
Vec3d dp24 = p4-p2; Vec3d dp24 = p4-p2;
double d1 = dp13.Length2(); double d1 = dp13.Length2();
double d2 = dp24.Length2(); double d2 = dp24.Length2();
// if(fabs(d.X()) <= eps && fabs(d.Y())<= eps) // if(fabs(d.X()) <= eps && fabs(d.Y())<= eps)
//if (d.Length2() < sqr(eps)) //if (d.Length2() < sqr(eps))
if (d.Length2() < sqr(eps)*d1 && d.Length2() < sqr(eps)*d2) if (d.Length2() < sqr(eps)*d1 && d.Length2() < sqr(eps)*d2)
{ {
//Solve Linear System //Solve Linear System
Vec2d sol;
Vec2d sol;
SolveLinearSystemLS (b, c, p-a, sol); SolveLinearSystemLS (b, c, p-a, sol);
lami[0] = sol.X(); lami[0] = sol.X();
lami[1] = sol.Y(); lami[1] = sol.Y();
if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
/* /*
lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/ 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());
@ -4347,66 +4350,103 @@ namespace netgen
*/ */
} }
else else
//if(fabs(dxb) <= eps) if(fabs(dxb) <= eps*fabs(dxc))
if(fabs(dxb) <= eps*fabs(dxc))
{ {
lami[1] = (dxp-dxa)/dxc; lami[1] = (dxp-dxa)/dxc;
//if(fabs(b.X()-d.X()*lami[1])>=eps) if(fabs(b.X()+d.X()*lami[1])>=fabs(b.Y()+d.Y()*lami[1]))
if(fabs(b.X()+d.X()*lami[1])>=fabs(b.Y()+d.Y()*lami[1]))
lami[0] = (p.X()-a.X() - c.X()*lami[1])/(b.X()+d.X()*lami[1]); lami[0] = (p.X()-a.X() - c.X()*lami[1])/(b.X()+d.X()*lami[1]);
else else
lami[0] = (p.Y()-a.Y() - c.Y()*lami[1])/(b.Y()+d.Y()*lami[1]); lami[0] = (p.Y()-a.Y() - c.Y()*lami[1])/(b.Y()+d.Y()*lami[1]);
if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
} }
else else
//if(fabs(dxc) <= eps) if(fabs(dxc) <= eps*fabs(dxb))
if(fabs(dxc) <= eps*fabs(dxb))
{ {
lami[0] = (dxp-dxa)/dxb; lami[0] = (dxp-dxa)/dxb;
//if(fabs(c.X()-d.X()*lami[0])>=eps) if(fabs(c.X()+d.X()*lami[0])>=fabs(c.Y()+d.Y()*lami[0]))
if(fabs(c.X()+d.X()*lami[0])>=fabs(c.Y()+d.Y()*lami[0]))
lami[1] = (p.X()-a.X() - b.X()*lami[0])/(c.X()+d.X()*lami[0]); lami[1] = (p.X()-a.X() - b.X()*lami[0])/(c.X()+d.X()*lami[0]);
else else
lami[1] = (p.Y()-a.Y() - b.Y()*lami[0])/(c.Y()+d.Y()*lami[0]); lami[1] = (p.Y()-a.Y() - b.Y()*lami[0])/(c.Y()+d.Y()*lami[0]);
if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
} }
else //Solve quadratic equation else //Solve quadratic equation
{ {
//if(fabs(d.X()) >= eps) c2 = -d.X()*dxb;
if(fabs(d.X()) >= fabs(d.Y())) c1 = b.X()*dxc - c.X()*dxb + d.X()*(dxp-dxa);
{ c0 = c.X()*(dxp-dxa) + (a.X()-p.X())*dxc;
c2 = d.X()*dxc; double rt = c1*c1 - 4*c2*c0;
//c1 = d.X()*dxc - c.X()*dxb - d.X()*(dxp-dxa);
c1 = b.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);
c1 = b.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; if (rt < 0.) return false;
lami[1] = (-c1 + sqrt(rt))/2/c2; lami[1] = (-c1 + sqrt(rt))/2/c2;
if(lami[1]<=1.+eps && lami[1]>=0.-eps) if(lami[1]<=1.+eps && lami[1]>=0.-eps)
{ {
lami[0] = (dxp - dxa -dxc*lami[1])/dxb; lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
if(lami[0]<=1.+eps && lami[0]>=0.-eps) if(lami[0]<=1.+eps && lami[0]>=0.-eps)
return true; return true;
} }
lami[1] = (-c1 - sqrt(rt))/2/c2;
lami[1] = (-c1 - sqrt(rt))/2/c2; lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
lami[0] = (dxp - dxa -dxc*lami[1])/dxb;
} if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
c2 = d.Y()*dxb;
c1 = b.Y()*dxc - c.Y()*dxb + d.Y()*(dxp-dxa);
c0 = c.Y()*(dxp -dxa) + (a.Y()-p.Y())*dxc;
rt = c1*c1 - 4*c2*c0;
if (rt < 0.) return false;
lami[1] = (-c1 + sqrt(rt))/2/c2;
if(lami[1]<=1.+eps && lami[1]>=0.-eps)
{
lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
if(lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
}
lami[1] = (-c1 - sqrt(rt))/2/c2;
lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
c2 = -d.X()*dxc;
c1 = -b.X()*dxc + c.X()*dxb + d.X()*(dxp-dxa);
c0 = b.X()*(dxp -dxa) + (a.X()-p.X())*dxb;
rt = c1*c1 - 4*c2*c0;
if (rt < 0.) return false;
lami[1] = (-c1 + sqrt(rt))/2/c2;
if(lami[1]<=1.+eps && lami[1]>=0.-eps)
{
lami[0] = (dxp - dxa -dxc*lami[1])/dxb;
if(lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
}
lami[1] = (-c1 - sqrt(rt))/2/c2;
lami[0] = (dxp - dxa -dxc*lami[1])/dxb;
if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
return true;
}
//cout << "lam0,1 = " << lami[0] << ", " << lami[1] << endl; //cout << "lam0,1 = " << lami[0] << ", " << lami[1] << endl;
if( lami[0] <= 1.+eps && lami[0] >= -eps && lami[1]<=1.+eps && lami[1]>=-eps) /*if( lami[0] <= 1.+eps && lami[0] >= -eps && lami[1]<=1.+eps && lami[1]>=-eps)
{ {
if(consider3D) if(consider3D)
{ {
@ -4419,7 +4459,7 @@ namespace netgen
} }
else else
return true; return true;
} }*/
return false; return false;