ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-26 04:40:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

vector/matrix access, small optimization in CalcTetBadness

Browse Source
This commit is contained in:
Joachim Schoeberl 2009-06-06 19:33:32 +00:00
parent 35bc5819db
commit 611a53801d
34 changed files with 4191 additions and 4314 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
configure.ac
View File

@ -1,4 +1,4 @@
AC_INIT([netgen],[4.9.9],[],[]) AC_INIT([netgen],[4.9.10-dev],[],[])
AM_INIT_AUTOMAKE([-Wall -Werror foreign]) AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_MACRO_DIR([m4])

22
libsrc/geom2d/spline.hpp
View File

@ -587,8 +587,6 @@ Vec<D> SplineSeg3<D> :: GetTangent (const double t) const
template<int D> template<int D>
void SplineSeg3<D> :: GetCoeff (Vector & u) const void SplineSeg3<D> :: GetCoeff (Vector & u) const
{ {
double t;
int i;
Point<D> p; Point<D> p;
DenseMatrix a(6, 6); DenseMatrix a(6, 6);
DenseMatrix ata(6, 6); DenseMatrix ata(6, 6);
@ -598,23 +596,23 @@ void SplineSeg3<D> :: GetCoeff (Vector & u) const
// ata.SetSymmetric(1); // ata.SetSymmetric(1);
t = 0; double t = 0;
for (i = 1; i <= 5; i++, t += 0.25) for (int i = 0; i < 5; i++, t += 0.25)
{ {
p = GetPoint (t); p = GetPoint (t);
a.Elem(i, 1) = p(0) * p(0); a(i, 0) = p(0) * p(0);
a.Elem(i, 2) = p(1) * p(1); a(i, 1) = p(1) * p(1);
a.Elem(i, 3) = p(0) * p(1); a(i, 2) = p(0) * p(1);
a.Elem(i, 4) = p(0); a(i, 3) = p(0);
a.Elem(i, 5) = p(1); a(i, 4) = p(1);
a.Elem(i, 6) = 1; a(i, 5) = 1;
} }
a.Elem(6, 1) = 1; a(5, 0) = 1;
CalcAtA (a, ata); CalcAtA (a, ata);
u = 0; u = 0;
u.Elem(6) = 1; u(5) = 1;
a.MultTrans (u, f); a.MultTrans (u, f);
ata.Solve (f, u); ata.Solve (f, u);
} }

40
libsrc/gprim/geomtest3d.cpp
View File

