From 50b81a747fb0ff57cce994c323584dd3809d9c19 Mon Sep 17 00:00:00 2001 From: william Date: Thu, 28 Aug 2014 10:54:00 +0100 Subject: [PATCH 01/15] BUG: mantis #1373: eigen value/vector calculation now handles repeated eigen values --- applications/test/tensor/Test-tensor.C | 116 +++++ .../primitives/Tensor/tensor/tensor.C | 417 ++++++------------ 2 files changed, 252 insertions(+), 281 deletions(-) diff --git a/applications/test/tensor/Test-tensor.C b/applications/test/tensor/Test-tensor.C index 6da37d33..e406ab7f 100644 --- a/applications/test/tensor/Test-tensor.C +++ b/applications/test/tensor/Test-tensor.C @@ -69,5 +69,121 @@ int main() Info<< (symm(t7) && t7) - (0.5*(t7 + t7.T()) && t7) << endl; Info<< (t7 && symm(t7)) - (t7 && 0.5*(t7 + t7.T())) << endl; + + /* + // Lots of awkward eigenvector tests ... + + tensor T_rand_real + ( + 0.9999996423721313, 0.3330855667591095, 0.6646450161933899, + 0.9745196104049683, 0.0369445420801640, 0.0846728682518005, + 0.6474838852882385, 0.1617118716239929, 0.2041363865137100 + ); + Debug(T_rand_real); + vector L_rand_real(eigenValues(T_rand_real)); + Debug(L_rand_real); + tensor U_rand_real(eigenVectors(T_rand_real)); + Debug(U_rand_real); + + Info << endl << endl; + + tensor T_rand_imag + ( + 0.8668024539947510, 0.1664607226848602, 0.8925783634185791, + 0.9126510620117188, 0.7408077120780945, 0.1499115079641342, + 0.0936608463525772, 0.7615650296211243, 0.8953040242195129 + ); + Debug(T_rand_imag); + vector L_rand_imag(eigenValues(T_rand_imag)); + Debug(L_rand_imag); + tensor U_rand_imag(eigenVectors(T_rand_imag)); + Debug(U_rand_imag); + + Info << endl << endl; + + tensor T_rand_symm + ( + 1.9999992847442627, 1.3076051771640778, 1.3121289014816284, + 1.3076051771640778, 0.0738890841603279, 0.2463847398757935, + 1.3121289014816284, 0.2463847398757935, 0.4082727730274200 + ); + Debug(T_rand_symm); + vector L_rand_symm(eigenValues(T_rand_symm)); + Debug(L_rand_symm); + tensor U_rand_symm(eigenVectors(T_rand_symm)); + Debug(U_rand_symm); + + Info << endl << endl; + + symmTensor T_rand_Symm + ( + 1.9999992847442627, 1.3076051771640778, 1.3121289014816284, + 0.0738890841603279, 0.2463847398757935, + 0.4082727730274200 + ); + Debug(T_rand_Symm); + vector L_rand_Symm(eigenValues(T_rand_Symm)); + Debug(L_rand_Symm); + tensor U_rand_Symm(eigenVectors(T_rand_Symm)); + Debug(U_rand_Symm); + + Info << endl << endl; + + tensor T_rand_diag + ( + 0.8668024539947510, 0, 0, + 0, 0.7408077120780945, 0, + 0, 0, 0.8953040242195129 + ); + Debug(T_rand_diag); + vector L_rand_diag(eigenValues(T_rand_diag)); + Debug(L_rand_diag); + tensor U_rand_diag(eigenVectors(T_rand_diag)); + Debug(U_rand_diag); + + Info << endl << endl; + + tensor T_repeated + ( + 0, 1, 1, + 1, 0, 1, + 1, 1, 0 + ); + Debug(T_repeated); + vector L_repeated(eigenValues(T_repeated)); + Debug(L_repeated); + tensor U_repeated(eigenVectors(T_repeated)); + Debug(U_repeated); + + Info << endl << endl; + + tensor T_repeated_zero + ( + 1, 1, 1, + 1, 1, 1, + 1, 1, 1 + ); + Debug(T_repeated_zero); + vector L_repeated_zero(eigenValues(T_repeated_zero)); + Debug(L_repeated_zero); + tensor U_repeated_zero(eigenVectors(T_repeated_zero)); + Debug(U_repeated_zero); + + Info << endl << endl; + + tensor T_triple + ( + 2, 0, 0, + 0, 2, 0, + 0, 0, 2 + ); + Debug(T_triple); + vector L_triple(eigenValues(T_triple)); + Debug(L_triple); + tensor U_triple(eigenVectors(T_triple)); + Debug(U_triple); + */ + + return 0; } diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C index 7e68450a..30f5ee24 100644 --- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C +++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C @@ -89,99 +89,80 @@ namespace Foam Foam::vector Foam::eigenValues(const tensor& t) { - scalar i = 0; - scalar ii = 0; - scalar iii = 0; + // The eigenvalues + scalar i, ii, iii; - if - ( - ( - mag(t.xy()) + mag(t.xz()) + mag(t.yx()) - + mag(t.yz()) + mag(t.zx()) + mag(t.zy()) - ) - < SMALL - ) + // Coefficients of the characteristic polynmial + // x^3 + a*x^2 + b*x + c = 0 + scalar a = + - t.xx() - t.yy() - t.zz(); + + scalar b = + t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz() + - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz(); + + scalar c = + - t.xx()*t.yy()*t.zz() + - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx() + + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx(); + + // Auxillary variables + scalar aBy3 = a/3; + + scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3 + scalar PPP = P*P*P; + + scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2 + scalar QQ = Q*Q; + + // Three identical roots + if (mag(P) < SMALL && mag(Q) < SMALL) { - // diagonal matrix - i = t.xx(); - ii = t.yy(); - iii = t.zz(); + return vector(- aBy3, - aBy3, - aBy3); } + + // Two identical roots and one distinct root + else if (mag(PPP/QQ - 1) < SMALL) + { + scalar sqrtP = sqrt(P); + scalar signQ = sign(Q); + + i = ii = signQ*sqrtP - aBy3; + iii = - 2*signQ*sqrtP - aBy3; + } + + // Three distinct roots + else if (PPP > QQ) + { + scalar sqrtP = sqrt(P); + scalar value = cos(acos(Q/sqrt(PPP))/3); + scalar delta = sqrt(3 - 3*value*value); + + i = - 2*sqrtP*value - aBy3; + ii = sqrtP*(value + delta) - aBy3; + iii = sqrtP*(value - delta) - aBy3; + } + + // One real root, two imaginary roots + // based on the above logic, PPP must be less than QQ else { - scalar a = -t.xx() - t.yy() - t.zz(); + WarningIn("eigenValues(const tensor&)") + << "complex eigenvalues detected for tensor: " << t + << endl; - scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz() - - t.xy()*t.yx() - t.xz()*t.zx() - t.yz()*t.zy(); - - scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.zx() - - t.xz()*t.yx()*t.zy() + t.xz()*t.yy()*t.zx() - + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy(); - - // If there is a zero root - if (mag(c) < 1e-100) + if (mag(P) < SMALL) { - scalar disc = sqr(a) - 4*b; - - if (disc >= -SMALL) - { - scalar q = -0.5*sqrt(max(0.0, disc)); - - i = 0; - ii = -0.5*a + q; - iii = -0.5*a - q; - } - else - { - FatalErrorIn("eigenValues(const tensor&)") - << "zero and complex eigenvalues in tensor: " << t - << abort(FatalError); - } + i = cbrt(QQ/2); } else { - scalar Q = (a*a - 3*b)/9; - scalar R = (2*a*a*a - 9*a*b + 27*c)/54; - - scalar R2 = sqr(R); - scalar Q3 = pow3(Q); - - // Three different real roots - if (R2 < Q3) - { - scalar sqrtQ = sqrt(Q); - scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ)))); - - scalar m2SqrtQ = -2*sqrtQ; - scalar aBy3 = a/3; - - i = m2SqrtQ*cos(theta/3) - aBy3; - ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3; - iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3; - } - else - { - scalar A = cbrt(R + sqrt(R2 - Q3)); - - // Three equal real roots - if (A < SMALL) - { - scalar root = -a/3; - return vector(root, root, root); - } - else - { - // Complex roots - WarningIn("eigenValues(const tensor&)") - << "complex eigenvalues detected for tensor: " << t - << endl; - - return vector::zero; - } - } + scalar w = cbrt(- Q - sqrt(QQ - PPP)); + i = w + P/w - aBy3; } - } + return vector(-VGREAT, i, VGREAT); + } // Sort the eigenvalues into ascending order if (i > ii) @@ -203,24 +184,35 @@ Foam::vector Foam::eigenValues(const tensor& t) } -Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda) +Foam::vector Foam::eigenVector +( + const tensor& t, + const scalar lambda +) { - if (mag(lambda) < SMALL) - { - return vector::zero; - } + // Constantly rotating direction ensures different eigenvectors are + // generated when called sequentially with a multiple eigenvalue + static vector direction(0,0,1); + vector oldDirection(direction); + scalar temp = direction[2]; + direction[2] = direction[1]; + direction[1] = direction[0]; + direction[0] = temp; - // Construct the matrix for the eigenvector problem + // Construct the linear system for this eigenvalue tensor A(t - lambda*I); - // Calculate the sub-determinants of the 3 components - scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy(); - scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx(); - scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx(); + // Determinants of the 2x2 sub-matrices used to find the eigenvectors + scalar sd0, sd1, sd2; + scalar magSd0, magSd1, magSd2; - scalar magSd0 = mag(sd0); - scalar magSd1 = mag(sd1); - scalar magSd2 = mag(sd2); + // Sub-determinants for a unique eivenvalue + sd0 = A.yy()*A.zz() - A.yz()*A.zy(); + sd1 = A.zz()*A.xx() - A.zx()*A.xz(); + sd2 = A.xx()*A.yy() - A.xy()*A.yx(); + magSd0 = mag(sd0); + magSd1 = mag(sd1); + magSd2 = mag(sd2); // Evaluate the eigenvector using the largest sub-determinant if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL) @@ -231,9 +223,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda) (A.yz()*A.zx() - A.zz()*A.yx())/sd0, (A.zy()*A.yx() - A.yy()*A.zx())/sd0 ); - ev /= mag(ev); - return ev; + return ev/mag(ev); } else if (magSd1 >= magSd2 && magSd1 > SMALL) { @@ -243,9 +234,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda) 1, (A.zx()*A.xy() - A.xx()*A.zy())/sd1 ); - ev /= mag(ev); - return ev; + return ev/mag(ev); } else if (magSd2 > SMALL) { @@ -255,14 +245,55 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda) (A.yx()*A.xz() - A.xx()*A.yz())/sd2, 1 ); - ev /= mag(ev); - return ev; + return ev/mag(ev); } - else + + // Sub-determinants for a repeated eigenvalue + sd0 = A.yy()*direction.z() - A.yz()*direction.y(); + sd1 = A.zz()*direction.x() - A.zx()*direction.z(); + sd2 = A.xx()*direction.y() - A.xy()*direction.x(); + magSd0 = mag(sd0); + magSd1 = mag(sd1); + magSd2 = mag(sd2); + + // Evaluate the eigenvector using the largest sub-determinant + if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL) { - return vector::zero; + vector ev + ( + 1, + (A.yz()*direction.x() - direction.z()*A.yx())/sd0, + (direction.y()*A.yx() - A.yy()*direction.x())/sd0 + ); + + return ev/mag(ev); } + else if (magSd1 >= magSd2 && magSd1 > SMALL) + { + vector ev + ( + (direction.z()*A.zy() - A.zz()*direction.y())/sd1, + 1, + (A.zx()*direction.y() - direction.x()*A.zy())/sd1 + ); + + return ev/mag(ev); + } + else if (magSd2 > SMALL) + { + vector ev + ( + (A.xy()*direction.z() - direction.y()*A.xz())/sd2, + (direction.x()*A.xz() - A.xx()*direction.z())/sd2, + 1 + ); + + return ev/mag(ev); + } + + // Triple eigenvalue + return oldDirection; } @@ -283,195 +314,19 @@ Foam::tensor Foam::eigenVectors(const tensor& t) Foam::vector Foam::eigenValues(const symmTensor& t) { - scalar i = 0; - scalar ii = 0; - scalar iii = 0; - - if - ( - ( - mag(t.xy()) + mag(t.xz()) + mag(t.xy()) - + mag(t.yz()) + mag(t.xz()) + mag(t.yz()) - ) - < SMALL - ) - { - // diagonal matrix - i = t.xx(); - ii = t.yy(); - iii = t.zz(); - } - else - { - scalar a = -t.xx() - t.yy() - t.zz(); - - scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz() - - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz(); - - scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz() - - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz() - + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz(); - - // If there is a zero root - if (mag(c) < 1e-100) - { - scalar disc = sqr(a) - 4*b; - - if (disc >= -SMALL) - { - scalar q = -0.5*sqrt(max(0.0, disc)); - - i = 0; - ii = -0.5*a + q; - iii = -0.5*a - q; - } - else - { - FatalErrorIn("eigenValues(const tensor&)") - << "zero and complex eigenvalues in tensor: " << t - << abort(FatalError); - } - } - else - { - scalar Q = (a*a - 3*b)/9; - scalar R = (2*a*a*a - 9*a*b + 27*c)/54; - - scalar R2 = sqr(R); - scalar Q3 = pow3(Q); - - // Three different real roots - if (R2 < Q3) - { - scalar sqrtQ = sqrt(Q); - scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ)))); - - scalar m2SqrtQ = -2*sqrtQ; - scalar aBy3 = a/3; - - i = m2SqrtQ*cos(theta/3) - aBy3; - ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3; - iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3; - } - else - { - scalar A = cbrt(R + sqrt(R2 - Q3)); - - // Three equal real roots - if (A < SMALL) - { - scalar root = -a/3; - return vector(root, root, root); - } - else - { - // Complex roots - WarningIn("eigenValues(const symmTensor&)") - << "complex eigenvalues detected for symmTensor: " << t - << endl; - - return vector::zero; - } - } - } - } - - - // Sort the eigenvalues into ascending order - if (i > ii) - { - Swap(i, ii); - } - - if (ii > iii) - { - Swap(ii, iii); - } - - if (i > ii) - { - Swap(i, ii); - } - - return vector(i, ii, iii); + return eigenValues(tensor(t)); } Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda) { - if (mag(lambda) < SMALL) - { - return vector::zero; - } - - // Construct the matrix for the eigenvector problem - symmTensor A(t - lambda*I); - - // Calculate the sub-determinants of the 3 components - scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz(); - scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz(); - scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy(); - - scalar magSd0 = mag(sd0); - scalar magSd1 = mag(sd1); - scalar magSd2 = mag(sd2); - - // Evaluate the eigenvector using the largest sub-determinant - if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL) - { - vector ev - ( - 1, - (A.yz()*A.xz() - A.zz()*A.xy())/sd0, - (A.yz()*A.xy() - A.yy()*A.xz())/sd0 - ); - ev /= mag(ev); - - return ev; - } - else if (magSd1 >= magSd2 && magSd1 > SMALL) - { - vector ev - ( - (A.xz()*A.yz() - A.zz()*A.xy())/sd1, - 1, - (A.xz()*A.xy() - A.xx()*A.yz())/sd1 - ); - ev /= mag(ev); - - return ev; - } - else if (magSd2 > SMALL) - { - vector ev - ( - (A.xy()*A.yz() - A.yy()*A.xz())/sd2, - (A.xy()*A.xz() - A.xx()*A.yz())/sd2, - 1 - ); - ev /= mag(ev); - - return ev; - } - else - { - return vector::zero; - } + return eigenVector(tensor(t), lambda); } Foam::tensor Foam::eigenVectors(const symmTensor& t) { - vector evals(eigenValues(t)); - - tensor evs - ( - eigenVector(t, evals.x()), - eigenVector(t, evals.y()), - eigenVector(t, evals.z()) - ); - - return evs; + return eigenVectors(tensor(t)); } From 68fcd7ad038ef62c7dfd7a15bca8078785248241 Mon Sep 17 00:00:00 2001 From: william Date: Thu, 28 Aug 2014 11:25:43 +0100 Subject: [PATCH 02/15] BUG: mantis #1373: quick return for diagonal matrices --- .../primitives/Tensor/tensor/tensor.C | 139 ++++++++++-------- 1 file changed, 79 insertions(+), 60 deletions(-) diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C index 30f5ee24..fa20ef6b 100644 --- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C +++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C @@ -92,76 +92,95 @@ Foam::vector Foam::eigenValues(const tensor& t) // The eigenvalues scalar i, ii, iii; - // Coefficients of the characteristic polynmial - // x^3 + a*x^2 + b*x + c = 0 - scalar a = - - t.xx() - t.yy() - t.zz(); - - scalar b = - t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz() - - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz(); - - scalar c = - - t.xx()*t.yy()*t.zz() - - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx() - + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx(); - - // Auxillary variables - scalar aBy3 = a/3; - - scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3 - scalar PPP = P*P*P; - - scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2 - scalar QQ = Q*Q; - - // Three identical roots - if (mag(P) < SMALL && mag(Q) < SMALL) + // diagonal matrix + if + ( + ( + mag(t.xy()) + mag(t.xz()) + mag(t.yx()) + + mag(t.yz()) + mag(t.zx()) + mag(t.zy()) + ) + < SMALL + ) { - return vector(- aBy3, - aBy3, - aBy3); + i = t.xx(); + ii = t.yy(); + iii = t.zz(); } - // Two identical roots and one distinct root - else if (mag(PPP/QQ - 1) < SMALL) - { - scalar sqrtP = sqrt(P); - scalar signQ = sign(Q); - - i = ii = signQ*sqrtP - aBy3; - iii = - 2*signQ*sqrtP - aBy3; - } - - // Three distinct roots - else if (PPP > QQ) - { - scalar sqrtP = sqrt(P); - scalar value = cos(acos(Q/sqrt(PPP))/3); - scalar delta = sqrt(3 - 3*value*value); - - i = - 2*sqrtP*value - aBy3; - ii = sqrtP*(value + delta) - aBy3; - iii = sqrtP*(value - delta) - aBy3; - } - - // One real root, two imaginary roots - // based on the above logic, PPP must be less than QQ + // non-diagonal matrix else { - WarningIn("eigenValues(const tensor&)") - << "complex eigenvalues detected for tensor: " << t - << endl; + // Coefficients of the characteristic polynmial + // x^3 + a*x^2 + b*x + c = 0 + scalar a = + - t.xx() - t.yy() - t.zz(); - if (mag(P) < SMALL) + scalar b = + t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz() + - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz(); + + scalar c = + - t.xx()*t.yy()*t.zz() + - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx() + + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx(); + + // Auxillary variables + scalar aBy3 = a/3; + + scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3 + scalar PPP = P*P*P; + + scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2 + scalar QQ = Q*Q; + + // Three identical roots + if (mag(P) < SMALL && mag(Q) < SMALL) { - i = cbrt(QQ/2); + return vector(- aBy3, - aBy3, - aBy3); } + + // Two identical roots and one distinct root + else if (mag(PPP/QQ - 1) < SMALL) + { + scalar sqrtP = sqrt(P); + scalar signQ = sign(Q); + + i = ii = signQ*sqrtP - aBy3; + iii = - 2*signQ*sqrtP - aBy3; + } + + // Three distinct roots + else if (PPP > QQ) + { + scalar sqrtP = sqrt(P); + scalar value = cos(acos(Q/sqrt(PPP))/3); + scalar delta = sqrt(3 - 3*value*value); + + i = - 2*sqrtP*value - aBy3; + ii = sqrtP*(value + delta) - aBy3; + iii = sqrtP*(value - delta) - aBy3; + } + + // One real root, two imaginary roots + // based on the above logic, PPP must be less than QQ else { - scalar w = cbrt(- Q - sqrt(QQ - PPP)); - i = w + P/w - aBy3; - } + WarningIn("eigenValues(const tensor&)") + << "complex eigenvalues detected for tensor: " << t + << endl; - return vector(-VGREAT, i, VGREAT); + if (mag(P) < SMALL) + { + i = cbrt(QQ/2); + } + else + { + scalar w = cbrt(- Q - sqrt(QQ - PPP)); + i = w + P/w - aBy3; + } + + return vector(-VGREAT, i, VGREAT); + } } // Sort the eigenvalues into ascending order @@ -192,7 +211,7 @@ Foam::vector Foam::eigenVector { // Constantly rotating direction ensures different eigenvectors are // generated when called sequentially with a multiple eigenvalue - static vector direction(0,0,1); + static vector direction(1,0,0); vector oldDirection(direction); scalar temp = direction[2]; direction[2] = direction[1]; From 5dc2851823511e5146d5c11edab6501e99e2b1e9 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 28 Aug 2014 15:32:38 +0100 Subject: [PATCH 03/15] ENH: Ensight: updated support for generic tensor --- src/sampling/sampledSurface/writers/ensight/ensightPTraits.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C b/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C index 8d3ed3be..1511e636 100644 --- a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C +++ b/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -40,7 +40,7 @@ const char* const Foam::ensightPTraits::typeName = "tensor symm"; const char* const Foam::ensightPTraits::typeName = - Foam::pTraits::typeName; + "tensor asym"; // ************************************************************************* // From ad188431c463e355675f2d45b99a1deec3c815a8 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 29 Aug 2014 10:37:26 +0100 Subject: [PATCH 04/15] ENH: polyMesh: avoid virtual mechanism --- src/OpenFOAM/meshes/polyMesh/polyMesh.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index 668a5953..945e727b 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -1179,7 +1179,9 @@ Foam::tmp Foam::polyMesh::movePoints if (debug && moveError) { - write(); + // Write mesh to ease debugging. Note we want to avoid calling + // e.g. fvMesh::write since meshPhi not yet complete. + polyMesh::write(); } return sweptVols; From a94eb42d39a7eaee697077a6c287afb3f23c5fd0 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 3 Sep 2014 11:45:53 +0100 Subject: [PATCH 05/15] ENH: pointToPointPlanarInterpolation: additional nearest only interpolation --- src/meshTools/Make/options | 2 + .../pointToPointPlanarInterpolation.C | 249 +++++++++++------- .../pointToPointPlanarInterpolation.H | 12 +- 3 files changed, 166 insertions(+), 97 deletions(-) diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options index 1b22dab2..1713152e 100644 --- a/src/meshTools/Make/options +++ b/src/meshTools/Make/options @@ -1,7 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/triSurface/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ -I$(LIB_SRC)/fileFormats/lnInclude LIB_LIBS = \ -ltriSurface \ + -lsurfMesh \ -lfileFormats diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C index 1f7de7c5..5d69ba1f 100644 --- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -29,8 +29,9 @@ License #include "vector2D.H" #include "triSurface.H" #include "triSurfaceTools.H" -#include "OFstream.H" +#include "OBJstream.H" #include "Time.H" +#include "matchPoints.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -139,114 +140,172 @@ void Foam::pointToPointPlanarInterpolation::calcWeights const pointField& destPoints ) { - tmp tlocalVertices - ( - referenceCS_.localPosition(sourcePoints) - ); - vectorField& localVertices = tlocalVertices(); - - const boundBox bb(localVertices, true); - const point bbMid(bb.midpoint()); - - if (debug) + if (nearestOnly_) { - Info<< "pointToPointPlanarInterpolation::readData :" - << " Perturbing points with " << perturb_ - << " fraction of a random position inside " << bb - << " to break any ties on regular meshes." - << nl << endl; - } - - Random rndGen(123456); - forAll(localVertices, i) - { - localVertices[i] += - perturb_ - *(rndGen.position(bb.min(), bb.max())-bbMid); - } - - // Determine triangulation - List localVertices2D(localVertices.size()); - forAll(localVertices, i) - { - localVertices2D[i][0] = localVertices[i][0]; - localVertices2D[i][1] = localVertices[i][1]; - } - - triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); - - tmp tlocalFaceCentres - ( - referenceCS_.localPosition + labelList destToSource; + bool fullMatch = matchPoints ( - destPoints - ) - ); - const pointField& localFaceCentres = tlocalFaceCentres(); + destPoints, + sourcePoints, + scalarField(destPoints.size(), GREAT), + true, // verbose + destToSource + ); - if (debug) - { - Pout<< "pointToPointPlanarInterpolation::readData :" - <<" Dumping triangulated surface to triangulation.stl" << endl; - s.write("triangulation.stl"); - - OFstream str("localFaceCentres.obj"); - Pout<< "readSamplePoints :" - << " Dumping face centres to " << str.name() << endl; - - forAll(localFaceCentres, i) + if (!fullMatch) { - const point& p = localFaceCentres[i]; - str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl; + FatalErrorIn("pointToPointPlanarInterpolation::calcWeights(..)") + << "Did not find a corresponding sourcePoint for every face" + << " centre" << exit(FatalError); + } + + nearestVertex_.setSize(destPoints.size()); + nearestVertexWeight_.setSize(destPoints.size()); + forAll(nearestVertex_, i) + { + nearestVertex_[i][0] = destToSource[i]; + nearestVertex_[i][1] = -1; + nearestVertex_[i][2] = -1; + + nearestVertexWeight_[i][0] = 1.0; + nearestVertexWeight_[i][1] = 0.0; + nearestVertexWeight_[i][2] = 0.0; + } + + if (debug) + { + forAll(destPoints, i) + { + label v0 = nearestVertex_[i][0]; + + Pout<< "For location " << destPoints[i] + << " sampling vertex " << v0 + << " at:" << sourcePoints[v0] + << " distance:" << mag(sourcePoints[v0]-destPoints[i]) + << endl; + } + + OBJstream str("destToSource.obj"); + Pout<< "pointToPointPlanarInterpolation::calcWeights :" + << " Dumping lines from face centres to original points to " + << str.name() << endl; + + forAll(destPoints, i) + { + label v0 = nearestVertex_[i][0]; + str.write(linePointRef(destPoints[i], sourcePoints[v0])); + } } } - - // Determine interpolation onto face centres. - triSurfaceTools::calcInterpolationWeights - ( - s, - localFaceCentres, // points to interpolate to - nearestVertex_, - nearestVertexWeight_ - ); - - if (debug) + else { - forAll(sourcePoints, i) + tmp tlocalVertices + ( + referenceCS_.localPosition(sourcePoints) + ); + vectorField& localVertices = tlocalVertices(); + + const boundBox bb(localVertices, true); + const point bbMid(bb.midpoint()); + + if (debug) { - Pout<< "source:" << i << " at:" << sourcePoints[i] - << " 2d:" << localVertices[i] - << endl; + Info<< "pointToPointPlanarInterpolation::calcWeights :" + << " Perturbing points with " << perturb_ + << " fraction of a random position inside " << bb + << " to break any ties on regular meshes." + << nl << endl; } - - forAll(destPoints, i) + Random rndGen(123456); + forAll(localVertices, i) { - label v0 = nearestVertex_[i][0]; - label v1 = nearestVertex_[i][1]; - label v2 = nearestVertex_[i][2]; + localVertices[i] += + perturb_ + *(rndGen.position(bb.min(), bb.max())-bbMid); + } - Pout<< "For location " << destPoints[i] - << " 2d:" << localFaceCentres[i] - << " sampling vertices" << nl - << " " << v0 - << " at:" << sourcePoints[v0] - << " weight:" << nearestVertexWeight_[i][0] << nl; + // Determine triangulation + List localVertices2D(localVertices.size()); + forAll(localVertices, i) + { + localVertices2D[i][0] = localVertices[i][0]; + localVertices2D[i][1] = localVertices[i][1]; + } - if (v1 != -1) + triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); + + tmp tlocalFaceCentres + ( + referenceCS_.localPosition + ( + destPoints + ) + ); + const pointField& localFaceCentres = tlocalFaceCentres(); + + if (debug) + { + Pout<< "pointToPointPlanarInterpolation::calcWeights :" + <<" Dumping triangulated surface to triangulation.stl" << endl; + s.write("triangulation.stl"); + + OBJstream str("localFaceCentres.obj"); + Pout<< "pointToPointPlanarInterpolation::calcWeights :" + << " Dumping face centres to " << str.name() << endl; + + forAll(localFaceCentres, i) { - Pout<< " " << v1 - << " at:" << sourcePoints[v1] - << " weight:" << nearestVertexWeight_[i][1] << nl; + str.write(localFaceCentres[i]); } - if (v2 != -1) + } + + // Determine interpolation onto face centres. + triSurfaceTools::calcInterpolationWeights + ( + s, + localFaceCentres, // points to interpolate to + nearestVertex_, + nearestVertexWeight_ + ); + + if (debug) + { + forAll(sourcePoints, i) { - Pout<< " " << v2 - << " at:" << sourcePoints[v2] - << " weight:" << nearestVertexWeight_[i][2] << nl; + Pout<< "source:" << i << " at:" << sourcePoints[i] + << " 2d:" << localVertices[i] + << endl; } - Pout<< endl; + forAll(destPoints, i) + { + label v0 = nearestVertex_[i][0]; + label v1 = nearestVertex_[i][1]; + label v2 = nearestVertex_[i][2]; + + Pout<< "For location " << destPoints[i] + << " 2d:" << localFaceCentres[i] + << " sampling vertices" << nl + << " " << v0 + << " at:" << sourcePoints[v0] + << " weight:" << nearestVertexWeight_[i][0] << nl; + + if (v1 != -1) + { + Pout<< " " << v1 + << " at:" << sourcePoints[v1] + << " weight:" << nearestVertexWeight_[i][1] << nl; + } + if (v2 != -1) + { + Pout<< " " << v2 + << " at:" << sourcePoints[v2] + << " weight:" << nearestVertexWeight_[i][2] << nl; + } + + Pout<< endl; + } } } } @@ -258,13 +317,14 @@ Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation ( const pointField& sourcePoints, const pointField& destPoints, - const scalar perturb + const scalar perturb, + const bool nearestOnly ) : perturb_(perturb), + nearestOnly_(nearestOnly), referenceCS_(calcCoordinateSystem(sourcePoints)), nPoints_(sourcePoints.size()) - { calcWeights(sourcePoints, destPoints); } @@ -279,6 +339,7 @@ Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation ) : perturb_(perturb), + nearestOnly_(false), referenceCS_(referenceCS), nPoints_(sourcePoints.size()) { diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H index e3428d80..4683dcd4 100644 --- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -56,6 +56,9 @@ class pointToPointPlanarInterpolation //- Perturbation factor const scalar perturb_; + //- Whether to use nearest point only (avoids triangulation, projection) + const bool nearestOnly_; + //- Coordinate system coordinateSystem referenceCS_; @@ -91,12 +94,15 @@ public: // Constructors //- Construct from 3D locations. Determines local coordinate system - // from sourcePoints and maps onto that. + // from sourcePoints and maps onto that. If nearestOnly skips any + // local coordinate system and triangulation and uses nearest vertex + // only pointToPointPlanarInterpolation ( const pointField& sourcePoints, const pointField& destPoints, - const scalar perturb + const scalar perturb, + const bool nearestOnly = false ); //- Construct from coordinate system and locations. From c5d310b1b96503e9a603208f2f2aaf41fdc52ac7 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 3 Sep 2014 11:54:30 +0100 Subject: [PATCH 06/15] ENH: mapFields: expose AMI mapping methods --- .../preProcessing/mapFields/mapFields.C | 64 +++- .../AMIInterpolation/AMIInterpolation.C | 191 +++++++---- .../AMIInterpolation/AMIInterpolation.H | 37 +- .../meshToMesh/meshToMesh.C | 322 +++++++++++------- .../meshToMesh/meshToMesh.H | 70 +++- .../meshToMesh/meshToMeshI.H | 7 + 6 files changed, 494 insertions(+), 197 deletions(-) diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C index 6a6e6e0e..5a293d79 100644 --- a/applications/utilities/preProcessing/mapFields/mapFields.C +++ b/applications/utilities/preProcessing/mapFields/mapFields.C @@ -41,7 +41,8 @@ void mapConsistentMesh ( const fvMesh& meshSource, const fvMesh& meshTarget, - const meshToMesh::interpolationMethod& mapMethod, + const word& mapMethod, + const word& AMIMapMethod, const bool subtract, const HashSet& selectedFields, const bool noLagrangian @@ -50,7 +51,7 @@ void mapConsistentMesh Info<< nl << "Consistently creating and mapping fields for time " << meshSource.time().timeName() << nl << endl; - meshToMesh interp(meshSource, meshTarget, mapMethod); + meshToMesh interp(meshSource, meshTarget, mapMethod, AMIMapMethod); if (subtract) { @@ -79,7 +80,8 @@ void mapSubMesh const fvMesh& meshTarget, const HashTable& patchMap, const wordList& cuttingPatches, - const meshToMesh::interpolationMethod& mapMethod, + const word& mapMethod, + const word& AMIMapMethod, const bool subtract, const HashSet& selectedFields, const bool noLagrangian @@ -93,6 +95,7 @@ void mapSubMesh meshSource, meshTarget, mapMethod, + AMIMapMethod, patchMap, cuttingPatches ); @@ -186,6 +189,12 @@ int main(int argc, char *argv[]) "word", "specify the mapping method (direct|mapNearest|cellVolumeWeight)" ); + argList::addOption + ( + "patchMapMethod", + "word", + "specify the patch mapping method (direct|mapNearest|faceAreaWeight)" + ); argList::addBoolOption ( "subtract", @@ -231,17 +240,50 @@ int main(int argc, char *argv[]) const bool consistent = args.optionFound("consistent"); - meshToMesh::interpolationMethod mapMethod = - meshToMesh::imCellVolumeWeight; - if (args.optionFound("mapMethod")) + word mapMethod = meshToMesh::interpolationMethodNames_ + [ + meshToMesh::imCellVolumeWeight + ]; + + if (args.optionReadIfPresent("mapMethod", mapMethod)) { - mapMethod = meshToMesh::interpolationMethodNames_[args["mapMethod"]]; - - Info<< "Mapping method: " - << meshToMesh::interpolationMethodNames_[mapMethod] << endl; + Info<< "Mapping method: " << mapMethod << endl; } + + word patchMapMethod; + if (meshToMesh::interpolationMethodNames_.found(mapMethod)) + { + // Lookup corresponding AMI method + meshToMesh::interpolationMethod method = + meshToMesh::interpolationMethodNames_[mapMethod]; + + patchMapMethod = AMIPatchToPatchInterpolation::interpolationMethodToWord + ( + meshToMesh::interpolationMethodAMI(method) + ); + } + + // Optionally override + if (args.optionFound("patchMapMethod")) + { + patchMapMethod = args["patchMapMethod"]; + + Info<< "Patch mapping method: " << patchMapMethod << endl; + } + + + if (patchMapMethod.empty()) + { + FatalErrorIn(args.executable()) + << "No valid patchMapMethod for method " << mapMethod + << ". Please supply one through the 'patchMapMethod' option" + << exit(FatalError); + } + + + const bool subtract = args.optionFound("subtract"); if (subtract) { @@ -314,6 +356,7 @@ int main(int argc, char *argv[]) meshSource, meshTarget, mapMethod, + patchMapMethod, subtract, selectedFields, noLagrangian @@ -328,6 +371,7 @@ int main(int argc, char *argv[]) patchMap, addProcessorPatches(meshTarget, cuttingPatches), mapMethod, + patchMapMethod, subtract, selectedFields, noLagrangian diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index ea79e069..58fb989e 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -536,66 +536,13 @@ void Foam::AMIInterpolation::agglomerate } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - template -Foam::AMIInterpolation::AMIInterpolation +void Foam::AMIInterpolation::constructFromSurface ( const SourcePatch& srcPatch, const TargetPatch& tgtPatch, - const faceAreaIntersect::triangulationMode& triMode, - const bool requireMatch, - const interpolationMethod& method, - const scalar lowWeightCorrection, - const bool reverseTarget + const autoPtr& surfPtr ) -: - method_(method), - reverseTarget_(reverseTarget), - requireMatch_(requireMatch), - singlePatchProc_(-999), - lowWeightCorrection_(lowWeightCorrection), - srcAddress_(), - srcWeights_(), - srcWeightsSum_(), - tgtAddress_(), - tgtWeights_(), - tgtWeightsSum_(), - triMode_(triMode), - srcMapPtr_(NULL), - tgtMapPtr_(NULL) -{ - update(srcPatch, tgtPatch); -} - - -template -Foam::AMIInterpolation::AMIInterpolation -( - const SourcePatch& srcPatch, - const TargetPatch& tgtPatch, - const autoPtr& surfPtr, - const faceAreaIntersect::triangulationMode& triMode, - const bool requireMatch, - const interpolationMethod& method, - const scalar lowWeightCorrection, - const bool reverseTarget -) -: - method_(method), - reverseTarget_(reverseTarget), - requireMatch_(requireMatch), - singlePatchProc_(-999), - lowWeightCorrection_(lowWeightCorrection), - srcAddress_(), - srcWeights_(), - srcWeightsSum_(), - tgtAddress_(), - tgtWeights_(), - tgtWeightsSum_(), - triMode_(triMode), - srcMapPtr_(NULL), - tgtMapPtr_(NULL) { if (surfPtr.valid()) { @@ -658,6 +605,134 @@ Foam::AMIInterpolation::AMIInterpolation } +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const interpolationMethod& method, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(interpolationMethodToWord(method)), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + update(srcPatch, tgtPatch); +} + + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const word& methodName, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(methodName), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + update(srcPatch, tgtPatch); +} + + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surfPtr, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const interpolationMethod& method, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(interpolationMethodToWord(method)), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + constructFromSurface(srcPatch, tgtPatch, surfPtr); +} + + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surfPtr, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const word& methodName, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(methodName), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + constructFromSurface(srcPatch, tgtPatch, surfPtr); +} + + template Foam::AMIInterpolation::AMIInterpolation ( @@ -666,7 +741,7 @@ Foam::AMIInterpolation::AMIInterpolation const labelList& targetRestrictAddressing ) : - method_(fineAMI.method_), + methodName_(fineAMI.methodName_), reverseTarget_(fineAMI.reverseTarget_), requireMatch_(fineAMI.requireMatch_), singlePatchProc_(fineAMI.singlePatchProc_), @@ -883,7 +958,7 @@ void Foam::AMIInterpolation::update ( AMIMethod::New ( - interpolationMethodToWord(method_), + methodName_, srcPatch, newTgtPatch, srcMagSf_, @@ -1000,7 +1075,7 @@ void Foam::AMIInterpolation::update ( AMIMethod::New ( - interpolationMethodToWord(method_), + methodName_, srcPatch, tgtPatch, srcMagSf_, diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H index 7dbacd9e..27d077f1 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H @@ -110,7 +110,7 @@ private: // Private data //- Interpolation method - interpolationMethod method_; + const word methodName_; //- Flag to indicate that the two patches are co-directional and // that the orientation of the target patch should be reversed @@ -250,7 +250,7 @@ private: ); - // Constructor helper + // Constructor helpers static void agglomerate ( @@ -269,6 +269,12 @@ private: autoPtr& tgtMap ); + void constructFromSurface + ( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surfPtr + ); public: @@ -286,6 +292,19 @@ public: const bool reverseTarget = false ); + //- Construct from components + AMIInterpolation + ( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch = true, + const word& methodName = + interpolationMethodToWord(imFaceAreaWeight), + const scalar lowWeightCorrection = -1, + const bool reverseTarget = false + ); + //- Construct from components, with projection surface AMIInterpolation ( @@ -299,6 +318,20 @@ public: const bool reverseTarget = false ); + //- Construct from components, with projection surface + AMIInterpolation + ( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surf, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch = true, + const word& methodName = + interpolationMethodToWord(imFaceAreaWeight), + const scalar lowWeightCorrection = -1, + const bool reverseTarget = false + ); + //- Construct from agglomeration of AMIInterpolation. Agglomeration // passed in as new coarse size and addressing from fine from coarse AMIInterpolation diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C index 2b2b7e97..91c998ce 100644 --- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C +++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C @@ -119,6 +119,7 @@ void Foam::meshToMesh::normaliseWeights void Foam::meshToMesh::calcAddressing ( + const word& methodName, const polyMesh& src, const polyMesh& tgt ) @@ -127,7 +128,7 @@ void Foam::meshToMesh::calcAddressing ( meshToMeshMethod::New ( - interpolationMethodNames_[method_], + methodName, src, tgt ) @@ -150,11 +151,11 @@ void Foam::meshToMesh::calcAddressing } -void Foam::meshToMesh::calculate() +void Foam::meshToMesh::calculate(const word& methodName) { Info<< "Creating mesh-to-mesh addressing for " << srcRegion_.name() << " and " << tgtRegion_.name() << " regions using " - << interpolationMethodNames_[method_] << endl; + << methodName << endl; singleMeshProc_ = calcDistribution(srcRegion_, tgtRegion_); @@ -241,7 +242,7 @@ void Foam::meshToMesh::calculate() } } - calcAddressing(srcRegion_, newTgt); + calcAddressing(methodName, srcRegion_, newTgt); // per source cell the target cell address in newTgt mesh forAll(srcToTgtCellAddr_, i) @@ -320,7 +321,7 @@ void Foam::meshToMesh::calculate() } else { - calcAddressing(srcRegion_, tgtRegion_); + calcAddressing(methodName, srcRegion_, tgtRegion_); normaliseWeights ( @@ -342,12 +343,9 @@ void Foam::meshToMesh::calculate() Foam::AMIPatchToPatchInterpolation::interpolationMethod -Foam::meshToMesh::interpolationMethodAMI -( - const interpolationMethod method -) const +Foam::meshToMesh::interpolationMethodAMI(const interpolationMethod method) { - switch (method_) + switch (method) { case imDirect: { @@ -374,7 +372,7 @@ Foam::meshToMesh::interpolationMethodAMI "const interpolationMethod method" ") const" ) - << "Unhandled enumeration " << method_ + << "Unhandled enumeration " << method << abort(FatalError); } } @@ -383,90 +381,66 @@ Foam::meshToMesh::interpolationMethodAMI } -const Foam::PtrList& -Foam::meshToMesh::patchAMIs() const +void Foam::meshToMesh::calculatePatchAMIs(const word& AMIMethodName) { - if (patchAMIs_.empty()) + if (!patchAMIs_.empty()) { - const word amiMethod = - AMIPatchToPatchInterpolation::interpolationMethodToWord - ( - interpolationMethodAMI(method_) - ); - - patchAMIs_.setSize(srcPatchID_.size()); - - forAll(srcPatchID_, i) - { - label srcPatchI = srcPatchID_[i]; - label tgtPatchI = tgtPatchID_[i]; - - const polyPatch& srcPP = srcRegion_.boundaryMesh()[srcPatchI]; - const polyPatch& tgtPP = tgtRegion_.boundaryMesh()[tgtPatchI]; - - Info<< "Creating AMI between source patch " << srcPP.name() - << " and target patch " << tgtPP.name() - << " using " << amiMethod - << endl; - - Info<< incrIndent; - - patchAMIs_.set - ( - i, - new AMIPatchToPatchInterpolation - ( - srcPP, - tgtPP, - faceAreaIntersect::tmMesh, - false, - interpolationMethodAMI(method_), - -1, - true // flip target patch since patch normals are aligned - ) - ); - - Info<< decrIndent; - } + FatalErrorIn("meshToMesh::calculatePatchAMIs()") + << "patch AMI already calculated" + << exit(FatalError); } - return patchAMIs_; + patchAMIs_.setSize(srcPatchID_.size()); + + forAll(srcPatchID_, i) + { + label srcPatchI = srcPatchID_[i]; + label tgtPatchI = tgtPatchID_[i]; + + const polyPatch& srcPP = srcRegion_.boundaryMesh()[srcPatchI]; + const polyPatch& tgtPP = tgtRegion_.boundaryMesh()[tgtPatchI]; + + Info<< "Creating AMI between source patch " << srcPP.name() + << " and target patch " << tgtPP.name() + << " using " << AMIMethodName + << endl; + + Info<< incrIndent; + + patchAMIs_.set + ( + i, + new AMIPatchToPatchInterpolation + ( + srcPP, + tgtPP, + faceAreaIntersect::tmMesh, + false, + AMIMethodName, + -1, + true // flip target patch since patch normals are aligned + ) + ); + + Info<< decrIndent; + } } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::meshToMesh::meshToMesh +void Foam::meshToMesh::constructNoCuttingPatches ( - const polyMesh& src, - const polyMesh& tgt, - const interpolationMethod& method, - bool interpAllPatches + const word& methodName, + const word& AMIMethodName, + const bool interpAllPatches ) -: - srcRegion_(src), - tgtRegion_(tgt), - srcPatchID_(), - tgtPatchID_(), - patchAMIs_(), - cuttingPatches_(), - srcToTgtCellAddr_(), - tgtToSrcCellAddr_(), - srcToTgtCellWght_(), - tgtToSrcCellWght_(), - method_(method), - V_(0.0), - singleMeshProc_(-1), - srcMapPtr_(NULL), - tgtMapPtr_(NULL) { if (interpAllPatches) { - const polyBoundaryMesh& srcBM = src.boundaryMesh(); - const polyBoundaryMesh& tgtBM = tgt.boundaryMesh(); + const polyBoundaryMesh& srcBM = srcRegion_.boundaryMesh(); + const polyBoundaryMesh& tgtBM = tgtRegion_.boundaryMesh(); - DynamicList