МЕТОДЫ КОМПАКТНОГО ХРАНЕНИЯ МАТРИЦЫ ЖЕСТКОСТИ

2 МЕТОДЫ КОМПАКТНОГО ХРАНЕНИЯ МАТРИЦЫ ЖЕСТКОСТИ

Матрица жесткости, получающаяся при применении МКЭ, обладает симметричной структурой, что позволяет в общем случае хранить только верхнюю треугольную часть матрицы. Однако для задач с большим количеством неизвестных это так же приводит к проблеме нехватки памяти. Предлагаемый в данной работе метод, позволяет хранить только ненулевые члены матрицы жесткости. Суть его заключается в следующем.

Первоначально, с целью выявления связей каждого узла с другими, производится анализ структуры дискретизации области на КЭ. Например, для КЭ - сетки, изображенной на рис. 1, соответствующая структура связей будет иметь вид:

№ узла	1	2	3	4	5	6	7
Связи	1, 2, 5, 6, 7	1, 2, 3, 6	2, 3, 4, 6	3, 4, 5, 6, 7	1, 4, 5, 7	1, 2, 3, 4, 6, 7	1, 4, 5, 6, 7

Т
огда, для хранения матрицы жесткости необходимо построчно запоминать информацию о коэффициентах, соответствующих узлам, с которыми связан данный узел. На рис. 2 приведены матрица жесткости и ее компактное представление для сетки изображенной на рис 1 [9].

Текст подпрограммы, реализующий предложенный алгоритм анализа структуры КЭ-разбиения тела, приведен в Приложении 1.

Данный способ компактного хранения матрицы жесткости позволяет легко его использовать совместно с каким-нибудь численным методом. Наиболее удобным для этой цели представляется использование вышеизложенного итерационного метода Ланцоша, так как на каждой итерации требуется только перемножать матрицу коэффициентов СЛАУ и заданный вектор. Следовательно, для использования предложенного метода компактного хранения СЛАУ необходимо построить прямое и обратное преобразование в первоначальную квадратную матрицу.

Пусть – элемент первоначальной квадратной матрицы размерностью , а - ее компактное представление. Тогда для обратного преобразования будут справедливы следующие соотношения:

, (*)

где m – количество степеней свободы (m=1,2,3).

Для прямого преобразования будут справедливы соотношения, обратные к соотношениям (*).

3 ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Для проверки предлагаемого метода компактного хранения матрицы жесткости была решена задача о контактном взаимодействии оболочечной конструкции и ложемента [12] (рис. 4).

Д
анная задача часто возникает на практике при транспортировке или хранении с горизонтальным расположением оси оболочечные конструкции устанавливаются на круговые опоры - ложементы. Взаимодействие подкрепленных оболочечных конструкций и ложементов осуществляется через опорные шпангоуты, протяженность которых вдоль оси оболочки соизмерима с шириной ложементов и много меньше радиуса оболочки и величины зоны контакта.

Данная задача решалась методом конечных элементов при помощи системы FORL [5]. Дискретная модель ложемента (в трехмерной постановке) представлена на Рис. 5.

При построении данной КЭ-модели было использовано 880 узлов и 2016 КЭ в форме тетраэдра. Полный размер матрицы жесткости для такой задачи составляет байт, что приблизительно равно 2,7 Мбайт оперативной памяти. Размер упакованного представления составил около 315 Кбайт.

Данная задача решалась на ЭВМ с процессором Pentium 166 и 32 МБ ОЗУ двумя способами – методом Гаусса и методом Ланцоша. Сопоставление результатов решения приведено в Таблице 1.

Таблица 1.

	Время решения (сек)
Метод Гаусса	280	2.2101	-2.4608	1.3756	-5.2501	1.7406	-2.3489
Метод Ланцоша	150	2.2137	-2.4669	1.3904	-5.2572	1.7433	-2.3883

Из Таблицы 1 легко видеть, что результаты решения СЛАУ методом Гаусса и методом Ланцоша хорошо согласуются между собой, при этом время решения вторым способом почти в два раза меньше, чем в случае использования метода Гаусса.

ВЫВОДЫ.

В данной работе были рассмотрены способы компактного хранения матрицы коэффициентов системы линейных алгебраических уравнений (СЛАУ) и методы ее решения. Разработан алгоритм компактного хранения матрицы жесткости, позволяющий в несколько раз (иногда более чем в десятки раз) сократить объем требуемой памяти для хранения такой матрицы жесткости.

Классические методы хранения, учитывающие симметричную и ленточную структуру матриц жесткости, возникающих при применении метода конечных элементов (МКЭ), как правило, не применимы при решении контактных задач, так как при их решении матрицы жесткости нескольких тел объединяются в одну общую матрицу, ширина ленты которой может стремиться к порядку системы.

Предложенная в работе методика компактного хранения матриц коэффициентов СЛАУ и использования метода Ланцоша позволили на примере решения контактных задач добиться существенной экономии процессорного времени и затрат оперативной памяти.

