| |||||
МЕНЮ
| Алгоритм компактного хранения и решения СЛАУ высокого порядкасоставляет [pic] байт, что приблизительно равно 2,7 Мбайт оперативной памяти. Размер упакованного представления составил около 315 Кбайт. Данная задача решалась на ЭВМ с процессором Pentium 166 и 32 МБ ОЗУ двумя способами – методом Гаусса и методом Ланцоша. Сопоставление результатов решения приведено в Таблице 1. Таблица 1. | |Время |[pic] |[pic] |[pic] |[pic] |[pic] |[pic] | | |решения | | | | | | | | |(сек) | | | | | | | |Метод |280 |2.2101 |-2.460|1.3756 |-5.2501|1.7406 |-2.3489| |Гаусса | | |8 | | | | | |Метод |150 |2.2137 |-2.466|1.3904 |-5.2572|1.7433 |-2.3883| |Ланцоша | | |9 | | | | | Из Таблицы 1 легко видеть, что результаты решения СЛАУ методом Гаусса и методом Ланцоша хорошо согласуются между собой, при этом время решения вторым способом почти в два раза меньше, чем в случае использования метода Гаусса. ВЫВОДЫ. В данной работе были рассмотрены способы компактного хранения матрицы коэффициентов системы линейных алгебраических уравнений (СЛАУ) и методы ее решения. Разработан алгоритм компактного хранения матрицы жесткости, позволяющий в несколько раз (иногда более чем в десятки раз) сократить объем требуемой памяти для хранения такой матрицы жесткости. Классические методы хранения, учитывающие симметричную и ленточную структуру матриц жесткости, возникающих при применении метода конечных элементов (МКЭ), как правило, не применимы при решении контактных задач, так как при их решении матрицы жесткости нескольких тел объединяются в одну общую матрицу, ширина ленты которой может стремиться к порядку системы. Предложенная в работе методика компактного хранения матриц коэффициентов СЛАУ и использования метода Ланцоша позволили на примере решения контактных задач добиться существенной экономии процессорного времени и затрат оперативной памяти. СПИСОК ССЫЛОК. 1. Зенкевич О., Морган К. Конечные методы и аппроксимация // М.: Мир, 1980 2. Зенкевич О., Метод конечных элементов // М.: Мир., 1975 3. Стрэнг Г., Фикс Дж. Теория метода конечных элементов // М.: Мир, 1977 4. Бахвалов Н.С.,Жидков Н.П., Кобельков Г.М. Численные методы // М.: наука, 1987 5. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления // М.:Наука, 1984 6. Бахвалов Н.С. Численные методы // М.: Наука, 1975 7. Годунов С.К. Решение систем линейных уравнений // Новосибирск: Наука, 1980 8. Гоменюк С.И., Толок В.А. Инструментальная система анализа задач механики деформируемого твердого тела // Приднiпровський науковий вiсник – 1997. – №4. 9. 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 10. А.Джордж, Дж. Лиу, Численное решение больших разреженных систем уравнений // Москва, Мир, 1984 11. D.J. Rose, “A graph theoretic study of the numerical solution of sparse positive definite system of linear equations” // New York, Academic Press, 1972 12. Мосаковский В.И., Гудрамович В.С., Макеев Е.М., Контактные задачи теории оболочек и стержней // М.:”Машиностроение”, 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&,Vec tor&, Matrix&,DWORD&,DWORD&), ReadCubeData(char*,char*,FILE*,FILE*,Vector&,Vector&,Vector< double>&, 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) { |
ИНТЕРЕСНОЕ | |||
|