diff --git a/Examples/ADMFD/elastic.cpp b/Examples/ADMFD/elastic.cpp index 6648d1a6e5f0778a5f077bef1bb690d3e0995e82..bd2dcad6cc2e2f9204fce5492a352ca751518d0e 100644 --- a/Examples/ADMFD/elastic.cpp +++ b/Examples/ADMFD/elastic.cpp @@ -455,6 +455,7 @@ int main(int argc,char ** argv) const rMatrix iBtB(viBtB,6,6); const rMatrix Curl(vCurl,9,9); const rMatrix I = rMatrix::Unit(3); + const rMatrix I9 = rMatrix::Unit(9); //PrintSV(B); std::cout << "B^T*B" << std::endl; @@ -492,10 +493,10 @@ int main(int argc,char ** argv) ttt = Timer(); //Assemble gradient matrix W on cells #if defined(USE_OMP) -#pragma omp parallel +//#pragma omp parallel #endif { - rMatrix N, R, L, T, K(9,9), C(6,6), W1, W2, W3, U,S,V, w, u, v; + rMatrix N, R, L, M, T, K(9,9), C(6,6), W1, W2, W3, U,S,V, w, u, v; rMatrix x(3,1), xf(3,1), n(3,1); double area; //area of the face double volume; //volume of the cell @@ -518,9 +519,11 @@ int main(int argc,char ** argv) //K += rMatrix::Unit(9)*1.0e-6*K.FrobeniusNorm(); //PrintSV(K); N.Resize(3*NF,9); //co-normals + //NQ.Resize(3*NF,9); //T.Resize(3*NF,9); //transversals R.Resize(3*NF,9); //directions L.Resize(3*NF,3*NF); + M.Resize(3*NF,3*NF); L.Zero(); //A.Resize(3*NF,3*NF); @@ -540,6 +543,32 @@ int main(int argc,char ** argv) N(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose()); + /* + Cell c2 = cell.Neighbour(faces[k]); + + if( c2.isValid() ) + { + KTensor(tag_C[c2],K2); + CTensor(tag_C[c2],C2); + //T1 = I.Kronecker(n.Transpose())*K*I.Kronecker(n); + T2 = I.Kronecker(n.Transpose())*K2*I.Kronecker(n); + //TagRealArray tag_G = m->GetTag("REFERENCE_GRADIENT"); + Q1 = I9 + I.Kronecker(n)*T2.Invert()*I.Kronecker(n.Transpose())*(K-K2); + //Q2 = I9 + I.Kronecker(n)*T1.Invert()*I.Kronecker(n.Transpose())*(K2-K); + //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*(I9 - B*Q*((B*Q).Transpose()*B*Q)*(B*Q).Transpose()); + Q1 = Q1.Root(); + //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*(I9 - (Q1.Transpose()*B)*((Q1.Transpose()*B).Transpose()*(Q1.Transpose()*B))*(Q1.Transpose()*B).Transpose()); + //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*Q1*(I9 - B*(B.Transpose()*B)*B.Transpose()); + //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*(I9 - B*(B.Transpose()*B)*B.Transpose())*Q1; + NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*Q1; + //NQ(3*k,3*(k+1),0,9).Zero(); + } + else + { + NQ(3*k,3*(k+1),0,9).Zero(); + } + */ + //std::cout << "I\otimes n^T" << std::endl; //I.Kronecker(n.Transpose()).Print(); @@ -642,30 +671,61 @@ int main(int argc,char ** argv) //W2 = L - (L*R*B)*((L*R*B).Transpose()*R*B).CholeskyInvert()*(L*R*B).Transpose(); //W2 = W1.Trace()*(rMatrix::Unit(3*NF) - (R*B)*((R*B).Transpose()*R*B).Invert()*(R*B).Transpose()); - if( true ) + if( false ) { double alpha = 1; double beta = alpha; //R = R*(rMatrix::Unit(9) + B*iBtB*B.Transpose())*0.5; - K += 1.0e-5*K.Trace()*(rMatrix::Unit(9) - B*iBtB*B.Transpose()); + //if( cell.GetElementType() == Element::Tet ) + //K += 1.0e-3*K.Trace()*(rMatrix::Unit(9) - B*iBtB*B.Transpose()); + + + W1 = (N*K+alpha*L*R)*((N*K+alpha*L*R).Transpose()*R).PseudoInvert(1.0e-11)*(N*K+alpha*L*R).Transpose(); + + //R += N*(rMatrix::Unit(9) - B*iBtB*B.Transpose())*K.Trace()*2*NF/volume; + + W2 = L - (1+beta)*(L*R)*((L*R).Transpose()*R).PseudoInvert(1.0e-11)*(L*R).Transpose(); - W1 = (N*K+alpha*L*R)*((N*K+alpha*L*R).Transpose()*R).Invert()*(N*K+alpha*L*R).Transpose(); - W2 = L - (1+beta)*(L*R)*((L*R).Transpose()*R).PseudoInvert()*(L*R).Transpose(); + if( cell.GetElementType() == Element::Tet ) + { + //R = R*B*iBtB; + //W2 += (L - (L*R)*((L*R).Transpose()*L*R).PseudoInvert(1.0e-11)*(L*R).Transpose()); + N = N*B; + R = R*B*iBtB;//*B.Transpose(); + + //W1 = (N*C)*((N*C).Transpose()*R).Invert()*(N*C).Transpose(); + W2 += 1.0e-5*(L - (L*R)*((L*R).Transpose()*R).PseudoInvert(1.0e-11)*(L*R).Transpose()); + } + +#pragma omp critical + { + std::cout << "W1: "; PrintSV(W1); + std::cout << "W2: "; PrintSV(W2); + std::cout << "S : "; PrintSV(W1+W2); + } } else { - //W2 = L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose(); + N = N*B; - R = R*B*iBtB; + R = R*B;//*iBtB; + + W1 = (N*C*(B.Transpose()*B))*((N*C*(B.Transpose()*B)).Transpose()*R).Invert()*(N*C*(B.Transpose()*B)).Transpose(); + W2 = L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose(); - W1 = (N*C)*((N*C).Transpose()*R).Invert()*(N*C).Transpose(); //double alpha = 1; //W1 = (N*K+alpha*L*R)*((N*K+alpha*L*R).Transpose()*R).Invert()*(N*K+alpha*L*R).Transpose() - alpha*(L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose(); +#pragma omp critical + { + std::cout << "W1: "; PrintSV(W1); + std::cout << "W2: "; PrintSV(W2); + std::cout << "S : "; PrintSV(W1+W2); + } //R = R*B*iBtB; - W2 += L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose(); + //W2 = L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose(); //W2 = 2*W1.Trace()/(3*NF) *(rMatrix::Unit(3*NF) - (R)*((R).Transpose()*R).Invert()*(R).Transpose()); /* #pragma omp critical @@ -1192,8 +1252,8 @@ int main(int argc,char ** argv) //R.GetJacobian().Save("A.mtx",&Text); - Solver S(Solver::INNER_MPTILU2); - //Solver S(Solver::INNER_MPTILUC); + //Solver S(Solver::INNER_MPTILU2); + Solver S(Solver::INNER_MPTILUC); //Solver S("superlu"); S.SetParameter("relative_tolerance", "1.0e-14"); S.SetParameter("absolute_tolerance", "1.0e-12"); diff --git a/Source/Headers/inmost_dense.h b/Source/Headers/inmost_dense.h index 5c32f79e0a9286ab74e305a76475a99910cb98b1..88c9b6d268ec2cabcfdcadf22fcff5839f289a0d 100644 --- a/Source/Headers/inmost_dense.h +++ b/Source/Headers/inmost_dense.h @@ -751,9 +751,25 @@ namespace INMOST /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. /// If *ierr == -1 on input, then prints out information in case of failure. /// In case of failure *ierr = 1, in case of no failure *ierr = 0. - /// @return A pair of pseudo-inverse matrix and boolean. If boolean is true, - /// then the matrix was inverted successfully. + /// @return A pseudo-inverse of the matrix. Matrix > PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0, int * ierr = NULL) const; + + /// Calcuate A^n, where n is some real value. + /// @param n Real value. + /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. + /// If *ierr == -1 on input, then prints out information in case of failure. + /// In case of failure *ierr = 1, in case of no failure *ierr = 0. + /// The error may be caused by error in SVD algorithm. + /// @return Matrix in power of n. + //Matrix > Power(INMOST_DATA_REAL_TYPE n = 1, int * ierr = NULL) const; + /// Calculate square root of A matrix by Babylonian method. + /// @param iter Number of iterations. + /// @param tol Convergence tolerance. + /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. + /// If *ierr == -1 on input, then prints out information in case of failure. + /// In case of failure *ierr = 1, in case of no failure *ierr = 0. + /// @return Square root of a matrix. + Matrix > Root(INMOST_DATA_ENUM_TYPE iter = 25, INMOST_DATA_REAL_TYPE tol = 1.0e-7, int * ierr = NULL) const; /// Solves the system of equations of the form A*X=B, with A and B matrices. /// Uses Moore-Penrose pseudo-inverse of the matrix A and calculates X = A^+*B. /// @param B Matrix at the right hand side. @@ -2571,7 +2587,58 @@ namespace INMOST if( ierr ) *ierr = 0; return ret; } - + template + Matrix > + AbstractMatrix::Root(INMOST_DATA_ENUM_TYPE iter, INMOST_DATA_REAL_TYPE tol, int *ierr) const + { + assert(Rows() == Cols()); + Matrix > ret(Cols(),Cols()); + Matrix > ret0(Cols(),Cols()); + ret.Zero(); + ret0.Zero(); + int k = 0; + for(k = 0; k < Cols(); ++k) ret(k,k) = ret0(k,k) = 1; + while(k < iter) + { + ret0 = ret; + ret = 0.5*(ret + (*this)*ret.Invert()); + if( (ret - ret0).FrobeniusNorm() < tol ) return ret; + } + if( ierr ) + { + if( *ierr == -1 ) std::cout << "Failed to find square root of matrix by Babylonian method" << std::endl; + *ierr = 1; + return ret; + } + return ret; + } + /* + template + Matrix > + AbstractMatrix::Power(INMOST_DATA_REAL_TYPE n, int * ierr) const + { + Matrix > ret(Cols(),Rows()); + Matrix > L,S,iL; + bool success = Eigensolver(L,S,iL); + if( !success ) + { + if( ierr ) + { + if( *ierr == -1 ) std::cout << "Failed to compute eigenvalue decomposition of the matrix" << std::endl; + *ierr = 1; + return ret; + } + else throw MatrixEigensolverFail; + } + for(INMOST_DATA_ENUM_TYPE k = 0; k < S.Cols(); ++k) S(k,k) = pow(S(k,k),n); + if( n >= 0 ) + ret = U*S*V.Transpose(); + else + ret = V*S*U.Transpose(); + if( ierr ) *ierr = 0; + return ret; + } + */ template template Matrix::type, pool_array_t::type> >