@ -23,12 +23,12 @@ IntersectTriangleLine (const Point<3> ** tri, const Point<3> ** line)
<< "tri = " << *tri[0] << ", " << *tri[1] << ", " << *tri[2] << endl << "tri = " << *tri[0] << ", " << *tri[1] << ", " << *tri[2] << endl
<< "line = " << *line[0] << ", " << *line[1] << endl; << "line = " << *line[0] << ", " << *line[1] << endl;
*/ */
for (i = 1; i <= 3; i++) for (i = 0; i < 3; i++)
{ {
a.Elem(i, 1) = -vl.X(i); a(i, 0) = -vl.X(i+1);
a.Elem(i, 2) = vt1.X(i); a(i, 1) = vt1.X(i+1);
a.Elem(i, 3) = vt2.X(i); a(i, 2) = vt2.X(i+1);
rs.Elem(i) = vrs.X(i); rs(i) = vrs.X(i+1);
} }
double det = a.Det(); double det = a.Det();
@ -70,12 +70,12 @@ IntersectTriangleLine (const Point<3> ** tri, const Point<3> ** line)
double eps = 1e-6; double eps = 1e-6;
if ( if (
(lami.Get(1) >= -eps && lami.Get(1) <= 1+eps && (lami(0) >= -eps && lami(0) <= 1+eps &&
lami.Get(2) >= -eps && lami.Get(3) >= -eps && lami(1) >= -eps && lami(2) >= -eps &&
lami.Get(2) + lami.Get(3) <= 1+eps) && ! lami(1) + lami(2) <= 1+eps) && !
(lami.Get(1) >= eps && lami.Get(1) <= 1-eps && (lami(0) >= eps && lami(0) <= 1-eps &&
lami.Get(2) >= eps && lami.Get(3) >= eps && lami(1) >= eps && lami(2) >= eps &&
lami.Get(2) + lami.Get(3) <= 1-eps) ) lami(1) + lami(2) <= 1-eps) )
{ {
@ -97,8 +97,8 @@ IntersectTriangleLine (const Point<3> ** tri, const Point<3> ** line)
} }
if (lami.Get(1) >= 0 && lami.Get(1) <= 1 && if (lami(0) >= 0 && lami(0) <= 1 &&
lami.Get(2) >= 0 && lami.Get(3) >= 0 && lami.Get(2) + lami.Get(3) <= 1) lami(1) >= 0 && lami(2) >= 0 && lami(1) + lami(2) <= 1)
{ {
return 1; return 1;
@ -860,12 +860,12 @@ int CalcTriangleCenter (const Point3d ** pts, Point3d & c)
Vec3d v1(*pts[0], *pts[1]); Vec3d v1(*pts[0], *pts[1]);
Vec3d v2(*pts[0], *pts[2]); Vec3d v2(*pts[0], *pts[2]);
rs.Elem(1) = v1 * v1; rs(0) = v1 * v1;
rs.Elem(2) = v2 * v2; rs(1) = v2 * v2;
a.Elem(1,1) = 2 * rs.Get(1); a(0,0) = 2 * rs(0);
a.Elem(1,2) = a.Elem(2,1) = 2 * (v1 * v2); a(0,1) = a(1,0) = 2 * (v1 * v2);
a.Elem(2,2) = 2 * rs.Get(2); a(1,1) = 2 * rs(1);
if (fabs (a.Det()) <= 1e-12 * h * h) if (fabs (a.Det()) <= 1e-12 * h * h)
{ {
@ -877,8 +877,8 @@ int CalcTriangleCenter (const Point3d ** pts, Point3d & c)
inva.Mult (rs, sol); inva.Mult (rs, sol);
c = *pts[0]; c = *pts[0];
v1 *= sol.Get(1); v1 *= sol(0);
v2 *= sol.Get(2); v2 *= sol(1);
c += v1; c += v1;
c += v2; c += v2;

16
libsrc/gprim/transform3d.cpp
View File

@ -90,21 +90,21 @@ void Transformation3d :: CalcInverse (Transformation3d & inv) const
static Vector b(3), sol(3); static Vector b(3), sol(3);
int i, j; int i, j;
for (i = 1; i <= 3; i++) for (int i = 0; i < 3; i++)
{ {
b.Elem(i) = offset[i-1]; b(i) = offset[i];
for (j = 1; j <= 3; j++) for (j = 0; j < 3; j++)
a.Elem(i, j) = lin[i-1][j-1]; a(i, j) = lin[i][j];
} }
::netgen::CalcInverse (a, inva); ::netgen::CalcInverse (a, inva);
inva.Mult (b, sol); inva.Mult (b, sol);
for (i = 1; i <= 3; i++) for (int i = 0; i < 3; i++)
{ {
inv.offset[i-1] = -sol.Get(i); inv.offset[i] = -sol(i);
for (j = 1; j <= 3; j++) for (j = 0; j < 3; j++)
inv.lin[i-1][j-1] = inva.Elem(i, j); inv.lin[i][j] = inva(i, j);
} }
} }

4
libsrc/linalg/Makefile.am
View File

@ -2,6 +2,8 @@ noinst_HEADERS = densemat.hpp linalg.hpp polynomial.hpp vector.hpp opti.hpp
AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include
METASOURCES = AUTO METASOURCES = AUTO
noinst_LTLIBRARIES = libla.la noinst_LTLIBRARIES = libla.la
libla_la_SOURCES = densemat.cpp polynomial.cpp vector.cpp bfgs.cpp linopt.cpp linsearch.cpp libla_la_SOURCES = densemat.cpp polynomial.cpp bfgs.cpp linopt.cpp linsearch.cpp
# vector.cpp
libla_la_LDFLAGS = -rdynamic libla_la_LDFLAGS = -rdynamic

57
libsrc/linalg/bfgs.cpp
View File

@ -1,9 +1,9 @@
/***************************************************************************/ /***************************************************************************/
/* */ /* */
/* Vorlesung Optimierung I, Gfrerer, WS94/95 */ /* Vorlesung Optimierung I, Gfrerer, WS94/95 */
/* BFGS-Verfahren zur Lösung freier nichtlinearer Optimierungsprobleme */ /* BFGS-Verfahren zur Lösung freier nichtlinearer Optimierungsprobleme */
/* */ /* */
/* Programmautor: Joachim Schöberl */ /* Programmautor: Joachim Schöberl */
/* Matrikelnummer: 9155284 */ /* Matrikelnummer: 9155284 */
/* */ /* */
/***************************************************************************/ /***************************************************************************/
@ -25,37 +25,36 @@ void Cholesky (const DenseMatrix & a,
double x; double x;
int i, j, k;
int n = a.Height(); int n = a.Height();
// (*testout) << "a = " << a << endl; // (*testout) << "a = " << a << endl;
l = a; l = a;
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
{ {
for (j = i; j <= n; j++) for (int j = i; j <= n; j++)
{ {
x = l.Get(i, j); x = l.Get(i, j);
for (k = 1; k < i; k++) for (int k = 1; k < i; k++)
x -= l.Get(i, k) * l.Get(j, k) * d.Get(k); x -= l.Get(i, k) * l.Get(j, k) * d(k-1);
if (i == j) if (i == j)
{ {
d.Elem(i) = x; d(i-1) = x;
} }
else else
{ {
l.Elem(j, i) = x / d.Get(k); l.Elem(j, i) = x / d(i-1);
} }
} }
} }
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
{ {
l.Elem(i, i) = 1; l.Elem(i, i) = 1;
for (j = i+1; j <= n; j++) for (int j = i+1; j <= n; j++)
l.Elem(i, j) = 0; l.Elem(i, j) = 0;
} }
@ -162,9 +161,7 @@ int LDLtUpdate (DenseMatrix & l, Vector & d, double a, const Vector & u)
// Rueckgabewert: 0 .. D bleibt positiv definit // Rueckgabewert: 0 .. D bleibt positiv definit
// 1 .. sonst // 1 .. sonst
int i, j, n; int n = l.Height();
n = l.Height();
Vector v(n); Vector v(n);
double t, told, xi; double t, told, xi;
@ -172,9 +169,9 @@ int LDLtUpdate (DenseMatrix & l, Vector & d, double a, const Vector & u)
told = 1; told = 1;
v = u; v = u;
for (j = 1; j <= n; j++) for (int j = 1; j <= n; j++)
{ {
t = told + a * sqr (v.Elem(j)) / d.Get(j); t = told + a * sqr (v(j-1)) / d(j-1);
if (t <= 0) if (t <= 0)
{ {
@ -182,14 +179,14 @@ int LDLtUpdate (DenseMatrix & l, Vector & d, double a, const Vector & u)
return 1; return 1;
} }
xi = a * v.Elem(j) / (d.Get(j) * t); xi = a * v(j-1) / (d(j-1) * t);
d.Elem(j) *= t / told; d(j-1) *= t / told;
for (i = j + 1; i <= n; i++) for (int i = j + 1; i <= n; i++)
{ {
v.Elem(i) -= v.Elem(j) * l.Elem(i, j); v(i-1) -= v(j-1) * l.Elem(i, j);
l.Elem(i, j) += xi * v.Elem(i); l.Elem(i, j) += xi * v(i-1);
} }
told = t; told = t;
@ -209,7 +206,7 @@ double BFGS (
{ {
int i, j, n = x.Size(); int n = x.Size();
long it; long it;
char a1crit, a3acrit; char a1crit, a3acrit;
@ -233,14 +230,14 @@ double BFGS (
int ifail; // o: 0 .. Erfolg int ifail; // o: 0 .. Erfolg
// -1 .. Unterschreitung von fmin // -1 .. Unterschreitung von fmin
// 1 .. kein Erfolg bei Liniensuche // 1 .. kein Erfolg bei Liniensuche
// 2 .. Überschreitung von itmax // 2 .. Überschreitung von itmax
typx = par.typx; typx = par.typx;
typf = par.typf; typf = par.typf;
l = 0; l = 0;
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
l.Elem(i, i) = 1; l.Elem(i, i) = 1;
f = fun.FuncGrad (x, g); f = fun.FuncGrad (x, g);
@ -255,10 +252,10 @@ double BFGS (
if (it % (5 * n) == 0) if (it % (5 * n) == 0)
{ {
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
d.Elem(i) = typf/ sqr (typx.Get(i)); // 1; d(i-1) = typf/ sqr (typx(i-1)); // 1;
for (i = 2; i <= n; i++) for (int i = 2; i <= n; i++)
for (j = 1; j < i; j++) for (int j = 1; j < i; j++)
l.Elem(i, j) = 0; l.Elem(i, j) = 0;
/* /*
@ -357,8 +354,8 @@ double BFGS (
hd = eps * max2 (typf, fabs (f)); hd = eps * max2 (typf, fabs (f));
a1crit = 1; a1crit = 1;
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
if ( fabs (g.Elem(i)) * max2 (typx.Elem(i), fabs (x.Elem(i))) > hd) if ( fabs (g(i-1)) * max2 (typx(i-1), fabs (x(i-1))) > hd)
a1crit = 0; a1crit = 0;

262
libsrc/linalg/densemat.cpp
View File

@ -76,7 +76,7 @@ namespace netgen
/* /*
DenseMatrix & DenseMatrix :: operator= (const BaseMatrix & m2) DenseMatrix & DenseMatrix :: operator= (const BaseMatrix & m2)
{ {
int i, j; int i, j;
@ -132,7 +132,7 @@ DenseMatrix & DenseMatrix :: operator= (const BaseMatrix & m2)
} }
DenseMatrix & DenseMatrix :: operator-= (const DenseMatrix & m2) DenseMatrix & DenseMatrix :: operator-= (const DenseMatrix & m2)
{ {
int i; int i;
double * p, * q; double * p, * q;
@ -164,15 +164,15 @@ DenseMatrix & DenseMatrix :: operator-= (const DenseMatrix & m2)
/* /*
double & DenseMatrix :: operator() (int i, int j) double & DenseMatrix :: operator() (int i, int j)
{ {
if (i >= 1 && j >= 1 && i <= height && j <= width) if (i >= 1 && j >= 1 && i <= height && j <= width)
return Elem(i,j); return Elem(i,j);
else (*myerr) << "DenseMatrix: index (" << i << "," << j << ") out of range (1.." else (*myerr) << "DenseMatrix: index (" << i << "," << j << ") out of range (1.."
<< height << ",1.." << width << ")\n"; << height << ",1.." << width << ")\n";
static double dummy = 0; static double dummy = 0;
return dummy; return dummy;
} }
double DenseMatrix :: operator() (int i, int j) const double DenseMatrix :: operator() (int i, int j) const
{ {
@ -186,7 +186,7 @@ double & DenseMatrix :: operator() (int i, int j)
} }
*/ */
DenseMatrix & DenseMatrix :: operator= (double v) DenseMatrix & DenseMatrix :: operator= (double v)
{ {
int i; int i;
double * p = data; double * p = data;
@ -200,7 +200,7 @@ DenseMatrix & DenseMatrix :: operator= (double v)
DenseMatrix & DenseMatrix :: operator*= (double v) DenseMatrix & DenseMatrix :: operator*= (double v)
{ {
int i; int i;
double * p = data; double * p = data;
@ -213,7 +213,7 @@ DenseMatrix & DenseMatrix :: operator*= (double v)
} }
double DenseMatrix :: Det () const double DenseMatrix :: Det () const
{ {
if (width != height) if (width != height)
{ {
@ -223,15 +223,14 @@ double DenseMatrix :: Det () const
switch (width) switch (width)
{ {
case 1: return Get(1, 1); case 1: return data[0];
case 2: return Get(1) * Get(4) - Get(2) * Get(3); case 2: return data[0] * data[3] - data[1] * data[2];
case 3: return data[0] * data[4] * data[8]
case 3: return Get(1) * Get(5) * Get(9) + data[1] * data[5] * data[6]
+ Get(2) * Get(6) * Get(7) + data[2] * data[3] * data[7]
+ Get(3) * Get(4) * Get(8) - data[0] * data[5] * data[7]
- Get(1) * Get(6) * Get(8) - data[1] * data[3] * data[8]
- Get(2) * Get(4) * Get(9) - data[2] * data[4] * data[6];
- Get(3) * Get(5) * Get(7);
default: default:
{ {
(*myerr) << "Matrix :: Det: general size not implemented (size=" << width << ")" << endl; (*myerr) << "Matrix :: Det: general size not implemented (size=" << width << ")" << endl;
@ -241,11 +240,9 @@ double DenseMatrix :: Det () const
} }
void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2) void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2)
{ {
// int i, j, k, n;
double det; double det;
// DenseMatrix m1 = hm1;
if (m1.width != m1.height) if (m1.width != m1.height)
{ {
@ -268,35 +265,35 @@ void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2)
return; return;
} }
det = 1e0 / det; det = 1.0 / det;
switch (m1.width) switch (m1.width)
{ {
case 1: case 1:
{ {
m2.Set(1, 1, det); m2(0,0) = det;
return; return;
} }
case 2: case 2:
{ {
m2.Set(1, 1, det * m1.Get(4)); m2(0,0) = det * m1(3);
m2.Set(2, 2, det * m1.Get(1)); m2(1,1) = det * m1(0);
m2.Set(1, 2, - det * m1.Get(2)); m2(0,1) = -det * m1(1);
m2.Set(2, 1, - det * m1.Get(3)); m2(1,0) = - det * m1(2);
return; return;
} }
case 3: case 3:
{ {
m2.Set(1, 1, det * (m1.Get(5) * m1.Get(9) - m1.Get(6) * m1.Get(8))); m2(0, 0) = det * (m1(4) * m1(8) - m1(5) * m1(7));
m2.Set(2, 1, -det * (m1.Get(4) * m1.Get(9) - m1.Get(6) * m1.Get(7))); m2(1, 0) = -det * (m1(3) * m1(8) - m1(5) * m1(6));
m2.Set(3, 1, det * (m1.Get(4) * m1.Get(8) - m1.Get(5) * m1.Get(7))); m2(2, 0) = det * (m1(3) * m1(7) - m1(4) * m1(6));
m2.Set(1, 2, -det * (m1.Get(2) * m1.Get(9) - m1.Get(3) * m1.Get(8))); m2(0, 1) = -det * (m1(1) * m1(8) - m1(2) * m1(7));
m2.Set(2, 2, det * (m1.Get(1) * m1.Get(9) - m1.Get(3) * m1.Get(7))); m2(1, 1) = det * (m1(0) * m1(8) - m1(2) * m1(6));
m2.Set(3, 2, -det * (m1.Get(1) * m1.Get(8) - m1.Get(2) * m1.Get(7))); m2(2, 1) = -det * (m1(0) * m1(7) - m1(1) * m1(6));
m2.Set(1, 3, det * (m1.Get(2) * m1.Get(6) - m1.Get(3) * m1.Get(5))); m2(0, 2) = det * (m1(1) * m1(5) - m1(2) * m1(4));
m2.Set(2, 3, -det * (m1.Get(1) * m1.Get(6) - m1.Get(3) * m1.Get(4))); m2(1, 2) = -det * (m1(0) * m1(5) - m1(2) * m1(3));
m2.Set(3, 3, det * (m1.Get(1) * m1.Get(5) - m1.Get(2) * m1.Get(4))); m2(2, 2) = det * (m1(0) * m1(4) - m1(1) * m1(3));
return; return;
} }
} }
@ -367,14 +364,14 @@ void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2)
// check: A = L L^t // check: A = L L^t
// for (i = 1; i <= n; i++) // for (i = 1; i <= n; i++)
// for (j = 1; j <= n; j++) // for (j = 1; j <= n; j++)
// { // {
// x = 0; // x = 0;
// for (k = 1; k <= i && k <= j; k++) // for (k = 1; k <= i && k <= j; k++)
// x += m2.Get(i, k) * m2.Get(j, k); // x += m2.Get(i, k) * m2.Get(j, k);
// (*testout) << "err " << i << "," << j << " = " << (m1.Get(i, j) - x) << endl; // (*testout) << "err " << i << "," << j << " = " << (m1.Get(i, j) - x) << endl;
// } // }
@ -405,15 +402,15 @@ void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2)
} }
} }
// (*testout) << "check L^-1" << endl; // (*testout) << "check L^-1" << endl;
// for (i = 1; i <= n; i++) // for (i = 1; i <= n; i++)
// for (j = 1; j <= n; j++) // for (j = 1; j <= n; j++)
// { // {
// x = 0; // x = 0;
// for (k = j; k <= i; k++) // for (k = j; k <= i; k++)
// x += hm.Get(i, k) * m2.Get(j, k); // x += hm.Get(i, k) * m2.Get(j, k);
// (*testout) << "i, j = " << i << "," << j << " x = " << x << endl; // (*testout) << "i, j = " << i << "," << j << " x = " << x << endl;
// } // }
// calc A^-1 = L^-T * L^-1 // calc A^-1 = L^-T * L^-1
@ -533,9 +530,9 @@ void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2)
for (i = 1; i <= n; i++) for (i = 1; i <= n; i++)
{ {
for (k = 1; k <= n; k++) for (k = 1; k <= n; k++)
hv.Elem(p.Get(k)) = m2.Get(i, k); hv(p.Get(k)-1) = m2.Get(i, k);
for (k = 1; k <= n; k++) for (k = 1; k <= n; k++)
m2.Elem(i, k) = hv.Get(k); m2.Elem(i, k) = hv(k-1);
} }
@ -614,12 +611,12 @@ void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2)
for (j = 1; j < i; j++) for (j = 1; j < i; j++)
m2.Elem(i, j) = m2.Elem(j, i); m2.Elem(i, j) = m2.Elem(j, i);
} }
*/ */
} }
} }
void CalcAAt (const DenseMatrix & a, DenseMatrix & m2) void CalcAAt (const DenseMatrix & a, DenseMatrix & m2)
{ {
int n1 = a.Height(); int n1 = a.Height();
int n2 = a.Width(); int n2 = a.Width();
@ -665,34 +662,7 @@ void CalcAAt (const DenseMatrix & a, DenseMatrix & m2)
#ifdef ABC void CalcAtA (const DenseMatrix & a, DenseMatrix & m2)
BaseMatrix * DenseMatrix :: InverseMatrix (const BitArray * /* inner */) const
{
if (Height() != Width())
{
(*myerr) << "BaseMatrix::InverseMatrix(): Matrix not symmetric" << endl;
return new DenseMatrix(1);
}
else
{
if (Symmetric())
{
(*mycout) << "Invmat not available" << endl;
BaseMatrix * invmat = NULL;
return invmat;
}
DenseMatrix * invmat = new DenseMatrix (Height());
CalcInverse (*this, *invmat);
return invmat;
}
}
#endif
void CalcAtA (const DenseMatrix & a, DenseMatrix & m2)
{ {
int n1 = a.Height(); int n1 = a.Height();
int n2 = a.Width(); int n2 = a.Width();
@ -717,10 +687,7 @@ void CalcAtA (const DenseMatrix & a, DenseMatrix & m2)
void CalcABt (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2)
void CalcABt (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2)
{ {
int n1 = a.Height(); int n1 = a.Height();
int n2 = a.Width(); int n2 = a.Width();
@ -759,7 +726,7 @@ void CalcABt (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2)
} }
void CalcAtB (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2) void CalcAtB (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2)
{ {
int n1 = a.Height(); int n1 = a.Height();
int n2 = a.Width(); int n2 = a.Width();
@ -805,7 +772,7 @@ void CalcAtB (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2)
DenseMatrix operator* (const DenseMatrix & m1, const DenseMatrix & m2) DenseMatrix operator* (const DenseMatrix & m1, const DenseMatrix & m2)
{ {
DenseMatrix temp (m1.Height(), m2.Width()); DenseMatrix temp (m1.Height(), m2.Width());
@ -825,7 +792,7 @@ DenseMatrix operator* (const DenseMatrix & m1, const DenseMatrix & m2)
} }
void Mult (const DenseMatrix & m1, const DenseMatrix & m2, DenseMatrix & m3) void Mult (const DenseMatrix & m1, const DenseMatrix & m2, DenseMatrix & m3)
{ {
double sum; double sum;
double *p1, *p1s, *p1sn, *p1snn, *p2, *p2s, *p2sn, *p3; double *p1, *p1s, *p1sn, *p1snn, *p2, *p2s, *p2sn, *p3;
@ -954,7 +921,7 @@ void Mult (const DenseMatrix & m1, const DenseMatrix & m2, DenseMatrix & m3)
DenseMatrix operator+ (const DenseMatrix & m1, const DenseMatrix & m2) DenseMatrix operator+ (const DenseMatrix & m1, const DenseMatrix & m2)
{ {
DenseMatrix temp (m1.Height(), m1.Width()); DenseMatrix temp (m1.Height(), m1.Width());
int i, j; int i, j;
@ -981,8 +948,8 @@ DenseMatrix operator+ (const DenseMatrix & m1, const DenseMatrix & m2)
void Transpose (const DenseMatrix & m1, DenseMatrix & m2) void Transpose (const DenseMatrix & m1, DenseMatrix & m2)
{ {
int w = m1.Width(); int w = m1.Width();
int h = m1.Height(); int h = m1.Height();
int i, j; int i, j;
@ -1000,11 +967,11 @@ void Transpose (const DenseMatrix & m1, DenseMatrix & m2)
pm1 += w; pm1 += w;
} }
} }
} }
/* /*
void DenseMatrix :: Mult (const Vector & v, Vector & prod) const void DenseMatrix :: Mult (const Vector & v, Vector & prod) const
{ {
double sum, val; double sum, val;
const double * mp, * sp; const double * mp, * sp;
@ -1019,7 +986,7 @@ void DenseMatrix :: Mult (const Vector & v, Vector & prod) const
if (prod.Size() != n) if (prod.Size() != n)
prod.SetSize (n); prod.SetSize (n);
#ifdef DEVELOP #ifdef DEVELOP
if (!n) if (!n)
{ {
cout << "DenseMatrix::Mult mheight = 0" << endl; cout << "DenseMatrix::Mult mheight = 0" << endl;
@ -1038,7 +1005,7 @@ void DenseMatrix :: Mult (const Vector & v, Vector & prod) const
(*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl; (*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl;
} }
else else
#endif #endif
{ {
if (Symmetric()) if (Symmetric())
{ {
@ -1087,10 +1054,10 @@ void DenseMatrix :: Mult (const Vector & v, Vector & prod) const
} }
} }
} }
*/ */
void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const
{ {
// const Vector & v = (const Vector&)bv; // .CastToVector(); // const Vector & v = (const Vector&)bv; // .CastToVector();
// Vector & prod = (Vector & )bprod; // .CastToVector(); // Vector & prod = (Vector & )bprod; // .CastToVector();
@ -1112,7 +1079,7 @@ void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const
prod.SetSize (w); prod.SetSize (w);
const double * pmat = &Get(1, 1); const double * pmat = &Get(1, 1);
const double * pv = &v.Get(1); const double * pv = &v(0);
prod = 0; prod = 0;
@ -1121,7 +1088,7 @@ void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const
double val = *pv; double val = *pv;
++pv; ++pv;
double * pprod = &prod.Elem(1); double * pprod = &prod(0);
for (j = w-1; j >= 0; --j, ++pmat, ++pprod) for (j = w-1; j >= 0; --j, ++pmat, ++pprod)
{ {
@ -1146,7 +1113,7 @@ void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const
} }
void DenseMatrix :: Residuum (const Vector & x, const Vector & b, void DenseMatrix :: Residuum (const Vector & x, const Vector & b,
Vector & res) const Vector & res) const
{ {
double sum; double sum;
@ -1166,26 +1133,25 @@ void DenseMatrix :: Residuum (const Vector & x, const Vector & b,
} }
else else
{ {
int i, j;
int h = Height(); int h = Height();
int w = Width(); int w = Width();
const double * mp = &Get(1, 1); const double * mp = &Get(1, 1);
for (i = 1; i <= h; i++) for (int i = 1; i <= h; i++)
{ {
sum = b.Get(i); sum = b(i-1);
const double * xp = &x.Get(1); const double * xp = &x(0);
for (j = 1; j <= w; ++j, ++mp, ++xp) for (int j = 1; j <= w; ++j, ++mp, ++xp)
sum -= *mp * *xp; sum -= *mp * *xp;
res.Elem(i) = sum; res(i-1) = sum;
} }
} }
} }
#ifdef ABC #ifdef ABC
double DenseMatrix :: EvaluateBilinearform (const Vector & hx) const double DenseMatrix :: EvaluateBilinearform (const Vector & hx) const
{ {
double sum = 0, hsum; double sum = 0, hsum;
// const Vector & hx = x.CastToVector(); // const Vector & hx = x.CastToVector();
@ -1208,12 +1174,12 @@ double DenseMatrix :: EvaluateBilinearform (const Vector & hx) const
} }
} }
// testout << "sum = " << sum << endl; // testout << "sum = " << sum << endl;
return sum; return sum;
} }
void DenseMatrix :: MultElementMatrix (const Array<int> & pnum, void DenseMatrix :: MultElementMatrix (const Array<int> & pnum,
const Vector & hx, Vector & hy) const Vector & hx, Vector & hy)
{ {
int i, j; int i, j;
@ -1239,7 +1205,7 @@ void DenseMatrix :: MultElementMatrix (const Array<int> & pnum,
} }
void DenseMatrix :: MultTransElementMatrix (const Array<int> & pnum, void DenseMatrix :: MultTransElementMatrix (const Array<int> & pnum,
const Vector & hx, Vector & hy) const Vector & hx, Vector & hy)
{ {
int i, j; int i, j;
@ -1256,14 +1222,14 @@ void DenseMatrix :: MultTransElementMatrix (const Array<int> & pnum,
#endif #endif
void DenseMatrix :: Solve (const Vector & v, Vector & sol) const void DenseMatrix :: Solve (const Vector & v, Vector & sol) const
{ {
DenseMatrix temp (*this); DenseMatrix temp (*this);
temp.SolveDestroy (v, sol); temp.SolveDestroy (v, sol);
} }
void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol) void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol)
{ {
double q; double q;
@ -1327,17 +1293,17 @@ void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol)
return; return;
} }
p.Elem(i) = 1 / sqrt(x); p(i-1) = 1 / sqrt(x);
} }
else else
{ {
Elem(j, i) = x * p.Get(i); Elem(j, i) = x * p(i-1);
} }
} }
} }
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
Elem(i, i) = 1 / p.Get(i); Elem(i, i) = 1 / p(i-1);
// A = L L^t // A = L L^t
// L stored in left-lower triangle // L stored in left-lower triangle
@ -1347,29 +1313,29 @@ void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol)
// Solve L sol = sol // Solve L sol = sol
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
{ {
double val = sol.Get(i); double val = sol(i-1);
const double * pij = &Get(i, 1); const double * pij = &Get(i, 1);
const double * solj = &sol.Get(1); const double * solj = &sol(0);
for (j = 1; j < i; j++, ++pij, ++solj) for (int j = 1; j < i; j++, ++pij, ++solj)
val -= *pij * *solj; val -= *pij * *solj;
// for (j = 1; j < i; j++) // for (j = 1; j < i; j++)
// val -= Get(i, j) * sol.Get(j); // val -= Get(i, j) * sol.Get(j);
sol.Elem(i) = val / Get(i, i); sol(i-1) = val / Get(i, i);
} }
// Solve L^t sol = sol // Solve L^t sol = sol
for (i = n; i >= 1; i--) for (int i = n; i >= 1; i--)
{ {
double val = sol.Get(i) / Get(i, i); double val = sol(i-1) / Get(i, i);
sol.Elem(i) = val; sol(i-1) = val;
double * solj = &sol.Elem(1); double * solj = &sol(0);
const double * pij = &Get(i, 1); const double * pij = &Get(i, 1);
for (j = 1; j < i; ++j, ++pij, ++solj) for (j = 1; j < i; ++j, ++pij, ++solj)
@ -1383,10 +1349,10 @@ void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol)
else else
{ {
// (*mycout) << "gauss" << endl; // (*mycout) << "gauss" << endl;
int i, j, k, n = Height(); int n = Height();
for (i = 1; i <= n; i++) for (int i = 1; i <= n; i++)
{ {
for (j = i+1; j <= n; j++) for (int j = i+1; j <= n; j++)
{ {
q = Get(j,i) / Get(i,i); q = Get(j,i) / Get(i,i);
if (q) if (q)
@ -1394,42 +1360,42 @@ void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol)
const double * pik = &Get(i, i+1); const double * pik = &Get(i, i+1);
double * pjk = &Elem(j, i+1); double * pjk = &Elem(j, i+1);
for (k = i+1; k <= n; ++k, ++pik, ++pjk) for (int k = i+1; k <= n; ++k, ++pik, ++pjk)
*pjk -= q * *pik; *pjk -= q * *pik;
// for (k = i+1; k <= Height(); k++) // for (k = i+1; k <= Height(); k++)
// Elem(j, k) -= q * Get(i,k); // Elem(j, k) -= q * Get(i,k);
sol.Elem(j) -= q * sol.Get(i); sol(j-1) -= q * sol(i-1);
} }
} }
} }
for (i = n; i >= 1; i--) for (int i = n; i >= 1; i--)
{ {
q = sol.Get(i); q = sol(i-1);
for (j = i+1; j <= n; j++) for (int j = i+1; j <= n; j++)
q -= Get(i,j) * sol.Get(j); q -= Get(i,j) * sol(j-1);
sol.Elem(i) = q / Get(i,i); sol(i-1) = q / Get(i,i);
} }
} }
} }
/* /*
BaseMatrix * DenseMatrix :: Copy () const BaseMatrix * DenseMatrix :: Copy () const
{ {
return new DenseMatrix (*this); return new DenseMatrix (*this);
} }
*/ */
ostream & operator<< (ostream & ost, const DenseMatrix & m) ostream & operator<< (ostream & ost, const DenseMatrix & m)
{ {
for (int i = 0; i < m.Height(); i++) for (int i = 0; i < m.Height(); i++)
{ {
for (int j = 0; j < m.Width(); j++) for (int j = 0; j < m.Width(); j++)
@ -1437,7 +1403,7 @@ ostream & operator<< (ostream & ost, const DenseMatrix & m)
ost << endl; ost << endl;
} }
return ost; return ost;
} }

