Алгоритм компактного хранения и решения СЛАУ высокого порядка

Алгоритм компактного хранения и решения СЛАУ высокого порядка ВВЕДЕНИЕ. Метод конечных элементов является численным методом для дифференциальных уравнений, встречающихся в физике 1. Возникновение этого метода связано с решением задач космических исследований 1950 г Впервые он был опубликован в работе Тернера, Клужа, Мартина и Топпа. Эта работа способствовала появлению других работ был опубликован ряд статей с применениями метода конечных элементов к задачам строительной механики и механики сплошных сред. Важный вклад в теоретическую разработку метода сделал в 1963 г. Мелош, который показал, что метод конечных элементов можно рассматривать как один из вариантов хорошо известного метода Рэлея-Ритца. В строительной механике метод конечных элементов минимизацией потенциальной энергии позволяет свести задачу к системе линейных уравнений равновесия 2,3. Одной из существующих трудностей, возникающих при численной реализации решения контактных задач теории упругости методом конечных элементов МКЭ, является решение систем линейных алгебраических уравнений СЛАУ большого порядка вида Большинство существующих методов решения таких систем разработаны в предположении того, что матрица A имеет ленточную структуру, причем ширина ленты , где n2 - порядок . Однако, при использовании МКЭ для численного решения контактных задач возможны случаи, когда ширина ленты 1 ОБЗОР МЕТОДОВ РЕШЕНИЯ СЛАУ, ВОЗНИКАЮЩИХ В МКЭ Основная идея метода конечных элементов состоит в том, что любую непрерывную величину, такую, как температура, давление и перемещение, можно аппроксимировать дискретной моделью, которая строится на множестве кусочно-непрерывных функций, определенных на конечном числе подобластей.

Кусочно-непрерывные функции определяются с помощью значений непрерывной величины в конечном числе точек рассматриваемой области 1,2,3. В общем случае непрерывная величина заранее неизвестна и нужно определить значения этой величины в некоторых внутренних точках области.

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

После этого можно перейти к общему случаю.

Итак, при построении конкретной модели непрерывной величины поступают следующим образом 1. В рассматриваемой области фиксируется конечное число точек. Эти точки называются узловыми точками или просто узлами. 2. Значение непрерывной величины в каждой узловой точке считается переменной, которая должна быть определена. 3. Область определения непрерывной величины разбивается на конечное число подобластей, называемых элементами.

Эти элементы имеют общие узловые точки и в совокупности аппроксимируют форму области. 4 .Непрерывная величина апроксимируется на каждом элементе функцией, которая определяется с помощью узловых значений этой величины.Для каждого элемента определяется своя функция, но функции подбираются таким образом, чтобы сохранялась непрерывность величины вдоль границ элемента.

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

Кроме того, существуют задачи с такой структурой матрицы, для которой прямые методы всегда предпочтительнее, чем итерационные. 1.1 Точные методы решения СЛАУ Рассмотрим ряд точных методов решения СЛАУ 4,5. Решение систем n-линейных уравнении с n-неизвестными по формулам Крамера. Пусть дана система линейных уравнений, в которой число уравнений равно числу неизвестных Предположим, что определитель системы d не равен нулю. Если теперь заменить последовательно в определителе столбцы коэффициентов при неизвестных хj столбцом свободных членов bj, то получатся соответственно n определителей d1 dn. Теорема Крамера.

Система n линейных уравнений с n неизвестными, определитель которой отличен от нуля, всегда совместна и имеет единственное решение, вычисляемое по формулам x1d1d x2d2d xn-1dn-1d xndnd Решение произвольных систем линейных уравнений.Пусть произвольная система линейных уравнений, где число уравнений системы не равно числу n неизвестных.

Предположим, что система 3 совместна и r minm,n, тогда в матрицах А и А найдутся r линейно независимых строк, а остальные m-r строк окажутся их линейными комбинациями. Перестановкой уравнений можно добиться того, что эти r линейно независимых строк займут первые r мест. Отсюда следует, что любое из последних m - r уравнений системы 3 можно представить как сумму первых r уравнений которые называются линейно независимыми или базисными, взятых с некоторыми коэффициентами.