СПИСОК ССЫЛОК.

Зенкевич О., Морган К. Конечные методы и аппроксимация // М.: Мир, 1980

Зенкевич О., Метод конечных элементов // М.: Мир., 1975

Стрэнг Г., Фикс Дж. Теория метода конечных элементов // М.: Мир, 1977

Бахвалов Н.С.,Жидков Н.П., Кобельков Г.М. Численные методы // М.: наука, 1987

Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления // М.:Наука, 1984

Бахвалов Н.С. Численные методы // М.: Наука, 1975

Годунов С.К. Решение систем линейных уравнений // Новосибирск: Наука, 1980

Гоменюк С.И., Толок В.А. Инструментальная система анализа задач механики деформируемого твердого тела // Приднiпровський науковий вiсник – 1997. – №4.

F.G. Gustavson, “Some basic techniques for solving sparse matrix algorithms”, // editer by D.J. Rose and R.A.Willoughby, Plenum Press, New York, 1972

А.Джордж, Дж. Лиу, Численное решение больших разреженных систем уравнений // Москва, Мир, 1984

D.J. Rose, “A graph theoretic study of the numerical solution of sparse positive definite system of linear equations” // New York, Academic Press, 1972

Мосаковский В.И., Гудрамович В.С., Макеев Е.М., Контактные задачи теории оболочек и стержней // М.:”Машиностроение”, 1978

ПРИЛОЖЕНИЕ 1

Исходный текст программы, реализующий анализ структуры КЭ-разбиения объекта.

#include

#include

#include

#include

#include

#include "matrix.h"

#define BASE3D_4 4

#define BASE3D_8 8

#define BASE3D_10 10

const double Eps = 1.0E-10;

DWORD CurrentType = BASE3D_10;

void PrintHeader(void)

{

printf("Command format: AConvert - [/Options]\n");

printf("Switch: -t10 - Tetraedr(10)\n");

printf(" -c8 - Cube(8)\n");

printf(" -s4 - Shell(4)\n");

printf(" -s8 - Shell(8)\n\n");

printf("Optins: /8 - convert Tetraedr(10)->8*Tetraedr(4)\n");

printf(" /6 - convert Cube(8)->6*Tetraedr(4)\n");

}

bool Output(char* fname,Vector& x,Vector& y,Vector& z, Matrix& tr, DWORD n,

DWORD NumNewPoints,DWORD ntr,Matrix& Bounds,DWORD CountBn)

{

char* Label = "NTRout";

int type = CurrentType,

type1 = (type==BASE3D_4 || type==BASE3D_10) ? 3 : 4;

DWORD NewSize,

i,

j;

ofstream Out;

if (type==BASE3D_4) type1 = 3;

else if (type==BASE3D_8) type1 = 4;

else if (type==BASE3D_10) type1 = 6;

Out.open(fname,ios::out | ios:: binary);

if (Out.bad()) return true;

Out.write((const char*)Label,6 * sizeof(char));

if (Out.fail()) return true;

Out.write((const char*)&type,sizeof(int));

if (Out.fail()) return true;

Out.write((const char*)&CountBn,sizeof(DWORD));

if (Out.fail())

{

Out.close();

return true;

}

Out.write((const char*)&(NewSize = n + NumNewPoints),sizeof(DWORD));

if (Out.fail()) return true;

Out.write((const char*)&(NumNewPoints),sizeof(DWORD));

if (Out.fail())

{

Out.close();

return true;

}

for (DWORD i = 0; i < n; i++)

{

Out.write((const char*)&x[i],sizeof(double));

Out.write((const char*)&y[i],sizeof(double));

Out.write((const char*)&z[i],sizeof(double));

if (Out.fail())

{

Out.close();

return true;

}

}

for (i = 0; i < NumNewPoints; i++)

{

Out.write((const char*)&x[n + i],sizeof(double));

Out.write((const char*)&y[n + i],sizeof(double));

if (Out.fail())

{

Out.close();

return true;

}

}

Out.write((const char*)&(ntr),sizeof(DWORD));

if (Out.fail())

{

Out.close();

return true;

}

for (i = 0; i < ntr; i++)

for (j = 0; j < (DWORD)type; j++)

{

DWORD out = tr[i][j];

Out.write((const char*)&out,sizeof(DWORD));

if (Out.fail())

{

Out.close();

return true;

}

}

for (i = 0; i < CountBn; i++)

for (j = 0; j < (DWORD)type1; j++)

{

DWORD out = Bounds[i][j];

Out.write((const char*)&out,sizeof(DWORD));

if (Out.fail())

{

Out.close();

return true;

}

}

{

//*********************

// Create Links

printf("Create links...\r");

Vector Link(n),

Size(n);

Matrix Links(n,n);

DWORD Count;

int type = CurrentType;

for (DWORD i = 0; i < n; i++)

{

for (DWORD j = 0; j < ntr; j++)

for (DWORD k = 0; k < (DWORD)type; k++)

if (tr[j][k] == i)

for (DWORD m = 0; m < (DWORD)type; m++) Link[tr[j][m]] = 1;

Count = 0;

for (DWORD m = 0; m < n; m++)

if (Link[m]) Count++;

Size[i] = Count;

Count = 0;

for (DWORD m = 0; m < n; m++)

if (Link[m])

Links[i][Count++] = m;

//Set zero

Link.ReSize(n);

}

// Output

//*********************

for (DWORD i = 0; i < n; i++)

{

DWORD Sz = Size[i];

Out.write((const char*)&Sz,sizeof(DWORD));

for (DWORD j = 0; j < Sz; j++)

Out.write((const char*)&(Links[i][j]),sizeof(DWORD));

}

//*********************

}

printf(" \r");

printf("Points: %ld\n",n);

printf("FE: %ld\n",ntr);

Out.close();

return false;

}

