 |
Итерационные методы решения СЛАУ
|
|
|
|
Метод итераций (метод последовательных приближений).
Приближенные методы решения систем линейных уравнений позволяют получать значения корней системы с заданной точностью в виде предела последовательности некоторых векторов. Процесс построения такой последовательности называется итерационным (повторяющимся).
Эффективность применения приближенных методов зависят от выбора начального вектора и быстроты сходимости процесса.
Рассмотрим метод итераций (метод последовательных приближений). Пусть дана система n линейных уравнений с n неизвестными:
Ах=b, (14)
Предполагая, что диагональные элементы aii 
0 (i = 2,..., n), выразим xi через первое уравнение систем x2 - через второе уравнение и т. д. В результате получим систему, эквивалентную системе (14):
Обозначим 
; 
, где i == 1, 2,...,n; j == 1,2,..., n. Тогда система (15) запишется таким образом в матричной форме
Решим систему (16) методом последовательных приближений. За нулевое приближение примем столбец свободных членов. Любое (k+1)-е приближение вычисляют по формуле
Если последовательность приближений x(0),...,x(k) имеет предел 
, то этот предел является решением системы (15), поскольку в силу свойства предела 
, т.е. 
[4,6].
Метод Зейделя.
Метод Зейделя представляет собой модификацию метода последовательных приближений. В методе Зейделя при вычислении (k+1)-го приближения неизвестного xi (i>1) учитываются уже найденные ранее (k+1)-е приближения неизвестных xi-1.
Пусть дана линейная система, приведенная к нормальному виду:
(17)
Выбираем произвольно начальные приближения неизвестных и подставляем в первое уравнение системы (17). Полученное первое приближение подставляем во второе уравнение системы и так далее до последнего уравнения. Аналогично строим вторые, третьи и т.д. итерации.
|
|
|
Таким образом, предполагая, что k-е приближения 
известны, методом Зейделя строим (k+1)-е приближение по следующим формулам:
где k=0,1,...,n
Метод Ланцоша.
Для решения СЛАУ высокого порядка (1), матрица, коэффициентов которой хранится в компактном нижеописанном виде, наиболее удобным итерационным методом является метод Ланцоша [4], схема которого имеет вид:
(18)
где
Преимуществом данного метода является его высокая скорость сходимости к точному решению. Кроме того, доказано, что он обладает свойством «квадратичного окончания», т.е. для положительно определенной матрицы можно гарантировано получить точное решение при количестве итераций 
. Размер требуемой памяти на каждой итерации не изменяется, т.к. не требует преобразование матрицы 
. В качестве критерия остановки данного итерационного процесса обычно используют соотношение
, (19)
где 
- заданная точность. В качестве другого критерия сходимости иногда удобнее использовать среднеквадратичную разность между решениями, полученными на соседних итерациях:
(20)
Среднеквадратичную разность необходимо контролировать при выполнении каждых k наперед заданных итераций.
Отдельно следует рассмотреть проблему выбора начального приближения 
. Доказывается, что при положительно определенной матрице , итерационный процесс (18) всегда сходится при любом выборе начального приближения. При решении контактных задач, когда для уточнения граничных условий в зоне предполагаемого контакта требуется большое количество решений СЛАУ вида (1), в качестве начального приближения для первого расчета используется правая часть системы (1), а для каждого последующего пересчета - решение, полученное на предыдущем. Такая схема позволяет значительно сократить количество итераций, необходимых для достижения заданной точности (19) или (20) [10,11].
|
|
|
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 легко видеть, что результаты решения СЛАУ методом Гаусса и методом Ланцоша хорошо согласуются между собой, при этом время решения вторым способом почти в два раза меньше, чем в случае использования метода Гаусса.
ВЫВОДЫ.
В данной работе были рассмотрены способы компактного хранения матрицы коэффициентов системы линейных алгебраических уравнений (СЛАУ) и методы ее решения. Разработан алгоритм компактного хранения матрицы жесткости, позволяющий в несколько раз (иногда более чем в десятки раз) сократить объем требуемой памяти для хранения такой матрицы жесткости.
Классические методы хранения, учитывающие симметричную и ленточную структуру матриц жесткости, возникающих при применении метода конечных элементов (МКЭ), как правило, не применимы при решении контактных задач, так как при их решении матрицы жесткости нескольких тел объединяются в одну общую матрицу, ширина ленты которой может стремиться к порядку системы.
|
|
|
Предложенная в работе методика компактного хранения матриц коэффициентов СЛАУ и использования метода Ланцоша позволили на примере решения контактных задач добиться существенной экономии процессорного времени и затрат оперативной памяти.
ПРИЛОЖЕНИЕ 1
Исходный текст программы, реализующий анализ структуры КЭ-разбиения объекта.
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <fstream.h>
#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 -<switch> <file1.in> <file2.in> <file3.out> [/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<double>& x,Vector<double>& y,Vector<double>& z, Matrix<DWORD>& tr, DWORD n,
DWORD NumNewPoints,DWORD ntr,Matrix<DWORD>& 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<DWORD> Link(n),
Size(n);
Matrix<DWORD> 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<double>& X,Vector<double>& Y,Vector<double>& Z, Matrix<DWORD>& FE,DWORD NumTr,Matrix<DWORD>& 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<DWORD> 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<double> X(1000),
Y(1000),
Z(1000);
DWORD NumPoints,
NumFE,
NumBounds = 0,
tmp;
Matrix<DWORD> FE(1000,10),
Bounds;
bool ReadTetraedrData(char*,char*,FILE*,FILE*,Vector<double>&,Vector<double>&,Vector<double>&,
Matrix<DWORD>&,DWORD&,DWORD&),
ReadCubeData(char*,char*,FILE*,FILE*,Vector<double>&,Vector<double>&,Vector<double>&,
Matrix<DWORD>&,DWORD&,DWORD&);
void Convert824(Matrix<DWORD>&,DWORD&),
Convert1024(Matrix<DWORD>&,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<DWORD>& FE,DWORD& NumFE)
{
Matrix<DWORD> 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<DWORD>& FE,DWORD& NumFE)
{
Matrix<DWORD> 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<double>& X,Vector<double>& Y,Vector<double>& Z,
Matrix<DWORD>& 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<double> 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<DWORD> 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<double>& X,Vector<double>& Y,Vector<double>& Z,
Matrix<DWORD>& 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<double> 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<DWORD> 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<double> Buffer;
public:
RVector(void) {}
~RVector() {}
RVector(DWORD Size) { Buffer.ReSize(Size); }
RVector(RVector& right) { Buffer = right.Buffer; }
RVector(Vector<double>& 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<double>& 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<double> Right;
Vector<double>* Array;
Vector<DWORD>* Links;
uint Dim;
DWORD Size;
public:
TSMatrix(void) { Size = 0; Dim = 0; Array = NULL; Links = NULL; }
TSMatrix(Vector<DWORD>*,DWORD,uint);
~TSMatrix(void) { if (Array) delete [] Array; }
Vector<double>& GetRight(void) { return Right; }
DWORD GetSize(void) { return Size; }
uint GetDim(void) { return Dim; }
Vector<double>& GetVector(DWORD i) { return Array[i]; }
Vector<DWORD>* GetLinks(void) { return Links; }
void SetLinks(Vector<DWORD>* l) { Links = l; }
void Add(Matrix<double>&,Vector<DWORD>&);
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<double>&);
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<double> 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<DWORD>* links, DWORD size, uint dim)
{
Dim = dim;
Links = links;
Size = size;
Right.ReSize(Dim * Size);
Array = new Vector<double>[Size];
for (DWORD i = 0; i < Size; i++)
Array[i].ReSize(Links[i].Size() * Dim * Dim);
}
void TSMatrix::Add(Matrix<double>& FEMatr,Vector<DWORD>& 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<double>& 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<double>[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));
}
Список литературы
Зенкевич О., Морган К. Конечные методы и аппроксимация // М.: Мир, 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
|
|
|