Тогда система эквивалентна следующей системе r уравнений с n неизвестными Предположим, что минор r-го порядка, составленный из коэффициентов при первых r неизвестных, отличен от нуля Мr 0, т. е. является базисным минором. В этом случае неизвестные, коэффициенты при которых составляют базисный минор, называются базисными неизвестными, а остальные n - r - свободными неизвестными.В каждом из уравнений системы 4 перенесем в правую часть все члены со свободными неизвестными xr1 xn. Тогда получим систему, которая содержит r уравнений с r базисными неизвестными.

Так как определитель этой системы есть базисный минор Mr то система имеет единственное решение относительно базисных неизвестных, которое можно найти по формулам Крамера. Давая свободным неизвестным произвольные числовые значения, получим общее решение исходной системы. Однородная система линейных уравнений.Пусть дана однородная система линейных уравнений n неизвестными Так как добавление столбца из нулей не изменяет ранга матрицы системы, то на основании теоремы Кронекера - Kaneлли эта система всегда совместна и имеет, по крайней мере, нулевое решение. Если определитель системы 5 отличен от нуля и число уравнений системы равно числу неизвестных, то по теореме Крамера нулевое решение является единственным.

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

Тогда приходим к системе, решая которую по формулам Крамера, выразим r базисных неизвестных x1 хr через n - r свободных неизвестных. Система 5 имеет бесчисленное множество решений. Среди этого множества есть решения, линейно независимые между собой.Фундаментальной системой решений называются n - r линейно независимых решений однородной системы уравнений.

Метод главных элементов. Пусть дана система n линейных уравнений с n неизвестными расширенная матрица системы 6 . Выберем ненулевой наибольший по модулю и не принадлежащий столбцу свободных членов элемент apq матрицы , который называется главным элементом, и вычислим множители mi-aiqapq для всех строк с номерами i p р - я строка, содержащая главный элемент, называется главной строкой.Далее к каждой неглавной i-й строке прибавим главную строку, умноженную на соответствующий множитель mi для этой строки.

В результате получим новую матрицу, все элементы q-го столбца которой, кроме apq, состоят из нулей. Отбросив этот столбец и главную p-ю получим новую матрицу, число строк и столбцов которой на единицу меньше.Повторяем те же операции с получившейся матрицей, после чего получаем новую матрицу и т.д. Таким образом, построим последовательность матриц, последняя из которых является двучленной матрицей-строкой главной строкой.

Для определения неизвестных xi объединяем в систему все главные строки, начиная с последней. Изложенный метод решения системы линейных уравнений с n неизвестными называется методом главных элементов. Необходимое условие его применения состоит том, что определитель матрицы не равен нулю 6,7. Схема Халецкого. Пусть система линейных уравнений дана в матричном виде. Axb 7 Где А - квадратная матрица порядка n, а x,b - векторы столбцы.Представим матрицу А в виде произведения нижней треугольной матрицы С и верхней треугольной матрицы В с единичной диагональю, т.е. АСВ, Где , Причем элементы сij и bij определяются по формулам , Уравнение 7 можно записать в следующем виде CBxb. 9 Произведение Bx матрицы B на вектор-столбец x является вектором-столбцом, который обозначим через y Bxy. 10 Тогда уравнение 9 перепишем в виде Cyb. 11 Здесь элементы сij известны, так как матрица А системы 7 считается уже разложенной на произведение двух треугольных матриц С и В. Перемножив матрицы в левой части равенства 11, получаем систему уравнений из которой получаем следующие формулы для определения неизвестных неизвестные yi удобно вычислять вместе с элементами bij. После того как все yi определены по формулам 12, подставляем их в уравнение 10. Так как коэффициенты bij определены 8, то значения неизвестных, начиная с последнего, вычисляем по следующим формулам К прямым методам, использующим свойство разреженности А, можно отнести алгоритм минимальной степени, алгоритм минимального дефицита, древовидное блочное разбиение для асимметричного разложения, методы вложенных или параллельных сечений и др. Метод Гаусса.

Пусть дана система Ax b где А матрица размерности m x m. В предположении, что , первое уравнение системы , делим на коэффициент , в результате получаем уравнение Затем из каждого из остальных уравнений вычитается первое уравнение, умноженное на соответствующий коэффициент . В результате эти уравнения преобразуются к виду первое неизвестное оказалось исключенным из всех уравнений, кроме первого.

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

Из -го уравнения системы 2 определяем , из -го уравнения определяем и т.д. до . Совокупность таких вычислений называют обратным ходом метода Гаусса.

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

Процесс построения такой последовательности называется итерационным повторяющимся.

