Commit 3cc9d06d authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Adapt library and examples for 64-bit integer option (USE_INT64=ON); adapt...

Adapt library and examples for 64-bit integer option (USE_INT64=ON); adapt library (but not yet examples) for 32-bit float option (USE_FP64=OFF); make gcc happy on library compilation  with -Wall -Wextra -pedantic; replace malloc/realloc/free with new/delete in container.hpp
parent b5445e4c
...@@ -40,6 +40,8 @@ option(USE_SOLVER_PETSC "Use PETSc solvers" OFF) ...@@ -40,6 +40,8 @@ option(USE_SOLVER_PETSC "Use PETSc solvers" OFF)
option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF) option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF)
option(USE_SOLVER_SUPERLU "Use SuperLU solver" OFF) option(USE_SOLVER_SUPERLU "Use SuperLU solver" OFF)
option(USE_SOLVER_SUPERLU_DOWNLOAD "Attempts to download SuperLU solver if library not found" OFF) option(USE_SOLVER_SUPERLU_DOWNLOAD "Attempts to download SuperLU solver if library not found" OFF)
option(USE_SOLVER_K3BIILU2 "Activate K3BIILU2 solver (requires USE_MPI and USE_SOLVER_METIS)")
option(USE_SOLVER_FCBIILU2 "Activate FCBIILU2 solver (requires fcbiilu2.cpp)")
option(USE_OPTIMIZER "Use Optimization toolkit" OFF) option(USE_OPTIMIZER "Use Optimization toolkit" OFF)
option(USE_OPTIMIZER_BAYESIAN "Compile Optimization module against Limbo library to provide Bayesian optimization algorithm" OFF) option(USE_OPTIMIZER_BAYESIAN "Compile Optimization module against Limbo library to provide Bayesian optimization algorithm" OFF)
#option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF) #option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF)
......
...@@ -48,6 +48,7 @@ int main(int argc,char ** argv) ...@@ -48,6 +48,7 @@ int main(int argc,char ** argv)
Mesh * m = new Mesh(); // Create an empty mesh Mesh * m = new Mesh(); // Create an empty mesh
double ttt = Timer(); double ttt = Timer();
bool repartition = false; bool repartition = false;
(void)repartition;
m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
if( m->GetProcessorRank() == 0 ) // If the current process is the master one if( m->GetProcessorRank() == 0 ) // If the current process is the master one
std::cout << argv[0] << std::endl; std::cout << argv[0] << std::endl;
......
project(ADMFD) project(ADMFD)
add_executable(MFDDIFF diffusion.cpp) add_executable(MFDDIFF diffusion.cpp)
add_executable(NMFDDIFF diffusion_nonlinear.cpp)
add_executable(MFDELAST elastic.cpp) add_executable(MFDELAST elastic.cpp)
target_link_libraries(MFDDIFF inmost) target_link_libraries(MFDDIFF inmost)
target_link_libraries(NMFDDIFF inmost)
target_link_libraries(MFDELAST inmost) target_link_libraries(MFDELAST inmost)
if(USE_SOLVER) if(USE_SOLVER)
...@@ -56,6 +58,58 @@ if(USE_MPI) ...@@ -56,6 +58,58 @@ if(USE_MPI)
endif() endif()
endif(USE_MPI) endif(USE_MPI)
if(USE_SOLVER)
if(USE_SOLVER_ANI)
message("linking NMFDDIFF with ani3d and BLAS")
target_link_libraries(NMFDDIFF ani3d ${BLAS_LIBRARIES})
if(BLAS_LINKER_FLAGS)
set_target_properties(NMFDDIFF PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}")
endif()
endif()
if(USE_SOLVER_PETSC)
message("linking NMFDDIFF with PETSc")
target_link_libraries(NMFDDIFF ${PETSC_LIBRARIES})
endif()
if(USE_SOLVER_TRILINOS)
message("linking NMFDDIFF with Trilinos")
target_link_libraries(NMFDDIFF ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_SOLVER_METIS)
message("linking NMFDDIFF with Metis")
target_link_libraries(NMFDDIFF ${METIS_LIBRARIES})
endif()
if(USE_SOLVER_MONDRIAAN)
message("linking NMFDDIFF with Mondriaan")
target_link_libraries(NMFDDIFF ${MONDRIAAN_LIBRARIES})
endif()
if(USE_SOLVER_SUPERLU)
message("linking NMFDDIFF with SuperLU")
target_link_libraries(NMFDDIFF ${SUPERLU_LIBRARIES})
endif()
endif()
if(USE_PARTITIONER)
if(USE_PARTITIONER_ZOLTAN)
message("linking NMFDDIFF with Zoltan")
target_link_libraries(NMFDDIFF ${ZOLTAN_LIBRARIES})
endif()
if(USE_PARTITIONER_PARMETIS)
message("linking NMFDDIFF with ParMETIS")
target_link_libraries(NMFDDIFF ${PARMETIS_LIBRARIES})
endif()
endif()
if(USE_MPI)
message("linking NMFDDIFF with MPI")
target_link_libraries(NMFDDIFF ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(NMFDDIFF PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
if(USE_SOLVER) if(USE_SOLVER)
if(USE_SOLVER_ANI) if(USE_SOLVER_ANI)
message("linking MFDELAST with ani3d and BLAS") message("linking MFDELAST with ani3d and BLAS")
...@@ -107,4 +161,5 @@ if(USE_MPI) ...@@ -107,4 +161,5 @@ if(USE_MPI)
endif(USE_MPI) endif(USE_MPI)
install(TARGETS MFDDIFF EXPORT inmost-targets RUNTIME DESTINATION bin) install(TARGETS MFDDIFF EXPORT inmost-targets RUNTIME DESTINATION bin)
install(TARGETS NMFDDIFF EXPORT inmost-targets RUNTIME DESTINATION bin)
install(TARGETS MFDELAST EXPORT inmost-targets RUNTIME DESTINATION bin) install(TARGETS MFDELAST EXPORT inmost-targets RUNTIME DESTINATION bin)
...@@ -205,18 +205,15 @@ int main(int argc,char ** argv) ...@@ -205,18 +205,15 @@ int main(int argc,char ** argv)
rMatrix NK, L, R, Areas; rMatrix NK, L, R, Areas;
rMatrix x(1,3), xf(1,3), n(1,3); rMatrix x(1,3), xf(1,3), n(1,3);
double area, dist; //area of the face double area, dist; //area of the face
double volume; //volume of the cell
#if defined(USE_OMP) #if defined(USE_OMP)
#pragma omp for #pragma omp for
#endif #endif
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{ {
Mesh * mesh = m;
Cell cell = m->CellByLocalID(q); Cell cell = m->CellByLocalID(q);
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
int NF = (int)faces.size(); //number of faces; int NF = (int)faces.size(); //number of faces;
rMatrix W(NF,NF); rMatrix W(NF,NF);
volume = cell->Volume(); //volume of the cell
cell->Barycenter(x.data()); cell->Barycenter(x.data());
//get permeability for the cell //get permeability for the cell
rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(), rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),
...@@ -349,7 +346,7 @@ int main(int argc,char ** argv) ...@@ -349,7 +346,7 @@ int main(int argc,char ** argv)
{ {
R.Clear(); //clean up the residual R.Clear(); //clean up the residual
double tttt = Timer(); double tttt = Timer();
int total = 0, dmp = 0; //~ int total = 0, dmp = 0;
#if defined(USE_OMP) #if defined(USE_OMP)
#pragma omp parallel #pragma omp parallel
#endif #endif
......
...@@ -108,8 +108,8 @@ int main(int argc,char ** argv) ...@@ -108,8 +108,8 @@ int main(int argc,char ** argv)
TagRealArray tag_PG; // Pressure gradient //~ TagRealArray tag_PG; // Pressure gradient
TagRealArray tag_WG; // matrix to reconstruct gradient //~ TagRealArray tag_WG; // matrix to reconstruct gradient
if( m->GetProcessorsNumber() > 1 ) //skip for one processor job if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
{ // Exchange ghost cells { // Exchange ghost cells
...@@ -145,8 +145,8 @@ int main(int argc,char ** argv) ...@@ -145,8 +145,8 @@ int main(int argc,char ** argv)
if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh? if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh?
tag_P = m->GetTag("PRESSURE"); //Get the pressure tag_P = m->GetTag("PRESSURE"); //Get the pressure
tag_PG = m->CreateTag("PRESSURE_GRADIENT",DATA_REAL,CELL,NONE,3); //~ tag_PG = m->CreateTag("PRESSURE_GRADIENT",DATA_REAL,CELL,NONE,3);
tag_WG = m->CreateTag("WGRAD",DATA_REAL,CELL,NONE); //~ tag_WG = m->CreateTag("WGRAD",DATA_REAL,CELL,NONE);
if( !tag_P.isValid() || !tag_P.isDefined(CELL) ) // Pressure was not initialized or was not defined on nodes if( !tag_P.isValid() || !tag_P.isDefined(CELL) ) // Pressure was not initialized or was not defined on nodes
{ {
...@@ -187,9 +187,8 @@ int main(int argc,char ** argv) ...@@ -187,9 +187,8 @@ int main(int argc,char ** argv)
#if defined(USE_OMP) #if defined(USE_OMP)
#pragma omp for #pragma omp for
#endif #endif
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) for( integer q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{ {
Mesh * mesh = m;
Cell cell = m->CellByLocalID(q); Cell cell = m->CellByLocalID(q);
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
int NF = (int)faces.size(); //number of faces; int NF = (int)faces.size(); //number of faces;
...@@ -201,28 +200,46 @@ int main(int argc,char ** argv) ...@@ -201,28 +200,46 @@ int main(int argc,char ** argv)
R.Resize(NF,3); //directions R.Resize(NF,3); //directions
L.Resize(NF,NF); //two-point transmissibility coefficient L.Resize(NF,NF); //two-point transmissibility coefficient
L.Zero(); L.Zero();
tag_WS[cell].resize(NF);
for(int k = 0; k < NF; ++k) //loop over faces for(int k = 0; k < NF; ++k) //loop over faces
{ {
area = faces[k].Area(); area = faces[k].Area();
faces[k].Barycenter(xf.data()); faces[k].Barycenter(xf.data());
faces[k].OrientedUnitNormal(cell->self(),n.data()); faces[k].OrientedUnitNormal(cell->self(),n.data());
dist = n.DotProduct(xf-x); dist = (n.DotProduct(xf-x));
if( dist < 0 ) std::cout << "negative dist! " << dist << " cell " << q << " face " << k << std::endl;
// assemble matrix of directions // assemble matrix of directions
R(k,k+1,0,3) = (xf-x); R(k,k+1,0,3) = (xf-x);
// assemble matrix of co-normals // assemble matrix of co-normals
N(k,k+1,0,3) = n*K*area; N(k,k+1,0,3) = n*K*area;
L(k,k) = n.DotProduct(n*K)*area/dist; L(k,k) = n.DotProduct(n*K)*area/dist;
tag_WS[cell][k] = L(k,k);
} //end of loop over faces } //end of loop over faces
double ca = 0, ba = ca; //parameters //~ double ca = 0, ba = ca; //parameters
tag_WA[cell].resize(NF*NF); tag_WA[cell].resize(NF*NF);
tag_WS[cell].resize(NF*NF);
tag_WG[cell].resize(3*NF); //~ tag_WG[cell].resize(3*NF);
//tag_WA(cell,NF,NF) = (N+ca*L*R)*((N+ca*L*R).Transpose()*R).Invert()*(N+ca*L*R).Transpose(); //tag_WA(cell,NF,NF) = (N+ca*L*R)*((N+ca*L*R).Transpose()*R).Invert()*(N+ca*L*R).Transpose();
//tag_WS(cell,NF,NF) = L - (1+ba)*(L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose(); //tag_WS(cell,NF,NF) = L - (1+ba)*(L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
tag_WA(cell,NF,NF) = N*(N.Transpose()*R).Invert()*N.Transpose(); tag_WA(cell,NF,NF) = N*(N.Transpose()*R).Invert()*N.Transpose() - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
tag_WS(cell,NF,NF) = L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose(); /*
tag_WG(cell,3,NF) = (N.Transpose()*R).Invert()*N.Transpose(); #pragma omp critical
{
std::cout << "cell " << cell.LocalID() << std::endl;
std::cout << "WA: " << std::endl;
tag_WA(cell,NF,NF).Print();
for(int k = 0; k < NF; ++k)
{
double rowsum = 0;
for(int q = 0; q < NF; ++q)
rowsum += tag_WA(cell,NF,NF)(k,q);
std::cout << k << " T " << L(k,k) << " rowsum " << rowsum+L(k,k) << " no T " << rowsum << std::endl;
}
}
*/
//~ tag_WS(cell,NF,1) = L;
//~ tag_WG(cell,3,NF) = (N.Transpose()*R).Invert()*N.Transpose();
} //end of loop over cells } //end of loop over cells
} }
std::cout << "Construct W matrix: " << Timer() - ttt << std::endl; std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
...@@ -277,18 +294,17 @@ int main(int argc,char ** argv) ...@@ -277,18 +294,17 @@ int main(int argc,char ** argv)
do do
{ {
Resid.Clear(); //clean up the residual Resid.Clear(); //clean up the residual
double tttt = Timer(); //~ double tttt = Timer();
int total = 0, dmp = 0;
#if defined(USE_OMP) #if defined(USE_OMP)
#pragma omp parallel #pragma omp parallel
#endif #endif
{ {
vMatrix pF; //vector of pressure differences on faces vMatrix pF; //vector of pressure differences on faces
vMatrix FLUX; //computed flux on faces vMatrix FLUX; //computed flux on faces
vMatrix N; variable u, v;
rMatrix R; variable muu, muv;
rMatrix x(1,3), xf(1,3), n(1,3), coefs, fluxes; //~ double muu,muv;
double area; //area of the face //~ double area; //area of the face
#if defined(USE_OMP) #if defined(USE_OMP)
#pragma omp for #pragma omp for
#endif #endif
...@@ -301,7 +317,7 @@ int main(int argc,char ** argv) ...@@ -301,7 +317,7 @@ int main(int argc,char ** argv)
for(int k = 0; k < NF; ++k) for(int k = 0; k < NF; ++k)
pF(k,0) = (P(faces[k]) - P(cell)); pF(k,0) = (P(faces[k]) - P(cell));
FLUX = tag_WA(cell,NF,NF)*pF; //approximation of fluxes on faces FLUX = tag_WA(cell,NF,NF)*pF; //approximation of fluxes on faces
tag_PG(cell,3,1) = tag_WG(cell,3,NF)*pF; //gradient of function //~ tag_PG(cell,3,1) = tag_WG(cell,3,NF)*pF; //gradient of function
for(int k = 0; k < NF; ++k) //loop over faces of current cell for(int k = 0; k < NF; ++k) //loop over faces of current cell
{ {
int ind = (faces[k].BackCell() == cell ? 0 : 1); int ind = (faces[k].BackCell() == cell ? 0 : 1);
...@@ -318,6 +334,7 @@ int main(int argc,char ** argv) ...@@ -318,6 +334,7 @@ int main(int argc,char ** argv)
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
int NF = (int)faces.size(); int NF = (int)faces.size();
//get permeability for the cell //get permeability for the cell
/*
cell->Barycenter(x.data()); cell->Barycenter(x.data());
rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(), rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),
cell->RealArrayDF(tag_K).size()); cell->RealArrayDF(tag_K).size());
...@@ -336,6 +353,7 @@ int main(int argc,char ** argv) ...@@ -336,6 +353,7 @@ int main(int argc,char ** argv)
N(k,k+1,0,3) = n*K*area; N(k,k+1,0,3) = n*K*area;
int ind = (faces[k].BackCell() == cell ? 0 : 1); int ind = (faces[k].BackCell() == cell ? 0 : 1);
FLUX(k,0); = tag_WF[faces[k]][ind] //restore fluxes
if( !faces[k].Boundary() ) if( !faces[k].Boundary() )
{ {
const double eps = 1.0e-6; const double eps = 1.0e-6;
...@@ -355,26 +373,35 @@ int main(int argc,char ** argv) ...@@ -355,26 +373,35 @@ int main(int argc,char ** argv)
coefs(k,0) = 1; coefs(k,0) = 1;
} }
} //end of loop over faces } //end of loop over faces
pF.Resize(NF,1); */
for(int k = 0; k < NF; ++k) for(int k = 0; k < NF; ++k)
pF(k,0) = (P(faces[k]) - P(cell));
if( false ) if( nit != 0 )
{ {
std::cout << "fluxes: " << std::endl; int ind = (faces[k].BackCell() == cell ? 0 : 1);
fluxes.Transpose().Print(); //FLUX(k,0) = tag_WF[faces[k]][ind];
std::cout << "coefs: " << std::endl; u = tag_WF[faces[k]][ind];
coefs.Transpose().Print(); if( !faces[k].Boundary() )
std::cout << " inv(N^T*R) " << std::endl; {
(N.Transpose()*R).Invert().Print(); v = tag_WF[faces[k]][1-ind];
std::cout << "Approximation W" << std::endl; //~ muu = fabs(get_value(v)) + 1.0e-5;
(N*(N.Transpose()*R).Invert()*N.Transpose()).Print(); //~ muv = fabs(get_value(u)) + 1.0e-5;
std::cout << "Approximation W (original)" << std::endl; //~ muu = get_value(v)*get_value(v) + 1.0e-5;
tag_WA(cell,NF,NF).Print(); //~ muv = get_value(u)*get_value(u) + 1.0e-5;
std::cout << "Stabilization W" << std::endl; //~ muu = 0.5;
tag_WS(cell,NF,NF).Print(); //~ muv = 0.5;
scanf("%*c"); muu = sqrt(v*v + 1.0e-4);
muv = sqrt(u*u + 1.0e-4);
//~ muu = variation(sqrt(v*v + 1.0e-10),0.25);
//~ muv = variation(sqrt(u*u + 1.0e-10),0.25);
//~ muu = v*v + 1.0e-2;
//~ muv = u*u + 1.0e-2;
//~ FLUX(k,0) = u*variation(muu/(muu+muv),0.1) - v*variation(muv/(muu+muv),0.1);
FLUX(k,0) = (u*muu - v*muv)/(muu+muv);
//~ FLUX(k,0) = 2*u*variation(muu/(muu+muv),0.5);
}
else FLUX(k,0) = u;
} }
FLUX = (N*(N.Transpose()*R).Invert()*N.Transpose() + tag_WS(cell,NF,NF))*pF; //approximation of fluxes on faces for(int k = 0; k < NF; ++k)
FLUX(k,0) += tag_WS[cell][k]*(P(faces[k]) - P(cell));
if( cell.GetStatus() != Element::Ghost ) if( cell.GetStatus() != Element::Ghost )
{ {
for(int k = 0; k < NF; ++k) //loop over faces of current cell for(int k = 0; k < NF; ++k) //loop over faces of current cell
...@@ -391,7 +418,10 @@ int main(int argc,char ** argv) ...@@ -391,7 +418,10 @@ int main(int argc,char ** argv)
Resid[index] -= BC[0]*P(faces[k]) + BC[1]*FLUX(k,0) - BC[2]; Resid[index] -= BC[0]*P(faces[k]) + BC[1]*FLUX(k,0) - BC[2];
} }
else else
{
Resid[index] -= FLUX(k,0); Resid[index] -= FLUX(k,0);
}
Locks.UnLock(index); Locks.UnLock(index);
} }
} //end of loop over cells } //end of loop over cells
...@@ -412,9 +442,9 @@ int main(int argc,char ** argv) ...@@ -412,9 +442,9 @@ int main(int argc,char ** argv)
//R[P.Index(m->BeginCell()->self())] = P[m->BeginCell()->self()]; //R[P.Index(m->BeginCell()->self())] = P[m->BeginCell()->self()];
} }
std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl; //~ std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
std::cout << "Nonlinear residual: " << Resid.Norm() << "\t\t" << std::endl; std::cout << "Nonlinear residual: " << Resid.Norm() << " " << std::endl;
if( Resid.Norm() < 1.0e-4 ) break; if( Resid.Norm() < 1.0e-4 ) break;
...@@ -425,14 +455,14 @@ int main(int argc,char ** argv) ...@@ -425,14 +455,14 @@ int main(int argc,char ** argv)
S.SetParameter("verbosity","1"); S.SetParameter("verbosity","1");
S.SetParameter("relative_tolerance", "1.0e-14"); S.SetParameter("relative_tolerance", "1.0e-14");
S.SetParameter("absolute_tolerance", "1.0e-12"); S.SetParameter("absolute_tolerance", "1.0e-12");
S.SetParameter("drop_tolerance", "1.0e-4"); S.SetParameter("drop_tolerance", "1.0e-3");
S.SetParameter("reuse_tolerance", "1.0e-5"); S.SetParameter("reuse_tolerance", "1.0e-4");
S.SetMatrix(Resid.GetJacobian()); S.SetMatrix(Resid.GetJacobian());
if( S.Solve(Resid.GetResidual(),Update) ) if( S.Solve(Resid.GetResidual(),Update) )
{ {
double alpha = 0.5; double alpha = 1;
#if defined(USE_OMP) #if defined(USE_OMP)
#pragma omp parallel for #pragma omp parallel for
#endif #endif
...@@ -524,6 +554,7 @@ int main(int argc,char ** argv) ...@@ -524,6 +554,7 @@ int main(int argc,char ** argv)
m->Save("out.pvtk"); m->Save("out.pvtk");
m->Save("solution.pmf"); m->Save("solution.pmf");
m->Save("solution.vtu");
delete m; //clean up the mesh delete m; //clean up the mesh
} }
......
...@@ -41,15 +41,15 @@ bool print_matrix = false; ...@@ -41,15 +41,15 @@ bool print_matrix = false;
template<typename T> template<typename T>
static void SVD2Eigen(const Matrix<T> & U, Matrix<T> & S, Matrix<T> & V) static void SVD2Eigen(const Matrix<T> & U, Matrix<T> & S, Matrix<T> & V)
{ {
for (int i = 0; i < V.Cols(); ++i) for (unsigned i = 0; i < V.Cols(); ++i)
{ {
double dot = 0.0; double dot = 0.0;
for (int j = 0; j < V.Rows(); ++j) for (unsigned j = 0; j < V.Rows(); ++j)
dot += get_value(U(j, i))*get_value(V(j, i)); dot += get_value(U(j, i))*get_value(V(j, i));
if (dot < 0.0) if (dot < 0.0)
{ {
S(i, i) *= -1; S(i, i) *= -1;
for (int j = 0; j < V.Rows(); ++j) for (unsigned j = 0; j < V.Rows(); ++j)
V(j, i) *= -1; V(j, i) *= -1;
} }
} }
...@@ -290,7 +290,7 @@ void PrintSV(const rMatrix & A) ...@@ -290,7 +290,7 @@ void PrintSV(const rMatrix & A)
A.SVD(U,S,V); A.SVD(U,S,V);
std::cout << "singular values:"; std::cout << "singular values:";
int cnt = 0; int cnt = 0;
for(int k = 0; k < A.Cols(); ++k) if( fabs(S(k,k)) > 1.0e-10 ) for(unsigned k = 0; k < A.Cols(); ++k) if( fabs(S(k,k)) > 1.0e-10 )
{ {
std::cout << " " << S(k,k); std::cout << " " << S(k,k);
cnt++; cnt++;
...@@ -304,10 +304,10 @@ void PrintRS(const rMatrix & A) ...@@ -304,10 +304,10 @@ void PrintRS(const rMatrix & A)
{ {
double sum; double sum;
std::cout << "row sum:"; std::cout << "row sum:";
for(int i = 0; i < A.Rows(); ++i) for(unsigned i = 0; i < A.Rows(); ++i)
{ {
sum = 0; sum = 0;
for(int j = 0; j < A.Cols(); ++j) for(unsigned j = 0; j < A.Cols(); ++j)
sum+= A(i,j); sum+= A(i,j);
std::cout << " " << sum; std::cout << " " << sum;
} }
...@@ -513,18 +513,17 @@ int main(int argc,char ** argv) ...@@ -513,18 +513,17 @@ int main(int argc,char ** argv)
rMatrix N, R, L, M(9,9), T, K(9,9), iK(9,9), C(6,6), Q, W1, W2, W3, W3s, U,S,V, w, u, v; rMatrix N, R, L, M(9,9), T, K(9,9), iK(9,9), C(6,6), Q, W1, W2, W3, W3s, U,S,V, w, u, v;
rMatrix x(3,1), xf(3,1), n(3,1); rMatrix x(3,1), xf(3,1), n(3,1);
double area; //area of the face double area; //area of the face
double volume; //volume of the cell //double volume; //volume of the cell
double dist; //distance from cell center to face double dist; //distance from cell center to face
#if defined(USE_OMP) #if defined(USE_OMP)
#pragma omp for #pragma omp for
#endif #endif
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{ {
Mesh * mesh = m;
Cell cell = m->CellByLocalID(q); Cell cell = m->CellByLocalID(q);
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
int NF = (int)faces.size(); //number of faces; int NF = (int)faces.size(); //number of faces;
volume = cell->Volume(); //volume of the cell //~ volume = cell->Volume(); //volume of the cell
cell->Barycenter(x.data());; cell->Barycenter(x.data());;
//get permeability for the cell //get permeability for the cell
KTensor(tag_C[cell],K); KTensor(tag_C[cell],K);
...@@ -628,7 +627,6 @@ int main(int argc,char ** argv) ...@@ -628,7 +627,6 @@ int main(int argc,char ** argv)
#endif #endif
for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) ) for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
{ {
Mesh * mesh = m;
Face f = m->FaceByLocalID(q); Face f = m->FaceByLocalID(q);
Storage::integer_array rows = f->IntegerArray(tag_Row); Storage::integer_array rows = f->IntegerArray(tag_Row);
...@@ -649,7 +647,7 @@ int main(int argc,char ** argv) ...@@ -649,7 +647,7 @@ int main(int argc,char ** argv)
tag_i_rhs(f,3,1).Zero(); tag_i_rhs(f,3,1).Zero();
real aF = 1; real aF = 1;
//f_rhs = i_rhs = 0; //f_rhs = i_rhs = 0;