bool Test(DWORD* a,DWORD* b)

{

bool result;

int NumPoints = 3;

if (CurrentType == BASE3D_8) NumPoints = 4;

else if (CurrentType == BASE3D_10) NumPoints = 6;

for (int i = 0; i < NumPoints; i++)

{

result = false;

for (int j = 0; j < NumPoints; j++)

if (b[j] == a[i])

{

result = true;

break;

}

if (result == false) return false;

}

return true;

}

void Convert(Vector& X,Vector& Y,Vector& Z, Matrix& FE,DWORD NumTr,Matrix& Bounds,DWORD& BnCount)

{

int cData8[6][5] = {{0,4,5,1,7},

{6,2,3,7,0},

{4,6,7,5,0},

{2,0,1,3,5},

{1,5,7,3,4},

{6,4,0,2,1}},

cData4[4][4] = {{0,1,2,3},

{1,3,2,0},

{3,0,2,1},

{0,3,1,2}},

cData10[4][7] = {{0,1,2,4,5,6,3},

{0,1,3,4,8,7,2},

{1,3,2,8,9,5,0},

{0,2,3,6,9,7,1}},

cData[6][7],

Data[6],

l,

Num1,

Num2,

m;

DWORD i,

j,

p[6],

pp[6],

Index;

Matrix BoundList(4 * NumTr,6);

double cx,

cy,

cz,

x1,

y1,

z1,

x2,

y2,

z2,

x3,

y3,

z3;

Bounds.ReSize(4 * NumTr,6);

switch (CurrentType)

{

case BASE3D_4:

Num1 = 4;

Num2 = 3;

for (l = 0; l < Num1; l++)

for (m = 0; m < Num2+1; m++)

cData[l][m] = cData4[l][m];

break;

case BASE3D_8:

Num1 = 6;

Num2 = 4;

for (l = 0; l < Num1; l++)

for (m = 0; m < Num2 + 1; m++)

cData[l][m] = cData8[l][m];

break;

case BASE3D_10:

Num1 = 4;

Num2 = 6;

for (l = 0; l < Num1; l++)

for (m = 0; m < Num2+1; m++)

cData[l][m] = cData10[l][m];

}

printf("Create bounds...\r");

for (i = 0; i < NumTr - 1; i++)

for (int j = 0; j < Num1; j++)

if (!BoundList[i][j])

{

for (l = 0; l < Num2; l++)

p[l] = FE[i][cData[j][l]];

for (DWORD k = i + 1; k < NumTr; k++)

for (int m = 0; m < Num1; m++)

if (!BoundList[k][m])

{

for (int l = 0; l < Num2; l++)

pp[l] = FE[k][cData[m][l]];

if (Test(p,pp))

BoundList[i][j] = BoundList[k][m] = 1;

}

}

for (i = 0; i < NumTr; i++)

for (j = 0; j < (DWORD)Num1; j++)

if (BoundList[i][j] == 0)

{

if (CurrentType == BASE3D_4)

{

cx = X[FE[i][cData[j][3]]];

cy = Y[FE[i][cData[j][3]]];

cz = Z[FE[i][cData[j][3]]];

}

else

if (CurrentType == BASE3D_10)

{

cx = X[FE[i][cData[j][6]]];

cy = Y[FE[i][cData[j][6]]];

cz = Z[FE[i][cData[j][6]]];

}

else

{

cx = X[FE[i][cData[j][4]]];

cy = Y[FE[i][cData[j][4]]];

cz = Z[FE[i][cData[j][4]]];

}

x1 = X[FE[i][cData[j][0]]];

y1 = Y[FE[i][cData[j][0]]];

z1 = Z[FE[i][cData[j][0]]];

x2 = X[FE[i][cData[j][1]]];

y2 = Y[FE[i][cData[j][1]]];

z2 = Z[FE[i][cData[j][1]]];

x3 = X[FE[i][cData[j][2]]];

y3 = Y[FE[i][cData[j][2]]];

z3 = Z[FE[i][cData[j][2]]];

for (l = 0; l < Num2; l++)

Data[l] = cData[j][l];

if ( ((cx-x1)*(y2-y1)*(z3-z1) + (cy-y1)*(z2-z1)*(x3-x1) + (y3-y1)*(cz-z1)*(x2-x1) -

(x3-x1)*(y2-y1)*(cz-z1) - (y3-y1)*(z2-z1)*(cx-x1) - (cy-y1)*(z3-z1)*(x2-x1)) > 0)

{

if (CurrentType == BASE3D_4)

{

Data[0] = cData[j][0];

Data[1] = cData[j][2];

Data[2] = cData[j][1];

}

else

if (CurrentType == BASE3D_10)

{

Data[0] = cData[j][0];

Data[1] = cData[j][2];

Data[2] = cData[j][1];

Data[3] = cData[j][5];

Data[5] = cData[j][3];

}

else

{

Data[0] = cData[j][0];

Data[1] = cData[j][3];

Data[2] = cData[j][2];

Data[3] = cData[j][1];

}

}

for (l = 0; l < Num2; l++)

Bounds[BnCount][l] = FE[i][Data[l]];

BnCount++;

}

}

