/usr/share/netgen/libsrc/meshing/findip2.hpp is in netgen-headers 4.9.13.dfsg-8build2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 | // find inner point
template <typename POINTArray, typename FACEArray>
inline int FindInnerPoint2 (POINTArray & points,
FACEArray & faces,
Point3d & p)
{
static int timer = NgProfiler::CreateTimer ("FindInnerPoint2");
NgProfiler::RegionTimer reg (timer);
Array<Vec3d> a;
Array<double> c;
Mat<3> m, inv;
Vec<3> rs, x, pmin;
int nf = faces.Size();
a.SetSize (nf);
c.SetSize (nf);
for (int i = 0; i < nf; i++)
{
Point3d p1 = points.Get(faces[i][0]);
a[i] = Cross (points.Get(faces[i][1]) - p1,
points.Get(faces[i][2]) - p1);
a[i] /= a[i].Length();
c[i] = - (a[i].X() * p1.X() + a[i].Y() * p1.Y() + a[i].Z() * p1.Z());
}
x = 0;
double hmax = 0;
for (int i = 0; i < nf; i++)
{
const Element2d & el = faces[i];
for (int j = 1; j <= 3; j++)
{
double hi = Dist (points.Get(el.PNumMod(j)),
points.Get(el.PNumMod(j+1)));
if (hi > hmax) hmax = hi;
}
}
double fmin = 0;
for (int i1 = 1; i1 <= nf; i1++)
for (int i2 = i1+1; i2 <= nf; i2++)
for (int i3 = i2+1; i3 <= nf; i3++)
for (int i4 = i3+1; i4 <= nf; i4++)
{
m(0, 0) = a.Get(i1).X() - a.Get(i2).X();
m(0, 1) = a.Get(i1).Y() - a.Get(i2).Y();
m(0, 2) = a.Get(i1).Z() - a.Get(i2).Z();
rs(0) = c.Get(i2) - c.Get(i1);
m(1, 0) = a.Get(i1).X() - a.Get(i3).X();
m(1, 1) = a.Get(i1).Y() - a.Get(i3).Y();
m(1, 2) = a.Get(i1).Z() - a.Get(i3).Z();
rs(1) = c.Get(i3) - c.Get(i1);
m(2, 0) = a.Get(i1).X() - a.Get(i4).X();
m(2, 1) = a.Get(i1).Y() - a.Get(i4).Y();
m(2, 2) = a.Get(i1).Z() - a.Get(i4).Z();
rs(2) = c.Get(i4) - c.Get(i1);
if (fabs (Det (m)) > 1e-10)
{
CalcInverse (m, inv);
x = inv * rs;
double f = -1e10;
for (int i = 0; i < nf; i++)
{
double hd =
x(0) * a[i].X() + x(1) * a[i].Y() + x(2) * a[i].Z() + c[i];
if (hd > f) f = hd;
if (hd > fmin) break;
}
if (f < fmin)
{
fmin = f;
pmin = x;
}
}
}
p = Point3d (pmin(0), pmin(1), pmin(2));
(*testout) << "fmin = " << fmin << endl;
return (fmin < -1e-3 * hmax);
}
|