Алгоритм компактного збереження і рішення СЛАУ високого порядка
Аби вирішити СЛАУ в МКЭ потрібно вибрати метод рішення. Остаточне рішення про застосування итерационных або отримання прямих методів рішення СЛАУ необхідно ухвалити з урахуванням аналізу структури досліджуваної математичної завдання. Прямі на методи вирішення СЛАУ вигідніше використовувати, якщо потрібно вирішувати багато однакових систем з різними правими частинами, чи якщо матриця, А є… Читать ещё >
Алгоритм компактного збереження і рішення СЛАУ високого порядка (реферат, курсовая, диплом, контрольная)
Алгоритм компактного збереження і рішення СЛАУ високого порядка.
Метод кінцевих елементів є численным методом для диференційних рівнянь, можна зустріти у фізиці [1]. Виникнення цього методу пов’язані з рішенням завдань космічних досліджень (1950 р.). Вперше він опублікований роботі Тернера, Клужа, Мартіна і Топпа. Ця робота сприяла появі інших робіт; було опубліковано ряд статей з застосуваннями методу кінцевих елементів до завдань будівельної механіки і механіки суцільних середовищ. Важливий внесок у теоретичну розробку методу зробив у 1963 р. Мелош, який довів, що метод кінцевих елементів можна розглядати, як одне із варіантів добре відомого методу Рэлея-Ритца. У будівельної механіці метод кінцевих елементів мінімізацією потенційної енергії дозволяє звести завдання до системи лінійних рівнянь рівноваги [2,3].
Однією з труднощів, які виникають за чисельної реалізації рішення контактних завдань теорії пружності методом кінцевих елементів (МКЭ), є рішення систем лінійних алгебраїчних рівнянь (СЛАУ) великого порядку вида.
[pic].
Більшість існуючих методів вирішення цих систем розроблено у припущенні те, що матриця A має стрічкову структуру, причому ширина стрічки [pic], де n2 — порядок[pic]. Проте, під час використання МКЭ для чисельного рішення контактних завдань можливі випадки, коли ширина стрічки [pic] [5].
1 ОГЛЯД МЕТОДІВ РІШЕННЯ СЛАУ, ВИНИКАЮЧИХ У МКЭ.
Основна ідея методу кінцевих елементів у тому, що будь-яку безперервну величину, таку, як температура, тиск переміщення, можна апроксимувати дискретної моделлю, побудована на безлічі кусочнобезперервних функцій, певних на кінцевому числі подобластей. Кусочнобезперервні функції визначаються з допомогою значень безупинної величини в кінцевому числі точок аналізованої області [1,2,3].
У випадку безперервна величина заздалегідь невідома і треба визначити значення цієї величини у деяких внутрішніх точках області. Дискретну модель, проте, дуже просто побудувати, якщо спочатку припустити, що числові значення цієї величини у кожному внутрішньої точці області відомі. Після цього можна можливість перейти до загальному випадку. Отже, при побудові конкретної моделі безупинної величини надходять наступним образом:
1. У області фіксується кінцеве число точок. Ці точки називаються вузловими точками чи навіть узлами.
2. Значення безупинної величини у кожному вузловий точці вважається перемінної, що має бути определена.
3. Область визначення безупинної величини розбивається на кінцеве число подобластей, званих елементами. Ці елементи мають загальні вузлові крапки й разом аппроксимируют форму области.
4 .Безперервна величина апроксимируется кожному елементі функцією, що визначається з допомогою вузлових значень цієї величини. До кожного елемента визначається своя функція, але функції підбираються в такий спосіб, щоб зберігалася безперервність величини уздовж кордонів элемента.
Аби вирішити СЛАУ в МКЭ потрібно вибрати метод рішення. Остаточне рішення про застосування итерационных або отримання прямих методів рішення СЛАУ необхідно ухвалити з урахуванням аналізу структури досліджуваної математичної завдання. Прямі на методи вирішення СЛАУ вигідніше використовувати, якщо потрібно вирішувати багато однакових систем з різними правими частинами, чи якщо матриця, А є положительно-определенной. З іншого боку, існують завдання з такою структурою матриці, на яку прямі методи завжди краще, ніж итерационные.
1. Точні на методи вирішення СЛАУ.
Розглянемо ряд точних методів рішення СЛАУ [4,5].
Рішення систем n-линейных рівнянні з n-неизвестными по формулам.
Крамера.
Нехай дана система лінійних рівнянь, у якій кількість рівнянь одно числу неизвестных:
[pic].
Припустимо, що визначник системи d не нульовий. Якщо тепер замінити послідовний у определителе стовпчики коефіцієнтів при невідомих хj стовпцем вільних членів bj, то вийдуть відповідно n визначників d1,…, dn.
Теорему Крамера. Система n лінійних рівнянь з n невідомими, визначник якої різниться від нуля, завжди совместна і має єдине рішення, вычисляемое по формулам:
x1=d1/d; x2=d2/d;…; xn-1=dn-1/d; xn=dn/d;
Рішення довільних систем лінійних уравнений.
Пусть.
[pic].
довільна система лінійних рівнянь, де число рівнянь системи не одно числу n невідомих. Припустимо, що систему (3) совместна і r[pic]min{m, n}, тоді матрицях Проте й, А знайдуться r лінійно незалежних рядків, інші ж m-r рядків виявляться їх лінійними комбінаціями. Перестановкою рівнянь можна домогтися, що це r лінійно незалежних рядків займуть перші r мест.
Звідси випливає, що будь-який з усіх m — r рівнянь системи (3) можна подати як суму перших r рівнянь (які називаються лінійно незалежними чи засадничими), узятих з декотрими коефіцієнтами. Тоді система еквівалентна наступній системі r рівнянь з n неизвестными.
[pic].
Припустимо, що мінор r-го порядку, складений із коефіцієнтів при перших r невідомих, різниться від нуля Мr [pic] 0, т. е. є базисним мінором. І тут невідомі, коефіцієнти у яких становлять базисний мінор, називаються засадничими невідомими, інші ж n — r — вільними неизвестными.
У кожному з рівнянь системи (4) перенесемо в праву частина усіх членів зі вільними невідомими xr+1,…, xn. Тоді одержимо систему, яка містить r рівнянь з r засадничими невідомими. Оскільки визначник цієї системи є базисний мінор Mr то система має єдине рішення щодо базисних невідомих, що можна знайти по формулам Крамера. Даючи вільним невідомим довільні числові значення, одержимо загальне рішення вихідної системы.
Однорідна система лінійних уравнений.
Нехай дана однорідна система лінійних рівнянь n неизвестными.
[pic].
Оскільки додавання шпальти з нулів не змінює рангу матриці системи, то, на підставі теореми Кронекера — Kaneлли цю систему завжди совместна і має, по крайнього заходу, нульовий рішення. Якщо визначник системи (5) різниться від нуля і кількість рівнянь системи одно числу невідомих, то теоремі Крамера нульовий рішення є единственным.
У разі, коли ранг матриці системи (5) менше ніж невідомих, т. е. r (А)< n, дана система крім нульового рішення матиме і ненульові рішення. Для перебування цих рішень на системі (5) виділяємо r лінійно незалежних рівнянь, інші відкидаємо. У виділених рівняннях у частині залишаємо r базисних невідомих, інші ж n — r вільних невідомих переносимо в праву частина. Тоді дійшли системі, вирішуючи яку по формулам Крамера, висловимо r базисних невідомих x1,…, хr через n — r вільних неизвестных.
Система (5) має незліченну кількість рішень. Серед цього безлічі має рішення, лінійно незалежні між собой.
Фундаментальною системою рішень називаються n — r лінійно незалежних рішень однорідної системи уравнений.
Метод головних элементов.
Нехай дана система n лінійних рівнянь з n неизвестными.
[pic].
розширена матриця системи (6) [pic]. Виберемо ненульовий найбільший по модулю і приналежний стовпцю вільних членів елемент apq матриці [pic], що називається головним елементом, і обчислимо множники mi=- aiq/apq всім рядків з номерами i[pic]p (р — я рядок, що містить головний елемент, називається головною строкой).
Далі до кожної неглавной i-го рядку додамо головну рядок, помноженій на відповідний множник mi; з цією строки.
Через війну одержимо нову матрицю, все елементи q-го шпальти якої, крім apq, складаються з нулей.
Відкинувши цей стовпець і головну p-ю одержимо нову матрицю, число рядків і шпальт чим одиницю менше. Повторюємо самі операції з получившейся матрицею, після чого отримуємо нову матрицю і т.д.
Отже, побудуємо послідовність матриць, останню з яких є двучленной матрицей-строкой (головною рядком). Для визначення невідомих xi об'єднуємо до системи все головні рядки, починаючи з последней.
Викладене метод рішення системи лінійних рівнянь з n невідомими називається методом головних елементів. Необхідна умова його застосування полягає тому, що визначник матриці не нульовий [6,7].
Схема Халецкого.
Нехай система лінійних рівнянь дана в матричному виде.
Ax=b (7).
Де, А — квадратна матриця порядку n, а x, b — вектори столбцы.
Уявімо матрицю На вигляді твори нижньої трикутною матриці З повагою та верхньої трикутною матриці У з одиничної діагоналлю, т. е.
А=СВ,.
Где.
[pic], [pic].
Причому елементи сij і bij визначаються по формулам:
[pic],.
[pic].
Рівняння (7) можна записати наступного виде:
CBx=b. (9).
Твір Bx матриці B на вектор-столбец x є векторомстовпцем, який позначимо через y:
Bx=y. (10).
Тоді рівняння (9) перепишемо в виде:
Cy=b. (11).
Тут елементи сij відомі, оскільки матриця, А системи (7) вважається вже розкладеною на твір двох трикутних матриць З повагою та В.
Перемноживши матриці у частині рівності (11), отримуємо систему рівнянь з якої отримуємо такі формули визначення неизвестных:
[pic].
невідомі yi зручно обраховувати разом із елементами bij.
Коли все yi визначено по формулам (12), підставляємо в рівняння (10).
Оскільки коефіцієнти bij визначено (8), то значення невідомих, починаючи з останнього, обчислюємо за такими формулам:
[pic].
До прямим методам, використовує властивість розрідженості А, можна віднести: алгоритм мінімальний ступінь, алгоритм мінімального дефіциту, древовидное блокове розбивка для асиметричного розкладання, методи вкладених чи паралельних перетинів і др.
Метод Гаусса.
Нехай дана система.
Ax = b.
де, А — матриця розмірності m x m.
У припущенні, що [pic], перше рівняння системы.
[pic], [pic].
ділимо на коефіцієнт [pic], внаслідок отримуємо уравнение.
[pic].
Потім з кожного з інших рівнянь віднімається перше рівняння, помножена на відповідний коефіцієнт [pic]. Через війну ці рівняння перетворюються до виду.
[pic] перше невідоме виявилося які з всіх рівнянь, крім першого. Далі в припущенні, що [pic], ділимо друге рівняння на коефіцієнт [pic] і виключаємо невідоме із усіх рівнянь, починаючи з другого тощо. Через війну послідовного винятку невідомих система рівнянь перетворюється на систему рівнянь із трикутною матрицей.
[pic].
Сукупність проведених обчислень називається прямим ходом методу Гаусса.
З [pic]-го рівняння системи (2) визначаємо [pic], з ([pic])-го рівняння визначаємо [pic] тощо. до [pic]. Сукупність таких обчислень називають зворотним ходом методу Гаусса.
Реалізація прямого методу Гаусса вимагає [pic] арифметичних операцій, а зворотного — [pic] арифметичних операций.
1.2 Итерационные на методи вирішення СЛАУ.
Метод ітерацій (метод послідовних приближений).
Наближені на методи вирішення систем лінійних рівнянь дозволяють отримувати значення коренів системи із заданою точністю як краю послідовності деяких векторів. Процес побудови такий послідовності називається итерационным (повторяющимся).
Ефективність застосування наближених методів залежить від вибору початкового вектора і швидкості збіжності процесса.
Розглянемо метод ітерацій (метод послідовних наближень). Нехай дана система n лінійних рівнянь з n неизвестными:
Ах=b, (14).
Припускаючи, що діагональні елементи aii [pic] 0 (і = 2, …, n), висловимо xi через перше рівняння систем x2 — через друге рівняння тощо. буд. Через війну одержимо систему, еквівалентну системі (14):
[pic].
Означимо [pic]; [pic], де і == 1, 2, …, n; j == 1,2,…, n. Тоді система (15) запишеться отже, у матричної форме.
[pic].
Вирішимо систему (16) методом послідовних наближень. За нульовий наближення приймемо стовпець вільних членів. Будь-яке (k+1)-е наближення обчислюють по формуле.
[pic].
Якщо послідовність наближень x (0),…, x (k) має межа [pic], цей межа розв’язує системи (15), що у силу властивості краю [pic], тобто. [pic] [4,6].
Метод Зейделя.
Метод Зейделя є модифікацію методу послідовних наближень. У методі Зейделя при обчисленні (k+1)-го наближення невідомого xi (i>1) враховуються вже знайдені раніше (k+1)-е наближення невідомих xi-1.
Нехай дана лінійна система, наведена до нормальному виду:
[pic] (17).
Вибираємо довільно початкові наближення невідомих і підставляємо в перше рівняння системи (17). Отримане перше наближення підставляємо у друге рівняння системи та таке інше аж до останнього рівняння. Аналогічно будуємо другі, треті тощо. итерации.
Отже, припускаючи, що k-е наближення [pic]известны, методом Зейделя будуємо (k+1)-е наближення за такими формулам:
[pic].
[pic].
де k=0,1,…, n.
[pic].
Метод Ланцоша.
Аби вирішити СЛАУ високого порядку (1), матриця, коефіцієнтів якої зберігається в компактному нижеописанном вигляді, найзручнішим итерационным методом є метод Ланцоша [4], схема якого має вид:
[pic] (18).
[pic].
где.
[pic].
Перевагою цього методу є його висока швидкість збіжності до точному рішенню. З іншого боку, доведено, що він має здатність «квадратичного закінчення», тобто. для позитивно певної матриці можна гарантовано отримати точне рішення у кількості ітерацій [pic]. Розмір необхідної пам’яті з кожної ітерації не змінюється, т.к. не вимагає перетворення матриці [pic]. Як критерію зупинки даного итерационного процесу зазвичай використовують соотношение.
[pic], (19).
де [pic]- задана точність. Як іншого критерію збіжності іноді зручніше використовувати среднеквадратичную різницю між рішеннями, отриманими на сусідніх итерациях:
[pic] (20).
Среднеквадратичную різницю необхідно контролювати і під час кожних k наперед заданих итераций.
Окремо слід подивитися на проблему вибору початкового наближення [pic]. Доводиться, що з позитивно певної матриці [pic], итерационный процес (18) завжди сходиться незалежно від виборі початкового наближення. За позитивного рішення контактних завдань, коли для уточнення граничних умов у зони очікуваного контакту потрібно дуже багато рішень СЛАУ виду (1), як початкового наближення на першому розрахунку використовується права частина системи (1), а кожного наступного перерахунку — рішення, отримане попередньому. Така схема дозволяє значно зменшити кількість ітерацій, необхідні досягнення заданої точності (19) чи (20) [10,11].
2 МЕТОДИ КОМПАКТНОГО ЗБЕРІГАННЯ МАТРИЦІ ЖЕСТКОСТИ.
Матриця жорсткості, получающаяся при застосуванні МКЭ, має симетричній структурою, що дозволяє у випадку зберігати лише верхню трикутну частина матриці. Проте задля завдань із велику кількість невідомих це так призводить до проблемі нестачі пам’яті. Запропонований в цій роботі метод, дозволяє зберігати лише ненульові члени матриці жорсткості. Суть його в следующем.
Спочатку, для виявлення зв’язків кожного вузла коїться з іншими, виробляється аналіз структури дискретизації області на КЭ. Наприклад, для КЭ — сітки, зображеною на рис. 1, відповідна структура зв’язків буде мати вид: |№ вузла |1 |2 |3 |4 |5 |6 |7 | |Зв'язки |1, 2, |1, 2, |2, 3, |3, 4, |1, 4, |1, 2, |1, 4, | | |5, 6, 7|3, 6 |4, 6 |5, 6, 7|5, 7 |3, 4, |5, 6, 7| | | | | | | |6, 7 | |.
Тоді, для зберігання матриці жорсткості необхідно через підрядник запам’ятовувати інформацію про коефіцієнти, відповідних вузлам, із якими пов’язаний даний вузол. На рис. 2 наведено матриця жорсткості і його компактне уявлення для сітки зображеною на рис. 1 [9].
Текст підпрограми, який реалізує запропонований алгоритм аналізу структури КЭ-разбиения тіла, приведено у Додатку 1.
Цей спосіб компактного зберігання матриці жорсткості дозволяє легко його використовувати що з якимось численным методом. Найбільш зручним цієї мети видається використання вищевикладеного итерационного методу Ланцоша, бо в кожної ітерації потрібно лише перемножать матрицю коефіцієнтів СЛАУ і поставлене вектор. Отже, від використання запропонованого методу компактного зберігання СЛАУ необхідно побудувати пряме і зворотне перетворення на початкову квадратну матрицу.
Нехай [pic]- елемент початкової квадратної матриці размерностью [pic], а [pic] - її компактне уявлення. Тоді протилежного перетворення будуть справедливі такі соотношения:
[pic], (*).
где m — кількість ступенів свободи (m=1,2,3).
Для прямого перетворення будуть справедливі співвідношення, зворотні до співвідношенням (*).
3 ЧИСЕЛЬНІ ЭКСПЕРИМЕНТЫ.
Для перевірки запропонованого методу компактного зберігання матриці жорсткості було вирішено завдання про контактному взаємодії оболочечной конструкції і ложемента [12] (рис. 4).
Це завдання часто виникає практично при транспортуванні або зберіганні з горизонтальним розташуванням осі оболонкові конструкції встановлюються на кругові опори — ложементи. Взаємодія підкріплених оболочечных конструкцій і ложементів здійснюється через опорні шпангоути, протяжність яких вздовж осі оболонки порівнянна з шириною ложементів і набагато меншою радіуса оболонки, та величини зони контакта.
Це завдання вирішувалася методом кінцевих елементів з допомогою системи FORL [5]. Дискретна модель ложемента (в тривимірної постановці) представлена на Рис. 5.
При побудові даної КЭ-модели було використане 880 вузлів і 2016 КЭ в формі тетраедра. Повний розмір матриці жорсткості такому масштабному завданню становить [pic] байт, що дорівнює 2,7 Мбайт оперативної пам’яті. Розмір упакованого уявлення становить близько 315 Кбайт.
Це завдання вирішувалася на ЕОМ з процесором Pentium 166 і 32 МБ ОЗУ двома шляхами — методом Гаусса і методом Ланцоша. Зіставлення результатів рішення наведено в Таблиці 1.
Таблиця 1. | |Час |[pic] |[pic] |[pic] |[pic] |[pic] |[pic] | | |рішення | | | | | | | | |(сік) | | | | | | | |Метод |280 |2.2101 |-2.460|1.3756 |-5.2501|1.7406 |-2.3489| |Гаусса | | |8 | | | | | |Метод |150 |2.2137 |-2.466|1.3904 |-5.2572|1.7433 |-2.3883| |Ланцоша | | |9 | | | | |.
З Таблиці 1 то зрозуміло, що результати рішення СЛАУ методом Гаусса і методом Ланцоша добре узгоджуються між собою, у своїй час вирішення другим способом майже двічі менше, ніж у випадку використання методу Гаусса.
ВЫВОДЫ.
У цьому роботі було розглянуто способи компактного зберігання матриці коефіцієнтів системи лінійних алгебраїчних рівнянь (СЛАУ) й методи його рішення. Розроблено алгоритм компактного зберігання матриці жорсткості, дозволяє у кілька разів (іноді - більш ніж у десятки раз) скоротити обсяг необхідної пам’яті для зберігання такої матриці жесткости.
Класичні методи зберігання, враховують симметричную і стрічкову структуру матриць жорсткості, які виникають за застосуванні методу кінцевих елементів (МКЭ), зазвичай, неспроможні під час вирішення контактних завдань, бо за розв’язанні матриці жорсткості кількох тіл об'єднують у одну загальну матрицю, ширина стрічки якій у змозі йти до порядку системы.
Запропонована у роботі методика компактного зберігання матриць коефіцієнтів СЛАУ та ефективного використання методу Ланцоша дозволили з прикладу рішення контактних завдань домогтися істотну економію процесорного часу й витрат оперативної памяти.
СПИСОК ССЫЛОК.
1. Зенкевич Про., Морган До. Кінцеві методи лікування й апроксимація // М.: Мир,.
2. Зенкевич Про., Метод кінцевих елементів // М.: Світ., 1975.
3. Стрэнг Р., Фікс Дж. Теорія методу кінцевих елементів // М.: Мир,.
4. Хвальків Н.С., Жидков Н. П., Кобельков Г. М. Чисельні методи // М.: наука, 1987.
5. Воєводін В.В., Кузнєцов Ю.О. Матриці і обчислення // М.:Наука, 1984.
6. Хвальків М.С. Чисельні методи // М.: Наука, 1975.
7. Годунов С. К. Рішення систем лінійних рівнянь // Новосибирск:
Наука, 1980.
8. Гоменюк С.І., Толок В. А. Інструментальна система аналізу завдань механіки деформируемого твердого тіла // Приднiпровський науковий вiсник — 1997. — № 4.
9. F.G. Gustavson, «Some basic techniques for solving sparse matrix algorithms», // editer by D.J. Rose and R.A.Willoughby, Plenum.
Press, New York, 1972.
10. А. Джордж, Дж. Ліу, Кількісна рішення великих розріджених систем рівнянь // Москва, Світ, 1984.
11. D.J. Rose, «A graph theoretic study of the numerical solution of sparse positive definite system of linear equations» // New York,.
Academic Press, 1972.
12. Мосаковский В.І., Гудрамович В. С., Макєєв О.М., Контактні завдання теорії оболонок і стрижнів // М.:"Машиностроение", 1978.
ДОДАТОК 1.
Вихідний текст програми, який реалізує аналіз структури КЭ-разбиения объекта.
#include.
#include.
#include.
#include.
#include.
#include «matrix.h «.
#define BASE3D4 4.
#define BASE3D8 8.
#define BASE3D10 10.
const double Eps = 1.0E-10;
DWORD CurrentType = BASE3D10;
void PrintHeader (void).
{ printf («Command format: AConvert — [/Options]n »); printf («Switch: -t10 — Tetraedr (10)n »); printf («-c8 — Cube (8)n »); printf («-s4 — Shell (4)n »); printf («-s8 — Shell (8)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& x, Vector& y, Vector& z, Matrix& tr, DWORD n,.
DWORD NumNewPoints, DWORD ntr, Matrix& Bounds, DWORD CountBn).
{ char* Label = «NTRout »; int type = CurrentType, type1 = (type==BASE3D4 || type==BASE3D10)? 3: 4;
DWORD NewSize, і, j; ofstream Out;
if (type==BASE3D4) type1 = 3; else if (type==BASE3D8) type1 = 4; else if (type==BASE3D10) 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 і = 0; і < 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 (і = 0; і < 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 (і = 0; і < ntr; і++) 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 (і = 0; і < CountBn; і++) for (j = 0; j < (DWORD)type1; j++).
{.
DWORD out = Bounds[i][j];
Out.write ((const char*)&out, sizeof (DWORD)); if (Out.fail ()).
{.
Out.close (); return true;
}.
}.
{.
//*********************.
// Create Links printf («Create links… r »);
Vector Link (n),.
Size (n);
Matrix Links (n, n);
DWORD Count; int type = CurrentType;
for (DWORD і = 0; і < n; i++).
{ for (DWORD j = 0; j < ntr; j++) for (DWORD k = 0; k < (DWORD)type; k++) if (tr[j][k] == і) 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 і = 0; і < 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 == BASE3D8) NumPoints = 4; else if (CurrentType == BASE3D10) NumPoints = 6;
for (int і = 0; і < NumPoints; i++).
{ result = false; for (int j = 0; j < NumPoints; j++) if (b[j] == a[i]).
{ result = true; break;
} if (result == false) return false;
} return true;
}.
void Convert (Vector& X, Vector& Y, Vector& Z, Matrix& FE, DWORD NumTr, Matrix& Bounds, DWORD& BnCount).
{ int cData8[6][5] = {{0,4,5,1,7},.
{6,2,3,7,0},.
{4,6,7,5,0},.
{2,0,1,3,5},.
{1,5,7,3,4},.
{6,4,0,2,1}}, cData4[4][4] = {{0,1,2,3},.
{1,3,2,0},.
{3,0,2,1},.
{0,3,1,2}}, cData10[4][7] = {{0,1,2,4,5,6,3},.
{0,1,3,4,8,7,2},.
{1,3,2,8,9,5,0},.
{0,2,3,6,9,7,1}}, cData[6][7],.
Data[6], l,.
Num1,.
Num2, m;
DWORD і, j, p[6], pp[6],.
Index;
Matrix BoundList (4 * NumTr, 6); double cx, cy, cz, x1, y1, z1, x2, y2, z2, x3, y3, z3;
Bounds.ReSize (4 * NumTr, 6); switch (CurrentType).
{ case BASE3D4:
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 BASE3D8:
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 BASE3D10:
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 (і = 0; і < NumTr — 1; і++) 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 = і + 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 (і = 0; і < NumTr; і++) for (j = 0; j < (DWORD)Num1; j++) if (BoundList[i][j] == 0).
{ if (CurrentType == BASE3D4).
{ cx = X[FE[i][cData[j][3]]]; cy = Y[FE[i][cData[j][3]]]; cz = Z[FE[i][cData[j][3]]];
} else if (CurrentType == BASE3D10).
{ 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) — (cyy1)*(z3-z1)*(x2-x1)) > 0).
{ if (CurrentType == BASE3D4).
{.
Data[0] = cData[j][0];
Data[1] = cData[j][2];
Data[2] = cData[j][1];
} else if (CurrentType == BASE3D10).
{.
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© 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 = BASE3D10; else if (!strcmp (sw, «-c8 »)).
CurrentType = BASE3D8; 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 X (1000),.
Y (1000),.
Z (1000);
DWORD NumPoints,.
NumFE,.
NumBounds = 0, tmp;
Matrix FE (1000,10),.
Bounds; bool ReadTetraedrData (char*, char*, FILE*, FILE*, Vector&, Vector&, Vec tor&,.
Matrix&, DWORD&, DWORD&),.
ReadCubeData (char*, char*, FILE*, FILE*, Vector&, Vector&, Vector< double>&,.
Matrix&, DWORD&, DWORD&); void Convert824(Matrix&, DWORD&),.
Convert1024(Matrix&, DWORD&);
if ((in1 = fopen (fn1, «r »)) == NULL).
{ printf («Unable open file %p.s », fn1); return false;
} if ((in2 = fopen (fn2, «r »)) == NULL).
{ printf («Unable open file %p.s », fn2); return false;
}.
if (CurrentType == BASE3D10).
{ 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 == BASE3D8).
{ if (!ReadCubeData (fn1,fn2,in1,in2,X, Y, Z, FE, NumPoints, NumFE)) return false; if (!strcmp (Op, «/6 »)).
{.
// Create 6*Tetraedr (4).
Convert824(FE, NumFE);
}.
Convert (X, Y, Z, FE, NumFE, Bounds, NumBounds); return! Output (fn3,X, Y, Z, FE, NumPoints, 0, NumFE, Bounds, NumBounds);
} return false;
}.
void Convert824(Matrix& FE, DWORD& NumFE).
{.
Matrix nFE (6 * NumFE, 4);
DWORD data[][4] = {.
{ 0,2,3,6 },.
{ 4,5,1,7 },.
{ 0,4,1,3 },.
{ 6,7,3,4 },.
{ 1,3,7,4 },.
{ 0,4,6,3 }.
},.
Current = 0;
for (DWORD і = 0; і < NumFE; і++) for (DWORD j = 0; j < 6; j++).
{ for (DWORD k = 0; k < 4; k++) nFE[Current][k] = FE[i][data[j][k]];
Current++;
}.
CurrentType = BASE3D4;
NumFE = Current;
FE = nFE;
}.
void Convert1024(Matrix& FE, DWORD& NumFE).
{.
Matrix nFE (8 * NumFE, 4);
DWORD data[][4] = {.
{ 3,7,8,9 },.
{ 0,4,7,6 },.
{ 2,5,9,6 },.
{ 7,9,8,6 },.
{ 4,8,5,1 },.
{ 4,5,8,6 },.
{ 7,8,4,6 },.
{ 8,9,5,6 }.
},.
Current = 0;
for (DWORD і = 0; і < NumFE; і++) for (DWORD j = 0; j < 8; j++).
{ for (DWORD k = 0; k < 4; k++) nFE[Current][k] = FE[i][data[j][k]];
Current++;
}.
CurrentType = BASE3D4;
NumFE = Current;
FE = nFE;
}.
bool ReadTetraedrData (char* fn1, char* fn2, FILE* in1, FILE* in2, Vector& X, Vector& Y, Vector& Z,.
Matrix& FE, DWORD& NumPoints, DWORD& NumFE).
{ double tx, ty, tz; char TextBuffer[1001];
DWORD Current = 0L, tmp;
while (!feof (in1)).
{ if (fgets (TextBuffer, 1000, in1) == NULL).
{ if (feof (in1)) break; printf («Unable read file %p.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 t1 = X, t2 = Y, t3 = Z;
X.ReSize (2 * X. Size ());
Y.ReSize (2 * X. Size ());
Z.ReSize (2 * X. Size ()); for (DWORD і = 0; і < 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 %p.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 %p.s », fn2); fclose (in2); return false;
} if (sscanf (TextBuffer, «%ld %ld » ,&FE[Current][8],&FE[Current][9]) ≠ 2).
{ printf («Unable read file %p.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 t = FE;
FE.ReSize (2 * FE. Size1(), 10); for (DWORD і = 0; і < Current; і++) 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 і = 0; і < NumFE; і++) 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& X, Vector& Y, Vector& Z,.
Matrix& FE, DWORD& NumPoints, DWORD& NumFE).
{ double tx, ty, tz; char TextBuffer[1001];
DWORD Current = 0L, tmp;
while (!feof (in1)).
{ if (fgets (TextBuffer, 1000, in1) == NULL).
{ if (feof (in1)) break; printf («Unable read file %p.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 t1 = X, t2 = Y, t3 = Z;
X.ReSize (2 * X. Size ());
Y.ReSize (2 * X. Size ());
Z.ReSize (2 * X. Size ()); for (DWORD і = 0; і < 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 %p.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 t = FE;
FE.ReSize (2 * FE. Size1(), 10); for (DWORD і = 0; і < Current; і++) 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 і = 0; і < NumFE; і++) for (DWORD j = 0; j < 10; j++).
FE[i][j]—; printf («r »); return true;}.
ДОДАТОК 2.
Вихідний текст програми, реалізує алгоритм компактного збереження і рішення СЛАУ високого порядка.
#include «matrix.h «.
class RVector.
{ private:
Vector Buffer; public:
RVector (void) {}.
~RVector () {}.
RVector (DWORD Size) { Buffer. ReSize (Size); }.
RVector (RVector& right) { Buffer = right. Buffer; }.
RVector (Vector& right) { Buffer = right; }.
DWORD Size (void) { return Buffer. Size (); } void ReSize (DWORD Size) { Buffer. ReSize (Size); } double& operator [] (DWORD і) { return Buffer[i]; }.
RVector& operator = (RVector& right) { Buffer = right. Buffer; return *this; }.
RVector& operator = (Vector& 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 Right;
Vector* Array;
Vector* Links; uint Dim;
DWORD Size; public:
TSMatrix (void) { Size = 0; Dim = 0; Array = NULL; Links = NULL; }.
TSMatrix (Vector*, DWORD, uint);
~TSMatrix (void) { if (Array) delete [] Array; }.
Vector& GetRight (void) { return Right; }.
DWORD GetSize (void) { return Size; } uint GetDim (void) { return Dim; }.
Vector& GetVector (DWORD і) { return Array[i]; }.
Vector* GetLinks (void) { return Links; } void SetLinks (Vector* l) { Links = l; } void Add (Matrix&, Vector&); 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&); 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 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 — і]; }.
};
//************************.
#include «smatrix.h «.
double Norm (RVector& Left, RVector& Right).
{ double Ret = 0;
for (DWORD і = 0; і < Left. Size (); i++).
Ret += Left[i] * Right[i]; return Ret;
}.
void RVector: Sub (RVector& Right).
{ for (DWORD і = 0; і < Size (); i++).
(*this)[i] -= Right[i];
}.
void RVector: Add (RVector& Right).
{ for (DWORD і = 0; і < Size (); i++).
(*this)[i] += Right[i];
}.
void RVector: Mul (double koff).
{ for (DWORD і = 0; і < Size (); i++).
(*this)[i] *= koff;
}.
void RVector: Sub (RVector& Right, double koff).
{ for (DWORD і = 0; і < Size (); i++).
(*this)[i] -= Right[i]*koff;
}.
TSMatrix:TSMatrix (Vector* links, DWORD size, uint dim).
{.
Dim = dim;
Links = links;
Size = size;
Right.ReSize (Dim * Size);
Array = new Vector[Size]; for (DWORD і = 0; і < Size; i++).
Array[i]. ReSize (Links[i].Size () * Dim * Dim);
}.
void TSMatrix: Add (Matrix& FEMatr, Vector& FE).
{ double Res;
DWORD RRow;
for (DWORD і = 0L; і < FE. Size (); і++) 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 і = 0L; і < FE. Size (); і++) 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 (і = 0; і < Links[I]. Size (); і++) if (Links[I][i] == J) return і; return DWORD (-1);
}.
void TSMatrix: Restore (Matrix& Matr).
{.
DWORD і, j,.
NRow,.
NPoint,.
NLink,.
Pos;
Matr.ReSize (Size * Dim, Size * Dim + 1); for (і = 0; і < Size; і++) for (j = 0; 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 (і = 0; і < Right. Size (); і++) 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, і; double koff = Array[Row][Col], val;
if (!Case).
Right[Dim * Index + Position] = Value; else.
{.
Right[Index * Dim + Position] = Value * koff; for (і = 0L; і < Size * Dim; і++) if (і ≠ 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 і, j,.
NRow,.
NPoint,.
NLink,.
Pos;
Res.ReSize (Arr.Size ()); for (і = 0; і < Size; і++) for (j = 0; 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 і = 0; і < 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 і = 0; і < 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[Size];
Right.ReSize (Size * Dim); for (DWORD і = 0; і < 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 і = 0; і < Size * Dim; i++).
In.read ((char*)&(Right[i]), sizeof (double));
}.
———————————;
Малюнок 2.
[pic].
[pic].
Рис. 5.
Малюнок 1.