21
libsrc/linalg/densemat.hpp
View File

@ -2,7 +2,7 @@
#define FILE_DENSEMAT #define FILE_DENSEMAT
/**************************************************************************/ /**************************************************************************/
/* File: densemat.hh */ /* File: densemat.hpp */
/* Author: Joachim Schoeberl */ /* Author: Joachim Schoeberl */
/* Date: 01. Oct. 94 */ /* Date: 01. Oct. 94 */
/**************************************************************************/ /**************************************************************************/
@ -12,11 +12,6 @@
*/ */
#include <assert.h>
class DenseMatrix class DenseMatrix
{ {
protected: protected:
@ -67,12 +62,8 @@ public:
#ifdef DEBUG #ifdef DEBUG
if (prod.Size() != height) if (prod.Size() != height)
{ {
cerr << "Mult: wrong vector size " << endl; (*myerr) << "Mult: wrong vector size " << endl;
assert (1);
// prod.SetSize (height);
} }
if (!height) if (!height)
{ {
cout << "DenseMatrix::Mult height = 0" << endl; cout << "DenseMatrix::Mult height = 0" << endl;
@ -94,13 +85,13 @@ public:
#endif #endif
{ {
mp = data; mp = data;
dp = &prod.Elem(1); dp = &prod(0);
for (int i = 1; i <= height; i++) for (int i = 0; i < height; i++)
{ {
sum = 0; sum = 0;
sp = &v.Get(1); sp = &v(0);
for (int j = 1; j <= width; j++) for (int j = 0; j < width; j++)
{ {
// sum += Get(i,j) * v.Get(j); // sum += Get(i,j) * v.Get(j);
sum += *mp * *sp; sum += *mp * *sp;

1
libsrc/linalg/linalg.hpp
View File

@ -10,7 +10,6 @@
/* /*
Data types for basic linear algebra Data types for basic linear algebra
more data types are found in linalgl.hpp
The basic concepts include the data types The basic concepts include the data types

12
libsrc/linalg/linopt.cpp
View File

@ -37,9 +37,9 @@ void LinearOptimize (const DenseMatrix & a, const Vector & b,
m.Elem(3, j) = a.Get(i3, j); m.Elem(3, j) = a.Get(i3, j);
} }
rs.Elem(1) = b.Get(i1); rs(0) = b(i1-1);
rs.Elem(2) = b.Get(i2); rs(1) = b(i2-1);
rs.Elem(3) = b.Get(i3); rs(2) = b(i3-1);
if (fabs (m.Det()) < 1e-12) continue; if (fabs (m.Det()) < 1e-12) continue;
@ -59,9 +59,9 @@ void LinearOptimize (const DenseMatrix & a, const Vector & b,
*/ */
double rmin = res.Elem(1); double rmin = res(0);
for (int hi = 2; hi <= res.Size(); hi++) for (int hi = 1; hi < res.Size(); hi++)
if (res.Elem(hi) < rmin) rmin = res.Elem(hi); if (res(hi) < rmin) rmin = res(hi);
if ( (f < fmin) && rmin >= -1e-8) if ( (f < fmin) && rmin >= -1e-8)
{ {

30
libsrc/linalg/linsearch.cpp
View File

@ -2,7 +2,7 @@
/* */ /* */
/* Problem: Liniensuche */ /* Problem: Liniensuche */
/* */ /* */
/* Programmautor: Joachim Schöberl */ /* Programmautor: Joachim Schöberl */
/* Matrikelnummer: 9155284 */ /* Matrikelnummer: 9155284 */
/* */ /* */
/* Algorithmus nach: */ /* Algorithmus nach: */
@ -98,36 +98,36 @@ void MinFunction :: ApproximateHesse (const Vector & x,
double eps = 1e-6; double eps = 1e-6;
double f, f11, f12, f21, f22; double f, f11, f12, f21, f22;
for (i = 1; i <= n; i++) for (i = 0; i < n; i++)
{ {
for (j = 1; j < i; j++) for (j = 0; j < i; j++)
{ {
hx = x; hx = x;
hx.Elem(i) = x.Get(i) + eps; hx(i) = x(i) + eps;
hx.Elem(j) = x.Get(j) + eps; hx(j) = x(j) + eps;
f11 = Func(hx); f11 = Func(hx);
hx.Elem(i) = x.Get(i) + eps; hx(i) = x(i) + eps;
hx.Elem(j) = x.Get(j) - eps; hx(j) = x(j) - eps;
f12 = Func(hx); f12 = Func(hx);
hx.Elem(i) = x.Get(i) - eps; hx(i) = x(i) - eps;
hx.Elem(j) = x.Get(j) + eps; hx(j) = x(j) + eps;
f21 = Func(hx); f21 = Func(hx);
hx.Elem(i) = x.Get(i) - eps; hx(i) = x(i) - eps;
hx.Elem(j) = x.Get(j) - eps; hx(j) = x(j) - eps;
f22 = Func(hx); f22 = Func(hx);
hesse.Elem(i, j) = hesse.Elem(j, i) = hesse(i, j) = hesse(j, i) =
(f11 + f22 - f12 - f21) / (2 * eps * eps); (f11 + f22 - f12 - f21) / (2 * eps * eps);
} }
hx = x; hx = x;
f = Func(x); f = Func(x);
hx.Elem(i) = x.Get(i) + eps; hx(i) = x(i) + eps;
f11 = Func(hx); f11 = Func(hx);
hx.Elem(i) = x.Get(i) - eps; hx(i) = x(i) - eps;
f22 = Func(hx); f22 = Func(hx);
hesse.Elem(i, i) = (f11 + f22 - 2 * f) / (eps * eps); hesse(i, i) = (f11 + f22 - 2 * f) / (eps * eps);
} }
// (*testout) << "hesse = " << hesse << endl; // (*testout) << "hesse = " << hesse << endl;
} }

18
libsrc/linalg/polynomial.cpp
View File

@ -152,24 +152,6 @@ MaxUnitSquare ()
if (hv > maxv) maxv = hv; if (hv > maxv) maxv = hv;
return maxv; return maxv;
// (*testout) << "maxv = " << maxv << " =~= ";
/*
double x, y;
maxv = -1e8;
double val;
for (x = 0; x <= 1.01; x += 0.1)
for (y = 0; y <= 1.01; y += 0.1)
{
val = Value (x, y);
if (val > maxv)
maxv = val;
}
// (*testout) << maxv << endl;
return maxv;
*/
}; };

230
libsrc/linalg/vector.cpp
View File

@ -6,34 +6,34 @@
namespace netgen namespace netgen
{ {
double BaseVector :: shit = 0; double BaseVector :: shit = 0;
// %FD Constructs a vector of length zero // %FD Constructs a vector of length zero
BaseVector :: BaseVector () BaseVector :: BaseVector ()
{ {
length = 0; length = 0;
} }
// %FD Constructs a vector of given length // %FD Constructs a vector of given length
BaseVector :: BaseVector ( BaseVector :: BaseVector (
INDEX alength // length of the vector INDEX alength // length of the vector
) )
{ {
length = alength; length = alength;
} }
// %FD Sets length of the vector, old vector will be destroyed // %FD Sets length of the vector, old vector will be destroyed
void void
BaseVector :: SetLength ( BaseVector :: SetLength (
INDEX alength // new length of the vector INDEX alength // new length of the vector
) )
{ {
length = alength; length = alength;
} }
// %FD Changes length of the vector, old values stay valid // %FD Changes length of the vector, old values stay valid
void void
BaseVector :: ChangeLength ( BaseVector :: ChangeLength (
INDEX alength // new length of the vector INDEX alength // new length of the vector
) )
{ {
@ -42,9 +42,9 @@ BaseVector :: ChangeLength (
// %FD { Write a vector with the help of the '<<' operator onto a stream } // %FD { Write a vector with the help of the '<<' operator onto a stream }
ostream & // stream for further use ostream & // stream for further use
operator<< ( operator<< (
ostream & s, // stream to write vector onto ostream & s, // stream to write vector onto
const BaseVector & v // vector to write const BaseVector & v // vector to write
) )
@ -53,10 +53,10 @@ operator<< (
} }
// %FD{ Divides every component of the vector by the scalar c. // %FD{ Divides every component of the vector by the scalar c.
// The function checks for division by zero } // The function checks for division by zero }
BaseVector & // result vector BaseVector & // result vector
BaseVector :: operator/= ( BaseVector :: operator/= (
double c // scalar to divide by double c // scalar to divide by
) )
{ {
@ -70,9 +70,9 @@ BaseVector :: operator/= (
} }
// %FD Creates a copy of the object // %FD Creates a copy of the object
BaseVector * // pointer to the new vector BaseVector * // pointer to the new vector
BaseVector :: Copy () const BaseVector :: Copy () const
{ {
cerr << "Base_vector::Copy called" << endl << flush; cerr << "Base_vector::Copy called" << endl << flush;
return NULL; return NULL;
@ -81,30 +81,30 @@ BaseVector :: Copy () const
void BaseVector :: GetElementVector (const Array<INDEX> & pnum, void BaseVector :: GetElementVector (const Array<INDEX> & pnum,
BaseVector & elvec) const BaseVector & elvec) const
{ {
int i; int i;
for (i = 1; i <= pnum.Size(); i++) for (i = 1; i <= pnum.Size(); i++)
elvec(i) = (*this)(pnum.Get(i)); elvec(i) = (*this)(pnum.Get(i));
} }
void BaseVector :: SetElementVector (const Array<INDEX> & pnum, void BaseVector :: SetElementVector (const Array<INDEX> & pnum,
const BaseVector & elvec) const BaseVector & elvec)
{ {
int i; int i;
for (i = 1; i <= pnum.Size(); i++) for (i = 1; i <= pnum.Size(); i++)
(*this)(pnum.Get(i)) = elvec(i); (*this)(pnum.Get(i)) = elvec(i);
} }
void BaseVector :: AddElementVector (const Array<INDEX> & pnum, void BaseVector :: AddElementVector (const Array<INDEX> & pnum,
const BaseVector & elvec) const BaseVector & elvec)
{ {
int i; int i;
for (i = 1; i <= pnum.Size(); i++) for (i = 1; i <= pnum.Size(); i++)
(*this)(pnum.Get(i)) += elvec(i); (*this)(pnum.Get(i)) += elvec(i);
} }
@ -116,39 +116,39 @@ void BaseVector :: AddElementVector (const Array<INDEX> & pnum,
TempVector :: ~TempVector () TempVector :: ~TempVector ()
{ {
delete vec; delete vec;
} }
TempVector BaseVector :: operator+ (const BaseVector & v2) const TempVector BaseVector :: operator+ (const BaseVector & v2) const
{ {
return (*Copy()) += v2; return (*Copy()) += v2;
} }
TempVector BaseVector :: operator- (const BaseVector & v2) const TempVector BaseVector :: operator- (const BaseVector & v2) const
{ {
return (*Copy()) -= v2; return (*Copy()) -= v2;
} }
TempVector BaseVector :: operator- () const TempVector BaseVector :: operator- () const
{ {
return (*Copy()) *= -1; return (*Copy()) *= -1;
} }
TempVector operator* (const BaseVector & v1, double scal) TempVector operator* (const BaseVector & v1, double scal)
{ {
return (*v1.Copy()) *= scal; return (*v1.Copy()) *= scal;
} }
TempVector operator/ (const BaseVector & v1, double scal) TempVector operator/ (const BaseVector & v1, double scal)
{ {
return (*v1.Copy()) /= scal; return (*v1.Copy()) /= scal;
} }
TempVector operator* (double scal, const BaseVector & v1) TempVector operator* (double scal, const BaseVector & v1)
{ {
return v1 * scal; return v1 * scal;
} }
@ -157,7 +157,7 @@ TempVector operator* (double scal, const BaseVector & v1)
BaseVector * TempVector :: Copy () const BaseVector * TempVector :: Copy () const
{ {
return vec->Copy(); return vec->Copy();
} }
@ -171,30 +171,30 @@ BaseVector * TempVector :: Copy () const
double Vector :: shit = 0; double Vector :: shit = 0;
class clVecpool class clVecpool
{ {
public: public:
Array<double *> vecs; Array<double *> vecs;
Array<INDEX> veclens; Array<INDEX> veclens;
~clVecpool(); ~clVecpool();
}; };
clVecpool :: ~clVecpool() clVecpool :: ~clVecpool()
{ {
int i; int i;
for (i = 1; i <= vecs.Size(); i++) for (i = 1; i <= vecs.Size(); i++)
delete vecs.Elem(i); delete vecs.Elem(i);
} }
static clVecpool vecpool; static clVecpool vecpool;
static double * NewDouble (INDEX len) static double * NewDouble (INDEX len)
{ {
if (len < 10) if (len < 10)
return new double[len]; return new double[len];
else else
@ -211,10 +211,10 @@ static double * NewDouble (INDEX len)
return new double[len]; return new double[len];
} }
} }
static void DeleteDouble (INDEX len, double * dp) static void DeleteDouble (INDEX len, double * dp)
{ {
if (len < 10) if (len < 10)
delete [] dp; delete [] dp;
else else
@ -222,16 +222,16 @@ static void DeleteDouble (INDEX len, double * dp)
vecpool.vecs.Append (dp); vecpool.vecs.Append (dp);
vecpool.veclens.Append (len); vecpool.veclens.Append (len);
} }
} }
Vector :: Vector () : BaseVector() Vector :: Vector () : BaseVector()
{ {
data = NULL; data = NULL;
} }
Vector :: Vector (INDEX alength) : BaseVector (alength) Vector :: Vector (INDEX alength) : BaseVector (alength)
{ {
if (length) if (length)
{ {
@ -249,7 +249,7 @@ Vector :: Vector (INDEX alength) : BaseVector (alength)
} }
Vector :: Vector (const Vector & v2) Vector :: Vector (const Vector & v2)
{ {
length = v2.length; length = v2.length;
@ -273,8 +273,8 @@ Vector :: Vector (const Vector & v2)
} }
Vector :: ~Vector () Vector :: ~Vector ()
{ {
// veclenfile << "~Vector delete: " << length << endl; // veclenfile << "~Vector delete: " << length << endl;
if (data) if (data)
{ {
@ -282,9 +282,9 @@ Vector :: ~Vector ()
// delete [] data; // delete [] data;
} }
} }
void Vector :: SetLength (INDEX alength) void Vector :: SetLength (INDEX alength)
{ {
if (length == alength) return; if (length == alength) return;
@ -307,8 +307,8 @@ void Vector :: SetLength (INDEX alength)
} }
} }
void Vector :: ChangeLength (INDEX alength) void Vector :: ChangeLength (INDEX alength)
{ {
(*mycout) << "Vector::ChangeLength called" << endl; (*mycout) << "Vector::ChangeLength called" << endl;
if (length == alength) return; if (length == alength) return;
@ -337,14 +337,14 @@ void Vector :: ChangeLength (INDEX alength)
length = alength; length = alength;
} }
/// NEW RM /// NEW RM
void Vector::SetBlockLength (INDEX /* blength */) void Vector::SetBlockLength (INDEX /* blength */)
{ {
MyError("BaseVector::SetBlockLength was called for a Vector"); MyError("BaseVector::SetBlockLength was called for a Vector");
} }
double & Vector :: operator() (INDEX i) double & Vector :: operator() (INDEX i)
{ {
if (i >= 1 && i <= length) return Elem(i); if (i >= 1 && i <= length) return Elem(i);
else (*myerr) << "\nindex " << i << " out of range (" else (*myerr) << "\nindex " << i << " out of range ("
@ -352,7 +352,7 @@ double & Vector :: operator() (INDEX i)
return shit; return shit;
} }
double Vector :: operator() (INDEX i) const double Vector :: operator() (INDEX i) const
{ {
if (i >= 1 && i <= length) return Get(i); if (i >= 1 && i <= length) return Get(i);
else (*myerr) << "\nindex " << i << " out of range (" else (*myerr) << "\nindex " << i << " out of range ("
@ -362,7 +362,7 @@ double Vector :: operator() (INDEX i) const
double Vector :: SupNorm () const double Vector :: SupNorm () const
{ {
double sup = 0; double sup = 0;
@ -373,7 +373,7 @@ double Vector :: SupNorm () const
return sup; return sup;
} }
double Vector :: L2Norm () const double Vector :: L2Norm () const
{ {
double sum = 0; double sum = 0;
@ -383,7 +383,7 @@ double Vector :: L2Norm () const
return sqrt (sum); return sqrt (sum);
} }
double Vector :: L1Norm () const double Vector :: L1Norm () const
{ {
double sum = 0; double sum = 0;
@ -393,7 +393,7 @@ double Vector :: L1Norm () const
return sum; return sum;
} }
double Vector :: Max () const double Vector :: Max () const
{ {
if (!Length()) return 0; if (!Length()) return 0;
double m = Get(1); double m = Get(1);
@ -402,7 +402,7 @@ double Vector :: Max () const
return m; return m;
} }
double Vector :: Min () const double Vector :: Min () const
{ {
if (!Length()) return 0; if (!Length()) return 0;
double m = Get(1); double m = Get(1);
@ -412,8 +412,8 @@ double Vector :: Min () const
} }
/* /*
ostream & operator<<(ostream & s, const Vector & v) ostream & operator<<(ostream & s, const Vector & v)
{ {
int w = s.width(); int w = s.width();
if (v.Length()) if (v.Length())
@ -434,9 +434,9 @@ ostream & operator<<(ostream & s, const Vector & v)
return s; return s;
} }
*/ */
ostream & Vector :: Print (ostream & s) const ostream & Vector :: Print (ostream & s) const
{ {
int w = s.width(); int w = s.width();
if (Length()) if (Length())
@ -460,7 +460,7 @@ ostream & Vector :: Print (ostream & s) const
BaseVector & Vector :: operator+= (const BaseVector & v2) BaseVector & Vector :: operator+= (const BaseVector & v2)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
@ -472,7 +472,7 @@ BaseVector & Vector :: operator+= (const BaseVector & v2)
return *this; return *this;
} }
BaseVector & Vector :: operator-= (const BaseVector & v2) BaseVector & Vector :: operator-= (const BaseVector & v2)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
@ -484,7 +484,7 @@ BaseVector & Vector :: operator-= (const BaseVector & v2)
return *this; return *this;
} }
BaseVector & Vector :: operator*= (double c) BaseVector & Vector :: operator*= (double c)
{ {
for (INDEX i = 1; i <= Length(); i++) for (INDEX i = 1; i <= Length(); i++)
Elem(i) *= c; Elem(i) *= c;
@ -493,7 +493,7 @@ BaseVector & Vector :: operator*= (double c)
BaseVector & Vector :: Add (double scal, const BaseVector & v2) BaseVector & Vector :: Add (double scal, const BaseVector & v2)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
@ -513,7 +513,7 @@ BaseVector & Vector :: Add (double scal, const BaseVector & v2)
return *this; return *this;
} }
BaseVector & Vector :: Add2 (double scal, const BaseVector & v2, BaseVector & Vector :: Add2 (double scal, const BaseVector & v2,
double scal3, const BaseVector & v3) double scal3, const BaseVector & v3)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
@ -536,7 +536,7 @@ BaseVector & Vector :: Add2 (double scal, const BaseVector & v2,
return *this; return *this;
} }
BaseVector & Vector :: Set (double scal, const BaseVector & v2) BaseVector & Vector :: Set (double scal, const BaseVector & v2)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
@ -557,9 +557,9 @@ BaseVector & Vector :: Set (double scal, const BaseVector & v2)
} }
BaseVector & Vector :: Set2 (double scal , const BaseVector & v2, BaseVector & Vector :: Set2 (double scal , const BaseVector & v2,
double scal3, const BaseVector & v3) double scal3, const BaseVector & v3)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
const Vector & hv3 = v3.CastToVector(); const Vector & hv3 = v3.CastToVector();
@ -578,11 +578,11 @@ BaseVector & Vector :: Set2 (double scal , const BaseVector & v2,
else else
(*myerr) << "Vector::Set: illegal dimension" << endl; (*myerr) << "Vector::Set: illegal dimension" << endl;
return *this; return *this;
} }
void Vector :: GetPart (int startpos, BaseVector & v2) const void Vector :: GetPart (int startpos, BaseVector & v2) const
{ {
Vector & hv2 = v2.CastToVector(); Vector & hv2 = v2.CastToVector();
if (Length() >= startpos + v2.Length() - 1) if (Length() >= startpos + v2.Length() - 1)
@ -594,12 +594,12 @@ void Vector :: GetPart (int startpos, BaseVector & v2) const
} }
else else
MyError ("Vector::GetPart: Vector to short"); MyError ("Vector::GetPart: Vector to short");
} }
// NEW RM // NEW RM
void Vector :: SetPart (int startpos, const BaseVector & v2) void Vector :: SetPart (int startpos, const BaseVector & v2)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
INDEX i; INDEX i;
INDEX n = v2.Length(); INDEX n = v2.Length();
@ -618,10 +618,10 @@ void Vector :: SetPart (int startpos, const BaseVector & v2)
} }
else else
MyError ("Vector::SetPart: Vector to short"); MyError ("Vector::SetPart: Vector to short");
} }
void Vector :: AddPart (int startpos, double val, const BaseVector & v2) void Vector :: AddPart (int startpos, double val, const BaseVector & v2)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
INDEX i; INDEX i;
INDEX n = v2.Length(); INDEX n = v2.Length();
@ -640,12 +640,12 @@ void Vector :: AddPart (int startpos, double val, const BaseVector & v2)
} }
else else
MyError ("Vector::AddPart: Vector to short"); MyError ("Vector::AddPart: Vector to short");
} }
double Vector :: operator* (const BaseVector & v2) const double Vector :: operator* (const BaseVector & v2) const
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
@ -666,7 +666,7 @@ double Vector :: operator* (const BaseVector & v2) const
return sum; return sum;
} }
void Vector :: SetRandom () void Vector :: SetRandom ()
{ {
INDEX i; INDEX i;
for (i = 1; i <= Length(); i++) for (i = 1; i <= Length(); i++)
@ -680,8 +680,8 @@ void Vector :: SetRandom ()
} }
/* /*
TempVector Vector :: operator- () const TempVector Vector :: operator- () const
{ {
Vector & sum = *(Vector*)Copy(); Vector & sum = *(Vector*)Copy();
@ -694,9 +694,9 @@ TempVector Vector :: operator- () const
(*myerr) << "operator+ (Vector, Vector): sum.Length() not ok" << endl; (*myerr) << "operator+ (Vector, Vector): sum.Length() not ok" << endl;
return sum; return sum;
} }
*/ */
BaseVector & Vector::operator= (const Vector & v2) BaseVector & Vector::operator= (const Vector & v2)
{ {
SetLength (v2.Length()); SetLength (v2.Length());
@ -710,7 +710,7 @@ BaseVector & Vector::operator= (const Vector & v2)
return *this; return *this;
} }
BaseVector & Vector::operator= (const BaseVector & v2) BaseVector & Vector::operator= (const BaseVector & v2)
{ {
const Vector & hv2 = v2.CastToVector(); const Vector & hv2 = v2.CastToVector();
@ -727,7 +727,7 @@ BaseVector & Vector::operator= (const BaseVector & v2)
} }
BaseVector & Vector::operator= (double scal) BaseVector & Vector::operator= (double scal)
{ {
if (!Length()) (*myerr) << "Vector::operator= (double) : data not allocated" if (!Length()) (*myerr) << "Vector::operator= (double) : data not allocated"
<< endl; << endl;
@ -739,13 +739,13 @@ BaseVector & Vector::operator= (double scal)
} }
BaseVector * Vector :: Copy () const BaseVector * Vector :: Copy () const
{ {
return new Vector (*this); return new Vector (*this);
} }
void Vector :: Swap (BaseVector & v2) void Vector :: Swap (BaseVector & v2)
{ {
Vector & hv2 = v2.CastToVector(); Vector & hv2 = v2.CastToVector();
swap (length, hv2.length); swap (length, hv2.length);
@ -755,32 +755,32 @@ void Vector :: Swap (BaseVector & v2)
void Vector :: GetElementVector (const Array<INDEX> & pnum, void Vector :: GetElementVector (const Array<INDEX> & pnum,
BaseVector & elvec) const BaseVector & elvec) const
{ {
int i; int i;
Vector & helvec = elvec.CastToVector(); Vector & helvec = elvec.CastToVector();
for (i = 1; i <= pnum.Size(); i++) for (i = 1; i <= pnum.Size(); i++)
helvec.Elem(i) = Get(pnum.Get(i)); helvec.Elem(i) = Get(pnum.Get(i));
} }
void Vector :: SetElementVector (const Array<INDEX> & pnum, void Vector :: SetElementVector (const Array<INDEX> & pnum,
const BaseVector & elvec) const BaseVector & elvec)
{ {
int i; int i;
const Vector & helvec = elvec.CastToVector(); const Vector & helvec = elvec.CastToVector();
for (i = 1; i <= pnum.Size(); i++) for (i = 1; i <= pnum.Size(); i++)
Elem(pnum.Get(i)) = helvec.Get(i); Elem(pnum.Get(i)) = helvec.Get(i);
} }
void Vector :: AddElementVector (const Array<INDEX> & pnum, void Vector :: AddElementVector (const Array<INDEX> & pnum,
const BaseVector & elvec) const BaseVector & elvec)
{ {
int i; int i;
const Vector & helvec = elvec.CastToVector(); const Vector & helvec = elvec.CastToVector();
for (i = 1; i <= pnum.Size(); i++) for (i = 1; i <= pnum.Size(); i++)
Elem(pnum.Get(i)) += helvec.Get(i); Elem(pnum.Get(i)) += helvec.Get(i);
} }
} }
#endif #endif

