Для прямого перетворення будуть слушні співвідношення, зворотні до співвідношень (2.1) – (2.4).
Розглянемо як приклад кінцево-елементну модель, зображену на рис. 2.3. З кожним вузлом сітки зв'язане деяке число вузлів, які можуть бути розбиті на дві множини: першу множину становлять вузли, номери яких менше i – номера розглядаємого вузла, а друга множина – вузли, номери яких більше i. У табл. 2.2 зображена структура матриці СЛАР, відповідна до вузлових зв'язків кінцево-елементної моделі, представленої на рис. 2.3.
У рівнянні (2.5) присутні суми добутків двох типів. Перша сума складається з добутків елементів aij
– матриці A
і елементів вектора невідомих
, у яких j<i
, а для добутків другої суми маємо j>i
. Для побудови добутків першої суми використовуються вектори
і
, приклад структури яких показаний на рис. 2.5. У векторі
зберігаються номери тих вузлів, які пов'язані з даним розглядаємим вузлом, і їх номера менше номера даного вузла. Наприклад, розглянемо вузол з номером 5 (див. рис. 2.3). Цей вузол входить в елементи з номерами 2 і 3. У силу цього він має зв'язок з вузлами 4, 2 і 1, номера яких менше 5. Номера цих вузлів записуються у вектор
у послідовності, визначеній номерами елементів і правилом обходу їх вузлів. На рис. 2.3 стрілками показаний початок і напрямок обходу вузлів елементів. Цим пояснюється порядок появи номерів 4, 2 і 1 у векторі
. У векторі
зберігаються номери індексів, що координують у векторі
розташування номера першого вузла всього масиву номерів вузлів, пов'язаних з розглядаємим вузлом. Наприклад, як видно на рис. 2.3, масив вузлів, пов'язаних з вузлом 5, починається з індексу 6. Останній індекс масиву визначається початковим індексом вузла, наступного за розглянутим мінус одиницю. Наприклад, для вузла 5 останній індекс масиву у векторі
рівний 8=9-1, де номер 9 визначає початок масиву вузлів для вузла 6 (див. рис. 2.3).
У тих рідких випадках, коли той або інший розглядаємий вузол кінцево-елементної моделі не має зв'язку з вузлами, номера яких менше або більше номера даного вузла, тоді у відповідні гнізда векторів
або
заносяться нулі. Наприклад, вузол 4 (див. рис. 2.3) не має зв'язку з вузлами, номера яких менше 4, це враховується записом нуля в гніздо 4 вектора
(рис. 2.5). Вузол 3 не має зв'язку з вузлами, номера яких більше 3 – цій обставині відповідає запис нуля в гніздо 3 вектора
. Виключення становить запис в 1-е гніздо вектора
, куди заноситься номер останнього гнізда вектора
. Для всіх інших вузлів ці параметри визначаються тривіально. Для кожної кінцево-елементної моделі інформаційні вектори будуються один раз – відразу ж після того, як закінчена побудова масиву, що описує кінцеві елементи сітки.
Тепер розглянемо побудову сум добутків, що входять у рівняння (2.5). Розгляд почнемо з другої суми, а саме:
алгоритмічно виконується складніше. При фіксованому номері i
визначаються номери першого n і останнього m (n<m) гнізд вектора
, у якому зберігаються номери j
вузлів, пов'язаних з вузлом i
, і для яких j<i
. Цим самим задаються компоненти xi
вектора
. Тепер необхідно визначити елементи
aij
=
aji
(i≠j)
матриці A
, що зберігаються у векторі
. Ця процедура виконується в такий спосіб. Послідовно з гнізд із n по m у векторі
беруться номери вузлів j (j<i)
. По цих номерах j
з вектора
визначаються номери першого k і останнього l (k<l) гнізд вектора
. Для кожного вузла j
у діапазоні гнізд (l–k) вектора
перебуває номер i
, оскільки, по-перше, вузол i
пов'язаний з кожним вузлом j
, а, по-друге, i>j
. Переглядаючи гнізда вектора
і порівнюючи номера вузлів, що перебувають у цих гніздах, з номером i
, можна визначити номер гнізда, у якому для даного вузла j
і його діапазону номерів гнізд (l–k) перебуває номер вузла i
. По номеру гнізда вузла i
у векторі
можна знайти у векторі
елемент aij
=aji
(i≠j, j<i)
.
Розглянемо приклад. Побудуємо суму (2.8) для вузла 5 кінцево-елементної моделі, зображеної на рис. 2.3. Із цим вузлом, як ми вже відзначали вище, зв'язані вузли 4, 2 і 1. Номера цих вузлів зберігаються у векторі
у гніздах з 4 по 6, причому їх номера менше 5. Таким чином, компоненти вектора
відомі: х4
, х2
і x1
. Визначимо елементи матриці A
. Вузлу 4 у векторі
відповідають гнізда з 8 по 12. У цих гніздах зберігаються номери вузлів, пов'язаних з номером 4, причому всі ці номери більше 4. Один із цих номерів – 5. Цей номер перебуває в гнізді вектора
з номером 12. У гнізді з цим же номером 12 у векторі
зберігається елемент a45
=a54
(j=4<i=5)
матриці А
. Аналогічно відшукуються елементи a25
(j=2<i=5)
і a15
(j=1<i=5)
. Таким чином, одержуємо вираз для
:
Із представленого вище алгоритму видно, що обчислювальні витрати на побудову сум добутків
і
різні. Друга сума будується значно швидше. Це пов'язано з тим, що при побудові першої суми необхідно проводити додатковий порівняльний аналіз номерів вузлів, розміщених у векторі
. При розгляді процедур алгоритму передбачалося, що всі необхідні дані розташовуються в оперативній пам'яті комп'ютера. Довжина одномірних масивів, у яких розміщаються вектори
,
і
, залежить від структури розглядаємої кінцево-елементної моделі. У табл. 2.3 представлені дані по розмірності цих векторів для трьох типів кінцево-елементних моделей, побудованих у результаті розбивки одиничного куба на 20-вузлові квадратичні елементи. Куб розбивався в наступних пропорціях: 5:5:5 (КЕМ № 1), 10:10:10 (КЕМ № 2) і 20:20:20 (КЕМ № 3).
Таблиця 2.3 – Розмірність векторів КЕМ
Номер КЕМ |
Число вузлів |
Число елементів |
|
|
|
1 |
756 |
125 |
16079 |
16071 |
16079 |
2 |
4961 |
1000 |
121709 |
121691 |
121709 |
3 |
35721 |
8000 |
946619 |
946619 |
946619 |
Порівняння даного алгоритму з алгоритмами, розглянутими в роботах [17–19], показує наступне. Для реалізації запропонованого алгоритму потрібно більше оперативної пам'яті, оскільки для координування ненульових елементів використовуються не два масиви (вектори номерів рядків і стовпців), а чотири. Додаткові вектори
і
, розмірності яких не перевищують числа діагональних елементів, дозволяють локалізувати зони пошуку елементів-співмножників при побудові сум (2.6) і (2.8). За рахунок цього досягається економія часових витрат на розв’язок СЛАР.
На закінчення відзначу, що викладений вище алгоритм формування сум у рівнянні (2.5) може бути повністю застосований і для побудови інших аналогічних матричних добутків, необхідних для реалізації ітераційних процесів.
3 Оптимізація обчислень
Задача вирішення систем рівнянь виду
досить поширена. Вона виникає при описі декількох процесів, кожний з яких описується лінійним рівнянням. Сюди відноситься вивчення процесів, що протікають у лінійних електричних схемах, апроксимація систем диференціальних рівнянь, зокрема, які описують процеси теплопровідності.
Процедура розв’язку великої системи рівнянь на одному процесорі приводить до значних витрат часу. Побудова багатопроцесорної обчислювальної системи не завжди є доступним і оптимальним рішенням. З підвищенням надійності і швидкодії сучасних комп'ютерних мереж з'явилася можливість розв’язку подібних завдань у розподіленім обчислювальнім середовищі [21,22].
Оскільки компоненти середовища фізично рознесені в просторі, то кожний процес має локальну пам'ять і використання глобальних структур даних неможливо. Така схема взаємодії припускає обмін повідомленнями між процесами.
Пропонується будувати взаємодію між процесами по наступному принципу. Є деякий процес, який одержує дані від користувача. Він проводить попередню обробку цих даних. Потім процес формує повідомлення і розсилає їх необхідному числу допоміжних процесів, які виконують безпосередні обчислення. Відповідь повертається першому процесу у вигляді повідомлень від допоміжних процесів. Тут відбувається аналіз результатів обчислень і, при необхідності, процедура повторюється. Такий принцип взаємодії процесів зветься «керуючий-робітники». Перший процес є керуючим, допоміжні процеси – це робітники [23].
Існує кілька технологій побудови розподілених обчислювальних систем. Одна з перших, технологія RPC (Remote Procedure Call), заснована на віддалених викликах процедур. Процесу дозволяється викликати процедури, реалізовані на іншому, віддаленому комп'ютері. При цьому параметри процедури передаються на віддалений комп'ютер, там здійснюється виконання процедури і результати вертаються в місце виклику.
Технологія COM (Component Object Model) – це об'єктна модель компонентів. Вона застосовується для зв'язку між компонентами і програмами, а також при описі API і взаємодії різних мов і середовищ програмування. Одним з головних аспектів технології COM є реалізація клієнт-серверних додатків за допомогою інтерфейсів. В COM-програму входять COM-сервер, COM-клієнт і COM-інтерфейс.
COM-сервер – програма або бібліотека, яка надає деякі служби програмі-клієнтові. COM-сервер містить один або кілька COM-об'єктів. Кожний COM-об'єкт реалізує одну або кілька служб. Кожна служба описується інтерфейсом.
COM-клієнт – програма, яка використовує служби COM-сервера, запитуючи при цьому описані інтерфейси.
COM-інтерфейс – абстрактний клас, що описує властивості і методи COM-об'єкта. Клієнт через інтерфейс об'єкта може використовувати його методи і властивості як свої власні.
Існує також технологія CORBA (Common Object Request Broker Architecture) – узагальнена архітектура брокера об'єктних запитів. Дана технологія також заснована на ключовому понятті інтерфейсу. Об'єкт надає службу. Клієнт звертається до неї, використовуючи інтерфейс.
У технологіях COM і CORBA реалізована ідея розподілених об'єктів (distributed objects). При цьому інтерфейс служить універсальним засобом взаємодії. Розподілені об'єкти розміщаються на різних комп'ютерах. Коли процес викликає метод, реалізація інтерфейсу на машині з процесом перетворює виклик методу в повідомлення, передане об'єкту. Об'єкт виконує запит і повертає результати. Реалізація інтерфейсу перетворює відповідне повідомлення в повертаєме значення, яке передається викликаючому процесу.
Реалізація розподіленої системи, що розв’язує систему рівнянь виду
, можлива засобами технології COM. Перенос моделі «керуючий-робітники» на COM означає наступне. Керуючий процес повинен бути реалізований як COM-клієнт, робочі процеси, що виконують обчислення, реалізуються як COM-сервера послуги, що надають свої послуги через інтерфейси. COM-клієнт приймає вихідну систему від користувача, виконує перетворення до виду, зручного для ітерації, ділить систему на блоки, розсилає частини вихідних даних COM-серверам, використовуючи їх інтерфейси. Коли відповідь отримана, клієнт виконує збирання результатів і аналіз точності отриманого розв’язку. При досягненні заданої точності процес обчислень припиняється. Важливо, що в цьому випадку COM-клієнт реалізується як одиночна програма, що використовує послуги декількох віддалених COM-серверів.
У свою чергу, COM-сервера містять COM-об'єкти, які і реалізують функціональність, що забезпечує розв’язок систем виду
ітераційними методами.
Слід зазначити, що при такій реалізації розв’язок окремих блоків рівнянь вихідної системи здійснюється на різних хостах розподіленої системи. Це дозволяє робити обчислення паралельно і незалежно, що значно скорочує час розв’язку в порівнянні з обчисленнями на однопроцесорній системі.
4. Чисельні експерименти
Була розв’язана задача з кількістю невідомих, що перевищує чотири мільйона та проведене порівняння ефективності прямих і ітераційних методів. Конфігурація тестового комп'ютера наведена в табл. 4.1.
Таблиця 4.1 – Конфігурація тестового комп'ютера
Мікропроцесор |
Intel Core 2 Duo E7200 (Socket LGA775, 3172 МГц) |
Оперативна пам'ять |
4*(Kingston ValueRAM KVR667D2N5/1G) (4 ГБ) |
Системна плата |
ASUS P5K (Socket LGA775, Intel P35+ICH9) |
ОС |
Microsoft Windows Vista Ultimate x64 SP1 |
4.1 Пружне деформування тонкостінної просторової рами
Була розглянута задача про пружне деформування тонкостінної просторової рами (рис. 4.1) при наступних геометричних параметрах: L=H=2м, L1
=1,98м, L2
=0,01м. Пружні властивості матеріалу рами: E=106
Па, ν=0,3. Зовнішнє навантаження P=102
Н.
Кінцево-елементна модель усієї конструкції містила 1449604 вузла і 4374320 скінчених елементів у формі тетраедра. СЛАР розв’язувалася з використанням першої схеми. Час розв’язку склав близько двох годин.
4.2 Контактна взаємодія оболонкової конструкції і ложемента
Була вирішена задача про контактну взаємодію оболонкової конструкції і ложемента (рис. 4.2).
Дана задача часто виникає на практиці. При транспортуванні або зберіганні з горизонтальним розташуванням осі оболонкові конструкції встановлюються на кругові опори – ложементи. Взаємодія підкріплених оболонкових конструкцій і ложементів здійснюється через опорні шпангоути, довжина яких уздовж вісі оболонки порівняна з шириною ложементів і багато менше радіуса оболонки і величини зони контакту.
Дискретна модель ложемента (у тривимірній постановці) представлена на рис. 4.3.
При побудові даної КЕ-моделі було використано 6689вузлів і 15323 КЕ у формі тетраедра. Розмір нестиснутої матриці жорсткості для такої задачі становить (6689*3)2
*8=3221475912 байт, що приблизно дорівнює 3ГБ оперативної пам'яті.
Задача розв’язуваласядвома способами – методом Гауса і методом Ланцоша. Порівняння результатів розв’язку наведено в табл. 4.2.
Таблиця 4.2 – Результати розв’язку
Метод Гауса |
Метод Ланцоша |
Час розв’язку, (сек.) |
3137 |
2194 |
З табл. 4.2 легко бачити, що час розв’язку методом Ланцоша значно менше, ніж у випадку використання методу Гауса.
Висновки
У даній роботі розглянуті способи компактного зберігання матриці коефіцієнтів системи лінійних алгебраїчних рівнянь і методи розв’язку СЛАР. Розглянуті алгоритми компактного зберігання матриці жорсткості дозволяють значно скоротити об'єм необхідної пам'яті для зберігання матриці жорсткості.
Класичні методи зберігання, що враховують симетричну і стрічкову структуру матриць жорсткості, що виникають при застосуванні методу скінчених елементів, як правило, не застосовні при рішенні контактних задач, тому що при їхньому рішенні матриці жорсткості декількох тіл поєднуються в одну загальну матрицю, ширина стрічки якої може прагнути до розмірності системи.
Викладені у роботі методи компактного зберігання матриці коефіцієнтів СЛАР дозволяють на прикладі розв’язку контактних задач досягти значної економії оперативної пам'яті. Це дозволяє розв’язувати системи з мільйонами невідомих на сучасних персональних комп'ютерах.
Вважаю, що для оптимальної продуктивності при розв’язку СЛАР високих ступенів треба використовувати 24 ГБ оперативної пам'яті стандарту DDRIII 1333, мікропроцесор Intel Core i7-965 ExtremeEdition, системну плату ASUS RampageIIExtreme, дві відеокарти AMD Radeon HD 4870X2 та твердотільний накопичувачOCZSummit. Потужність такої конфігурації може скласти близько 5 терафлопс, що досить непогано. Для пришвидшення розрахунків вважаю доцільним розгін мікропроцесора, використання у програмі спеціалізованих наборів інструкцій процесорів по роботі з числами з плаваючою крапкою та використання у програмі ресурсів відеопроцесора відеокарти, за умови оптимізації програми під такі нововведення. Також значно пришвидшити розрахунки допоможе використання більшої кількості процесорів (або процесорних ядер).
Наведеніалгоритмирозв’язку СЛАР високих ступенів у задачах механіки можуть ефективно застосовуватися для аналізу тривимірних задач високої розмірності, що особливо актуально в задачах сучасного машинобудування.
Перелік посилань
1. Зенкевич О., Морган К. Конечные методы и аппроксимация. – М.: Мир, 1980. – 479 с.
2. Зенкевич О. Метод конечных элементов. – М.: Мир, 1975. – 217 с.
3. Стрэнг Г., Фикс Дж. Теория метода конечных элементов. – М.: Мир, 1977. – 346 с.
4. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. – М.: Наука, 1987. – 632 с.
5. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. – М.: Наука, 1984. – 176 с.
6. Бахвалов Н.С. Численные методы. – М.: Наука, 1975. – 153 с.
7. Годунов С.К. Решение систем линейных уравнений. – Новосибирск: Наука, 1980. – 81 с.
8. Gustavson F.G. Some basic techniques for solving sparse matrix algorithms. – New York: Plenum Press, 1972. – 215 p.
9. Джордж А., Лиу Дж. Численное решение больших разрежённых систем уравнений. – М.: Мир, 1984. – 333 с.
10. Rose D.J. A graph theoretic study of the numerical solution of sparse positive definite system of linear equations. – New York: Academic Press, 1972. – 473 p.
11. Мосаковский В.И., Гудрамович В.С., Макеев Е.М. Контактные задачи теории оболочек и стержней. – М.: Машиностроение, 1978. – 718 с.
12. Рвачев В.Л., Шевченко А.Н. Проблемно-ориентированные языки и системы для инженерных расчётов. – К.: Техніка, 1988. – 198 с.
13. Голуб Дж., Ван Лоун Ч. Матричные вычисления. – М.: Мир, 1999. – 548 с.
14. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. – М.: Наука, 1986. – 288 с.
15. Тьюарсон Р. Разрежённые матрицы. – М.: Мир, 1977. – 191 с.
16. Брамеллер А., Аллан Р., Хэмэм Я. Слабозаполненные матрицы. – М.: Энергия, 1979. – 192 с.
17. Эстербю О., Златев 3. Прямые методы для разрежённых матриц. – М.: Мир, 1987. – 120 с.
18. Писсанецки С. Технология разреженных матриц. – М.: Мир, 1988. – 410 с.
19. Сабоннадьер Ж.-К., Кулон Ж.-Л. Метод конечных элементов и САПР. – М.: Мир, 1989. – 192 с.
20. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. – М.: Наука, 1978. – 592 с.
21. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. – М.: Мир, 1991. – 347 с.
22. Валях Е. Последовательно-параллельные вычисления. – М.: Мир, 1985. – 216 с.
23. Таненбаум Э., Ван Стеен М. Распределённые системы. Принципы и парадигмы. – СПб.: Питер, 2003. – 481 с.
Додаток А
Вихідний текст програми
розв
’
язку СЛАР з розрідженою матрицею
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <fstream.h>
#include "matrix.h"
#include "smatrix.h"
#define BASE3D_4 4
#define BASE3D_8 8
#define BASE3D_10 10
const double Eps=1.0E-10;
DWORD CurrentType=BASE3D_10;
class RVector
{
private:
Vector<double> Buffer;
public:
RVector() {}
~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;
int Dim;
DWORD Size;
public:
TSMatrix(void) { Size=0; Dim=0; Array=NULL; Links=NULL; }
TSMatrix(Vector<DWORD> *, DWORD, int);
~TSMatrix(void) { if (Array) delete [] Array; }
Vector<double> &GetRight(void) { return Right; }
DWORD GetSize(void) { return Size; }
int 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]; }
};
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, int 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[INT(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);
NPoint=(j-NRow*(Array[i].Size()/Dim))/Dim;
NLink=j%Dim;
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));
}
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;
}
}
}
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; 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; 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; 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;
}
}
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;
if (argc<5 || argc>6) 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);
return;
}
if (argc==6)
{
op=argv[5];
if (strcmp(op, "/8") && strcmp(op, "/6"))
{
printf("Unknown options %s\n\n", op);
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;
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")) 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")) 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);
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;
}
}
}
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);
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]--;
return true;
}