void main(int argc,char** argv)

{

char *input1,

*input2,

*input3,

*op = "",

*sw;

bool CreateFile(char*,char*,char*,char*);

printf("ANSYS->FORL file convertor. ZSU(c) 1998.\n\n");

if (argc < 5 || argc > 6)

{

PrintHeader();

return;

}

sw = argv[1];

input1 = argv[2];

input2 = argv[3];

input3 = argv[4];

if (!strcmp(sw,"-t10"))

CurrentType = BASE3D_10;

else

if (!strcmp(sw,"-c8"))

CurrentType = BASE3D_8;

else

{

printf("Unknown switch %s\n\n",sw);

PrintHeader();

return;

}

if (argc == 6)

{

op = argv[5];

if (strcmp(op,"/8") && strcmp(op,"/6"))

{

printf("Unknown options %s\n\n",op);

PrintHeader();

return;

}

}

if (CreateFile(input1,input2,input3,op))

printf("OK\n");

}

bool CreateFile(char* fn1,char* fn2,char* fn3,char* Op)

{

FILE *in1,

*in2,

*in3;

Vector X(1000),

Y(1000),

Z(1000);

DWORD NumPoints,

NumFE,

NumBounds = 0,

tmp;

Matrix FE(1000,10),

Bounds;

bool ReadTetraedrData(char*,char*,FILE*,FILE*,Vector&,Vector&,Vector&,

Matrix&,DWORD&,DWORD&),

ReadCubeData(char*,char*,FILE*,FILE*,Vector&,Vector&,Vector&,

Matrix&,DWORD&,DWORD&);

void Convert824(Matrix&,DWORD&),

Convert1024(Matrix&,DWORD&);

if ((in1 = fopen(fn1,"r")) == NULL)

{

printf("Unable open file %s",fn1);

return false;

}

if ((in2 = fopen(fn2,"r")) == NULL)

{

printf("Unable open file %s",fn2);

return false;

}

if (CurrentType == BASE3D_10)

{

if (!ReadTetraedrData(fn1,fn2,in1,in2,X,Y,Z,FE,NumPoints,NumFE)) return false;

if (!strcmp(Op,"/8"))

{

// Create 8*Tetraedr(4)

Convert1024(FE,NumFE);

}

Convert(X,Y,Z,FE,NumFE,Bounds,NumBounds);

return !Output(fn3,X,Y,Z,FE,NumPoints,0,NumFE,Bounds,NumBounds);

}

if (CurrentType == BASE3D_8)

{

if (!ReadCubeData(fn1,fn2,in1,in2,X,Y,Z,FE,NumPoints,NumFE)) return false;

if (!strcmp(Op,"/6"))

{

// Create 6*Tetraedr(4)

Convert824(FE,NumFE);

}

Convert(X,Y,Z,FE,NumFE,Bounds,NumBounds);

return !Output(fn3,X,Y,Z,FE,NumPoints,0,NumFE,Bounds,NumBounds);

}

return false;

}

void Convert824(Matrix& FE,DWORD& NumFE)

{

Matrix nFE(6 * NumFE,4);

DWORD data[][4] = {

{ 0,2,3,6 },

{ 4,5,1,7 },

{ 0,4,1,3 },

{ 6,7,3,4 },

{ 1,3,7,4 },

{ 0,4,6,3 }

},

Current = 0;

for (DWORD i = 0; i < NumFE; i++)

for (DWORD j = 0; j < 6; j++)

{

for (DWORD k = 0; k < 4; k++)

nFE[Current][k] = FE[i][data[j][k]];

Current++;

}

CurrentType = BASE3D_4;

NumFE = Current;

FE = nFE;

}

