From bdb4477e567fe372e6eb4d311953af1595b0512e Mon Sep 17 00:00:00 2001 From: Kirill Terekhov Date: Tue, 13 Nov 2018 16:38:05 +0300 Subject: [PATCH] Add elasticity example and fix a bug Elasticity example in Examples/ADMFD/elasticity.cpp Bug fix - when assigning a matrix with different storage container, the number of rows and columns is changed after the assignment itself, which leads to incorrect behavior of opeartor(i,j) of the matrix. --- Examples/ADMFD/CMakeLists.txt | 100 ++- Examples/ADMFD/{main.cpp => diffusion.cpp} | 0 Examples/ADMFD/elastic.cpp | 940 +++++++++++++++++++++ Source/Headers/inmost_dense.h | 65 +- 4 files changed, 1052 insertions(+), 53 deletions(-) rename Examples/ADMFD/{main.cpp => diffusion.cpp} (100%) create mode 100644 Examples/ADMFD/elastic.cpp diff --git a/Examples/ADMFD/CMakeLists.txt b/Examples/ADMFD/CMakeLists.txt index ceb3331..4346734 100644 --- a/Examples/ADMFD/CMakeLists.txt +++ b/Examples/ADMFD/CMakeLists.txt @@ -1,58 +1,110 @@ project(ADMFD) -set(SOURCE main.cpp) -add_executable(ADMFD ${SOURCE}) +add_executable(MFDDIFF diffusion.cpp) +add_executable(MFDELAST elastic.cpp) -target_link_libraries(ADMFD inmost) +target_link_libraries(MFDDIFF inmost) +target_link_libraries(MFDELAST inmost) if(USE_SOLVER) if(USE_SOLVER_ANI) - message("linking ADMFD with ani3d and BLAS") - target_link_libraries(ADMFD ani3d ${BLAS_LIBRARIES}) + message("linking MFDDIFF with ani3d and BLAS") + target_link_libraries(MFDDIFF ani3d ${BLAS_LIBRARIES}) if(BLAS_LINKER_FLAGS) - set_target_properties(ADMFD PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}") + set_target_properties(MFDDIFF PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}") endif() endif() if(USE_SOLVER_PETSC) - message("linking ADMFD with PETSc") - target_link_libraries(ADMFD ${PETSC_LIBRARIES}) + message("linking MFDDIFF with PETSc") + target_link_libraries(MFDDIFF ${PETSC_LIBRARIES}) endif() if(USE_SOLVER_TRILINOS) - message("linking ADMFD with Trilinos") - target_link_libraries(ADMFD ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES}) + message("linking MFDDIFF with Trilinos") + target_link_libraries(MFDDIFF ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES}) endif() if(USE_SOLVER_METIS) - message("linking ADMFD with Metis") - target_link_libraries(ADMFD ${METIS_LIBRARIES}) + message("linking MFDDIFF with Metis") + target_link_libraries(MFDDIFF ${METIS_LIBRARIES}) endif() if(USE_SOLVER_MONDRIAAN) - message("linking ADMFD with Mondriaan") - target_link_libraries(ADMFD ${MONDRIAAN_LIBRARIES}) + message("linking MFDDIFF with Mondriaan") + target_link_libraries(MFDDIFF ${MONDRIAAN_LIBRARIES}) endif() if(USE_SOLVER_SUPERLU) - message("linking ADMFD with SuperLU") - target_link_libraries(ADMFD ${SUPERLU_LIBRARIES}) + message("linking MFDDIFF with SuperLU") + target_link_libraries(MFDDIFF ${SUPERLU_LIBRARIES}) endif() endif() if(USE_PARTITIONER) if(USE_PARTITIONER_ZOLTAN) - message("linking ADMFD with Zoltan") - target_link_libraries(ADMFD ${ZOLTAN_LIBRARIES}) + message("linking MFDDIFF with Zoltan") + target_link_libraries(MFDDIFF ${ZOLTAN_LIBRARIES}) endif() if(USE_PARTITIONER_PARMETIS) - message("linking ADMFD with ParMETIS") - target_link_libraries(ADMFD ${PARMETIS_LIBRARIES}) + message("linking MFDDIFF with ParMETIS") + target_link_libraries(MFDDIFF ${PARMETIS_LIBRARIES}) endif() endif() if(USE_MPI) - message("linking ADMFD with MPI") - target_link_libraries(ADMFD ${MPI_LIBRARIES}) + message("linking MFDDIFF with MPI") + target_link_libraries(MFDDIFF ${MPI_LIBRARIES}) if(MPI_LINK_FLAGS) - set_target_properties(ADMFD PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + set_target_properties(MFDDIFF PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") endif() endif(USE_MPI) -install(TARGETS ADMFD EXPORT inmost-targets RUNTIME DESTINATION bin) +if(USE_SOLVER) + if(USE_SOLVER_ANI) + message("linking MFDELAST with ani3d and BLAS") + target_link_libraries(MFDELAST ani3d ${BLAS_LIBRARIES}) + if(BLAS_LINKER_FLAGS) + set_target_properties(MFDELAST PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}") + endif() + endif() + if(USE_SOLVER_PETSC) + message("linking MFDELAST with PETSc") + target_link_libraries(MFDELAST ${PETSC_LIBRARIES}) + endif() + if(USE_SOLVER_TRILINOS) + message("linking MFDELAST with Trilinos") + target_link_libraries(MFDELAST ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES}) + endif() + if(USE_SOLVER_METIS) + message("linking MFDELAST with Metis") + target_link_libraries(MFDELAST ${METIS_LIBRARIES}) + endif() + if(USE_SOLVER_MONDRIAAN) + message("linking MFDELAST with Mondriaan") + target_link_libraries(MFDELAST ${MONDRIAAN_LIBRARIES}) + endif() + if(USE_SOLVER_SUPERLU) + message("linking MFDELAST with SuperLU") + target_link_libraries(MFDELAST ${SUPERLU_LIBRARIES}) + endif() +endif() + + +if(USE_PARTITIONER) + if(USE_PARTITIONER_ZOLTAN) + message("linking MFDELAST with Zoltan") + target_link_libraries(MFDELAST ${ZOLTAN_LIBRARIES}) + endif() + if(USE_PARTITIONER_PARMETIS) + message("linking MFDELAST with ParMETIS") + target_link_libraries(MFDELAST ${PARMETIS_LIBRARIES}) + endif() +endif() + +if(USE_MPI) + message("linking MFDELAST with MPI") + target_link_libraries(MFDELAST ${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(MFDELAST PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) + +install(TARGETS MFDDIFF EXPORT inmost-targets RUNTIME DESTINATION bin) +install(TARGETS MFDELAST EXPORT inmost-targets RUNTIME DESTINATION bin) diff --git a/Examples/ADMFD/main.cpp b/Examples/ADMFD/diffusion.cpp similarity index 100% rename from Examples/ADMFD/main.cpp rename to Examples/ADMFD/diffusion.cpp diff --git a/Examples/ADMFD/elastic.cpp b/Examples/ADMFD/elastic.cpp new file mode 100644 index 0000000..4f5e5f2 --- /dev/null +++ b/Examples/ADMFD/elastic.cpp @@ -0,0 +1,940 @@ +#include "inmost.h" +#include +#include +#include +#include + + +using namespace INMOST; + +#ifndef M_PI +#define M_PI 3.141592653589 +#endif + +#if defined(USE_MPI) +#define BARRIER MPI_Barrier(MPI_COMM_WORLD); +#else +#define BARRIER +#endif + +//shortcuts +typedef Storage::bulk bulk; +typedef Storage::real real; +typedef Storage::integer integer; +typedef Storage::enumerator enumerator; +typedef Storage::real_array real_array; +typedef Storage::var_array var_array; + +bool print_niter = false; //save file on nonlinear iterations + +//If true, then tag ELASTIC_TENSOR contains stiffness tensor C, +//determining strain to stress relation C*\varepsilon = \sigma, otherwise +//the tag contains compliance tensor S, determining stress to strain relation +//S*\sigma = \varepsilon. Here \varepsilon = \frac{G+G^T}{2}, +//with G - gradient of displacement matrix +bool inverse_tensor = true; + +//indices for tensor +const int IC11 = 0; +const int IC12 = 1; +const int IC13 = 2; +const int IC14 = 3; +const int IC15 = 4; +const int IC16 = 5; +const int IC22 = 6; +const int IC23 = 7; +const int IC24 = 8; +const int IC25 = 9; +const int IC26 = 10; +const int IC33 = 11; +const int IC34 = 12; +const int IC35 = 13; +const int IC36 = 14; +const int IC44 = 15; +const int IC45 = 16; +const int IC46 = 17; +const int IC55 = 18; +const int IC56 = 19; +const int IC66 = 20; + + +/* + + stored tensor + C_1 C_2 C_3 C_4 C_5 C_6 + C_7 C_8 C_9 C_10 C_11 + C_12 C_13 C_14 C_15 + C_16 C_17 C_18 + C_19 C_20 + C_21 + + + voigt notation tensor + + c_11 c_12 c_13 c_14 c_15 c_16 + c_12 c_22 c_23 c_24 c_25 c_26 + c_13 c_23 c_33 c_34 c_35 c_36 + c_14 c_24 c_34 c_44 c_45 c_46 + c_15 c_25 c_35 c_45 c_55 c_56 + c_16 c_26 c_36 c_46 c_56 c_66 + + co-normal tensors + + | K_1 \vec{n} \cdot \nabla u + K_4 \vec{n} \cdot \nabla v + K_5 \vec{n} \cdot \nabla w | + F = | K_4^T \vec{n} \cdot \nabla u + K_2 \vec{n} \cdot \nabla v + K_6 \vec{n} \cdot \nabla w | + | K_5^T \vec{n} \cdot \nabla u + K_6^T \vec{n} \cdot \nabla v + K_3 \vec{n} \cdot \nabla w | + + tensors representations + + | c_11 c_16 c_15 | + K_1 = | c_16 c_66 c_56 | + | c_15 c_56 c_55 | + + | c_66 c_26 c_46 | + K_2 = | c_26 c_22 c_24 | + | c_46 c_24 c_44 | + + | c_55 c_45 c_35 | + K_3 = | c_45 c_44 c_34 | + | c_35 c_34 c_33 | + + | c_16 c_66 c_56 | + K_4 = | c_12 c_26 c_25 | + | c_14 c_46 c_45 | + + | c_15 c_56 c_55 | + K_5 = | c_14 c_46 c_45 | + | c_13 c_36 c_35 | + + | c_56 c_25 c_45 | + K_6 = | c_46 c_24 c_44 | + | c_36 c_23 c_34 | + + two-point flux approximation + + F = T (\vec{u}_2 - \vec{u}_1) + + here the transimissibility matrix is + + T = C'_1 C'_2 (d_1 C'_2 + d_2 C'_1)^{-1} + + | k_1 k_4 k_5 | + C' = | k_4 k_2 k_6 | + | k_5 k_6 k_3 | + + with + + k_i = \vec{n} \cdot K_i \vec{n} = \vec{n} \cdot K_i^T \vec{n} + + */ + + +void CTensor(const real_array & Cv, rMatrix & C) +{ + C = rMatrix::FromTensor(Cv.data(),Cv.size(),6); +} + +void KTensor(const real_array & Cv, rMatrix & K) +{ + int si = 0, sj = 0; + rMatrix C = rMatrix::FromTensor(Cv.data(),Cv.size(),6); + + K(0+si,0+sj) = C(0,0); K(0+si,1+sj) = C(0,5); K(0+si,2+sj) = C(0,4); + K(1+si,0+sj) = C(0,5); K(1+si,1+sj) = C(5,5); K(1+si,2+sj) = C(4,5); + K(2+si,0+sj) = C(0,4); K(2+si,1+sj) = C(4,5); K(2+si,2+sj) = C(4,4); + + si = sj = 3; + + K(0+si,0+sj) = C(5,5); K(0+si,1+sj) = C(1,5); K(0+si,2+sj) = C(3,5); + K(1+si,0+sj) = C(1,5); K(1+si,1+sj) = C(1,1); K(1+si,2+sj) = C(1,3); + K(2+si,0+sj) = C(3,5); K(2+si,1+sj) = C(1,3); K(2+si,2+sj) = C(3,3); + + si = sj = 6; + + K(0+si,0+sj) = C(4,4); K(0+si,1+sj) = C(3,4); K(0+si,2+sj) = C(2,4); + K(1+si,0+sj) = C(3,4); K(1+si,1+sj) = C(3,3); K(1+si,2+sj) = C(2,3); + K(2+si,0+sj) = C(2,4); K(2+si,1+sj) = C(2,3); K(2+si,2+sj) = C(2,2); + + si = 0, sj = 3; + + K(0+si,0+sj) = C(0,5); K(0+si,1+sj) = C(0,1); K(0+si,2+sj) = C(0,3); + K(1+si,0+sj) = C(5,5); K(1+si,1+sj) = C(1,5); K(1+si,2+sj) = C(3,5); + K(2+si,0+sj) = C(4,5); K(2+si,1+sj) = C(1,4); K(2+si,2+sj) = C(3,4); + + si = 3, sj = 0; + + K(0+si,0+sj) = C(0,5); K(0+si,1+sj) = C(5,5); K(0+si,2+sj) = C(4,5); + K(1+si,0+sj) = C(0,1); K(1+si,1+sj) = C(1,5); K(1+si,2+sj) = C(1,4); + K(2+si,0+sj) = C(0,3); K(2+si,1+sj) = C(3,5); K(2+si,2+sj) = C(3,4); + + si = 0, sj = 6; + + K(0+si,0+sj) = C(0,4); K(0+si,1+sj) = C(0,3); K(0+si,2+sj) = C(0,2); + K(1+si,0+sj) = C(4,5); K(1+si,1+sj) = C(3,5); K(1+si,2+sj) = C(2,5); + K(2+si,0+sj) = C(4,4); K(2+si,1+sj) = C(3,4); K(2+si,2+sj) = C(2,4); + + si = 6, sj = 0; + + K(0+si,0+sj) = C(0,4); K(0+si,1+sj) = C(4,5); K(0+si,2+sj) = C(4,4); + K(1+si,0+sj) = C(0,3); K(1+si,1+sj) = C(3,5); K(1+si,2+sj) = C(3,4); + K(2+si,0+sj) = C(0,2); K(2+si,1+sj) = C(2,5); K(2+si,2+sj) = C(2,4); + + si = 3, sj = 6; + + K(0+si,0+sj) = C(4,5); K(0+si,1+sj) = C(3,5); K(0+si,2+sj) = C(2,5); + K(1+si,0+sj) = C(1,4); K(1+si,1+sj) = C(1,3); K(1+si,2+sj) = C(1,2); + K(2+si,0+sj) = C(3,4); K(2+si,1+sj) = C(3,3); K(2+si,2+sj) = C(2,3); + + si = 6, sj = 3; + + K(0+si,0+sj) = C(4,5); K(0+si,1+sj) = C(1,4); K(0+si,2+sj) = C(3,4); + K(1+si,0+sj) = C(3,5); K(1+si,1+sj) = C(1,3); K(1+si,2+sj) = C(3,3); + K(2+si,0+sj) = C(2,5); K(2+si,1+sj) = C(1,2); K(2+si,2+sj) = C(2,3); +} + +void GetBC(const real_array & bc, const rMatrix & n, rMatrix & Ra, rMatrix & Rb, rMatrix & r) +{ + const rMatrix I = rMatrix::Unit(3); + double alpha_perp, alpha_parallel, beta_perp, beta_parallel; + + if( bc.size() == 6 ) + { + double alpha, beta, proj; + alpha = bc[0]; + beta = bc[1]; + proj = bc[2]; + r(0,0) = bc[3]; + r(1,0) = bc[4]; + r(2,0) = bc[5]; + if( proj ) + { + alpha_perp = alpha; + beta_parallel = beta; + alpha_parallel = beta_perp = 0; + } + else + { + alpha_perp = alpha_parallel = alpha; + beta_perp = beta_parallel = beta; + } + } + else if( bc.size() == 7 ) + { + alpha_perp = bc[0]; + beta_perp = bc[1]; + alpha_parallel = bc[2]; + beta_parallel = bc[3]; + r(0,0) = bc[4]; + r(1,0) = bc[5]; + r(2,0) = bc[6]; + } + Ra = alpha_parallel*I + (alpha_perp-alpha_parallel)*n*n.Transpose(); + Rb = beta_parallel*I + (beta_perp-beta_parallel)*n*n.Transpose(); +} + +void PrintSV(const rMatrix & A) +{ + rMatrix U,S,V; + A.SVD(U,S,V); + std::cout << "singular values:"; + int cnt = 0; + for(int k = 0; k < A.Rows(); ++k) if( fabs(S(k,k)) > 1.0e-10 ) + { + std::cout << " " << S(k,k); + cnt++; + } + else + std::cout << " 0"; + std::cout << " count: " << cnt << "/" << A.Rows() << std::endl; +} + +void PrintRS(const rMatrix & A) +{ + double sum; + std::cout << "row sum:"; + for(int i = 0; i < A.Rows(); ++i) + { + sum = 0; + for(int j = 0; j < A.Cols(); ++j) + sum+= A(i,j); + std::cout << " " << sum; + } + std::cout << std::endl; +} + +int main(int argc,char ** argv) +{ + Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity +#if defined(USE_PARTITIONER) + Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity +#endif + if( argc > 1 ) + { + + double ttt; // Variable used to measure timing + bool repartition = false; // Is it required to redistribute the mesh? + Mesh * m = new Mesh(); // Create an empty mesh + { // Load the mesh + ttt = Timer(); + m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh + if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl; + if( m->isParallelFileFormat(argv[1]) ) //The format is + { + m->Load(argv[1]); // Load mesh from the parallel file format + repartition = true; // Ask to repartition the mesh + } + else if( m->GetProcessorRank() == 0 ) m->Load(argv[1]); // Load mesh from the serial file format + BARRIER + if( m->GetProcessorRank() == 0 ) std::cout << "Load the mesh: " << Timer()-ttt << std::endl; + } + + +#if defined(USE_PARTITIONER) + if (m->GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner + { + { // Compute mesh partitioning + ttt = Timer(); + Partitioner p(m); //Create Partitioning object + p.SetMethod(Partitioner::INNER_KMEANS,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner + p.Evaluate(); // Compute the partitioner and store new processor ID in the mesh + BARRIER + if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl; + } + + { //Distribute the mesh + ttt = Timer(); + m->Redistribute(); // Redistribute the mesh data + m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring + BARRIER + if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl; + } + } +#endif + + { // prepare geometrical data on the mesh + ttt = Timer(); + Mesh::GeomParam table; + table[CENTROID] = CELL | FACE; //Compute averaged center of mass + table[NORMAL] = FACE; //Compute normals + table[ORIENTATION] = FACE; //Check and fix normal orientation + table[MEASURE] = CELL | FACE; //Compute volumes and areas + //table[BARYCENTER] = CELL | FACE; //Compute volumetric center of mass + m->PrepareGeometricData(table); //Ask to precompute the data + BARRIER + if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl; + } + + // data tags for + TagRealArray tag_UVW;// Displacement + TagRealArray tag_C; // Elasticity tensor in Voigt format + TagRealArray tag_F; // Forcing term + TagRealArray tag_BC; // Boundary conditions + TagRealArray tag_W; // Gradient matrix + TagRealArray tag_FLUX; // Flux (for error) + + if( m->GetProcessorsNumber() > 1 ) //skip for one processor job + { // Exchange ghost cells + ttt = Timer(); + m->ExchangeGhost(1,FACE); // Produce layer of ghost cells + BARRIER + if( m->GetProcessorRank() == 0 ) std::cout << "Exchange ghost: " << Timer()-ttt << std::endl; + } + + std::cout << "Cells " << m->TotalNumberOf(CELL) << std::endl; + std::cout << "Faces " << m->TotalNumberOf(FACE) << std::endl; + std::cout << "Edges " << m->TotalNumberOf(EDGE) << std::endl; + std::cout << "Nodes " << m->TotalNumberOf(NODE) << std::endl; + + const double vB[] = + { + 1,0,0,0,0,0, + 0,0,0,0,0,1, + 0,0,0,0,1,0, + 0,0,0,0,0,1, + 0,1,0,0,0,0, + 0,0,0,1,0,0, + 0,0,0,0,1,0, + 0,0,0,1,0,0, + 0,0,1,0,0,0 + }; + const rMatrix B(vB,9,6); + const rMatrix I = rMatrix::Unit(3); + + (B.Transpose()*B).Print(); + + { //initialize data + if( m->HaveTag("ELASTIC_TENSOR") ) // is elasticity tensor already defined on the mesh? + tag_C = m->GetTag("ELASTIC_TENSOR"); // get the elasticity tensor + else + { + std::cout << "No ELASTIC_TENSOR on mesh" << std::endl; + return -1; + } + + if( m->HaveTag("FORCE") ) //Is there force on the mesh? + { + tag_F = m->GetTag("FORCE"); //initial force + assert(tag_F.isDefined(CELL)); //assuming it was defined on cells + } // end of force + + if( m->HaveTag("BOUNDARY_CONDITION") ) //Is there boundary condition on the mesh? + { + tag_BC = m->GetTag("BOUNDARY_CONDITION"); + //initialize unknowns at boundary + } + tag_W = m->CreateTag("W",DATA_REAL,CELL,NONE); + + ttt = Timer(); + //Assemble gradient matrix W on cells +#if defined(USE_OMP) +#pragma omp parallel +#endif + { + rMatrix N, R, L, T, K(9,9), C(6,6), U,S,V, W1, W2; + rMatrix x(3,1), xf(3,1), n(3,1); + double area; //area of the face + double volume; //volume of the cell + double dist; //distance from cell center to face +#if defined(USE_OMP) +#pragma omp for +#endif + for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) + { + Mesh * mesh = m; + Cell cell = m->CellByLocalID(q); + ElementArray faces = cell->getFaces(); //obtain faces of the cell + int NF = (int)faces.size(); //number of faces; + volume = cell->Volume(); //volume of the cell + cell->Centroid(x.data()); + //get permeability for the cell + KTensor(tag_C[cell],K); + //CTensor(tag_C[cell],C); + if( !inverse_tensor ) K = K.Invert(); + //K += rMatrix::Unit(9)*1.0e-6*K.FrobeniusNorm(); + //PrintSV(K); + N.Resize(3*NF,9); //co-normals + //T.Resize(3*NF,9); //transversals + R.Resize(3*NF,9); //directions + L.Resize(3*NF,3*NF); + L.Zero(); + //A.Resize(3*NF,3*NF); + //A.Zero(); + for(int k = 0; k < NF; ++k) //loop over faces + { + area = faces[k].Area(); + faces[k].Centroid(xf.data()); + faces[k].OrientedUnitNormal(cell->self(),n.data()); + dist = n.DotProduct(xf-x); + // assemble matrix of directions + R(3*k,3*(k+1),0,9) = I.Kronecker((xf-x).Transpose()); + // assemble matrix of co-normals + + L(3*k,3*(k+1),3*k,3*(k+1)) = (area/dist)*I.Kronecker(n.Transpose())*K*I.Kronecker(n); + + N(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose()); + + //std::cout << "I\otimes n^T" << std::endl; + //I.Kronecker(n.Transpose()).Print(); + + //std::cout << "I\otimes n^T B" << std::endl; + + //(I.Kronecker(n.Transpose())*B).Print(); + + //T(3*k,3*(k+1),0,9) = N(3*k,3*(k+1),0,9)*K - L(3*k,3*(k+1),3*k,3*(k+1))*I.Kronecker(n.Transpose()); + //A(3*k,3*(k+1),3*k,3*(k+1)) = I*area; + //NK(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose()); + } //end of loop over faces + tag_W[cell].resize(9*NF*NF); + //int ierr = -1; + //tag_W(cell,3*NF,3*NF) = N*(N.Transpose()*R).PseudoInvert(1.0e-9)*N.Transpose() //stability part + // + (rMatrix::Unit(3*NF) - R*(R.Transpose()*R).PseudoInvert(1.0e-9)*R.Transpose())* + //(2.0/(static_cast(NF)*volume)*(N*K.PseudoInvert(1.0e-9)*N.Transpose()).Trace()); + + //+ (rMatrix::Unit(3*NF) - R*(N.Transpose()*R).PseudoInvert(1.0e-12)*N.Transpose())* + // (2.0/(static_cast(NF)*volume)*(N*K.Invert()*N.Transpose()).Trace()); + + //tag_W(cell,3*NF,3*NF) = N*(N.Transpose()*R).PseudoInvert(1.0e-9)*N.Transpose() + (L - L*R*((L*R).Transpose()*R).PseudoInvert(1.0e-9)*(L*R).Transpose()); + + + //tag_W(cell,3*NF,3*NF) = L + (N*K-L*R)*((N*K-L*R).Transpose()*R).PseudoInvert(1.0e-9)*(N*K-L*R).Transpose() + (L - L*R*((L*R).Transpose()*R).PseudoInvert(1.0e-9)*(L*R).Transpose()); + + //tag_W(cell,3*NF,3*NF) = (N*K)*((N*K).Transpose()*R).PseudoInvert(1.0e-9)*(N*K).Transpose() + (L - L*R*((L*R).Transpose()*R).PseudoInvert(1.0e-9)*(L*R).Transpose()); + + //W1.Resize(3*NF,3*NF); + //W2.Resize(3*NF,3*NF); + + W1 = (N*K)*((N*K).Transpose()*R).PseudoInvert(1.0e-7)*(N*K).Transpose(); + W2 = L - (L*R)*((L*R).Transpose()*R).PseudoInvert(1.0e-7)*(L*R).Transpose(); + + tag_W(cell,3*NF,3*NF) = W1+W2; + //W2 = W1.Trace()*(rMatrix::Unit(3*NF) - R*(R.Transpose()*R).PseudoInvert(1.0e-7)*R.Transpose()); + //W2 = W1.Trace()*(rMatrix::Unit(3*NF) - R*((N*K).Transpose()*R).PseudoInvert(1.0e-7)*(N*K).Transpose()); + /* + std::cout << "L " << std::endl; + L.Print(); + + std::cout << "N*K*N^T" << std::endl; + (N*K*N.Transpose()).Print(); + + std::cout << "N*R^T" << std::endl; + (N*R.Transpose()).Print(); + */ + //L = N*K*N.Transpose()*(N*R.Transpose()).PseudoInvert(1.0e-9); + //std::cout << "L2 " << std::endl; + // L.Print(); + + //PrintSV(W1); + //PrintSV(W2); + //PrintSV(W1+W2); + + //tag_W(cell,3*NF,3*NF) = W1+W2; + + //std::cout << "W1(" << W1.Rows() << "," << W1.Cols() << ")" << std::endl; + //W1.Print(); + + //std::cout << "W2(" << W2.Rows() << "," << W2.Cols() << ")" << std::endl; + //W2.Print(); + /* + tag_W(cell,3*NF,3*NF) = + L + //+ (N*K - L*R)*((N*K-L*R).Transpose()*R).PseudoInvert(1.0e-12)*(N*K-L*R).Transpose() + + (N*K)*((N*K).Transpose()*R).PseudoInvert(1.0e-7)*(N*K).Transpose() + //+ L.Trace()*(rMatrix::Unit(3*NF) - R*((R).Transpose()*R).PseudoInvert(1.0e-14)*(R).Transpose()) + - (L*R)*((L*R).Transpose()*R).PseudoInvert(1.0e-7)*(L*R).Transpose() + ; + */ + //(tag_W(cell,3*NF,3*NF) - W1-W2).Print(); + //tag_W(cell,3*NF,3*NF) = W1+W2; + + + //std::cout << " (NK)^T*R "; + //PrintSV((N*K).Transpose()*R); + + //std::cout << " (LR)^T*R "; + //PrintSV((L*R).Transpose()*R); + + + //std::cout << " vol mat: " << std::endl; + //((N).Transpose()*R).PseudoInvert(1.0e-12).Print(); + + //std::cout << " grad mat: ";// << std::endl; + //PrintSV(((N).Transpose()*R).PseudoInvert(1.0e-12)*(N).Transpose()); + //(((N).Transpose()*R).PseudoInvert(1.0e-12)*(N).Transpose()).Print(); + /* + std::cout << " W1 "; + PrintSV(L); + + L.Print(); + + std::cout << " W2 "; + PrintSV((N*K - L*R)*((N*K).Transpose()*R).PseudoInvert(1.0e-12)*(N*K).Transpose()); + + ((N*K - L*R)*((N*K).Transpose()*R).PseudoInvert(1.0e-12)*(N*K).Transpose()).Print(); + + std::cout << " W "; + PrintSV(tag_W(cell,3*NF,3*NF)); + tag_W(cell,3*NF,3*NF).Print(); + */ + //Convergent on tetra + //tag_W(cell,3*NF,3*NF) = + //L + //+ (N*K-L*R)*((N*K).Transpose()*R).PseudoInvert(1.0e-14)*(N*K).Transpose() + //+ (N*K-L*R)*((N*K-L*R).Transpose()*R).PseudoInvert(1.0e-14)*(N*K-L*R).Transpose() + //+ (L - L*R*((L*R).Transpose()*R).PseudoInvert(1.0e-14)*(L*R).Transpose()) + //+ L.Trace()*(rMatrix::Unit(3*NF) - R*((R).Transpose()*R).PseudoInvert(1.0e-14)*(R).Transpose()) + ; + + //tag_W(cell,3*NF,3*NF) = (N*K)*((N*K).Transpose()*R).PseudoInvert(1.0e-14)*(N*K).Transpose() + (L - L*R*((L*R).Transpose()*R).PseudoInvert(1.0e-14)*(L*R).Transpose()); + + + //tag_W(cell,3*NF,3*NF) = N*K*((N*K).Transpose()*R).PseudoInvert(1.0e-9)*(N*K).Transpose() + (L - L*R*((L*R).Transpose()*R).PseudoInvert(1.0e-9)*(L*R).Transpose()); + + //tag_W(cell,3*NF,3*NF) = L + (N*K-L*R)*(R.Transpose()*R).PseudoInvert(1.0e-9)*R.Transpose(); + //tag_W(cell,3*NF,3*NF) = (N*K)*(N.Transpose()*R).PseudoInvert(1.0e-9)*N.Transpose(); + + //std::cout << "K*V" << std::endl; + //(K).Print(); + //std::cout << "N^T*R" << std::endl; + //(N.Transpose()*R/volume).Print(); + + //std::cout << "R^T*R" << std::endl; + //(R.Transpose()*R).Print(); + //std::cout << "(R^T*R)^{-1}" << std::endl; + //(R.Transpose()*R).PseudoInvert(1.0e-12).Print(); + //std::cout << "W:" < 1.0e-3 ) + { + std::cout << "error: " << std::endl; + (N*K - tag_W(cell,3*NF,3*NF)*R).Print(); + } + + //std::cout << "L "; + //PrintSV(L); + + //std::cout << "W2 "; + //PrintSV((N*K-L*R)*((N*K-L*R).Transpose()*R).PseudoInvert(1.0e-9)*(N*K-L*R).Transpose()); + + //std::cout << "R "; + //PrintSV(rMatrix::Unit(3*NF) - R*(R.Transpose()*R).PseudoInvert(1.0e-9)*R.Transpose()); + + //std::cout << "consistency "; + //PrintSV(NK*(NK.Transpose()*R).PseudoInvert(1.0e-12,&ierr)*NK.Transpose()); + //std::cout << "stability "; + //PrintSV((1.5/(static_cast(NF)*volume)*(NK*K.Invert()*NK.Transpose()).Trace())*(rMatrix::Unit(3*NF) - R*(R.Transpose()*R).PseudoInvert(1.0e-12)*R.Transpose())); + + + //std::cout << "K "; + //PrintSV(K); + if( false ) + { + std::cout << "total "; + PrintSV(tag_W(cell,3*NF,3*NF)); + } + + + //std::cout << "Check:" << std::endl; + //(NK - tag_W(cell,3*NF,3*NF)*R).Print(); + //std::cout << "R^T*R:" << std::endl; + //(R.Transpose()*R).Print(); + //std::cout << "NK^T*R:" << std::endl; + //(NK.Transpose()*R).Print(); + /* + if( ierr ) + { + std::cout << "K:" << std::endl; + K.Print(); + std::cout << "NK:" << std::endl; + NK.Print(); + std::cout << "R:" << std::endl; + R.Print(); + std::cout << "NK^T*R:" << std::endl; + (NK.Transpose()*R).Print(); + std::cout << "sym? " << tag_W(cell,3*NF,3*NF).isSymmetric() << std::endl; + std::cout << "ierr " << ierr << std::endl; + } + */ + } //end of loop over cells + } + std::cout << "Construct W matrix: " << Timer() - ttt << std::endl; + + if( m->HaveTag("UVW") ) //Is there a displacement on the mesh? + tag_UVW = m->GetTag("UVW"); //Get the pressure + + if( !tag_UVW.isValid() || !tag_UVW.isDefined(CELL) ) // Pressure was not initialized or was not defined on nodes + { + //srand(1); // Randomization + tag_UVW = m->CreateTag("UVW",DATA_REAL,CELL|FACE,NONE,3); // Create a new tag for the displacement + for(Mesh::iteratorElement e = m->BeginElement(CELL|FACE); e != m->EndElement(); ++e) //Loop over mesh cells + tag_UVW(*e,3,1).Zero();//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1] + } + + if( !tag_UVW.isDefined(FACE) ) + { + tag_UVW = m->CreateTag("UVW",DATA_REAL,FACE,NONE,3); + for(Mesh::iteratorElement e = m->BeginElement(FACE); e != m->EndElement(); ++e) //Loop over mesh cells + tag_UVW(*e,3,1).Zero();//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1] + } + m->ExchangeData(tag_UVW,CELL|FACE,0); //Synchronize initial solution with boundary unknowns + + + tag_FLUX = m->CreateTag("FLUX",DATA_REAL,FACE,NONE,3); + } //end of initialize data + + + std::cout << "Initialization done" << std::endl; + + + integer nit = 0; + ttt = Timer(); + + { //Main loop for problem solution + Automatizator aut; // declare class to help manage unknowns + Automatizator::MakeCurrent(&aut); + BlockEntry UVW(CELL|FACE); + UVW.AddTag(tag_UVW); + aut.RegisterEntry(UVW); + aut.EnumerateEntries(); //enumerate all primary variables + std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl; + + Residual R("",aut.GetFirstIndex(),aut.GetLastIndex()); + Sparse::LockService Locks(aut.GetFirstIndex(),aut.GetLastIndex()); + Sparse::AnnotationService Text(aut.GetFirstIndex(),aut.GetLastIndex()); + Sparse::Vector Update ("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for update + {//Annotate matrix + for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) + { + Cell cell = m->CellByLocalID(q); + if( cell.GetStatus() != Element::Ghost ) + { + Text.SetAnnotation(UVW.Index(cell,0),"Divergence, U"); + Text.SetAnnotation(UVW.Index(cell,1),"Divergence, V"); + Text.SetAnnotation(UVW.Index(cell,2),"Divergence, W"); + } + } + for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) ) + { + Face face = m->FaceByLocalID(q); + if( face.GetStatus() != Element::Ghost ) + { + if( tag_BC.isValid() && face.HaveData(tag_BC) ) + { + Text.SetAnnotation(UVW.Index(face,0),"Boundary condition, U"); + Text.SetAnnotation(UVW.Index(face,1),"Boundary condition, V"); + Text.SetAnnotation(UVW.Index(face,2),"Boundary condition, W"); + } + else + { + Text.SetAnnotation(UVW.Index(face,0),"Flux continuity, U"); + Text.SetAnnotation(UVW.Index(face,1),"Flux continuity, V"); + Text.SetAnnotation(UVW.Index(face,2),"Flux continuity, W"); + } + } + } + } + + std::cout << "Matrix was annotated" << std::endl; + + do + { + R.Clear(); //clean up the residual + double tttt = Timer(); + int total = 0, dmp = 0; +#if defined(USE_OMP) +#pragma omp parallel +#endif + { + + vMatrix U,T; //vector of pressure differences and fluxes on faces + rMatrix Ra(3,3), Rb(3,3), r(3,1), n(3,1); +#if defined(USE_OMP) +#pragma omp for +#endif + for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) ) + { + Face face = m->FaceByLocalID(q); + tag_FLUX(face,3,1).Zero(); + } +#if defined(USE_OMP) +#pragma omp for +#endif + for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells + { + Cell cell = m->CellByLocalID(q); + ElementArray faces = cell->getFaces(); //obtain faces of the cell + int NF = (int)faces.size(); + + raMatrix W = raMatrixMake(tag_W[cell].data(),3*NF,3*NF); //Matrix for gradient + + U.Resize(3*NF,1); + T.Resize(3*NF,1); + + for(int k = 0; k < NF; ++k) + U(3*k,3*(k+1),0,1) = (UVW[faces[k]] - UVW[cell]); + T = W*U; //fluxes on faces + if( cell.GetStatus() != Element::Ghost ) + { + for(int k = 0; k < NF; ++k) //loop over faces of current cell + R[UVW.Index(cell)] += T(3*k,3*(k+1),0,1); + } + for(int k = 0; k < NF; ++k) //loop over faces of current cell + { + if( faces[k].GetStatus() == Element::Ghost ) continue; + int index = UVW.Index(faces[k],0); + Locks.Lock(index); + if( tag_BC.isValid() && faces[k].HaveData(tag_BC) ) + { + faces[k]->UnitNormal(n.data()); + GetBC(faces[k].RealArray(tag_BC),n,Ra,Rb,r); + if( Rb.FrobeniusNorm() ) Rb.Print(); + R[UVW.Index(faces[k])] -= Ra*UVW[faces[k]] + Rb*T(3*k,3*(k+1),0,1)/faces[k].Area() - r; + + + tag_FLUX(faces[k],3,1) += rMatrix(T(3*k,3*(k+1),0,1))/faces[k].Area(); + } + else + { + R[UVW.Index(faces[k])] -= T(3*k,3*(k+1),0,1); + + tag_FLUX(faces[k],3,1) += 0.5*rMatrix(T(3*k,3*(k+1),0,1))/faces[k].Area()*(faces[k].FaceOrientedOutside(cell)?1:-1); + } + Locks.UnLock(index); + } + } //end of loop over cells + + + if( tag_F.isValid() ) + { +#if defined(USE_OMP) +#pragma omp for +#endif + for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) + { + Cell cell = m->CellByLocalID(q); + if( cell.GetStatus() == Element::Ghost ) continue; + if( cell->HaveData(tag_F) ) R[UVW.Index(cell)] += tag_F(cell,3,1)*cell->Volume(); + } + } + } + std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl; + + std::cout << "Nonlinear residual: " << R.Norm() << "\t\t" << std::endl; + + if( R.Norm() < 1.0e-4 ) break; + tttt = Timer(); + + 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"); + S.SetParameter("drop_tolerance", "1.0e-2"); + S.SetParameter("reuse_tolerance", "1.0e-4"); + + S.SetMatrix(R.GetJacobian()); + + if( S.Solve(R.GetResidual(),Update) ) + { +#if defined(USE_OMP) +#pragma omp parallel for +#endif + for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) + { + Cell cell = m->CellByLocalID(q); + if( cell->GetStatus() == Element::Ghost ) continue; + tag_UVW(cell,3,1) -= Update[UVW.Index(cell)]; + } +#if defined(USE_OMP) +#pragma omp parallel for +#endif + for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) ) + { + Face face = m->FaceByLocalID(q); + if( face->GetStatus() == Element::Ghost ) continue; + tag_UVW(face,3,1) -= Update[UVW.Index(face)]; + } + m->ExchangeData(tag_UVW, CELL|FACE, 0); + + if( print_niter ) + { + std::stringstream str; + str << "iter" << nit; + if( m->GetProcessorsNumber() == 1 ) + str << ".vtk"; + else + str << ".pvtk"; + m->Save(str.str()); + } + } + else + { + std::cout << "Unable to solve: " << S.ReturnReason() << std::endl; + break; + } + std::cout << "solved in " << Timer() - tttt << "\t\t\t" << std::endl; + ++nit; + } while( R.Norm() > 1.0e-4 && nit < 10); //check the residual norm + } + std::cout << "Solved problem in " << Timer() - ttt << " seconds with " << nit << " iterations " << std::endl; + + if( m->HaveTag("REFERENCE_SOLUTION") ) + { + TagReal tag_E = m->CreateTag("ERROR",DATA_REAL,CELL,NONE,1); + TagRealArray tag_R = m->GetTag("REFERENCE_SOLUTION"); + real C, L2, volume; + C = L2 = volume = 0.0; + for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) + { + Cell cell = m->CellByLocalID(q); + real err = (tag_UVW(cell,3,1) - tag_R(cell,3,1)).FrobeniusNorm(); + real vol = cell->Volume(); + if( C < fabs(err) ) C = fabs(err); + L2 += err*err*vol; + volume += vol; + tag_E[cell] = err; + } + C = m->AggregateMax(C); + L2 = m->Integrate(L2); + volume = m->Integrate(volume); + L2 = sqrt(L2/volume); + if( m->GetProcessorRank() == 0 ) std::cout << "Error on cells, C-norm " << C << " L2-norm " << L2 << std::endl; + C = L2 = volume = 0.0; + if( tag_R.isDefined(FACE) ) + { + tag_E = m->CreateTag("ERROR",DATA_REAL,FACE,NONE,1); + for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) ) + { + Face face = m->FaceByLocalID(q); + real err = (tag_UVW(face,3,1) - tag_R(face,3,1)).FrobeniusNorm(); + real vol = (face->BackCell()->Volume() + (face->FrontCell().isValid() ? face->FrontCell()->Volume() : 0))*0.5; + if( C < fabs(err) ) C = fabs(err); + L2 += err*err*vol; + volume += vol; + tag_E[face] = err; + } + C = m->AggregateMax(C); + L2 = m->Integrate(L2); + volume = m->Integrate(volume); + L2 = sqrt(L2/volume); + if( m->GetProcessorRank() == 0 ) std::cout << "Error on faces, C-norm " << C << " L2-norm " << L2 << std::endl; + } + else std::cout << "Reference solution was not defined on faces" << std::endl; + } + + if( m->HaveTag("REFERENCE_FLUX") ) + { + TagReal tag_E = m->CreateTag("ERROR_FLUX",DATA_REAL,FACE,NONE,1); + TagRealArray tag_R = m->GetTag("REFERENCE_FLUX"); + real C, L2, volume; + C = L2 = volume = 0.0; + for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) ) + { + Face face = m->FaceByLocalID(q); + real err = (tag_FLUX(face,3,1) - tag_R(face,3,1)).FrobeniusNorm(); + real vol = face->BackCell()->Volume()/(real)face->BackCell().nbAdjElements(FACE); + if( face->FrontCell().isValid() ) + vol += face->FrontCell()->Volume()/(real)face->FrontCell().nbAdjElements(FACE); + if( C < fabs(err) ) C = fabs(err); + L2 += err*err*vol; + volume += vol; + tag_E[face] = err; + } + C = m->AggregateMax(C); + L2 = m->Integrate(L2); + volume = m->Integrate(volume); + L2 = sqrt(L2/volume); + if( m->GetProcessorRank() == 0 ) std::cout << "Error of fluxes, C-norm " << C << " L2-norm " << L2 << std::endl; + } + + if( m->GetProcessorsNumber() == 1 ) + m->Save("out.vtk"); + else + m->Save("out.pvtk"); + m->Save("out.pmf"); + + delete m; //clean up the mesh + } + else + { + std::cout << argv[0] << " mesh_file" << std::endl; + } + +#if defined(USE_PARTITIONER) + Partitioner::Finalize(); // Finalize the partitioner activity +#endif + Solver::Finalize(); // Finalize solver and close MPI activity + return 0; +} diff --git a/Source/Headers/inmost_dense.h b/Source/Headers/inmost_dense.h index 52b2d42..5c32f79 100644 --- a/Source/Headers/inmost_dense.h +++ b/Source/Headers/inmost_dense.h @@ -699,13 +699,13 @@ namespace INMOST } /// Check if the matrix is symmetric. /// @return Returns true if the matrix is symmetric, otherwise false. - bool isSymmetric() const + bool isSymmetric(double eps = 1.0e-7) const { if( Rows() != Cols() ) return false; for(enumerator k = 0; k < Rows(); ++k) { for(enumerator l = k+1; l < Rows(); ++l) - if( fabs((*this)(k,l)-(*this)(l,k)) > 1.0e-7 ) + if( fabs((*this)(k,l)-(*this)(l,k)) > eps ) return false; } return true; @@ -910,7 +910,7 @@ namespace INMOST /// @param i Column index. /// @param j Row index. /// @return Reference to element. - Var & operator()(enumerator i, enumerator j) + __INLINE Var & operator()(enumerator i, enumerator j) { assert(i >= 0 && i < n); assert(j >= 0 && j < n); @@ -922,7 +922,7 @@ namespace INMOST /// @param i Column index. /// @param j Row index. /// @return Reference to constant element. - const Var & operator()(enumerator i, enumerator j) const + __INLINE const Var & operator()(enumerator i, enumerator j) const { assert(i >= 0 && i < n); assert(j >= 0 && j < n); @@ -932,23 +932,23 @@ namespace INMOST /// Return raw pointer to matrix data, stored in row-wise format. /// @return Pointer to data. - Var * data() {return space.data();} + __INLINE Var * data() {return space.data();} /// Return raw pointer to matrix data without right of change, /// stored in row-wise format. /// @return Pointer to constant data. - const Var * data() const {return space.data();} + __INLINE const Var * data() const {return space.data();} /// Obtain number of rows. /// @return Number of rows. - enumerator Rows() const {return n;} + __INLINE enumerator Rows() const {return n;} /// Obtain number of columns. /// @return Number of columns. - enumerator Cols() const {return n;} + __INLINE enumerator Cols() const {return n;} /// Obtain number of rows. /// @return Reference to number of rows. - enumerator & Rows() {return n;} + __INLINE enumerator & Rows() {return n;} /// Obtain number of rows. /// @return Reference to number of columns. - enumerator & Cols() {return n;} + __INLINE enumerator & Cols() {return n;} /// Convert values in array into square matrix. /// Supports the following representation, depending on the size /// of input array and size of side of final tensors' matrix: @@ -1348,18 +1348,18 @@ namespace INMOST { if( Cols()*Rows() != other.Cols()*other.Rows() ) space.resize(other.Cols()*other.Rows()); + n = other.Rows(); + m = other.Cols(); for(enumerator i = 0; i < other.Rows(); ++i) for(enumerator j = 0; j < other.Cols(); ++j) assign((*this)(i,j),other(i,j)); - n = other.Rows(); - m = other.Cols(); return *this; } /// Access element of the matrix by row and column indices. /// @param i Column index. /// @param j Row index. /// @return Reference to element. - Var & operator()(enumerator i, enumerator j) + __INLINE Var & operator()(enumerator i, enumerator j) { assert(i >= 0 && i < n); assert(j >= 0 && j < m); @@ -1371,7 +1371,7 @@ namespace INMOST /// @param i Column index. /// @param j Row index. /// @return Reference to constant element. - const Var & operator()(enumerator i, enumerator j) const + __INLINE const Var & operator()(enumerator i, enumerator j) const { assert(i >= 0 && i < n); assert(j >= 0 && j < m); @@ -1381,23 +1381,23 @@ namespace INMOST /// Return raw pointer to matrix data, stored in row-wise format. /// @return Pointer to data. - Var * data() {return space.data();} + __INLINE Var * data() {return space.data();} /// Return raw pointer to matrix data without right of change, /// stored in row-wise format. /// @return Pointer to constant data. - const Var * data() const {return space.data();} + __INLINE const Var * data() const {return space.data();} /// Obtain number of rows. /// @return Number of rows. - enumerator Rows() const {return n;} + __INLINE enumerator Rows() const {return n;} /// Obtain number of columns. /// @return Number of columns. - enumerator Cols() const {return m;} + __INLINE enumerator Cols() const {return m;} /// Obtain number of rows. /// @return Reference to number of rows. - enumerator & Rows() {return n;} + __INLINE enumerator & Rows() {return n;} /// Obtain number of rows. /// @return Reference to number of columns. - enumerator & Cols() {return m;} + __INLINE enumerator & Cols() {return m;} /// Construct row permutation matrix from array of new positions for rows. /// Row permutation matrix multiplies matrix from left. /// Column permutation matrix is obtained by transposition and is multiplied @@ -1793,10 +1793,10 @@ namespace INMOST public: /// Number of rows in submatrix. /// @return Number of rows. - enumerator Rows() const {return erow-brow;} + __INLINE enumerator Rows() const {return erow-brow;} /// Number of columns in submatrix. /// @return Number of columns. - enumerator Cols() const {return ecol-bcol;} + __INLINE enumerator Cols() const {return ecol-bcol;} /// Create submatrix for a matrix. /// @param rM Reference to the matrix that stores elements. /// @param first_row First row in the matrix. @@ -1836,7 +1836,7 @@ namespace INMOST /// @param i Column index. /// @param j Row index. /// @return Reference to element. - Var & operator()(enumerator i, enumerator j) + __INLINE Var & operator()(enumerator i, enumerator j) { assert(i >= 0 && i < Rows()); assert(j >= 0 && j < Cols()); @@ -1848,7 +1848,7 @@ namespace INMOST /// @param i Column index. /// @param j Row index. /// @return Reference to constant element. - const Var & operator()(enumerator i, enumerator j) const + __INLINE const Var & operator()(enumerator i, enumerator j) const { assert(i >= 0 && i < Rows()); assert(j >= 0 && j < Cols()); @@ -1894,10 +1894,10 @@ namespace INMOST public: /// Number of rows in submatrix. /// @return Number of rows. - enumerator Rows() const {return erow-brow;} + __INLINE enumerator Rows() const {return erow-brow;} /// Number of columns in submatrix. /// @return Number of columns. - enumerator Cols() const {return ecol-bcol;} + __INLINE enumerator Cols() const {return ecol-bcol;} /// Create submatrix for a matrix. /// @param rM Reference to the matrix that stores elements. /// @param first_row First row in the matrix. @@ -1912,7 +1912,7 @@ namespace INMOST /// @param i Column index. /// @param j Row index. /// @return Reference to element. - Var & operator()(enumerator i, enumerator j) + __INLINE Var & operator()(enumerator i, enumerator j) { assert(i >= 0 && i < Rows()); assert(j >= 0 && j < Cols()); @@ -1924,7 +1924,7 @@ namespace INMOST /// @param i Column index. /// @param j Row index. /// @return Reference to constant element. - const Var & operator()(enumerator i, enumerator j) const + __INLINE const Var & operator()(enumerator i, enumerator j) const { assert(i >= 0 && i < Rows()); assert(j >= 0 && j < Cols()); @@ -2463,7 +2463,14 @@ namespace INMOST { if( ierr ) { - if( *ierr == -1 ) std::cout << "Failed to invert matrix" << std::endl; + if( *ierr == -1 ) + { + std::cout << "Failed to invert matrix diag " << get_value(AtA(i,i)) << std::endl; + std::cout << "rhs:"; + for(enumerator k = 0; k < l; k++) + std::cout << " " << get_value(AtB(i,k)); + std::cout << std::endl; + } *ierr = i+1; } else throw MatrixSolveFail; -- 2.26.2