Эффективность применения приближенных методов зависят от выбора начального вектора и быстроты сходимости процесса.Рассмотрим метод итераций метод последовательных приближений. Пусть дана система n линейных уравнений с n неизвестными Ахb, 14 Предполагая, что диагональные элементы aii 0 i 2, n, выразим xi через первое уравнение систем x2 - через второе уравнение и т. д. В результате получим систему, эквивалентную системе 14 Обозначим , где i 1, 2, n j 1,2 n. Тогда система 15 запишется таким образом в матричной форме Решим систему 16 методом последовательных приближений. За нулевое приближение примем столбец свободных членов.

Любое k1-е приближение вычисляют по формуле Если последовательность приближений x0 xk имеет предел , то этот предел является решением системы 15, поскольку в силу свойства предела , т.е. 4,6. Метод Зейделя. Метод Зейделя представляет собой модификацию метода последовательных приближений.В методе Зейделя при вычислении k1-го приближения неизвестного xi i 1 учитываются уже найденные ранее k1-е приближения неизвестных xi-1. Пусть дана линейная система, приведенная к нормальному виду 17 Выбираем произвольно начальные приближения неизвестных и подставляем в первое уравнение системы 17. Полученное первое приближение подставляем во второе уравнение системы и так далее до последнего уравнения.

Аналогично строим вторые, третьи и т.д. итерации.Таким образом, предполагая, что k-е приближения известны, методом Зейделя строим k1-е приближение по следующим формулам где k0,1 n Метод Ланцоша.

Для решения СЛАУ высокого порядка 1, матрица, коэффициентов которой хранится в компактном нижеописанном виде, наиболее удобным итерационным методом является метод Ланцоша 4, схема которого имеет вид 18 где Преимуществом данного метода является его высокая скорость сходимости к точному решению.Кроме того, доказано, что он обладает свойством квадратичного окончания, т.е. для положительно определенной матрицы можно гарантировано получить точное решение при количестве итераций . Размер требуемой памяти на каждой итерации не изменяется, т.к. не требует преобразование матрицы . В качестве критерия остановки данного итерационного процесса обычно используют соотношение , 19 где - заданная точность.

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

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

Такая схема позволяет значительно сократить количество итераций, необходимых для достижения заданной точности 19 или 20 10,11. 2 МЕТОДЫ КОМПАКТНОГО ХРАНЕНИЯ МАТРИЦЫ ЖЕСТКОСТИ Матрица жесткости, получающаяся при применении МКЭ, обладает симметричной структурой, что позволяет в общем случае хранить только верхнюю треугольную часть матрицы.Однако для задач с большим количеством неизвестных это так же приводит к проблеме нехватки памяти.

Предлагаемый в данной работе метод, позволяет хранить только ненулевые члены матрицы жесткости. Суть его заключается в следующем.Первоначально, с целью выявления связей каждого узла с другими, производится анализ структуры дискретизации области на КЭ. Например, для КЭ - сетки, изображенной на рис. 1, соответствующая структура связей будет иметь вид узла1234567Связи1, 2, 5, 6, 71, 2, 3, 62, 3, 4, 63, 4, 5, 6, 71, 4, 5, 71, 2, 3, 4, 6, 71, 4, 5, 6, 7 Тогда, для хранения матрицы жесткости необходимо построчно запоминать информацию о коэффициентах, соответствующих узлам, с которыми связан данный узел. На рис. 2 приведены матрица жесткости и ее компактное представление для сетки изображенной на рис 1 9. Текст подпрограммы, реализующий предложенный алгоритм анализа структуры КЭ-разбиения тела, приведен в Приложении 1. Данный способ компактного хранения матрицы жесткости позволяет легко его использовать совместно с каким-нибудь численным методом.

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

Следовательно, для использования предложенного метода компактного хранения СЛАУ необходимо построить прямое и обратное преобразование в первоначальную квадратную матрицу. Пусть элемент первоначальной квадратной матрицы размерностью , а - ее компактное представление.Тогда для обратного преобразования будут справедливы следующие соотношения , где m количество степеней свободы m1,2,3. Для прямого преобразования будут справедливы соотношения, обратные к соотношениям . 3 ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ Для проверки предлагаемого метода компактного хранения матрицы жесткости была решена задача о контактном взаимодействии оболочечной конструкции и ложемента 12 рис. 4. Данная задача часто возникает на практике при транспортировке или хранении с горизонтальным расположением оси оболочечные конструкции устанавливаются на круговые опоры - ложементы.

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

