Commit 8fa00d6a authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Solver updates

Added support of Trilinos: AztecOO, ML, Ifpack, Belos. Added Petsc
divergence reasons. Described parameters, some parameters unified with
Trilinos.
parent d7217797
...@@ -59,7 +59,8 @@ option(COMPILE_TESTS "Compile some tests" OFF) ...@@ -59,7 +59,8 @@ option(COMPILE_TESTS "Compile some tests" OFF)
option(USE_PARTITIONER_PARMETIS "Use ParMetis partitioner" OFF) option(USE_PARTITIONER_PARMETIS "Use ParMetis partitioner" OFF)
option(USE_PARTITIONER_ZOLTAN "Use Zoltan partitioner" OFF) option(USE_PARTITIONER_ZOLTAN "Use Zoltan partitioner" OFF)
option(USE_SOLVER_PETSC "Use PETSc solver" OFF) option(USE_SOLVER_PETSC "Use PETSc solvers" OFF)
option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF)
option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF) option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF)
option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF) option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF)
option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates for automatic differentiation" OFF) option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates for automatic differentiation" OFF)
...@@ -142,6 +143,40 @@ if(USE_SOLVER_PETSC) ...@@ -142,6 +143,40 @@ if(USE_SOLVER_PETSC)
endif() endif()
if(USE_SOLVER_TRILINOS)
set(CMAKE_PREFIX_PATH ${TRILINOS_PATH} ${CMAKE_PREFIX_PATH})
find_package(Trilinos PATHS ${TRILINOS_PATH}/lib/cmake/Trilinos ${TRILINOS_PATH})
if(NOT Trilinos_FOUND)
set(USE_SOLVER_TRILINOS OFF)
message("Trilinos NOT FOUND")
else()
MESSAGE("\nFound Trilinos! Here are the details: ")
MESSAGE(" Trilinos_DIR = ${Trilinos_DIR}")
MESSAGE(" Trilinos_VERSION = ${Trilinos_VERSION}")
MESSAGE(" Trilinos_PACKAGE_LIST = ${Trilinos_PACKAGE_LIST}")
MESSAGE(" Trilinos_LIBRARIES = ${Trilinos_LIBRARIES}")
MESSAGE(" Trilinos_INCLUDE_DIRS = ${Trilinos_INCLUDE_DIRS}")
MESSAGE(" Trilinos_LIBRARY_DIRS = ${Trilinos_LIBRARY_DIRS}")
MESSAGE(" Trilinos_TPL_LIST = ${Trilinos_TPL_LIST}")
MESSAGE(" Trilinos_TPL_INCLUDE_DIRS = ${Trilinos_TPL_INCLUDE_DIRS}")
MESSAGE(" Trilinos_TPL_LIBRARIES = ${Trilinos_TPL_LIBRARIES}")
MESSAGE(" Trilinos_TPL_LIBRARY_DIRS = ${Trilinos_TPL_LIBRARY_DIRS}")
MESSAGE(" Trilinos_BUILD_SHARED_LIBS = ${Trilinos_BUILD_SHARED_LIBS}")
MESSAGE("End of Trilinos details\n")
include_directories(${Trilinos_INCLUDE_DIRS})
include_directories(${Trilinos_TPL_INCLUDE_DIRS})
link_directories(${Trilinos_LIBRARY_DIRS})
link_directories(${Trilinos_TPL_LIBRARY_DIRS})
set(USE_SOLVER_TRILINOS ON)
message("Trilinos FOUND")
endif()
endif()
if(USE_AUTODIFF_OPENCL) if(USE_AUTODIFF_OPENCL)
find_package(OpenCL) find_package(OpenCL)
if(OPENCL_FOUND) if(OPENCL_FOUND)
...@@ -167,6 +202,14 @@ if(USE_AUTODIFF_EXPRESSION_TEMPLATES) ...@@ -167,6 +202,14 @@ if(USE_AUTODIFF_EXPRESSION_TEMPLATES)
endif() endif()
if(MSVC)
if(USE_SOLVER_TRILINOS)
message("Putting workaround for Visual Studio that allow to use Trilinos Release libraries in Debug mode")
message("Note that this workaround may affect your debugging experience, you may want to debug without Trilinos")
add_definitions(-D_ITERATOR_DEBUG_LEVEL=0)
endif()
endif()
configure_file("inmost_options_cmake.h" "${PROJECT_BINARY_DIR}/inmost_options.h") configure_file("inmost_options_cmake.h" "${PROJECT_BINARY_DIR}/inmost_options.h")
include_directories("${PROJECT_BINARY_DIR}") include_directories("${PROJECT_BINARY_DIR}")
......
...@@ -5,4 +5,4 @@ add_subdirectory(MatSolve) ...@@ -5,4 +5,4 @@ add_subdirectory(MatSolve)
add_subdirectory(GridGen) add_subdirectory(GridGen)
add_subdirectory(FVDiscr) add_subdirectory(FVDiscr)
#add_subdirectory(OctreeCutcell) #add_subdirectory(OctreeCutcell)
#add_subdirectory(Solver) add_subdirectory(Solver)
...@@ -4,28 +4,40 @@ find_package(OpenGL) ...@@ -4,28 +4,40 @@ find_package(OpenGL)
find_package(GLUT) find_package(GLUT)
if(OPENGL_FOUND) if(OPENGL_FOUND)
if(GLUT_FOUND) if(GLUT_FOUND)
include_directories(${OPENGL_INCLUDE_DIR}) message("linking DrawMatrix with GLUT and OpenGL")
include_directories(${GLUT_INCLUDE_DIR}) include_directories(${OPENGL_INCLUDE_DIR})
add_executable(DrawMatrix ${SOURCE}) include_directories(${GLUT_INCLUDE_DIR})
target_link_libraries(DrawMatrix inmost ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES}) add_executable(DrawMatrix ${SOURCE})
if(USE_MPI) target_link_libraries(DrawMatrix inmost ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES})
target_link_libraries(DrawMatrix ${MPI_LIBRARIES}) if(USE_MPI)
if(MPI_LINK_FLAGS) target_link_libraries(DrawMatrix ${MPI_LIBRARIES})
set_target_properties(DrawMatrix PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") if(MPI_LINK_FLAGS)
endif() set_target_properties(DrawMatrix PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif(USE_MPI) endif()
if(USE_SOLVER_ANI3D) endif(USE_MPI)
target_link_libraries(DrawMatrix ani3d) if(USE_SOLVER)
endif() if(USE_SOLVER_ANI3D)
if(USE_SOLVER_PETSC) message("linking DrawMatrix with ANI3d")
add_definitions(${PETSC_DEFINITIONS}) target_link_libraries(DrawMatrix ani3d)
target_link_libraries(DrawMatrix ${PETSC_LIBRARIES}) endif()
endif() if(USE_SOLVER_PETSC)
else(GLUT_FOUND) message("linking DrawMatrix with PETSc")
message("GLUT not found") add_definitions(${PETSC_DEFINITIONS})
endif(GLUT_FOUND) target_link_libraries(DrawMatrix ${PETSC_LIBRARIES})
endif()
if(USE_SOLVER_TRILINOS)
message("linking DrawMatrix with Trilinos")
link_directories(${Trilinos_LIBRARY_DIRS})
link_directories(${Trilinos_TPL_LIBRARY_DIRS})
target_link_libraries(DrawMatrix ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
endif()
else(GLUT_FOUND)
message("GLUT not found, not building DrawMatrix")
endif(GLUT_FOUND)
else(OPENGL_FOUND) else(OPENGL_FOUND)
message("OpenGL not found") message("OpenGL not found, not building DrawMatrix")
endif(OPENGL_FOUND) endif(OPENGL_FOUND)
//g++ main.cpp rotate.cpp -L/usr/X11R6/lib -lX11 -lXi -lXmu -lGL -lglut -lGLU ../../inmost.a -O5 //g++ main.cpp rotate.cpp -L/usr/X11R6/lib -lX11 -lXi -lXmu -lGL -lglut -lGLU ../../inmost.a -O5
// press space - explode mesh to see connection // press space - explode mesh to see connection
#define _CRT_SECURE_NO_WARNINGS
#include "../../inmost.h" #include "../../inmost.h"
#include "my_glut.h" #include "my_glut.h"
#include <iostream> #include <iostream>
...@@ -140,7 +141,7 @@ class Reorder_ARMS ...@@ -140,7 +141,7 @@ class Reorder_ARMS
for (INMOST_DATA_ENUM_TYPE k = 0; k < neighbours[block_entries[j]].size(); ++k) for (INMOST_DATA_ENUM_TYPE k = 0; k < neighbours[block_entries[j]].size(); ++k)
A[neighbours[block_entries[j]][k]].SetMarker(); A[neighbours[block_entries[j]][k]].SetMarker();
} }
sblock += block_entries.size(); sblock += static_cast<int>(block_entries.size());
block_ends.push_back(sblock); block_ends.push_back(sblock);
nblocks++; nblocks++;
stack.clear(); stack.clear();
...@@ -168,7 +169,7 @@ class Reorder_ARMS ...@@ -168,7 +169,7 @@ class Reorder_ARMS
} }
INMOST_DATA_ENUM_TYPE GetOrder(INMOST_DATA_ENUM_TYPE ind) INMOST_DATA_ENUM_TYPE GetOrder(INMOST_DATA_ENUM_TYPE ind)
{ {
return neighbours[ind].size(); return static_cast<Storage::enumerator>(neighbours[ind].size());
INMOST_DATA_ENUM_TYPE ret = 0, k; INMOST_DATA_ENUM_TYPE ret = 0, k;
for (k = 0; k < neighbours[ind].size(); k++) for (k = 0; k < neighbours[ind].size(); k++)
if (!A[neighbours[ind][k]].GetMarker()) if (!A[neighbours[ind][k]].GetMarker())
...@@ -306,11 +307,11 @@ public: ...@@ -306,11 +307,11 @@ public:
if (qi == qj && qi != block_ends.end()) if (qi == qj && qi != block_ends.end())
{ {
p = qi - block_ends.begin(); p = qi - block_ends.begin();
if (p < nblocks) if (p < static_cast<unsigned>(nblocks))
nnz[p] += 1.0; nnz[p] += 1.0;
} }
} }
for (INMOST_DATA_ENUM_TYPE k = 0; k < nblocks; ++k) for (INMOST_DATA_ENUM_TYPE k = 0; k < static_cast<unsigned>(nblocks); ++k)
{ {
nnz[k] /= static_cast<INMOST_DATA_REAL_TYPE>((block_ends[k + 1] - block_ends[k])*(block_ends[k + 1] - block_ends[k])); nnz[k] /= static_cast<INMOST_DATA_REAL_TYPE>((block_ends[k + 1] - block_ends[k])*(block_ends[k + 1] - block_ends[k]));
nnzall += nnz[k]; nnzall += nnz[k];
......
...@@ -5,22 +5,32 @@ add_executable(FVDiscr ${SOURCE}) ...@@ -5,22 +5,32 @@ add_executable(FVDiscr ${SOURCE})
target_link_libraries(FVDiscr inmost) target_link_libraries(FVDiscr inmost)
if(USE_SOLVER_ANI) if(USE_SOLVER)
message("linking FVDiscr with ani3d") if(USE_SOLVER_ANI)
target_link_libraries(FVDiscr ani3d) message("linking FVDiscr with ani3d")
target_link_libraries(FVDiscr ani3d)
endif()
if(USE_SOLVER_PETSC)
message("linking FVDiscr with PETSc")
link_directories(${PETSC_LIBRARY_DIRS})
target_link_libraries(FVDiscr ${PETSC_LIBRARIES})
endif()
if(USE_SOLVER_TRILINOS)
message("linking FVDiscr with Trilinos")
link_directories(${Trilinos_LIBRARY_DIRS})
link_directories(${Trilinos_TPL_LIBRARY_DIRS})
target_link_libraries(FVDiscr ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
endif() endif()
if(USE_SOLVER_PETSC)
message("linking FVDiscr with PETSc")
link_directories(${PETSC_LIBRARY_DIRS})
target_link_libraries(FVDiscr ${PETSC_LIBRARIES})
endif()
if(USE_PARTITIONER) if(USE_PARTITIONER)
if(USE_PARTITIONER_ZOLTAN) if(USE_PARTITIONER_ZOLTAN)
message("linking FVDiscr with Zoltan")
target_link_libraries(FVDiscr ${ZOLTAN_LIBRARIES}) target_link_libraries(FVDiscr ${ZOLTAN_LIBRARIES})
endif() endif()
if(USE_PARTITIONER_PARMETIS) if(USE_PARTITIONER_PARMETIS)
message("linking FVDiscr with ParMETIS")
target_link_libraries(FVDiscr ${PARMETIS_LIBRARIES}) target_link_libraries(FVDiscr ${PARMETIS_LIBRARIES})
endif() endif()
endif() endif()
...@@ -32,3 +42,5 @@ if(USE_MPI) ...@@ -32,3 +42,5 @@ if(USE_MPI)
set_target_properties(FVDiscr PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") set_target_properties(FVDiscr PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif() endif()
endif(USE_MPI) endif(USE_MPI)
...@@ -13,12 +13,20 @@ if(USE_MPI) ...@@ -13,12 +13,20 @@ if(USE_MPI)
endif() endif()
endif(USE_MPI) endif(USE_MPI)
if(USE_SOLVER_ANI) if(USE_SOLVER)
message("linking MatSolve with ani3d") if(USE_SOLVER_ANI)
target_link_libraries(MatSolve ani3d) message("linking MatSolve with ani3d")
endif() target_link_libraries(MatSolve ani3d)
if(USE_SOLVER_PETSC) endif()
message("linking MatSolve with PETSc") if(USE_SOLVER_PETSC)
link_directories(${PETSC_LIBRARY_DIRS}) message("linking MatSolve with PETSc")
target_link_libraries(MatSolve ${PETSC_LIBRARIES}) link_directories(${PETSC_LIBRARY_DIRS})
target_link_libraries(MatSolve ${PETSC_LIBRARIES})
endif()
if(USE_SOLVER_TRILINOS)
message("linking MatSolve with Trilinos")
link_directories(${Trilinos_LIBRARY_DIRS})
link_directories(${Trilinos_TPL_LIBRARY_DIRS})
target_link_libraries(MatSolve ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
endif() endif()
...@@ -7,27 +7,21 @@ find_package(OpenGL) ...@@ -7,27 +7,21 @@ find_package(OpenGL)
find_package(GLUT) find_package(GLUT)
if(OPENGL_FOUND) if(OPENGL_FOUND)
if(GLUT_FOUND) if(GLUT_FOUND)
include_directories(${OPENGL_INCLUDE_DIR}) message("linking OldDrawGrid with GLUT and OpenGL")
include_directories(${GLUT_INCLUDE_DIR}) include_directories(${OPENGL_INCLUDE_DIR})
add_executable(OldDrawGrid ${SOURCE}) include_directories(${GLUT_INCLUDE_DIR})
target_link_libraries(OldDrawGrid inmost ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES}) add_executable(OldDrawGrid ${SOURCE})
if(USE_MPI) target_link_libraries(OldDrawGrid inmost ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES})
target_link_libraries(OldDrawGrid ${MPI_LIBRARIES}) if(USE_MPI)
if(MPI_LINK_FLAGS) target_link_libraries(OldDrawGrid ${MPI_LIBRARIES})
set_target_properties(OldDrawGrid PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") if(MPI_LINK_FLAGS)
endif() set_target_properties(OldDrawGrid PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif(USE_MPI) endif()
if(USE_SOLVER_ANI3D) endif(USE_MPI)
target_link_libraries(OldDrawGrid ani3d) else(GLUT_FOUND)
endif() message("GLUT not found, not building OldDrawGrid")
if(USE_SOLVER_PETSC) endif(GLUT_FOUND)
add_definitions(${PETSC_DEFINITIONS})
target_link_libraries(OldDrawGrid ${PETSC_LIBRARIES})
endif()
else(GLUT_FOUND)
message("GLUT not found")
endif(GLUT_FOUND)
else(OPENGL_FOUND) else(OPENGL_FOUND)
message("OpenGL not found") message("OpenGL not found, not building OldDrawGrid")
endif(OPENGL_FOUND) endif(OPENGL_FOUND)
...@@ -14,11 +14,19 @@ if(USE_SOLVER_PETSC) ...@@ -14,11 +14,19 @@ if(USE_SOLVER_PETSC)
link_directories(${PETSC_LIBRARY_DIRS}) link_directories(${PETSC_LIBRARY_DIRS})
target_link_libraries(Solver ${PETSC_LIBRARIES}) target_link_libraries(Solver ${PETSC_LIBRARIES})
endif() endif()
if(USE_SOLVER_TRILINOS)
message("linking Solver with Trilinos")
link_directories(${Trilinos_LIBRARY_DIRS})
link_directories(${Trilinos_TPL_LIBRARY_DIRS})
target_link_libraries(Solver ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_PARTITIONER) if(USE_PARTITIONER)
if(USE_PARTITIONER_ZOLTAN) if(USE_PARTITIONER_ZOLTAN)
message("linking Solver with Zoltan")
target_link_libraries(Solver ${ZOLTAN_LIBRARIES}) target_link_libraries(Solver ${ZOLTAN_LIBRARIES})
endif() endif()
if(USE_PARTITIONER_PARMETIS) if(USE_PARTITIONER_PARMETIS)
message("linking Solver with ParMETIS")
target_link_libraries(Solver ${PARMETIS_LIBRARIES}) target_link_libraries(Solver ${PARMETIS_LIBRARIES})
endif() endif()
endif() endif()
......
...@@ -212,7 +212,7 @@ int main(int argc,char ** argv) ...@@ -212,7 +212,7 @@ int main(int argc,char ** argv)
if( !repartition && m->HaveTag("CM_INDEX") ) if( !repartition && m->HaveTag("CM_INDEX") )
{ {
id = m->GetTag("CM_INDEX"); id = m->GetTag("CM_INDEX");
m->ExchangeData(id,CELL); m->ExchangeData(id,CELL,0);
} }
else else
{ {
...@@ -263,7 +263,7 @@ int main(int argc,char ** argv) ...@@ -263,7 +263,7 @@ int main(int argc,char ** argv)
Solver::Matrix A; Solver::Matrix A;
Solver::Vector x,b; Solver::Vector x,b;
std::map<GeometricData,ElementType> table; Mesh::GeomParam table;
table[MEASURE] = CELL | FACE; table[MEASURE] = CELL | FACE;
table[CENTROID] = CELL | FACE; table[CENTROID] = CELL | FACE;
...@@ -293,10 +293,10 @@ int main(int argc,char ** argv) ...@@ -293,10 +293,10 @@ int main(int argc,char ** argv)
{ {
//~ std::cout << face->LocalID() << " / " << m->NumberOfFaces() << std::endl; //~ std::cout << face->LocalID() << " / " << m->NumberOfFaces() << std::endl;
Element::Status s1,s2; Element::Status s1,s2;
Cell * r1 = face->BackCell(); Cell r1 = face->BackCell();
Cell * r2 = face->FrontCell(); Cell r2 = face->FrontCell();
if( ((r1 == NULL || (s1 = r1->GetStatus()) == Element::Ghost)?0:1)+ if( ((!r1.isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1)+
((r2 == NULL || (s2 = r2->GetStatus()) == Element::Ghost)?0:1) == 0) continue; ((!r2.isValid() || (s2 = r2->GetStatus()) == Element::Ghost)?0:1) == 0) continue;
Storage::real f_nrm[3], r1_cnt[3], r2_cnt[3], f_cnt[3], d1[3], d2[3], T1, T2, Coef; Storage::real f_nrm[3], r1_cnt[3], r2_cnt[3], f_cnt[3], d1[3], d2[3], T1, T2, Coef;
Storage::real f_area = face->Area(); Storage::real f_area = face->Area();
Storage::real vol1 = r1->Volume(), vol2; Storage::real vol1 = r1->Volume(), vol2;
...@@ -308,7 +308,7 @@ int main(int argc,char ** argv) ...@@ -308,7 +308,7 @@ int main(int argc,char ** argv)
f_nrm[2] /= f_area; f_nrm[2] /= f_area;
r1->Barycenter(r1_cnt); r1->Barycenter(r1_cnt);
face->Barycenter(f_cnt); face->Barycenter(f_cnt);
if( r2 == NULL ) //boundary condition if( !r2.isValid() ) //boundary condition
{ {
if( r1->Integer(mat) == 2 ) //direchlet boundary for second material if( r1->Integer(mat) == 2 ) //direchlet boundary for second material
{ {
...@@ -433,7 +433,7 @@ int main(int argc,char ** argv) ...@@ -433,7 +433,7 @@ int main(int argc,char ** argv)
if( m->GetProcessorRank() == 0 ) std::cout << "Retrive data: " << Timer()-ttt << std::endl; if( m->GetProcessorRank() == 0 ) std::cout << "Retrive data: " << Timer()-ttt << std::endl;
ttt = Timer(); ttt = Timer();
m->ExchangeData(phi,CELL); m->ExchangeData(phi,CELL,0);
BARRIER BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Exchange phi: " << Timer()-ttt << std::endl; if( m->GetProcessorRank() == 0 ) std::cout << "Exchange phi: " << Timer()-ttt << std::endl;
...@@ -444,19 +444,19 @@ int main(int argc,char ** argv) ...@@ -444,19 +444,19 @@ int main(int argc,char ** argv)
m->SetParallelStrategy(s); m->SetParallelStrategy(s);
ttt = Timer(); ttt = Timer();
m->ExchangeData(tensor_K,CELL); m->ExchangeData(tensor_K,CELL,0);
BARRIER BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Exchange dense fixed data: " << Timer()-ttt << std::endl; if( m->GetProcessorRank() == 0 ) std::cout << "Exchange dense fixed data: " << Timer()-ttt << std::endl;
ttt = Timer(); ttt = Timer();
m->ExchangeData(m->ProcessorsTag(),CELL); m->ExchangeData(m->ProcessorsTag(),CELL,0);
BARRIER BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Exchange dense variable data: " << Timer()-ttt << std::endl; if( m->GetProcessorRank() == 0 ) std::cout << "Exchange dense variable data: " << Timer()-ttt << std::endl;
ttt = Timer(); ttt = Timer();
m->ExchangeData(test,CELL); m->ExchangeData(test,CELL,0);
BARRIER BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Exchange sparse fixed data: " << Timer()-ttt << std::endl; if( m->GetProcessorRank() == 0 ) std::cout << "Exchange sparse fixed data: " << Timer()-ttt << std::endl;
} }
......
//g++ main.cpp -lpetsc -L/usr/X11R6/lib -lX11 -ldmumps -lmumps_common -lmpi_f77 -lscalapack -lpord -lblacs -lparmetis -lmetis -lmpi -lHYPRE -lgfortran -lblas -llapack ../../mspp.a
//mpicxx main.cpp -lpetsc -L/usr/X11R6/lib -lX11 -ldmumps -lmumps_common -lmpi_f77 -lscalapack -lpord -lblacs -lparmetis -lmetis -lmpi -lHYPRE -lgfortran -lblas -llapack ../../mspp.a -lzoltan -lparmetis -lmetis
// run ./a.out grids/rezultMesh.vtk MATERIALS -ksp_monitor -ksp_view -mat_type aij -vec_type standard
#include "../../mspp.h"
#include <stdio.h>
using namespace MSPP;
#define BARRIER MPI_Barrier(MPI_COMM_WORLD);
void make_vec(Storage::real v1[3], Storage::real v2[3], Storage::real out[3])
{
out[0] = v1[0] - v2[0];
out[1] = v1[1] - v2[1];
out[2] = v1[2] - v2[2];
}
Storage::real dot_prod(Storage::real v1[3], Storage::real v2[3])
{
return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
}
void tensor_prod(Storage::real K[9], Storage::real v[3], Storage::real out[3])
{
out[0] = v[0]*K[0] + v[1]*K[1] + v[2]*K[2];
out[1] = v[0]*K[3] + v[1]*K[4] + v[2]*K[5];
out[2] = v[0]*K[6] + v[1]*K[7] + v[2]*K[8];
}
Storage::real transmissibility(Storage::real vec[3], Storage::real K_tensor[9], Storage::real normal_face[3])
{
Storage::real Kn[3];
tensor_prod(K_tensor,normal_face,Kn);
return dot_prod(vec,Kn);
}
Storage::real tensor_K_mat_const(Storage::integer mat, int comp)
{
if( comp == 0 ) //X
{
switch(mat)
{
case 0: return 850;
case 1: return 950;
case 2: return 1050;
case 3: return 450;
case 4: return 650;
case 5: return 550;
default: return 1000;
}
}
else if( comp == 1 ) //Y
{
switch(mat)
{
case 0: return 750;
case 1: return 450;
case 2: return 950;
case 3: return 350;
case 4: return 650;
case 5: return 750;
default: return 1000;
}
}
else if( comp == 2 ) //Z
{
switch(mat)
{
case 0: return 12;
case 1: return 11;
case 2: return 13;
case 3: return 15;
case 4: return 14;
case 5: return 13;
default: return 10;
}
}
else return 0;
}
Storage::real tensor_K_mat_var(Storage::integer mat, int comp)
{
if( comp == 0 ) //X
{
switch(mat)
{
case 0: return 150;
case 1: return 50;
case 2: return 80;
case 3: return 90;
case 4: return 100;
case 5: return 20;
default: return 10;
}
}
else if( comp == 1 ) //Y
{
switch(mat)
{
case 0: return 150;
case 1: return 50;
case 2: return 80;
case 3: return 90;
case 4: return 100;
case 5: return 20;
default: return 10;
}
}
else if( comp == 2 ) //Z
{
switch(mat)
{
case 0: return 2;
case 1: return 5;
case 2: return 3;