8
libsrc/linalg/vector.hpp
View File

@ -35,9 +35,9 @@ public:
double & operator() (int i) { return data[i]; } double & operator() (int i) { return data[i]; }
const double & operator() (int i) const { return data[i]; } const double & operator() (int i) const { return data[i]; }
double & Elem (int i) { return data[i-1]; } // double & Elem (int i) { return data[i-1]; }
const double & Get (int i) const { return data[i-1]; } // const double & Get (int i) const { return data[i-1]; }
void Set (int i, double val) { data[i-1] = val; } // void Set (int i, double val) { data[i-1] = val; }
FlatVector & operator*= (double scal) FlatVector & operator*= (double scal)
{ {
@ -92,7 +92,7 @@ public:
{ delete [] data; } { delete [] data; }
Vector & operator= (const FlatVector & v) Vector & operator= (const FlatVector & v)
{ memcpy (data, &v.Get(1), s*sizeof(double)); return *this; } { memcpy (data, &v(0), s*sizeof(double)); return *this; }
Vector & operator= (double scal) Vector & operator= (double scal)
{ {

28
libsrc/meshing/adfront3.cpp
View File

@ -793,25 +793,25 @@ bool AdFront3 :: Inside (const Point<3> & p) const
v1 = p2 - p1; v1 = p2 - p1;
v2 = p3 - p1; v2 = p3 - p1;
a.Elem(1, 1) = v1.X(); a(0, 0) = v1.X();
a.Elem(2, 1) = v1.Y(); a(1, 0) = v1.Y();
a.Elem(3, 1) = v1.Z(); a(2, 0) = v1.Z();
a.Elem(1, 2) = v2.X(); a(0, 1) = v2.X();
a.Elem(2, 2) = v2.Y(); a(1, 1) = v2.Y();
a.Elem(3, 2) = v2.Z(); a(2, 1) = v2.Z();
a.Elem(1, 3) = -n.X(); a(0, 2) = -n.X();
a.Elem(2, 3) = -n.Y(); a(1, 2) = -n.Y();
a.Elem(3, 3) = -n.Z(); a(2, 2) = -n.Z();
b.Elem(1) = p(0) - p1(0); b(0) = p(0) - p1(0);
b.Elem(2) = p(1) - p1(1); b(1) = p(1) - p1(1);
b.Elem(3) = p(2) - p1(2); b(2) = p(2) - p1(2);
CalcInverse (a, ainv); CalcInverse (a, ainv);
ainv.Mult (b, u); ainv.Mult (b, u);
if (u.Elem(1) >= 0 && u.Elem(2) >= 0 && u.Elem(1)+u.Elem(2) <= 1 && if (u(0) >= 0 && u(1) >= 0 && u(0)+u(1) <= 1 &&
u.Elem(3) > 0) u(2) > 0)
{ {
cnt++; cnt++;
} }

29
libsrc/meshing/bisect.cpp
View File

@ -1091,26 +1091,23 @@ namespace netgen
T_MPRISMS & mprisms, T_MPRISMS & mprisms,
const Mesh & mesh) const Mesh & mesh)
{ {
int i, j, k;
int step;
int marked = 0; int marked = 0;
int np = mesh.GetNP(); int np = mesh.GetNP();
Vector hv(np); Vector hv(np);
for (i = 1; i <= np; i++) for (int i = 0; i < np; i++)
hv.Elem(i) = mesh.GetH (mesh.Point(i)); hv(i) = mesh.GetH (mesh.Point(i+1));
double hfac = 1; double hfac = 1;
for (step = 1; step <= 2; step++) for (int step = 1; step <= 2; step++)
{ {
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
{ {
double h = 0; double h = 0;
for (j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
for (k = j+1; k < 4; k++) for (int k = j+1; k < 4; k++)
{ {
const Point<3> & p1 = mesh.Point (mtets.Get(i).pnums[j]); const Point<3> & p1 = mesh.Point (mtets.Get(i).pnums[j]);
const Point<3> & p2 = mesh.Point (mtets.Get(i).pnums[k]); const Point<3> & p2 = mesh.Point (mtets.Get(i).pnums[k]);
@ -1120,9 +1117,9 @@ namespace netgen
h = sqrt (h); h = sqrt (h);
double hshould = 1e10; double hshould = 1e10;
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
{ {
double hi = hv.Get (mtets.Get(i).pnums[j]); double hi = hv (mtets.Get(i).pnums[j]-1);
if (hi < hshould) if (hi < hshould)
hshould = hi; hshould = hi;
} }
@ -1145,12 +1142,12 @@ namespace netgen
} }
} }
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
{ {
double h = 0; double h = 0;
for (j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
for (k = j+1; k < 3; k++) for (int k = j+1; k < 3; k++)
{ {
const Point<3> & p1 = mesh.Point (mprisms.Get(i).pnums[j]); const Point<3> & p1 = mesh.Point (mprisms.Get(i).pnums[j]);
const Point<3> & p2 = mesh.Point (mprisms.Get(i).pnums[k]); const Point<3> & p2 = mesh.Point (mprisms.Get(i).pnums[k]);
@ -1160,9 +1157,9 @@ namespace netgen
h = sqrt (h); h = sqrt (h);
double hshould = 1e10; double hshould = 1e10;
for (j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {
double hi = hv.Get (mprisms.Get(i).pnums[j]); double hi = hv (mprisms.Get(i).pnums[j]-1);
if (hi < hshould) if (hi < hshould)
hshould = hi; hshould = hi;
} }

36
libsrc/meshing/hprefinement.cpp
View File

@ -117,7 +117,7 @@ namespace netgen
} }
void HPRefElement :: SetType( HPREF_ELEMENT_TYPE t) void HPRefElement :: SetType (HPREF_ELEMENT_TYPE t)
{ {
type = t; type = t;
switch(type) switch(type)
@ -129,11 +129,17 @@ namespace netgen
case HP_PRISM: np=6; break; case HP_PRISM: np=6; break;
case HP_PYRAMID: np=5; break; case HP_PYRAMID: np=5; break;
case HP_HEX: np=8; break; case HP_HEX: np=8; break;
default:
cerr << "HPRefElement: illegal type " << type << endl;
throw NgException ("HPRefElement::SetType: illegal type");
} }
for(int k=0;k<8;k++)
for(int k = 0; k < 8;k++)
{ {
pnums[k]=0; pnums[k]=0;
for(int l=0;l<3;l++) param[k][l]=0.; for(int l = 0; l < 3; l++)
param[k][l]=0.;
} }
} }
@ -587,6 +593,11 @@ namespace netgen
case PYRAMID: case PYRAMID:
hpel.type = HP_PYRAMID; hpel.type = HP_PYRAMID;
break; break;
default:
cerr << "HPRefElement: illegal elementtype (1) " << mesh[i].GetType() << endl;
throw NgException ("HPRefElement: illegal elementtype (1)");
} }
elements.Append(hpel); elements.Append(hpel);
} }
@ -700,7 +711,7 @@ namespace netgen
HPRef_Struct * hprs = Get_HPRef_Struct (el.type); HPRef_Struct * hprs = Get_HPRef_Struct (el.type);
int newlevel = el.levelx + 1; int newlevel = el.levelx + 1;
int oldnp(0); int oldnp = 0;
switch (hprs->geom) switch (hprs->geom)
{ {
case HP_SEGM: oldnp = 2; break; case HP_SEGM: oldnp = 2; break;
@ -710,6 +721,10 @@ namespace netgen
case HP_PYRAMID: oldnp = 5; break; case HP_PYRAMID: oldnp = 5; break;
case HP_PRISM: oldnp = 6; break; case HP_PRISM: oldnp = 6; break;
case HP_HEX: oldnp = 8; break; case HP_HEX: oldnp = 8; break;
default:
cerr << "HPRefElement: illegal type (3) " << hprs->geom << endl;
throw NgException ("HPRefElement::SetType: illegal type (3)");
} }
@ -830,6 +845,8 @@ namespace netgen
case HP_PRISM: newel.np=6; break; case HP_PRISM: newel.np=6; break;
case HP_TET: newel.np=4; break; case HP_TET: newel.np=4; break;
case HP_PYRAMID: newel.np=5; break; case HP_PYRAMID: newel.np=5; break;
default:
throw NgException (string("hprefinement.cpp: illegal type"));
} }
for (int k = 0; k < newel.np; k++) for (int k = 0; k < newel.np; k++)
@ -915,6 +932,11 @@ namespace netgen
case HP_PRISM: newel.np=6; break; case HP_PRISM: newel.np=6; break;
case HP_TET: newel.np=4; break; case HP_TET: newel.np=4; break;
case HP_PYRAMID: newel.np=5; break; case HP_PYRAMID: newel.np=5; break;
default:
cerr << "HPRefElement: illegal type (4) " << hprsnew->geom << endl;
throw NgException ("HPRefElement: illegal type (4)");
} }
newel.type = hprs->neweltypes[j]; newel.type = hprs->neweltypes[j];
for (int k = 0; k < 8; k++) for (int k = 0; k < 8; k++)
@ -1477,6 +1499,12 @@ namespace netgen
ord_dir[2] = 2; ord_dir[2] = 2;
ned = 8; ned = 8;
break; break;
default:
cerr << "HPRefElement: illegal elementtype (2) " << mesh[i].GetType() << endl;
throw NgException ("HPRefElement: illegal elementtype (2)");
} }
for (int j=0;j<ned;j++) for (int j=0;j<ned;j++)

