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

ВВЕДЕНИЕ.

Метод конечных элементов является численным методом длядифференциальных уравнений, встречающихся в физике [1]. Возникновение этогометода связано с решением задач космических исследований (1950 г.). Впервые онбыл опубликован в работе Тернера, Клужа, Мартина и Топпа. Эта работа способствовалапоявлению других работ; был опубликован ряд статей  с применениями методаконечных элементов к задачам строительной механики и механики сплошных сред.Важный вклад в теоретическую разработку метода сделал в 1963 г. Мелош, которыйпоказал, что метод конечных элементов можно рассматривать как один из вариантовхорошо известного метода Рэлея-Ритца. В строительной механике метод конечныхэлементов минимизацией потенциальной  энергии позволяет свести задачу к системелинейных уравнений равновесия [2,3].

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

/>

Большинство существующих методов решения таких системразработаны в предположении того, что матрица A имеет ленточную структуру,причем ширина ленты />, где n2 — порядок/>. Однако, прииспользовании МКЭ для численного решения контактных задач возможны случаи,когда ширина ленты /> [5].

1 ОБЗОР МЕТОДОВ РЕШЕНИЯ СЛАУ, ВОЗНИКАЮЩИХВ МКЭ

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

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

1. В рассматриваемой области фиксируется конечное числоточек. Эти точки называются узловыми точками или просто узлами.

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

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

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

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

Точные методырешения СЛАУ

Рассмотрим ряд точных методов решения СЛАУ [4,5].

Решение систем n-линейных уравнении с n-неизвестными поформулам Крамера.

Пусть дана система линейных уравнений, в которой числоуравнений равно числу неизвестных:

/>

Предположим, что определитель системы d не равен нулю.Если теперь заменить последовательно в определителе столбцы коэффициентов принеизвестных хj столбцом свободных членов bj, то получатся соответственно nопределителей d1,...,dn.

Теорема Крамера. Система n линейных уравнений с nнеизвестными, определитель которой отличен от нуля, всегда совместна и имеетединственное решение, вычисляемое по формулам:

x1=d1/d; x2=d2/d;....; xn-1=dn-1/d; xn=dn/d;

Решение произвольных систем линейных уравнений.

Пусть

/>

произвольная система линейных уравнений, где числоуравнений системы не равно числу n неизвестных. Предположим, что система (3)совместна и r/>min{m,n},тогда в матрицах А и А найдутся r линейно независимых строк, а остальные m-rстрок окажутся их линейными комбинациями. Перестановкой уравнений можнодобиться того, что эти r линейно независимых строк займут первые r мест.

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

/>

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

В каждом из уравнений системы (4) перенесем в правуючасть все члены со свободными неизвестными xr+1,..., 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=-aiq/apq для всех строк сномерами i/>p(р — я строка, содержащая главный элемент, называется главной строкой).

Далее к каждой неглавной i-й строке прибавим главнуюстроку, умноженную на соответствующий множитель mi; для этой строки.

В результате получим новую матрицу, все элементы q-гостолбца которой, кроме apq, состоят из нулей.

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

Таким образом, построим последовательность матриц,последняя из которых является двучленной матрицей-строкой (главной строкой).Для определения неизвестных xi объединяем в систему все главные строки, начинаяс последней.

Изложенный метод решения системы линейных уравнений с nнеизвестными называется методом главных элементов. Необходимое условие егоприменения состоит том, что определитель матрицы не равен нулю [6,7].

Схема Халецкого.

Пусть система линейных уравнений дана в матричном виде.

Ax=b     (7)

Где А — квадратная матрица порядка n, а x,b — векторыстолбцы.

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

А=СВ,

Где

/>,  />

Причем элементы сij  и bij определяются по формулам:

/>,

/>

Уравнение (7) можно записать в следующем виде:

CBx=b.          (9)

Произведение Bx матрицы B на вектор-столбец x являетсявектором-столбцом, который обозначим через y:

Bx=y.        (10)

Тогда уравнение (9) перепишем  в виде:

Cy=b.     (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) методом последовательных приближений.За нулевое приближение примем столбец свободных членов. Любое (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 /> <td/> />
Тогда, для хранения матрицыжесткости необходимо построчно запоминать информацию о коэффициентах,соответствующих  узлам, с которыми связан данный узел.  На рис. 2 приведены  матрица жесткости и ее компактное представление для сетки изображенной на рис 1[9].

/>

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

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

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

/>,                                                                       (*)

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

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

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

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

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

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

/>

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

Даннаязадача решалась на ЭВМ с процессором Pentium166  и 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 — convertTetraedr(10)->8*Tetraedr(4)\n»);

  printf("        /6 — convertCube(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,DWORDntr,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((constchar*)&type,sizeof(int));

  if (Out.fail()) return true;

  Out.write((constchar*)&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((constchar*)&(NumNewPoints),sizeof(DWORD));

  if (Out.fail())

   {

     Out.close();

     return true;

   }

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

   {

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

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

     Out.write((constchar*)&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((constchar*)&(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((constchar*)&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((constchar*)&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((constchar*)&Sz,sizeof(DWORD));

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

      Out.write((constchar*)&(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 fileconvertor. 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)    { returnBuffer.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 doubleNorm(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(DWORDi) { 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, DWORDK, 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,DWORDIndex2,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,DWORDIndex2,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,doublekoff)

{

  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);

}

voidTSMatrix::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,DWORDPosition,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(DWORDIndex,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 matrixalgorithms”, // 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 ofsparse positive definite system of linear equations” // New York, AcademicPress, 1972

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

Для подготовки данной работы были использованы материалыс сайта www.ed.vseved.ru/

еще рефераты
Еще работы по математике