Commit 2af8a709 authored by Dmitry Bagaev's avatar Dmitry Bagaev
Browse files

Merge branch 'INMOST-DEV/master' into fixfc2

parents 6b1fd49b 86e6d1f4
...@@ -15,4 +15,5 @@ if(USE_AUTODIFF AND USE_SOLVER AND USE_MESH) ...@@ -15,4 +15,5 @@ if(USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
add_subdirectory(ADFVDiscr) add_subdirectory(ADFVDiscr)
add_subdirectory(ADMFD) add_subdirectory(ADMFD)
endif(USE_AUTODIFF AND USE_SOLVER AND USE_MESH) endif(USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
#add_subdirectory(OctreeCutcell) #add_subdirectory(OctreeCutcell)
\ No newline at end of file add_subdirectory(Octree)
\ No newline at end of file
CXX=mpicxx
CC=mpicxx
#links for cluster
PETSC_DIR=/data4t/terekhov/Packages/petsc-3.3-p3
PETSC_ARCH=arch-linux2-c-opt
INCPATH = -I../../ #-I/data4t/terekhov/program/include -I${PETSC_DIR}/${PETSC_ARCH}/include -I${PETSC_DIR}/include
OPTFLAGS= -O0 -g -Wall
#OPTFLAGS= -O3
#MYLIBS=../../mspp.a
MYLIBS= ../../msppd.a
CXXFLAGS=$(OPTFLAGS) $(INCPATH)
LDFLAGS=$(OPTFLAGS)
PETSC = -L${PETSC_DIR}/${PETSC_ARCH}/lib/
#-lpetsc -L/usr/X11R6/lib -lX11 -ldmumps -lmumps_common -lmpi_f77 -lscalapack -lpord -lblacs -lparmetis -lmetis -lflapack -lfblas -lmpi -lHYPRE -lsuperlu_dist_2.4
#PETSC+=-lpetsc -L/usr/X11R6/lib64 -lX11 -lpthread -lmpigc4 -lstdc++ -ldmumps -lmumps_common -lscalapack -lpord -lblacs -lparmetis -lmetis -lmpi -lHYPRE -lflapack -lfblas -L/opt/intel/impi/3.2.2.006/lib64 -L/opt/intel/mpi-rt/3.2.2 -lmpigc4 -lmpigf -lmpi ../../ILU2/lib/libilu-2.3.a -lgfortran
PETSC =-lpetsc -L/usr/X11R6/lib -lX11 -lparmetis -lmetis -lmpi_f77 -lflapack -lfblas -lgfortran
ILU2 = ../../ILU2/lib/libilu-2.3.a -lgfortran
LDLIBS=$(MYLIBS) $(PETSC) $(ILU2)
targets=main
all: $(targets)
main: main.o
clean:
rm -f $(targets) *.o
project(Octree)
set(SOURCE main.cpp
rotate.cpp
rotate.h
octgrid.cpp
octgrid.h
my_glut.h)
add_executable(Octree ${SOURCE})
target_link_libraries(Octree inmost lapack blas)
if(USE_MPI)
message("linking with mpi")
target_link_libraries(Octree ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
message("adding link flags")
set_target_properties(Octree PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
find_package(OpenGL)
find_package(GLUT)
if(OPENGL_FOUND)
if(GLUT_FOUND)
include_directories(${OPENGL_INCLUDE_DIR})
include_directories(${GLUT_INCLUDE_DIR})
add_definitions(-D__GRAPHICS__)
target_link_libraries(Octree ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES})
else(GLUT_FOUND)
message("GLUT not found")
endif(GLUT_FOUND)
else(OPENGL_FOUND)
message("OpenGL not found")
endif(OPENGL_FOUND)
if(USE_PARTITIONER_ZOLTAN)
message("linking Octree with Zoltan")
target_link_libraries(Octree ${ZOLTAN_LIBRARIES})
endif()
if(USE_PARTITIONER_PARMETIS)
message("linking Octree with ParMETIS")
target_link_libraries(Octree ${PARMETIS_LIBRARIES})
endif()
This diff is collapsed.
#if defined(__GRAPHICS__)
#if defined (__APPLE__) || defined(MAXOSX)
#include <GLUT/glut.h>
#endif
#if defined(_WIN32)
//#include <GL/glut.h>
#include <windows.h>
#include <GL/glut.h>
//#include "glut.h"
#pragma comment(lib,"glut32.lib")
#endif
#if defined(__linux__)
#include <GL/glut.h>
#endif
#endif
This diff is collapsed.
#ifndef _OCTGRID_H
#define _OCTGRID_H
#include "inmost.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define DIM 3
using namespace INMOST;
#define TSNH {printf("(%s:%d) this should not happen!\n",__FILE__,__LINE__); exit(-1);}
#define WAITNL {int ret; char c; printf("press Enter\n"); ret = scanf("%c",&c);}
#define VBUFFER 131072
#define CBUFFER 131072
/// Struct conatains tags for one cell
struct Cell_Tags
{
/// Center of cell (3 coordinates)
Tag center;
/// Sides of cell (3 doubles)
Tag side;
/// Number of cell in z dimension (up/down).
Tag floor;
/// Only for child cell. Octree division means 8 children. This tag contains id among other children.
/// Numeration: 1 2 above them 5 6
/// 0 3 4 7
Tag chld_num;
/// Child id for parent
Tag par_chld_nums;
/// Level of partition
Tag level;
/// Cell should be split
Tag to_split;
Tag busy;
Tag leaf;
Tag vertexes;
Tag children;
Tag parent;
Tag is_valid;
/// Number of processor before redistribute
Tag proc;
Tag i;
};
/// Main object, contains all information about mesh
struct grid
{
/// Inmost mesh
Mesh* mesh;
Cell_Tags c_tags;
void (*transformation)(double xyz[3]);
void (*rev_transformation)(double xyz[3]);
int (*cell_should_split)(struct grid * g, Cell cell, int level);
int (*cell_should_unite)(struct grid * g, Cell cell);
};
/// Often after redistribution brother cells splits to different processor
/// This function corrects this division. Puts all children to one processor
void correct_brothers(struct grid* g, int size, int rank, int type);
/// Debug. Print information about all cells
void print_all_about_cells(struct grid * g);
/// Initialize grid. Create tags and other data
void gridInit(struct grid * g, int n[3]);
/// Split cells using "cell_should_split" rule
void gridRefine(struct grid * g);
/// Unite cells using "cell_should_unite" rule
void gridCoarse(struct grid * g);
/// Unite and coarse cells
void gridAMR(struct grid * g, int action);
/// Default method of transformation
void default_transformation(double xyz[3]);
/// Default method of unite
int default_cell_should_unite(struct grid * g, int cell);
/// Default method of split
int default_cell_should_split(struct grid * g, int cell);
#endif
#if defined(__GRAPHICS__)
/*********************************
Implementation of quaternion-based
object rotation
Functions
clickmotion - call when holded mouse moved
click - call when user clicks
motion - call when mouse just moves
quatinit - flush rotation value
rotate - multiply GL matrix
Dependency: rotate.h,
Standard: math.h
Specific: glut.h
**********************************/
#include "my_glut.h"
#include "rotate.h"
#include "math.h"
#include <stdio.h>
struct quaternion
{
double x,y,z,w;
};
struct vector
{
double x,y,z;
};
// Rotation
struct quaternion q;
struct vector drag, onclick;
double rotate_mx,rotate_my;
extern int width, height;
//
void rotate_clickmotion(int nmx, int nmy) // Mouse
{
struct vector n;
double norm,length,t;
rotate_mx = 2.*(nmx/(double)width - 0.5);
rotate_my = 2.*(0.5 - nmy/(double)height);
norm = rotate_mx*rotate_mx + rotate_my*rotate_my;
if( norm > 1.0 )
{
length = sqrt(norm);
drag.x = rotate_mx/length;
drag.y = rotate_my/length;
drag.z = 0.0;
}
else
{
drag.x = rotate_mx;
drag.y = rotate_my;
drag.z = sqrt(1.0-norm);
}
n.x = drag.y*onclick.z - drag.z*onclick.y;
n.y = drag.z*onclick.x - drag.x*onclick.z;
n.z = drag.x*onclick.y - drag.y*onclick.x;
if ( n.x*n.x + n.y*n.y + n.z*n.z > 10e-7 )
{
t = drag.x*onclick.x + drag.y*onclick.y + drag.z*onclick.z;
q.x = + q.x*t + q.y*n.z - q.z*n.y + q.w*n.x;
q.y = - q.x*n.z + q.y*t + q.z*n.x + q.w*n.y;
q.z = + q.x*n.y - q.y*n.x + q.z*t + q.w*n.z;
q.w = - q.x*n.x - q.y*n.y - q.z*n.z + q.w*t;
onclick.x = drag.x;
onclick.y = drag.y;
onclick.z = drag.z;
}
glutPostRedisplay();
}
void rotate_motion(int nmx, int nmy) // Mouse
{
rotate_mx = 2.*(nmx/(double)width - 0.5);
rotate_my = 2.*(0.5 - nmy/(double)height);
}
void rotate_click(int b, int s, int nmx, int nmy) // Mouse
{
(void ) s;
double norm,length;
switch(b)
{
case GLUT_LEFT_BUTTON:
rotate_mx = 2.*(nmx/(double)width - 0.5);
rotate_my = 2.*(0.5 - nmy/(double)height);
norm = rotate_mx*rotate_mx + rotate_my*rotate_my;
if( norm > 1.0 )
{
length = sqrt(norm);
drag.x = rotate_mx/length;
drag.y = rotate_my/length;
drag.z = 0.0;
}
else
{
drag.x = rotate_mx;
drag.y = rotate_my;
drag.z = sqrt(1.0-norm);
}
onclick.x = drag.x;
onclick.y = drag.y;
onclick.z = drag.z;
break;
}
glutPostRedisplay();
}
void quatinit()
{
q.x = 0.0;
q.y = 0.0;
q.z = 0.0;
q.w = 1.0;
}
void rotate()
{
double rot[16];
rot[ 0] = (q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z);
rot[ 1] = 2.*(q.x*q.y - q.w*q.z);
rot[ 2] = 2.*(q.x*q.z + q.w*q.y);
rot[ 3] = 0.0;
rot[ 4] = 2.*(q.x*q.y + q.w*q.z);
rot[ 5] = (q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z);
rot[ 6] = 2.*(q.y*q.z - q.w*q.x);
rot[ 7] = 0.0;
rot[ 8] = 2.*(q.x*q.z - q.w*q.y);
rot[ 9] = 2.*(q.y*q.z + q.w*q.x);
rot[10] = (q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z);
rot[11] = 0.0;
rot[12] = 0.0;
rot[13] = 0.0;
rot[14] = 0.0;
rot[15] = (q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
glMultMatrixd(rot);
}
#endif
#if defined(__GRAPHICS__)
/*********************************
Implementation of quaternion-based
object rotation
Functions
clickmotion - call when holded mouse moved
click - call when user clicks
motion - call when mouse just moves
quatinit - flush rotation value
rotate - multiply GL matrix
Dependency: rotate.h,
Standard: math.h
Specific: glut.h
**********************************/
#ifndef _ROTATE_H
#define _ROTATE_H
void rotate_clickmotion(int nmx, int nmy);
void rotate_motion(int nmx, int nmy);
void rotate_click(int b, int s, int nmx, int nmy);
void quatinit();
void rotate();
#endif
#endif
...@@ -220,8 +220,10 @@ namespace INMOST ...@@ -220,8 +220,10 @@ namespace INMOST
friend class Storage; friend class Storage;
friend class Mesh; friend class Mesh;
}; };
class TagManager //implemented in tag.cpp class TagManager //implemented in tag.cpp
{ {
protected: protected:
TagManager(); TagManager();
...@@ -607,6 +609,123 @@ namespace INMOST ...@@ -607,6 +609,123 @@ namespace INMOST
__INLINE ElementSet getAsSet () const; __INLINE ElementSet getAsSet () const;
friend class Mesh; friend class Mesh;
}; };
//////////////////////////////////////////////////////////////////////
/// Helper classes for class Tag //
//////////////////////////////////////////////////////////////////////
class TagReal : public Tag
{
public:
TagReal() : Tag() {}
TagReal(const TagReal & b) : Tag(b) {}
TagReal(const Tag & b) : Tag(b) {}
TagReal & operator = (TagReal const & b) {Tag::operator =(b); return *this;}
TagReal & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::real & operator [](const Storage & arg) const {return arg.Real(*static_cast<const Tag*>(this));}
};
class TagInteger : public Tag
{
public:
TagInteger() : Tag() {}
TagInteger(const TagInteger & b) : Tag(b) {}
TagInteger(const Tag & b) : Tag(b) {}
TagInteger & operator = (TagInteger const & b) {Tag::operator =(b); return *this;}
TagInteger & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::integer & operator [](const Storage & arg) const {return arg.Integer(*static_cast<const Tag*>(this));}
};
class TagBulk : public Tag
{
public:
TagBulk() : Tag() {}
TagBulk(const TagBulk & b) : Tag(b) {}
TagBulk(const Tag & b) : Tag(b) {}
TagBulk & operator = (TagBulk const & b) {Tag::operator =(b); return *this;}
TagBulk & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::bulk & operator [](const Storage & arg) const {return arg.Bulk(*static_cast<const Tag*>(this));}
};
class TagReference : public Tag
{
public:
TagReference() : Tag() {}
TagReference(const TagReference & b) : Tag(b) {}
TagReference(const Tag & b) : Tag(b) {}
TagReference & operator = (TagReference const & b) {Tag::operator =(b); return *this;}
TagReference & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::reference & operator [](const Storage & arg) const {return arg.Reference(*static_cast<const Tag*>(this));}
};
class TagRealArray : public Tag
{
public:
TagRealArray() : Tag() {}
TagRealArray(const TagRealArray & b) : Tag(b) {}
TagRealArray(const Tag & b) : Tag(b) {}
TagRealArray & operator = (TagRealArray const & b) {Tag::operator =(b); return *this;}
TagRealArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::real_array operator [](const Storage & arg) const {return arg.RealArray(*static_cast<const Tag*>(this));}
};
class TagIntegerArray : public Tag
{
public:
TagIntegerArray() : Tag() {}
TagIntegerArray(const TagIntegerArray & b) : Tag(b) {}
TagIntegerArray(const Tag & b) : Tag(b) {}
TagIntegerArray & operator = (TagIntegerArray const & b) {Tag::operator =(b); return *this;}
TagIntegerArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::integer_array operator [](const Storage & arg) const {return arg.IntegerArray(*static_cast<const Tag*>(this));}
};
class TagBulkArray : public Tag
{
public:
TagBulkArray() : Tag() {}
TagBulkArray(const TagBulkArray & b) : Tag(b) {}
TagBulkArray(const Tag & b) : Tag(b) {}
TagBulkArray & operator = (TagBulkArray const & b) {Tag::operator =(b); return *this;}
TagBulkArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::bulk_array operator [](const Storage & arg) const {return arg.BulkArray(*static_cast<const Tag*>(this));}
};
class TagReferenceArray : public Tag
{
public:
TagReferenceArray() : Tag() {}
TagReferenceArray(const TagReferenceArray & b) : Tag(b) {}
TagReferenceArray(const Tag & b) : Tag(b) {}
TagReferenceArray & operator = (TagReferenceArray const & b) {Tag::operator =(b); return *this;}
TagReferenceArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::reference_array operator [](const Storage & arg) const {return arg.ReferenceArray(*static_cast<const Tag*>(this));}
};
#if defined(USE_AUTODIFF)
class TagVariable : public Tag
{
public:
TagVariable() : Tag() {}
TagVariable(const TagVariable & b) : Tag(b) {}
TagVariable(const Tag & b) : Tag(b) {}
TagVariable & operator = (TagVariable const & b) {Tag::operator =(b); return *this;}
TagVariable & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::var & operator [](const Storage & arg) const {return arg.Variable(*static_cast<const Tag*>(this));}
};
class TagVariableArray : public Tag
{
public:
TagVariableArray() : Tag() {}
TagVariableArray(const TagVariableArray & b) : Tag(b) {}
TagVariableArray(const Tag & b) : Tag(b) {}
TagVariableArray & operator = (TagVariableArray const & b) {Tag::operator =(b); return *this;}
TagVariableArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::var_array operator [](const Storage & arg) const {return arg.VariableArray(*static_cast<const Tag*>(this));}
};
#endif //USE_AUTODIFF
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
/// Inline functions for class Tag // /// Inline functions for class Tag //
......
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
#endif #endif
#include <iomanip> #include <iomanip>
// Matrix with n columns and m rows // Matrix with n rows and m columns
// __m__ // __m__
// | | // | |
// n| | // n| |
...@@ -611,6 +611,19 @@ namespace INMOST ...@@ -611,6 +611,19 @@ namespace INMOST
} }
return ret; return ret;
} }
/// Inverts matrix using Crout-LU decomposition with full pivoting for
/// maximum element. Works well for non-singular matrices, for singular
/// matrices look see Matrix::PseudoInvert.
/// @param print_fail Verbose output for singular matrices.
/// @return A pair of inverse matrix and boolean. If boolean is true,
/// then the matrix was inverted successfully.
/// @see Matrix::PseudoInvert.
/// \todo
/// 1. Introduce Matrix::Solve, that uses the same algorithm with
/// right hand side, a matrix. Then Matrix::Invert is Matrix::Solve
/// with unit matrix.
/// 2. Maximum product transversal + block pivoting instead of pivoting
/// by maxium element.
std::pair<Matrix,bool> Invert(bool print_fail = false) const std::pair<Matrix,bool> Invert(bool print_fail = false) const
{ {
std::pair<Matrix,bool> ret = std::make_pair(Matrix(m,n),true); std::pair<Matrix,bool> ret = std::make_pair(Matrix(m,n),true);
...@@ -793,8 +806,27 @@ namespace INMOST ...@@ -793,8 +806,27 @@ namespace INMOST
for(enumerator i = 0; i < n*m; ++i) ret += space[i]*space[i]; for(enumerator i = 0; i < n*m; ++i) ret += space[i]*space[i];
return sqrt(ret); return sqrt(ret);
} }
/// Convert values in array into matrix. /// Convert values in array into square matrix.
/// Depending on /// Supports the following representation, depending on the size
/// of input array and size of side of final tensors' matrix:
///
/// representation | (array size, tensor size)
///
/// scalar | (1,1), (1,2), (1,3), (1,6)
///
/// diagonal | (2,2), (2,2), (3,3), (6,6)
///
/// symmetric | (3,2), (6,3), (21,6)
///
/// full | (4,2), (9,3), (36,6)
///
/// For full matrix elements in array are enumerated row by row.
/// For symmetric matrix elements in array are enumerated row by
/// row starting from diagonal.
/// @param K Array of elements to be converted into tensor.
/// @param size Size of the input array.
/// @param matsize Size of the final tensor.
/// @return Matrix of the tensor of size matsize by matsize.
static Matrix<Var> FromTensor(Var * K, enumerator size, enumerator matsize = 3) static Matrix<Var> FromTensor(Var * K, enumerator size, enumerator matsize = 3)
{ {
Matrix<Var> Kc(matsize,matsize); Matrix<Var> Kc(matsize,matsize);
...@@ -909,43 +941,8 @@ namespace INMOST ...@@ -909,43 +941,8 @@ namespace INMOST
} }
return Kc; return Kc;
} }
static Matrix<Var> FromElasticTensor(Var * K, enumerator size) ///Retrive vector in matrix form from array.
{ ///
Matrix<Var> Kc(6,6);
switch(size)
{
case 1: //scalar permeability tensor
Kc(0,0) = Kc(1,1) = Kc(2,2) = Kc(3,3) = Kc(4,4) = Kc(5,5) = K[0];