166
libsrc/meshing/improve3.cpp
View File

@ -482,18 +482,18 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
pnew = Center (p1, p2); pnew = Center (p1, p2);
Vector px(3); Vector px(3);
px.Elem(1) = pnew.X(); px(0) = pnew.X();
px.Elem(2) = pnew.Y(); px(1) = pnew.Y();
px.Elem(3) = pnew.Z(); px(2) = pnew.Z();
if (elerrs[ei] > 0.1 * badmax) if (elerrs[ei] > 0.1 * badmax)
BFGS (px, pf, par); BFGS (px, pf, par);
bad2 = pf.Func (px); bad2 = pf.Func (px);
pnew.X() = px.Get(1); pnew.X() = px(0);
pnew.Y() = px.Get(2); pnew.Y() = px(1);
pnew.Z() = px.Get(3); pnew.Z() = px(2);
int hpinew = mesh.AddPoint (pnew); int hpinew = mesh.AddPoint (pnew);
@ -584,11 +584,6 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
const BitArray * working_elements) const BitArray * working_elements)
{ {
int j, k, l;
ElementIndex ei;
SurfaceElementIndex sei;
PointIndex pi1(0), pi2(0), pi3(0), pi4(0), pi5(0), pi6(0); PointIndex pi1(0), pi2(0), pi3(0), pi4(0), pi5(0), pi6(0);
int cnt = 0; int cnt = 0;
@ -600,7 +595,6 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
int np = mesh.GetNP(); int np = mesh.GetNP();
int ne = mesh.GetNE(); int ne = mesh.GetNE();
//int nse = mesh.GetNSE();
// contains at least all elements at node // contains at least all elements at node
TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np); TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np);
@ -614,16 +608,6 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
multithread.task = "Swap Improve"; multithread.task = "Swap Improve";
// mesh.CalcSurfacesOfNode (); // mesh.CalcSurfacesOfNode ();
/*
for (i = 1; i <= GetNE(); i++)
if (volelements.Get(i).PNum(1))
if (!LegalTet (volelements.Get(i)))
{
cout << "detected illegal tet, 1" << endl;
(*testout) << "detected illegal tet1: " << i << endl;
}
*/
INDEX_3_HASHTABLE<int> faces(mesh.GetNOpenElements()/3 + 2); INDEX_3_HASHTABLE<int> faces(mesh.GetNOpenElements()/3 + 2);
if (goal == OPT_CONFORM) if (goal == OPT_CONFORM)
@ -645,32 +629,15 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
} }
// find elements on node // find elements on node
for (ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
for (j = 0; j < mesh[ei].GetNP(); j++) for (int j = 0; j < mesh[ei].GetNP(); j++)
elementsonnode.Add (mesh[ei][j], ei); elementsonnode.Add (mesh[ei][j], ei);
/*
BitArray illegaltet(GetNE());
MarkIllegalElements();
if (goal == OPT_QUALITY || goal == OPT_LEGAL)
{
int cntill = 0;
for (i = 1; i <= GetNE(); i++)
{
// if (!LegalTet (volelements.Get(i)))
if (VolumeElement(i).flags.illegal)
{
cntill++;
illegaltet.Set (i);
}
}
// (*mycout) << cntill << " illegal tets" << endl;
}
*/
INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5); // INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5);
INDEX_2_CLOSED_HASHTABLE<int> edgeused(12 * ne + 5);
for (ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
{ {
if (multithread.terminate) if (multithread.terminate)
break; break;
@ -695,7 +662,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
// int onlybedges = 1; // int onlybedges = 1;
for (j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {
// loop over edges // loop over edges
@ -726,7 +693,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
edgeused.Set (i2, 1); edgeused.Set (i2, 1);
hasbothpoints.SetSize (0); hasbothpoints.SetSize (0);
for (k = 0; k < elementsonnode[pi1].Size(); k++) for (int k = 0; k < elementsonnode[pi1].Size(); k++)
{ {
bool has1 = 0, has2 = 0; bool has1 = 0, has2 = 0;
ElementIndex elnr = elementsonnode[pi1][k]; ElementIndex elnr = elementsonnode[pi1][k];
@ -734,7 +701,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
if (elem.IsDeleted()) continue; if (elem.IsDeleted()) continue;
for (l = 0; l < elem.GetNP(); l++) for (int l = 0; l < elem.GetNP(); l++)
{ {
if (elem[l] == pi1) has1 = 1; if (elem[l] == pi1) has1 = 1;
if (elem[l] == pi2) has2 = 1; if (elem[l] == pi2) has2 = 1;
@ -742,7 +709,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
if (has1 && has2) if (has1 && has2)
{ // only once { // only once
for (l = 0; l < hasbothpoints.Size(); l++) for (int l = 0; l < hasbothpoints.Size(); l++)
if (hasbothpoints[l] == elnr) if (hasbothpoints[l] == elnr)
has1 = 0; has1 = 0;
@ -752,7 +719,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
} }
bool puretet = 1; bool puretet = 1;
for (k = 0; k < hasbothpoints.Size(); k++) for (int k = 0; k < hasbothpoints.Size(); k++)
if (mesh[hasbothpoints[k]].GetType () != TET) if (mesh[hasbothpoints[k]].GetType () != TET)
puretet = 0; puretet = 0;
if (!puretet) if (!puretet)
@ -763,7 +730,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
if ( nsuround == 3 ) if ( nsuround == 3 )
{ {
Element & elem = mesh[hasbothpoints[0]]; Element & elem = mesh[hasbothpoints[0]];
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elem[l] != pi1 && elem[l] != pi2) if (elem[l] != pi1 && elem[l] != pi2)
{ {
pi4 = pi3; pi4 = pi3;
@ -784,16 +751,16 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
} }
pi5 = 0; pi5 = 0;
for (k = 1; k < 3; k++) for (int k = 1; k < 3; k++)
{ {
const Element & elemk = mesh[hasbothpoints[k]]; const Element & elemk = mesh[hasbothpoints[k]];
bool has1 = 0; bool has1 = 0;
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elemk[l] == pi4) if (elemk[l] == pi4)
has1 = 1; has1 = 1;
if (has1) if (has1)
{ {
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elemk[l] != pi1 && elemk[l] != pi2 && elemk[l] != pi4) if (elemk[l] != pi1 && elemk[l] != pi2 && elemk[l] != pi4)
pi5 = elemk[l]; pi5 = elemk[l];
} }
@ -905,21 +872,20 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
el22.flags.illegal_valid = 0; el22.flags.illegal_valid = 0;
mesh[hasbothpoints[0]] = el21; mesh[hasbothpoints[0]] = el21;
mesh[hasbothpoints[1]] = el22; mesh[hasbothpoints[1]] = el22;
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
mesh[hasbothpoints[2]][l] = 0; mesh[hasbothpoints[2]][l] = 0;
mesh[hasbothpoints[2]].Delete(); mesh[hasbothpoints[2]].Delete();
for (k = 0; k < 2; k++) for (int k = 0; k < 2; k++)
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]); elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]);
} }
} }
if (nsuround == 4) if (nsuround == 4)
{ {
const Element & elem1 = mesh[hasbothpoints[0]]; const Element & elem1 = mesh[hasbothpoints[0]];
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elem1[l] != pi1 && elem1[l] != pi2) if (elem1[l] != pi1 && elem1[l] != pi2)
{ {
pi4 = pi3; pi4 = pi3;
@ -938,32 +904,32 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
} }
pi5 = 0; pi5 = 0;
for (k = 1; k < 4; k++) for (int k = 1; k < 4; k++)
{ {
const Element & elem = mesh[hasbothpoints[k]]; const Element & elem = mesh[hasbothpoints[k]];
bool has1 = 0; bool has1 = 0;
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elem[l] == pi4) if (elem[l] == pi4)
has1 = 1; has1 = 1;
if (has1) if (has1)
{ {
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elem[l] != pi1 && elem[l] != pi2 && elem[l] != pi4) if (elem[l] != pi1 && elem[l] != pi2 && elem[l] != pi4)
pi5 = elem[l]; pi5 = elem[l];
} }
} }
pi6 = 0; pi6 = 0;
for (k = 1; k < 4; k++) for (int k = 1; k < 4; k++)
{ {
const Element & elem = mesh[hasbothpoints[k]]; const Element & elem = mesh[hasbothpoints[k]];
bool has1 = 0; bool has1 = 0;
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elem[l] == pi3) if (elem[l] == pi3)
has1 = 1; has1 = 1;
if (has1) if (has1)
{ {
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elem[l] != pi1 && elem[l] != pi2 && elem[l] != pi3) if (elem[l] != pi1 && elem[l] != pi2 && elem[l] != pi3)
pi6 = elem[l]; pi6 = elem[l];
} }
@ -1183,8 +1149,8 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
mesh[hasbothpoints[2]] = el3; mesh[hasbothpoints[2]] = el3;
mesh[hasbothpoints[3]] = el4; mesh[hasbothpoints[3]] = el4;
for (k = 0; k < 4; k++) for (int k = 0; k < 4; k++)
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]); elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]);
} }
else if (swap3) else if (swap3)
@ -1216,8 +1182,8 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
mesh[hasbothpoints[3]] = el4b; mesh[hasbothpoints[3]] = el4b;
for (k = 0; k < 4; k++) for (int k = 0; k < 4; k++)
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]); elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]);
} }
} }
@ -1231,7 +1197,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
Element & elem = mesh[hasbothpoints[0]]; Element & elem = mesh[hasbothpoints[0]];
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elem[l] != pi1 && elem[l] != pi2) if (elem[l] != pi1 && elem[l] != pi2)
{ {
pi4 = pi3; pi4 = pi3;
@ -1259,12 +1225,12 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
tetused = 0; tetused = 0;
tetused[0] = 1; tetused[0] = 1;
for (l = 2; l < nsuround; l++) for (int l = 2; l < nsuround; l++)
{ {
int oldpi = suroundpts[l-1]; int oldpi = suroundpts[l-1];
int newpi = 0; int newpi = 0;
for (k = 0; k < nsuround && !newpi; k++) for (int k = 0; k < nsuround && !newpi; k++)
if (!tetused[k]) if (!tetused[k])
{ {
const Element & nel = mesh[hasbothpoints[k]]; const Element & nel = mesh[hasbothpoints[k]];
@ -1284,7 +1250,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
bad1 = 0; bad1 = 0;
for (k = 0; k < nsuround; k++) for (int k = 0; k < nsuround; k++)
{ {
hel[0] = pi1; hel[0] = pi1;
hel[1] = pi2; hel[1] = pi2;
@ -1303,11 +1269,11 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
int confedge = -1; int confedge = -1;
double badopt = bad1; double badopt = bad1;
for (l = 0; l < nsuround; l++) for (int l = 0; l < nsuround; l++)
{ {
bad2 = 0; bad2 = 0;
for (k = l+1; k <= nsuround + l - 2; k++) for (int k = l+1; k <= nsuround + l - 2; k++)
{ {
hel[0] = suroundpts[l]; hel[0] = suroundpts[l];
hel[1] = suroundpts[k % nsuround]; hel[1] = suroundpts[k % nsuround];
@ -1344,7 +1310,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
(bad2 <= 100 * bad1 && bad2 <= 1e18) || (bad2 <= 100 * bad1 && bad2 <= 1e18) ||
(bad2 <= 1e8); (bad2 <= 1e8);
for (k = l+1; k <= nsuround + l - 2; k++) for (int k = l+1; k <= nsuround + l - 2; k++)
{ {
INDEX_3 hi3(suroundpts[l], INDEX_3 hi3(suroundpts[l],
suroundpts[k % nsuround], suroundpts[k % nsuround],
@ -1359,7 +1325,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
} }
} }
for (k = l+2; k <= nsuround+l-2; k++) for (int k = l+2; k <= nsuround+l-2; k++)
{ {
if (mesh.BoundaryEdge (suroundpts[l], if (mesh.BoundaryEdge (suroundpts[l],
suroundpts[k % nsuround])) suroundpts[k % nsuround]))
@ -1385,7 +1351,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
// (*mycout) << nsuround << "->" << 2 * (nsuround-2) << " " << flush; // (*mycout) << nsuround << "->" << 2 * (nsuround-2) << " " << flush;
cnt++; cnt++;
for (k = bestl+1; k <= nsuround + bestl - 2; k++) for (int k = bestl+1; k <= nsuround + bestl - 2; k++)
{ {
int k1; int k1;
@ -1420,7 +1386,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
elementsonnode.Add (hel[k1], mesh.GetNE()-1); elementsonnode.Add (hel[k1], mesh.GetNE()-1);
} }
for (k = 0; k < nsuround; k++) for (int k = 0; k < nsuround; k++)
{ {
Element & rel = mesh[hasbothpoints[k]]; Element & rel = mesh[hasbothpoints[k]];
/* /*
@ -1513,8 +1479,6 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
} }
} }
ElementIndex ei;
SurfaceElementIndex sei;
PointIndex pi1, pi2, pi3, pi4, pi5, pi6; PointIndex pi1, pi2, pi3, pi4, pi5, pi6;
PointIndex pi1other, pi2other; PointIndex pi1other, pi2other;
@ -1548,11 +1512,11 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
// find elements on node // find elements on node
for (ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
for (int j = 0; j < mesh[ei].GetNP(); j++) for (int j = 0; j < mesh[ei].GetNP(); j++)
elementsonnode.Add (mesh[ei][j], ei); elementsonnode.Add (mesh[ei][j], ei);
for (sei = 0; sei < nse; sei++) for (SurfaceElementIndex sei = 0; sei < nse; sei++)
for(int j=0; j<mesh[sei].GetNP(); j++) for(int j=0; j<mesh[sei].GetNP(); j++)
{ {
surfaceelementsonnode.Add(mesh[sei][j], sei); surfaceelementsonnode.Add(mesh[sei][j], sei);
@ -1563,9 +1527,10 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
bool periodic; bool periodic;
int idnum(-1); int idnum(-1);
INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5); // INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5);
INDEX_2_CLOSED_HASHTABLE<int> edgeused(12 * ne + 5);
for (ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
{ {
if (multithread.terminate) if (multithread.terminate)
break; break;
@ -2354,18 +2319,12 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
*/ */
void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
{ {
int j, k, l;
ElementIndex ei, eli1, eli2, elnr;
SurfaceElementIndex sei;
PointIndex pi1(0), pi2(0), pi3(0), pi4(0), pi5(0); PointIndex pi1(0), pi2(0), pi3(0), pi4(0), pi5(0);
int cnt = 0;
Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
int cnt = 0;
double bad1, bad2; double bad1, bad2;
int np = mesh.GetNP(); int np = mesh.GetNP();
@ -2383,7 +2342,6 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
// TestOk(); // TestOk();
/* /*
CalcSurfacesOfNode (); CalcSurfacesOfNode ();
for (i = 1; i <= GetNE(); i++) for (i = 1; i <= GetNE(); i++)
@ -2404,15 +2362,15 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
// find elements on node // find elements on node
for (ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
for (j = 0; j < mesh[ei].GetNP(); j++) for (int j = 0; j < mesh[ei].GetNP(); j++)
elementsonnode.Add (mesh[ei][j], ei); elementsonnode.Add (mesh[ei][j], ei);
for (sei = 0; sei < nse; sei++) for (SurfaceElementIndex sei = 0; sei < nse; sei++)
for (j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
belementsonnode.Add (mesh[sei][j], sei); belementsonnode.Add (mesh[sei][j], sei);
for (eli1 = 0; eli1 < ne; eli1++) for (ElementIndex eli1 = 0; eli1 < ne; eli1++)
{ {
if (multithread.terminate) if (multithread.terminate)
break; break;
@ -2431,7 +2389,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
// cout << "eli = " << eli1 << endl; // cout << "eli = " << eli1 << endl;
// (*testout) << "swapimp2, eli = " << eli1 << "; el = " << mesh[eli1] << endl; // (*testout) << "swapimp2, eli = " << eli1 << "; el = " << mesh[eli1] << endl;
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
{ {
// loop over faces // loop over faces
@ -2463,13 +2421,13 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
bool bface = 0; bool bface = 0;
for (k = 0; k < belementsonnode[pi1].Size(); k++) for (int k = 0; k < belementsonnode[pi1].Size(); k++)
{ {
const Element2d & bel = const Element2d & bel =
mesh[belementsonnode[pi1][k]]; mesh[belementsonnode[pi1][k]];
bool bface1 = 1; bool bface1 = 1;
for (l = 0; l < 3; l++) for (int l = 0; l < 3; l++)
if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3) if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3)
{ {
bface1 = 0; bface1 = 0;
@ -2487,9 +2445,9 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
FlatArray<ElementIndex> row = elementsonnode[pi1]; FlatArray<ElementIndex> row = elementsonnode[pi1];
for (k = 0; k < row.Size(); k++) for (int k = 0; k < row.Size(); k++)
{ {
eli2 = row[k]; ElementIndex eli2 = row[k];
// cout << "\rei1 = " << eli1 << ", pi1 = " << pi1 << ", k = " << k << ", ei2 = " << eli2 // cout << "\rei1 = " << eli1 << ", pi1 = " << pi1 << ", k = " << k << ", ei2 = " << eli2
// << ", getne = " << mesh.GetNE(); // << ", getne = " << mesh.GetNE();
@ -2502,7 +2460,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
continue; continue;
int comnodes=0; int comnodes=0;
for (l = 1; l <= 4; l++) for (int l = 1; l <= 4; l++)
if (elem2.PNum(l) == pi1 || elem2.PNum(l) == pi2 || if (elem2.PNum(l) == pi1 || elem2.PNum(l) == pi2 ||
elem2.PNum(l) == pi3) elem2.PNum(l) == pi3)
{ {
@ -2586,7 +2544,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
} }
*/ */
// cout << "neli = " << neli << endl; // cout << "neli = " << neli << endl;
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
{ {
elementsonnode.Add (el31[l], eli1); elementsonnode.Add (el31[l], eli1);
elementsonnode.Add (el32[l], eli2); elementsonnode.Add (el32[l], eli2);

14
libsrc/meshing/improve3.hpp
View File

@ -20,8 +20,18 @@ public:
extern double CalcBad (const Mesh::T_POINTS & points, const Element & elem, inline double
double h); CalcBad (const Mesh::T_POINTS & points, const Element & elem,
double h)
{
if (elem.GetType() == TET)
return CalcTetBadness (points[elem[0]], points[elem[1]],
points[elem[2]], points[elem[3]], h);
return 0;
}
extern double CalcTotalBad (const Mesh::T_POINTS & points, extern double CalcTotalBad (const Mesh::T_POINTS & points,
const Mesh::T_VOLELEMENTS & elements); const Mesh::T_VOLELEMENTS & elements);

6
libsrc/meshing/meshtool.cpp
View File

@ -202,7 +202,7 @@ namespace netgen
Vec3d v2 (p1, p3); Vec3d v2 (p1, p3);
Vec3d v3 (p1, p4); Vec3d v3 (p1, p4);
vol = -Determinant (v1, v2, v3) / 6; vol = Determinant (v1, v2, v3) * (-0.166666666666666);
ll1 = v1.Length2(); ll1 = v1.Length2();
ll2 = v2.Length2(); ll2 = v2.Length2();
@ -276,7 +276,7 @@ namespace netgen
Vec3d v5 (*pp2, *pp4); Vec3d v5 (*pp2, *pp4);
Vec3d v6 (*pp3, *pp4); Vec3d v6 (*pp3, *pp4);
vol = -Determinant (v1, v2, v3) / 6; vol = Determinant (v1, v2, v3) * (-0.166666666666666);
Vec3d gradvol; Vec3d gradvol;
Cross (v5, v4, gradvol); Cross (v5, v4, gradvol);
@ -344,10 +344,8 @@ namespace netgen
graderr += (1/(h*h) - h*h/(ll1*ll1)) * gradll1; graderr += (1/(h*h) - h*h/(ll1*ll1)) * gradll1;
graderr += (1/(h*h) - h*h/(ll2*ll2)) * gradll2; graderr += (1/(h*h) - h*h/(ll2*ll2)) * gradll2;
graderr += (1/(h*h) - h*h/(ll3*ll3)) * gradll3; graderr += (1/(h*h) - h*h/(ll3*ll3)) * gradll3;
cout << "?";
} }
double errpow; double errpow;
if (teterrpow == 2) if (teterrpow == 2)
{ {

625
libsrc/meshing/meshtype.cpp
View File

File diff suppressed because it is too large Load Diff

8
libsrc/meshing/netrule2.cpp
View File

@ -128,10 +128,10 @@ int netrule :: IsLineInFreeZone2 (const Point2d & p1, const Point2d & p2) const
{ {
int left, right, allleft, allright; int left, right, allleft, allright;
if (p1.X() > fzmaxx && p2.X() > fzmaxx || if ( (p1.X() > fzmaxx && p2.X() > fzmaxx) ||
p1.X() < fzminx && p2.X() < fzminx || (p1.X() < fzminx && p2.X() < fzminx) ||
p1.Y() > fzmaxy && p2.Y() > fzmaxy || (p1.Y() > fzmaxy && p2.Y() > fzmaxy) ||
p1.Y() < fzminy && p2.Y() < fzminy) return 0; (p1.Y() < fzminy && p2.Y() < fzminy) ) return 0;
for (int i = 1; i <= transfreezone.Size(); i++) for (int i = 1; i <= transfreezone.Size(); i++)
{ {

6
libsrc/meshing/netrule3.cpp
View File

@ -61,7 +61,7 @@ void vnetrule :: SetFreeZoneTransformation (const Vector & allp, int tolclass)
for (i = 1; i <= 3; i++) for (i = 1; i <= 3; i++)
{ {
for (j = 1; j <= np; j++) for (j = 1; j <= np; j++)
vp.Elem(j) = allp.Get(i+3*j-3); vp(j-1) = allp(i+3*j-3-1);
oldutofreezone->Mult (vp, vfp1); oldutofreezone->Mult (vp, vfp1);
oldutofreezonelimit->Mult (vp, vfp2); oldutofreezonelimit->Mult (vp, vfp2);
@ -70,7 +70,7 @@ void vnetrule :: SetFreeZoneTransformation (const Vector & allp, int tolclass)
vfp1.Add (lam2, vfp2); vfp1.Add (lam2, vfp2);
for (j = 1; j <= nfp; j++) for (j = 1; j <= nfp; j++)
transfreezone.Elem(j).X(i) = vfp1.Elem(j); transfreezone.Elem(j).X(i) = vfp1(j-1);
} }
// MARK(setfz2); // MARK(setfz2);
@ -814,7 +814,7 @@ int vnetrule :: IsTriangleInFreeSet (const Point3d & p1, const Point3d & p2,
if (it > minit) if (it > minit)
(*testout) << "act1,2,3 = " << act1 << act2 << act3 << endl; (*testout) << "act1,2,3 = " << act1 << act2 << act3 << endl;
if (act1 && act2 || act1 && act3 || act2 && act3) return 0; if ( (act1 && act2) || (act1 && act3) || (act2 && act3) ) return 0;
if (act1) if (act1)
{ {

10
libsrc/meshing/parser3.cpp
View File

@ -706,13 +706,13 @@ void vnetrule :: LoadRule (istream & ist)
if (quality < 100) if (quality < 100)
{ {
for (i = 1; i <= 3; i++) for (int i = 1; i <= 3; i++)
{ {
for (j = 1; j <= points.Size(); j++) for (int j = 1; j <= points.Size(); j++)
vp.Elem(j) = points.Get(j).X(i); vp(j-1) = points.Get(j).X(i);
oldutofreezone->Mult(vp, vfp); oldutofreezone->Mult(vp, vfp);
for (j = 1; j <= freezone.Size(); j++) for (int j = 1; j <= freezone.Size(); j++)
freezone.Elem(j).X(i) = vfp.Get(j); freezone.Elem(j).X(i) = vfp(j-1);
} }
// for (i = 1; i <= freezone.Size(); i++) // for (i = 1; i <= freezone.Size(); i++)
// (*testout) << "freepoint: " << freezone.Get(i) << endl; // (*testout) << "freepoint: " << freezone.Get(i) << endl;

10
libsrc/meshing/ruler2.cpp
View File

@ -397,11 +397,11 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
{ {
oldu.SetSize (2 * rule->GetNOldP()); oldu.SetSize (2 * rule->GetNOldP());
for (i = 1; i <= rule->GetNOldP(); i++) for (int i = 1; i <= rule->GetNOldP(); i++)
{ {
Vec2d ui(rule->GetPoint(i), lpoints.Get(pmap.Get(i))); Vec2d ui(rule->GetPoint(i), lpoints.Get(pmap.Get(i)));
oldu.Set (2*i-1, ui.X()); oldu (2*i-2) = ui.X();
oldu.Set (2*i , ui.Y()); oldu (2*i-1) = ui.Y();
} }
rule -> SetFreeZoneTransformation (oldu, tolerance); rule -> SetFreeZoneTransformation (oldu, tolerance);
@ -513,8 +513,8 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
for (i = oldnp + 1; i <= rule->GetNP(); i++) for (i = oldnp + 1; i <= rule->GetNP(); i++)
{ {
np = rule->GetPoint(i); np = rule->GetPoint(i);
np.X() += newu.Elem (2 * (i-oldnp) - 1); np.X() += newu (2 * (i-oldnp) - 2);
np.Y() += newu.Elem (2 * (i-oldnp)); np.Y() += newu (2 * (i-oldnp) - 1);
pmap.Elem(i) = lpoints.Append (np); pmap.Elem(i) = lpoints.Append (np);
} }

24
libsrc/meshing/ruler3.cpp
View File

@ -603,13 +603,13 @@ int Meshing3 :: ApplyRules
{ {
const Point3d & lp = lpoints.Get(pmap.Get(i)); const Point3d & lp = lpoints.Get(pmap.Get(i));
const Point3d & rp = rule->GetPoint(i); const Point3d & rp = rule->GetPoint(i);
oldu.Set (3*i-2, lp.X()-rp.X()); oldu (3*i-3) = lp.X()-rp.X();
oldu.Set (3*i-1, lp.Y()-rp.Y()); oldu (3*i-2) = lp.Y()-rp.Y();
oldu.Set (3*i , lp.Z()-rp.Z()); oldu (3*i-1) = lp.Z()-rp.Z();
allp.Set (3*i-2, lp.X()); allp (3*i-3) = lp.X();
allp.Set (3*i-1, lp.Y()); allp (3*i-2) = lp.Y();
allp.Set (3*i , lp.Z()); allp (3*i-1) = lp.Z();
} }
if (rule->GetNP() > rule->GetNOldP()) if (rule->GetNP() > rule->GetNOldP())
@ -623,9 +623,9 @@ int Meshing3 :: ApplyRules
for (i = rule->GetNOldP()+1; i <= rule->GetNP(); i++) for (i = rule->GetNOldP()+1; i <= rule->GetNP(); i++)
{ {
const Point3d & rp = rule->GetPoint(i); const Point3d & rp = rule->GetPoint(i);
allp.Set (3*i-2, rp.X() + newu.Get(3*i-2 - idiff)); allp (3*i-3) = rp.X() + newu(3*i-3 - idiff);
allp.Set (3*i-1, rp.Y() + newu.Get(3*i-1 - idiff)); allp (3*i-2) = rp.Y() + newu(3*i-2 - idiff);
allp.Set (3*i , rp.Z() + newu.Get(3*i - idiff)); allp (3*i-1) = rp.Z() + newu(3*i-1 - idiff);
} }
rule->SetFreeZoneTransformation (allp, rule->SetFreeZoneTransformation (allp,
@ -865,9 +865,9 @@ int Meshing3 :: ApplyRules
for (i = oldnp + 1; i <= rule->GetNP(); i++) for (i = oldnp + 1; i <= rule->GetNP(); i++)
{ {
np = rule->GetPoint(i); np = rule->GetPoint(i);
np.X() += newu.Elem (3 * (i-oldnp) - 2); np.X() += newu (3 * (i-oldnp) - 3);
np.Y() += newu.Elem (3 * (i-oldnp) - 1); np.Y() += newu (3 * (i-oldnp) - 2);
np.Z() += newu.Elem (3 * (i-oldnp)); np.Z() += newu (3 * (i-oldnp) - 1);
pmap.Elem(i) = lpoints.Append (np); pmap.Elem(i) = lpoints.Append (np);
} }

6
libsrc/meshing/smoothing2.5.cpp
View File

@ -219,9 +219,9 @@ namespace netgen
while (loci <= 5 && !moveisok) while (loci <= 5 && !moveisok)
{ {
loci ++; loci ++;
mesh[pi](0) = origp(0) + x.Get(1)*fact; mesh[pi](0) = origp(0) + x(0)*fact;
mesh[pi](1) = origp(1) + x.Get(2)*fact; mesh[pi](1) = origp(1) + x(1)*fact;
mesh[pi](2) = origp(2) + x.Get(3)*fact; mesh[pi](2) = origp(2) + x(2)*fact;
fact = fact/2.; fact = fact/2.;

10
libsrc/meshing/smoothing2.cpp
View File

@ -441,7 +441,7 @@ namespace netgen
vgrad = 0.0; vgrad = 0.0;
badness = 0; badness = 0;
pp1 = sp1 + x.Get(1) * t1; pp1 = sp1 + x(0) * t1;
meshthis -> ProjectPoint2 (surfi, surfi2, pp1); meshthis -> ProjectPoint2 (surfi, surfi2, pp1);
for (j = 0; j < locelements.Size(); j++) for (j = 0; j < locelements.Size(); j++)
@ -540,7 +540,7 @@ namespace netgen
pts2d.Elem(pi) = Point2d (t1 * (mesh.Point(pi) - sp1), pts2d.Elem(pi) = Point2d (t1 * (mesh.Point(pi) - sp1),
t2 * (mesh.Point(pi) - sp1)); t2 * (mesh.Point(pi) - sp1));
} }
pts2d.Elem(gpi) = Point2d (x.Get(1), x.Get(2)); pts2d.Elem(gpi) = Point2d (x(0), x(1));
for (int k = 1; k <= 2; k++) for (int k = 1; k <= 2; k++)
@ -553,7 +553,7 @@ namespace netgen
hbad = bel. hbad = bel.
CalcJacobianBadnessDirDeriv (pts2d, lpi, vdir, hderiv); CalcJacobianBadnessDirDeriv (pts2d, lpi, vdir, hderiv);
grad.Elem(k) += hderiv; grad(k-1) += hderiv;
if (k == 1) if (k == 1)
badness += hbad; badness += hbad;
} }
@ -590,7 +590,7 @@ namespace netgen
// pp1 = sp1; // pp1 = sp1;
// pp1.Add2 (x.Get(1), t1, x.Get(2), t2); // pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
pp1 = sp1 + x.Get(1) * t1 + x.Get(2) * t2; pp1 = sp1 + x(0) * t1 + x(1) * t2;
static Array<Point2d> pts2d; static Array<Point2d> pts2d;
pts2d.SetSize(mesh.GetNP()); pts2d.SetSize(mesh.GetNP());
@ -611,7 +611,7 @@ namespace netgen
pts2d.Elem(pi) = Point2d (t1 * (mesh.Point(pi) - sp1), pts2d.Elem(pi) = Point2d (t1 * (mesh.Point(pi) - sp1),
t2 * (mesh.Point(pi) - sp1)); t2 * (mesh.Point(pi) - sp1));
} }
pts2d.Elem(gpi) = Point2d (x.Get(1), x.Get(2)); pts2d.Elem(gpi) = Point2d (x(0), x(1));
vdir = Vec2d (dir(0), dir(1)); vdir = Vec2d (dir(0), dir(1));

159
libsrc/meshing/smoothing3.cpp
View File

@ -149,15 +149,15 @@ namespace netgen
static double eps = 1e-6; static double eps = 1e-6;
hx = x; hx = x;
for (int i = 1; i <= 3; i++) for (int i = 0; i < 3; i++)
{ {
hx.Elem(i) = x.Get(i) + eps * h; hx(i) = x(i) + eps * h;
double fr = Func (hx); double fr = Func (hx);
hx.Elem(i) = x.Get(i) - eps * h; hx(i) = x(i) - eps * h;
double fl = Func (hx); double fl = Func (hx);
hx.Elem(i) = x.Get(i); hx(i) = x(i);
g.Elem(i) = (fr - fl) / (2 * eps * h); g(i) = (fr - fl) / (2 * eps * h);
} }
return Func(x); return Func(x);
@ -246,17 +246,17 @@ namespace netgen
static Vector res; static Vector res;
res.SetSize (m.Height()); res.SetSize (m.Height());
for (i = 1;i <= 3; i++) for (i = 0;i < 3; i++)
hv.Elem(i) = vp.Get(i); hv(i) = vp(i);
hv.Elem(4) = 1; hv(3) = 1;
m.Mult (hv, res); m.Mult (hv, res);
for (i = 1; i <= res.Size(); i++) for (i = 1; i <= res.Size(); i++)
{ {
if (res.Get(i) < 1e-10) if (res(i-1) < 1e-10)
badness += 1e24; badness += 1e24;
else else
badness += 1 / res.Get(i); badness += 1 / res(i-1);
} }
return badness; return badness;
@ -269,15 +269,15 @@ namespace netgen
static double eps = 1e-6; static double eps = 1e-6;
hx = x; hx = x;
for (int i = 1; i <= 3; i++) for (int i = 0; i < 3; i++)
{ {
hx.Elem(i) = x.Get(i) + eps * h; hx(i) = x(i) + eps * h;
double fr = Func (hx); double fr = Func (hx);
hx.Elem(i) = x.Get(i) - eps * h; hx(i) = x(i) - eps * h;
double fl = Func (hx); double fl = Func (hx);
hx.Elem(i) = x.Get(i); hx(i) = x(i);
g.Elem(i) = (fr - fl) / (2 * eps * h); g(i) = (fr - fl) / (2 * eps * h);
} }
return Func(x); return Func(x);
@ -374,42 +374,12 @@ namespace netgen
double PointFunction :: PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const double PointFunction :: PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const
{ {
double f;//, delta = h * 1e-6; double f = 0;
// Point<3> hpp; // f = PointFunctionValue (pp);
f = PointFunctionValue (pp);
/*
hpp = pp;
hpp.X() = pp.X() + delta;
fr = PointFunctionValue (hpp);
hpp.X() = pp.X() - delta;
fl = PointFunctionValue (hpp);
grad.Elem(1) = (fr - fl) / (2 * delta);
hpp = pp;
hpp.Y() = pp.Y() + delta;
fr = PointFunctionValue (hpp);
hpp.Y() = pp.Y() - delta;
fl = PointFunctionValue (hpp);
grad.Elem(2) = (fr - fl) / (2 * delta);
hpp = pp;
hpp.Z() = pp.Z() + delta;
fr = PointFunctionValue (hpp);
hpp.Z() = pp.Z() - delta;
fl = PointFunctionValue (hpp);
grad.Elem(3) = (fr - fl) / (2 * delta);
*/
// new gradient calculation
// double badness;
Point<3> hp; Point<3> hp;
Vec<3> vgradi, vgrad(0,0,0); Vec<3> vgradi, vgrad(0,0,0);
// badness = 0;
hp = points[actpind]; hp = points[actpind];
points[actpind] = Point<3> (pp); points[actpind] = Point<3> (pp);
@ -421,7 +391,7 @@ namespace netgen
for (int k = 1; k <= 4; k++) for (int k = 1; k <= 4; k++)
if (el.PNum(k) == actpind) if (el.PNum(k) == actpind)
{ {
CalcTetBadnessGrad (points[el.PNum(1)], f += CalcTetBadnessGrad (points[el.PNum(1)],
points[el.PNum(2)], points[el.PNum(2)],
points[el.PNum(3)], points[el.PNum(3)],
points[el.PNum(4)], -1, k, vgradi); points[el.PNum(4)], -1, k, vgradi);
@ -606,19 +576,19 @@ namespace netgen
static Vector di; static Vector di;
int n = m.Height(); int n = m.Height();
p4.Elem(1) = pp(0); p4(0) = pp(0);
p4.Elem(2) = pp(1); p4(1) = pp(1);
p4.Elem(3) = pp(2); p4(2) = pp(2);
p4.Elem(4) = 1; p4(3) = 1;
di.SetSize (n); di.SetSize (n);
m.Mult (p4, di); m.Mult (p4, di);
double sum = 0; double sum = 0;
for (int i = 1; i <= n; i++) for (int i = 0; i < n; i++)
{ {
if (di.Get(i) > 0) if (di(i) > 0)
sum += 1 / di.Get(i); sum += 1 / di(i);
else else
return 1e16; return 1e16;
} }
@ -635,25 +605,25 @@ namespace netgen
int n = m.Height(); int n = m.Height();
p4.Elem(1) = pp(0); p4(0) = pp(0);
p4.Elem(2) = pp(1); p4(1) = pp(1);
p4.Elem(3) = pp(2); p4(2) = pp(2);
p4.Elem(4) = 1; p4(3) = 1;
di.SetSize (n); di.SetSize (n);
m.Mult (p4, di); m.Mult (p4, di);
double sum = 0; double sum = 0;
grad = 0; grad = 0;
for (int i = 1; i <= n; i++) for (int i = 0; i < n; i++)
{ {
if (di.Get(i) > 0) if (di(i) > 0)
{ {
double idi = 1 / di.Get(i); double idi = 1 / di(i);
sum += idi; sum += idi;
grad(0) -= idi * idi * m.Get(i, 1); grad(0) -= idi * idi * m(i, 0);
grad(1) -= idi * idi * m.Get(i, 2); grad(1) -= idi * idi * m(i, 1);
grad(2) -= idi * idi * m.Get(i, 3); grad(2) -= idi * idi * m(i, 2);
} }
else else
{ {
@ -774,9 +744,9 @@ namespace netgen
} }
hx = x; hx = x;
hx.Elem(i) = x.Get(i) + eps; hx(i-1) = x(i-1) + eps;
f11 = Func(hx); f11 = Func(hx);
hx.Elem(i) = x.Get(i) - eps; hx(i-1) = x(i-1) - eps;
f22 = Func(hx); f22 = Func(hx);
hesse.Elem(i, i) = (f11 + f22 - 2 * f) / (eps * eps) + 1e-12; hesse.Elem(i, i) = (f11 + f22 - 2 * f) / (eps * eps) + 1e-12;
@ -978,17 +948,6 @@ double Opti3EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
double CalcBad (const Mesh::T_POINTS & points, const Element & elem,
double h)
{
if (elem.GetType() == TET)
return CalcTetBadness (points[elem.PNum(1)],
points[elem.PNum(2)],
points[elem.PNum(3)],
points[elem.PNum(4)], h);
return 0;
}
extern double teterrpow; extern double teterrpow;
double CalcTotalBad (const Mesh::T_POINTS & points, double CalcTotalBad (const Mesh::T_POINTS & points,
@ -1100,10 +1059,10 @@ double JacobianPointFunction :: Func (const Vector & v) const
Point<3> hp = points.Elem(actpind); Point<3> hp = points.Elem(actpind);
points.Elem(actpind) = hp + Vec<3> (v.Get(1), v.Get(2), v.Get(3)); points.Elem(actpind) = hp + Vec<3> (v(0), v(1), v(2));
if(onplane) if(onplane)
points.Elem(actpind) -= (v.Get(1)*nv(0)+v.Get(2)*nv(1)+v.Get(3)*nv(2)) * nv; points.Elem(actpind) -= (v(0)*nv(0)+v(1)*nv(1)+v(2)*nv(2)) * nv;
for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++) for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++)
@ -1129,10 +1088,10 @@ FuncGrad (const Vector & x, Vector & g) const
double badness = 0;//, hbad; double badness = 0;//, hbad;
Point<3> hp = points.Elem(actpind); Point<3> hp = points.Elem(actpind);
points.Elem(actpind) = hp + Vec<3> (x.Get(1), x.Get(2), x.Get(3)); points.Elem(actpind) = hp + Vec<3> (x(0), x(1), x(2));
if(onplane) if(onplane)
points.Elem(actpind) -= (x.Get(1)*nv(0)+x.Get(2)*nv(1)+x.Get(3)*nv(2)) * nv; points.Elem(actpind) -= (x(0)*nv(0)+x(1)*nv(1)+x(2)*nv(2)) * nv;
Vec<3> hderiv; Vec<3> hderiv;
//Vec3d vdir; //Vec3d vdir;
@ -1154,7 +1113,7 @@ FuncGrad (const Vector & x, Vector & g) const
CalcJacobianBadnessGradient (points, lpi, hderiv); CalcJacobianBadnessGradient (points, lpi, hderiv);
for(k=0; k<3; k++) for(k=0; k<3; k++)
g.Elem(k+1) += hderiv(k); g(k) += hderiv(k);
/* /*
for (k = 1; k <= 3; k++) for (k = 1; k <= 3; k++)
@ -1174,10 +1133,10 @@ FuncGrad (const Vector & x, Vector & g) const
if(onplane) if(onplane)
{ {
double scal = nv(0)*g.Get(1) + nv(1)*g.Get(2) + nv(2)*g.Get(3); double scal = nv(0)*g(0) + nv(1)*g(1) + nv(2)*g(2);
g.Elem(1) -= scal*nv(0); g(0) -= scal*nv(0);
g.Elem(2) -= scal*nv(1); g(1) -= scal*nv(1);
g.Elem(3) -= scal*nv(2); g(2) -= scal*nv(2);
} }
//(*testout) << "g = " << g << endl; //(*testout) << "g = " << g << endl;
@ -1197,14 +1156,14 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
double badness = 0; double badness = 0;
Point<3> hp = points.Elem(actpind); Point<3> hp = points.Elem(actpind);
points.Elem(actpind) = Point<3> (hp + Vec3d (x.Get(1), x.Get(2), x.Get(3))); points.Elem(actpind) = Point<3> (hp + Vec3d (x(0), x(1), x(2)));
if(onplane) if(onplane)
points.Elem(actpind) -= (Vec3d (x.Get(1), x.Get(2), x.Get(3))*nv) * nv; points.Elem(actpind) -= (Vec3d (x(0), x(1), x(2))*nv) * nv;
double hderiv; double hderiv;
deriv = 0; deriv = 0;
Vec<3> vdir(dir.Get(1), dir.Get(2), dir.Get(3)); Vec<3> vdir(dir(0), dir(1), dir(2));
if(onplane) if(onplane)
{ {
@ -1569,9 +1528,9 @@ void Mesh :: ImproveMesh (OPTIMIZEGOAL goal)
//*testout << "start BFGS, pok" << endl; //*testout << "start BFGS, pok" << endl;
BFGS (x, freeminf, par); BFGS (x, freeminf, par);
//*testout << "BFGS complete, pok" << endl; //*testout << "BFGS complete, pok" << endl;
points[i](0) += x.Get(1); points[i](0) += x(0);
points[i](1) += x.Get(2); points[i](1) += x(1);
points[i](2) += x.Get(3); points[i](2) += x(2);
} }
} }
PrintDot ('\n'); PrintDot ('\n');
@ -1694,9 +1653,9 @@ void Mesh :: ImproveMeshJacobian (OPTIMIZEGOAL goal, const BitArray * usepoint)
//*testout << "start BFGS, Jacobian" << endl; //*testout << "start BFGS, Jacobian" << endl;
BFGS (x, pf, par); BFGS (x, pf, par);
//*testout << "end BFGS, Jacobian" << endl; //*testout << "end BFGS, Jacobian" << endl;
points.Elem(i)(0) += x.Get(1); points.Elem(i)(0) += x(0);
points.Elem(i)(1) += x.Get(2); points.Elem(i)(1) += x(1);
points.Elem(i)(2) += x.Get(3); points.Elem(i)(2) += x(2);
} }
else else
{ {
@ -1873,12 +1832,12 @@ void Mesh :: ImproveMeshJacobianOnSurface (const BitArray & usepoint,
BFGS (x, pf_sum, par); BFGS (x, pf_sum, par);
for(j=1; j<=3; j++) for(j=0; j<3; j++)
points.Elem(i)(j-1) += x.Get(j);// - scal*nv[i-1].X(j); points.Elem(i)(j) += x(j);// - scal*nv[i-1].X(j);
if(brother != -1) if(brother != -1)
for(j=1; j<=3; j++) for(j=0; j<3; j++)
points.Elem(brother)(j-1) += x.Get(j);// - scal*nv[brother-1].X(j); points.Elem(brother)(j) += x(j);// - scal*nv[brother-1].X(j);
} }

14
libsrc/stlgeom/stlgeom.cpp
View File

@ -183,7 +183,7 @@ void STLGeometry :: SmoothNormals()
for (k = 0; k < 3; k++) for (k = 0; k < 3; k++)
hv.Elem(k+1) = ngeom(k); hv(k) = ngeom(k);
hm.Mult (hv, hv2); hm.Mult (hv, hv2);
/* /*
@ -230,7 +230,7 @@ void STLGeometry :: SmoothNormals()
} }
for (k = 0; k < 3; k++) for (k = 0; k < 3; k++)
hv.Elem(k+1) = nnb(k); hv(k) = nnb(k);
hm.Mult (hv, hv2); hm.Mult (hv, hv2);
/* /*
@ -245,7 +245,7 @@ void STLGeometry :: SmoothNormals()
} }
m.Solve (rhs, sol); m.Solve (rhs, sol);
Vec3d newn(sol.Get(1), sol.Get(2), sol.Get(3)); Vec3d newn(sol(0), sol(1), sol(2));
newn /= (newn.Length() + 1e-24); newn /= (newn.Length() + 1e-24);
GetTriangle(i).SetNormal(newn); GetTriangle(i).SetNormal(newn);
@ -1562,9 +1562,9 @@ void STLGeometry :: NeighbourAnglesOfSelectedTrig()
for (i = 1; i <= NONeighbourTrigs(st); i++) for (i = 1; i <= NONeighbourTrigs(st); i++)
{ {
PrintMessage(1," triangle ", NeighbourTrig(st,i), ": angle = ", PrintMessage(1," triangle ", NeighbourTrig(st,i), ": angle = ",
180./M_PI*GetAngle(st, NeighbourTrig(st,i)), "°", 180./M_PI*GetAngle(st, NeighbourTrig(st,i)), "°",
", calculated = ", 180./M_PI*Angle(GetTriangle(st).GeomNormal(points), ", calculated = ", 180./M_PI*Angle(GetTriangle(st).GeomNormal(points),
GetTriangle(NeighbourTrig(st,i)).GeomNormal(points)), "°"); GetTriangle(NeighbourTrig(st,i)).GeomNormal(points)), "°");
} }
} }
} }
@ -2601,7 +2601,7 @@ void STLGeometry :: AddFaceEdges()
{ {
PrintFnStart("Add starting edges for faces"); PrintFnStart("Add starting edges for faces");
//für Kugel eine STLLine hinzufügen (Vorteil: verfeinerbar, unabhängig von Auflösung der Geometrie!!!): //für Kugel eine STLLine hinzufügen (Vorteil: verfeinerbar, unabhängig von Auflösung der Geometrie!!!):
//Grenze von 1. gefundener chart //Grenze von 1. gefundener chart
Array<int> edgecnt; Array<int> edgecnt;
@ -2676,7 +2676,7 @@ void STLGeometry :: LinkEdges()
//worked edges //worked edges
Array<int> we(GetNE()); Array<int> we(GetNE());
//setlineendpoints; wenn 180°, dann keine endpunkte //setlineendpoints; wenn 180°, dann keine endpunkte
//nur punkte mit 2 edges kommen in frage, da bei mehr oder weniger punkten ohnehin ein meshpoint hinkommt //nur punkte mit 2 edges kommen in frage, da bei mehr oder weniger punkten ohnehin ein meshpoint hinkommt
Vec3d v1,v2; Vec3d v1,v2;

13
libsrc/visualization/vsmesh.cpp
View File

@ -3,7 +3,6 @@
#include <mystdlib.h> #include <mystdlib.h>
#include <myadt.hpp> #include <myadt.hpp>
#include <meshing.hpp> #include <meshing.hpp>
#include <csg.hpp> #include <csg.hpp>
@ -983,9 +982,9 @@ namespace netgen
for (int i = 1; i <= locms.Size(); i++) for (int i = 1; i <= locms.Size(); i++)
{ {
Point3d p = mesh->Point(i); Point3d p = mesh->Point(i);
locms.Elem(i) = mesh->GetH (p); locms(i-1) = mesh->GetH (p);
if (locms.Elem(i) > maxh) maxh = locms.Elem(i); if (locms(i-1) > maxh) maxh = locms(i-1);
if (locms.Elem(i) < minh) minh = locms.Elem(i); if (locms(i-1) < minh) minh = locms(i-1);
} }
if (!locms.Size()) if (!locms.Size())
{ minh = 1; maxh = 10; } { minh = 1; maxh = 10; }
@ -1145,11 +1144,11 @@ namespace netgen
if (vispar.colormeshsize) if (vispar.colormeshsize)
{ {
SetOpenGlColor (locms.Get(el[0]), minh, maxh, 1); SetOpenGlColor (locms(el[0]-1), minh, maxh, 1);
glVertex3dv (lp0); glVertex3dv (lp0);
SetOpenGlColor (locms.Get(el[1]), minh, maxh, 1); SetOpenGlColor (locms(el[1]-1), minh, maxh, 1);
glVertex3dv (lp1); glVertex3dv (lp1);
SetOpenGlColor (locms.Get(el[2]), minh, maxh, 1); SetOpenGlColor (locms(el[2]-1), minh, maxh, 1);
glVertex3dv (lp2); glVertex3dv (lp2);
} }
else else

10
ng/nginterface.cpp
View File

@ -149,10 +149,8 @@ using namespace netgen;
void Ng_LoadGeometry (const char * filename) void Ng_LoadGeometry (const char * filename)
{ {
// if (printmessage_importance>0)
// cout << "CALLED NG LOAD GEOMETRY" << endl;
if (printmessage_importance>0)
cout << "CALLED NG LOAD GEOMETRY" << endl;
geometry.Reset (new CSGeometry ()); geometry.Reset (new CSGeometry ());
geometry2d.Reset (); geometry2d.Reset ();
@ -171,6 +169,8 @@ void Ng_LoadGeometry (const char * filename)
if (strcmp(filename,"")==0) if (strcmp(filename,"")==0)
return; return;
PrintMessage (1, "Load geometry from file ", filename);
ifstream infile (filename); ifstream infile (filename);
if ((strcmp (&filename[strlen(filename)-3], "geo") == 0) || if ((strcmp (&filename[strlen(filename)-3], "geo") == 0) ||
@ -1936,8 +1936,6 @@ void Ng_SetSolutionData (Ng_SolutionData * soldata)
// vssolution.ClearSolutionData (); // vssolution.ClearSolutionData ();
VisualSceneSolution::SolData * vss = new VisualSceneSolution::SolData; VisualSceneSolution::SolData * vss = new VisualSceneSolution::SolData;
cout << "Add solution " << soldata->name << ", type = " << soldata->soltype << endl;
vss->name = new char[strlen (soldata->name)+1]; vss->name = new char[strlen (soldata->name)+1];
strcpy (vss->name, soldata->name); strcpy (vss->name, soldata->name);

2
nglib/Makefile.am
View File

@ -1,5 +1,7 @@
include_HEADERS = nglib.h include_HEADERS = nglib.h
dist_pkgdata_DATA = cube.surf
AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include $(OCCFLAGS) $(MPI_INCLUDES) AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include $(OCCFLAGS) $(MPI_INCLUDES)
lib_LTLIBRARIES = libnglib.la lib_LTLIBRARIES = libnglib.la

22
nglib/cube.surf Normal file
View File

@ -0,0 +1,22 @@
8
0 0 0
1 0 0
1 1 1
1 0 1
0 1 1
0 0 1
0 1 0
1 1 0
12
2 1 7
8 2 7
6 1 2
4 6 2
4 3 5
5 6 4
8 3 4
8 4 2
5 3 8
7 5 8
1 6 5
7 1 5