void Convert1024(Matrix& FE,DWORD& NumFE)

{

Matrix nFE(8 * NumFE,4);

DWORD data[][4] = {

{ 3,7,8,9 },

{ 0,4,7,6 },

{ 2,5,9,6 },

{ 7,9,8,6 },

{ 4,8,5,1 },

{ 4,5,8,6 },

{ 7,8,4,6 },

{ 8,9,5,6 }

},

Current = 0;

for (DWORD i = 0; i < NumFE; i++)

for (DWORD j = 0; j < 8; j++)

{

for (DWORD k = 0; k < 4; k++)

nFE[Current][k] = FE[i][data[j][k]];

Current++;

}

CurrentType = BASE3D_4;

NumFE = Current;

FE = nFE;

}

bool ReadTetraedrData(char* fn1,char* fn2,FILE* in1,FILE* in2,Vector& X,Vector& Y,Vector& Z,

Matrix& FE,DWORD& NumPoints,DWORD& NumFE)

{

double tx,

ty,

tz;

char TextBuffer[1001];

DWORD Current = 0L,

tmp;

while (!feof(in1))

{

if (fgets(TextBuffer,1000,in1) == NULL)

{

if (feof(in1)) break;

printf("Unable read file %s",fn1);

fclose(in1);

fclose(in2);

return false;

}

if (sscanf(TextBuffer,"%ld %lf %lf %lf", &NumPoints,&tx,&ty,&tz) != 4) continue;

X[Current] = tx;

Y[Current] = ty;

Z[Current] = tz;

if (++Current == 999)

{

Vector t1 = X,

t2 = Y,

t3 = Z;

X.ReSize(2 * X.Size());

Y.ReSize(2 * X.Size());

Z.ReSize(2 * X.Size());

for (DWORD i = 0; i < Current; i++)

{

X[i] = t1[i];

Y[i] = t2[i];

Z[i] = t3[i];

}

}

if (Current % 100 == 0)

printf("Line: %ld\r",Current);

}

fclose(in1);

printf(" \r");

NumPoints = Current;

Current = 0L;

while (!feof(in2))

{

if (fgets(TextBuffer,1000,in2) == NULL)

{

if (feof(in2)) break;

printf("Unable read file %s",fn2);

fclose(in2);

return false;

}

if (sscanf(TextBuffer,"%d %d %d %d %d %ld %ld %ld %ld %ld %ld %ld %ld",

&tmp,&tmp,&tmp,&tmp,&tmp,

&FE[Current][0],&FE[Current][1],&FE[Current][2],&FE[Current][3],

&FE[Current][4],&FE[Current][5],&FE[Current][6],&FE[Current][7]) != 13) continue;

if (fgets(TextBuffer,1000,in2) == NULL)

{

printf("Unable read file %s",fn2);

fclose(in2);

return false;

}

if (sscanf(TextBuffer,"%ld %ld",&FE[Current][8],&FE[Current][9]) != 2)

{

printf("Unable read file %s",fn2);

fclose(in2);

return false;

}

{

if (fabs((tx = 0.5*(X[FE[Current][0] - 1] + X[FE[Current][1] - 1])) - X[FE[Current][4] - 1]) > Eps)

X[FE[Current][4] - 1] = tx;

if (fabs((ty = 0.5*(Y[FE[Current][0] - 1] + Y[FE[Current][1] - 1])) - Y[FE[Current][4] - 1]) > Eps)

Y[FE[Current][4] - 1] = ty;

if (fabs((tz = 0.5*(Z[FE[Current][0] - 1] + Z[FE[Current][1] - 1])) - Z[FE[Current][4] - 1]) > Eps)

Z[FE[Current][4] - 1] = tz;

if (fabs((tx = 0.5*(X[FE[Current][2] - 1] + X[FE[Current][1] - 1])) - X[FE[Current][5] - 1]) > Eps)

X[FE[Current][5] - 1] = tx;

if (fabs((ty = 0.5*(Y[FE[Current][2] - 1] + Y[FE[Current][1] - 1])) - Y[FE[Current][5] - 1]) > Eps)

Y[FE[Current][5] - 1] = ty;

if (fabs((tz = 0.5*(Z[FE[Current][2] - 1] + Z[FE[Current][1] - 1])) - Z[FE[Current][5] - 1]) > Eps)

Z[FE[Current][5] - 1] = tz;

if (fabs((tx = 0.5*(X[FE[Current][0] - 1] + X[FE[Current][2] - 1])) - X[FE[Current][6] - 1]) > Eps)

X[FE[Current][6] - 1] = tx;

if (fabs((ty = 0.5*(Y[FE[Current][0] - 1] + Y[FE[Current][2] - 1])) - Y[FE[Current][6] - 1]) > Eps)

Y[FE[Current][6] - 1] = ty;

if (fabs((tz = 0.5*(Z[FE[Current][0] - 1] + Z[FE[Current][2] - 1])) - Z[FE[Current][6] - 1]) > Eps)

Z[FE[Current][6] - 1] = tz;

if (fabs((tx = 0.5*(X[FE[Current][0] - 1] + X[FE[Current][3] - 1])) - X[FE[Current][7] - 1]) > Eps)

X[FE[Current][7] - 1] = tx;

if (fabs((ty = 0.5*(Y[FE[Current][0] - 1] + Y[FE[Current][3] - 1])) - Y[FE[Current][7] - 1]) > Eps)

Y[FE[Current][7] - 1] = ty;

if (fabs((tz = 0.5*(Z[FE[Current][0] - 1] + Z[FE[Current][3] - 1])) - Z[FE[Current][7] - 1]) > Eps)

Z[FE[Current][7] - 1] = tz;

if (fabs((tx = 0.5*(X[FE[Current][3] - 1] + X[FE[Current][1] - 1])) - X[FE[Current][8] - 1]) > Eps)

X[FE[Current][8] - 1] = tx;

if (fabs((ty = 0.5*(Y[FE[Current][3] - 1] + Y[FE[Current][1] - 1])) - Y[FE[Current][8] - 1]) > Eps)

Y[FE[Current][8] - 1] = ty;

if (fabs((tz = 0.5*(Z[FE[Current][3] - 1] + Z[FE[Current][1] - 1])) - Z[FE[Current][8] - 1]) > Eps)

Z[FE[Current][8] - 1] = tz;

if (fabs((tx = 0.5*(X[FE[Current][3] - 1] + X[FE[Current][2] - 1])) - X[FE[Current][9] - 1]) > Eps)

X[FE[Current][9] - 1] = tx;

if (fabs((ty = 0.5*(Y[FE[Current][3] - 1] + Y[FE[Current][2] - 1])) - Y[FE[Current][9] - 1]) > Eps)

Y[FE[Current][9] - 1] = ty;

if (fabs((tz = 0.5*(Z[FE[Current][3] - 1] + Z[FE[Current][2] - 1])) - Z[FE[Current][9] - 1]) > Eps)

Z[FE[Current][9] - 1] = tz;

}

if (++Current == 999)

{

Matrix t = FE;

FE.ReSize(2 * FE.Size1(),10);

for (DWORD i = 0; i < Current; i++)

for (DWORD j = 0; j < 10; j++)

FE[i][j] = t[i][j];

}

if (Current % 100 == 0)

printf("Line: %ld\r",Current);

}

NumFE = Current;

for (DWORD i = 0; i < NumFE; i++)

for (DWORD j = 0; j < 10; j++)

FE[i][j]--;

printf(" \r");

return true;

}

