Загрузка...
Учебная работа. Реферат: Алгоритм компактного хранения и решения СЛАУ высокого порядка
способ конечных частей является численным способом для
дифференциальных уравнений, встречающихся в физике [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)
совместна и rmin{m,n},
тогда в матрицах А и А найдутся r линейно независящих строк, а другие m-r
строк окажутся их линейными комбинациями. Перестановкой уравнений можно
достигнуть того, что эти r линейно независящих строк займут 1-ые 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 для всех строк с
номерами ip
(р — я строчка, содержащая основной элемент, именуется главной строчкой).
Дальше к каждой неглавной 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) считается уже разложенной на произведение 2-ух треугольных матриц С и В.
Перемножив матрицы в левой части равенства (11),
получаем систему уравнений из которой получаем последующие формулы для
определения неведомых:
неведомые yi комфортно вычислять вкупе с элементами bij.
Опосля того как все yi определены по формулам (12),
подставляем их в уравнение(10).
Потому что коэффициенты bij определены (8), то значения
неведомых, начиная с крайнего, вычисляем по последующим формулам:
К прямым способам, использующим свойство разреженности А,
можно отнести: метод малой степени, метод малого недостатка,
древовидное блочное разбиение для асимметричного разложения, способы вложенных
либо параллельных сечений и др.
способ Гаусса.
Пусть дана система
Ax
= b
где А – матрица размерности m x m.
В предположении, что , 1-ое уравнение системы
,
делим на коэффициент , в итоге получаем уравнение
Потом из всякого из других уравнений вычитается
1-ое уравнение, умноженное на соответственный коэффициент . В итоге эти
уравнения преобразуются к виду
1-ое неведомое оказалось исключенным из всех
уравнений, не считая первого. Дальше в предположении, что , делим 2-ое уравнение на
коэффициент и
исключаем неведомое из всех уравнений, начиная со второго и т.д. В итоге
поочередного исключения неведомых система уравнений преобразуется в
систему уравнений с треугольной матрицей
совокупа
проведенных вычислений именуется прямым ходом способа Гаусса.
Из -го уравнения системы (2)
определяем ,
из ()-го
уравнения определяем и т.д. до . совокупа таковых вычислений
именуют оборотным ходом способа Гаусса.
Реализация прямого способа Гаусса просит арифметических операций,
а оборотного — арифметических операций.
1.2. Итерационные
способы решения СЛАУ
Способ итераций (способ поочередных приближений).
Приближенные способы решения систем линейных уравнений
разрешают получать значения корней системы с данной точностью в виде предела
последовательности неких векторов. процесс построения таковой
последовательности именуется итерационным (циклическим).
Эффективность внедрения приближенных способов зависят от
выбора исходного вектора и быстроты сходимости процесса.
Разглядим способ итераций (способ поочередных
приближений). Пусть дана система n линейных уравнений с n неведомыми:
Ах=b, (14)
Предполагая, что диагональные элементы aii 0 (i = 2, …,
n), выразим xi через 1-ое уравнение систем x2 — через 2-ое уравнение и т.
д. В итоге получим систему, эквивалентную системе (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)
Избираем произвольно исходные приближения неведомых и
подставляем в 1-ое уравнение системы (17). Приобретенное 1-ое приближение
подставляем во 2-ое уравнение системы и так дальше до крайнего уравнения.
Аналогично строим 2-ые, третьи и т.д. итерации.
Таковым образом, предполагая, что k-е приближения известны,
способом Зейделя строим (k+1)-е приближение по последующим формулам:
где k=0,1,…,n
способ Ланцоша.
Для решения СЛАУ высочайшего порядка (1), матрица,
коэффициентов которой хранится в малогабаритном нижеописанном виде, более
комфортным итерационным способом является способ Ланцоша [4], схема которого имеет
вид:
(18)
где
Преимуществом данного способа является его высочайшая
скорость сходимости к четкому решению. Не считая того, подтверждено, что он владеет
свойством «квадратичного окончания», т.е. для положительно определенной матрицы
можно гарантировано получить четкое решение при количестве итераций . Размер
требуемой памяти на каждой итерации не меняется, т.к. не просит
преобразование матрицы . В качестве аспекта остановки
данного итерационного процесса обычно употребляют соотношение
, (19)
где — данная точность. В качестве
другого аспекта сходимости время от времени удобнее применять среднеквадратичную
разность меж решениями, приобретенными на примыкающих итерациях:
(20)
Среднеквадратичную разность нужно надзирать
при выполнении каждых k наперед данных итераций.
Раздельно следует разглядеть делему выбора исходного
приближения .
Доказывается, что при положительно определенной матрице , итерационный процесс
(18) постоянно сходится при любом выборе исходного приближения. При решении
контактных задач, когда для уточнения граничных критерий в зоне предполагаемого
контакта требуется огромное количество решений СЛАУ вида (1), в качестве
исходного приближения для первого расчета употребляется правая часть системы
(1), а для всякого следующего пересчета — решение, приобретенное на прошлом.
Таковая схема дозволяет существенно уменьшить количество итераций, нужных
для заслуги данной точности (19) либо (20) [10,11].
2 способы
КОМПАКТНОГО ХРАНЕНИЯ МАТРИЦЫ ЖЕСТКОСТИ
Матрица
жесткости, получающаяся при применении МКЭ, владеет симметричной структурой,
что дозволяет в общем случае хранить лишь верхнюю треугольную часть матрицы.
Но для задач с огромным количеством неведомых это так же приводит к
дилемме нехватки памяти. Предлагаемый в данной работе способ, дозволяет хранить
лишь ненулевые члены матрицы жесткости. Сущность его заключается в последующем.
Сначало,
с целью выявления связей всякого узла с иными, делается анализ структуры
дискретизации области на КЭ. к примеру, для КЭ — сетки, изображенной на рис. 1,
соответственная структура связей будет иметь вид:
№ узла
1
2
3
4
5
6
7
Связи
1, 2, 5, 6, 7
1, 2, 3, 6
2, 3, 4, 6
3, 4, 5, 6, 7
1, 4, 5, 7
1, 2, 3, 4, 6, 7
1, 4, 5, 6, 7
Тогда, для хранения матрицы
жесткости нужно построчно запоминать информацию о коэффициентах,
соответственных узлам, с которыми связан данный узел. На рис. 2 приведены
матрица жесткости и ее малогабаритное текст
подпрограммы, реализующий предложенный метод анализа структуры КЭ-разбиения
тела, приведен в Приложении 1.
Данный
метод малогабаритного хранения матрицы жесткости дозволяет просто его применять
вместе с каким-либо численным способом. Более комфортным для данной цели
представляется внедрение вышеизложенного итерационного способа Ланцоша, так
как на каждой итерации требуется лишь перемножать матрицу коэффициентов СЛАУ
и данный вектор. Как следует, для использования предложенного способа
малогабаритного хранения СЛАУ нужно выстроить прямое и оборотное преобразование
в первоначальную квадратную матрицу.
Пусть
– элемент
начальной квадратной матрицы размерностью , а — ее малогабаритное количество степеней свободы (m=1,2,3).
Для
прямого преобразования будут справедливы соотношения, оборотные к соотношениям
.
3 ЧИСЛЕННЫЕ
ЭКСПЕРИМЕНТЫ
Для
проверки предлагаемого способа малогабаритного хранения матрицы жесткости была
решена задачка о контактном содействии оболочечной конструкции и ложемента
[12] (рис. 4).
Данная задачка нередко возникает на
практике при транспортировке либо хранении с горизонтальным расположением оси
оболочечные конструкции инсталлируются на радиальные опоры — ложементы.
Взаимодействие подкрепленных оболочечных конструкций и ложементов
осуществляется через опорные шпангоуты, протяженность которых вдоль оси
оболочки соизмерима с шириной ложементов и много меньше радиуса оболочки и
величины зоны контакта.
Данная
задачка решалась способом конечных частей с помощью системы FORL [5]. Дискретная модель ложемента (в трехмерной постановке)
представлена на Рис. 5.
При
построении данной КЭ-модели было применено 880 узлов и 2016 КЭ в форме
тетраэдра. Полный размер матрицы жесткости для таковой задачки составляет б, что
примерно равно 2,7 Мбайт оперативки. Размер упакованного
представления составил около 315 Кбайт.
Данная
задачка решалась на ЭВМ с микропроцессором Pentium
166 и 32 МБ ОЗУ 2-мя методами – способом Гаусса и способом Ланцоша.
Сравнение результатов решения приведено в Таблице 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 просто созидать, что результаты решения СЛАУ способом Гаусса и способом
Ланцоша отлично согласуются меж собой, при всем этом время решения вторым методом
практически вдвое меньше, чем в случае использования способа Гаусса.
ВЫВОДЫ.
В
данной работе были рассмотрены методы малогабаритного хранения матрицы
коэффициентов системы линейных алгебраических уравнений (СЛАУ) и способы ее
решения. Разработан метод малогабаритного хранения матрицы жесткости,
позволяющий в несколько раз (время от времени наиболее чем в 10-ки раз) уменьшить размер
требуемой памяти для хранения таковой матрицы жесткости.
Традиционные
способы хранения, учитывающие симметричную и ленточную структуру матриц
жесткости, возникающих при применении способа конечных частей (МКЭ), как
правило, не применимы при решении контактных задач, потому что при их решении
матрицы жесткости нескольких тел соединяются воединыжды в одну общую матрицу, ширина
ленты которой может стремиться к порядку системы.
Предложенная
в работе методика малогабаритного хранения матриц коэффициентов СЛАУ и
использования способа Ланцоша дозволили на примере решения контактных задач
достигнуть значимой экономии процессорного времени и издержек оперативной
памяти.
ПРИЛОЖЕНИЕ 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)nn»);
printf(«Optins: /8 — convert
Tetraedr(10)->8*Tetraedr(4)n»);
printf(» /6 — convert
Cube(8)->6*Tetraedr(4)n»);
}
bool Output(char*
fname,Vector<double>& x,Vector<double>&
y,Vector<double>& z, Matrix<DWORD>& tr, DWORD n,
DWORD NumNewPoints,DWORD
ntr,Matrix<DWORD>& Bounds,DWORD CountBn)
{
char* Label = «NTRout»;
int type = CurrentType,
type1 = (type==BASE3D_4 ||
type==BASE3D_10) ? 3 : 4;
DWORD NewSize,
i,
j;
ofstream Out;
if (type==BASE3D_4) type1 = 3;
else if (type==BASE3D_8) type1 = 4;
else if (type==BASE3D_10) type1 = 6;
Out.open(fname,ios::out | ios:: binary);
if (Out.bad()) return true;
Out.write((const char*)Label,6 *
sizeof(char));
if (Out.fail()) return true;
Out.write((const
char*)&type,sizeof(int));
if (Out.fail()) return true;
Out.write((const
char*)&CountBn,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
Out.write((const char*)&(NewSize = n +
NumNewPoints),sizeof(DWORD));
if (Out.fail()) return true;
Out.write((const
char*)&(NumNewPoints),sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
for (DWORD i = 0; i < n; i++)
{
Out.write((const
char*)&x[i],sizeof(double));
Out.write((const
char*)&y[i],sizeof(double));
Out.write((const
char*)&z[i],sizeof(double));
if (Out.fail())
{
Out.close();
return true;
}
}
for (i = 0; i < NumNewPoints; i++)
{
Out.write((const char*)&x[n +
i],sizeof(double));
Out.write((const char*)&y[n +
i],sizeof(double));
if (Out.fail())
{
Out.close();
return true;
}
}
Out.write((const
char*)&(ntr),sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
for (i = 0; i < ntr; i++)
for (j = 0; j < (DWORD)type; j++)
{
DWORD out = tr[i][j];
Out.write((const
char*)&out,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
}
for (i = 0; i < CountBn; i++)
for (j = 0; j < (DWORD)type1; j++)
{
DWORD out = Bounds[i][j];
Out.write((const
char*)&out,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
}
{
//*********************
// Create Links
printf(«Create links…r»);
Vector<DWORD> Link(n),
Size(n);
Matrix<DWORD> Links(n,n);
DWORD Count;
int type = CurrentType;
for (DWORD i = 0; i < n; i++)
{
for (DWORD j = 0; j < ntr; j++)
for (DWORD k = 0; k < (DWORD)type;
k++)
if (tr[j][k] == i)
for (DWORD m = 0; m <
(DWORD)type; m++) Link[tr[j][m]] = 1;
Count = 0;
for (DWORD m = 0; m < n; m++)
if (Link[m]) Count++;
Size[i] = Count;
Count = 0;
for (DWORD m = 0; m < n; m++)
if (Link[m])
Links[i][Count++] = m;
//Set zero
Link.ReSize(n);
}
// Output
//*********************
for (DWORD i = 0; i < n; i++)
{
DWORD Sz = Size[i];
Out.write((const
char*)&Sz,sizeof(DWORD));
for (DWORD j = 0; j < Sz; j++)
Out.write((const
char*)&(Links[i][j]),sizeof(DWORD));
}
//*********************
}
printf(»
r»);
printf(«Points: %ldn»,n);
printf(«FE: %ldn»,ntr);
Out.close();
return false;
}
bool Test(DWORD* a,DWORD* b)
{
bool result;
int NumPoints = 3;
if (CurrentType == BASE3D_8) NumPoints =
4;
else if (CurrentType == BASE3D_10)
NumPoints = 6;
for (int i = 0; i < NumPoints; i++)
{
result = false;
for (int j = 0; j < NumPoints; j++)
if (b[j] == a[i])
{
result = true;
break;
}
if (result == false) return false;
}
return true;
}
void Convert(Vector<double>&
X,Vector<double>& Y,Vector<double>& Z,
Matrix<DWORD>& FE,DWORD NumTr,Matrix<DWORD>&
Bounds,DWORD& BnCount)
{
int cData8[6][5] = {{0,4,5,1,7},
{6,2,3,7,0},
{4,6,7,5,0},
{2,0,1,3,5},
{1,5,7,3,4},
{6,4,0,2,1}},
cData4[4][4] = {{0,1,2,3},
{1,3,2,0},
{3,0,2,1},
{0,3,1,2}},
cData10[4][7] =
{{0,1,2,4,5,6,3},
{0,1,3,4,8,7,2},
{1,3,2,8,9,5,0},
{0,2,3,6,9,7,1}},
cData[6][7],
Data[6],
l,
Num1,
Num2,
m;
DWORD i,
j,
p[6],
pp[6],
Index;
Matrix<DWORD> BoundList(4 * NumTr,6);
double cx,
cy,
cz,
x1,
y1,
z1,
x2,
y2,
z2,
x3,
y3,
z3;
Bounds.ReSize(4 * NumTr,6);
switch (CurrentType)
{
case BASE3D_4:
Num1 = 4;
Num2 = 3;
for (l = 0; l < Num1; l++)
for (m = 0; m < Num2+1; m++)
cData[l][m] = cData4[l][m];
break;
case BASE3D_8:
Num1 = 6;
Num2 = 4;
for (l = 0; l < Num1; l++)
for (m = 0; m < Num2 + 1; m++)
cData[l][m] = cData8[l][m];
break;
case BASE3D_10:
Num1 = 4;
Num2 = 6;
for (l = 0; l < Num1; l++)
for (m = 0; m < Num2+1; m++)
cData[l][m] = cData10[l][m];
}
printf(«Create bounds…r»);
for (i = 0; i < NumTr — 1; i++)
for (int j = 0; j < Num1; j++)
if (!BoundList[i][j])
{
for (l = 0; l < Num2; l++)
p[l] = FE[i][cData[j][l]];
for (DWORD k = i + 1; k < NumTr;
k++)
for (int m = 0; m < Num1; m++)
if (!BoundList[k][m])
{
for (int l = 0; l < Num2; l++)
pp[l] = FE[k][cData[m][l]];
if (Test(p,pp))
BoundList[i][j] = BoundList[k][m]
= 1;
}
}
for (i = 0; i < NumTr; i++)
for (j = 0; j < (DWORD)Num1; j++)
if (BoundList[i][j] == 0)
{
if (CurrentType == BASE3D_4)
{
cx = X[FE[i][cData[j][3]]];
cy = Y[FE[i][cData[j][3]]];
cz = Z[FE[i][cData[j][3]]];
}
else
if (CurrentType == BASE3D_10)
{
cx = X[FE[i][cData[j][6]]];
cy = Y[FE[i][cData[j][6]]];
cz = Z[FE[i][cData[j][6]]];
}
else
{
cx = X[FE[i][cData[j][4]]];
cy = Y[FE[i][cData[j][4]]];
cz = Z[FE[i][cData[j][4]]];
}
x1 = X[FE[i][cData[j][0]]];
y1 = Y[FE[i][cData[j][0]]];
z1 = Z[FE[i][cData[j][0]]];
x2 = X[FE[i][cData[j][1]]];
y2 = Y[FE[i][cData[j][1]]];
z2 = Z[FE[i][cData[j][1]]];
x3 = X[FE[i][cData[j][2]]];
y3 = Y[FE[i][cData[j][2]]];
z3 = Z[FE[i][cData[j][2]]];
for (l = 0; l < Num2; l++)
Data[l] = cData[j][l];
if ( ((cx-x1)*(y2-y1)*(z3-z1) +
(cy-y1)*(z2-z1)*(x3-x1) + (y3-y1)*(cz-z1)*(x2-x1) —
(x3-x1)*(y2-y1)*(cz-z1) -
(y3-y1)*(z2-z1)*(cx-x1) — (cy-y1)*(z3-z1)*(x2-x1)) > 0)
{
if (CurrentType == BASE3D_4)
{
Data[0] = cData[j][0];
Data[1] = cData[j][2];
Data[2] = cData[j][1];
}
else
if (CurrentType == BASE3D_10)
{
Data[0] = cData[j][0];
Data[1] = cData[j][2];
Data[2] = cData[j][1];
Data[3] = cData[j][5];
Data[5] = cData[j][3];
}
else
{
Data[0] = cData[j][0];
Data[1] = cData[j][3];
Data[2] = cData[j][2];
Data[3] = cData[j][1];
}
}
for (l = 0; l < Num2; l++)
Bounds[BnCount][l] = FE[i][Data[l]];
BnCount++;
}
}
void main(int argc,char** argv)
{
char *input1,
*input2,
*input3,
*op = «»,
*sw;
bool CreateFile(char*,char*,char*,char*);
printf(«ANSYS->FORL file
convertor. ZSU(c) 1998.nn»);
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
%snn»,sw);
PrintHeader();
return;
}
if (argc == 6)
{
op = argv[5];
if (strcmp(op,»/8″)
&& strcmp(op,»/6″))
{
printf(«Unknown options
%snn»,op);
PrintHeader();
return;
}
}
if (CreateFile(input1,input2,input3,op))
printf(«OKn»);
}
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:
%ldr»,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:
%ldr»,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:
%ldr»,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:
%ldr»,Current);}
NumFE = Current;
for (DWORD i = 0; i < NumFE; i++)
for (DWORD j = 0; j < 10; j++)
FE[i][j]—;
printf(»
r»);
return true;}ПРИЛОЖЕНИЕ 2.
Начальный текст программки, реализующей метод
малогабаритного хранения и решения СЛАУ высочайшего порядка.
#include «matrix.h»
class RVector
{
private:
Vector<double> Buffer;
public:
RVector(void) {}
~RVector() {}
RVector(DWORD Size) {
Buffer.ReSize(Size); }
RVector(RVector& right) { Buffer =
right.Buffer; }
RVector(Vector<double>& right)
{ Buffer = right; }
DWORD Size(void) { return
Buffer.Size(); }
void ReSize(DWORD Size) {
Buffer.ReSize(Size); }
double& operator [] (DWORD i) {
return Buffer[i]; }
RVector& operator = (RVector&
right) { Buffer = right.Buffer; return *this; }
RVector& operator =
(Vector<double>& right) { Buffer = right; return *this; }
void Sub(RVector&);
void Sub(RVector&,double);
void Mul(double);
void Add(RVector&);
friend double
Norm(RVector&,RVector&);
};
class TSMatrix
{
private:
Vector<double> Right;
Vector<double>* Array;
Vector<DWORD>* Links;
uint Dim;
DWORD Size;
public:
TSMatrix(void) { Size = 0; Dim = 0; Array
= NULL; Links = NULL; }
TSMatrix(Vector<DWORD>*,DWORD,uint);
~TSMatrix(void) { if (Array) delete []
Array; }
Vector<double>& GetRight(void)
{ return Right; }
DWORD GetSize(void) { return Size; }
uint GetDim(void) { return Dim; }
Vector<double>& GetVector(DWORD
i) { return Array[i]; }
Vector<DWORD>* GetLinks(void) {
return Links; }
void SetLinks(Vector<DWORD>* l) {
Links = l; }
void
Add(Matrix<double>&,Vector<DWORD>&);
void Add(DWORD I, DWORD L, DWORD J, DWORD
K, double v)
{
DWORD Row = I,
Col = L * Links[I].Size() * Dim +
Find(I,J) * Dim + K;
Array[Row][Col] += v;
}
void Add(DWORD I, double v)
{
Right[I] += v;
}
DWORD Find(DWORD,DWORD);
void Restore(Matrix<double>&);
void Set(DWORD,DWORD,double,bool);
void Set(DWORD Index1,DWORD
Index2,double value)
{
DWORD I = Index1 / Dim,
L = Index1 % Dim,
J = Index2 / Dim,
K = Index2 % Dim,
Pos = Find(I,J),
Row = I,
Col;
if (Pos == DWORD(-1)) return;
Col = L * Links[I].Size() * Dim +
Find(I,J) * Dim + K;
Array[Row][Col] = value;
}
bool Get(DWORD Index1,DWORD
Index2,double& value)
{
DWORD I = Index1 / Dim,
L = Index1 % Dim,
J = Index2 / Dim,
K = Index2 % Dim,
Pos = Find(I,J),
Row = I,
Col;
value = 0;
if (Pos == DWORD(-1)) return false;
Col = L * Links[I].Size() * Dim +
Find(I,J) * Dim + K;
value = Array[Row][Col];
return true;
}
void Mul(RVector&,RVector&);
double Mul(DWORD,RVector&);
void write(ofstream&);
void read(ifstream&);
};
class RMatrix
{
private:
Vector<double> Buffer;
DWORD size;
public:
RMatrix(DWORD sz) { size = sz;
Buffer.ReSize(size*(size + 1)*0.5); }
~RMatrix() {}
DWORD Size(void) { return size; }
double& Get(DWORD i,DWORD j) {
return Buffer[(2*size + 1 — i)*0.5*i + j — i]; }
};
//************************
#include «smatrix.h»
double Norm(RVector& Left,RVector&
Right)
{
double Ret = 0;
for (DWORD i = 0; i < Left.Size(); i++)
Ret += Left[i] * Right[i];
return Ret;
}
void RVector::Sub(RVector& Right)
{
for (DWORD i = 0; i < Size(); i++)
(*this)[i] -= Right[i];
}
void RVector::Add(RVector& Right)
{
for (DWORD i = 0; i < Size(); i++)
(*this)[i] += Right[i];
}
void RVector::Mul(double koff)
{
for (DWORD i = 0; i < Size(); i++)
(*this)[i] *= koff;
}
void RVector::Sub(RVector& Right,double
koff)
{
for (DWORD i = 0; i < Size(); i++)
(*this)[i] -= Right[i]*koff;
}
TSMatrix::TSMatrix(Vector<DWORD>*
links, DWORD size, uint dim)
{
Dim = dim;
Links = links;
Size = size;
Right.ReSize(Dim * Size);
Array = new Vector<double>[Size];
for (DWORD i = 0; i < Size; i++)
Array[i].ReSize(Links[i].Size() * Dim *
Dim);
}
void TSMatrix::Add(Matrix<double>&
FEMatr,Vector<DWORD>& FE)
{
double Res;
DWORD RRow;
for (DWORD i = 0L; i < FE.Size(); i++)
for (DWORD l = 0L; l < Dim; l++)
for (DWORD j = 0L; j < FE.Size();
j++)
for (DWORD k = 0L; k < Dim; k++)
{
Res = FEMatr[i * Dim +
l][j * Dim + k];
if (Res) Add(FE[i],l,FE[j],k,Res);
}
for (DWORD i = 0L; i < FE.Size(); i++)
for (DWORD l = 0L; l < Dim; l++)
{
RRow = FE[UINT(i %
(FE.Size()))] * Dim + l;
Res = FEMatr[i * Dim +
l][FEMatr.Size1()];
if (Res) Add(RRow,Res);
}
}
DWORD TSMatrix::Find(DWORD I,DWORD J)
{
DWORD i;
for (i = 0; i < Links[I].Size(); i++)
if (Links[I][i] == J) return i;
return DWORD(-1);
}
void
TSMatrix::Restore(Matrix<double>& Matr)
{
DWORD i,
j,
NRow,
NPoint,
NLink,
Pos;
Matr.ReSize(Size * Dim,Size * Dim + 1);
for (i = 0; i < Size; i++)
for (j = 0; j < Array[i].Size(); j++)
{
NRow = j / (Array[i].Size() /
Dim); // Number of row
NPoint = (j — NRow * (Array[i].Size()
/ Dim)) / Dim; // Number of points
NLink = j %
Dim; // Number of link
Pos = Links[i][NPoint];
Matr[i * Dim + NRow][Pos * Dim +
NLink] = Array[i][j];
}
for (i = 0; i < Right.Size(); i++)
Matr[i][Matr.Size1()] = Right[i];
}
void TSMatrix::Set(DWORD Index,DWORD
Position,double Value,bool Case)
{
DWORD Row = Index,
Col = Position *
Links[Index].Size() * Dim + Find(Index,Index) * Dim + Position,
i;
double koff = Array[Row][Col],
val;
if (!Case)
Right[Dim * Index + Position] = Value;
else
{
Right[Index * Dim + Position] = Value *
koff;
for (i = 0L; i < Size * Dim; i++)
if (i != Index * Dim + Position)
{
Set(Index * Dim + Position,i,0);
Set(i,Index * Dim + Position,0);
if (Get(i,Index * Dim +
Position,val))
Right[i] -= val * Value;
}
}
}
void TSMatrix::Mul(RVector&
Arr,RVector& Res)
{
DWORD i,
j,
NRow,
NPoint,
NLink,
Pos;
Res.ReSize(Arr.Size());
for (i = 0; i < Size; i++)
for (j = 0; j < Array[i].Size(); j++)
{
NRow = j / (Array[i].Size() / Dim);
NPoint = (j — NRow * (Array[i].Size()
/ Dim)) / Dim;
NLink = j % Dim;
Pos = Links[i][NPoint];
Res[i * Dim + NRow] += Arr[Pos * Dim +
NLink] * Array[i][j];
}
}
double TSMatrix::Mul(DWORD
Index,RVector& Arr)
{
DWORD j,
I = Index / Dim,
L = Index % Dim,
Start = L * (Array[I].Size() /
Dim),
Stop = Start + (Array[I].Size() /
Dim),
NRow,
NPoint,
NLink,
Pos;
double Res = 0;
for (j = Start; j < Stop; j++)
{
NRow = j / (Array[I].Size() / Dim);
NPoint = (j — NRow * (Array[I].Size()
/ Dim)) / Dim;
NLink = j % Dim;
Pos = Links[I][NPoint];
Res += Arr[Pos * Dim + NLink] *
Array[I][j];
}
return Res;
}
void TSMatrix::write(ofstream& Out)
{
DWORD ColSize;
Out.write((char*)&(Dim),sizeof(DWORD));
Out.write((char*)&(Size),sizeof(DWORD));
for (DWORD i = 0; i < Size; i++)
{
ColSize = Array[i].Size();
Out.write((char*)&(ColSize),sizeof(DWORD));
for (DWORD j = 0; j < ColSize; j++)
Out.write((char*)&(Array[i][j]),sizeof(double));
}
for (DWORD i = 0; i < Size * Dim; i++)
Out.write((char*)&(Right[i]),sizeof(double));
}
void TSMatrix::read(ifstream& In)
{
DWORD ColSize;
In.read((char*)&(Dim),sizeof(DWORD));
In.read((char*)&(Size),sizeof(DWORD));
if (Array) delete [] Array;
Array = new Vector<double>[Size];
Right.ReSize(Size * Dim);
for (DWORD i = 0; i < Size; i++)
{
In.read((char*)&(ColSize),sizeof(DWORD));
Array[i].ReSize(ColSize);
for (DWORD j = 0; j < ColSize; j++)
In.read((char*)&(Array[i][j]),sizeof(double));
}
for (DWORD i = 0; i < Size * Dim; i++)
In.read((char*)&(Right[i]),sizeof(double));
}
Перечень литературы
Зенкевич
О., Морган К. Конечные способы и аппроксимация // М.: мир, 1980
Зенкевич
О., Способ конечных частей // М.: мир., 1975
Стрэнг
Г., Фикс Дж. Теория способа конечных частей // М.: мир, 1977
Бахвалов
Н.С.,Жидков Н.П., Кобельков Г.М. Численные способы // М.: наука, 1987
Воеводин
В.В., Кузнецов Ю.А. Матрицы и вычисления // М.:Наука, 1984
Бахвалов
Н.С. Численные способы // М.: Наука, 1975
Годунов
С.К. Решение систем линейных уравнений // Новосибирск: Наука, 1980
Гоменюк
С.И., Толок В.А. Инструментальная система анализа задач механики деформируемого
твердого тела // Приднiпровський
науковий вiсник – 1997. – №4.
F.G. Gustavson, “Some basic techniques for solving sparse matrix
algorithms”, // editer by D.J. Rose and R.A.Willoughby, Plenum Press, New York,
1972
А.Джордж,
Дж. Лиу, Численное решение огромных разреженных систем уравнений // Москва, мир,
1984
D.J. Rose, “A graph theoretic study of the numerical solution of
sparse positive definite system of linear equations” // New York, Academic
Press, 1972
Мосаковский
В.И., Гудрамович В.С., Макеев Е.М., Контактные задачки теории оболочек и
стержней // М.:”Машиностроение”, 1978
]]>