diff --git a/configure.ac b/configure.ac index 94dca563..2caa0bee 100644 --- a/configure.ac +++ b/configure.ac @@ -1,4 +1,4 @@ -AC_INIT([netgen],[4.9.9],[],[]) +AC_INIT([netgen],[4.9.10-dev],[],[]) AM_INIT_AUTOMAKE([-Wall -Werror foreign]) AC_CONFIG_MACRO_DIR([m4]) diff --git a/libsrc/geom2d/spline.hpp b/libsrc/geom2d/spline.hpp index 5ac5018f..15c6f18d 100644 --- a/libsrc/geom2d/spline.hpp +++ b/libsrc/geom2d/spline.hpp @@ -587,8 +587,6 @@ Vec SplineSeg3 :: GetTangent (const double t) const template void SplineSeg3 :: GetCoeff (Vector & u) const { - double t; - int i; Point p; DenseMatrix a(6, 6); DenseMatrix ata(6, 6); @@ -598,23 +596,23 @@ void SplineSeg3 :: GetCoeff (Vector & u) const // ata.SetSymmetric(1); - t = 0; - for (i = 1; i <= 5; i++, t += 0.25) + double t = 0; + for (int i = 0; i < 5; i++, t += 0.25) { p = GetPoint (t); - a.Elem(i, 1) = p(0) * p(0); - a.Elem(i, 2) = p(1) * p(1); - a.Elem(i, 3) = p(0) * p(1); - a.Elem(i, 4) = p(0); - a.Elem(i, 5) = p(1); - a.Elem(i, 6) = 1; + a(i, 0) = p(0) * p(0); + a(i, 1) = p(1) * p(1); + a(i, 2) = p(0) * p(1); + a(i, 3) = p(0); + a(i, 4) = p(1); + a(i, 5) = 1; } - a.Elem(6, 1) = 1; + a(5, 0) = 1; CalcAtA (a, ata); u = 0; - u.Elem(6) = 1; + u(5) = 1; a.MultTrans (u, f); ata.Solve (f, u); } diff --git a/libsrc/gprim/geomtest3d.cpp b/libsrc/gprim/geomtest3d.cpp index 97501d7e..bb4000b3 100644 --- a/libsrc/gprim/geomtest3d.cpp +++ b/libsrc/gprim/geomtest3d.cpp @@ -23,12 +23,12 @@ IntersectTriangleLine (const Point<3> ** tri, const Point<3> ** line) << "tri = " << *tri[0] << ", " << *tri[1] << ", " << *tri[2] << 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.Elem(i, 2) = vt1.X(i); - a.Elem(i, 3) = vt2.X(i); - rs.Elem(i) = vrs.X(i); + a(i, 0) = -vl.X(i+1); + a(i, 1) = vt1.X(i+1); + a(i, 2) = vt2.X(i+1); + rs(i) = vrs.X(i+1); } double det = a.Det(); @@ -70,12 +70,12 @@ IntersectTriangleLine (const Point<3> ** tri, const Point<3> ** line) double eps = 1e-6; if ( - (lami.Get(1) >= -eps && lami.Get(1) <= 1+eps && - lami.Get(2) >= -eps && lami.Get(3) >= -eps && - lami.Get(2) + lami.Get(3) <= 1+eps) && ! - (lami.Get(1) >= eps && lami.Get(1) <= 1-eps && - lami.Get(2) >= eps && lami.Get(3) >= eps && - lami.Get(2) + lami.Get(3) <= 1-eps) ) + (lami(0) >= -eps && lami(0) <= 1+eps && + lami(1) >= -eps && lami(2) >= -eps && + lami(1) + lami(2) <= 1+eps) && ! + (lami(0) >= eps && lami(0) <= 1-eps && + lami(1) >= eps && lami(2) >= 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 && - lami.Get(2) >= 0 && lami.Get(3) >= 0 && lami.Get(2) + lami.Get(3) <= 1) + if (lami(0) >= 0 && lami(0) <= 1 && + lami(1) >= 0 && lami(2) >= 0 && lami(1) + lami(2) <= 1) { return 1; @@ -860,12 +860,12 @@ int CalcTriangleCenter (const Point3d ** pts, Point3d & c) Vec3d v1(*pts[0], *pts[1]); Vec3d v2(*pts[0], *pts[2]); - rs.Elem(1) = v1 * v1; - rs.Elem(2) = v2 * v2; + rs(0) = v1 * v1; + rs(1) = v2 * v2; - a.Elem(1,1) = 2 * rs.Get(1); - a.Elem(1,2) = a.Elem(2,1) = 2 * (v1 * v2); - a.Elem(2,2) = 2 * rs.Get(2); + a(0,0) = 2 * rs(0); + a(0,1) = a(1,0) = 2 * (v1 * v2); + a(1,1) = 2 * rs(1); if (fabs (a.Det()) <= 1e-12 * h * h) { @@ -877,8 +877,8 @@ int CalcTriangleCenter (const Point3d ** pts, Point3d & c) inva.Mult (rs, sol); c = *pts[0]; - v1 *= sol.Get(1); - v2 *= sol.Get(2); + v1 *= sol(0); + v2 *= sol(1); c += v1; c += v2; diff --git a/libsrc/gprim/transform3d.cpp b/libsrc/gprim/transform3d.cpp index ea62fffe..c3a6efe8 100644 --- a/libsrc/gprim/transform3d.cpp +++ b/libsrc/gprim/transform3d.cpp @@ -90,21 +90,21 @@ void Transformation3d :: CalcInverse (Transformation3d & inv) const static Vector b(3), sol(3); int i, j; - for (i = 1; i <= 3; i++) + for (int i = 0; i < 3; i++) { - b.Elem(i) = offset[i-1]; - for (j = 1; j <= 3; j++) - a.Elem(i, j) = lin[i-1][j-1]; + b(i) = offset[i]; + for (j = 0; j < 3; j++) + a(i, j) = lin[i][j]; } ::netgen::CalcInverse (a, inva); inva.Mult (b, sol); - for (i = 1; i <= 3; i++) + for (int i = 0; i < 3; i++) { - inv.offset[i-1] = -sol.Get(i); - for (j = 1; j <= 3; j++) - inv.lin[i-1][j-1] = inva.Elem(i, j); + inv.offset[i] = -sol(i); + for (j = 0; j < 3; j++) + inv.lin[i][j] = inva(i, j); } } diff --git a/libsrc/linalg/Makefile.am b/libsrc/linalg/Makefile.am index d39fa91b..9503b31e 100644 --- a/libsrc/linalg/Makefile.am +++ b/libsrc/linalg/Makefile.am @@ -2,6 +2,8 @@ noinst_HEADERS = densemat.hpp linalg.hpp polynomial.hpp vector.hpp opti.hpp AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include METASOURCES = AUTO 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 diff --git a/libsrc/linalg/bfgs.cpp b/libsrc/linalg/bfgs.cpp index d106b6f1..e0f40d68 100644 --- a/libsrc/linalg/bfgs.cpp +++ b/libsrc/linalg/bfgs.cpp @@ -1,9 +1,9 @@ /***************************************************************************/ /* */ /* 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 */ /* */ /***************************************************************************/ @@ -25,37 +25,36 @@ void Cholesky (const DenseMatrix & a, double x; - int i, j, k; int n = a.Height(); // (*testout) << "a = " << a << endl; 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); - for (k = 1; k < i; k++) - x -= l.Get(i, k) * l.Get(j, k) * d.Get(k); - + for (int k = 1; k < i; k++) + x -= l.Get(i, k) * l.Get(j, k) * d(k-1); + if (i == j) { - d.Elem(i) = x; + d(i-1) = x; } 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; - for (j = i+1; j <= n; j++) + for (int j = i+1; j <= n; j++) 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 // 1 .. sonst - int i, j, n; - - n = l.Height(); + int n = l.Height(); Vector v(n); double t, told, xi; @@ -172,9 +169,9 @@ int LDLtUpdate (DenseMatrix & l, Vector & d, double a, const Vector & u) told = 1; 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) { @@ -182,14 +179,14 @@ int LDLtUpdate (DenseMatrix & l, Vector & d, double a, const Vector & u) 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); - l.Elem(i, j) += xi * v.Elem(i); + v(i-1) -= v(j-1) * l.Elem(i, j); + l.Elem(i, j) += xi * v(i-1); } told = t; @@ -209,7 +206,7 @@ double BFGS ( { - int i, j, n = x.Size(); + int n = x.Size(); long it; char a1crit, a3acrit; @@ -233,14 +230,14 @@ double BFGS ( int ifail; // o: 0 .. Erfolg // -1 .. Unterschreitung von fmin // 1 .. kein Erfolg bei Liniensuche - // 2 .. Überschreitung von itmax + // 2 .. Ãœberschreitung von itmax typx = par.typx; typf = par.typf; l = 0; - for (i = 1; i <= n; i++) + for (int i = 1; i <= n; i++) l.Elem(i, i) = 1; f = fun.FuncGrad (x, g); @@ -255,10 +252,10 @@ double BFGS ( if (it % (5 * n) == 0) { - for (i = 1; i <= n; i++) - d.Elem(i) = typf/ sqr (typx.Get(i)); // 1; - for (i = 2; i <= n; i++) - for (j = 1; j < i; j++) + for (int i = 1; i <= n; i++) + d(i-1) = typf/ sqr (typx(i-1)); // 1; + for (int i = 2; i <= n; i++) + for (int j = 1; j < i; j++) l.Elem(i, j) = 0; /* @@ -357,8 +354,8 @@ double BFGS ( hd = eps * max2 (typf, fabs (f)); a1crit = 1; - for (i = 1; i <= n; i++) - if ( fabs (g.Elem(i)) * max2 (typx.Elem(i), fabs (x.Elem(i))) > hd) + for (int i = 1; i <= n; i++) + if ( fabs (g(i-1)) * max2 (typx(i-1), fabs (x(i-1))) > hd) a1crit = 0; diff --git a/libsrc/linalg/densemat.cpp b/libsrc/linalg/densemat.cpp index 56b6c158..768c7d5c 100644 --- a/libsrc/linalg/densemat.cpp +++ b/libsrc/linalg/densemat.cpp @@ -27,21 +27,21 @@ namespace netgen } /* - DenseMatrix :: DenseMatrix (int h, int w, const double * d) + DenseMatrix :: DenseMatrix (int h, int w, const double * d) : BaseMatrix (h, w) - { - int size = h * w; - int i; - - if (size) { - data = new double[size]; - for (i = 0; i < size; i++) - data[i] = d[i]; + int size = h * w; + int i; + + if (size) + { + data = new double[size]; + for (i = 0; i < size; i++) + data[i] = d[i]; } - else + else data = NULL; - } + } */ DenseMatrix :: DenseMatrix (const DenseMatrix & m2) @@ -76,21 +76,21 @@ namespace netgen /* -DenseMatrix & DenseMatrix :: operator= (const BaseMatrix & m2) - { - int i, j; + DenseMatrix & DenseMatrix :: operator= (const BaseMatrix & m2) + { + int i, j; - SetSize (m2.Height(), m2.Width()); + SetSize (m2.Height(), m2.Width()); - if (data) + if (data) for (i = 1; i <= Height(); i++) - for (j = 1; j <= Width(); j++) - Set (i, j, m2(i, j)); - else + for (j = 1; j <= Width(); j++) + Set (i, j, m2(i, j)); + else (*myerr) << "DenseMatrix::Operator=: Matrix not allocated" << endl; - return *this; - } + return *this; + } */ @@ -109,610 +109,661 @@ DenseMatrix & DenseMatrix :: operator= (const BaseMatrix & m2) double * p, * q; if (Height() != m2.Height() || Width() != m2.Width()) - { - (*myerr) << "DenseMatrix::Operator+=: Sizes don't fit" << endl; - return *this; - } + { + (*myerr) << "DenseMatrix::Operator+=: Sizes don't fit" << endl; + return *this; + } if (data) { p = data; q = m2.data; for (i = Width() * Height(); i > 0; i--) - { - *p += *q; - p++; - q++; + { + *p += *q; + p++; + q++; + } } - } - else - (*myerr) << "DenseMatrix::Operator+=: Matrix not allocated" << endl; + else + (*myerr) << "DenseMatrix::Operator+=: Matrix not allocated" << endl; - return *this; + return *this; } -DenseMatrix & DenseMatrix :: operator-= (const DenseMatrix & m2) + DenseMatrix & DenseMatrix :: operator-= (const DenseMatrix & m2) { - int i; - double * p, * q; + int i; + double * p, * q; - if (Height() != m2.Height() || Width() != m2.Width()) - { - (*myerr) << "DenseMatrix::Operator-=: Sizes don't fit" << endl; - return *this; - } - - if (data) - { - p = data; - q = m2.data; - for (i = Width() * Height(); i > 0; i--) + if (Height() != m2.Height() || Width() != m2.Width()) { - *p -= *q; - p++; - q++; + (*myerr) << "DenseMatrix::Operator-=: Sizes don't fit" << endl; + return *this; } - } - else - (*myerr) << "DenseMatrix::Operator-=: Matrix not allocated" << endl; - return *this; + if (data) + { + p = data; + q = m2.data; + for (i = Width() * Height(); i > 0; i--) + { + *p -= *q; + p++; + q++; + } + } + else + (*myerr) << "DenseMatrix::Operator-=: Matrix not allocated" << endl; + + return *this; } /* -double & DenseMatrix :: operator() (int i, int j) -{ - if (i >= 1 && j >= 1 && i <= height && j <= width) - return Elem(i,j); - else (*myerr) << "DenseMatrix: index (" << i << "," << j << ") out of range (1.." - << height << ",1.." << width << ")\n"; - static double dummy = 0; - return dummy; -} - - double DenseMatrix :: operator() (int i, int j) const - { + double & DenseMatrix :: operator() (int i, int j) + { if (i >= 1 && j >= 1 && i <= height && j <= width) - return Get(i,j); + return Elem(i,j); else (*myerr) << "DenseMatrix: index (" << i << "," << j << ") out of range (1.." - << height << ",1.." << width << ")\n"; + << height << ",1.." << width << ")\n"; + static double dummy = 0; + return dummy; + } + + double DenseMatrix :: operator() (int i, int j) const + { + if (i >= 1 && j >= 1 && i <= height && j <= width) + return Get(i,j); + else (*myerr) << "DenseMatrix: index (" << i << "," << j << ") out of range (1.." + << height << ",1.." << width << ")\n"; static double dummy = 0; return dummy; - } + } */ -DenseMatrix & DenseMatrix :: operator= (double v) + DenseMatrix & DenseMatrix :: operator= (double v) { - int i; - double * p = data; + int i; + double * p = data; - if (data) - for (i = width*height; i > 0; i--, p++) - *p = v; + if (data) + for (i = width*height; i > 0; i--, p++) + *p = v; - return *this; + return *this; } -DenseMatrix & DenseMatrix :: operator*= (double v) + DenseMatrix & DenseMatrix :: operator*= (double v) { - int i; - double * p = data; + int i; + double * p = data; - if (data) - for (i = width*height; i > 0; i--, p++) - *p *= v; + if (data) + for (i = width*height; i > 0; i--, p++) + *p *= v; - return *this; + return *this; } -double DenseMatrix :: Det () const + double DenseMatrix :: Det () const { - if (width != height) - { - (*myerr) << "DenseMatrix :: Det: width != height" << endl; - return 0; - } - - switch (width) - { - case 1: return Get(1, 1); - case 2: return Get(1) * Get(4) - Get(2) * Get(3); - - case 3: return Get(1) * Get(5) * Get(9) - + Get(2) * Get(6) * Get(7) - + Get(3) * Get(4) * Get(8) - - Get(1) * Get(6) * Get(8) - - Get(2) * Get(4) * Get(9) - - Get(3) * Get(5) * Get(7); - default: + if (width != height) { - (*myerr) << "Matrix :: Det: general size not implemented (size=" << width << ")" << endl; - return 0; + (*myerr) << "DenseMatrix :: Det: width != height" << endl; + return 0; + } + + switch (width) + { + case 1: return data[0]; + case 2: return data[0] * data[3] - data[1] * data[2]; + case 3: return data[0] * data[4] * data[8] + + data[1] * data[5] * data[6] + + data[2] * data[3] * data[7] + - data[0] * data[5] * data[7] + - data[1] * data[3] * data[8] + - data[2] * data[4] * data[6]; + default: + { + (*myerr) << "Matrix :: Det: general size not implemented (size=" << width << ")" << endl; + return 0; + } } - } } -void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2) + void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2) { - // int i, j, k, n; - double det; - // DenseMatrix m1 = hm1; + double det; - if (m1.width != m1.height) - { - (*myerr) << "CalcInverse: matrix not symmetric" << endl; - return; - } - if (m1.width != m2.width || m1.height != m2.height) - { - (*myerr) << "CalcInverse: dim(m2) != dim(m1)" << endl; - return; - } - - - if (m1.Width() <= 3) - { - det = m1.Det(); - if (det == 0) + if (m1.width != m1.height) { - (*myerr) << "CalcInverse: Matrix singular" << endl; - return; + (*myerr) << "CalcInverse: matrix not symmetric" << endl; + return; + } + if (m1.width != m2.width || m1.height != m2.height) + { + (*myerr) << "CalcInverse: dim(m2) != dim(m1)" << endl; + return; } - det = 1e0 / det; - switch (m1.width) + + if (m1.Width() <= 3) { - case 1: - { - m2.Set(1, 1, det); - return; - } - case 2: - { - m2.Set(1, 1, det * m1.Get(4)); - m2.Set(2, 2, det * m1.Get(1)); - m2.Set(1, 2, - det * m1.Get(2)); - m2.Set(2, 1, - det * m1.Get(3)); - return; - } - case 3: - { - m2.Set(1, 1, det * (m1.Get(5) * m1.Get(9) - m1.Get(6) * m1.Get(8))); - m2.Set(2, 1, -det * (m1.Get(4) * m1.Get(9) - m1.Get(6) * m1.Get(7))); - m2.Set(3, 1, det * (m1.Get(4) * m1.Get(8) - m1.Get(5) * m1.Get(7))); + det = m1.Det(); + if (det == 0) + { + (*myerr) << "CalcInverse: Matrix singular" << endl; + return; + } - m2.Set(1, 2, -det * (m1.Get(2) * m1.Get(9) - m1.Get(3) * m1.Get(8))); - m2.Set(2, 2, det * (m1.Get(1) * m1.Get(9) - m1.Get(3) * m1.Get(7))); - m2.Set(3, 2, -det * (m1.Get(1) * m1.Get(8) - m1.Get(2) * m1.Get(7))); + det = 1.0 / det; + switch (m1.width) + { + case 1: + { + m2(0,0) = det; + return; + } + case 2: + { + m2(0,0) = det * m1(3); + m2(1,1) = det * m1(0); + m2(0,1) = -det * m1(1); + m2(1,0) = - det * m1(2); + return; + } + case 3: + { + m2(0, 0) = det * (m1(4) * m1(8) - m1(5) * m1(7)); + m2(1, 0) = -det * (m1(3) * m1(8) - m1(5) * m1(6)); + m2(2, 0) = det * (m1(3) * m1(7) - m1(4) * m1(6)); - m2.Set(1, 3, det * (m1.Get(2) * m1.Get(6) - m1.Get(3) * m1.Get(5))); - m2.Set(2, 3, -det * (m1.Get(1) * m1.Get(6) - m1.Get(3) * m1.Get(4))); - m2.Set(3, 3, det * (m1.Get(1) * m1.Get(5) - m1.Get(2) * m1.Get(4))); - return; - } + m2(0, 1) = -det * (m1(1) * m1(8) - m1(2) * m1(7)); + m2(1, 1) = det * (m1(0) * m1(8) - m1(2) * m1(6)); + m2(2, 1) = -det * (m1(0) * m1(7) - m1(1) * m1(6)); + + m2(0, 2) = det * (m1(1) * m1(5) - m1(2) * m1(4)); + m2(1, 2) = -det * (m1(0) * m1(5) - m1(2) * m1(3)); + m2(2, 2) = det * (m1(0) * m1(4) - m1(1) * m1(3)); + return; + } + } } - } - else - { - int i, j, k, n; - n = m1.Height(); + else + { + int i, j, k, n; + n = m1.Height(); #ifdef CHOL - int dots = (n > 200); + int dots = (n > 200); - // Cholesky + // Cholesky - double x; - Vector p(n); + double x; + Vector p(n); - m2 = m1; - /* - m2.SetSymmetric(); - if (!m2.Symmetric()) - cerr << "m should be symmetric for Cholesky" << endl; - */ + m2 = m1; + /* + m2.SetSymmetric(); + if (!m2.Symmetric()) + cerr << "m should be symmetric for Cholesky" << endl; + */ - for (i = 1; i <= n; i++) - for (j = 1; j < i; j++) - m2.Elem(j, i) = m2.Get(i, j); + for (i = 1; i <= n; i++) + for (j = 1; j < i; j++) + m2.Elem(j, i) = m2.Get(i, j); - for (i = 1; i <= n; i++) - { - if (dots && i % 10 == 0) - (*mycout) << "." << flush; + for (i = 1; i <= n; i++) + { + if (dots && i % 10 == 0) + (*mycout) << "." << flush; - for (j = i; j <= n; j++) - { - x = m2.Get(i, j); + for (j = i; j <= n; j++) + { + x = m2.Get(i, j); - const double * pik = &m2.Get(i, 1); - const double * pjk = &m2.Get(j, 1); + const double * pik = &m2.Get(i, 1); + const double * pjk = &m2.Get(j, 1); - for (k = i-2; k >= 0; --k, ++pik, ++pjk) - x -= (*pik) * (*pjk); + for (k = i-2; k >= 0; --k, ++pik, ++pjk) + x -= (*pik) * (*pjk); - // for (k = i-1; k >= 1; --k) - // x -= m2.Get(j, k) * m2.Get(i, k); + // for (k = i-1; k >= 1; --k) + // x -= m2.Get(j, k) * m2.Get(i, k); - if (i == j) - { - if (x <= 0) - { - cerr << "Matrix indefinite 1" << endl; - return; - } + if (i == j) + { + if (x <= 0) + { + cerr << "Matrix indefinite 1" << endl; + return; + } - p.Elem(i) = 1 / sqrt(x); - } - else - { - m2.Elem(j, i) = x * p.Get(i); - } - } - } + p.Elem(i) = 1 / sqrt(x); + } + else + { + m2.Elem(j, i) = x * p.Get(i); + } + } + } - for (i = 1; i <= n; i++) - m2.Elem(i, i) = 1 / p.Get(i); + for (i = 1; i <= n; i++) + m2.Elem(i, i) = 1 / p.Get(i); - // check: A = L L^t + // check: A = L L^t -// for (i = 1; i <= n; i++) -// for (j = 1; j <= n; j++) -// { -// x = 0; -// for (k = 1; k <= i && k <= j; k++) -// x += m2.Get(i, k) * m2.Get(j, k); -// (*testout) << "err " << i << "," << j << " = " << (m1.Get(i, j) - x) << endl; -// } + // for (i = 1; i <= n; i++) + // for (j = 1; j <= n; j++) + // { + // x = 0; + // for (k = 1; k <= i && k <= j; k++) + // x += m2.Get(i, k) * m2.Get(j, k); + // (*testout) << "err " << i << "," << j << " = " << (m1.Get(i, j) - x) << endl; + // } - // calc L^{-1}, store upper triangle + // calc L^{-1}, store upper triangle - // DenseMatrix hm(n); - // hm = m2; + // DenseMatrix hm(n); + // hm = m2; - for (i = 1; i <= n; i++) - { - if (dots && i % 10 == 0) - (*mycout) << "+" << flush; + for (i = 1; i <= n; i++) + { + if (dots && i % 10 == 0) + (*mycout) << "+" << flush; - for (j = i; j <= n; j++) - { - x = 0; - if (j == i) x = 1; + for (j = i; j <= n; j++) + { + x = 0; + if (j == i) x = 1; - const double * pjk = &m2.Get(j, i); - const double * pik = &m2.Get(i, i); - for (k = i; k < j; k++, ++pjk, ++pik) - x -= *pik * *pjk; + const double * pjk = &m2.Get(j, i); + const double * pik = &m2.Get(i, i); + for (k = i; k < j; k++, ++pjk, ++pik) + x -= *pik * *pjk; - // for (k = i; k < j; k++) - // x -= m2.Get(j, k) * m2.Get(i, k); + // for (k = i; k < j; k++) + // x -= m2.Get(j, k) * m2.Get(i, k); - m2.Elem(i, j) = x / m2.Get(j, j); - } - } + m2.Elem(i, j) = x / m2.Get(j, j); + } + } -// (*testout) << "check L^-1" << endl; -// for (i = 1; i <= n; i++) -// for (j = 1; j <= n; j++) -// { -// x = 0; -// for (k = j; k <= i; k++) -// x += hm.Get(i, k) * m2.Get(j, k); -// (*testout) << "i, j = " << i << "," << j << " x = " << x << endl; -// } + // (*testout) << "check L^-1" << endl; + // for (i = 1; i <= n; i++) + // for (j = 1; j <= n; j++) + // { + // x = 0; + // for (k = j; k <= i; k++) + // x += hm.Get(i, k) * m2.Get(j, k); + // (*testout) << "i, j = " << i << "," << j << " x = " << x << endl; + // } - // calc A^-1 = L^-T * L^-1 + // calc A^-1 = L^-T * L^-1 - for (i = 1; i <= n; i++) - { - if (dots && i % 10 == 0) - (*mycout) << "-" << flush; + for (i = 1; i <= n; i++) + { + if (dots && i % 10 == 0) + (*mycout) << "-" << flush; - for (j = 1; j <= i; j++) - { - x = 0; - k = i; - if (j > i) k = j; + for (j = 1; j <= i; j++) + { + x = 0; + k = i; + if (j > i) k = j; - const double * pik = &m2.Get(i, k); - const double * pjk = &m2.Get(j, k); + const double * pik = &m2.Get(i, k); + const double * pjk = &m2.Get(j, k); - for ( ; k <= n; ++k, ++pik, ++pjk) - x += *pik * *pjk; - // for ( ; k <= n; k++) - // x += m2.Get(i, k) * m2.Get(j, k); + for ( ; k <= n; ++k, ++pik, ++pjk) + x += *pik * *pjk; + // for ( ; k <= n; k++) + // x += m2.Get(i, k) * m2.Get(j, k); - m2.Elem(i, j) = x; - } - } + m2.Elem(i, j) = x; + } + } - for (i = 1; i <= n; i++) - for (j = 1; j < i; j++) - m2.Elem(j, i) = m2.Get(i, j); + for (i = 1; i <= n; i++) + for (j = 1; j < i; j++) + m2.Elem(j, i) = m2.Get(i, j); - if (dots) (*mycout) << endl; + if (dots) (*mycout) << endl; #endif - // Gauss - Jordan - algorithm + // Gauss - Jordan - algorithm - int r, hi; - double max, hr; + int r, hi; + double max, hr; - Array p(n); // pivot-permutation - Vector hv(n); + Array p(n); // pivot-permutation + Vector hv(n); - m2 = m1; + m2 = m1; - /* - if (m2.Symmetric()) - for (i = 1; i <= n; i++) - for (j = 1; j < i; j++) - m2.Elem(j, i) = m2.Get(i, j); - */ + /* + if (m2.Symmetric()) + for (i = 1; i <= n; i++) + for (j = 1; j < i; j++) + m2.Elem(j, i) = m2.Get(i, j); + */ - // Algorithm of Stoer, Einf. i. d. Num. Math, S 145 + // Algorithm of Stoer, Einf. i. d. Num. Math, S 145 - for (j = 1; j <= n; j++) - p.Set(j, j); + for (j = 1; j <= n; j++) + p.Set(j, j); - for (j = 1; j <= n; j++) - { - // pivot search + for (j = 1; j <= n; j++) + { + // pivot search - max = fabs(m2.Get(j, j)); - r = j; + max = fabs(m2.Get(j, j)); + r = j; - for (i = j+1; i <= n ;i++) - if (fabs (m2.Get(i, j)) > max) - { - r = i; - max = fabs (m2.Get(i, j)); - } + for (i = j+1; i <= n ;i++) + if (fabs (m2.Get(i, j)) > max) + { + r = i; + max = fabs (m2.Get(i, j)); + } - if (max < 1e-20) - { - cerr << "Inverse matrix: matrix singular" << endl; - return; - } + if (max < 1e-20) + { + cerr << "Inverse matrix: matrix singular" << endl; + return; + } - r = j; + r = j; - // exchange rows - if (r > j) - { - for (k = 1; k <= n; k++) - { - hr = m2.Get(j, k); - m2.Elem(j, k) = m2.Get(r, k); - m2.Elem(r, k) = hr; - } - hi = p.Get(j); - p.Elem(j) = p.Get(r); - p.Elem(r) = hi; - } + // exchange rows + if (r > j) + { + for (k = 1; k <= n; k++) + { + hr = m2.Get(j, k); + m2.Elem(j, k) = m2.Get(r, k); + m2.Elem(r, k) = hr; + } + hi = p.Get(j); + p.Elem(j) = p.Get(r); + p.Elem(r) = hi; + } - // transformation + // transformation - hr = 1 / m2.Get(j, j); - for (i = 1; i <= n; i++) - m2.Elem(i, j) *= hr; - m2.Elem(j, j) = hr; + hr = 1 / m2.Get(j, j); + for (i = 1; i <= n; i++) + m2.Elem(i, j) *= hr; + m2.Elem(j, j) = hr; - for (k = 1; k <= n; k++) - if (k != j) - { - for (i = 1; i <= n; i++) - if (i != j) - m2.Elem(i, k) -= m2.Elem(i, j) * m2.Elem(j, k); - m2.Elem(j, k) *= -hr; - } - } + for (k = 1; k <= n; k++) + if (k != j) + { + for (i = 1; i <= n; i++) + if (i != j) + m2.Elem(i, k) -= m2.Elem(i, j) * m2.Elem(j, k); + m2.Elem(j, k) *= -hr; + } + } - // col exchange + // col exchange - for (i = 1; i <= n; i++) - { - for (k = 1; k <= n; k++) - hv.Elem(p.Get(k)) = m2.Get(i, k); - for (k = 1; k <= n; k++) - m2.Elem(i, k) = hv.Get(k); - } + for (i = 1; i <= n; i++) + { + for (k = 1; k <= n; k++) + hv(p.Get(k)-1) = m2.Get(i, k); + for (k = 1; k <= n; k++) + m2.Elem(i, k) = hv(k-1); + } - /* - if (m1.Symmetric()) - for (i = 1; i <= n; i++) - for (j = 1; j < i; j++) + /* + if (m1.Symmetric()) + for (i = 1; i <= n; i++) + for (j = 1; j < i; j++) m1.Elem(j, i) = m1.Get(i, j); - m2 = 0; + m2 = 0; - for (i = 1; i <= n; i++) - m2.Elem(i, i) = 1; + for (i = 1; i <= n; i++) + m2.Elem(i, i) = 1; - for (i = 1; i <= n; i++) - { - // (*mycout) << '.' << flush; - q = m1.Get(i, i); - for (k = 1; k <= n; k++) - { - m1.Elem(i, k) /= q; - m2.Elem(i, k) /= q; - } + for (i = 1; i <= n; i++) + { + // (*mycout) << '.' << flush; + q = m1.Get(i, i); + for (k = 1; k <= n; k++) + { + m1.Elem(i, k) /= q; + m2.Elem(i, k) /= q; + } - for (j = i+1; j <= n; j++) - { - q = m1.Elem(j, i); + for (j = i+1; j <= n; j++) + { + q = m1.Elem(j, i); - double * m1pi = &m1.Elem(i, i); - double * m1pj = &m1.Elem(j, i); + double * m1pi = &m1.Elem(i, i); + double * m1pj = &m1.Elem(j, i); - for (k = n; k >= i; --k, ++m1pi, ++m1pj) - *m1pj -= q * (*m1pi); + for (k = n; k >= i; --k, ++m1pi, ++m1pj) + *m1pj -= q * (*m1pi); - double * m2pi = &m2.Elem(i, 1); - double * m2pj = &m2.Elem(j, 1); + double * m2pi = &m2.Elem(i, 1); + double * m2pj = &m2.Elem(j, 1); - for (k = i; k > 0; --k, ++m2pi, ++m2pj) - *m2pj -= q * (*m2pi); + for (k = i; k > 0; --k, ++m2pi, ++m2pj) + *m2pj -= q * (*m2pi); - // for (k = 1; k <= n; k++) - // { - // m1.Elem(j, k) -= q * m1.Elem(i, k); - // m2.Elem(j, k) -= q * m2.Elem(i, k); - // } + // for (k = 1; k <= n; k++) + // { + // m1.Elem(j, k) -= q * m1.Elem(i, k); + // m2.Elem(j, k) -= q * m2.Elem(i, k); + // } - } - } + } + } - for (i = n; i >= 1; i--) - { - // (*mycout) << "+" << flush; - for (j = 1; j < i; j++) + for (i = n; i >= 1; i--) + { + // (*mycout) << "+" << flush; + for (j = 1; j < i; j++) { - q = m1.Elem(j, i); + q = m1.Elem(j, i); - double * m2pi = &m2.Elem(i, 1); - double * m2pj = &m2.Elem(j, 1); + double * m2pi = &m2.Elem(i, 1); + double * m2pj = &m2.Elem(j, 1); - for (k = n; k > 0; --k, ++m2pi, ++m2pj) - *m2pj -= q * (*m2pi); + for (k = n; k > 0; --k, ++m2pi, ++m2pj) + *m2pj -= q * (*m2pi); - // for (k = 1; k <= n; k++) - // { - // m1.Elem(j, k) -= q * m1.Elem(i, k); - // m2.Elem(j, k) -= q * m2.Elem(i, k); - // } + // for (k = 1; k <= n; k++) + // { + // m1.Elem(j, k) -= q * m1.Elem(i, k); + // m2.Elem(j, k) -= q * m2.Elem(i, k); + // } } - } + } - if (m2.Symmetric()) - { - for (i = 1; i <= n; i++) + if (m2.Symmetric()) + { + for (i = 1; i <= n; i++) 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 n2 = a.Width(); - int i, j, k; - double sum; - const double *p, *q, *p0; + int n1 = a.Height(); + int n2 = a.Width(); + int i, j, k; + double sum; + const double *p, *q, *p0; - if (m2.Height() != n1 || m2.Width() != n1) - { - (*myerr) << "CalcAAt: sizes don't fit" << endl; - return; - } - - for (i = 1; i <= n1; i++) - { - sum = 0; - p = &a.ConstElem(i, 1); - for (k = 1; k <= n2; k++) + if (m2.Height() != n1 || m2.Width() != n1) { - sum += *p * *p; - p++; + (*myerr) << "CalcAAt: sizes don't fit" << endl; + return; } - m2.Set(i, i, sum); - p0 = &a.ConstElem(i, 1); - q = a.data; - for (j = 1; j < i; j++) + for (i = 1; i <= n1; i++) { - sum = 0; - p = p0; + sum = 0; + p = &a.ConstElem(i, 1); + for (k = 1; k <= n2; k++) + { + sum += *p * *p; + p++; + } + m2.Set(i, i, sum); - for (k = 1; k <= n2; k++) + p0 = &a.ConstElem(i, 1); + q = a.data; + for (j = 1; j < i; j++) + { + sum = 0; + p = p0; + + for (k = 1; k <= n2; k++) + { + sum += *p * *q; + p++; + q++; + } + m2.Set(i, j, sum); + m2.Set(j, i, sum); + } + } + } + + + + void CalcAtA (const DenseMatrix & a, DenseMatrix & m2) + { + int n1 = a.Height(); + int n2 = a.Width(); + int i, j, k; + double sum; + + if (m2.Height() != n2 || m2.Width() != n2) + { + (*myerr) << "CalcAtA: sizes don't fit" << endl; + return; + } + + for (i = 1; i <= n2; i++) + for (j = 1; j <= n2; j++) { - sum += *p * *q; - p++; - q++; + sum = 0; + for (k = 1; k <= n1; k++) + sum += a.Get(k, i) * a.Get(k, j); + m2.Elem(i, j) = sum; } - m2.Set(i, j, sum); - m2.Set(j, i, sum); - } - } } -#ifdef ABC -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) + void CalcABt (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2) { - int n1 = a.Height(); - int n2 = a.Width(); - int i, j, k; - double sum; + int n1 = a.Height(); + int n2 = a.Width(); + int n3 = b.Height(); + int i, j, k; + double sum; - if (m2.Height() != n2 || m2.Width() != n2) - { - (*myerr) << "CalcAtA: sizes don't fit" << endl; - return; - } + if (m2.Height() != n1 || m2.Width() != n3 || b.Width() != n2) + { + (*myerr) << "CalcABt: sizes don't fit" << endl; + return; + } - for (i = 1; i <= n2; i++) - for (j = 1; j <= n2; j++) + double * pm2 = &m2.Elem(1, 1); + const double * pa1 = &a.Get(1, 1); + + for (i = 1; i <= n1; i++) + { + const double * pb = &b.Get(1, 1); + for (j = 1; j <= n3; j++) + { + sum = 0; + const double * pa = pa1; + + for (k = 1; k <= n2; k++) + { + sum += *pa * *pb; + pa++; pb++; + } + + *pm2 = sum; + pm2++; + } + pa1 += n2; + } + } + + + void CalcAtB (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2) + { + int n1 = a.Height(); + int n2 = a.Width(); + int n3 = b.Width(); + int i, j, k; + + if (m2.Height() != n2 || m2.Width() != n3 || b.Height() != n1) + { + (*myerr) << "CalcAtB: sizes don't fit" << endl; + return; + } + + for (i = 1; i <= n2 * n3; i++) + m2.data[i-1] = 0; + + for (i = 1; i <= n1; i++) + for (j = 1; j <= n2; j++) + { + const double va = a.Get(i, j); + double * pm2 = &m2.Elem(j, 1); + const double * pb = &b.Get(i, 1); + + for (k = 1; k <= n3; ++k, ++pm2, ++pb) + *pm2 += va * *pb; + // for (k = 1; k <= n3; k++) + // m2.Elem(j, k) += va * b.Get(i, k); + } + /* + for (i = 1; i <= n2; i++) + for (j = 1; j <= n3; j++) { sum = 0; for (k = 1; k <= n1; k++) - sum += a.Get(k, i) * a.Get(k, j); + sum += a.Get(k, i) * b.Get(k, j); m2.Elem(i, j) = sum; } + */ } @@ -720,391 +771,307 @@ void CalcAtA (const DenseMatrix & a, DenseMatrix & m2) -void CalcABt (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2) + + DenseMatrix operator* (const DenseMatrix & m1, const DenseMatrix & m2) { - int n1 = a.Height(); - int n2 = a.Width(); - int n3 = b.Height(); - int i, j, k; - double sum; + DenseMatrix temp (m1.Height(), m2.Width()); - if (m2.Height() != n1 || m2.Width() != n3 || b.Width() != n2) - { - (*myerr) << "CalcABt: sizes don't fit" << endl; - return; - } - - double * pm2 = &m2.Elem(1, 1); - const double * pa1 = &a.Get(1, 1); - - for (i = 1; i <= n1; i++) - { - const double * pb = &b.Get(1, 1); - for (j = 1; j <= n3; j++) - { - sum = 0; - const double * pa = pa1; - - for (k = 1; k <= n2; k++) - { - sum += *pa * *pb; - pa++; pb++; - } - - *pm2 = sum; - pm2++; - } - pa1 += n2; - } - } - - -void CalcAtB (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2) - { - int n1 = a.Height(); - int n2 = a.Width(); - int n3 = b.Width(); - int i, j, k; - - if (m2.Height() != n2 || m2.Width() != n3 || b.Height() != n1) - { - (*myerr) << "CalcAtB: sizes don't fit" << endl; - return; - } - - for (i = 1; i <= n2 * n3; i++) - m2.data[i-1] = 0; - - for (i = 1; i <= n1; i++) - for (j = 1; j <= n2; j++) + if (m1.Width() != m2.Height()) { - const double va = a.Get(i, j); - double * pm2 = &m2.Elem(j, 1); - const double * pb = &b.Get(i, 1); - - for (k = 1; k <= n3; ++k, ++pm2, ++pb) - *pm2 += va * *pb; - // for (k = 1; k <= n3; k++) - // m2.Elem(j, k) += va * b.Get(i, k); + (*myerr) << "DenseMatrix :: operator*: Matrix Size does not fit" << endl; } - /* - for (i = 1; i <= n2; i++) - for (j = 1; j <= n3; j++) + else if (temp.Height() != m1.Height()) { - sum = 0; - for (k = 1; k <= n1; k++) - sum += a.Get(k, i) * b.Get(k, j); - m2.Elem(i, j) = sum; + (*myerr) << "DenseMatrix :: operator*: temp not allocated" << endl; } - */ + else + { + Mult (m1, m2, temp); + } + return temp; } - - - - - -DenseMatrix operator* (const DenseMatrix & m1, const DenseMatrix & m2) + void Mult (const DenseMatrix & m1, const DenseMatrix & m2, DenseMatrix & m3) { - DenseMatrix temp (m1.Height(), m2.Width()); + double sum; + double *p1, *p1s, *p1sn, *p1snn, *p2, *p2s, *p2sn, *p3; - if (m1.Width() != m2.Height()) - { - (*myerr) << "DenseMatrix :: operator*: Matrix Size does not fit" << endl; - } - else if (temp.Height() != m1.Height()) - { - (*myerr) << "DenseMatrix :: operator*: temp not allocated" << endl; - } - else - { - Mult (m1, m2, temp); - } - return temp; - } + if (m1.Width() != m2.Height() || m1.Height() != m3.Height() || + m2.Width() != m3.Width() ) + { + (*myerr) << "DenseMatrix :: Mult: Matrix Size does not fit" << endl; + (*myerr) << "m1: " << m1.Height() << " x " << m1.Width() << endl; + (*myerr) << "m2: " << m2.Height() << " x " << m2.Width() << endl; + (*myerr) << "m3: " << m3.Height() << " x " << m3.Width() << endl; + return; + } + /* + else if (m1.Symmetric() || m2.Symmetric() || m3.Symmetric()) + { + (*myerr) << "DenseMatrix :: Mult: not implemented for symmetric matrices" << endl; + return; + } + */ + else + { + // int i, j, k; + int n1 = m1.Height(); + int n2 = m2.Width(); + int n3 = m1.Width(); + /* + for (i = n1 * n2-1; i >= 0; --i) + m3.data[i] = 0; -void Mult (const DenseMatrix & m1, const DenseMatrix & m2, DenseMatrix & m3) - { - double sum; - double *p1, *p1s, *p1sn, *p1snn, *p2, *p2s, *p2sn, *p3; - - if (m1.Width() != m2.Height() || m1.Height() != m3.Height() || - m2.Width() != m3.Width() ) - { - (*myerr) << "DenseMatrix :: Mult: Matrix Size does not fit" << endl; - (*myerr) << "m1: " << m1.Height() << " x " << m1.Width() << endl; - (*myerr) << "m2: " << m2.Height() << " x " << m2.Width() << endl; - (*myerr) << "m3: " << m3.Height() << " x " << m3.Width() << endl; - return; - } - /* - else if (m1.Symmetric() || m2.Symmetric() || m3.Symmetric()) - { - (*myerr) << "DenseMatrix :: Mult: not implemented for symmetric matrices" << endl; - return; - } - */ - else - { - // int i, j, k; - int n1 = m1.Height(); - int n2 = m2.Width(); - int n3 = m1.Width(); - - /* - for (i = n1 * n2-1; i >= 0; --i) - m3.data[i] = 0; - - const double * pm1 = &m1.Get(1, 1); - for (i = 1; i <= n1; i++) - { + const double * pm1 = &m1.Get(1, 1); + for (i = 1; i <= n1; i++) + { const double * pm2 = &m2.Get(1, 1); double * pm3i = &m3.Elem(i, 1); for (j = 1; j <= n3; j++) - { - const double vm1 = *pm1; - ++pm1; - // const double vm1 = m1.Get(i, j); - double * pm3 = pm3i; - // const double * pm2 = &m2.Get(j, 1); + { + const double vm1 = *pm1; + ++pm1; + // const double vm1 = m1.Get(i, j); + double * pm3 = pm3i; + // const double * pm2 = &m2.Get(j, 1); - for (k = 0; k < n2; k++) - { - *pm3 += vm1 * *pm2; - ++pm2; - ++pm3; - } + for (k = 0; k < n2; k++) + { + *pm3 += vm1 * *pm2; + ++pm2; + ++pm3; + } - // for (k = 1; k <= n2; k++) - // m3.Elem(i, k) += m1.Get(i, j) * m2.Get(j, k); - } - } + // for (k = 1; k <= n2; k++) + // m3.Elem(i, k) += m1.Get(i, j) * m2.Get(j, k); + } + } */ - /* - for (i = 1; i <= n1; i++) - for (j = 1; j <= n2; j++) + /* + for (i = 1; i <= n1; i++) + for (j = 1; j <= n2; j++) { - sum = 0; - for (k = 1; k <= n3; k++) - sum += m1.Get(i, k) * m2.Get(k, j); - m3.Set(i, j, sum); + sum = 0; + for (k = 1; k <= n3; k++) + sum += m1.Get(i, k) * m2.Get(k, j); + m3.Set(i, j, sum); } - */ + */ - /* - for (i = 1; i <= n1; i++) - { + /* + for (i = 1; i <= n1; i++) + { const double pm1i = &m1.Get(i, 1); const double pm2j = &m2.Get(1, 1); for (j = 1; j <= n2; j++) - { - double sum = 0; - const double * pm1 = pm1i; - const double * pm2 = pm2j; - pm2j++; + { + double sum = 0; + const double * pm1 = pm1i; + const double * pm2 = pm2j; + pm2j++; - for (k = 1; k <= n3; k++) - { - sum += *pm1 * *pm2; - ++pm1; - pm2 += n2; - } + for (k = 1; k <= n3; k++) + { + sum += *pm1 * *pm2; + ++pm1; + pm2 += n2; + } - m3.Set (i, j, sum); - } - } + m3.Set (i, j, sum); + } + } */ - p3 = m3.data; - p1s = m1.data; - p2sn = m2.data + n2; - p1snn = p1s + n1 * n3; + p3 = m3.data; + p1s = m1.data; + p2sn = m2.data + n2; + p1snn = p1s + n1 * n3; - while (p1s != p1snn) - { - p1sn = p1s + n3; - p2s = m2.data; + while (p1s != p1snn) + { + p1sn = p1s + n3; + p2s = m2.data; - while (p2s != p2sn) - { - sum = 0; - p1 = p1s; - p2 = p2s; - p2s++; + while (p2s != p2sn) + { + sum = 0; + p1 = p1s; + p2 = p2s; + p2s++; - while (p1 != p1sn) - { - sum += *p1 * *p2; - p1++; - p2 += n2; - } - *p3++ = sum; - } - p1s = p1sn; - } - } + while (p1 != p1sn) + { + sum += *p1 * *p2; + p1++; + p2 += n2; + } + *p3++ = sum; + } + p1s = p1sn; + } + } } -DenseMatrix operator+ (const DenseMatrix & m1, const DenseMatrix & m2) + DenseMatrix operator+ (const DenseMatrix & m1, const DenseMatrix & m2) { - DenseMatrix temp (m1.Height(), m1.Width()); - int i, j; + DenseMatrix temp (m1.Height(), m1.Width()); + int i, j; - if (m1.Width() != m2.Width() || m1.Height() != m2.Height()) - { - (*myerr) << "BaseMatrix :: operator+: Matrix Size does not fit" << endl; - } - else if (temp.Height() != m1.Height()) - { - (*myerr) << "BaseMatrix :: operator+: temp not allocated" << endl; - } - else - { - for (i = 1; i <= m1.Height(); i++) - for (j = 1; j <= m1.Width(); j++) - { - temp.Set(i, j, m1.Get(i, j) + m2.Get(i, j)); - } - } - return temp; + if (m1.Width() != m2.Width() || m1.Height() != m2.Height()) + { + (*myerr) << "BaseMatrix :: operator+: Matrix Size does not fit" << endl; + } + else if (temp.Height() != m1.Height()) + { + (*myerr) << "BaseMatrix :: operator+: temp not allocated" << endl; + } + else + { + for (i = 1; i <= m1.Height(); i++) + for (j = 1; j <= m1.Width(); j++) + { + temp.Set(i, j, m1.Get(i, j) + m2.Get(i, j)); + } + } + return temp; } -void Transpose (const DenseMatrix & m1, DenseMatrix & m2) -{ - int w = m1.Width(); - int h = m1.Height(); - int i, j; - - m2.SetSize (w, h); - - double * pm2 = &m2.Elem(1, 1); - for (j = 1; j <= w; j++) - { - const double * pm1 = &m1.Get(1, j); - for (i = 1; i <= h; i++) - { - *pm2 = *pm1; - pm2 ++; - pm1 += w; - } - } -} - - -/* -void DenseMatrix :: Mult (const Vector & v, Vector & prod) const + void Transpose (const DenseMatrix & m1, DenseMatrix & m2) { - double sum, val; - const double * mp, * sp; - double * dp; - // const Vector & v = bv.CastToVector(); - // Vector & prod = bprod.CastToVector(); - + int w = m1.Width(); + int h = m1.Height(); + int i, j; - int n = Height(); - int m = Width(); + m2.SetSize (w, h); - if (prod.Size() != n) - prod.SetSize (n); - -#ifdef DEVELOP - if (!n) - { - cout << "DenseMatrix::Mult mheight = 0" << endl; - } - if (!m) - { - cout << "DenseMatrix::Mult mwidth = 0" << endl; - } - - if (m != v.Size()) - { - (*myerr) << "\nMatrix and Vector don't fit" << endl; - } - else if (Height() != prod.Size()) - { - (*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl; - } - else -#endif - { - if (Symmetric()) - { - int i, j; - - - for (i = 1; i <= n; i++) - { - sp = &v.Get(1); - dp = &prod.Elem(1); - mp = &Get(i, 1); - - val = v.Get(i); - sum = Get(i, i) * val; - - for (j = 1; j < i; ++j, ++mp, ++sp, ++dp) - { - sum += *mp * *sp; - *dp += val * *mp; - } - - prod.Elem(i) = sum; - } - } - else - { - mp = data; - dp = &prod.Elem(1); - for (int i = 1; i <= n; i++) - { - sum = 0; - sp = &v.Get(1); - - for (int j = 1; j <= m; j++) - { - // sum += Get(i,j) * v.Get(j); - sum += *mp * *sp; - mp++; - sp++; - } - - // prod.Set (i, sum); - *dp = sum; - dp++; - } - } - } + double * pm2 = &m2.Elem(1, 1); + for (j = 1; j <= w; j++) + { + const double * pm1 = &m1.Get(1, j); + for (i = 1; i <= h; i++) + { + *pm2 = *pm1; + pm2 ++; + pm1 += w; + } + } } -*/ -void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const -{ - // const Vector & v = (const Vector&)bv; // .CastToVector(); - // Vector & prod = (Vector & )bprod; // .CastToVector(); /* - if (Height() != v.Size()) + void DenseMatrix :: Mult (const Vector & v, Vector & prod) const + { + double sum, val; + const double * mp, * sp; + double * dp; + // const Vector & v = bv.CastToVector(); + // Vector & prod = bprod.CastToVector(); + + + int n = Height(); + int m = Width(); + + if (prod.Size() != n) + prod.SetSize (n); + + #ifdef DEVELOP + if (!n) + { + cout << "DenseMatrix::Mult mheight = 0" << endl; + } + if (!m) + { + cout << "DenseMatrix::Mult mwidth = 0" << endl; + } + + if (m != v.Size()) { (*myerr) << "\nMatrix and Vector don't fit" << endl; } - else if (Width() != prod.Size()) + else if (Height() != prod.Size()) { (*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl; } - else + else + #endif + { + if (Symmetric()) + { + int i, j; + + + for (i = 1; i <= n; i++) + { + sp = &v.Get(1); + dp = &prod.Elem(1); + mp = &Get(i, 1); + + val = v.Get(i); + sum = Get(i, i) * val; + + for (j = 1; j < i; ++j, ++mp, ++sp, ++dp) + { + sum += *mp * *sp; + *dp += val * *mp; + } + + prod.Elem(i) = sum; + } + } + else + { + mp = data; + dp = &prod.Elem(1); + for (int i = 1; i <= n; i++) + { + sum = 0; + sp = &v.Get(1); + + for (int j = 1; j <= m; j++) + { + // sum += Get(i,j) * v.Get(j); + sum += *mp * *sp; + mp++; + sp++; + } + + // prod.Set (i, sum); + *dp = sum; + dp++; + } + } + } + } */ + + void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const + { + // const Vector & v = (const Vector&)bv; // .CastToVector(); + // Vector & prod = (Vector & )bprod; // .CastToVector(); + + /* + if (Height() != v.Size()) + { + (*myerr) << "\nMatrix and Vector don't fit" << endl; + } + else if (Width() != prod.Size()) + { + (*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl; + } + else + */ { int i, j; int w = Width(), h = Height(); @@ -1112,7 +1079,7 @@ void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const prod.SetSize (w); const double * pmat = &Get(1, 1); - const double * pv = &v.Get(1); + const double * pv = &v(0); prod = 0; @@ -1121,7 +1088,7 @@ void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const double val = *pv; ++pv; - double * pprod = &prod.Elem(1); + double * pprod = &prod(0); for (j = w-1; j >= 0; --j, ++pmat, ++pprod) { @@ -1130,314 +1097,313 @@ void DenseMatrix :: MultTrans (const Vector & v, Vector & prod) const } /* - double sum; + double sum; - for (i = 1; i <= Width(); i++) + for (i = 1; i <= Width(); i++) { - sum = 0; + sum = 0; - for (int j = 1; j <= Height(); j++) - sum += Get(j, i) * v.Get(j); + for (int j = 1; j <= Height(); j++) + sum += Get(j, i) * v.Get(j); - prod.Set (i, sum); + prod.Set (i, sum); } */ } } -void DenseMatrix :: Residuum (const Vector & x, const Vector & b, - Vector & res) const + void DenseMatrix :: Residuum (const Vector & x, const Vector & b, + Vector & res) const { - double sum; - // const Vector & x = bx.CastToVector(); - // const Vector & b = bb.CastToVector(); - // Vector & res = bres.CastToVector(); + double sum; + // const Vector & x = bx.CastToVector(); + // const Vector & b = bb.CastToVector(); + // Vector & res = bres.CastToVector(); - res.SetSize (Height()); + res.SetSize (Height()); - if (Width() != x.Size() || Height() != b.Size()) - { - (*myerr) << "\nMatrix and Vector don't fit" << endl; - } - else if (Height() != res.Size()) - { - (*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl; - } - else - { - int i, j; - int h = Height(); - int w = Width(); - const double * mp = &Get(1, 1); + if (Width() != x.Size() || Height() != b.Size()) + { + (*myerr) << "\nMatrix and Vector don't fit" << endl; + } + else if (Height() != res.Size()) + { + (*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl; + } + else + { + int h = Height(); + int w = Width(); + const double * mp = &Get(1, 1); - for (i = 1; i <= h; i++) - { - sum = b.Get(i); - const double * xp = &x.Get(1); + for (int i = 1; i <= h; i++) + { + sum = b(i-1); + const double * xp = &x(0); - for (j = 1; j <= w; ++j, ++mp, ++xp) - sum -= *mp * *xp; + for (int j = 1; j <= w; ++j, ++mp, ++xp) + sum -= *mp * *xp; - res.Elem(i) = sum; - } - } + res(i-1) = sum; + } + } } #ifdef ABC -double DenseMatrix :: EvaluateBilinearform (const Vector & hx) const + double DenseMatrix :: EvaluateBilinearform (const Vector & hx) const { - double sum = 0, hsum; - // const Vector & hx = x.CastToVector(); - int i, j; + double sum = 0, hsum; + // const Vector & hx = x.CastToVector(); + int i, j; - if (Width() != hx.Size() || Height() != hx.Size()) - { - (*myerr) << "Matrix::EvaluateBilinearForm: sizes don't fit" << endl; - } - else - { - for (i = 1; i <= Height(); i++) + if (Width() != hx.Size() || Height() != hx.Size()) { - hsum = 0; - for (j = 1; j <= Height(); j++) - { - hsum += Get(i, j) * hx.Get(j); - } - sum += hsum * hx.Get(i); + (*myerr) << "Matrix::EvaluateBilinearForm: sizes don't fit" << endl; + } + else + { + for (i = 1; i <= Height(); i++) + { + hsum = 0; + for (j = 1; j <= Height(); j++) + { + hsum += Get(i, j) * hx.Get(j); + } + sum += hsum * hx.Get(i); + } } - } -// testout << "sum = " << sum << endl; - return sum; + // testout << "sum = " << sum << endl; + return sum; } -void DenseMatrix :: MultElementMatrix (const Array & pnum, - const Vector & hx, Vector & hy) + void DenseMatrix :: MultElementMatrix (const Array & pnum, + const Vector & hx, Vector & hy) { - int i, j; - // const Vector & hx = x.CastToVector(); - // Vector & hy = y.CastToVector(); + int i, j; + // const Vector & hx = x.CastToVector(); + // Vector & hy = y.CastToVector(); - if (Symmetric()) - { - for (i = 1; i <= Height(); i++) + if (Symmetric()) { - for (j = 1; j < i; j++) - { - hy.Elem(pnum.Get(i)) += Get(i, j) * hx.Get(pnum.Get(j)); - hy.Elem(pnum.Get(j)) += Get(i, j) * hx.Get(pnum.Get(i)); - } - hy.Elem(pnum.Get(j)) += Get(i, i) * hx.Get(pnum.Get(i)); + for (i = 1; i <= Height(); i++) + { + for (j = 1; j < i; j++) + { + hy.Elem(pnum.Get(i)) += Get(i, j) * hx.Get(pnum.Get(j)); + hy.Elem(pnum.Get(j)) += Get(i, j) * hx.Get(pnum.Get(i)); + } + hy.Elem(pnum.Get(j)) += Get(i, i) * hx.Get(pnum.Get(i)); + } } - } - else - for (i = 1; i <= Height(); i++) - for (j = 1; j <= Width(); j++) - hy.Elem(pnum.Get(i)) += Get(i, j) * hx.Get(pnum.Get(j)); + else + for (i = 1; i <= Height(); i++) + for (j = 1; j <= Width(); j++) + hy.Elem(pnum.Get(i)) += Get(i, j) * hx.Get(pnum.Get(j)); } -void DenseMatrix :: MultTransElementMatrix (const Array & pnum, - const Vector & hx, Vector & hy) + void DenseMatrix :: MultTransElementMatrix (const Array & pnum, + const Vector & hx, Vector & hy) { - int i, j; - // const Vector & hx = x.CastToVector(); - // Vector & hy = y.CastToVector(); + int i, j; + // const Vector & hx = x.CastToVector(); + // Vector & hy = y.CastToVector(); - if (Symmetric()) - MultElementMatrix (pnum, hx, hy); - else - for (i = 1; i <= Height(); i++) - for (j = 1; j <= Width(); j++) - hy.Elem(pnum.Get(i)) += Get(j, i) * hx.Get(pnum.Get(j)); + if (Symmetric()) + MultElementMatrix (pnum, hx, hy); + else + for (i = 1; i <= Height(); i++) + for (j = 1; j <= Width(); j++) + hy.Elem(pnum.Get(i)) += Get(j, i) * hx.Get(pnum.Get(j)); } #endif -void DenseMatrix :: Solve (const Vector & v, Vector & sol) const -{ - DenseMatrix temp (*this); - temp.SolveDestroy (v, sol); -} - - -void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol) + void DenseMatrix :: Solve (const Vector & v, Vector & sol) const { - double q; - - if (Width() != Height()) - { - (*myerr) << "SolveDestroy: Matrix not square"; - return; - } - if (Width() != v.Size()) - { - (*myerr) << "SolveDestroy: Matrix and Vector don't fit"; - return; - } - - sol = v; - if (Height() != sol.Size()) - { - (*myerr) << "SolveDestroy: Solution Vector not ok"; - return; - } - - - if (0 /* Symmetric() */) - { - - // Cholesky factorization - - int i, j, k, n; - n = Height(); - - // Cholesky - - double x; - Vector p(n); - - for (i = 1; i <= n; i++) - for (j = 1; j < i; j++) - Elem(j, i) = Get(i, j); - - for (i = 1; i <= n; i++) - { - // (*mycout) << "." << flush; - for (j = i; j <= n; j++) - { - x = Get(i, j); - - const double * pik = &Get(i, 1); - const double * pjk = &Get(j, 1); - - for (k = i-2; k >= 0; --k, ++pik, ++pjk) - x -= (*pik) * (*pjk); - - // for (k = i-1; k >= 1; --k) - // x -= Get(j, k) * Get(i, k); - - if (i == j) - { - if (x <= 0) - { - cerr << "Matrix indefinite" << endl; - return; - } - - p.Elem(i) = 1 / sqrt(x); - } - else - { - Elem(j, i) = x * p.Get(i); - } - } - } - - for (i = 1; i <= n; i++) - Elem(i, i) = 1 / p.Get(i); - - // A = L L^t - // L stored in left-lower triangle - - - sol = v; - - // Solve L sol = sol - - for (i = 1; i <= n; i++) - { - double val = sol.Get(i); - - const double * pij = &Get(i, 1); - const double * solj = &sol.Get(1); - - for (j = 1; j < i; j++, ++pij, ++solj) - val -= *pij * *solj; - // for (j = 1; j < i; j++) - // val -= Get(i, j) * sol.Get(j); - - sol.Elem(i) = val / Get(i, i); - } - - // Solve L^t sol = sol - - for (i = n; i >= 1; i--) - { - double val = sol.Get(i) / Get(i, i); - sol.Elem(i) = val; - - double * solj = &sol.Elem(1); - const double * pij = &Get(i, 1); - - for (j = 1; j < i; ++j, ++pij, ++solj) - *solj -= val * *pij; - // for (j = 1; j < i; j++) - // sol.Elem(j) -= Get(i, j) * val; - } - - - } - else - { - // (*mycout) << "gauss" << endl; - int i, j, k, n = Height(); - for (i = 1; i <= n; i++) - { - for (j = i+1; j <= n; j++) - { - q = Get(j,i) / Get(i,i); - if (q) - { - const double * pik = &Get(i, i+1); - double * pjk = &Elem(j, i+1); - - for (k = i+1; k <= n; ++k, ++pik, ++pjk) - *pjk -= q * *pik; - - // for (k = i+1; k <= Height(); k++) - // Elem(j, k) -= q * Get(i,k); - - - sol.Elem(j) -= q * sol.Get(i); - } - } - } - - for (i = n; i >= 1; i--) - { - q = sol.Get(i); - for (j = i+1; j <= n; j++) - q -= Get(i,j) * sol.Get(j); - - sol.Elem(i) = q / Get(i,i); - } - } + DenseMatrix temp (*this); + temp.SolveDestroy (v, sol); } -/* -BaseMatrix * DenseMatrix :: Copy () const + void DenseMatrix :: SolveDestroy (const Vector & v, Vector & sol) { - return new DenseMatrix (*this); + double q; + + if (Width() != Height()) + { + (*myerr) << "SolveDestroy: Matrix not square"; + return; + } + if (Width() != v.Size()) + { + (*myerr) << "SolveDestroy: Matrix and Vector don't fit"; + return; + } + + sol = v; + if (Height() != sol.Size()) + { + (*myerr) << "SolveDestroy: Solution Vector not ok"; + return; + } + + + if (0 /* Symmetric() */) + { + + // Cholesky factorization + + int i, j, k, n; + n = Height(); + + // Cholesky + + double x; + Vector p(n); + + for (i = 1; i <= n; i++) + for (j = 1; j < i; j++) + Elem(j, i) = Get(i, j); + + for (i = 1; i <= n; i++) + { + // (*mycout) << "." << flush; + for (j = i; j <= n; j++) + { + x = Get(i, j); + + const double * pik = &Get(i, 1); + const double * pjk = &Get(j, 1); + + for (k = i-2; k >= 0; --k, ++pik, ++pjk) + x -= (*pik) * (*pjk); + + // for (k = i-1; k >= 1; --k) + // x -= Get(j, k) * Get(i, k); + + if (i == j) + { + if (x <= 0) + { + cerr << "Matrix indefinite" << endl; + return; + } + + p(i-1) = 1 / sqrt(x); + } + else + { + Elem(j, i) = x * p(i-1); + } + } + } + + for (int i = 1; i <= n; i++) + Elem(i, i) = 1 / p(i-1); + + // A = L L^t + // L stored in left-lower triangle + + + sol = v; + + // Solve L sol = sol + + for (int i = 1; i <= n; i++) + { + double val = sol(i-1); + + const double * pij = &Get(i, 1); + const double * solj = &sol(0); + + for (int j = 1; j < i; j++, ++pij, ++solj) + val -= *pij * *solj; + // for (j = 1; j < i; j++) + // val -= Get(i, j) * sol.Get(j); + + sol(i-1) = val / Get(i, i); + } + + // Solve L^t sol = sol + + for (int i = n; i >= 1; i--) + { + double val = sol(i-1) / Get(i, i); + sol(i-1) = val; + + double * solj = &sol(0); + const double * pij = &Get(i, 1); + + for (j = 1; j < i; ++j, ++pij, ++solj) + *solj -= val * *pij; + // for (j = 1; j < i; j++) + // sol.Elem(j) -= Get(i, j) * val; + } + + + } + else + { + // (*mycout) << "gauss" << endl; + int n = Height(); + for (int i = 1; i <= n; i++) + { + for (int j = i+1; j <= n; j++) + { + q = Get(j,i) / Get(i,i); + if (q) + { + const double * pik = &Get(i, i+1); + double * pjk = &Elem(j, i+1); + + for (int k = i+1; k <= n; ++k, ++pik, ++pjk) + *pjk -= q * *pik; + + // for (k = i+1; k <= Height(); k++) + // Elem(j, k) -= q * Get(i,k); + + + sol(j-1) -= q * sol(i-1); + } + } + } + + for (int i = n; i >= 1; i--) + { + q = sol(i-1); + for (int j = i+1; j <= n; j++) + q -= Get(i,j) * sol(j-1); + + sol(i-1) = q / Get(i,i); + } + } } -*/ - - -ostream & operator<< (ostream & ost, const DenseMatrix & m) -{ - for (int i = 0; i < m.Height(); i++) + /* + BaseMatrix * DenseMatrix :: Copy () const { - for (int j = 0; j < m.Width(); j++) - ost << m.Get(i+1,j+1) << " "; - ost << endl; + return new DenseMatrix (*this); } - return ost; -} + */ + + + + + ostream & operator<< (ostream & ost, const DenseMatrix & m) + { + for (int i = 0; i < m.Height(); i++) + { + for (int j = 0; j < m.Width(); j++) + ost << m.Get(i+1,j+1) << " "; + ost << endl; + } + return ost; + } diff --git a/libsrc/linalg/densemat.hpp b/libsrc/linalg/densemat.hpp index 9daec518..a6b1e042 100644 --- a/libsrc/linalg/densemat.hpp +++ b/libsrc/linalg/densemat.hpp @@ -2,7 +2,7 @@ #define FILE_DENSEMAT /**************************************************************************/ -/* File: densemat.hh */ +/* File: densemat.hpp */ /* Author: Joachim Schoeberl */ /* Date: 01. Oct. 94 */ /**************************************************************************/ @@ -12,11 +12,6 @@ */ -#include - - - - class DenseMatrix { protected: @@ -67,12 +62,8 @@ public: #ifdef DEBUG if (prod.Size() != height) { - cerr << "Mult: wrong vector size " << endl; - assert (1); - // prod.SetSize (height); + (*myerr) << "Mult: wrong vector size " << endl; } - - if (!height) { cout << "DenseMatrix::Mult height = 0" << endl; @@ -94,13 +85,13 @@ public: #endif { mp = data; - dp = &prod.Elem(1); - for (int i = 1; i <= height; i++) + dp = &prod(0); + for (int i = 0; i < height; i++) { 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 += *mp * *sp; diff --git a/libsrc/linalg/linalg.hpp b/libsrc/linalg/linalg.hpp index 7be9b4ad..95d0c823 100644 --- a/libsrc/linalg/linalg.hpp +++ b/libsrc/linalg/linalg.hpp @@ -10,7 +10,6 @@ /* Data types for basic linear algebra - more data types are found in linalgl.hpp The basic concepts include the data types diff --git a/libsrc/linalg/linopt.cpp b/libsrc/linalg/linopt.cpp index a5d381e6..dc5b53fa 100644 --- a/libsrc/linalg/linopt.cpp +++ b/libsrc/linalg/linopt.cpp @@ -37,9 +37,9 @@ void LinearOptimize (const DenseMatrix & a, const Vector & b, m.Elem(3, j) = a.Get(i3, j); } - rs.Elem(1) = b.Get(i1); - rs.Elem(2) = b.Get(i2); - rs.Elem(3) = b.Get(i3); + rs(0) = b(i1-1); + rs(1) = b(i2-1); + rs(2) = b(i3-1); if (fabs (m.Det()) < 1e-12) continue; @@ -59,9 +59,9 @@ void LinearOptimize (const DenseMatrix & a, const Vector & b, */ - double rmin = res.Elem(1); - for (int hi = 2; hi <= res.Size(); hi++) - if (res.Elem(hi) < rmin) rmin = res.Elem(hi); + double rmin = res(0); + for (int hi = 1; hi < res.Size(); hi++) + if (res(hi) < rmin) rmin = res(hi); if ( (f < fmin) && rmin >= -1e-8) { diff --git a/libsrc/linalg/linsearch.cpp b/libsrc/linalg/linsearch.cpp index 4c297fd9..dcdb6644 100644 --- a/libsrc/linalg/linsearch.cpp +++ b/libsrc/linalg/linsearch.cpp @@ -2,7 +2,7 @@ /* */ /* Problem: Liniensuche */ /* */ -/* Programmautor: Joachim Schöberl */ +/* Programmautor: Joachim Schöberl */ /* Matrikelnummer: 9155284 */ /* */ /* Algorithmus nach: */ @@ -98,36 +98,36 @@ void MinFunction :: ApproximateHesse (const Vector & x, double eps = 1e-6; 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.Elem(i) = x.Get(i) + eps; - hx.Elem(j) = x.Get(j) + eps; + hx(i) = x(i) + eps; + hx(j) = x(j) + eps; f11 = Func(hx); - hx.Elem(i) = x.Get(i) + eps; - hx.Elem(j) = x.Get(j) - eps; + hx(i) = x(i) + eps; + hx(j) = x(j) - eps; f12 = Func(hx); - hx.Elem(i) = x.Get(i) - eps; - hx.Elem(j) = x.Get(j) + eps; + hx(i) = x(i) - eps; + hx(j) = x(j) + eps; f21 = Func(hx); - hx.Elem(i) = x.Get(i) - eps; - hx.Elem(j) = x.Get(j) - eps; + hx(i) = x(i) - eps; + hx(j) = x(j) - eps; f22 = Func(hx); - hesse.Elem(i, j) = hesse.Elem(j, i) = + hesse(i, j) = hesse(j, i) = (f11 + f22 - f12 - f21) / (2 * eps * eps); } hx = x; f = Func(x); - hx.Elem(i) = x.Get(i) + eps; + hx(i) = x(i) + eps; f11 = Func(hx); - hx.Elem(i) = x.Get(i) - eps; + hx(i) = x(i) - eps; 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; } diff --git a/libsrc/linalg/polynomial.cpp b/libsrc/linalg/polynomial.cpp index f947e0a1..cc515aac 100644 --- a/libsrc/linalg/polynomial.cpp +++ b/libsrc/linalg/polynomial.cpp @@ -152,24 +152,6 @@ MaxUnitSquare () if (hv > maxv) maxv = hv; 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; - */ }; diff --git a/libsrc/linalg/vector.cpp b/libsrc/linalg/vector.cpp index 56c2059f..859147e4 100644 --- a/libsrc/linalg/vector.cpp +++ b/libsrc/linalg/vector.cpp @@ -6,160 +6,104 @@ namespace netgen { -double BaseVector :: shit = 0; + double BaseVector :: shit = 0; -// %FD Constructs a vector of length zero -BaseVector :: BaseVector () + // %FD Constructs a vector of length zero + BaseVector :: BaseVector () { - length = 0; + length = 0; } -// %FD Constructs a vector of given length -BaseVector :: BaseVector ( - INDEX alength // length of the vector - ) + // %FD Constructs a vector of given length + BaseVector :: BaseVector ( + INDEX alength // length of the vector + ) { - length = alength; + length = alength; } -// %FD Sets length of the vector, old vector will be destroyed -void -BaseVector :: SetLength ( - INDEX alength // new length of the vector - ) + // %FD Sets length of the vector, old vector will be destroyed + void + BaseVector :: SetLength ( + INDEX alength // new length of the vector + ) { - length = alength; + length = alength; } -// %FD Changes length of the vector, old values stay valid -void -BaseVector :: ChangeLength ( - INDEX alength // new length of the vector - ) + // %FD Changes length of the vector, old values stay valid + void + BaseVector :: ChangeLength ( + INDEX alength // new length of the vector + ) { - length = alength; + length = alength; } -// %FD { Write a vector with the help of the '<<' operator onto a stream } -ostream & // stream for further use -operator<< ( - ostream & s, // stream to write vector onto - const BaseVector & v // vector to write - ) + // %FD { Write a vector with the help of the '<<' operator onto a stream } + ostream & // stream for further use + operator<< ( + ostream & s, // stream to write vector onto + const BaseVector & v // vector to write + ) { - return v.Print (s); + return v.Print (s); } -// %FD{ Divides every component of the vector by the scalar c. -// The function checks for division by zero } -BaseVector & // result vector -BaseVector :: operator/= ( - double c // scalar to divide by - ) + // %FD{ Divides every component of the vector by the scalar c. + // The function checks for division by zero } + BaseVector & // result vector + BaseVector :: operator/= ( + double c // scalar to divide by + ) { - if (c) - return (*this) *= (1/c); - else - { - (*myerr) << "operator/=: division by zero" << endl; - return *this; - } + if (c) + return (*this) *= (1/c); + else + { + (*myerr) << "operator/=: division by zero" << endl; + return *this; + } } -// %FD Creates a copy of the object -BaseVector * // pointer to the new vector -BaseVector :: Copy () const + // %FD Creates a copy of the object + BaseVector * // pointer to the new vector + BaseVector :: Copy () const { - cerr << "Base_vector::Copy called" << endl << flush; - return NULL; + cerr << "Base_vector::Copy called" << endl << flush; + return NULL; } -void BaseVector :: GetElementVector (const Array & pnum, - BaseVector & elvec) const -{ - int i; - for (i = 1; i <= pnum.Size(); i++) - elvec(i) = (*this)(pnum.Get(i)); -} - -void BaseVector :: SetElementVector (const Array & pnum, - const BaseVector & elvec) -{ - int i; - for (i = 1; i <= pnum.Size(); i++) - (*this)(pnum.Get(i)) = elvec(i); -} - - -void BaseVector :: AddElementVector (const Array & pnum, - const BaseVector & elvec) -{ - int i; - for (i = 1; i <= pnum.Size(); i++) - (*this)(pnum.Get(i)) += elvec(i); -} - - - - - - - - - - - -TempVector :: ~TempVector () + void BaseVector :: GetElementVector (const Array & pnum, + BaseVector & elvec) const { - delete vec; + int i; + for (i = 1; i <= pnum.Size(); i++) + elvec(i) = (*this)(pnum.Get(i)); } -TempVector BaseVector :: operator+ (const BaseVector & v2) const + void BaseVector :: SetElementVector (const Array & pnum, + const BaseVector & elvec) { - return (*Copy()) += v2; - } - -TempVector BaseVector :: operator- (const BaseVector & v2) const - { - return (*Copy()) -= v2; - } - -TempVector BaseVector :: operator- () const - { - return (*Copy()) *= -1; + int i; + for (i = 1; i <= pnum.Size(); i++) + (*this)(pnum.Get(i)) = elvec(i); } -TempVector operator* (const BaseVector & v1, double scal) + void BaseVector :: AddElementVector (const Array & pnum, + const BaseVector & elvec) { - return (*v1.Copy()) *= scal; - } - -TempVector operator/ (const BaseVector & v1, double scal) - { - return (*v1.Copy()) /= scal; - } - - -TempVector operator* (double scal, const BaseVector & v1) - { - return v1 * scal; - } - - - - - -BaseVector * TempVector :: Copy () const - { - return vec->Copy(); + int i; + for (i = 1; i <= pnum.Size(); i++) + (*this)(pnum.Get(i)) += elvec(i); } @@ -171,616 +115,672 @@ BaseVector * TempVector :: Copy () const -double Vector :: shit = 0; -class clVecpool -{ -public: - Array vecs; - Array veclens; + TempVector :: ~TempVector () + { + delete vec; + } - ~clVecpool(); -}; + TempVector BaseVector :: operator+ (const BaseVector & v2) const + { + return (*Copy()) += v2; + } -clVecpool :: ~clVecpool() -{ - int i; - for (i = 1; i <= vecs.Size(); i++) - delete vecs.Elem(i); -} + TempVector BaseVector :: operator- (const BaseVector & v2) const + { + return (*Copy()) -= v2; + } -static clVecpool vecpool; + TempVector BaseVector :: operator- () const + { + return (*Copy()) *= -1; + } + + + TempVector operator* (const BaseVector & v1, double scal) + { + return (*v1.Copy()) *= scal; + } + + TempVector operator/ (const BaseVector & v1, double scal) + { + return (*v1.Copy()) /= scal; + } + + + TempVector operator* (double scal, const BaseVector & v1) + { + return v1 * scal; + } -static double * NewDouble (INDEX len) -{ - if (len < 10) - return new double[len]; - else - { - int i; - for (i = 1; i <= vecpool.veclens.Size(); i++) - if (vecpool.veclens.Get(i) == len) - { - double * hvec = vecpool.vecs.Get(i); - vecpool.vecs.DeleteElement(i); - vecpool.veclens.DeleteElement(i); - return hvec; - } + + BaseVector * TempVector :: Copy () const + { + return vec->Copy(); + } + + + + + + + + + + + double Vector :: shit = 0; + + class clVecpool + { + public: + Array vecs; + Array veclens; + + ~clVecpool(); + }; + + clVecpool :: ~clVecpool() + { + int i; + for (i = 1; i <= vecs.Size(); i++) + delete vecs.Elem(i); + } + + static clVecpool vecpool; + + + + static double * NewDouble (INDEX len) + { + if (len < 10) return new double[len]; - } -} + else + { + int i; + for (i = 1; i <= vecpool.veclens.Size(); i++) + if (vecpool.veclens.Get(i) == len) + { + double * hvec = vecpool.vecs.Get(i); + vecpool.vecs.DeleteElement(i); + vecpool.veclens.DeleteElement(i); + return hvec; + } -static void DeleteDouble (INDEX len, double * dp) -{ - if (len < 10) - delete [] dp; - else - { - vecpool.vecs.Append (dp); - vecpool.veclens.Append (len); - } -} - - - -Vector :: Vector () : BaseVector() - { - data = NULL; + return new double[len]; + } } -Vector :: Vector (INDEX alength) : BaseVector (alength) + static void DeleteDouble (INDEX len, double * dp) { - if (length) - { - // data = new double[length]; - data = NewDouble (length); + if (len < 10) + delete [] dp; + else + { + vecpool.vecs.Append (dp); + vecpool.veclens.Append (len); + } + } + + + + Vector :: Vector () : BaseVector() + { + data = NULL; + } + + Vector :: Vector (INDEX alength) : BaseVector (alength) + { + if (length) + { + // data = new double[length]; + data = NewDouble (length); + + if (!data) + { + length = 0; + (*myerr) << "Vector not allocated" << endl; + } + } + else + data = NULL; + } + + + Vector :: Vector (const Vector & v2) + { + length = v2.length; + + if (length) + { + // data = new double[length]; + data = NewDouble (length); + + if (data) + { + memcpy (data, v2.data, length * sizeof (double)); + } + else + { + length = 0; + (*myerr) << "Vector::Vector : Vector not allocated" << endl; + } + } + else + data = NULL; + } + + + Vector :: ~Vector () + { + // veclenfile << "~Vector delete: " << length << endl; + if (data) + { + DeleteDouble (length, data); + // delete [] data; + } + + } + + void Vector :: SetLength (INDEX alength) + { + if (length == alength) return; + + if (data) + { + DeleteDouble (length, data); + // delete [] data; + } + data = NULL; + length = alength; + + if (length == 0) return; + // data = new double[length]; + data = NewDouble (length); if (!data) { - length = 0; - (*myerr) << "Vector not allocated" << endl; + length = 0; + (*myerr) << "Vector::SetLength: Vector not allocated" << endl; } - } - else - data = NULL; } - -Vector :: Vector (const Vector & v2) + void Vector :: ChangeLength (INDEX alength) { - length = v2.length; - - if (length) - { - // data = new double[length]; - data = NewDouble (length); - - if (data) - { - memcpy (data, v2.data, length * sizeof (double)); - } - else - { - length = 0; - (*myerr) << "Vector::Vector : Vector not allocated" << endl; - } - } - else - data = NULL; - } - - -Vector :: ~Vector () -{ - // veclenfile << "~Vector delete: " << length << endl; - if (data) - { - DeleteDouble (length, data); - // delete [] data; - } - -} - -void Vector :: SetLength (INDEX alength) - { - if (length == alength) return; - - if (data) - { - DeleteDouble (length, data); - // delete [] data; - } - data = NULL; - length = alength; - - if (length == 0) return; - // data = new double[length]; - data = NewDouble (length); - - if (!data) - { - length = 0; - (*myerr) << "Vector::SetLength: Vector not allocated" << endl; - } - } - -void Vector :: ChangeLength (INDEX alength) -{ - (*mycout) << "Vector::ChangeLength called" << endl; - if (length == alength) return; + (*mycout) << "Vector::ChangeLength called" << endl; + if (length == alength) return; - if (alength == 0) - { - // delete [] data; - DeleteDouble (length, data); - length = 0; - return; - } + if (alength == 0) + { + // delete [] data; + DeleteDouble (length, data); + length = 0; + return; + } - double * olddata = data; + double * olddata = data; + + data = NewDouble (alength); + // data = new double[alength]; + if (!data) + { + length = 0; + (*myerr) << "Vector::SetLength: Vector not allocated" << endl; + delete [] olddata; + } + + memcpy (data, olddata, min2(alength, length)); - data = NewDouble (alength); - // data = new double[alength]; - if (!data) - { - length = 0; - (*myerr) << "Vector::SetLength: Vector not allocated" << endl; delete [] olddata; - } - - memcpy (data, olddata, min2(alength, length)); - - delete [] olddata; - length = alength; + length = alength; } -/// NEW RM -void Vector::SetBlockLength (INDEX /* blength */) -{ - MyError("BaseVector::SetBlockLength was called for a Vector"); -} - - -double & Vector :: operator() (INDEX i) + /// NEW RM + void Vector::SetBlockLength (INDEX /* blength */) { - if (i >= 1 && i <= length) return Elem(i); - else (*myerr) << "\nindex " << i << " out of range (" - << 1 << "," << Length() << ")\n"; - return shit; + MyError("BaseVector::SetBlockLength was called for a Vector"); } -double Vector :: operator() (INDEX i) const + + double & Vector :: operator() (INDEX i) { - if (i >= 1 && i <= length) return Get(i); - else (*myerr) << "\nindex " << i << " out of range (" - << 1 << "," << Length() << ")\n" << flush; - return shit; + if (i >= 1 && i <= length) return Elem(i); + else (*myerr) << "\nindex " << i << " out of range (" + << 1 << "," << Length() << ")\n"; + return shit; + } + + double Vector :: operator() (INDEX i) const + { + if (i >= 1 && i <= length) return Get(i); + else (*myerr) << "\nindex " << i << " out of range (" + << 1 << "," << Length() << ")\n" << flush; + return shit; } -double Vector :: SupNorm () const + double Vector :: SupNorm () const { - double sup = 0; + double sup = 0; - for (INDEX i = 1; i <= Length(); i++) - if (fabs (Get(i)) > sup) - sup = fabs(Get(i)); + for (INDEX i = 1; i <= Length(); i++) + if (fabs (Get(i)) > sup) + sup = fabs(Get(i)); - return sup; + return sup; } -double Vector :: L2Norm () const + double Vector :: L2Norm () const { - double sum = 0; + double sum = 0; - for (INDEX i = 1; i <= Length(); i++) - sum += Get(i) * Get(i); + for (INDEX i = 1; i <= Length(); i++) + sum += Get(i) * Get(i); - return sqrt (sum); + return sqrt (sum); } -double Vector :: L1Norm () const + double Vector :: L1Norm () const { - double sum = 0; + double sum = 0; - for (INDEX i = 1; i <= Length(); i++) - sum += fabs (Get(i)); + for (INDEX i = 1; i <= Length(); i++) + sum += fabs (Get(i)); - return sum; + return sum; } -double Vector :: Max () const + double Vector :: Max () const { - if (!Length()) return 0; - double m = Get(1); - for (INDEX i = 2; i <= Length(); i++) - if (Get(i) > m) m = Get(i); - return m; + if (!Length()) return 0; + double m = Get(1); + for (INDEX i = 2; i <= Length(); i++) + if (Get(i) > m) m = Get(i); + return m; } -double Vector :: Min () const + double Vector :: Min () const { - if (!Length()) return 0; - double m = Get(1); - for (INDEX i = 2; i <= Length(); i++) - if (Get(i) < m) m = Get(i); - return m; + if (!Length()) return 0; + double m = Get(1); + for (INDEX i = 2; i <= Length(); i++) + if (Get(i) < m) m = Get(i); + return m; } -/* -ostream & operator<<(ostream & s, const Vector & v) - { - int w = s.width(); - if (v.Length()) + /* + ostream & operator<<(ostream & s, const Vector & v) + { + int w = s.width(); + if (v.Length()) { s.width(0); s << '('; for (INDEX i = 1; i < v.Length(); i++) - { - s.width(w); - s << v.Get(i) << ","; - if (i % 8 == 0) s << endl << ' '; - } + { + s.width(w); + s << v.Get(i) << ","; + if (i % 8 == 0) s << endl << ' '; + } s.width(w); s << v.Get(v.Length()) << ')'; } - else + else s << "(Vector not allocated)"; - return s; - } -*/ - -ostream & Vector :: Print (ostream & s) const - { - int w = s.width(); - if (Length()) - { - s.width(0); - s << '('; - for (INDEX i = 1; i < Length(); i++) - { - s.width(w); - s << Get(i) << ","; - if (i % 8 == 0) s << endl << ' '; - } - s.width(w); - s << Get(Length()) << ')'; + return s; } - else - s << "(Vector not allocated)"; + */ - return s; + ostream & Vector :: Print (ostream & s) const + { + int w = s.width(); + if (Length()) + { + s.width(0); + s << '('; + for (INDEX i = 1; i < Length(); i++) + { + s.width(w); + s << Get(i) << ","; + if (i % 8 == 0) s << endl << ' '; + } + s.width(w); + s << Get(Length()) << ')'; + } + else + s << "(Vector not allocated)"; + + return s; } -BaseVector & Vector :: operator+= (const BaseVector & v2) + BaseVector & Vector :: operator+= (const BaseVector & v2) { - const Vector & hv2 = v2.CastToVector(); + const Vector & hv2 = v2.CastToVector(); - if (Length() == hv2.Length()) + if (Length() == hv2.Length()) + for (INDEX i = 1; i <= Length(); i++) + Elem (i) += hv2.Get(i); + else + (*myerr) << "operator+= illegal dimension" << endl; + return *this; + } + + BaseVector & Vector :: operator-= (const BaseVector & v2) + { + const Vector & hv2 = v2.CastToVector(); + + if (Length() == hv2.Length()) + for (INDEX i = 1; i <= Length(); i++) + Elem (i) -= hv2.Get(i); + else + (*myerr) << "operator-= illegal dimension" << endl; + return *this; + } + + BaseVector & Vector :: operator*= (double c) + { for (INDEX i = 1; i <= Length(); i++) - Elem (i) += hv2.Get(i); - else - (*myerr) << "operator+= illegal dimension" << endl; - return *this; - } - -BaseVector & Vector :: operator-= (const BaseVector & v2) - { - const Vector & hv2 = v2.CastToVector(); - - if (Length() == hv2.Length()) - for (INDEX i = 1; i <= Length(); i++) - Elem (i) -= hv2.Get(i); - else - (*myerr) << "operator-= illegal dimension" << endl; - return *this; - } - -BaseVector & Vector :: operator*= (double c) - { - for (INDEX i = 1; i <= Length(); i++) - Elem(i) *= c; - return *this; + Elem(i) *= c; + return *this; } -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(); - if (Length() == hv2.Length()) - { - double * p1 = data; - double * p2 = hv2.data; - - for (INDEX i = Length(); i > 0; i--) + if (Length() == hv2.Length()) { - (*p1) += scal * (*p2); - p1++; p2++; + double * p1 = data; + double * p2 = hv2.data; + + for (INDEX i = Length(); i > 0; i--) + { + (*p1) += scal * (*p2); + p1++; p2++; + } } - } - else - (*myerr) << "Vector::Add: illegal dimension" << endl; - return *this; + else + (*myerr) << "Vector::Add: illegal dimension" << endl; + return *this; } -BaseVector & Vector :: Add2 (double scal, const BaseVector & v2, - double scal3, const BaseVector & v3) + BaseVector & Vector :: Add2 (double scal, const BaseVector & v2, + double scal3, const BaseVector & v3) { - const Vector & hv2 = v2.CastToVector(); - const Vector & hv3 = v3.CastToVector(); + const Vector & hv2 = v2.CastToVector(); + const Vector & hv3 = v3.CastToVector(); - if (Length() == hv2.Length()) - { - double * p1 = data; - double * p2 = hv2.data; - double * p3 = hv3.data; - - for (INDEX i = Length(); i > 0; i--) + if (Length() == hv2.Length()) { - (*p1) += (scal * (*p2) + scal3 * (*p3)); - p1++; p2++; p3++; + double * p1 = data; + double * p2 = hv2.data; + double * p3 = hv3.data; + + for (INDEX i = Length(); i > 0; i--) + { + (*p1) += (scal * (*p2) + scal3 * (*p3)); + p1++; p2++; p3++; + } } - } - else - (*myerr) << "Vector::Add: illegal dimension" << endl; - return *this; + else + (*myerr) << "Vector::Add: illegal dimension" << endl; + 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(); - if (Length() == hv2.Length()) - { - double * p1 = data; - double * p2 = hv2.data; - - for (INDEX i = Length(); i > 0; i--) + if (Length() == hv2.Length()) { - (*p1) = scal * (*p2); - p1++; p2++; + double * p1 = data; + double * p2 = hv2.data; + + for (INDEX i = Length(); i > 0; i--) + { + (*p1) = scal * (*p2); + p1++; p2++; + } } - } - else - (*myerr) << "Vector::Set: illegal dimension" << endl; - return *this; + else + (*myerr) << "Vector::Set: illegal dimension" << endl; + return *this; } -BaseVector & Vector :: Set2 (double scal , const BaseVector & v2, - double scal3, const BaseVector & v3) -{ - const Vector & hv2 = v2.CastToVector(); - const Vector & hv3 = v3.CastToVector(); + BaseVector & Vector :: Set2 (double scal , const BaseVector & v2, + double scal3, const BaseVector & v3) + { + const Vector & hv2 = v2.CastToVector(); + const Vector & hv3 = v3.CastToVector(); - if (Length() == hv2.Length() && Length() == hv3.Length()) - { - double * p1 = data; - double * p2 = hv2.data; - double * p3 = hv3.data; + if (Length() == hv2.Length() && Length() == hv3.Length()) + { + double * p1 = data; + double * p2 = hv2.data; + double * p3 = hv3.data; - for (INDEX i = Length(); i > 0; i--) - { - (*p1) = scal * (*p2) + scal3 * (*p3); - p1++; p2++; p3++; - } - } - else - (*myerr) << "Vector::Set: illegal dimension" << endl; - return *this; -} - - -void Vector :: GetPart (int startpos, BaseVector & v2) const -{ - Vector & hv2 = v2.CastToVector(); - - if (Length() >= startpos + v2.Length() - 1) - { - const double * p1 = &Get(startpos); - double * p2 = &hv2.Elem(1); - - memcpy (p2, p1, hv2.Length() * sizeof(double)); - } - else - MyError ("Vector::GetPart: Vector to short"); -} - - -// NEW RM -void Vector :: SetPart (int startpos, const BaseVector & v2) -{ - const Vector & hv2 = v2.CastToVector(); - INDEX i; - INDEX n = v2.Length(); - - if (Length() >= startpos + n - 1) - { - double * p1 = &Elem(startpos); - const double * p2 = &hv2.Get(1); - - for (i = 1; i <= n; i++) - { - (*p1) = (*p2); - p1++; - p2++; - } - } - else - MyError ("Vector::SetPart: Vector to short"); -} - -void Vector :: AddPart (int startpos, double val, const BaseVector & v2) -{ - const Vector & hv2 = v2.CastToVector(); - INDEX i; - INDEX n = v2.Length(); - - if (Length() >= startpos + n - 1) - { - double * p1 = &Elem(startpos); - const double * p2 = &hv2.Get(1); - - for (i = 1; i <= n; i++) - { - (*p1) += val * (*p2); - p1++; - p2++; - } - } - else - MyError ("Vector::AddPart: Vector to short"); -} - - - - -double Vector :: operator* (const BaseVector & v2) const - { - const Vector & hv2 = v2.CastToVector(); - - double sum = 0; - double * p1 = data; - double * p2 = hv2.data; - - if (Length() == hv2.Length()) - { - for (INDEX i = Length(); i > 0; i--) - { - sum += (*p1) * (*p2); - p1++; p2++; + for (INDEX i = Length(); i > 0; i--) + { + (*p1) = scal * (*p2) + scal3 * (*p3); + p1++; p2++; p3++; + } } - } - else - (*myerr) << "Scalarproduct illegal dimension" << endl; - return sum; + else + (*myerr) << "Vector::Set: illegal dimension" << endl; + return *this; } -void Vector :: SetRandom () - { - INDEX i; - for (i = 1; i <= Length(); i++) - Elem(i) = rand (); - double l2 = L2Norm(); - if (l2 > 0) - (*this) /= l2; + void Vector :: GetPart (int startpos, BaseVector & v2) const + { + Vector & hv2 = v2.CastToVector(); + + if (Length() >= startpos + v2.Length() - 1) + { + const double * p1 = &Get(startpos); + double * p2 = &hv2.Elem(1); + + memcpy (p2, p1, hv2.Length() * sizeof(double)); + } + else + MyError ("Vector::GetPart: Vector to short"); + } + + + // NEW RM + void Vector :: SetPart (int startpos, const BaseVector & v2) + { + const Vector & hv2 = v2.CastToVector(); + INDEX i; + INDEX n = v2.Length(); + + if (Length() >= startpos + n - 1) + { + double * p1 = &Elem(startpos); + const double * p2 = &hv2.Get(1); + + for (i = 1; i <= n; i++) + { + (*p1) = (*p2); + p1++; + p2++; + } + } + else + MyError ("Vector::SetPart: Vector to short"); + } + + void Vector :: AddPart (int startpos, double val, const BaseVector & v2) + { + const Vector & hv2 = v2.CastToVector(); + INDEX i; + INDEX n = v2.Length(); + + if (Length() >= startpos + n - 1) + { + double * p1 = &Elem(startpos); + const double * p2 = &hv2.Get(1); + + for (i = 1; i <= n; i++) + { + (*p1) += val * (*p2); + p1++; + p2++; + } + } + else + MyError ("Vector::AddPart: Vector to short"); + } + + + + + double Vector :: operator* (const BaseVector & v2) const + { + const Vector & hv2 = v2.CastToVector(); + + double sum = 0; + double * p1 = data; + double * p2 = hv2.data; + + if (Length() == hv2.Length()) + { + for (INDEX i = Length(); i > 0; i--) + { + sum += (*p1) * (*p2); + p1++; p2++; + } + } + else + (*myerr) << "Scalarproduct illegal dimension" << endl; + return sum; + } + + void Vector :: SetRandom () + { + INDEX i; + for (i = 1; i <= Length(); i++) + Elem(i) = rand (); + + double l2 = L2Norm(); + if (l2 > 0) + (*this) /= l2; // Elem(i) = 1.0 / double(i); // Elem(i) = drand48(); } -/* -TempVector Vector :: operator- () const - { - Vector & sum = *(Vector*)Copy(); + /* + TempVector Vector :: operator- () const + { + Vector & sum = *(Vector*)Copy(); - if (sum.Length () == Length()) + if (sum.Length () == Length()) { for (INDEX i = 1; i <= Length(); i++) - sum.Set (i, Get(i)); + sum.Set (i, Get(i)); } - else + else (*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()); - if (data == v2.data) return *this; + if (data == v2.data) return *this; - if (v2.Length() == Length()) - memcpy (data, v2.data, sizeof (double) * Length()); - else - (*myerr) << "Vector::operator=: not allocated" << endl; + if (v2.Length() == Length()) + memcpy (data, v2.data, sizeof (double) * Length()); + else + (*myerr) << "Vector::operator=: not allocated" << endl; - 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(); - SetLength (hv2.Length()); + SetLength (hv2.Length()); - if (data == hv2.data) return *this; + if (data == hv2.data) return *this; - if (hv2.Length() == Length()) - memcpy (data, hv2.data, sizeof (double) * Length()); - else - (*myerr) << "Vector::operator=: not allocated" << endl; + if (hv2.Length() == Length()) + memcpy (data, hv2.data, sizeof (double) * Length()); + else + (*myerr) << "Vector::operator=: not allocated" << endl; - return *this; + return *this; } -BaseVector & Vector::operator= (double scal) + BaseVector & Vector::operator= (double scal) { - if (!Length()) (*myerr) << "Vector::operator= (double) : data not allocated" - << endl; + if (!Length()) (*myerr) << "Vector::operator= (double) : data not allocated" + << endl; - for (INDEX i = 1; i <= Length(); i++) - Set (i, scal); + for (INDEX i = 1; i <= Length(); i++) + Set (i, scal); - return *this; + return *this; } -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(); - swap (length, hv2.length); - swap (data, hv2.data); + Vector & hv2 = v2.CastToVector(); + swap (length, hv2.length); + swap (data, hv2.data); } -void Vector :: GetElementVector (const Array & pnum, - BaseVector & elvec) const -{ - int i; - Vector & helvec = elvec.CastToVector(); - for (i = 1; i <= pnum.Size(); i++) - helvec.Elem(i) = Get(pnum.Get(i)); -} + void Vector :: GetElementVector (const Array & pnum, + BaseVector & elvec) const + { + int i; + Vector & helvec = elvec.CastToVector(); + for (i = 1; i <= pnum.Size(); i++) + helvec.Elem(i) = Get(pnum.Get(i)); + } -void Vector :: SetElementVector (const Array & pnum, - const BaseVector & elvec) -{ - int i; - const Vector & helvec = elvec.CastToVector(); - for (i = 1; i <= pnum.Size(); i++) - Elem(pnum.Get(i)) = helvec.Get(i); -} + void Vector :: SetElementVector (const Array & pnum, + const BaseVector & elvec) + { + int i; + const Vector & helvec = elvec.CastToVector(); + for (i = 1; i <= pnum.Size(); i++) + Elem(pnum.Get(i)) = helvec.Get(i); + } -void Vector :: AddElementVector (const Array & pnum, - const BaseVector & elvec) -{ - int i; - const Vector & helvec = elvec.CastToVector(); - for (i = 1; i <= pnum.Size(); i++) - Elem(pnum.Get(i)) += helvec.Get(i); -} + void Vector :: AddElementVector (const Array & pnum, + const BaseVector & elvec) + { + int i; + const Vector & helvec = elvec.CastToVector(); + for (i = 1; i <= pnum.Size(); i++) + Elem(pnum.Get(i)) += helvec.Get(i); + } } #endif diff --git a/libsrc/linalg/vector.hpp b/libsrc/linalg/vector.hpp index 90ee657f..fbd38df9 100644 --- a/libsrc/linalg/vector.hpp +++ b/libsrc/linalg/vector.hpp @@ -35,9 +35,9 @@ public: double & operator() (int i) { return data[i]; } const double & operator() (int i) const { return data[i]; } - double & Elem (int i) { return data[i-1]; } - const double & Get (int i) const { return data[i-1]; } - void Set (int i, double val) { data[i-1] = val; } + // double & Elem (int i) { return data[i-1]; } + // const double & Get (int i) const { return data[i-1]; } + // void Set (int i, double val) { data[i-1] = val; } FlatVector & operator*= (double scal) { @@ -92,7 +92,7 @@ public: { delete [] data; } 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) { diff --git a/libsrc/meshing/adfront3.cpp b/libsrc/meshing/adfront3.cpp index 85c4abf8..c5dbb8d4 100644 --- a/libsrc/meshing/adfront3.cpp +++ b/libsrc/meshing/adfront3.cpp @@ -793,25 +793,25 @@ bool AdFront3 :: Inside (const Point<3> & p) const v1 = p2 - p1; v2 = p3 - p1; - a.Elem(1, 1) = v1.X(); - a.Elem(2, 1) = v1.Y(); - a.Elem(3, 1) = v1.Z(); - a.Elem(1, 2) = v2.X(); - a.Elem(2, 2) = v2.Y(); - a.Elem(3, 2) = v2.Z(); - a.Elem(1, 3) = -n.X(); - a.Elem(2, 3) = -n.Y(); - a.Elem(3, 3) = -n.Z(); + a(0, 0) = v1.X(); + a(1, 0) = v1.Y(); + a(2, 0) = v1.Z(); + a(0, 1) = v2.X(); + a(1, 1) = v2.Y(); + a(2, 1) = v2.Z(); + a(0, 2) = -n.X(); + a(1, 2) = -n.Y(); + a(2, 2) = -n.Z(); - b.Elem(1) = p(0) - p1(0); - b.Elem(2) = p(1) - p1(1); - b.Elem(3) = p(2) - p1(2); + b(0) = p(0) - p1(0); + b(1) = p(1) - p1(1); + b(2) = p(2) - p1(2); CalcInverse (a, ainv); ainv.Mult (b, u); - if (u.Elem(1) >= 0 && u.Elem(2) >= 0 && u.Elem(1)+u.Elem(2) <= 1 && - u.Elem(3) > 0) + if (u(0) >= 0 && u(1) >= 0 && u(0)+u(1) <= 1 && + u(2) > 0) { cnt++; } diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 163996eb..f347f1e3 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -1091,26 +1091,23 @@ namespace netgen T_MPRISMS & mprisms, const Mesh & mesh) { - int i, j, k; - int step; - int marked = 0; int np = mesh.GetNP(); Vector hv(np); - for (i = 1; i <= np; i++) - hv.Elem(i) = mesh.GetH (mesh.Point(i)); + for (int i = 0; i < np; i++) + hv(i) = mesh.GetH (mesh.Point(i+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; - for (j = 0; j < 3; j++) - for (k = j+1; k < 4; k++) + for (int j = 0; j < 3; j++) + for (int k = j+1; k < 4; k++) { const Point<3> & p1 = mesh.Point (mtets.Get(i).pnums[j]); const Point<3> & p2 = mesh.Point (mtets.Get(i).pnums[k]); @@ -1120,9 +1117,9 @@ namespace netgen h = sqrt (h); 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) 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; - for (j = 0; j < 2; j++) - for (k = j+1; k < 3; k++) + for (int j = 0; j < 2; j++) + for (int k = j+1; k < 3; k++) { const Point<3> & p1 = mesh.Point (mprisms.Get(i).pnums[j]); const Point<3> & p2 = mesh.Point (mprisms.Get(i).pnums[k]); @@ -1160,9 +1157,9 @@ namespace netgen h = sqrt (h); 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) hshould = hi; } diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 5ea6b239..cc22d1bf 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -117,7 +117,7 @@ namespace netgen } - void HPRefElement :: SetType( HPREF_ELEMENT_TYPE t) + void HPRefElement :: SetType (HPREF_ELEMENT_TYPE t) { type = t; switch(type) @@ -129,11 +129,17 @@ namespace netgen case HP_PRISM: np=6; break; case HP_PYRAMID: np=5; 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; - 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: hpel.type = HP_PYRAMID; break; + + default: + cerr << "HPRefElement: illegal elementtype (1) " << mesh[i].GetType() << endl; + throw NgException ("HPRefElement: illegal elementtype (1)"); + } elements.Append(hpel); } @@ -700,7 +711,7 @@ namespace netgen HPRef_Struct * hprs = Get_HPRef_Struct (el.type); int newlevel = el.levelx + 1; - int oldnp(0); + int oldnp = 0; switch (hprs->geom) { case HP_SEGM: oldnp = 2; break; @@ -710,6 +721,10 @@ namespace netgen case HP_PYRAMID: oldnp = 5; break; case HP_PRISM: oldnp = 6; 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_TET: newel.np=4; break; case HP_PYRAMID: newel.np=5; break; + default: + throw NgException (string("hprefinement.cpp: illegal type")); } for (int k = 0; k < newel.np; k++) @@ -915,6 +932,11 @@ namespace netgen case HP_PRISM: newel.np=6; break; case HP_TET: newel.np=4; 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]; for (int k = 0; k < 8; k++) @@ -1477,6 +1499,12 @@ namespace netgen ord_dir[2] = 2; ned = 8; break; + + + default: + cerr << "HPRefElement: illegal elementtype (2) " << mesh[i].GetType() << endl; + throw NgException ("HPRefElement: illegal elementtype (2)"); + } for (int j=0;j 0.1 * badmax) BFGS (px, pf, par); bad2 = pf.Func (px); - pnew.X() = px.Get(1); - pnew.Y() = px.Get(2); - pnew.Z() = px.Get(3); + pnew.X() = px(0); + pnew.Y() = px(1); + pnew.Z() = px(2); int hpinew = mesh.AddPoint (pnew); @@ -584,11 +584,6 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, 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); int cnt = 0; @@ -600,7 +595,6 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, int np = mesh.GetNP(); int ne = mesh.GetNE(); - //int nse = mesh.GetNSE(); // contains at least all elements at node TABLE elementsonnode(np); @@ -614,16 +608,6 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, multithread.task = "Swap Improve"; // 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 faces(mesh.GetNOpenElements()/3 + 2); if (goal == OPT_CONFORM) @@ -645,32 +629,15 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, } // find elements on node - for (ei = 0; ei < ne; ei++) - for (j = 0; j < mesh[ei].GetNP(); j++) + for (ElementIndex ei = 0; ei < ne; ei++) + for (int j = 0; j < mesh[ei].GetNP(); j++) 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 edgeused(2 * ne + 5); - - for (ei = 0; ei < ne; ei++) + // INDEX_2_HASHTABLE edgeused(2 * ne + 5); + INDEX_2_CLOSED_HASHTABLE edgeused(12 * ne + 5); + + for (ElementIndex ei = 0; ei < ne; ei++) { if (multithread.terminate) break; @@ -695,7 +662,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, // int onlybedges = 1; - for (j = 0; j < 6; j++) + for (int j = 0; j < 6; j++) { // loop over edges @@ -726,7 +693,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, edgeused.Set (i2, 1); 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; ElementIndex elnr = elementsonnode[pi1][k]; @@ -734,7 +701,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, 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] == pi2) has2 = 1; @@ -742,7 +709,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, if (has1 && has2) { // only once - for (l = 0; l < hasbothpoints.Size(); l++) + for (int l = 0; l < hasbothpoints.Size(); l++) if (hasbothpoints[l] == elnr) has1 = 0; @@ -752,7 +719,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, } bool puretet = 1; - for (k = 0; k < hasbothpoints.Size(); k++) + for (int k = 0; k < hasbothpoints.Size(); k++) if (mesh[hasbothpoints[k]].GetType () != TET) puretet = 0; if (!puretet) @@ -763,7 +730,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, if ( nsuround == 3 ) { 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) { pi4 = pi3; @@ -784,16 +751,16 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, } pi5 = 0; - for (k = 1; k < 3; k++) + for (int k = 1; k < 3; k++) { const Element & elemk = mesh[hasbothpoints[k]]; bool has1 = 0; - for (l = 0; l < 4; l++) + for (int l = 0; l < 4; l++) if (elemk[l] == pi4) has1 = 1; 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) pi5 = elemk[l]; } @@ -905,21 +872,20 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, el22.flags.illegal_valid = 0; mesh[hasbothpoints[0]] = el21; 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]].Delete(); - for (k = 0; k < 2; k++) - for (l = 0; l < 4; l++) + for (int k = 0; k < 2; k++) + for (int l = 0; l < 4; l++) elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]); } - } - - + } + if (nsuround == 4) { 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) { pi4 = pi3; @@ -938,32 +904,32 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, } pi5 = 0; - for (k = 1; k < 4; k++) + for (int k = 1; k < 4; k++) { const Element & elem = mesh[hasbothpoints[k]]; bool has1 = 0; - for (l = 0; l < 4; l++) + for (int l = 0; l < 4; l++) if (elem[l] == pi4) has1 = 1; 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) pi5 = elem[l]; } } pi6 = 0; - for (k = 1; k < 4; k++) + for (int k = 1; k < 4; k++) { const Element & elem = mesh[hasbothpoints[k]]; bool has1 = 0; - for (l = 0; l < 4; l++) + for (int l = 0; l < 4; l++) if (elem[l] == pi3) has1 = 1; 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) pi6 = elem[l]; } @@ -1183,8 +1149,8 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, mesh[hasbothpoints[2]] = el3; mesh[hasbothpoints[3]] = el4; - for (k = 0; k < 4; k++) - for (l = 0; l < 4; l++) + for (int k = 0; k < 4; k++) + for (int l = 0; l < 4; l++) elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]); } else if (swap3) @@ -1216,8 +1182,8 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, mesh[hasbothpoints[3]] = el4b; - for (k = 0; k < 4; k++) - for (l = 0; l < 4; l++) + for (int k = 0; k < 4; k++) + for (int l = 0; l < 4; l++) elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]); } } @@ -1231,7 +1197,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, 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) { pi4 = pi3; @@ -1259,12 +1225,12 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, tetused = 0; tetused[0] = 1; - for (l = 2; l < nsuround; l++) + for (int l = 2; l < nsuround; l++) { int oldpi = suroundpts[l-1]; int newpi = 0; - for (k = 0; k < nsuround && !newpi; k++) + for (int k = 0; k < nsuround && !newpi; k++) if (!tetused[k]) { const Element & nel = mesh[hasbothpoints[k]]; @@ -1284,7 +1250,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, bad1 = 0; - for (k = 0; k < nsuround; k++) + for (int k = 0; k < nsuround; k++) { hel[0] = pi1; hel[1] = pi2; @@ -1303,11 +1269,11 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, int confedge = -1; double badopt = bad1; - for (l = 0; l < nsuround; l++) + for (int l = 0; l < nsuround; l++) { 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[1] = suroundpts[k % nsuround]; @@ -1344,7 +1310,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, (bad2 <= 100 * bad1 && bad2 <= 1e18) || (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], 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], suroundpts[k % nsuround])) @@ -1385,7 +1351,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, // (*mycout) << nsuround << "->" << 2 * (nsuround-2) << " " << flush; cnt++; - for (k = bestl+1; k <= nsuround + bestl - 2; k++) + for (int k = bestl+1; k <= nsuround + bestl - 2; k++) { int k1; @@ -1420,7 +1386,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, 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]]; /* @@ -1513,8 +1479,6 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal, } } - ElementIndex ei; - SurfaceElementIndex sei; PointIndex pi1, pi2, pi3, pi4, pi5, pi6; PointIndex pi1other, pi2other; @@ -1548,11 +1512,11 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal, // 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++) elementsonnode.Add (mesh[ei][j], ei); - for (sei = 0; sei < nse; sei++) + for (SurfaceElementIndex sei = 0; sei < nse; sei++) for(int j=0; j edgeused(2 * ne + 5); + // INDEX_2_HASHTABLE edgeused(2 * ne + 5); + INDEX_2_CLOSED_HASHTABLE edgeused(12 * ne + 5); - for (ei = 0; ei < ne; ei++) + for (ElementIndex ei = 0; ei < ne; ei++) { if (multithread.terminate) break; @@ -2354,18 +2319,12 @@ void MeshOptimize3d :: SwapImproveSurface (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); - - int cnt = 0; - Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); + int cnt = 0; double bad1, bad2; int np = mesh.GetNP(); @@ -2383,7 +2342,6 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) // TestOk(); - /* CalcSurfacesOfNode (); for (i = 1; i <= GetNE(); i++) @@ -2404,15 +2362,15 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) // find elements on node - for (ei = 0; ei < ne; ei++) - for (j = 0; j < mesh[ei].GetNP(); j++) + for (ElementIndex ei = 0; ei < ne; ei++) + for (int j = 0; j < mesh[ei].GetNP(); j++) elementsonnode.Add (mesh[ei][j], ei); - for (sei = 0; sei < nse; sei++) - for (j = 0; j < 3; j++) + for (SurfaceElementIndex sei = 0; sei < nse; sei++) + for (int j = 0; j < 3; j++) belementsonnode.Add (mesh[sei][j], sei); - for (eli1 = 0; eli1 < ne; eli1++) + for (ElementIndex eli1 = 0; eli1 < ne; eli1++) { if (multithread.terminate) break; @@ -2431,7 +2389,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) // cout << "eli = " << 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 @@ -2463,13 +2421,13 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) bool bface = 0; - for (k = 0; k < belementsonnode[pi1].Size(); k++) + for (int k = 0; k < belementsonnode[pi1].Size(); k++) { const Element2d & bel = mesh[belementsonnode[pi1][k]]; 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) { bface1 = 0; @@ -2487,9 +2445,9 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) FlatArray 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 // << ", getne = " << mesh.GetNE(); @@ -2502,7 +2460,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) continue; 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 || elem2.PNum(l) == pi3) { @@ -2586,7 +2544,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) } */ // cout << "neli = " << neli << endl; - for (l = 0; l < 4; l++) + for (int l = 0; l < 4; l++) { elementsonnode.Add (el31[l], eli1); elementsonnode.Add (el32[l], eli2); diff --git a/libsrc/meshing/improve3.hpp b/libsrc/meshing/improve3.hpp index 395fd5a0..100ef54a 100644 --- a/libsrc/meshing/improve3.hpp +++ b/libsrc/meshing/improve3.hpp @@ -20,8 +20,18 @@ public: -extern double CalcBad (const Mesh::T_POINTS & points, const Element & elem, - double h); +inline double +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, const Mesh::T_VOLELEMENTS & elements); diff --git a/libsrc/meshing/meshtool.cpp b/libsrc/meshing/meshtool.cpp index 9ebc560e..b69e9c26 100644 --- a/libsrc/meshing/meshtool.cpp +++ b/libsrc/meshing/meshtool.cpp @@ -202,7 +202,7 @@ namespace netgen Vec3d v2 (p1, p3); Vec3d v3 (p1, p4); - vol = -Determinant (v1, v2, v3) / 6; + vol = Determinant (v1, v2, v3) * (-0.166666666666666); ll1 = v1.Length2(); ll2 = v2.Length2(); @@ -276,7 +276,7 @@ namespace netgen Vec3d v5 (*pp2, *pp4); Vec3d v6 (*pp3, *pp4); - vol = -Determinant (v1, v2, v3) / 6; + vol = Determinant (v1, v2, v3) * (-0.166666666666666); Vec3d 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/(ll2*ll2)) * gradll2; graderr += (1/(h*h) - h*h/(ll3*ll3)) * gradll3; - cout << "?"; } - double errpow; if (teterrpow == 2) { diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index ac910f32..7c95d957 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -12,13 +12,6 @@ namespace netgen return s; } - /* - MultiPointGeomInfo :: MultiPointGeomInfo() - { - cnt = 0; - } - */ - int MultiPointGeomInfo :: AddPointGeomInfo (const PointGeomInfo & gi) { @@ -37,22 +30,6 @@ namespace netgen } - /* - void MultiPointGeomInfo :: - Init () - { - cnt = 0; - } - - void MultiPointGeomInfo :: - DeleteAll () - { - cnt = 0; - } - */ - - - Segment :: Segment() { pnums[0] = -1; @@ -74,13 +51,13 @@ namespace netgen pnums[2] = -1; meshdocval = 0; /* - geominfo[0].trignum=-1; - geominfo[1].trignum=-1; + geominfo[0].trignum=-1; + geominfo[1].trignum=-1; - epgeominfo[0].edgenr = 1; - epgeominfo[0].dist = 0; - epgeominfo[1].edgenr = 1; - epgeominfo[1].dist = 0; + epgeominfo[0].edgenr = 1; + epgeominfo[0].dist = 0; + epgeominfo[1].edgenr = 1; + epgeominfo[1].dist = 0; */ bcname = 0; @@ -88,20 +65,20 @@ namespace netgen Segment::Segment (const Segment & other) : - edgenr(other.edgenr), - singedge_left(other.singedge_left), - singedge_right(other.singedge_right), - seginfo(other.seginfo), - si(other.si), - domin(other.domin), - domout(other.domout), - tlosurf(other.tlosurf), - geominfo(), - surfnr1(other.surfnr1), - surfnr2(other.surfnr2), - epgeominfo(), - meshdocval(other.meshdocval), - hp_elnr(other.hp_elnr) + edgenr(other.edgenr), + singedge_left(other.singedge_left), + singedge_right(other.singedge_right), + seginfo(other.seginfo), + si(other.si), + domin(other.domin), + domout(other.domout), + tlosurf(other.tlosurf), + geominfo(), + surfnr1(other.surfnr1), + surfnr2(other.surfnr2), + epgeominfo(), + meshdocval(other.meshdocval), + hp_elnr(other.hp_elnr) { for (int j = 0; j < 3; j++) pnums[j] = other.pnums[j]; @@ -217,7 +194,7 @@ namespace netgen refflag = 1; strongrefflag = false; #ifdef PARALLEL - isghost = 0; + isghost = 0; #endif } @@ -225,686 +202,686 @@ namespace netgen Element2d :: Element2d (int pi1, int pi2, int pi3) -{ - pnum[0] = pi1; - pnum[1] = pi2; - pnum[2] = pi3; - np = 3; - typ = TRIG; - pnum[3] = 0; - pnum[4] = 0; - pnum[5] = 0; + { + pnum[0] = pi1; + pnum[1] = pi2; + pnum[2] = pi3; + np = 3; + typ = TRIG; + pnum[3] = 0; + pnum[4] = 0; + pnum[5] = 0; - for (int i = 0; i < ELEMENT2D_MAXPOINTS; i++) - geominfo[i].trignum = 0; - index = 0; - badel = 0; - refflag = 1; - strongrefflag = false; - deleted = 0; - orderx = ordery = 1; + for (int i = 0; i < ELEMENT2D_MAXPOINTS; i++) + geominfo[i].trignum = 0; + index = 0; + badel = 0; + refflag = 1; + strongrefflag = false; + deleted = 0; + orderx = ordery = 1; #ifdef PARALLEL - isghost = 0; + isghost = 0; #endif -} + } -Element2d :: Element2d (int pi1, int pi2, int pi3, int pi4) -{ - pnum[0] = pi1; - pnum[1] = pi2; - pnum[2] = pi3; - pnum[3] = pi4; - np = 4; - typ = QUAD; + Element2d :: Element2d (int pi1, int pi2, int pi3, int pi4) + { + pnum[0] = pi1; + pnum[1] = pi2; + pnum[2] = pi3; + pnum[3] = pi4; + np = 4; + typ = QUAD; - pnum[4] = 0; - pnum[5] = 0; + pnum[4] = 0; + pnum[5] = 0; - for (int i = 0; i < ELEMENT2D_MAXPOINTS; i++) - geominfo[i].trignum = 0; - index = 0; - badel = 0; - refflag = 1; - strongrefflag = false; - deleted = 0; - orderx = ordery = 1; + for (int i = 0; i < ELEMENT2D_MAXPOINTS; i++) + geominfo[i].trignum = 0; + index = 0; + badel = 0; + refflag = 1; + strongrefflag = false; + deleted = 0; + orderx = ordery = 1; #ifdef PARALLEL - isghost = 0; + isghost = 0; #endif -} + } -/* -void Element2d :: SetType (ELEMENT_TYPE atyp) -{ - typ = atyp; - switch (typ) + /* + void Element2d :: SetType (ELEMENT_TYPE atyp) + { + typ = atyp; + switch (typ) { case TRIG: np = 3; break; case QUAD: np = 4; break; case TRIG6: np = 6; break; case QUAD6: np = 6; break; default: - PrintSysError ("Element2d::SetType, illegal type ", typ); + PrintSysError ("Element2d::SetType, illegal type ", typ); } -} -*/ - - -void Element2d :: GetBox (const T_POINTS & points, Box3d & box) const -{ - box.SetPoint (points.Get(pnum[0])); - for (unsigned i = 1; i < np; i++) - box.AddPoint (points.Get(pnum[i])); -} - -bool Element2d :: operator==(const Element2d & el2) const -{ - bool retval = (el2.GetNP() == np); - for(int i= 0; retval && i ipdtrig; -Array ipdquad; - - -int Element2d :: GetNIP () const -{ - int nip; - switch (np) - { - case 3: nip = 1; break; - case 4: nip = 4; break; - default: nip = 0; break; - } - return nip; -} - -void Element2d :: -GetIntegrationPoint (int ip, Point2d & p, double & weight) const -{ - static double eltriqp[1][3] = - { - { 1.0/3.0, 1.0/3.0, 0.5 } - }; - - static double elquadqp[4][3] = - { - { 0, 0, 0.25 }, - { 0, 1, 0.25 }, - { 1, 0, 0.25 }, - { 1, 1, 0.25 } - }; - - double * pp = 0; - switch (typ) - { - case TRIG: pp = &eltriqp[0][0]; break; - case QUAD: pp = &elquadqp[ip-1][0]; break; - default: - PrintSysError ("Element2d::GetIntegrationPoint, illegal type ", typ); - } - - p.X() = pp[0]; - p.Y() = pp[1]; - weight = pp[2]; -} - -void Element2d :: -GetTransformation (int ip, const Array & points, - DenseMatrix & trans) const -{ - int np = GetNP(); - static DenseMatrix pmat(2, np), dshape(2, np); - pmat.SetSize (2, np); - dshape.SetSize (2, np); - - Point2d p; - double w; - - GetPointMatrix (points, pmat); - GetIntegrationPoint (ip, p, w); - GetDShape (p, dshape); - - CalcABt (pmat, dshape, trans); - - /* - (*testout) << "p = " << p << endl - << "pmat = " << pmat << endl - << "dshape = " << dshape << endl - << "tans = " << trans << endl; */ -} - -void Element2d :: -GetTransformation (int ip, class DenseMatrix & pmat, - class DenseMatrix & trans) const -{ - int np = GetNP(); - -#ifdef DEBUG - if (pmat.Width() != np || pmat.Height() != 2) - { - (*testout) << "GetTransofrmation: pmat doesn't fit" << endl; - return; - } -#endif - - ComputeIntegrationPointData (); - DenseMatrix * dshapep; - switch (typ) - { - case TRIG: dshapep = &ipdtrig.Get(ip)->dshape; break; - case QUAD: dshapep = &ipdquad.Get(ip)->dshape; break; - default: - PrintSysError ("Element2d::GetTransformation, illegal type ", typ); - } - - CalcABt (pmat, *dshapep, trans); -} + void Element2d :: GetBox (const T_POINTS & points, Box3d & box) const + { + box.SetPoint (points.Get(pnum[0])); + for (unsigned i = 1; i < np; i++) + box.AddPoint (points.Get(pnum[i])); + } -void Element2d :: GetShape (const Point2d & p, Vector & shape) const -{ - if (shape.Size() != GetNP()) - { - cerr << "Element::GetShape: Length not fitting" << endl; - return; - } + bool Element2d :: operator==(const Element2d & el2) const + { + bool retval = (el2.GetNP() == np); + for(int i= 0; retval && i & p, FlatVector & shape) const -{ - switch (typ) - { - case TRIG: + void Element2d :: Invert2() + { + switch (typ) { - shape(0) = p(0); - shape(1) = p(1); - shape(2) = 1-p(0)-p(1); - break; + case TRIG: + { + Swap (pnum[1], pnum[2]); + break; + } + case QUAD: + { + Swap (pnum[0], pnum[3]); + Swap (pnum[1], pnum[2]); + break; + } + default: + { + cerr << "Element2d::Invert2, illegal element type " << int(typ) << endl; + } } + } - case QUAD: + int Element2d::HasFace(const Element2d & el) const + { + //nur für tets!!! hannes + for (int i = 1; i <= 3; i++) { - shape(0) = (1-p(0))*(1-p(1)); - shape(1) = p(0) *(1-p(1)); - shape(2) = p(0) * p(1) ; - shape(3) = (1-p(0))* p(1) ; - break; + if (PNumMod(i) == el[0] && + PNumMod(i+1) == el[1] && + PNumMod(i+2) == el[2]) + { + return 1; + } } - } -} + return 0; + } - - - - - - - - -void Element2d :: -GetDShape (const Point2d & p, DenseMatrix & dshape) const -{ -#ifdef DEBUG - if (dshape.Height() != 2 || dshape.Width() != np) - { - PrintSysError ("Element::DShape: Sizes don't fit"); - return; - } -#endif - - switch (typ) - { - case TRIG: - dshape.Elem(1, 1) = -1; - dshape.Elem(1, 2) = 1; - dshape.Elem(1, 3) = 0; - dshape.Elem(2, 1) = -1; - dshape.Elem(2, 2) = 0; - dshape.Elem(2, 3) = 1; - break; - case QUAD: - dshape.Elem(1, 1) = -(1-p.Y()); - dshape.Elem(1, 2) = (1-p.Y()); - dshape.Elem(1, 3) = p.Y(); - dshape.Elem(1, 4) = -p.Y(); - dshape.Elem(2, 1) = -(1-p.X()); - dshape.Elem(2, 2) = -p.X(); - dshape.Elem(2, 3) = p.X(); - dshape.Elem(2, 4) = (1-p.X()); - break; - - default: - PrintSysError ("Element2d::GetDShape, illegal type ", typ); - } -} - - - - -void Element2d :: -GetDShapeNew (const Point<2> & p, MatrixFixWidth<2> & dshape) const -{ - switch (typ) - { - case TRIG: + void Element2d :: NormalizeNumbering2 () + { + if (GetNP() == 3) { - dshape = 0; - dshape(0,0) = 1; - dshape(1,1) = 1; - dshape(2,0) = -1; - dshape(2,1) = -1; - break; + if (PNum(1) < PNum(2) && PNum(1) < PNum(3)) + return; + else + { + if (PNum(2) < PNum(3)) + { + PointIndex pi1 = PNum(2); + PNum(2) = PNum(3); + PNum(3) = PNum(1); + PNum(1) = pi1; + } + else + { + PointIndex pi1 = PNum(3); + PNum(3) = PNum(2); + PNum(2) = PNum(1); + PNum(1) = pi1; + } + } } - case QUAD: + else { - dshape(0,0) = -(1-p(1)); - dshape(0,1) = -(1-p(0)); - - dshape(1,0) = (1-p(1)); - dshape(1,1) = -p(0); - - dshape(2,0) = p(1); - dshape(2,1) = p(0); - - dshape(3,0) = -p(1); - dshape(3,1) = (1-p(0)); - break; - } - } -} - - - - - -void Element2d :: -GetPointMatrix (const Array & points, - DenseMatrix & pmat) const -{ - int np = GetNP(); - -#ifdef DEBUG - if (pmat.Width() != np || pmat.Height() != 2) - { - cerr << "Element::GetPointMatrix: sizes don't fit" << endl; - return; - } -#endif - - for (int i = 1; i <= np; i++) - { - const Point2d & p = points.Get(PNum(i)); - pmat.Elem(1, i) = p.X(); - pmat.Elem(2, i) = p.Y(); - } -} - - - - - -double Element2d :: CalcJacobianBadness (const Array & points) const -{ - int i, j; - int nip = GetNIP(); - static DenseMatrix trans(2,2); - static DenseMatrix pmat; - - pmat.SetSize (2, GetNP()); - GetPointMatrix (points, pmat); - - double err = 0; - for (i = 1; i <= nip; i++) - { - GetTransformation (i, pmat, trans); - - // Frobenius norm - double frob = 0; - for (j = 1; j <= 4; j++) - frob += sqr (trans.Get(j)); - frob = sqrt (frob); - frob /= 2; - - double det = trans.Det(); - - if (det <= 0) - err += 1e12; - else - err += frob * frob / det; - } - - err /= nip; - return err; -} - - - -static const int qip_table[4][4] = - { { 0, 1, 0, 3 }, - { 0, 1, 1, 2 }, - { 3, 2, 0, 3 }, - { 3, 2, 1, 2 } - }; - -double Element2d :: -CalcJacobianBadnessDirDeriv (const Array & points, - int pi, Vec2d & dir, double & dd) const -{ - if (typ == QUAD) - { - Mat<2,2> trans, dtrans; - Mat<2,4> vmat, pmat; + int mini = 1; + for (int i = 2; i <= GetNP(); i++) + if (PNum(i) < PNum(mini)) mini = i; - for (int j = 0; j < 4; j++) - { - const Point2d & p = points.Get( (*this)[j] ); - pmat(0, j) = p.X(); - pmat(1, j) = p.Y(); - } + Element2d hel = (*this); + for (int i = 1; i <= GetNP(); i++) + PNum(i) = hel.PNumMod (i+mini-1); + } + } - vmat = 0.0; - vmat(0, pi-1) = dir.X(); - vmat(1, pi-1) = dir.Y(); + + + + Array ipdtrig; + Array ipdquad; + + + int Element2d :: GetNIP () const + { + int nip; + switch (np) + { + case 3: nip = 1; break; + case 4: nip = 4; break; + default: nip = 0; break; + } + return nip; + } + + void Element2d :: + GetIntegrationPoint (int ip, Point2d & p, double & weight) const + { + static double eltriqp[1][3] = + { + { 1.0/3.0, 1.0/3.0, 0.5 } + }; + + static double elquadqp[4][3] = + { + { 0, 0, 0.25 }, + { 0, 1, 0.25 }, + { 1, 0, 0.25 }, + { 1, 1, 0.25 } + }; + + double * pp = 0; + switch (typ) + { + case TRIG: pp = &eltriqp[0][0]; break; + case QUAD: pp = &elquadqp[ip-1][0]; break; + default: + PrintSysError ("Element2d::GetIntegrationPoint, illegal type ", typ); + } + + p.X() = pp[0]; + p.Y() = pp[1]; + weight = pp[2]; + } + + void Element2d :: + GetTransformation (int ip, const Array & points, + DenseMatrix & trans) const + { + int np = GetNP(); + static DenseMatrix pmat(2, np), dshape(2, np); + pmat.SetSize (2, np); + dshape.SetSize (2, np); + + Point2d p; + double w; + + GetPointMatrix (points, pmat); + GetIntegrationPoint (ip, p, w); + GetDShape (p, dshape); + + CalcABt (pmat, dshape, trans); + + /* + (*testout) << "p = " << p << endl + << "pmat = " << pmat << endl + << "dshape = " << dshape << endl + << "tans = " << trans << endl; + */ + } + + void Element2d :: + GetTransformation (int ip, class DenseMatrix & pmat, + class DenseMatrix & trans) const + { + // int np = GetNP(); + +#ifdef DEBUG + if (pmat.Width() != np || pmat.Height() != 2) + { + (*testout) << "GetTransofrmation: pmat doesn't fit" << endl; + return; + } +#endif + + ComputeIntegrationPointData (); + DenseMatrix * dshapep = NULL; + switch (typ) + { + case TRIG: dshapep = &ipdtrig.Get(ip)->dshape; break; + case QUAD: dshapep = &ipdquad.Get(ip)->dshape; break; + default: + PrintSysError ("Element2d::GetTransformation, illegal type ", typ); + } + + CalcABt (pmat, *dshapep, trans); + } + + + + void Element2d :: GetShape (const Point2d & p, Vector & shape) const + { + if (shape.Size() != GetNP()) + { + cerr << "Element::GetShape: Length not fitting" << endl; + return; + } + + switch (typ) + { + case TRIG: + shape(0) = 1 - p.X() - p.Y(); + shape(1) = p.X(); + shape(2) = p.Y(); + break; + case QUAD: + shape(0) = (1-p.X()) * (1-p.Y()); + shape(1) = p.X() * (1-p.Y()); + shape(2) = p.X() * p.Y(); + shape(3) = (1-p.X()) * p.Y(); + break; + default: + PrintSysError ("Element2d::GetShape, illegal type ", typ); + } + } + + + + void Element2d :: GetShapeNew (const Point<2> & p, FlatVector & shape) const + { + switch (typ) + { + case TRIG: + { + shape(0) = p(0); + shape(1) = p(1); + shape(2) = 1-p(0)-p(1); + break; + } + + case QUAD: + { + shape(0) = (1-p(0))*(1-p(1)); + shape(1) = p(0) *(1-p(1)); + shape(2) = p(0) * p(1) ; + shape(3) = (1-p(0))* p(1) ; + break; + } + } + } + + + + + + + + + + void Element2d :: + GetDShape (const Point2d & p, DenseMatrix & dshape) const + { +#ifdef DEBUG + if (dshape.Height() != 2 || dshape.Width() != np) + { + PrintSysError ("Element::DShape: Sizes don't fit"); + return; + } +#endif + + switch (typ) + { + case TRIG: + dshape.Elem(1, 1) = -1; + dshape.Elem(1, 2) = 1; + dshape.Elem(1, 3) = 0; + dshape.Elem(2, 1) = -1; + dshape.Elem(2, 2) = 0; + dshape.Elem(2, 3) = 1; + break; + case QUAD: + dshape.Elem(1, 1) = -(1-p.Y()); + dshape.Elem(1, 2) = (1-p.Y()); + dshape.Elem(1, 3) = p.Y(); + dshape.Elem(1, 4) = -p.Y(); + dshape.Elem(2, 1) = -(1-p.X()); + dshape.Elem(2, 2) = -p.X(); + dshape.Elem(2, 3) = p.X(); + dshape.Elem(2, 4) = (1-p.X()); + break; + + default: + PrintSysError ("Element2d::GetDShape, illegal type ", typ); + } + } + + + + + void Element2d :: + GetDShapeNew (const Point<2> & p, MatrixFixWidth<2> & dshape) const + { + switch (typ) + { + case TRIG: + { + dshape = 0; + dshape(0,0) = 1; + dshape(1,1) = 1; + dshape(2,0) = -1; + dshape(2,1) = -1; + break; + } + case QUAD: + { + dshape(0,0) = -(1-p(1)); + dshape(0,1) = -(1-p(0)); + + dshape(1,0) = (1-p(1)); + dshape(1,1) = -p(0); + + dshape(2,0) = p(1); + dshape(2,1) = p(0); + + dshape(3,0) = -p(1); + dshape(3,1) = (1-p(0)); + break; + } + } + } + + + + + + void Element2d :: + GetPointMatrix (const Array & points, + DenseMatrix & pmat) const + { + int np = GetNP(); + +#ifdef DEBUG + if (pmat.Width() != np || pmat.Height() != 2) + { + cerr << "Element::GetPointMatrix: sizes don't fit" << endl; + return; + } +#endif + + for (int i = 1; i <= np; i++) + { + const Point2d & p = points.Get(PNum(i)); + pmat.Elem(1, i) = p.X(); + pmat.Elem(2, i) = p.Y(); + } + } + + + + + + double Element2d :: CalcJacobianBadness (const Array & points) const + { + int i, j; + int nip = GetNIP(); + static DenseMatrix trans(2,2); + static DenseMatrix pmat; + + pmat.SetSize (2, GetNP()); + GetPointMatrix (points, pmat); + + double err = 0; + for (i = 1; i <= nip; i++) + { + GetTransformation (i, pmat, trans); + + // Frobenius norm + double frob = 0; + for (j = 1; j <= 4; j++) + frob += sqr (trans.Get(j)); + frob = sqrt (frob); + frob /= 2; + + double det = trans.Det(); + + if (det <= 0) + err += 1e12; + else + err += frob * frob / det; + } + + err /= nip; + return err; + } + + + + static const int qip_table[4][4] = + { { 0, 1, 0, 3 }, + { 0, 1, 1, 2 }, + { 3, 2, 0, 3 }, + { 3, 2, 1, 2 } + }; + + double Element2d :: + CalcJacobianBadnessDirDeriv (const Array & points, + int pi, Vec2d & dir, double & dd) const + { + if (typ == QUAD) + { + Mat<2,2> trans, dtrans; + Mat<2,4> vmat, pmat; - double err = 0; - dd = 0; + for (int j = 0; j < 4; j++) + { + const Point2d & p = points.Get( (*this)[j] ); + pmat(0, j) = p.X(); + pmat(1, j) = p.Y(); + } - for (int i = 0; i < 4; i++) - { - int ix1 = qip_table[i][0]; - int ix2 = qip_table[i][1]; - int iy1 = qip_table[i][2]; - int iy2 = qip_table[i][3]; + vmat = 0.0; + vmat(0, pi-1) = dir.X(); + vmat(1, pi-1) = dir.Y(); + + double err = 0; + dd = 0; + + for (int i = 0; i < 4; i++) + { + int ix1 = qip_table[i][0]; + int ix2 = qip_table[i][1]; + int iy1 = qip_table[i][2]; + int iy2 = qip_table[i][3]; - trans(0,0) = pmat(0, ix2) - pmat(0,ix1); - trans(1,0) = pmat(1, ix2) - pmat(1,ix1); - trans(0,1) = pmat(0, iy2) - pmat(0,iy1); - trans(1,1) = pmat(1, iy2) - pmat(1,iy1); + trans(0,0) = pmat(0, ix2) - pmat(0,ix1); + trans(1,0) = pmat(1, ix2) - pmat(1,ix1); + trans(0,1) = pmat(0, iy2) - pmat(0,iy1); + trans(1,1) = pmat(1, iy2) - pmat(1,iy1); - double det = trans(0,0)*trans(1,1)-trans(1,0)*trans(0,1); + double det = trans(0,0)*trans(1,1)-trans(1,0)*trans(0,1); - if (det <= 0) - { - dd = 0; - return 1e12; - } + if (det <= 0) + { + dd = 0; + return 1e12; + } - dtrans(0,0) = vmat(0, ix2) - vmat(0,ix1); - dtrans(1,0) = vmat(1, ix2) - vmat(1,ix1); - dtrans(0,1) = vmat(0, iy2) - vmat(0,iy1); - dtrans(1,1) = vmat(1, iy2) - vmat(1,iy1); + dtrans(0,0) = vmat(0, ix2) - vmat(0,ix1); + dtrans(1,0) = vmat(1, ix2) - vmat(1,ix1); + dtrans(0,1) = vmat(0, iy2) - vmat(0,iy1); + dtrans(1,1) = vmat(1, iy2) - vmat(1,iy1); - // Frobenius norm - double frob = 0; - for (int j = 0; j < 4; j++) - frob += sqr (trans(j)); - frob = sqrt (frob); + // Frobenius norm + double frob = 0; + for (int j = 0; j < 4; j++) + frob += sqr (trans(j)); + frob = sqrt (frob); - double dfrob = 0; - for (int j = 0; j < 4; j++) - dfrob += trans(j) * dtrans(j); - dfrob = dfrob / frob; + double dfrob = 0; + for (int j = 0; j < 4; j++) + dfrob += trans(j) * dtrans(j); + dfrob = dfrob / frob; - frob /= 2; - dfrob /= 2; + frob /= 2; + dfrob /= 2; - // ddet = \sum_j det (m_j) with m_j = trans, except col j = dtrans - double ddet - = dtrans(0,0) * trans(1,1) - trans(0,1) * dtrans(1,0) - + trans(0,0) * dtrans(1,1) - dtrans(0,1) * trans(1,0); + // ddet = \sum_j det (m_j) with m_j = trans, except col j = dtrans + double ddet + = dtrans(0,0) * trans(1,1) - trans(0,1) * dtrans(1,0) + + trans(0,0) * dtrans(1,1) - dtrans(0,1) * trans(1,0); - err += frob * frob / det; - dd += (2 * frob * dfrob * det - frob * frob * ddet) / (det * det); - } + err += frob * frob / det; + dd += (2 * frob * dfrob * det - frob * frob * ddet) / (det * det); + } - err /= 4; - dd /= 4; - return err; - } + err /= 4; + dd /= 4; + return err; + } - int nip = GetNIP(); - static DenseMatrix trans(2,2), dtrans(2,2); - static DenseMatrix pmat, vmat; + int nip = GetNIP(); + static DenseMatrix trans(2,2), dtrans(2,2); + static DenseMatrix pmat, vmat; - pmat.SetSize (2, GetNP()); - vmat.SetSize (2, GetNP()); + pmat.SetSize (2, GetNP()); + vmat.SetSize (2, GetNP()); - GetPointMatrix (points, pmat); + GetPointMatrix (points, pmat); - vmat = 0.0; - vmat.Elem(1, pi) = dir.X(); - vmat.Elem(2, pi) = dir.Y(); + vmat = 0.0; + vmat.Elem(1, pi) = dir.X(); + vmat.Elem(2, pi) = dir.Y(); - double err = 0; - dd = 0; + double err = 0; + dd = 0; - for (int i = 1; i <= nip; i++) - { - GetTransformation (i, pmat, trans); - GetTransformation (i, vmat, dtrans); + for (int i = 1; i <= nip; i++) + { + GetTransformation (i, pmat, trans); + GetTransformation (i, vmat, dtrans); - // Frobenius norm - double frob = 0; - for (int j = 1; j <= 4; j++) - frob += sqr (trans.Get(j)); - frob = sqrt (frob); + // Frobenius norm + double frob = 0; + for (int j = 1; j <= 4; j++) + frob += sqr (trans.Get(j)); + frob = sqrt (frob); - double dfrob = 0; - for (int j = 1; j <= 4; j++) - dfrob += trans.Get(j) * dtrans.Get(j); - dfrob = dfrob / frob; + double dfrob = 0; + for (int j = 1; j <= 4; j++) + dfrob += trans.Get(j) * dtrans.Get(j); + dfrob = dfrob / frob; - frob /= 2; - dfrob /= 2; + frob /= 2; + dfrob /= 2; - double det = trans(0,0)*trans(1,1)-trans(1,0)*trans(0,1); + double det = trans(0,0)*trans(1,1)-trans(1,0)*trans(0,1); - // ddet = \sum_j det (m_j) with m_j = trans, except col j = dtrans - double ddet - = dtrans(0,0) * trans(1,1) - trans(0,1) * dtrans(1,0) - + trans(0,0) * dtrans(1,1) - dtrans(0,1) * trans(1,0); + // ddet = \sum_j det (m_j) with m_j = trans, except col j = dtrans + double ddet + = dtrans(0,0) * trans(1,1) - trans(0,1) * dtrans(1,0) + + trans(0,0) * dtrans(1,1) - dtrans(0,1) * trans(1,0); - if (det <= 0) - err += 1e12; - else - { - err += frob * frob / det; - dd += (2 * frob * dfrob * det - frob * frob * ddet) / (det * det); - } - } + if (det <= 0) + err += 1e12; + else + { + err += frob * frob / det; + dd += (2 * frob * dfrob * det - frob * frob * ddet) / (det * det); + } + } - err /= nip; - dd /= nip; - return err; -} + err /= nip; + dd /= nip; + return err; + } -double Element2d :: -CalcJacobianBadness (const T_POINTS & points, const Vec<3> & n) const -{ - int i, j; - int nip = GetNIP(); - static DenseMatrix trans(2,2); - static DenseMatrix pmat; + double Element2d :: + CalcJacobianBadness (const T_POINTS & points, const Vec<3> & n) const + { + int i, j; + int nip = GetNIP(); + static DenseMatrix trans(2,2); + static DenseMatrix pmat; - pmat.SetSize (2, GetNP()); + pmat.SetSize (2, GetNP()); - Vec<3> t1, t2; - t1 = n.GetNormal(); - t2 = Cross (n, t1); + Vec<3> t1, t2; + t1 = n.GetNormal(); + t2 = Cross (n, t1); - for (i = 1; i <= GetNP(); i++) - { - Point3d p = points.Get(PNum(i)); - pmat.Elem(1, i) = p.X() * t1(0) + p.Y() * t1(1) + p.Z() * t1(2); - pmat.Elem(2, i) = p.X() * t2(0) + p.Y() * t2(1) + p.Z() * t2(2); - } + for (i = 1; i <= GetNP(); i++) + { + Point3d p = points.Get(PNum(i)); + pmat.Elem(1, i) = p.X() * t1(0) + p.Y() * t1(1) + p.Z() * t1(2); + pmat.Elem(2, i) = p.X() * t2(0) + p.Y() * t2(1) + p.Z() * t2(2); + } - double err = 0; - for (i = 1; i <= nip; i++) - { - GetTransformation (i, pmat, trans); + double err = 0; + for (i = 1; i <= nip; i++) + { + GetTransformation (i, pmat, trans); - // Frobenius norm - double frob = 0; - for (j = 1; j <= 4; j++) - frob += sqr (trans.Get(j)); - frob = sqrt (frob); - frob /= 2; + // Frobenius norm + double frob = 0; + for (j = 1; j <= 4; j++) + frob += sqr (trans.Get(j)); + frob = sqrt (frob); + frob /= 2; - double det = trans.Det(); - if (det <= 0) - err += 1e12; - else - err += frob * frob / det; - } + double det = trans.Det(); + if (det <= 0) + err += 1e12; + else + err += frob * frob / det; + } - err /= nip; - return err; -} + err /= nip; + return err; + } -void Element2d :: ComputeIntegrationPointData () const -{ - switch (np) - { - case 3: if (ipdtrig.Size()) return; break; - case 4: if (ipdquad.Size()) return; break; - } + void Element2d :: ComputeIntegrationPointData () const + { + switch (np) + { + case 3: if (ipdtrig.Size()) return; break; + case 4: if (ipdquad.Size()) return; break; + } - for (int i = 1; i <= GetNIP(); i++) - { - IntegrationPointData * ipd = new IntegrationPointData; - Point2d hp; - GetIntegrationPoint (i, hp, ipd->weight); - ipd->p(0) = hp.X(); - ipd->p(1) = hp.Y(); - ipd->p(2) = 0; + for (int i = 1; i <= GetNIP(); i++) + { + IntegrationPointData * ipd = new IntegrationPointData; + Point2d hp; + GetIntegrationPoint (i, hp, ipd->weight); + ipd->p(0) = hp.X(); + ipd->p(1) = hp.Y(); + ipd->p(2) = 0; - ipd->shape.SetSize(GetNP()); - ipd->dshape.SetSize(2, GetNP()); + ipd->shape.SetSize(GetNP()); + ipd->dshape.SetSize(2, GetNP()); - GetShape (hp, ipd->shape); - GetDShape (hp, ipd->dshape); + GetShape (hp, ipd->shape); + GetDShape (hp, ipd->dshape); - switch (np) - { - case 3: ipdtrig.Append (ipd); break; - case 4: ipdquad.Append (ipd); break; - } - } -} + switch (np) + { + case 3: ipdtrig.Append (ipd); break; + case 4: ipdquad.Append (ipd); break; + } + } + } @@ -915,86 +892,86 @@ void Element2d :: ComputeIntegrationPointData () const -ostream & operator<<(ostream & s, const Element2d & el) -{ - s << "np = " << el.GetNP(); - for (int j = 1; j <= el.GetNP(); j++) - s << " " << el.PNum(j); - return s; -} + ostream & operator<<(ostream & s, const Element2d & el) + { + s << "np = " << el.GetNP(); + for (int j = 1; j <= el.GetNP(); j++) + s << " " << el.PNum(j); + return s; + } -ostream & operator<<(ostream & s, const Element & el) -{ - s << "np = " << el.GetNP(); - for (int j = 0; j < el.GetNP(); j++) - s << " " << int(el[j]); - return s; -} + ostream & operator<<(ostream & s, const Element & el) + { + s << "np = " << el.GetNP(); + for (int j = 0; j < el.GetNP(); j++) + s << " " << int(el[j]); + return s; + } -Element :: Element () -{ - typ = TET; - np = 4; - for (int i = 0; i < ELEMENT_MAXPOINTS; i++) - pnum[i] = 0; - index = 0; - flags.marked = 1; - flags.badel = 0; - flags.reverse = 0; - flags.illegal = 0; - flags.illegal_valid = 0; - flags.badness_valid = 0; - flags.refflag = 1; - flags.strongrefflag = false; - flags.deleted = 0; - flags.fixed = 0; - orderx = ordery = orderz = 1; + Element :: Element () + { + typ = TET; + np = 4; + for (int i = 0; i < ELEMENT_MAXPOINTS; i++) + pnum[i] = 0; + index = 0; + flags.marked = 1; + flags.badel = 0; + flags.reverse = 0; + flags.illegal = 0; + flags.illegal_valid = 0; + flags.badness_valid = 0; + flags.refflag = 1; + flags.strongrefflag = false; + flags.deleted = 0; + flags.fixed = 0; + orderx = ordery = orderz = 1; #ifdef PARALLEL - partitionNumber = -1; - isghost = 0; + partitionNumber = -1; + isghost = 0; #endif -} + } -Element :: Element (int anp) -{ - np = anp; - int i; - for (i = 0; i < ELEMENT_MAXPOINTS; i++) - pnum[i] = 0; - index = 0; - flags.marked = 1; - flags.badel = 0; - flags.reverse = 0; - flags.illegal = 0; - flags.illegal_valid = 0; - flags.badness_valid = 0; - flags.refflag = 1; - flags.strongrefflag = false; - flags.deleted = 0; - flags.fixed = 0; + Element :: Element (int anp) + { + np = anp; + int i; + for (i = 0; i < ELEMENT_MAXPOINTS; i++) + pnum[i] = 0; + index = 0; + flags.marked = 1; + flags.badel = 0; + flags.reverse = 0; + flags.illegal = 0; + flags.illegal_valid = 0; + flags.badness_valid = 0; + flags.refflag = 1; + flags.strongrefflag = false; + flags.deleted = 0; + flags.fixed = 0; - switch (np) - { - case 4: typ = TET; break; - case 5: typ = PYRAMID; break; - case 6: typ = PRISM; break; - case 8: typ = HEX; break; - case 10: typ = TET10; break; - default: cerr << "Element::Element: unknown element with " << np << " points" << endl; - } - orderx = ordery = orderz = 1; + switch (np) + { + case 4: typ = TET; break; + case 5: typ = PYRAMID; break; + case 6: typ = PRISM; break; + case 8: typ = HEX; break; + case 10: typ = TET10; break; + default: cerr << "Element::Element: unknown element with " << np << " points" << endl; + } + orderx = ordery = orderz = 1; #ifdef PARALLEL - isghost = 0; + isghost = 0; #endif -} + } -void Element :: SetOrder (const int aorder) + void Element :: SetOrder (const int aorder) { orderx = aorder; ordery = aorder; @@ -1002,428 +979,428 @@ void Element :: SetOrder (const int aorder) } -void Element :: SetOrder (const int ox, const int oy, const int oz) -{ - orderx = ox; - ordery = oy; - orderz = oz; -} + void Element :: SetOrder (const int ox, const int oy, const int oz) + { + orderx = ox; + ordery = oy; + orderz = oz; + } -Element :: Element (ELEMENT_TYPE type) -{ - SetType (type); + Element :: Element (ELEMENT_TYPE type) + { + SetType (type); - int i; - for (i = 0; i < ELEMENT_MAXPOINTS; i++) - pnum[i] = 0; - index = 0; - flags.marked = 1; - flags.badel = 0; - flags.reverse = 0; - flags.illegal = 0; - flags.illegal_valid = 0; - flags.badness_valid = 0; - flags.refflag = 1; - flags.strongrefflag = false; - flags.deleted = 0; - flags.fixed = 0; - orderx = ordery = orderz = 1; + int i; + for (i = 0; i < ELEMENT_MAXPOINTS; i++) + pnum[i] = 0; + index = 0; + flags.marked = 1; + flags.badel = 0; + flags.reverse = 0; + flags.illegal = 0; + flags.illegal_valid = 0; + flags.badness_valid = 0; + flags.refflag = 1; + flags.strongrefflag = false; + flags.deleted = 0; + flags.fixed = 0; + orderx = ordery = orderz = 1; #ifdef PARALLEL - isghost = 0; + isghost = 0; #endif -} + } -Element & Element :: operator= (const Element & el2) -{ - typ = el2.typ; - np = el2.np; - for (int i = 0; i < ELEMENT_MAXPOINTS; i++) - pnum[i] = el2.pnum[i]; - index = el2.index; - flags = el2.flags; - orderx = el2.orderx; - ordery = el2.ordery; - orderz = el2.orderz; - hp_elnr = el2.hp_elnr; - flags = el2.flags; - return *this; -} - - - -void Element :: SetNP (int anp) -{ - np = anp; - switch (np) - { - case 4: typ = TET; break; - case 5: typ = PYRAMID; break; - case 6: typ = PRISM; break; - case 8: typ = HEX; break; - case 10: typ = TET10; break; - // - default: break; - cerr << "Element::SetNP unknown element with " << np << " points" << endl; - } -} - - - -void Element :: SetType (ELEMENT_TYPE atyp) -{ - typ = atyp; - switch (atyp) - { - case TET: np = 4; break; - case PYRAMID: np = 5; break; - case PRISM: np = 6; break; - case HEX: np = 8; break; - case TET10: np = 10; break; - case PRISM12: np = 12; break; - - default: break; - cerr << "Element::SetType unknown type " << int(typ) << endl; - } -} - - - -void Element :: Invert() -{ - switch (GetNP()) - { - case 4: - { - Swap (PNum(3), PNum(4)); - break; - } - case 5: - { - Swap (PNum(1), PNum(4)); - Swap (PNum(2), PNum(3)); - break; - } - case 6: - { - Swap (PNum(1), PNum(4)); - Swap (PNum(2), PNum(5)); - Swap (PNum(3), PNum(6)); - break; - } - } -} - - -void Element :: Print (ostream & ost) const -{ - ost << np << " Points: "; - for (int i = 1; i <= np; i++) - ost << pnum[i-1] << " " << endl; -} - -void Element :: GetBox (const T_POINTS & points, Box3d & box) const -{ - box.SetPoint (points.Get(PNum(1))); - box.AddPoint (points.Get(PNum(2))); - box.AddPoint (points.Get(PNum(3))); - box.AddPoint (points.Get(PNum(4))); -} - -double Element :: Volume (const T_POINTS & points) const -{ - Vec<3> v1 = points.Get(PNum(2)) - points.Get(PNum(1)); - Vec<3> v2 = points.Get(PNum(3)) - points.Get(PNum(1)); - Vec<3> v3 = points.Get(PNum(4)) - points.Get(PNum(1)); - - return -(Cross (v1, v2) * v3) / 6; -} - - -void Element :: GetFace2 (int i, Element2d & face) const -{ - static const int tetfaces[][5] = - { { 3, 2, 3, 4, 0 }, - { 3, 3, 1, 4, 0 }, - { 3, 1, 2, 4, 0 }, - { 3, 2, 1, 3, 0 } }; - - static const int tet10faces[][7] = - { { 3, 2, 3, 4, 10, 9, 8 }, - { 3, 3, 1, 4, 7, 10, 6 }, - { 3, 1, 2, 4, 9, 7, 5 }, - { 3, 2, 1, 3, 6, 8, 5 } }; - - static const int pyramidfaces[][5] = - { { 4, 1, 4, 3, 2 }, - { 3, 1, 2, 5, 0 }, - { 3, 2, 3, 5, 0 }, - { 3, 3, 4, 5, 0 }, - { 3, 4, 1, 5, 0 } }; - - static const int prismfaces[][5] = + Element & Element :: operator= (const Element & el2) { - { 3, 1, 3, 2, 0 }, - { 3, 4, 5, 6, 0 }, - { 4, 1, 2, 5, 4 }, - { 4, 2, 3, 6, 5 }, - { 4, 3, 1, 4, 6 } - }; - - switch (np) - { - case 4: // tet - { - face.SetType(TRIG); - for (int j = 1; j <= 3; j++) - face.PNum(j) = PNum(tetfaces[i-1][j]); - break; - } - - case 10: // tet10 - { - face.SetType(TRIG6); - for (int j = 1; j <= 6; j++) - face.PNum(j) = PNum(tet10faces[i-1][j]); - break; - } - - case 5: // pyramid - { - // face.SetNP(pyramidfaces[i-1][0]); - face.SetType ( (i == 1) ? QUAD : TRIG); - for (int j = 1; j <= face.GetNP(); j++) - face.PNum(j) = PNum(pyramidfaces[i-1][j]); - break; - } - case 6: // prism - { - // face.SetNP(prismfaces[i-1][0]); - face.SetType ( (i >= 3) ? QUAD : TRIG); - for (int j = 1; j <= face.GetNP(); j++) - face.PNum(j) = PNum(prismfaces[i-1][j]); - break; - } - } -} + typ = el2.typ; + np = el2.np; + for (int i = 0; i < ELEMENT_MAXPOINTS; i++) + pnum[i] = el2.pnum[i]; + index = el2.index; + flags = el2.flags; + orderx = el2.orderx; + ordery = el2.ordery; + orderz = el2.orderz; + hp_elnr = el2.hp_elnr; + flags = el2.flags; + return *this; + } -void Element :: GetTets (Array & locels) const -{ - GetTetsLocal (locels); - int i, j; - for (i = 1; i <= locels.Size(); i++) - for (j = 1; j <= 4; j++) - locels.Elem(i).PNum(j) = PNum ( locels.Elem(i).PNum(j) ); -} + void Element :: SetNP (int anp) + { + np = anp; + switch (np) + { + case 4: typ = TET; break; + case 5: typ = PYRAMID; break; + case 6: typ = PRISM; break; + case 8: typ = HEX; break; + case 10: typ = TET10; break; + // + default: break; + cerr << "Element::SetNP unknown element with " << np << " points" << endl; + } + } -void Element :: GetTetsLocal (Array & locels) const -{ - int i, j; - locels.SetSize(0); - switch (GetType()) - { - case TET: - { - int linels[1][4] = - { { 1, 2, 3, 4 }, - }; - for (i = 0; i < 1; i++) - { - Element tet(4); - for (j = 1; j <= 4; j++) - tet.PNum(j) = linels[i][j-1]; - locels.Append (tet); - } - break; - } - case TET10: - { - int linels[8][4] = - { { 1, 5, 6, 7 }, - { 5, 2, 8, 9 }, - { 6, 8, 3, 10 }, - { 7, 9, 10, 4 }, - { 5, 6, 7, 9 }, - { 5, 6, 9, 8 }, - { 6, 7, 9, 10 }, - { 6, 8, 10, 9 } }; - for (i = 0; i < 8; i++) - { - Element tet(4); - for (j = 1; j <= 4; j++) - tet.PNum(j) = linels[i][j-1]; - locels.Append (tet); - } - break; - } - case PYRAMID: - { - int linels[2][4] = - { { 1, 2, 3, 5 }, - { 1, 3, 4, 5 } }; - for (i = 0; i < 2; i++) - { - Element tet(4); - for (j = 1; j <= 4; j++) - tet.PNum(j) = linels[i][j-1]; - locels.Append (tet); - } - break; - } - case PRISM: - case PRISM12: - { - int linels[3][4] = - { { 1, 2, 3, 4 }, - { 4, 2, 3, 5 }, - { 6, 5, 4, 3 } - }; - for (i = 0; i < 3; i++) - { - Element tet(4); - for (j = 0; j < 4; j++) - tet[j] = linels[i][j]; - locels.Append (tet); - } - break; - } - case HEX: - { - int linels[6][4] = - { { 1, 7, 2, 3 }, - { 1, 7, 3, 4 }, - { 1, 7, 4, 8 }, - { 1, 7, 8, 5 }, - { 1, 7, 5, 6 }, - { 1, 7, 6, 2 } - }; - for (i = 0; i < 6; i++) - { - Element tet(4); - for (j = 0; j < 4; j++) - tet[j] = linels[i][j]; - locels.Append (tet); - } - break; - } - default: - { - cerr << "GetTetsLocal not implemented for el with " << GetNP() << " nodes" << endl; - } - } -} -bool Element :: operator==(const Element & el2) const -{ - bool retval = (el2.GetNP() == np); - for(int i= 0; retval && i v1 = points.Get(PNum(2)) - points.Get(PNum(1)); + Vec<3> v2 = points.Get(PNum(3)) - points.Get(PNum(1)); + Vec<3> v3 = points.Get(PNum(4)) - points.Get(PNum(1)); + + return -(Cross (v1, v2) * v3) / 6; + } + + + void Element :: GetFace2 (int i, Element2d & face) const + { + static const int tetfaces[][5] = + { { 3, 2, 3, 4, 0 }, + { 3, 3, 1, 4, 0 }, + { 3, 1, 2, 4, 0 }, + { 3, 2, 1, 3, 0 } }; + + static const int tet10faces[][7] = + { { 3, 2, 3, 4, 10, 9, 8 }, + { 3, 3, 1, 4, 7, 10, 6 }, + { 3, 1, 2, 4, 9, 7, 5 }, + { 3, 2, 1, 3, 6, 8, 5 } }; + + static const int pyramidfaces[][5] = + { { 4, 1, 4, 3, 2 }, + { 3, 1, 2, 5, 0 }, + { 3, 2, 3, 5, 0 }, + { 3, 3, 4, 5, 0 }, + { 3, 4, 1, 5, 0 } }; + + static const int prismfaces[][5] = + { + { 3, 1, 3, 2, 0 }, + { 3, 4, 5, 6, 0 }, + { 4, 1, 2, 5, 4 }, + { 4, 2, 3, 6, 5 }, + { 4, 3, 1, 4, 6 } + }; + + switch (np) + { + case 4: // tet + { + face.SetType(TRIG); + for (int j = 1; j <= 3; j++) + face.PNum(j) = PNum(tetfaces[i-1][j]); + break; + } + + case 10: // tet10 + { + face.SetType(TRIG6); + for (int j = 1; j <= 6; j++) + face.PNum(j) = PNum(tet10faces[i-1][j]); + break; + } + + case 5: // pyramid + { + // face.SetNP(pyramidfaces[i-1][0]); + face.SetType ( (i == 1) ? QUAD : TRIG); + for (int j = 1; j <= face.GetNP(); j++) + face.PNum(j) = PNum(pyramidfaces[i-1][j]); + break; + } + case 6: // prism + { + // face.SetNP(prismfaces[i-1][0]); + face.SetType ( (i >= 3) ? QUAD : TRIG); + for (int j = 1; j <= face.GetNP(); j++) + face.PNum(j) = PNum(prismfaces[i-1][j]); + break; + } + } + } + + + + void Element :: GetTets (Array & locels) const + { + GetTetsLocal (locels); + int i, j; + for (i = 1; i <= locels.Size(); i++) + for (j = 1; j <= 4; j++) + locels.Elem(i).PNum(j) = PNum ( locels.Elem(i).PNum(j) ); + } + + void Element :: GetTetsLocal (Array & locels) const + { + int i, j; + locels.SetSize(0); + switch (GetType()) + { + case TET: + { + int linels[1][4] = + { { 1, 2, 3, 4 }, + }; + for (i = 0; i < 1; i++) + { + Element tet(4); + for (j = 1; j <= 4; j++) + tet.PNum(j) = linels[i][j-1]; + locels.Append (tet); + } + break; + } + case TET10: + { + int linels[8][4] = + { { 1, 5, 6, 7 }, + { 5, 2, 8, 9 }, + { 6, 8, 3, 10 }, + { 7, 9, 10, 4 }, + { 5, 6, 7, 9 }, + { 5, 6, 9, 8 }, + { 6, 7, 9, 10 }, + { 6, 8, 10, 9 } }; + for (i = 0; i < 8; i++) + { + Element tet(4); + for (j = 1; j <= 4; j++) + tet.PNum(j) = linels[i][j-1]; + locels.Append (tet); + } + break; + } + case PYRAMID: + { + int linels[2][4] = + { { 1, 2, 3, 5 }, + { 1, 3, 4, 5 } }; + for (i = 0; i < 2; i++) + { + Element tet(4); + for (j = 1; j <= 4; j++) + tet.PNum(j) = linels[i][j-1]; + locels.Append (tet); + } + break; + } + case PRISM: + case PRISM12: + { + int linels[3][4] = + { { 1, 2, 3, 4 }, + { 4, 2, 3, 5 }, + { 6, 5, 4, 3 } + }; + for (i = 0; i < 3; i++) + { + Element tet(4); + for (j = 0; j < 4; j++) + tet[j] = linels[i][j]; + locels.Append (tet); + } + break; + } + case HEX: + { + int linels[6][4] = + { { 1, 7, 2, 3 }, + { 1, 7, 3, 4 }, + { 1, 7, 4, 8 }, + { 1, 7, 8, 5 }, + { 1, 7, 5, 6 }, + { 1, 7, 6, 2 } + }; + for (i = 0; i < 6; i++) + { + Element tet(4); + for (j = 0; j < 4; j++) + tet[j] = linels[i][j]; + locels.Append (tet); + } + break; + } + default: + { + cerr << "GetTetsLocal not implemented for el with " << GetNP() << " nodes" << endl; + } + } + } + + bool Element :: operator==(const Element & el2) const + { + bool retval = (el2.GetNP() == np); + for(int i= 0; retval && i & points) const -{ - const static double tetpoints[4][3] = - { { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 }}; + void Element :: GetNodesLocal (Array & points) const + { + const static double tetpoints[4][3] = + { { 0, 0, 0 }, + { 1, 0, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 }}; - const static double prismpoints[6][3] = - { { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 }, - { 1, 0, 1 }, - { 0, 1, 1 } }; + const static double prismpoints[6][3] = + { { 0, 0, 0 }, + { 1, 0, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 }, + { 1, 0, 1 }, + { 0, 1, 1 } }; - const static double pyramidpoints[6][3] = - { { 0, 0, 0 }, - { 1, 0, 0 }, - { 1, 1, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 } }; + const static double pyramidpoints[6][3] = + { { 0, 0, 0 }, + { 1, 0, 0 }, + { 1, 1, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 } }; - const static double tet10points[10][3] = - { { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 }, - { 0.5, 0, 0 }, - { 0, 0.5, 0 }, - { 0, 0, 0.5 }, - { 0.5, 0.5, 0 }, - { 0.5, 0, 0.5 }, - { 0, 0.5, 0.5 } }; + const static double tet10points[10][3] = + { { 0, 0, 0 }, + { 1, 0, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 }, + { 0.5, 0, 0 }, + { 0, 0.5, 0 }, + { 0, 0, 0.5 }, + { 0.5, 0.5, 0 }, + { 0.5, 0, 0.5 }, + { 0, 0.5, 0.5 } }; - const static double hexpoints[8][3] = - { - { 0, 0, 0 }, - { 1, 0, 0 }, - { 1, 1, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 }, - { 1, 0, 1 }, - { 1, 1, 1 }, - { 0, 1, 1 } - }; + const static double hexpoints[8][3] = + { + { 0, 0, 0 }, + { 1, 0, 0 }, + { 1, 1, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 }, + { 1, 0, 1 }, + { 1, 1, 1 }, + { 0, 1, 1 } + }; - int np, i; - const double (*pp)[3]; - switch (GetType()) - { - case TET: + int np, i; + const double (*pp)[3]; + switch (GetType()) { - np = 4; - pp = tetpoints; - break; + case TET: + { + np = 4; + pp = tetpoints; + break; + } + case PRISM: + case PRISM12: + { + np = 6; + pp = prismpoints; + break; + } + case TET10: + { + np = 10; + pp = tet10points; + break; + } + case PYRAMID: + { + np = 5; + pp = pyramidpoints; + break; + } + case HEX: + { + np = 8; + pp = hexpoints; + break; + } + default: + { + cout << "GetNodesLocal not impelemented for element " << GetType() << endl; + np = 0; + } } - case PRISM: - case PRISM12: - { - np = 6; - pp = prismpoints; - break; - } - case TET10: - { - np = 10; - pp = tet10points; - break; - } - case PYRAMID: - { - np = 5; - pp = pyramidpoints; - break; - } - case HEX: - { - np = 8; - pp = hexpoints; - break; - } - default: - { - cout << "GetNodesLocal not impelemented for element " << GetType() << endl; - np = 0; - } - } - points.SetSize(0); - for (i = 0; i < np; i++) - points.Append (Point3d (pp[i][0], pp[i][1], pp[i][2])); -} + points.SetSize(0); + for (i = 0; i < np; i++) + points.Append (Point3d (pp[i][0], pp[i][1], pp[i][2])); + } #endif @@ -1431,1245 +1408,1239 @@ void Element :: GetNodesLocal (Array & points) const -void Element :: GetNodesLocalNew (Array > & points) const -{ - const static double tetpoints[4][3] = - { - { 1, 0, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 }, - { 0, 0, 0 } - }; - - const static double prismpoints[6][3] = - { - { 1, 0, 0 }, - { 0, 1, 0 }, - { 0, 0, 0 }, - { 1, 0, 1 }, - { 0, 1, 1 }, - { 0, 0, 1 } - }; - - const static double pyramidpoints[6][3] = - { { 0, 0, 0 }, - { 1, 0, 0 }, - { 1, 1, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 } }; - - const static double tet10points[10][3] = - { { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 }, - { 0.5, 0, 0 }, - { 0, 0.5, 0 }, - { 0, 0, 0.5 }, - { 0.5, 0.5, 0 }, - { 0.5, 0, 0.5 }, - { 0, 0.5, 0.5 } }; - - const static double hexpoints[8][3] = - { - { 0, 0, 0 }, - { 1, 0, 0 }, - { 1, 1, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 }, - { 1, 0, 1 }, - { 1, 1, 1 }, - { 0, 1, 1 } - }; - - - - int np, i; - const double (*pp)[3]; - switch (GetType()) - { - case TET: - { - np = 4; - pp = tetpoints; - break; - } - case PRISM: - case PRISM12: - { - np = 6; - pp = prismpoints; - break; - } - case TET10: - { - np = 10; - pp = tet10points; - break; - } - case PYRAMID: - { - np = 5; - pp = pyramidpoints; - break; - } - case HEX: - { - np = 8; - pp = hexpoints; - break; - } - default: - { - cout << "GetNodesLocal not impelemented for element " << GetType() << endl; - np = 0; - } - } - - points.SetSize(0); - for (i = 0; i < np; i++) - points.Append (Point<3> (pp[i][0], pp[i][1], pp[i][2])); -} - - - - - - - - - - - - - - - - - -void Element :: GetSurfaceTriangles (Array & surftrigs) const -{ - static int tet4trigs[][3] = - { { 2, 3, 4 }, - { 3, 1, 4 }, - { 1, 2, 4 }, - { 2, 1, 3 } }; - - static int tet10trigs[][3] = - { { 2, 8, 9 }, { 3, 10, 8}, { 4, 9, 10 }, { 9, 8, 10 }, - { 3, 6, 10 }, { 1, 7, 6 }, { 4, 10, 7 }, { 6, 7, 10 }, - { 1, 5, 7 }, { 2, 9, 5 }, { 4, 7, 9 }, { 5, 9, 7 }, - { 1, 6, 5 }, { 2, 5, 8 }, { 3, 8, 6 }, { 5, 6, 8 } - }; - - static int pyramidtrigs[][3] = + void Element :: GetNodesLocalNew (Array > & points) const { - { 1, 3, 2 }, - { 1, 4, 3 }, - { 1, 2, 5 }, - { 2, 3, 5 }, - { 3, 4, 5 }, - { 4, 1, 5 } - }; - - static int prismtrigs[][3] = - { - { 1, 3, 2 }, - { 4, 5, 6 }, - { 1, 2, 4 }, - { 4, 2, 5 }, - { 2, 3, 5 }, - { 5, 3, 6 }, - { 3, 1, 6 }, - { 6, 1, 4 } - }; + const static double tetpoints[4][3] = + { + { 1, 0, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 }, + { 0, 0, 0 } + }; - static int hextrigs[][3] = - { - { 1, 3, 2 }, - { 1, 4, 3 }, - { 5, 6, 7 }, - { 5, 7, 8 }, - { 1, 2, 6 }, - { 1, 6, 5 }, - { 2, 3, 7 }, - { 2, 7, 6 }, - { 3, 4, 8 }, - { 3, 8, 7 }, - { 4, 1, 8 }, - { 1, 5, 8 } - }; + const static double prismpoints[6][3] = + { + { 1, 0, 0 }, + { 0, 1, 0 }, + { 0, 0, 0 }, + { 1, 0, 1 }, + { 0, 1, 1 }, + { 0, 0, 1 } + }; + + const static double pyramidpoints[6][3] = + { { 0, 0, 0 }, + { 1, 0, 0 }, + { 1, 1, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 } }; + + const static double tet10points[10][3] = + { { 0, 0, 0 }, + { 1, 0, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 }, + { 0.5, 0, 0 }, + { 0, 0.5, 0 }, + { 0, 0, 0.5 }, + { 0.5, 0.5, 0 }, + { 0.5, 0, 0.5 }, + { 0, 0.5, 0.5 } }; - int j; - - int nf; - int (*fp)[3]; - - switch (GetType()) - { - case TET: - { - nf = 4; - fp = tet4trigs; - break; - } - case PYRAMID: - { - nf = 6; - fp = pyramidtrigs; - break; - } - case PRISM: - case PRISM12: - { - nf = 8; - fp = prismtrigs; - break; - } - case TET10: - { - nf = 16; - fp = tet10trigs; - break; - } - case HEX: - { - nf = 12; - fp = hextrigs; - break; - } - default: - { - nf = 0; - fp = NULL; - } - } + const static double hexpoints[8][3] = + { + { 0, 0, 0 }, + { 1, 0, 0 }, + { 1, 1, 0 }, + { 0, 1, 0 }, + { 0, 0, 1 }, + { 1, 0, 1 }, + { 1, 1, 1 }, + { 0, 1, 1 } + }; + - surftrigs.SetSize (nf); - for (j = 0; j < nf; j++) - { - surftrigs.Elem(j+1) = Element2d(TRIG); - surftrigs.Elem(j+1).PNum(1) = fp[j][0]; - surftrigs.Elem(j+1).PNum(2) = fp[j][1]; - surftrigs.Elem(j+1).PNum(3) = fp[j][2]; - } -} + int np, i; + const double (*pp)[3]; + switch (GetType()) + { + case TET: + { + np = 4; + pp = tetpoints; + break; + } + case PRISM: + case PRISM12: + { + np = 6; + pp = prismpoints; + break; + } + case TET10: + { + np = 10; + pp = tet10points; + break; + } + case PYRAMID: + { + np = 5; + pp = pyramidpoints; + break; + } + case HEX: + { + np = 8; + pp = hexpoints; + break; + } + default: + { + cout << "GetNodesLocal not impelemented for element " << GetType() << endl; + np = 0; + } + } + + points.SetSize(0); + for (i = 0; i < np; i++) + points.Append (Point<3> (pp[i][0], pp[i][1], pp[i][2])); + } -Array< AutoPtr < IntegrationPointData > > ipdtet; -Array< AutoPtr < IntegrationPointData > > ipdtet10; -int Element :: GetNIP () const -{ - int nip; - switch (typ) - { - case TET: nip = 1; break; - case TET10: nip = 8; break; - default: nip = 0; break; - } - return nip; -} -void Element :: -GetIntegrationPoint (int ip, Point<3> & p, double & weight) const -{ - static double eltetqp[1][4] = + + + + + + + + + void Element :: GetSurfaceTriangles (Array & surftrigs) const { - { 0.25, 0.25, 0.25, 1.0/6.0 } - }; + static int tet4trigs[][3] = + { { 2, 3, 4 }, + { 3, 1, 4 }, + { 1, 2, 4 }, + { 2, 1, 3 } }; - static double eltet10qp[8][4] = + static int tet10trigs[][3] = + { { 2, 8, 9 }, { 3, 10, 8}, { 4, 9, 10 }, { 9, 8, 10 }, + { 3, 6, 10 }, { 1, 7, 6 }, { 4, 10, 7 }, { 6, 7, 10 }, + { 1, 5, 7 }, { 2, 9, 5 }, { 4, 7, 9 }, { 5, 9, 7 }, + { 1, 6, 5 }, { 2, 5, 8 }, { 3, 8, 6 }, { 5, 6, 8 } + }; + + static int pyramidtrigs[][3] = + { + { 1, 3, 2 }, + { 1, 4, 3 }, + { 1, 2, 5 }, + { 2, 3, 5 }, + { 3, 4, 5 }, + { 4, 1, 5 } + }; + + static int prismtrigs[][3] = + { + { 1, 3, 2 }, + { 4, 5, 6 }, + { 1, 2, 4 }, + { 4, 2, 5 }, + { 2, 3, 5 }, + { 5, 3, 6 }, + { 3, 1, 6 }, + { 6, 1, 4 } + }; + + static int hextrigs[][3] = + { + { 1, 3, 2 }, + { 1, 4, 3 }, + { 5, 6, 7 }, + { 5, 7, 8 }, + { 1, 2, 6 }, + { 1, 6, 5 }, + { 2, 3, 7 }, + { 2, 7, 6 }, + { 3, 4, 8 }, + { 3, 8, 7 }, + { 4, 1, 8 }, + { 1, 5, 8 } + }; + + int j; + + int nf; + int (*fp)[3]; + + switch (GetType()) + { + case TET: + { + nf = 4; + fp = tet4trigs; + break; + } + case PYRAMID: + { + nf = 6; + fp = pyramidtrigs; + break; + } + case PRISM: + case PRISM12: + { + nf = 8; + fp = prismtrigs; + break; + } + case TET10: + { + nf = 16; + fp = tet10trigs; + break; + } + case HEX: + { + nf = 12; + fp = hextrigs; + break; + } + default: + { + nf = 0; + fp = NULL; + } + } + + + surftrigs.SetSize (nf); + for (j = 0; j < nf; j++) + { + surftrigs.Elem(j+1) = Element2d(TRIG); + surftrigs.Elem(j+1).PNum(1) = fp[j][0]; + surftrigs.Elem(j+1).PNum(2) = fp[j][1]; + surftrigs.Elem(j+1).PNum(3) = fp[j][2]; + } + } + + + + + + Array< AutoPtr < IntegrationPointData > > ipdtet; + Array< AutoPtr < IntegrationPointData > > ipdtet10; + + + + int Element :: GetNIP () const { - { 0.585410196624969, 0.138196601125011, 0.138196601125011, 1.0/24.0 }, - { 0.138196601125011, 0.585410196624969, 0.138196601125011, 1.0/24.0 }, - { 0.138196601125011, 0.138196601125011, 0.585410196624969, 1.0/24.0 }, - { 0.138196601125011, 0.138196601125011, 0.138196601125011, 1.0/24.0 }, - { 1, 0, 0, 1 }, - { 0, 1, 0, 1 }, - { 0, 0, 1, 1 }, - { 0, 0, 0, 1 }, - }; + int nip; + switch (typ) + { + case TET: nip = 1; break; + case TET10: nip = 8; break; + default: nip = 0; break; + } + return nip; + } + + void Element :: + GetIntegrationPoint (int ip, Point<3> & p, double & weight) const + { + static double eltetqp[1][4] = + { + { 0.25, 0.25, 0.25, 1.0/6.0 } + }; + + static double eltet10qp[8][4] = + { + { 0.585410196624969, 0.138196601125011, 0.138196601125011, 1.0/24.0 }, + { 0.138196601125011, 0.585410196624969, 0.138196601125011, 1.0/24.0 }, + { 0.138196601125011, 0.138196601125011, 0.585410196624969, 1.0/24.0 }, + { 0.138196601125011, 0.138196601125011, 0.138196601125011, 1.0/24.0 }, + { 1, 0, 0, 1 }, + { 0, 1, 0, 1 }, + { 0, 0, 1, 1 }, + { 0, 0, 0, 1 }, + }; - double * pp; - switch (typ) - { - case TET: pp = &eltetqp[0][0]; break; - case TET10: pp = &eltet10qp[ip-1][0]; break; - } + double * pp = NULL; + switch (typ) + { + case TET: pp = &eltetqp[0][0]; break; + case TET10: pp = &eltet10qp[ip-1][0]; break; + } - p(0) = pp[0]; - p(1) = pp[1]; - p(2) = pp[2]; - weight = pp[3]; -} + p(0) = pp[0]; + p(1) = pp[1]; + p(2) = pp[2]; + weight = pp[3]; + } -void Element :: -GetTransformation (int ip, const T_POINTS & points, - DenseMatrix & trans) const -{ - int np = GetNP(); - static DenseMatrix pmat(3, np), dshape(3, np); - pmat.SetSize (3, np); - dshape.SetSize (3, np); + void Element :: + GetTransformation (int ip, const T_POINTS & points, + DenseMatrix & trans) const + { + int np = GetNP(); + static DenseMatrix pmat(3, np), dshape(3, np); + pmat.SetSize (3, np); + dshape.SetSize (3, np); - Point<3> p; - double w; + Point<3> p; + double w; - GetPointMatrix (points, pmat); - GetIntegrationPoint (ip, p, w); - GetDShape (p, dshape); + GetPointMatrix (points, pmat); + GetIntegrationPoint (ip, p, w); + GetDShape (p, dshape); - CalcABt (pmat, dshape, trans); + CalcABt (pmat, dshape, trans); - /* - (*testout) << "p = " << p << endl - << "pmat = " << pmat << endl - << "dshape = " << dshape << endl - << "tans = " << trans << endl; - */ -} + /* + (*testout) << "p = " << p << endl + << "pmat = " << pmat << endl + << "dshape = " << dshape << endl + << "tans = " << trans << endl; + */ + } -void Element :: -GetTransformation (int ip, class DenseMatrix & pmat, - class DenseMatrix & trans) const -{ - int np = GetNP(); + void Element :: + GetTransformation (int ip, class DenseMatrix & pmat, + class DenseMatrix & trans) const + { + int np = GetNP(); - if (pmat.Width() != np || pmat.Height() != 3) - { - (*testout) << "GetTransofrmation: pmat doesn't fit" << endl; - return; - } + if (pmat.Width() != np || pmat.Height() != 3) + { + (*testout) << "GetTransofrmation: pmat doesn't fit" << endl; + return; + } - ComputeIntegrationPointData (); - DenseMatrix * dshapep; - switch (GetType()) - { - case TET: dshapep = &ipdtet.Get(ip)->dshape; break; - case TET10: dshapep = &ipdtet10.Get(ip)->dshape; break; - default: - PrintSysError ("Element::GetTransformation, illegal type ", int(typ)); - } + ComputeIntegrationPointData (); + DenseMatrix * dshapep = 0; + switch (GetType()) + { + case TET: dshapep = &ipdtet.Get(ip)->dshape; break; + case TET10: dshapep = &ipdtet10.Get(ip)->dshape; break; + default: + PrintSysError ("Element::GetTransformation, illegal type ", int(typ)); + } - CalcABt (pmat, *dshapep, trans); -} + CalcABt (pmat, *dshapep, trans); + } -void Element :: GetShape (const Point<3> & hp, Vector & shape) const -{ - Point3d p = hp; + void Element :: GetShape (const Point<3> & hp, Vector & shape) const + { + Point3d p = hp; - if (shape.Size() != GetNP()) - { + if (shape.Size() != GetNP()) + { + cerr << "Element::GetShape: Length not fitting" << endl; + return; + } + + switch (typ) + { + case TET: + { + shape(0) = 1 - p.X() - p.Y() - p.Z(); + shape(1) = p.X(); + shape(2) = p.Y(); + shape(3) = p.Z(); + break; + } + case TET10: + { + double lam1 = 1 - p.X() - p.Y() - p.Z(); + double lam2 = p.X(); + double lam3 = p.Y(); + double lam4 = p.Z(); + + shape(4) = 4 * lam1 * lam2; + shape(5) = 4 * lam1 * lam3; + shape(6) = 4 * lam1 * lam4; + shape(7) = 4 * lam2 * lam3; + shape(8) = 4 * lam2 * lam4; + shape(9) = 4 * lam3 * lam4; + + shape(0) = lam1 - 0.5 * (shape(4) + shape(5) + shape(6)); + shape(1) = lam2 - 0.5 * (shape(4) + shape(7) + shape(8)); + shape(2) = lam3 - 0.5 * (shape(5) + shape(7) + shape(9)); + shape(3) = lam4 - 0.5 * (shape(6) + shape(8) + shape(9)); + break; + } + + case PRISM: + { + Point<3> hp = p; + shape(0) = hp(0) * (1-hp(2)); + shape(1) = hp(1) * (1-hp(2)); + shape(2) = (1-hp(0)-hp(1)) * (1-hp(2)); + shape(3) = hp(0) * hp(2); + shape(4) = hp(1) * hp(2); + shape(5) = (1-hp(0)-hp(1)) * hp(2); + break; + } + case HEX: + { + Point<3> hp = p; + shape(0) = (1-hp(0))*(1-hp(1))*(1-hp(2)); + shape(1) = ( hp(0))*(1-hp(1))*(1-hp(2)); + shape(2) = ( hp(0))*( hp(1))*(1-hp(2)); + shape(3) = (1-hp(0))*( hp(1))*(1-hp(2)); + shape(4) = (1-hp(0))*(1-hp(1))*( hp(2)); + shape(5) = ( hp(0))*(1-hp(1))*( hp(2)); + shape(6) = ( hp(0))*( hp(1))*( hp(2)); + shape(7) = (1-hp(0))*( hp(1))*( hp(2)); + break; + } + } + } + + + + void Element :: GetShapeNew (const Point<3> & p, FlatVector & shape) const + { + /* + if (shape.Size() < GetNP()) + { cerr << "Element::GetShape: Length not fitting" << endl; return; - } - - switch (typ) - { - case TET: - { - shape.Elem(1) = 1 - p.X() - p.Y() - p.Z(); - shape.Elem(2) = p.X(); - shape.Elem(3) = p.Y(); - shape.Elem(4) = p.Z(); - break; } - case TET10: + */ + + switch (typ) { - double lam1 = 1 - p.X() - p.Y() - p.Z(); - double lam2 = p.X(); - double lam3 = p.Y(); - double lam4 = p.Z(); + case TET: + { + shape(0) = p(0); + shape(1) = p(1); + shape(2) = p(2); + shape(3) = 1-p(0)-p(1)-p(2); + break; + } + + case TET10: + { + double lam1 = p(0); + double lam2 = p(1); + double lam3 = p(2); + double lam4 = 1-p(0)-p(1)-p(2); - shape.Elem(5) = 4 * lam1 * lam2; - shape.Elem(6) = 4 * lam1 * lam3; - shape.Elem(7) = 4 * lam1 * lam4; - shape.Elem(8) = 4 * lam2 * lam3; - shape.Elem(9) = 4 * lam2 * lam4; - shape.Elem(10) = 4 * lam3 * lam4; + shape(0) = 2 * lam1 * (lam1-0.5); + shape(1) = 2 * lam2 * (lam2-0.5); + shape(2) = 2 * lam3 * (lam3-0.5); + shape(3) = 2 * lam4 * (lam4-0.5); + + shape(4) = 4 * lam1 * lam2; + shape(5) = 4 * lam1 * lam3; + shape(6) = 4 * lam1 * lam4; + shape(7) = 4 * lam2 * lam3; + shape(8) = 4 * lam2 * lam4; + shape(9) = 4 * lam3 * lam4; - shape.Elem(1) = lam1 - - 0.5 * (shape.Elem(5) + shape.Elem(6) + shape.Elem(7)); - shape.Elem(2) = lam2 - - 0.5 * (shape.Elem(5) + shape.Elem(8) + shape.Elem(9)); - shape.Elem(3) = lam3 - - 0.5 * (shape.Elem(6) + shape.Elem(8) + shape.Elem(10)); - shape.Elem(4) = lam4 - - 0.5 * (shape.Elem(7) + shape.Elem(9) + shape.Elem(10)); - break; - } + break; + } - case PRISM: + + case PYRAMID: + { + double noz = 1-p(2); + if (noz == 0.0) noz = 1e-10; + + double xi = p(0) / noz; + double eta = p(1) / noz; + shape(0) = (1-xi)*(1-eta) * (noz); + shape(1) = ( xi)*(1-eta) * (noz); + shape(2) = ( xi)*( eta) * (noz); + shape(3) = (1-xi)*( eta) * (noz); + shape(4) = p(2); + break; + } + + case PRISM: + { + shape(0) = p(0) * (1-p(2)); + shape(1) = p(1) * (1-p(2)); + shape(2) = (1-p(0)-p(1)) * (1-p(2)); + shape(3) = p(0) * p(2); + shape(4) = p(1) * p(2); + shape(5) = (1-p(0)-p(1)) * p(2); + break; + } + case HEX: + { + shape(0) = (1-p(0))*(1-p(1))*(1-p(2)); + shape(1) = ( p(0))*(1-p(1))*(1-p(2)); + shape(2) = ( p(0))*( p(1))*(1-p(2)); + shape(3) = (1-p(0))*( p(1))*(1-p(2)); + shape(4) = (1-p(0))*(1-p(1))*( p(2)); + shape(5) = ( p(0))*(1-p(1))*( p(2)); + shape(6) = ( p(0))*( p(1))*( p(2)); + shape(7) = (1-p(0))*( p(1))*( p(2)); + break; + } + } + } + + + + + void Element :: + GetDShape (const Point<3> & hp, DenseMatrix & dshape) const + { + Point3d p = hp; + + int np = GetNP(); + if (dshape.Height() != 3 || dshape.Width() != np) { - Point<3> hp = p; - shape(0) = hp(0) * (1-hp(2)); - shape(1) = hp(1) * (1-hp(2)); - shape(2) = (1-hp(0)-hp(1)) * (1-hp(2)); - shape(3) = hp(0) * hp(2); - shape(4) = hp(1) * hp(2); - shape(5) = (1-hp(0)-hp(1)) * hp(2); - break; + cerr << "Element::DShape: Sizes don't fit" << endl; + return; } - case HEX: + + double eps = 1e-6; + Vector shaper(np), shapel(np); + + for (int i = 1; i <= 3; i++) { - Point<3> hp = p; - shape(0) = (1-hp(0))*(1-hp(1))*(1-hp(2)); - shape(1) = ( hp(0))*(1-hp(1))*(1-hp(2)); - shape(2) = ( hp(0))*( hp(1))*(1-hp(2)); - shape(3) = (1-hp(0))*( hp(1))*(1-hp(2)); - shape(4) = (1-hp(0))*(1-hp(1))*( hp(2)); - shape(5) = ( hp(0))*(1-hp(1))*( hp(2)); - shape(6) = ( hp(0))*( hp(1))*( hp(2)); - shape(7) = (1-hp(0))*( hp(1))*( hp(2)); - break; - } - } -} - - - -void Element :: GetShapeNew (const Point<3> & p, FlatVector & shape) const -{ - /* - if (shape.Size() < GetNP()) - { - cerr << "Element::GetShape: Length not fitting" << endl; - return; - } - */ - - switch (typ) - { - case TET: - { - shape(0) = p(0); - shape(1) = p(1); - shape(2) = p(2); - shape(3) = 1-p(0)-p(1)-p(2); - break; - } - - case TET10: - { - double lam1 = p(0); - double lam2 = p(1); - double lam3 = p(2); - double lam4 = 1-p(0)-p(1)-p(2); - - shape(0) = 2 * lam1 * (lam1-0.5); - shape(1) = 2 * lam2 * (lam2-0.5); - shape(2) = 2 * lam3 * (lam3-0.5); - shape(3) = 2 * lam4 * (lam4-0.5); - - shape(4) = 4 * lam1 * lam2; - shape(5) = 4 * lam1 * lam3; - shape(6) = 4 * lam1 * lam4; - shape(7) = 4 * lam2 * lam3; - shape(8) = 4 * lam2 * lam4; - shape(9) = 4 * lam3 * lam4; - - break; - } - - - case PYRAMID: - { - double noz = 1-p(2); - if (noz == 0.0) noz = 1e-10; - - double xi = p(0) / noz; - double eta = p(1) / noz; - shape(0) = (1-xi)*(1-eta) * (noz); - shape(1) = ( xi)*(1-eta) * (noz); - shape(2) = ( xi)*( eta) * (noz); - shape(3) = (1-xi)*( eta) * (noz); - shape(4) = p(2); - break; - } - - case PRISM: - { - shape(0) = p(0) * (1-p(2)); - shape(1) = p(1) * (1-p(2)); - shape(2) = (1-p(0)-p(1)) * (1-p(2)); - shape(3) = p(0) * p(2); - shape(4) = p(1) * p(2); - shape(5) = (1-p(0)-p(1)) * p(2); - break; - } - case HEX: - { - shape(0) = (1-p(0))*(1-p(1))*(1-p(2)); - shape(1) = ( p(0))*(1-p(1))*(1-p(2)); - shape(2) = ( p(0))*( p(1))*(1-p(2)); - shape(3) = (1-p(0))*( p(1))*(1-p(2)); - shape(4) = (1-p(0))*(1-p(1))*( p(2)); - shape(5) = ( p(0))*(1-p(1))*( p(2)); - shape(6) = ( p(0))*( p(1))*( p(2)); - shape(7) = (1-p(0))*( p(1))*( p(2)); - break; - } - } -} - - - - -void Element :: -GetDShape (const Point<3> & hp, DenseMatrix & dshape) const -{ - Point3d p = hp; - - int np = GetNP(); - if (dshape.Height() != 3 || dshape.Width() != np) - { - cerr << "Element::DShape: Sizes don't fit" << endl; - return; - } - - int i, j; - double eps = 1e-6; - Vector shaper(np), shapel(np); - - for (i = 1; i <= 3; i++) - { - Point3d pr(p), pl(p); - pr.X(i) += eps; - pl.X(i) -= eps; + Point3d pr(p), pl(p); + pr.X(i) += eps; + pl.X(i) -= eps; - GetShape (pr, shaper); - GetShape (pl, shapel); - for (j = 1; j <= np; j++) - dshape.Elem(i, j) = (shaper.Get(j) - shapel.Get(j)) / (2 * eps); - } -} - - - -void Element :: -GetDShapeNew (const Point<3> & p, MatrixFixWidth<3> & dshape) const -{ - switch (typ) - { - case TET: - { - dshape = 0; - dshape(0,0) = 1; - dshape(1,1) = 1; - dshape(2,2) = 1; - dshape(3,0) = -1; - dshape(3,1) = -1; - dshape(3,2) = -1; - break; + GetShape (pr, shaper); + GetShape (pl, shapel); + for (int j = 0; j < np; j++) + dshape(i-1, j) = (shaper(j) - shapel(j)) / (2 * eps); } - case PRISM: - { - dshape = 0; - dshape(0,0) = 1-p(2); - dshape(0,2) = -p(0); - dshape(1,1) = 1-p(2); - dshape(1,2) = -p(1); - dshape(2,0) = -(1-p(2)); - dshape(2,1) = -(1-p(2)); - dshape(2,2) = -(1-p(0)-p(1)); + } - dshape(3,0) = p(2); - dshape(3,2) = p(0); - dshape(4,1) = p(2); - dshape(4,2) = p(1); - dshape(5,0) = -p(2); - dshape(5,1) = -p(2); - dshape(5,2) = 1-p(0)-p(1); - break; - } - default: + + void Element :: + GetDShapeNew (const Point<3> & p, MatrixFixWidth<3> & dshape) const + { + switch (typ) { - int np = GetNP(); - double eps = 1e-6; - Vector shaper(np), shapel(np); + case TET: + { + dshape = 0; + dshape(0,0) = 1; + dshape(1,1) = 1; + dshape(2,2) = 1; + dshape(3,0) = -1; + dshape(3,1) = -1; + dshape(3,2) = -1; + break; + } + case PRISM: + { + dshape = 0; + dshape(0,0) = 1-p(2); + dshape(0,2) = -p(0); + dshape(1,1) = 1-p(2); + dshape(1,2) = -p(1); + dshape(2,0) = -(1-p(2)); + dshape(2,1) = -(1-p(2)); + dshape(2,2) = -(1-p(0)-p(1)); + + dshape(3,0) = p(2); + dshape(3,2) = p(0); + dshape(4,1) = p(2); + dshape(4,2) = p(1); + dshape(5,0) = -p(2); + dshape(5,1) = -p(2); + dshape(5,2) = 1-p(0)-p(1); + break; + } + + default: + { + int np = GetNP(); + double eps = 1e-6; + Vector shaper(np), shapel(np); - for (int i = 1; i <= 3; i++) - { - Point3d pr(p), pl(p); - pr.X(i) += eps; - pl.X(i) -= eps; + for (int i = 1; i <= 3; i++) + { + Point3d pr(p), pl(p); + pr.X(i) += eps; + pl.X(i) -= eps; - GetShapeNew (pr, shaper); - GetShapeNew (pl, shapel); - for (int j = 1; j <= np; j++) - dshape.Elem(j, i) = (shaper.Get(j) - shapel.Get(j)) / (2 * eps); - } + GetShapeNew (pr, shaper); + GetShapeNew (pl, shapel); + for (int j = 0; j < np; j++) + dshape(j, i-1) = (shaper(j) - shapel(j)) / (2 * eps); + } + } } - } -} + } -void Element :: -GetPointMatrix (const T_POINTS & points, - DenseMatrix & pmat) const -{ - int np = GetNP(); - /* - if (pmat.Width() != np || pmat.Height() != 3) - { + void Element :: + GetPointMatrix (const T_POINTS & points, + DenseMatrix & pmat) const + { + int np = GetNP(); + /* + if (pmat.Width() != np || pmat.Height() != 3) + { cerr << "Element::GetPointMatrix: sizes don't fit" << endl; return; + } + */ + for (int i = 1; i <= np; i++) + { + const Point3d & p = points.Get(PNum(i)); + pmat.Elem(1, i) = p.X(); + pmat.Elem(2, i) = p.Y(); + pmat.Elem(3, i) = p.Z(); + } + } + + + + + + + double Element :: CalcJacobianBadness (const T_POINTS & points) const + { + int i, j; + int nip = GetNIP(); + static DenseMatrix trans(3,3); + static DenseMatrix pmat; + + pmat.SetSize (3, GetNP()); + GetPointMatrix (points, pmat); + + double err = 0; + for (i = 1; i <= nip; i++) + { + GetTransformation (i, pmat, trans); + + // Frobenius norm + double frob = 0; + for (j = 1; j <= 9; j++) + frob += sqr (trans.Get(j)); + frob = sqrt (frob); + frob /= 3; + + double det = -trans.Det(); + + if (det <= 0) + err += 1e12; + else + err += frob * frob * frob / det; + } + + err /= nip; + return err; + } + + double Element :: + CalcJacobianBadnessDirDeriv (const T_POINTS & points, + int pi, Vec<3> & dir, double & dd) const + { + int i, j, k; + int nip = GetNIP(); + static DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); + static DenseMatrix pmat, vmat; + + pmat.SetSize (3, GetNP()); + vmat.SetSize (3, GetNP()); + + GetPointMatrix (points, pmat); + + for (i = 1; i <= np; i++) + for (j = 1; j <= 3; j++) + vmat.Elem(j, i) = 0; + for (j = 1; j <= 3; j++) + vmat.Elem(j, pi) = dir(j-1); + + + + double err = 0; + dd = 0; + + for (i = 1; i <= nip; i++) + { + GetTransformation (i, pmat, trans); + GetTransformation (i, vmat, dtrans); + + + // Frobenius norm + double frob = 0; + for (j = 1; j <= 9; j++) + frob += sqr (trans.Get(j)); + frob = sqrt (frob); + + double dfrob = 0; + for (j = 1; j <= 9; j++) + dfrob += trans.Get(j) * dtrans.Get(j); + dfrob = dfrob / frob; + + frob /= 3; + dfrob /= 3; + + + double det = trans.Det(); + double ddet = 0; + + for (j = 1; j <= 3; j++) + { + hmat = trans; + for (k = 1; k <= 3; k++) + hmat.Elem(k, j) = dtrans.Get(k, j); + ddet += hmat.Det(); + } + + + det *= -1; + ddet *= -1; + + + if (det <= 0) + err += 1e12; + else + { + err += frob * frob * frob / det; + dd += (3 * frob * frob * dfrob * det - frob * frob * frob * ddet) / (det * det); + } + } + + err /= nip; + dd /= nip; + return err; + } + + double Element :: + CalcJacobianBadnessGradient (const T_POINTS & points, + int pi, Vec<3> & grad) const + { + int nip = GetNIP(); + static DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); + static DenseMatrix pmat, vmat; + + pmat.SetSize (3, GetNP()); + vmat.SetSize (3, GetNP()); + + GetPointMatrix (points, pmat); + + for (int i = 1; i <= np; i++) + for (int j = 1; j <= 3; j++) + vmat.Elem(j, i) = 0; + for (int j = 1; j <= 3; j++) + vmat.Elem(j, pi) = 1.; + + + double err = 0; + + double dfrob[3]; + + grad = 0; + + for (int i = 1; i <= nip; i++) + { + GetTransformation (i, pmat, trans); + GetTransformation (i, vmat, dtrans); + + // Frobenius norm + double frob = 0; + for (int j = 1; j <= 9; j++) + frob += sqr (trans.Get(j)); + frob = sqrt (frob); + + for(int k = 0; k<3; k++) + { + dfrob[k] = 0; + for (int j = 1; j <= 3; j++) + dfrob[k] += trans.Get(k+1,j) * dtrans.Get(k+1,j); + dfrob[k] = dfrob[k] / (3.*frob); + } + + frob /= 3; + + double det = trans.Det(); + double ddet[3]; // = 0; + + for(int k=1; k<=3; k++) + { + int km1 = (k > 1) ? (k-1) : 3; + int kp1 = (k < 3) ? (k+1) : 1; + ddet[k-1] = 0; + for(int j=1; j<=3; j++) + { + int jm1 = (j > 1) ? (j-1) : 3; + int jp1 = (j < 3) ? (j+1) : 1; + + ddet[k-1] += (-1.)* dtrans.Get(k,j) * ( trans.Get(km1,jm1)*trans.Get(kp1,jp1) - + trans.Get(km1,jp1)*trans.Get(kp1,jm1) ); + } + } + + + det *= -1; + + if (det <= 0) + err += 1e12; + else + { + err += frob * frob * frob / det; + double fac = (frob * frob)/(det * det); + for(int j=0; j<3; j++) + grad(j) += fac * (3 * dfrob[j] * det - frob * ddet[j]); + } + } + + err /= nip; + grad *= 1./nip; + return err; + } + + + + + + + void Element :: ComputeIntegrationPointData () const + { + switch (GetType()) + { + case TET: if (ipdtet.Size()) return; break; + case TET10: if (ipdtet10.Size()) return; break; + default: + PrintSysError ("Element::ComputeIntegrationPoint, illegal type ", int(typ)); + } + + switch (GetType()) + { + case TET: ipdtet.SetSize(GetNIP()); break; + case TET10: ipdtet10.SetSize(GetNIP()); break; + default: + PrintSysError ("Element::ComputeIntegrationPoint, illegal type2 ", int(typ)); + } + + + for (int i = 1; i <= GetNIP(); i++) + { + IntegrationPointData * ipd = new IntegrationPointData; + GetIntegrationPoint (i, ipd->p, ipd->weight); + ipd->shape.SetSize(GetNP()); + ipd->dshape.SetSize(3, GetNP()); + + GetShape (ipd->p, ipd->shape); + GetDShape (ipd->p, ipd->dshape); + + switch (GetType()) + { + case TET: ipdtet.Elem(i).Reset(ipd); break; + case TET10: ipdtet10.Elem(i).Reset(ipd); break; + default: + PrintSysError ("Element::ComputeIntegrationPoint(2), illegal type ", int(typ)); + } + } + } + + + + + + + + + FaceDescriptor :: FaceDescriptor() + { + surfnr = domin = domout = bcprop = 0; + domin_singular = domout_singular = 0.; + tlosurf = -1; + bcname = 0; + firstelement = -1; + } + + FaceDescriptor :: FaceDescriptor(const FaceDescriptor& other) + : surfnr(other.surfnr), domin(other.domin), domout(other.domout), + tlosurf(other.tlosurf), bcprop(other.bcprop), bcname(other.bcname), + domin_singular(other.domin_singular), domout_singular(other.domout_singular) + { + firstelement = -1; + } + + FaceDescriptor :: + FaceDescriptor(int surfnri, int domini, int domouti, int tlosurfi) + { + surfnr = surfnri; + domin = domini; + domout = domouti; + tlosurf = tlosurfi; + bcprop = surfnri; + domin_singular = domout_singular = 0.; + bcname = 0; + firstelement = -1; + } + + FaceDescriptor :: FaceDescriptor(const Segment & seg) + { + surfnr = seg.si; + domin = seg.domin+1; + domout = seg.domout+1; + tlosurf = seg.tlosurf+1; + bcprop = 0; + domin_singular = domout_singular = 0.; + bcname = 0; + firstelement = -1; + } + + int FaceDescriptor :: SegmentFits (const Segment & seg) + { + return + surfnr == seg.si && + domin == seg.domin+1 && + domout == seg.domout+1 && + tlosurf == seg.tlosurf+1; + } + + + string FaceDescriptor :: GetBCName () const + { + if ( bcname ) + return *bcname; + else + return "default"; + + } + + /* + void FaceDescriptor :: SetBCName (string * bcn) + { + bcname = bcn; } */ - for (int i = 1; i <= np; i++) - { - const Point3d & p = points.Get(PNum(i)); - pmat.Elem(1, i) = p.X(); - pmat.Elem(2, i) = p.Y(); - pmat.Elem(3, i) = p.Z(); - } -} + + + ostream & operator<<(ostream & s, const FaceDescriptor & fd) + { + s << "surfnr = " << fd.surfnr + << ", domin = " << fd.domin + << ", domout = " << fd.domout + << ", tlosurf = " << fd.tlosurf + << ", bcprop = " << fd.bcprop + << ", domin_sing = " << fd.domin_singular + << ", domout_sing = " << fd.domout_singular; + return s; + } -double Element :: CalcJacobianBadness (const T_POINTS & points) const -{ - int i, j; - int nip = GetNIP(); - static DenseMatrix trans(3,3); - static DenseMatrix pmat; + Identifications :: Identifications (Mesh & amesh) + : mesh(amesh) + { + identifiedpoints = new INDEX_2_HASHTABLE(100); + identifiedpoints_nr = new INDEX_3_HASHTABLE(100); + maxidentnr = 0; + } + + Identifications :: ~Identifications () + { + delete identifiedpoints; + delete identifiedpoints_nr; + } + + void Identifications :: Delete () + { + delete identifiedpoints; + identifiedpoints = new INDEX_2_HASHTABLE(100); + delete identifiedpoints_nr; + identifiedpoints_nr = new INDEX_3_HASHTABLE(100); + maxidentnr = 0; + } + + void Identifications :: Add (PointIndex pi1, PointIndex pi2, int identnr) + { + // (*testout) << "Identification::Add, pi1 = " << pi1 << ", pi2 = " << pi2 << ", identnr = " << identnr << endl; + INDEX_2 pair (pi1, pi2); + identifiedpoints->Set (pair, identnr); + + INDEX_3 tripl (pi1, pi2, identnr); + identifiedpoints_nr->Set (tripl, 1); + + if (identnr > maxidentnr) maxidentnr = identnr; + + if (identnr+1 > idpoints_table.Size()) + idpoints_table.ChangeSize (identnr+1); + idpoints_table.Add (identnr, pair); - pmat.SetSize (3, GetNP()); - GetPointMatrix (points, pmat); + // timestamp = NextTimeStamp(); + } + + int Identifications :: Get (PointIndex pi1, PointIndex pi2) const + { + INDEX_2 pair(pi1, pi2); + if (identifiedpoints->Used (pair)) + return identifiedpoints->Get(pair); + else + return 0; + } + + bool Identifications :: Get (PointIndex pi1, PointIndex pi2, int nr) const + { + INDEX_3 tripl(pi1, pi2, nr); + if (identifiedpoints_nr->Used (tripl)) + return 1; + else + return 0; + } + + + + int Identifications :: GetSymmetric (PointIndex pi1, PointIndex pi2) const + { + INDEX_2 pair(pi1, pi2); + if (identifiedpoints->Used (pair)) + return identifiedpoints->Get(pair); + + pair = INDEX_2 (pi2, pi1); + if (identifiedpoints->Used (pair)) + return identifiedpoints->Get(pair); - double err = 0; - for (i = 1; i <= nip; i++) - { - GetTransformation (i, pmat, trans); - - // Frobenius norm - double frob = 0; - for (j = 1; j <= 9; j++) - frob += sqr (trans.Get(j)); - frob = sqrt (frob); - frob /= 3; - - double det = -trans.Det(); - - if (det <= 0) - err += 1e12; - else - err += frob * frob * frob / det; - } - - err /= nip; - return err; -} - -double Element :: -CalcJacobianBadnessDirDeriv (const T_POINTS & points, - int pi, Vec<3> & dir, double & dd) const -{ - int i, j, k, l; - int nip = GetNIP(); - static DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); - static DenseMatrix pmat, vmat; - - pmat.SetSize (3, GetNP()); - vmat.SetSize (3, GetNP()); - - GetPointMatrix (points, pmat); - - for (i = 1; i <= np; i++) - for (j = 1; j <= 3; j++) - vmat.Elem(j, i) = 0; - for (j = 1; j <= 3; j++) - vmat.Elem(j, pi) = dir(j-1); - - - - double err = 0; - dd = 0; - - for (i = 1; i <= nip; i++) - { - GetTransformation (i, pmat, trans); - GetTransformation (i, vmat, dtrans); - - - // Frobenius norm - double frob = 0; - for (j = 1; j <= 9; j++) - frob += sqr (trans.Get(j)); - frob = sqrt (frob); - - double dfrob = 0; - for (j = 1; j <= 9; j++) - dfrob += trans.Get(j) * dtrans.Get(j); - dfrob = dfrob / frob; - - frob /= 3; - dfrob /= 3; - - - double det = trans.Det(); - double ddet = 0; - - for (j = 1; j <= 3; j++) - { - hmat = trans; - for (k = 1; k <= 3; k++) - hmat.Elem(k, j) = dtrans.Get(k, j); - ddet += hmat.Det(); - } - - - det *= -1; - ddet *= -1; - - - if (det <= 0) - err += 1e12; - else - { - err += frob * frob * frob / det; - dd += (3 * frob * frob * dfrob * det - frob * frob * frob * ddet) / (det * det); - } - } - - err /= nip; - dd /= nip; - return err; -} - -double Element :: -CalcJacobianBadnessGradient (const T_POINTS & points, - int pi, Vec<3> & grad) const -{ - int i, j, k, l; - int nip = GetNIP(); - static DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); - static DenseMatrix pmat, vmat; - - pmat.SetSize (3, GetNP()); - vmat.SetSize (3, GetNP()); - - GetPointMatrix (points, pmat); - - for (i = 1; i <= np; i++) - for (j = 1; j <= 3; j++) - vmat.Elem(j, i) = 0; - for (j = 1; j <= 3; j++) - vmat.Elem(j, pi) = 1.; - - - double err = 0; - - double dfrob[3]; - - grad = 0; - - for (i = 1; i <= nip; i++) - { - GetTransformation (i, pmat, trans); - GetTransformation (i, vmat, dtrans); - - // Frobenius norm - double frob = 0; - for (j = 1; j <= 9; j++) - frob += sqr (trans.Get(j)); - frob = sqrt (frob); - - for(k = 0; k<3; k++) - { - dfrob[k] = 0; - for (j = 1; j <= 3; j++) - dfrob[k] += trans.Get(k+1,j) * dtrans.Get(k+1,j); - dfrob[k] = dfrob[k] / (3.*frob); - } - - frob /= 3; - - double det = trans.Det(); - double ddet[3]; // = 0; - - for(k=1; k<=3; k++) - { - int km1 = (k > 1) ? (k-1) : 3; - int kp1 = (k < 3) ? (k+1) : 1; - ddet[k-1] = 0; - for(j=1; j<=3; j++) - { - int jm1 = (j > 1) ? (j-1) : 3; - int jp1 = (j < 3) ? (j+1) : 1; - - ddet[k-1] += (-1.)* dtrans.Get(k,j) * ( trans.Get(km1,jm1)*trans.Get(kp1,jp1) - - trans.Get(km1,jp1)*trans.Get(kp1,jm1) ); - } - } - - - det *= -1; - - if (det <= 0) - err += 1e12; - else - { - err += frob * frob * frob / det; - double fac = (frob * frob)/(det * det); - for(j=0; j<3; j++) - grad(j) += fac * (3 * dfrob[j] * det - frob * ddet[j]); - } - } - - err /= nip; - grad *= 1./nip; - return err; -} - - - - - - -void Element :: ComputeIntegrationPointData () const -{ - switch (GetType()) - { - case TET: if (ipdtet.Size()) return; break; - case TET10: if (ipdtet10.Size()) return; break; - default: - PrintSysError ("Element::ComputeIntegrationPoint, illegal type ", int(typ)); - } - - switch (GetType()) - { - case TET: ipdtet.SetSize(GetNIP()); break; - case TET10: ipdtet10.SetSize(GetNIP()); break; - default: - PrintSysError ("Element::ComputeIntegrationPoint, illegal type2 ", int(typ)); - } - - - for (int i = 1; i <= GetNIP(); i++) - { - IntegrationPointData * ipd = new IntegrationPointData; - GetIntegrationPoint (i, ipd->p, ipd->weight); - ipd->shape.SetSize(GetNP()); - ipd->dshape.SetSize(3, GetNP()); - - GetShape (ipd->p, ipd->shape); - GetDShape (ipd->p, ipd->dshape); - - switch (GetType()) - { - case TET: ipdtet.Elem(i).Reset(ipd); break; - case TET10: ipdtet10.Elem(i).Reset(ipd); break; - default: - PrintSysError ("Element::ComputeIntegrationPoint(2), illegal type ", int(typ)); - } - } -} - - - - - - - - -FaceDescriptor :: FaceDescriptor() -{ - surfnr = domin = domout = bcprop = 0; - domin_singular = domout_singular = 0.; - tlosurf = -1; - bcname = 0; - firstelement = -1; -} - -FaceDescriptor :: FaceDescriptor(const FaceDescriptor& other) - : surfnr(other.surfnr), domin(other.domin), domout(other.domout), - tlosurf(other.tlosurf), bcprop(other.bcprop), bcname(other.bcname), - domin_singular(other.domin_singular), domout_singular(other.domout_singular) -{ - firstelement = -1; -} - -FaceDescriptor :: -FaceDescriptor(int surfnri, int domini, int domouti, int tlosurfi) -{ - surfnr = surfnri; - domin = domini; - domout = domouti; - tlosurf = tlosurfi; - bcprop = surfnri; - domin_singular = domout_singular = 0.; - bcname = 0; - firstelement = -1; -} - -FaceDescriptor :: FaceDescriptor(const Segment & seg) -{ - surfnr = seg.si; - domin = seg.domin+1; - domout = seg.domout+1; - tlosurf = seg.tlosurf+1; - bcprop = 0; - domin_singular = domout_singular = 0.; - bcname = 0; - firstelement = -1; -} - -int FaceDescriptor :: SegmentFits (const Segment & seg) -{ - return - surfnr == seg.si && - domin == seg.domin+1 && - domout == seg.domout+1 && - tlosurf == seg.tlosurf+1; -} - - -string FaceDescriptor :: GetBCName () const -{ - if ( bcname ) - return *bcname; - else - return "default"; - -} - -/* -void FaceDescriptor :: SetBCName (string * bcn) -{ - bcname = bcn; -} -*/ - - -ostream & operator<<(ostream & s, const FaceDescriptor & fd) -{ - s << "surfnr = " << fd.surfnr - << ", domin = " << fd.domin - << ", domout = " << fd.domout - << ", tlosurf = " << fd.tlosurf - << ", bcprop = " << fd.bcprop - << ", domin_sing = " << fd.domin_singular - << ", domout_sing = " << fd.domout_singular; - return s; -} - - - - - - -Identifications :: Identifications (Mesh & amesh) - : mesh(amesh) -{ - identifiedpoints = new INDEX_2_HASHTABLE(100); - identifiedpoints_nr = new INDEX_3_HASHTABLE(100); - maxidentnr = 0; -} - -Identifications :: ~Identifications () -{ - delete identifiedpoints; - delete identifiedpoints_nr; -} - -void Identifications :: Delete () -{ - delete identifiedpoints; - identifiedpoints = new INDEX_2_HASHTABLE(100); - delete identifiedpoints_nr; - identifiedpoints_nr = new INDEX_3_HASHTABLE(100); - maxidentnr = 0; -} - -void Identifications :: Add (PointIndex pi1, PointIndex pi2, int identnr) -{ - // (*testout) << "Identification::Add, pi1 = " << pi1 << ", pi2 = " << pi2 << ", identnr = " << identnr << endl; - INDEX_2 pair (pi1, pi2); - identifiedpoints->Set (pair, identnr); - - INDEX_3 tripl (pi1, pi2, identnr); - identifiedpoints_nr->Set (tripl, 1); - - if (identnr > maxidentnr) maxidentnr = identnr; - - if (identnr+1 > idpoints_table.Size()) - idpoints_table.ChangeSize (identnr+1); - idpoints_table.Add (identnr, pair); - - // timestamp = NextTimeStamp(); -} - -int Identifications :: Get (PointIndex pi1, PointIndex pi2) const -{ - INDEX_2 pair(pi1, pi2); - if (identifiedpoints->Used (pair)) - return identifiedpoints->Get(pair); - else return 0; -} - -bool Identifications :: Get (PointIndex pi1, PointIndex pi2, int nr) const -{ - INDEX_3 tripl(pi1, pi2, nr); - if (identifiedpoints_nr->Used (tripl)) - return 1; - else - return 0; -} + } + void Identifications :: GetMap (int identnr, Array & identmap, bool symmetric) const + { + identmap.SetSize (mesh.GetNP()); + identmap = 0; -int Identifications :: GetSymmetric (PointIndex pi1, PointIndex pi2) const -{ - INDEX_2 pair(pi1, pi2); - if (identifiedpoints->Used (pair)) - return identifiedpoints->Get(pair); + if (identnr) + for (int i = 0; i < idpoints_table[identnr].Size(); i++) + { + INDEX_2 pair = idpoints_table[identnr][i]; + identmap[pair.I1()] = pair.I2(); + if(symmetric) + identmap[pair.I2()] = pair.I1(); + } - pair = INDEX_2 (pi2, pi1); - if (identifiedpoints->Used (pair)) - return identifiedpoints->Get(pair); - - return 0; -} - - -void Identifications :: GetMap (int identnr, Array & identmap, bool symmetric) const -{ - identmap.SetSize (mesh.GetNP()); - identmap = 0; - - if (identnr) - for (int i = 0; i < idpoints_table[identnr].Size(); i++) + else { - INDEX_2 pair = idpoints_table[identnr][i]; - identmap[pair.I1()] = pair.I2(); - if(symmetric) - identmap[pair.I2()] = pair.I1(); + cout << "getmap, identnr = " << identnr << endl; + + for (int i = 1; i <= identifiedpoints_nr->GetNBags(); i++) + for (int j = 1; j <= identifiedpoints_nr->GetBagSize(i); j++) + { + INDEX_3 i3; + int dummy; + identifiedpoints_nr->GetData (i, j, i3, dummy); + + if (i3.I3() == identnr || !identnr) + { + identmap.Elem(i3.I1()) = i3.I2(); + if(symmetric) + identmap.Elem(i3.I2()) = i3.I1(); + } + } } - else - { - cout << "getmap, identnr = " << identnr << endl; - - for (int i = 1; i <= identifiedpoints_nr->GetNBags(); i++) - for (int j = 1; j <= identifiedpoints_nr->GetBagSize(i); j++) - { - INDEX_3 i3; - int dummy; - identifiedpoints_nr->GetData (i, j, i3, dummy); - - if (i3.I3() == identnr || !identnr) - { - identmap.Elem(i3.I1()) = i3.I2(); - if(symmetric) - identmap.Elem(i3.I2()) = i3.I1(); - } - } - } - -} + } -void Identifications :: GetPairs (int identnr, - Array & identpairs) const -{ - identpairs.SetSize(0); + void Identifications :: GetPairs (int identnr, + Array & identpairs) const + { + identpairs.SetSize(0); - if (identnr == 0) + if (identnr == 0) + for (int i = 1; i <= identifiedpoints->GetNBags(); i++) + for (int j = 1; j <= identifiedpoints->GetBagSize(i); j++) + { + INDEX_2 i2; + int nr; + identifiedpoints->GetData (i, j, i2, nr); + identpairs.Append (i2); + } + else + for (int i = 1; i <= identifiedpoints_nr->GetNBags(); i++) + for (int j = 1; j <= identifiedpoints_nr->GetBagSize(i); j++) + { + INDEX_3 i3; + int dummy; + identifiedpoints_nr->GetData (i, j, i3 , dummy); + + if (i3.I3() == identnr) + identpairs.Append (INDEX_2(i3.I1(), i3.I2())); + } + } + + + void Identifications :: SetMaxPointNr (int maxpnum) + { for (int i = 1; i <= identifiedpoints->GetNBags(); i++) for (int j = 1; j <= identifiedpoints->GetBagSize(i); j++) - { - INDEX_2 i2; - int nr; - identifiedpoints->GetData (i, j, i2, nr); - identpairs.Append (i2); - } - else - for (int i = 1; i <= identifiedpoints_nr->GetNBags(); i++) - for (int j = 1; j <= identifiedpoints_nr->GetBagSize(i); j++) - { - INDEX_3 i3; - int dummy; - identifiedpoints_nr->GetData (i, j, i3 , dummy); - - if (i3.I3() == identnr) - identpairs.Append (INDEX_2(i3.I1(), i3.I2())); - } -} - - -void Identifications :: SetMaxPointNr (int maxpnum) -{ - for (int i = 1; i <= identifiedpoints->GetNBags(); i++) - for (int j = 1; j <= identifiedpoints->GetBagSize(i); j++) - { - INDEX_2 i2; - int nr; - identifiedpoints->GetData (i, j, i2, nr); + { + INDEX_2 i2; + int nr; + identifiedpoints->GetData (i, j, i2, nr); - if (i2.I1() > maxpnum || i2.I2() > maxpnum) - { - i2.I1() = i2.I2() = -1; - identifiedpoints->SetData (i, j, i2, -1); - } - } -} + if (i2.I1() > maxpnum || i2.I2() > maxpnum) + { + i2.I1() = i2.I2() = -1; + identifiedpoints->SetData (i, j, i2, -1); + } + } + } -void Identifications :: Print (ostream & ost) const -{ - ost << "Identifications:" << endl; - ost << "pairs: " << endl << *identifiedpoints << endl; - ost << "pairs and nr: " << endl << *identifiedpoints_nr << endl; - ost << "table: " << endl << idpoints_table << endl; -} + void Identifications :: Print (ostream & ost) const + { + ost << "Identifications:" << endl; + ost << "pairs: " << endl << *identifiedpoints << endl; + ost << "pairs and nr: " << endl << *identifiedpoints_nr << endl; + ost << "table: " << endl << idpoints_table << endl; + } -MeshingParameters :: MeshingParameters () -{ - optimize3d = "cmdmustm"; - //optimize3d = "cmdmstm"; - optsteps3d = 3; - optimize2d = "smsmsmSmSmSm"; - optsteps2d = 3; - opterrpow = 2; - blockfill = 1; - filldist = 0.1; - safety = 5; - relinnersafety = 3; - uselocalh = 1; - grading = 0.3; - delaunay = 1; - maxh = 1e10; - minh = 0; - meshsizefilename = NULL; - startinsurface = 0; - checkoverlap = 1; - checkoverlappingboundary = 1; - checkchartboundary = 1; - curvaturesafety = 2; - segmentsperedge = 1; - parthread = 0; + MeshingParameters :: MeshingParameters () + { + optimize3d = "cmdmustm"; + //optimize3d = "cmdmstm"; + optsteps3d = 3; + optimize2d = "smsmsmSmSmSm"; + optsteps2d = 3; + opterrpow = 2; + blockfill = 1; + filldist = 0.1; + safety = 5; + relinnersafety = 3; + uselocalh = 1; + grading = 0.3; + delaunay = 1; + maxh = 1e10; + minh = 0; + meshsizefilename = NULL; + startinsurface = 0; + checkoverlap = 1; + checkoverlappingboundary = 1; + checkchartboundary = 1; + curvaturesafety = 2; + segmentsperedge = 1; + parthread = 0; - elsizeweight = 0.2; - giveuptol2d = 200; - giveuptol = 10; - maxoutersteps = 10; - starshapeclass = 5; - baseelnp = 0; - sloppy = 1; + elsizeweight = 0.2; + giveuptol2d = 200; + giveuptol = 10; + maxoutersteps = 10; + starshapeclass = 5; + baseelnp = 0; + sloppy = 1; - badellimit = 175; - check_impossible = 0; - secondorder = 0; -} + badellimit = 175; + check_impossible = 0; + secondorder = 0; + } -void MeshingParameters :: Print (ostream & ost) const -{ - ost << "Meshing parameters: " << endl - << "optimize3d = " << optimize3d << endl - << "optsteps3d = " << optsteps3d << endl - << " optimize2d = " << optimize2d << endl - << " optsteps2d = " << optsteps2d << endl - << " opterrpow = " << opterrpow << endl - << " blockfill = " << blockfill << endl - << " filldist = " << filldist << endl - << " safety = " << safety << endl - << " relinnersafety = " << relinnersafety << endl - << " uselocalh = " << uselocalh << endl - << " grading = " << grading << endl - << " delaunay = " << delaunay << endl - << " maxh = " << maxh << endl; - if(meshsizefilename) - ost << " meshsizefilename = " << meshsizefilename << endl; - else - ost << " meshsizefilename = NULL" << endl; - ost << " startinsurface = " << startinsurface << endl - << " checkoverlap = " << checkoverlap << endl - << " checkchartboundary = " << checkchartboundary << endl - << " curvaturesafety = " << curvaturesafety << endl - << " segmentsperedge = " << segmentsperedge << endl - << " parthread = " << parthread << endl - << " elsizeweight = " << elsizeweight << endl - << " giveuptol2d = " << giveuptol2d << endl - << " giveuptol = " << giveuptol << endl - << " maxoutersteps = " << maxoutersteps << endl - << " starshapeclass = " << starshapeclass << endl - << " baseelnp = " << baseelnp << endl - << " sloppy = " << sloppy << endl - << " badellimit = " << badellimit << endl - << " secondorder = " << secondorder << endl - << " elementorder = " << elementorder << endl - << " quad = " << quad << endl - << " inverttets = " << inverttets << endl - << " inverttrigs = " << inverttrigs << endl; -} + void MeshingParameters :: Print (ostream & ost) const + { + ost << "Meshing parameters: " << endl + << "optimize3d = " << optimize3d << endl + << "optsteps3d = " << optsteps3d << endl + << " optimize2d = " << optimize2d << endl + << " optsteps2d = " << optsteps2d << endl + << " opterrpow = " << opterrpow << endl + << " blockfill = " << blockfill << endl + << " filldist = " << filldist << endl + << " safety = " << safety << endl + << " relinnersafety = " << relinnersafety << endl + << " uselocalh = " << uselocalh << endl + << " grading = " << grading << endl + << " delaunay = " << delaunay << endl + << " maxh = " << maxh << endl; + if(meshsizefilename) + ost << " meshsizefilename = " << meshsizefilename << endl; + else + ost << " meshsizefilename = NULL" << endl; + ost << " startinsurface = " << startinsurface << endl + << " checkoverlap = " << checkoverlap << endl + << " checkchartboundary = " << checkchartboundary << endl + << " curvaturesafety = " << curvaturesafety << endl + << " segmentsperedge = " << segmentsperedge << endl + << " parthread = " << parthread << endl + << " elsizeweight = " << elsizeweight << endl + << " giveuptol2d = " << giveuptol2d << endl + << " giveuptol = " << giveuptol << endl + << " maxoutersteps = " << maxoutersteps << endl + << " starshapeclass = " << starshapeclass << endl + << " baseelnp = " << baseelnp << endl + << " sloppy = " << sloppy << endl + << " badellimit = " << badellimit << endl + << " secondorder = " << secondorder << endl + << " elementorder = " << elementorder << endl + << " quad = " << quad << endl + << " inverttets = " << inverttets << endl + << " inverttrigs = " << inverttrigs << endl; + } -void MeshingParameters :: CopyFrom(const MeshingParameters & other) -{ - //strcpy(optimize3d,other.optimize3d); - optimize3d = other.optimize3d; - optsteps3d = other.optsteps3d; - //strcpy(optimize2d,other.optimize2d); - optimize2d = other.optimize2d; - optsteps2d = other.optsteps2d; - opterrpow = other.opterrpow; - blockfill = other.blockfill; - filldist = other.filldist; - safety = other.safety; - relinnersafety = other.relinnersafety; - uselocalh = other.uselocalh; - grading = other.grading; - delaunay = other.delaunay; - maxh = other.maxh; - //strcpy(const_cast(meshsizefilename), other.meshsizefilename); - //const_cast(meshsizefilename) = other.meshsizefilename; //??? - startinsurface = other.startinsurface; - checkoverlap = other.checkoverlap; - checkoverlappingboundary = other.checkoverlappingboundary; - checkchartboundary = other.checkchartboundary; - curvaturesafety = other.curvaturesafety; - segmentsperedge = other.segmentsperedge; - parthread = other.parthread; - elsizeweight = other.elsizeweight; - giveuptol2d = other.giveuptol2d; - giveuptol = other.giveuptol; - maxoutersteps = other.maxoutersteps; - starshapeclass = other.starshapeclass; - baseelnp = other.baseelnp; - sloppy = other.sloppy; - badellimit = other.badellimit; - secondorder = other.secondorder; - elementorder = other.elementorder; - quad = other.quad; - inverttets = other.inverttets; - inverttrigs = other.inverttrigs; -} + void MeshingParameters :: CopyFrom(const MeshingParameters & other) + { + //strcpy(optimize3d,other.optimize3d); + optimize3d = other.optimize3d; + optsteps3d = other.optsteps3d; + //strcpy(optimize2d,other.optimize2d); + optimize2d = other.optimize2d; + optsteps2d = other.optsteps2d; + opterrpow = other.opterrpow; + blockfill = other.blockfill; + filldist = other.filldist; + safety = other.safety; + relinnersafety = other.relinnersafety; + uselocalh = other.uselocalh; + grading = other.grading; + delaunay = other.delaunay; + maxh = other.maxh; + //strcpy(const_cast(meshsizefilename), other.meshsizefilename); + //const_cast(meshsizefilename) = other.meshsizefilename; //??? + startinsurface = other.startinsurface; + checkoverlap = other.checkoverlap; + checkoverlappingboundary = other.checkoverlappingboundary; + checkchartboundary = other.checkchartboundary; + curvaturesafety = other.curvaturesafety; + segmentsperedge = other.segmentsperedge; + parthread = other.parthread; + elsizeweight = other.elsizeweight; + giveuptol2d = other.giveuptol2d; + giveuptol = other.giveuptol; + maxoutersteps = other.maxoutersteps; + starshapeclass = other.starshapeclass; + baseelnp = other.baseelnp; + sloppy = other.sloppy; + badellimit = other.badellimit; + secondorder = other.secondorder; + elementorder = other.elementorder; + quad = other.quad; + inverttets = other.inverttets; + inverttrigs = other.inverttrigs; + } -DebugParameters :: DebugParameters () -{ - slowchecks = 0; - haltsuccess = 0; - haltnosuccess = 0; - haltlargequalclass = 0; - haltsegment = 0; - haltsegmentp1 = 0; - haltsegmentp2 = 0; -}; + DebugParameters :: DebugParameters () + { + slowchecks = 0; + haltsuccess = 0; + haltnosuccess = 0; + haltlargequalclass = 0; + haltsegment = 0; + haltsegmentp1 = 0; + haltsegmentp2 = 0; + }; diff --git a/libsrc/meshing/netrule2.cpp b/libsrc/meshing/netrule2.cpp index 02e19746..d760d0e9 100644 --- a/libsrc/meshing/netrule2.cpp +++ b/libsrc/meshing/netrule2.cpp @@ -128,10 +128,10 @@ int netrule :: IsLineInFreeZone2 (const Point2d & p1, const Point2d & p2) const { int left, right, allleft, allright; - if (p1.X() > fzmaxx && p2.X() > fzmaxx || - p1.X() < fzminx && p2.X() < fzminx || - p1.Y() > fzmaxy && p2.Y() > fzmaxy || - p1.Y() < fzminy && p2.Y() < fzminy) return 0; + if ( (p1.X() > fzmaxx && p2.X() > fzmaxx) || + (p1.X() < fzminx && p2.X() < fzminx) || + (p1.Y() > fzmaxy && p2.Y() > fzmaxy) || + (p1.Y() < fzminy && p2.Y() < fzminy) ) return 0; for (int i = 1; i <= transfreezone.Size(); i++) { diff --git a/libsrc/meshing/netrule3.cpp b/libsrc/meshing/netrule3.cpp index 50c44ae8..03771409 100644 --- a/libsrc/meshing/netrule3.cpp +++ b/libsrc/meshing/netrule3.cpp @@ -61,7 +61,7 @@ void vnetrule :: SetFreeZoneTransformation (const Vector & allp, int tolclass) for (i = 1; i <= 3; i++) { 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); oldutofreezonelimit->Mult (vp, vfp2); @@ -70,7 +70,7 @@ void vnetrule :: SetFreeZoneTransformation (const Vector & allp, int tolclass) vfp1.Add (lam2, vfp2); for (j = 1; j <= nfp; j++) - transfreezone.Elem(j).X(i) = vfp1.Elem(j); + transfreezone.Elem(j).X(i) = vfp1(j-1); } // MARK(setfz2); @@ -709,7 +709,7 @@ int vnetrule :: IsTriangleInFreeSet (const Point3d & p1, const Point3d & p2, minit = 1000; fold = 1E10; - + while (1) { @@ -814,7 +814,7 @@ int vnetrule :: IsTriangleInFreeSet (const Point3d & p1, const Point3d & p2, if (it > minit) (*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) { diff --git a/libsrc/meshing/parser3.cpp b/libsrc/meshing/parser3.cpp index 15284012..c6ccba4a 100644 --- a/libsrc/meshing/parser3.cpp +++ b/libsrc/meshing/parser3.cpp @@ -706,13 +706,13 @@ void vnetrule :: LoadRule (istream & ist) if (quality < 100) { - for (i = 1; i <= 3; i++) + for (int i = 1; i <= 3; i++) { - for (j = 1; j <= points.Size(); j++) - vp.Elem(j) = points.Get(j).X(i); + for (int j = 1; j <= points.Size(); j++) + vp(j-1) = points.Get(j).X(i); oldutofreezone->Mult(vp, vfp); - for (j = 1; j <= freezone.Size(); j++) - freezone.Elem(j).X(i) = vfp.Get(j); + for (int j = 1; j <= freezone.Size(); j++) + freezone.Elem(j).X(i) = vfp(j-1); } // for (i = 1; i <= freezone.Size(); i++) // (*testout) << "freepoint: " << freezone.Get(i) << endl; diff --git a/libsrc/meshing/ruler2.cpp b/libsrc/meshing/ruler2.cpp index 6b52a7cf..f65e2444 100644 --- a/libsrc/meshing/ruler2.cpp +++ b/libsrc/meshing/ruler2.cpp @@ -397,11 +397,11 @@ int Meshing2 ::ApplyRules (Array & lpoints, { 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))); - oldu.Set (2*i-1, ui.X()); - oldu.Set (2*i , ui.Y()); + oldu (2*i-2) = ui.X(); + oldu (2*i-1) = ui.Y(); } rule -> SetFreeZoneTransformation (oldu, tolerance); @@ -513,8 +513,8 @@ int Meshing2 ::ApplyRules (Array & lpoints, for (i = oldnp + 1; i <= rule->GetNP(); i++) { np = rule->GetPoint(i); - np.X() += newu.Elem (2 * (i-oldnp) - 1); - np.Y() += newu.Elem (2 * (i-oldnp)); + np.X() += newu (2 * (i-oldnp) - 2); + np.Y() += newu (2 * (i-oldnp) - 1); pmap.Elem(i) = lpoints.Append (np); } diff --git a/libsrc/meshing/ruler3.cpp b/libsrc/meshing/ruler3.cpp index 647cf15b..69a63914 100644 --- a/libsrc/meshing/ruler3.cpp +++ b/libsrc/meshing/ruler3.cpp @@ -603,13 +603,13 @@ int Meshing3 :: ApplyRules { const Point3d & lp = lpoints.Get(pmap.Get(i)); const Point3d & rp = rule->GetPoint(i); - oldu.Set (3*i-2, lp.X()-rp.X()); - oldu.Set (3*i-1, lp.Y()-rp.Y()); - oldu.Set (3*i , lp.Z()-rp.Z()); + oldu (3*i-3) = lp.X()-rp.X(); + oldu (3*i-2) = lp.Y()-rp.Y(); + oldu (3*i-1) = lp.Z()-rp.Z(); - allp.Set (3*i-2, lp.X()); - allp.Set (3*i-1, lp.Y()); - allp.Set (3*i , lp.Z()); + allp (3*i-3) = lp.X(); + allp (3*i-2) = lp.Y(); + allp (3*i-1) = lp.Z(); } if (rule->GetNP() > rule->GetNOldP()) @@ -623,9 +623,9 @@ int Meshing3 :: ApplyRules for (i = rule->GetNOldP()+1; i <= rule->GetNP(); i++) { const Point3d & rp = rule->GetPoint(i); - allp.Set (3*i-2, rp.X() + newu.Get(3*i-2 - idiff)); - allp.Set (3*i-1, rp.Y() + newu.Get(3*i-1 - idiff)); - allp.Set (3*i , rp.Z() + newu.Get(3*i - idiff)); + allp (3*i-3) = rp.X() + newu(3*i-3 - idiff); + allp (3*i-2) = rp.Y() + newu(3*i-2 - idiff); + allp (3*i-1) = rp.Z() + newu(3*i-1 - idiff); } rule->SetFreeZoneTransformation (allp, @@ -865,9 +865,9 @@ int Meshing3 :: ApplyRules for (i = oldnp + 1; i <= rule->GetNP(); i++) { np = rule->GetPoint(i); - np.X() += newu.Elem (3 * (i-oldnp) - 2); - np.Y() += newu.Elem (3 * (i-oldnp) - 1); - np.Z() += newu.Elem (3 * (i-oldnp)); + np.X() += newu (3 * (i-oldnp) - 3); + np.Y() += newu (3 * (i-oldnp) - 2); + np.Z() += newu (3 * (i-oldnp) - 1); pmap.Elem(i) = lpoints.Append (np); } diff --git a/libsrc/meshing/smoothing2.5.cpp b/libsrc/meshing/smoothing2.5.cpp index 52e94617..c9de8e21 100644 --- a/libsrc/meshing/smoothing2.5.cpp +++ b/libsrc/meshing/smoothing2.5.cpp @@ -219,9 +219,9 @@ namespace netgen while (loci <= 5 && !moveisok) { loci ++; - mesh[pi](0) = origp(0) + x.Get(1)*fact; - mesh[pi](1) = origp(1) + x.Get(2)*fact; - mesh[pi](2) = origp(2) + x.Get(3)*fact; + mesh[pi](0) = origp(0) + x(0)*fact; + mesh[pi](1) = origp(1) + x(1)*fact; + mesh[pi](2) = origp(2) + x(2)*fact; fact = fact/2.; diff --git a/libsrc/meshing/smoothing2.cpp b/libsrc/meshing/smoothing2.cpp index 1395d046..1794b840 100644 --- a/libsrc/meshing/smoothing2.cpp +++ b/libsrc/meshing/smoothing2.cpp @@ -441,7 +441,7 @@ namespace netgen vgrad = 0.0; badness = 0; - pp1 = sp1 + x.Get(1) * t1; + pp1 = sp1 + x(0) * t1; meshthis -> ProjectPoint2 (surfi, surfi2, pp1); for (j = 0; j < locelements.Size(); j++) @@ -540,7 +540,7 @@ namespace netgen pts2d.Elem(pi) = Point2d (t1 * (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++) @@ -552,8 +552,8 @@ namespace netgen hbad = bel. CalcJacobianBadnessDirDeriv (pts2d, lpi, vdir, hderiv); - - grad.Elem(k) += hderiv; + + grad(k-1) += hderiv; if (k == 1) badness += hbad; } @@ -590,7 +590,7 @@ namespace netgen // pp1 = sp1; // 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 pts2d; pts2d.SetSize(mesh.GetNP()); @@ -611,7 +611,7 @@ namespace netgen pts2d.Elem(pi) = Point2d (t1 * (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)); diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 7a77de4e..49131e44 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -149,15 +149,15 @@ namespace netgen static double eps = 1e-6; 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); - hx.Elem(i) = x.Get(i) - eps * h; + hx(i) = x(i) - eps * h; 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); @@ -246,17 +246,17 @@ namespace netgen static Vector res; res.SetSize (m.Height()); - for (i = 1;i <= 3; i++) - hv.Elem(i) = vp.Get(i); - hv.Elem(4) = 1; + for (i = 0;i < 3; i++) + hv(i) = vp(i); + hv(3) = 1; m.Mult (hv, res); for (i = 1; i <= res.Size(); i++) { - if (res.Get(i) < 1e-10) + if (res(i-1) < 1e-10) badness += 1e24; else - badness += 1 / res.Get(i); + badness += 1 / res(i-1); } return badness; @@ -269,15 +269,15 @@ namespace netgen static double eps = 1e-6; 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); - hx.Elem(i) = x.Get(i) - eps * h; + hx(i) = x(i) - eps * h; 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); @@ -374,42 +374,12 @@ namespace netgen double PointFunction :: PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const { - double f;//, delta = h * 1e-6; - // Point<3> hpp; + double f = 0; + // 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; Vec<3> vgradi, vgrad(0,0,0); - // badness = 0; - hp = points[actpind]; points[actpind] = Point<3> (pp); @@ -421,10 +391,10 @@ namespace netgen for (int k = 1; k <= 4; k++) if (el.PNum(k) == actpind) { - CalcTetBadnessGrad (points[el.PNum(1)], - points[el.PNum(2)], - points[el.PNum(3)], - points[el.PNum(4)], -1, k, vgradi); + f += CalcTetBadnessGrad (points[el.PNum(1)], + points[el.PNum(2)], + points[el.PNum(3)], + points[el.PNum(4)], -1, k, vgradi); vgrad += vgradi; } @@ -606,19 +576,19 @@ namespace netgen static Vector di; int n = m.Height(); - p4.Elem(1) = pp(0); - p4.Elem(2) = pp(1); - p4.Elem(3) = pp(2); - p4.Elem(4) = 1; + p4(0) = pp(0); + p4(1) = pp(1); + p4(2) = pp(2); + p4(3) = 1; di.SetSize (n); m.Mult (p4, di); double sum = 0; - for (int i = 1; i <= n; i++) + for (int i = 0; i < n; i++) { - if (di.Get(i) > 0) - sum += 1 / di.Get(i); + if (di(i) > 0) + sum += 1 / di(i); else return 1e16; } @@ -635,25 +605,25 @@ namespace netgen int n = m.Height(); - p4.Elem(1) = pp(0); - p4.Elem(2) = pp(1); - p4.Elem(3) = pp(2); - p4.Elem(4) = 1; + p4(0) = pp(0); + p4(1) = pp(1); + p4(2) = pp(2); + p4(3) = 1; di.SetSize (n); m.Mult (p4, di); double sum = 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; - grad(0) -= idi * idi * m.Get(i, 1); - grad(1) -= idi * idi * m.Get(i, 2); - grad(2) -= idi * idi * m.Get(i, 3); + grad(0) -= idi * idi * m(i, 0); + grad(1) -= idi * idi * m(i, 1); + grad(2) -= idi * idi * m(i, 2); } else { @@ -774,9 +744,9 @@ namespace netgen } hx = x; - hx.Elem(i) = x.Get(i) + eps; + hx(i-1) = x(i-1) + eps; f11 = Func(hx); - hx.Elem(i) = x.Get(i) - eps; + hx(i-1) = x(i-1) - eps; f22 = Func(hx); 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; double CalcTotalBad (const Mesh::T_POINTS & points, @@ -1100,10 +1059,10 @@ double JacobianPointFunction :: Func (const Vector & v) const 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) - 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++) @@ -1129,10 +1088,10 @@ FuncGrad (const Vector & x, Vector & g) const double badness = 0;//, hbad; 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) - 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; //Vec3d vdir; @@ -1154,7 +1113,7 @@ FuncGrad (const Vector & x, Vector & g) const CalcJacobianBadnessGradient (points, lpi, hderiv); for(k=0; k<3; k++) - g.Elem(k+1) += hderiv(k); + g(k) += hderiv(k); /* for (k = 1; k <= 3; k++) @@ -1174,10 +1133,10 @@ FuncGrad (const Vector & x, Vector & g) const if(onplane) { - double scal = nv(0)*g.Get(1) + nv(1)*g.Get(2) + nv(2)*g.Get(3); - g.Elem(1) -= scal*nv(0); - g.Elem(2) -= scal*nv(1); - g.Elem(3) -= scal*nv(2); + double scal = nv(0)*g(0) + nv(1)*g(1) + nv(2)*g(2); + g(0) -= scal*nv(0); + g(1) -= scal*nv(1); + g(2) -= scal*nv(2); } //(*testout) << "g = " << g << endl; @@ -1197,14 +1156,14 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const double badness = 0; 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) - 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; 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) { @@ -1569,9 +1528,9 @@ void Mesh :: ImproveMesh (OPTIMIZEGOAL goal) //*testout << "start BFGS, pok" << endl; BFGS (x, freeminf, par); //*testout << "BFGS complete, pok" << endl; - points[i](0) += x.Get(1); - points[i](1) += x.Get(2); - points[i](2) += x.Get(3); + points[i](0) += x(0); + points[i](1) += x(1); + points[i](2) += x(2); } } PrintDot ('\n'); @@ -1694,9 +1653,9 @@ void Mesh :: ImproveMeshJacobian (OPTIMIZEGOAL goal, const BitArray * usepoint) //*testout << "start BFGS, Jacobian" << endl; BFGS (x, pf, par); //*testout << "end BFGS, Jacobian" << endl; - points.Elem(i)(0) += x.Get(1); - points.Elem(i)(1) += x.Get(2); - points.Elem(i)(2) += x.Get(3); + points.Elem(i)(0) += x(0); + points.Elem(i)(1) += x(1); + points.Elem(i)(2) += x(2); } else { @@ -1873,12 +1832,12 @@ void Mesh :: ImproveMeshJacobianOnSurface (const BitArray & usepoint, BFGS (x, pf_sum, par); - for(j=1; j<=3; j++) - points.Elem(i)(j-1) += x.Get(j);// - scal*nv[i-1].X(j); + for(j=0; j<3; j++) + points.Elem(i)(j) += x(j);// - scal*nv[i-1].X(j); if(brother != -1) - for(j=1; j<=3; j++) - points.Elem(brother)(j-1) += x.Get(j);// - scal*nv[brother-1].X(j); + for(j=0; j<3; j++) + points.Elem(brother)(j) += x(j);// - scal*nv[brother-1].X(j); } diff --git a/libsrc/stlgeom/stlgeom.cpp b/libsrc/stlgeom/stlgeom.cpp index 2e055c7d..228ff5b4 100644 --- a/libsrc/stlgeom/stlgeom.cpp +++ b/libsrc/stlgeom/stlgeom.cpp @@ -183,7 +183,7 @@ void STLGeometry :: SmoothNormals() for (k = 0; k < 3; k++) - hv.Elem(k+1) = ngeom(k); + hv(k) = ngeom(k); hm.Mult (hv, hv2); /* @@ -230,7 +230,7 @@ void STLGeometry :: SmoothNormals() } for (k = 0; k < 3; k++) - hv.Elem(k+1) = nnb(k); + hv(k) = nnb(k); hm.Mult (hv, hv2); /* @@ -245,7 +245,7 @@ void STLGeometry :: SmoothNormals() } 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); GetTriangle(i).SetNormal(newn); @@ -1562,9 +1562,9 @@ void STLGeometry :: NeighbourAnglesOfSelectedTrig() for (i = 1; i <= NONeighbourTrigs(st); i++) { 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), - 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"); - //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 Array edgecnt; @@ -2676,7 +2676,7 @@ void STLGeometry :: LinkEdges() //worked edges Array 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 Vec3d v1,v2; diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index 1d4b26f7..13ced0ef 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -3,7 +3,6 @@ #include - #include #include #include @@ -983,9 +982,9 @@ namespace netgen for (int i = 1; i <= locms.Size(); i++) { Point3d p = mesh->Point(i); - locms.Elem(i) = mesh->GetH (p); - if (locms.Elem(i) > maxh) maxh = locms.Elem(i); - if (locms.Elem(i) < minh) minh = locms.Elem(i); + locms(i-1) = mesh->GetH (p); + if (locms(i-1) > maxh) maxh = locms(i-1); + if (locms(i-1) < minh) minh = locms(i-1); } if (!locms.Size()) { minh = 1; maxh = 10; } @@ -1145,11 +1144,11 @@ namespace netgen if (vispar.colormeshsize) { - SetOpenGlColor (locms.Get(el[0]), minh, maxh, 1); + SetOpenGlColor (locms(el[0]-1), minh, maxh, 1); glVertex3dv (lp0); - SetOpenGlColor (locms.Get(el[1]), minh, maxh, 1); + SetOpenGlColor (locms(el[1]-1), minh, maxh, 1); glVertex3dv (lp1); - SetOpenGlColor (locms.Get(el[2]), minh, maxh, 1); + SetOpenGlColor (locms(el[2]-1), minh, maxh, 1); glVertex3dv (lp2); } else diff --git a/ng/nginterface.cpp b/ng/nginterface.cpp index 0b81fbed..c086788b 100644 --- a/ng/nginterface.cpp +++ b/ng/nginterface.cpp @@ -149,10 +149,8 @@ using namespace netgen; 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 ()); geometry2d.Reset (); @@ -171,6 +169,8 @@ void Ng_LoadGeometry (const char * filename) if (strcmp(filename,"")==0) return; + PrintMessage (1, "Load geometry from file ", filename); + ifstream infile (filename); if ((strcmp (&filename[strlen(filename)-3], "geo") == 0) || @@ -1936,8 +1936,6 @@ void Ng_SetSolutionData (Ng_SolutionData * soldata) // vssolution.ClearSolutionData (); VisualSceneSolution::SolData * vss = new VisualSceneSolution::SolData; - cout << "Add solution " << soldata->name << ", type = " << soldata->soltype << endl; - vss->name = new char[strlen (soldata->name)+1]; strcpy (vss->name, soldata->name); diff --git a/nglib/Makefile.am b/nglib/Makefile.am index 002916fd..c46af5a7 100644 --- a/nglib/Makefile.am +++ b/nglib/Makefile.am @@ -1,5 +1,7 @@ include_HEADERS = nglib.h +dist_pkgdata_DATA = cube.surf + AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include $(OCCFLAGS) $(MPI_INCLUDES) lib_LTLIBRARIES = libnglib.la diff --git a/nglib/cube.surf b/nglib/cube.surf new file mode 100644 index 00000000..1fee858b --- /dev/null +++ b/nglib/cube.surf @@ -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