bool ReadCubeData(char* fn1,char*fn2,FILE* in1,FILE* in2,Vector& X,Vector& Y,Vector& Z,

Matrix& FE,DWORD& NumPoints,DWORD& NumFE)

{

double tx,

ty,

tz;

char TextBuffer[1001];

DWORD Current = 0L,

tmp;

while (!feof(in1))

{

if (fgets(TextBuffer,1000,in1) == NULL)

{

if (feof(in1)) break;

printf("Unable read file %s",fn1);

fclose(in1);

fclose(in2);

return false;

}

if (sscanf(TextBuffer,"%ld %lf %lf %lf", &NumPoints,&tx,&ty,&tz) != 4) continue;

X[Current] = tx;

Y[Current] = ty;

Z[Current] = tz;

if (++Current == 999)

{

Vector t1 = X,

t2 = Y,

t3 = Z;

X.ReSize(2 * X.Size());

Y.ReSize(2 * X.Size());

Z.ReSize(2 * X.Size());

for (DWORD i = 0; i < Current; i++)

{

X[i] = t1[i];

Y[i] = t2[i];

Z[i] = t3[i];

}

}

if (Current % 100 == 0)

printf("Line: %ld\r",Current);

}

fclose(in1);

printf(" \r");

NumPoints = Current;

Current = 0L;

while (!feof(in2))

{

if (fgets(TextBuffer,1000,in2) == NULL)

{

if (feof(in2)) break;

printf("Unable read file %s",fn2);

fclose(in2);

return false;

}

if (sscanf(TextBuffer,"%d %d %d %d %d %ld %ld %ld %ld %ld %ld %ld %ld",

&tmp,&tmp,&tmp,&tmp,&tmp,

&FE[Current][0],&FE[Current][1],&FE[Current][3],&FE[Current][2],

&FE[Current][4],&FE[Current][5],&FE[Current][7],&FE[Current][6]) != 13) continue;

if (++Current == 999)

{

Matrix t = FE;

FE.ReSize(2 * FE.Size1(),10);

for (DWORD i = 0; i < Current; i++)

for (DWORD j = 0; j < 10; j++)

FE[i][j] = t[i][j];

}

if (Current % 100 == 0)

printf("Line: %ld\r",Current);}

NumFE = Current;

for (DWORD i = 0; i < NumFE; i++)

for (DWORD j = 0; j < 10; j++)

FE[i][j]--;

printf(" \r");

return true;}

