Розробка алгоритму та програми чисельного розвязку систем лінійних алгебраїчних рівнянь з розрідженою (стр. 8 из 9)
Була розглянута задача про пружне деформування тонкостінної просторової рами (рис. 4.1) при наступних геометричних параметрах: L=H=2м, L1=1,98м, L2=0,01м. Пружні властивості матеріалу рами: E=106Па, ν=0,3. Зовнішнє навантаження P=102Н.
Кінцево-елементна модель усієї конструкції містила 1449604 вузла і 4374320 скінчених елементів у формі тетраедра. СЛАР розв’язувалася з використанням першої схеми. Час розв’язку склав близько двох годин.
4.2 Контактна взаємодія оболонкової конструкції і ложемента
Була вирішена задача про контактну взаємодію оболонкової конструкції і ложемента (рис. 4.2).
Дана задача часто виникає на практиці. При транспортуванні або зберіганні з горизонтальним розташуванням осі оболонкові конструкції встановлюються на кругові опори – ложементи. Взаємодія підкріплених оболонкових конструкцій і ложементів здійснюється через опорні шпангоути, довжина яких уздовж вісі оболонки порівняна з шириною ложементів і багато менше радіуса оболонки і величини зони контакту.
Дискретна модель ложемента (у тривимірній постановці) представлена на рис. 4.3.
При побудові даної КЕ-моделі було використано 6689вузлів і 15323 КЕ у формі тетраедра. Розмір нестиснутої матриці жорсткості для такої задачі становить (6689*3)2*8=3221475912 байт, що приблизно дорівнює 3ГБ оперативної пам'яті.
Задача розв’язуваласядвома способами – методом Гауса і методом Ланцоша. Порівняння результатів розв’язку наведено в табл. 4.2.
Таблиця 4.2 – Результати розв’язку
| Метод Гауса | Метод Ланцоша | |
| Час розв’язку, (сек.) | 3137 | 2194 |
З табл. 4.2 легко бачити, що час розв’язку методом Ланцоша значно менше, ніж у випадку використання методу Гауса.
У даній роботі розглянуті способи компактного зберігання матриці коефіцієнтів системи лінійних алгебраїчних рівнянь і методи розв’язку СЛАР. Розглянуті алгоритми компактного зберігання матриці жорсткості дозволяють значно скоротити об'єм необхідної пам'яті для зберігання матриці жорсткості.
Класичні методи зберігання, що враховують симетричну і стрічкову структуру матриць жорсткості, що виникають при застосуванні методу скінчених елементів, як правило, не застосовні при рішенні контактних задач, тому що при їхньому рішенні матриці жорсткості декількох тіл поєднуються в одну загальну матрицю, ширина стрічки якої може прагнути до розмірності системи.
Викладені у роботі методи компактного зберігання матриці коефіцієнтів СЛАР дозволяють на прикладі розв’язку контактних задач досягти значної економії оперативної пам'яті. Це дозволяє розв’язувати системи з мільйонами невідомих на сучасних персональних комп'ютерах.
Вважаю, що для оптимальної продуктивності при розв’язку СЛАР високих ступенів треба використовувати 24 ГБ оперативної пам'яті стандарту DDRIII 1333, мікропроцесор Intel Core i7-965 ExtremeEdition, системну плату ASUS RampageIIExtreme, дві відеокарти AMD Radeon HD 4870X2 та твердотільний накопичувачOCZSummit. Потужність такої конфігурації може скласти близько 5 терафлопс, що досить непогано. Для пришвидшення розрахунків вважаю доцільним розгін мікропроцесора, використання у програмі спеціалізованих наборів інструкцій процесорів по роботі з числами з плаваючою крапкою та використання у програмі ресурсів відеопроцесора відеокарти, за умови оптимізації програми під такі нововведення. Також значно пришвидшити розрахунки допоможе використання більшої кількості процесорів (або процесорних ядер).
Наведеніалгоритмирозв’язку СЛАР високих ступенів у задачах механіки можуть ефективно застосовуватися для аналізу тривимірних задач високої розмірності, що особливо актуально в задачах сучасного машинобудування.
1. Зенкевич О., Морган К. Конечные методы и аппроксимация. – М.: Мир, 1980. – 479 с.
2. Зенкевич О. Метод конечных элементов. – М.: Мир, 1975. – 217 с.
3. Стрэнг Г., Фикс Дж. Теория метода конечных элементов. – М.: Мир, 1977. – 346 с.
4. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. – М.: Наука, 1987. – 632 с.
5. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. – М.: Наука, 1984. – 176 с.
6. Бахвалов Н.С. Численные методы. – М.: Наука, 1975. – 153 с.
7. Годунов С.К. Решение систем линейных уравнений. – Новосибирск: Наука, 1980. – 81 с.
8. Gustavson F.G. Some basic techniques for solving sparse matrix algorithms. – New York: Plenum Press, 1972. – 215 p.
9. Джордж А., Лиу Дж. Численное решение больших разрежённых систем уравнений. – М.: Мир, 1984. – 333 с.
10. Rose D.J. A graph theoretic study of the numerical solution of sparse positive definite system of linear equations. – New York: Academic Press, 1972. – 473 p.
11. Мосаковский В.И., Гудрамович В.С., Макеев Е.М. Контактные задачи теории оболочек и стержней. – М.: Машиностроение, 1978. – 718 с.
12. Рвачев В.Л., Шевченко А.Н. Проблемно-ориентированные языки и системы для инженерных расчётов. – К.: Техніка, 1988. – 198 с.
13. Голуб Дж., Ван Лоун Ч. Матричные вычисления. – М.: Мир, 1999. – 548 с.
14. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. – М.: Наука, 1986. – 288 с.
15. Тьюарсон Р. Разрежённые матрицы. – М.: Мир, 1977. – 191 с.
16. Брамеллер А., Аллан Р., Хэмэм Я. Слабозаполненные матрицы. – М.: Энергия, 1979. – 192 с.
17. Эстербю О., Златев 3. Прямые методы для разрежённых матриц. – М.: Мир, 1987. – 120 с.
18. Писсанецки С. Технология разреженных матриц. – М.: Мир, 1988. – 410 с.
19. Сабоннадьер Ж.-К., Кулон Ж.-Л. Метод конечных элементов и САПР. – М.: Мир, 1989. – 192 с.
20. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. – М.: Наука, 1978. – 592 с.
21. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. – М.: Мир, 1991. – 347 с.
22. Валях Е. Последовательно-параллельные вычисления. – М.: Мир, 1985. – 216 с.
23. Таненбаум Э., Ван Стеен М. Распределённые системы. Принципы и парадигмы. – СПб.: Питер, 2003. – 481 с.
Вихідний текст програми розв’язку СЛАР з розрідженою матрицею
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <fstream.h>
#include "matrix.h"
#include "smatrix.h"
#define BASE3D_4 4
#define BASE3D_8 8
#define BASE3D_10 10
const double Eps=1.0E-10;
DWORD CurrentType=BASE3D_10;
class RVector
{
private:
Vector<double> Buffer;
public:
RVector() {}
~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;
int Dim;
DWORD Size;
public:
TSMatrix(void) { Size=0; Dim=0; Array=NULL; Links=NULL; }
TSMatrix(Vector<DWORD> *, DWORD, int);
~TSMatrix(void) { if (Array) delete [] Array; }
Vector<double> &GetRight(void) { return Right; }
DWORD GetSize(void) { return Size; }
int 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]; }
};
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, int 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[INT(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);
NPoint=(j-NRow*(Array[i].Size()/Dim))/Dim;
NLink=j%Dim;
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;