Данная задача решалась методом конечных элементов при помощи системы FORL 5. Дискретная модель ложемента в трехмерной постановке представлена на Рис. 5. При построении данной КЭ-модели было использовано 880 узлов и 2016 КЭ в форме тетраэдра.

Полный размер матрицы жесткости для такой задачи составляет байт, что приблизительно равно 2,7 Мбайт оперативной памяти. Размер упакованного представления составил около 315 Кбайт. Данная задача решалась на ЭВМ с процессором Pentium 166 и 32 МБ ОЗУ двумя способами методом Гаусса и методом Ланцоша.Сопоставление результатов решения приведено в Таблице 1. Таблица 1. Время решения сек Метод Гаусса2802.2101-2.46081.3756-5.25011.740 6-2.3489Метод Ланцоша1502.2137-2.46691.3904-5.25721.74 33-2.3883 Из Таблицы 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 math.h include stdio.h include string.h include stdlib.h include fstream.h include matrix.h define BASE3D4 4 define BASE3D8 8 define BASE3D10 10 const double Eps 1.0E-10 DWORD CurrentType BASE3D10 void PrintHeadervoid printfCommand format AConvert - switch file1.in file2.in file3.out Optionsn printfSwitch -t10 - Tetraedr10n printf -c8 - Cube8n printf -s4 - Shell4n printf -s8 - Shell8nn printfOptins 8 - convert Tetraedr10- 8Tetraedr4n printf 6 - convert Cube8- 6Tetraedr4n bool Outputchar 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 typeBASE3D4 typeBASE3D10 3 4 DWORD NewSize, i, j ofstream Out if typeBASE3D4 type1 3 else if typeBASE3D8 type1 4 else if typeBASE3D10 type1 6 Out.openfname,iosout ios binary if Out.bad return true Out.writeconst charLabel,6 sizeofchar if Out.fail return true Out.writeconst chartype,sizeofint if Out.fail return true Out.writeconst charCountBn,sizeofDWORD if Out.fail Out.close return true Out.writeconst charNewSize n NumNewPoints,sizeofDWORD if Out.fail return true Out.writeconst charNumNewPoints,sizeofDWORD if Out.fail Out.close return true for DWORD i 0 i n i Out.writeconst charxi,sizeofdouble Out.writeconst charyi,sizeofdouble Out.writeconst charzi,sizeofdouble if Out.fail Out.close return true for i 0 i NumNewPoints i Out.writeconst charxn i,sizeofdouble Out.writeconst charyn i,sizeofdouble if Out.fail Out.close return true Out.writeconst charntr,sizeofDWORD if Out.fail Out.close return true for i 0 i ntr i for j 0 j DWORDtype j DWORD out trij Out.writeconst charout,sizeofDWORD if Out.fail Out.close return true for i 0 i CountBn i for j 0 j DWORDtype1 j DWORD out Boundsij Out.writeconst charout,sizeofDWORD if Out.fail Out.close return true Create Links printfCreate links r Vector DWORD Linkn, Sizen Matrix DWORD Linksn,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 DWORDtype k if trjk i for DWORD m 0 m DWORDtype m Linktrjm 1 Count 0 for DWORD m 0 m n m if Linkm Count Sizei Count Count 0 for DWORD m 0 m n m if Linkm LinksiCount m Set zero Link.ReSizen Output for DWORD i 0 i n i DWORD Sz Sizei Out.writeconst charSz,sizeofDWORD for DWORD j 0 j Sz j Out.writeconst charLinksij,sizeofDWORD printf r printfPoints ldn,n printfFE ldn,ntr Out.close return false bool TestDWORD a,DWORD b bool result int NumPoints 3 if CurrentType BASE3D8 NumPoints 4 else if CurrentType BASE3D10 NumPoints 6 for int i 0 i NumPoints i result false for int j 0 j NumPoints j if bj ai result true break if result false return false return true void ConvertVector double X,Vector double Y,Vector double Z, Matrix DWORD FE,DWORD NumTr,Matrix DWORD Bounds,DWORD BnCount int cData865 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, cData444 0,1,2,3, 1,3,2,0, 3,0,2,1, 0,3,1,2, cData1047 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, cData67, Data6, l, Num1, Num2, m DWORD i, j, p6, pp6, Index Matrix DWORD BoundList4 NumTr,6 double cx, cy, cz, x1, y1, z1, x2, y2, z2, x3, y3, z3 Bounds.ReSize4 NumTr,6 switch CurrentType case BASE3D4 Num1 4 Num2 3 for l 0 l Num1 l for m 0 m Num21 m cDatalm cData4lm break case BASE3D8 Num1 6 Num2 4 for l 0 l Num1 l for m 0 m Num2 1 m cDatalm cData8lm break case BASE3D10 Num1 4 Num2 6 for l 0 l Num1 l for m 0 m Num21 m cDatalm cData10lm printfCreate bounds r for i 0 i NumTr - 1 i for int j 0 j Num1 j if BoundListij for l 0 l Num2 l pl FEicDatajl for DWORD k i 1 k NumTr k for int m 0 m Num1 m if BoundListkm for int l 0 l Num2 l ppl FEkcDataml if Testp,pp BoundListij BoundListkm 1 for i 0 i NumTr i for j 0 j DWORDNum1 j if BoundListij 0 if CurrentType BASE3D4 cx XFEicDataj3 cy YFEicDataj3 cz ZFEicDataj3 else if CurrentType BASE3D10 cx XFEicDataj6 cy YFEicDataj6 cz ZFEicDataj6 else cx XFEicDataj4 cy YFEicDataj4 cz ZFEicDataj4 x1 XFEicDataj0 y1 YFEicDataj0 z1 ZFEicDataj0 x2 XFEicDataj1 y2 YFEicDataj1 z2 ZFEicDataj1 x3 XFEicDataj2 y3 YFEicDataj2 z3 ZFEicDataj2 for l 0 l Num2 l Datal cDatajl if cx-x1y2-y1z3-z1 cy-y1z2-z1x3-x1 y3-y1cz-z1x2-x1 - x3-x1y2-y1cz-z1 - y3-y1z2-z1cx-x1 - cy-y1z3-z1x2-x1 0 if CurrentType BASE3D4 Data0 cDataj0 Data1 cDataj2 Data2 cDataj1 else if CurrentType BASE3D10 Data0 cDataj0 Data1 cDataj2 Data2 cDataj1 Data3 cDataj5 Data5 cDataj3 else Data0 cDataj0 Data1 cDataj3 Data2 cDataj2 Data3 cDataj1 for l 0 l Num2 l BoundsBnCountl FEiDatal BnCount void mainint argc,char argv char input1, input2, input3, op , sw bool CreateFilechar,char,char,char printfANSYS- FORL file convertor.