ПРИЛОЖЕНИЕ 2.

Исходный текст программы, реализующей алгоритм компактного хранения и решения СЛАУ высокого порядка.

#include "matrix.h"

class RVector

{

private:

Vector Buffer;

public:

RVector(void) {}

~RVector() {}

RVector(DWORD Size) { Buffer.ReSize(Size); }

RVector(RVector& right) { Buffer = right.Buffer; }

RVector(Vector& right) { Buffer = right; }

DWORD Size(void) { return Buffer.Size(); }

void ReSize(DWORD Size) { Buffer.ReSize(Size); }

double& operator [] (DWORD i) { return Buffer[i]; }

RVector& operator = (RVector& right) { Buffer = right.Buffer; return *this; }

RVector& operator = (Vector& right) { Buffer = right; return *this; }

void Sub(RVector&);

void Sub(RVector&,double);

void Mul(double);

void Add(RVector&);

friend double Norm(RVector&,RVector&);

};

class TSMatrix

{

private:

Vector Right;

Vector* Array;

Vector* Links;

uint Dim;

DWORD Size;

public:

TSMatrix(void) { Size = 0; Dim = 0; Array = NULL; Links = NULL; }

TSMatrix(Vector*,DWORD,uint);

~TSMatrix(void) { if (Array) delete [] Array; }

Vector& GetRight(void) { return Right; }

DWORD GetSize(void) { return Size; }

uint GetDim(void) { return Dim; }

Vector& GetVector(DWORD i) { return Array[i]; }

Vector* GetLinks(void) { return Links; }

void SetLinks(Vector* l) { Links = l; }

void Add(Matrix&,Vector&);

void Add(DWORD I, DWORD L, DWORD J, DWORD K, double v)

{

DWORD Row = I,

Col = L * Links[I].Size() * Dim + Find(I,J) * Dim + K;

Array[Row][Col] += v;

}

void Add(DWORD I, double v)

{

Right[I] += v;

}

DWORD Find(DWORD,DWORD);

void Restore(Matrix&);

void Set(DWORD,DWORD,double,bool);

void Set(DWORD Index1,DWORD Index2,double value)

{

DWORD I = Index1 / Dim,

L = Index1 % Dim,

J = Index2 / Dim,

K = Index2 % Dim,

Pos = Find(I,J),

Row = I,

Col;

if (Pos == DWORD(-1)) return;

Col = L * Links[I].Size() * Dim + Find(I,J) * Dim + K;

Array[Row][Col] = value;

}

bool Get(DWORD Index1,DWORD Index2,double& value)

{

DWORD I = Index1 / Dim,

L = Index1 % Dim,

J = Index2 / Dim,

K = Index2 % Dim,

Pos = Find(I,J),

Row = I,

Col;

value = 0;

if (Pos == DWORD(-1)) return false;

Col = L * Links[I].Size() * Dim + Find(I,J) * Dim + K;

value = Array[Row][Col];

return true;

}

void Mul(RVector&,RVector&);

double Mul(DWORD,RVector&);

void write(ofstream&);

void read(ifstream&);

};

class RMatrix

{

private:

Vector Buffer;

DWORD size;

public:

RMatrix(DWORD sz) { size = sz; Buffer.ReSize(size*(size + 1)*0.5); }

~RMatrix() {}

DWORD Size(void) { return size; }

double& Get(DWORD i,DWORD j) { return Buffer[(2*size + 1 - i)*0.5*i + j - i]; }

};

//************************

#include "smatrix.h"

double Norm(RVector& Left,RVector& Right)

{

double Ret = 0;

for (DWORD i = 0; i < Left.Size(); i++)

Ret += Left[i] * Right[i];

return Ret;

}

void RVector::Sub(RVector& Right)

{

for (DWORD i = 0; i < Size(); i++)

(*this)[i] -= Right[i];

}

void RVector::Add(RVector& Right)

{

for (DWORD i = 0; i < Size(); i++)

(*this)[i] += Right[i];

}

void RVector::Mul(double koff)

{

for (DWORD i = 0; i < Size(); i++)

(*this)[i] *= koff;

}

void RVector::Sub(RVector& Right,double koff)

{

for (DWORD i = 0; i < Size(); i++)

(*this)[i] -= Right[i]*koff;

}

TSMatrix::TSMatrix(Vector* links, DWORD size, uint dim)

{

Dim = dim;

Links = links;

Size = size;

Right.ReSize(Dim * Size);

Array = new Vector[Size];

for (DWORD i = 0; i < Size; i++)

Array[i].ReSize(Links[i].Size() * Dim * Dim);

}

void TSMatrix::Add(Matrix& FEMatr,Vector& FE)

{

double Res;

DWORD RRow;

for (DWORD i = 0L; i < FE.Size(); i++)

for (DWORD l = 0L; l < Dim; l++)

for (DWORD j = 0L; j < FE.Size(); j++)

for (DWORD k = 0L; k < Dim; k++)

{

Res = FEMatr[i * Dim + l][j * Dim + k];

if (Res) Add(FE[i],l,FE[j],k,Res);

}

for (DWORD i = 0L; i < FE.Size(); i++)

for (DWORD l = 0L; l < Dim; l++)

{

RRow = FE[UINT(i % (FE.Size()))] * Dim + l;

Res = FEMatr[i * Dim + l][FEMatr.Size1()];

if (Res) Add(RRow,Res);

}

}

DWORD TSMatrix::Find(DWORD I,DWORD J)

{

DWORD i;

for (i = 0; i < Links[I].Size(); i++)

if (Links[I][i] == J) return i;

return DWORD(-1);

}

void TSMatrix::Restore(Matrix& Matr)

{

DWORD i,

j,

NRow,

NPoint,

NLink,

Pos;

Matr.ReSize(Size * Dim,Size * Dim + 1);

for (i = 0; i < Size; i++)

for (j = 0; j < Array[i].Size(); j++)

{

NRow = j / (Array[i].Size() / Dim); // Number of row

NPoint = (j - NRow * (Array[i].Size() / Dim)) / Dim; // Number of points

NLink = j % Dim; // Number of link

Pos = Links[i][NPoint];

Matr[i * Dim + NRow][Pos * Dim + NLink] = Array[i][j];

}

for (i = 0; i < Right.Size(); i++) Matr[i][Matr.Size1()] = Right[i];

}

void TSMatrix::Set(DWORD Index,DWORD Position,double Value,bool Case)

{

DWORD Row = Index,

Col = Position * Links[Index].Size() * Dim + Find(Index,Index) * Dim + Position,

i;

double koff = Array[Row][Col],

val;

if (!Case)

Right[Dim * Index + Position] = Value;

else

{

Right[Index * Dim + Position] = Value * koff;

for (i = 0L; i < Size * Dim; i++)

if (i != Index * Dim + Position)

{

Set(Index * Dim + Position,i,0);

Set(i,Index * Dim + Position,0);

if (Get(i,Index * Dim + Position,val))

Right[i] -= val * Value;

}

}

}

void TSMatrix::Mul(RVector& Arr,RVector& Res)

{

DWORD i,

j,

NRow,

NPoint,

NLink,

Pos;

Res.ReSize(Arr.Size());

for (i = 0; i < Size; i++)

for (j = 0; j < Array[i].Size(); j++)

{

NRow = j / (Array[i].Size() / Dim);

NPoint = (j - NRow * (Array[i].Size() / Dim)) / Dim;

NLink = j % Dim;

Pos = Links[i][NPoint];

Res[i * Dim + NRow] += Arr[Pos * Dim + NLink] * Array[i][j];

}

}

double TSMatrix::Mul(DWORD Index,RVector& Arr)

{

DWORD j,

I = Index / Dim,

L = Index % Dim,

Start = L * (Array[I].Size() / Dim),

Stop = Start + (Array[I].Size() / Dim),

NRow,

NPoint,

NLink,

Pos;

double Res = 0;

for (j = Start; j < Stop; j++)

{

NRow = j / (Array[I].Size() / Dim);

NPoint = (j - NRow * (Array[I].Size() / Dim)) / Dim;

NLink = j % Dim;

Pos = Links[I][NPoint];

Res += Arr[Pos * Dim + NLink] * Array[I][j];

}

return Res;

}

void TSMatrix::write(ofstream& Out)

{

DWORD ColSize;

Out.write((char*)&(Dim),sizeof(DWORD));

Out.write((char*)&(Size),sizeof(DWORD));

for (DWORD i = 0; i < Size; i++)

{

ColSize = Array[i].Size();

Out.write((char*)&(ColSize),sizeof(DWORD));

for (DWORD j = 0; j < ColSize; j++)

Out.write((char*)&(Array[i][j]),sizeof(double));

}

for (DWORD i = 0; i < Size * Dim; i++)

Out.write((char*)&(Right[i]),sizeof(double));

}

void TSMatrix::read(ifstream& In)

{

DWORD ColSize;

In.read((char*)&(Dim),sizeof(DWORD));

In.read((char*)&(Size),sizeof(DWORD));

if (Array) delete [] Array;

Array = new Vector[Size];

Right.ReSize(Size * Dim);

for (DWORD i = 0; i < Size; i++)

{

In.read((char*)&(ColSize),sizeof(DWORD));

Array[i].ReSize(ColSize);

for (DWORD j = 0; j < ColSize; j++)

In.read((char*)&(Array[i][j]),sizeof(double));

}

for (DWORD i = 0; i < Size * Dim; i++)

In.read((char*)&(Right[i]),sizeof(double));