ZSUc 1998.nn if argc 5 argc 6 PrintHeader return sw argv1 input1 argv2 input2 argv3 input3 argv4 if strcmpsw t10 CurrentType BASE3D10 else if strcmpsw c8 CurrentType BASE3D8 else printfUnknown switch snn,sw PrintHeader return if argc 6 op argv5 if strcmpop,8 strcmpop,6 printfUnknown options snn,op PrintHeader return if CreateFileinput1,input2,input3,op printfOKn bool CreateFilechar fn1,char fn2,char fn3,char Op FILE in1, in2, in3 Vector double X1000, Y1000, Z1000 DWORD NumPoints, NumFE, NumBounds 0, tmp Matrix DWORD FE1000,10, Bounds bool ReadTetraedrDatachar,char,FILE,FILE,Vect or double ,Vector double ,Vector double , Matrix DWORD ,DWORD,DWORD, ReadCubeDatachar,char,FILE,FILE,Vector double ,Vector double ,Vector double , Matrix DWORD ,DWORD, DWORD void Convert824Matrix DWORD ,DWORD, Convert1024Matrix DWORD ,DWORD if in1 fopenfn1,r NULL printfUnable open file s,fn1 return false if in2 fopenfn2,r NULL printfUnable open file s,fn2 return false if CurrentType BASE3D10 if ReadTetraedrDatafn1,fn2,in1,in2,X,Y,Z,FE ,NumPoints,NumFE return false if strcmpOp,8 Create 8Tetraedr4 Convert1024FE,NumFE ConvertX,Y,Z,FE,NumFE,Bounds,NumBounds return Outputfn3,X,Y,Z,FE,NumPoints,0,NumFE,Bou nds,NumBounds if CurrentType BASE3D8 if ReadCubeDatafn1,fn2,in1,in2,X,Y,Z,FE,Num Points,NumFE return false if strcmpOp,6 Create 6Tetraedr4 Convert824FE,NumFE ConvertX,Y,Z,FE, NumFE,Bounds,NumBounds return Outputfn3,X,Y,Z,FE,NumPoints,0,NumFE,Bou nds,NumBounds return false void Convert824Matrix DWORD FE,DWORD NumFE Matrix DWORD nFE6 NumFE,4 DWORD data4 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 nFECurrentk FEidatajk Current CurrentType BASE3D4 NumFE Current FE nFE void Convert1024Matrix DWORD FE,DWORD NumFE Matrix DWORD nFE8 NumFE,4 DWORD data4 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 nFECurrentk FEidatajk Current CurrentType BASE3D4 NumFE Current FE nFE bool ReadTetraedrDatachar 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 TextBuffer1001 DWORD Current 0L, tmp while feofin1 if fgetsTextBuffer,1000,in1 NULL if feofin1 break printfUnable read file s,fn1 fclosein1 fclosein2 return false if sscanfTextBuffer,ld lf lf lf, NumPoints,tx,ty,tz 4 continue XCurrent tx YCurrent ty ZCurrent tz if Current 999 Vector double t1 X, t2 Y, t3 Z X.ReSize2 X.Size Y.ReSize2 X.Size Z.ReSize2 X.Size for DWORD i 0 i Current i Xi t1i Yi t2i Zi t3i if Current 100 0 printfLine ldr,Current fclosein1 printf r NumPoints Current Current 0L while feofin2 if fgetsTextBuffer,1000,in2 NULL if feofin2 break printfUnable read file s,fn2 fclosein2 return false if sscanfTextBuffer,d d d d d ld ld ld ld ld ld ld ld, tmp,tmp,tmp,tmp,tmp, FECurrent0,FECurrent1,FECurrent2,FECurre nt3, FECurrent4,FECurrent5,FECurrent6,FECurre nt7 13 continue if fgetsTextBuffer,1000,in2 NULL printfUnable read file s,fn2 fclosein2 return false if sscanfTextBuffer,ld ld,FECurrent8,FECurrent9 2 printfUnable read file s,fn2 fclosein2 return false if fabstx 0.5XFECurrent0 - 1 XFECurrent1 - 1 - XFECurrent4 - 1 Eps XFECurrent4 - 1 tx if fabsty 0.5YFECurrent0 - 1 YFECurrent1 - 1 - YFECurrent4 - 1 Eps YFECurrent4 - 1 ty if fabstz 0.5ZFECurrent0 - 1 ZFECurrent1 - 1 - ZFECurrent4 - 1 Eps ZFECurrent4 - 1 tz if fabstx 0.5XFECurrent2 - 1 XFECurrent1 - 1 - XFECurrent5 - 1 Eps XFECurrent5 - 1 tx if fabsty 0.5YFECurrent2 - 1 YFECurrent1 - 1 - YFECurrent5 - 1 Eps YFECurrent5 - 1 ty if fabstz 0.5ZFECurrent2 - 1 ZFECurrent1 - 1 - ZFECurrent5 - 1 Eps ZFECurrent5 - 1 tz if fabstx 0.5XFECurrent0 - 1 XFECurrent2 - 1 - XFECurrent6 - 1 Eps XFECurrent6 - 1 tx if fabsty 0.5YFECurrent0 - 1 YFECurrent2 - 1 - YFECurrent6 - 1 Eps YFECurrent6 - 1 ty if fabstz 0.5ZFECurrent0 - 1 ZFECurrent2 - 1 - ZFECurrent6 - 1 Eps ZFECurrent6 - 1 tz if fabstx 0.5XFECurrent0 - 1 XFECurrent3 - 1 - XFECurrent7 - 1 Eps XFECurrent7 - 1 tx if fabsty 0.5YFECurrent0 - 1 YFECurrent3 - 1 - YFECurrent7 - 1 Eps YFECurrent7 - 1 ty if fabstz 0.5ZFECurrent0 - 1 ZFECurrent3 - 1 - ZFECurrent7 - 1 Eps ZFECurrent7 - 1 tz if fabstx 0.5XFECurrent3 - 1 XFECurrent1 - 1 - XFECurrent8 - 1 Eps XFECurrent8 - 1 tx if fabsty 0.5YFECurrent3 - 1 YFECurrent1 - 1 - YFECurrent8 - 1 Eps YFECurrent8 - 1 ty if fabstz 0.5ZFECurrent3 - 1 ZFECurrent1 - 1 - ZFECurrent8 - 1 Eps ZFECurrent8 - 1 tz if fabstx 0.5XFECurrent3 - 1 XFECurrent2 - 1 - XFECurrent9 - 1 Eps XFECurrent9 - 1 tx if fabsty 0.5YFECurrent3 - 1 YFECurrent2 - 1 - YFECurrent9 - 1 Eps YFECurrent9 - 1 ty if fabstz 0.5ZFECurrent3 - 1 ZFECurrent2 - 1 - ZFECurrent9 - 1 Eps ZFECurrent9 - 1 tz if Current 999 Matrix DWORD t FE FE.ReSize2 FE.Size1,10 for DWORD i 0 i Current i for DWORD j 0 j 10 j FEij tij if Current 100 0 printfLine ldr,Current NumFE Current for DWORD i 0 i NumFE i for DWORD j 0 j 10 j FEij printf r return true bool ReadCubeDatachar fn1,charfn2,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 TextBuffer1001 DWORD Current 0L, tmp while feofin1 if fgetsTextBuffer,1000,in1 NULL if feofin1 break printfUnable read file s,fn1 fclosein1 fclosein2 return false if sscanfTextBuffer,ld lf lf lf, NumPoints,tx,ty,tz 4 continue XCurrent tx YCurrent ty ZCurrent tz if Current 999 Vector double t1 X, t2 Y, t3 Z X.ReSize2 X.Size Y.ReSize2 X.Size Z.ReSize2 X.Size for DWORD i 0 i Current i Xi t1i Yi t2i Zi t3i if Current 100 0 printfLine ldr,Current fclosein1 printf r NumPoints Current Current 0L while feofin2 if fgetsTextBuffer,1000,in2 NULL if feofin2 break printfUnable read file s,fn2 fclosein2 return false if sscanfTextBuffer,d d d d d ld ld ld ld ld ld ld ld, tmp,tmp,tmp,tmp,tmp, FECurrent0,FECurrent1,FECurrent3,FECurre nt2, FECurrent4,FECurrent5,FECurrent7,FECurre nt6 13 continue if Current 999 Matrix DWORD t FE FE.ReSize2 FE.Size1,10 for DWORD i 0 i Current i for DWORD j 0 j 10 j FEij tij if Current 100 0 printfLine ldr,Current NumFE Current for DWORD i 0 i NumFE i for DWORD j 0 j 10 j FEij printf r return true ПРИЛОЖЕНИЕ 2. Исходный текст программы, реализующей алгоритм компактного хранения и решения СЛАУ высокого порядка. include matrix.h class RVector private Vector double Buffer public RVectorvoid RVector RVectorDWORD Size Buffer.ReSizeSize RVectorRVector right Buffer right.Buffer RVectorVector double right Buffer right DWORD Sizevoid return Buffer.Size void ReSizeDWORD Size Buffer.ReSizeSize double operator DWORD i return Bufferi RVector operator RVector right Buffer right.Buffer return this RVector operator Vector double right Buffer right return this void SubRVector void SubRVector,double void Muldouble void AddRVector friend double NormRVector,RVector class TSMatrix private Vector double Right Vector double Array Vector DWORD Links uint Dim DWORD Size public TSMatrixvoid Size 0 Dim 0 Array NULL Links NULL TSMatrixVector DWORD ,DWORD,uint TSMatrixvoid if Array delete Array Vector double GetRightvoid return Right DWORD GetSizevoid return Size uint GetDimvoid return Dim Vector double GetVectorDWORD i return Arrayi Vector DWORD GetLinksvoid return Links void SetLinksVector DWORD l Links l void AddMatrix double ,Vector DWORD void AddDWORD I, DWORD L, DWORD J, DWORD K, double v DWORD Row I, Col L LinksI.Size Dim FindI,J Dim K ArrayRowCol v void AddDWORD I, double v RightI v DWORD FindDWORD,DWORD void RestoreMatrix double void SetDWORD,DWORD,double,bool void SetDWORD Index1,DWORD Index2,double value DWORD I Index1 Dim, L Index1 Dim, J Index2 Dim, K Index2 Dim, Pos FindI,J, Row I, Col if Pos DWORD-1 return Col L LinksI.Size Dim FindI,J Dim K ArrayRowCol value bool GetDWORD Index1,DWORD Index2,double value DWORD I Index1 Dim, L Index1 Dim, J Index2 Dim, K Index2 Dim, Pos FindI,J, Row I, Col value 0 if Pos DWORD-1 return false Col L LinksI.Size Dim FindI,J Dim K value ArrayRowCol return true void MulRVector,RVector double MulDWORD,RVector void writeofstream void readifstream class RMatrix private Vector double Buffer DWORD size public RMatrixDWORD sz size sz Buffer.ReSizesizesize 10.5 RMatrix DWORD Sizevoid return size double GetDWORD i, DWORD j return Buffer2size 1 - i0.5i j - i include smatrix.h double NormRVector Left,RVector Right double Ret 0 for DWORD i 0 i Left.Size i Ret Lefti Righti return Ret void RVectorSubRVector Right for DWORD i 0 i Size i thisi - Righti void RVectorAddRVector Right for DWORD i 0 i Size i thisi Righti void RVectorMuldouble koff for DWORD i 0 i Size i thisi koff void RVectorSubRVector Right,double koff for DWORD i 0 i Size i thisi - Rightikoff TSMatrixTSMatrixVector DWORD links, DWORD size, uint dim Dim dim Links links Size size Right.ReSizeDim Size Array new Vector double Size for DWORD i 0 i Size i Arrayi.ReSizeLinksi.Size Dim Dim void TSMatrixAddMatrix 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 FEMatri Dim lj Dim k if Res AddFEi,l,FEj,k,Res for DWORD i 0L i FE.Size i for DWORD l 0L l Dim l RRow FEUINTi FE.Size Dim l Res FEMatri Dim lFEMatr.Size1 if Res AddRRow,Res DWORD TSMatrixFindDWORD I,DWORD J DWORD i for i 0 i LinksI.Size i if LinksIi J return i return DWORD-1 void TSMatrixRestoreMatrix double Matr DWORD i, j, NRow, NPoint, NLink, Pos Matr.ReSizeSize Dim,Size Dim 1 for i 0 i Size i for j 0 j Arrayi.Size j NRow j Arrayi.Size Dim Number of row NPoint j - NRow Arrayi.Size Dim Dim Number of points NLink j Dim Number of link Pos LinksiNPoint Matri Dim NRowPos Dim NLink Arrayij for i 0 i Right.Size i MatriMatr.Size1 Righti void TSMatrixSetDWORD Index,DWORD Position,double Value,bool Case DWORD Row Index, Col Position LinksIndex.Size Dim FindIndex,Index Dim Position, i double koff ArrayRowCol, val if Case RightDim Index Position Value else RightIndex Dim Position Value koff for i 0L i Size Dim i if i Index Dim Position SetIndex Dim Position,i,0 Seti,Index Dim Position,0 if Geti,Index Dim Position,val Righti - val Value void TSMatrixMulRVector Arr,RVector Res DWORD i, j, NRow, NPoint, NLink, Pos Res.ReSizeArr.Size for i 0 i Size i for j 0 j Arrayi.Size j NRow j Arrayi.Size Dim NPoint j - NRow Arrayi.Size Dim Dim NLink j Dim Pos LinksiNPoint Resi Dim NRow ArrPos Dim NLink Arrayij double TSMatrixMulDWORD Index,RVector Arr DWORD j, I Index Dim, L Index Dim, Start L ArrayI.Size Dim, Stop Start ArrayI.Size Dim, NRow, NPoint, NLink, Pos double Res 0 for j Start j Stop j NRow j ArrayI.Size Dim NPoint j - NRow ArrayI.Size Dim Dim NLink j Dim Pos LinksINPoint Res ArrPos Dim NLink ArrayIj return Res void TSMatrixwriteofstream Out DWORD ColSize Out.writecharDim,sizeofDWORD Out.writecharSize,sizeofDWORD for DWORD i 0 i Size i ColSize Arrayi.Size Out.writecharColSize,sizeofDWORD for DWORD j 0 j ColSize j Out.writecharArrayij,sizeofdouble for DWORD i 0 i Size Dim i Out.writecharRighti,sizeofdouble void TSMatrixreadifstream In DWORD ColSize In.readcharDim,sizeofDWORD In.readcharSize,sizeofDWORD if Array delete Array Array new Vector double Size Right.ReSizeSize Dim for DWORD i 0 i Size i In.readcharColSize,sizeofDWORD Arrayi.ReSizeColSize for DWORD j 0 j ColSize j In.readcharArrayij,sizeofdouble for DWORD i 0 i Size Dim i In.readcharRighti,sizeofdouble.