Параллельные алгоритмы вычисления матриц жесткости конечных элементов высокого порядка многомерных задач матфизики Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 1, № 1, 1996

ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ ВЫЧИСЛЕНИЯ МАТРИЦ ЖЕСТКОСТИ КОНЕЧНЫХ ЭЛЕМЕНТОВ ВЫСОКОГО ПОРЯДКА МНОГОМЕРНЫХ ЗАДАЧ МАТФИЗИКИ*

А. Д. Матвеев Вычислительный центр СО РАН в г. Красноярске, Россия

Ю. В. Немировский Институт теоретической и прикладной механики СО РАН,

Новосибирск, Россия

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

Процедура построения матриц жесткости конечных элементов высокого порядка многомерных задач матфизики имеет матричную формулировку и сводится к вычислению интегралов от функций, построенных в местной системе координат или в Ь-координатах [1, 2]. Интегралы определяются численно, и эта процедура связана с выполнением большого объема вычислений.

Использование высокопроизводительной вычислительной техники — параллельных компьютеров при многочисленных исследованиях дает максимальную выгоду [3-5]. Поэтому в последнее время известные методы решения задач механики модифицируются с точки зрения выделения в них параллельных алгоритмов.

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

Краткая суть данного подхода состоит в том, что коэффициенты матрицы жесткости элементов определяются из условия минимизации квадратичного функционала краевой задачи, представленного в скалярной форме (в стандартной процедуре МКЭ функционал энергии имеет векторно-матричный вид), при определенном выборе структуры матрицы функций формы конечного элемента. Коэффициенты матрицы жесткости в данном случае выражаются в явном виде через группу интегралов.

Достоинства рассматриваемого подхода заключаются в том, что, во-первых, если переменные коэффициенты квазигармонического уравнения или модули упругости анизотропного неоднородного тела являются полиномами, то все интегралы в предложенных

*© А. Д. Матвеев, Ю. В. Немировский, 1996.

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

Вначале рассмотрим построение параллельных алгоритмов вычисления матриц жесткости конечных элементов для трехмерной задачи упругости. Не теряя общности суждений и для удобства изложения будем считать, что конечный элемент Уе второго порядка имеет форму прямоугольной призмы и расположен в декартовой системе координат ХУ Z (рис. 1), — местная система координат.

Рис. 1.

Пусть для анизотропного неоднородного конечного элемента Уе выполняются соотношения Гука и Коши [6]:

_ ди дх

дш дг

\Т

м = те},

дь ди дь

еу ду’ ^Ху ду + дх’

(1)

Ъ

ди дш дг дх'

1уг

дь д'ш дг + ду’

(2)

где {а} = {(ГхОуохТуХТхгТху}Т , {е} = {£х£уег^(уг 1хг 1хуУ — векторы напряжений и дефор-

маций; и,ь,ш — перемещения элемента Уе, [Б] — матрица модулей упругости С^ элемента Уе, С^ = С^г, С^ = С^(х,у,г), г,г] = 1, ..., 6; Т — транспонирование, функции Су удовлетворяют условиям эллиптичности [7].

Пусть функции С^ в области элемента Уe являются полиномами вида

е

г

т

Сг^ гкфк (г%к — действительные числа),

к=1

где _

фк = хакувкг_к (ак,вк, % — целые числа). (4)

Используя известные функции формы N1, ... , Ып элемента Уе, записанные в местной системе координат [1, 2], представим их в глобальной системе координат ХУZ, которые имеют вид

п

Na ^ ' таэ У] (га5 — действительные числа), (5)

3=1

Уз = хаэувэг1* (а],в],^3 — целые числа). (6)

Следуя МКЭ [1], аппроксимирующие функции перемещений и,ь,/ш элемента Уе определим соотношением

и } = Ш6}, Ш = ( . ) , (7)

ш ) У 5зп )

где структуру матрицы функций формы [N1 элемента Уе представим в форме [3]

[ы ]

N1 ...Ып 0 ... 0 0 ... 0

0 ... 0 N1 ...Ып 0 ... 0

0 ... 0 0 ... 0 N1 ...Ып

(8)

Согласно (7), (8), вектор {5} узловых неизвестных элемента Уе имеет следующую структуру: 5]^, ..., 5п есть узловые значения перемещений и, 5п+1, ... , 52п — перемещений V, 52п+1, ... , 53п — перемещений ш. Функционал энергии деформаций элемента Уе для представления (1) имеет вид

№е = J ^ 2 С11£х + С12ЄУ + С13ЄХ + С14^ух + С15^хг + С167ху^ +

Уе

+Єу(1 с^у + ^ + аиъ. + б'2^ + сЖ7ху) +

' 1

+Є.2 ^2 С33£, + С347у, + С35Іхх + С361ху^ +

+7у^ 2 С44іух + б45 7х, + С461х^ +

+7х^ ^ С55Іхг + б5б7ху^ + 1ху ^ 6667х^ ^У\ (9)

здесь Уе область элемента в системе координат ХУZ.

Из условия ~д5~ = 0, і = 1, ... , 3п, используя при этом в (9) представления (2), (7), (8),

т.е. учитывая структуру вектора {5}, получаем формулы для вычисления коэффициентов Кав (а, в = 1, ... , 3п) верхней треугольной части матрицы жесткости [К] элемента УЄ:

К,, = / {Сиаав + С\5(рав + р*,) + С,^ + ^) + Ст(и + Ы+

Уе

+ С55?«в + С66Ьав} ЛУ,

Ка+2п, в+2п — / {С33Яа/3 + С34(/ав + /ва) + С35(Рав + Рва) +

Уе

+ С45(^«в + ^ва) + С44Ьав + С55аав}dУ,

Ка+п, в+п / {С22Ьав + С24(/ав + /ва) + С26(^«в + ^в«) +

Уе

+ С46 (рав + рва) + С44?«в + С66аав }dУ,

здесь а = 1, ... , п, в = а, ..., п;

Ка+п,в+2п = / {С23/ав + С24Ьав + С25^0а +

Уе

+ С34?«в + С36Рав + С45Рва + С46^«в + С44/ва + С56аав "}^У-,

Ка,в+2п = / {С13рав + С14^«в + С15аав+

Уе

+ С35^«в + С36/ав + С45/ва + С46Ьав + С55рва + С56^ва}^У Ка,в+п = / {С12&ав + С14рав + С25Ува + С26Ьав + С45^ав+

Уе

+ С16аав + С56 рва + С46/ав + ^6^^^^

здесь а, в = 1,

п;

ду ду

дг дг

а = 1, ..., п, в = а, ..., п,

(10)

(11)

(12)

= дNа дNв = д^ д^ . = дNа дNв

рав с\ с\ , &ав о. гл , /ав г\ г\ ,

дх дг дх ду ду дг

а, в = 1, . . . , п.

Подставляя (5) в (12), (13) и учитывая (3), формулы (10), (11) представим в виде

Кав = С11ав + ... + С55ав + С66ав, Ка+п в+п = С02ав + ... + С44ав + ^

(13)

1а+п, в+п С22ав + ... + С44ав + С66ав

Ка+2п,в+2п = С%ав + ... + С44ав + С55ав,

а = 1, ..., п в = а, ..., п,

Ка+п,в+2п = С23ав + С24ав + ... + С56ав,

К-а в+2п = СГ3ав + С14ав + ... + СР

(14)

1а, в+2п = С13ав + С14ав Ка, в+п = С12ав + Ср4ав + ... + С<

55ва

Г

46ав,

где обозначено

С11а@ = С11аа@ЗУ\ • • • , С/6а@ = Сж/а/З■,

Уе

Уе

п п т

С11а/3 = Е Е Е Гаі Гві Гк1С11ізк, і=1 3=1 к=1

п п т

к 55ізк">

і=1 3=1 к=1

а = 1, • • •, п в = а,

, п,

п п т

С/3ав = Е Е Е Гаі Гв3 ГТСІзі3к, і=1 3=1 к=1

(16)

здесь а, в = 1,

п п т

с /

С46ав

„467^ /

121212ГРі Гаі ГкС 4біік,

і=1 3=1 к=1

п,

(17)

7=7/

С 11ізк = УХУ3хфк<ІУ, •••, С4бізк = УІуіфк<ІУ,

Уе

Уе

уХ = дУз/дх> уі = дУз/дУ> уі = дУз/дг^

(18)

Итак, в силу (14) — (17) коэффициенты матрицы жесткости элемента Уе в явном виде могут быть выражены через группу интегралов (18), для нахождения которых в силу (12), (13) достаточно вычислить 6 типов интегралов:

здесь і = 1,

31 = у уХуХфкзу, 3 = у

Уе Уе

, п; і = і, •••, п и

УІ УІ фк ІЇУ, З3

УІ УІ фк ,

Уе

34 = уХуУ фк ^, 35 = Угу УІ фк dv, 36 = уХуІ фк dv,

Уе

Уе

Уе

зДесь і,і = 1, • • • , п- В самом деле, например, С11ізк = С15ізк, Ср5ізк = СШ3к и т- д-

Пусть интегралы (18), т. е. 31, •••, 36, определяются численно. Заметим, что для элементов Уе формы тетраэдра интегралы 31 следует определять в Ь-координатах, для Уе (рис. 1) — в местной системе координат ^пС -

Рассмотрим один из вариантов параллельного алгоритма вычисления интегралов 3і, который выполняется на параллельном компьютере с шестью процессорами. Работа этого алгоритма показана на рис. 2. Основная программа расщепляется на шесть параллельных ветвей (вычислений): Ь1, •••, Ь6, которые реализуются соответственно на процессорах: Р1, • • • , Р6. Процессоры Рі работают одновременно и независимо друг от друга. Блок А

Рис. 2.

содержит операторы ввода и подготовки исходных данных, оператор распараллеливания вычислений (ветвей) и оператор распределения данных. Вначале выполняются операторы блока A, которые расположены в основной программе. Затем операторы ветвей Ll, ..., L6 соответственно вычисляют группы интегралов Jl, ..., J6 на процессорах Pl, ..., P6.

Блок B содержит операторы, которые вычисляют коэффициенты матрицы жесткости элемента по формулам (14) — (15) и расположены в программе. Эти операторы выполняются после завершения параллельных вычислений Li.

Для однородного элемента Ye (Cj = const, т. е. фk = const) общее число интегралов равно N = 3n(n + 1)/2 + 3n2. Следовательно, максимальное число параллельных ветвей равно N. Тогда в этом случае время выполнения параллельного алгоритма на параллельном компьютере, имеющем N процессоров, в t/N раз меньше, чем на моноЭВМ (t — время вычисления N интегралов на моноЭВМ). Например, для квадратичного элемента Ye(n = 20): N = 1830.

Отметим, что в блоках A,B, используя принцип потока данных [5], можно выделить ряд систем параллельных вычислений, имеющих структуру графов [4].

Расчеты на моноЭВМ показывают, что при численном интегрировании время построения матрицы жесткости однородного квадратичного элемента Ye формы куба по формулам (10) — (13) в 5 раз меньше времени вычисления [K] по матричной формуле МКЭ [1, 2], т.е. по формуле

Тг

[K]= [B]T [D][B]dY,

Ve

где [Б] — матрица, связывающая {е} и {5} при одинаковом выборе числа точек интегрирования в области V.

Для элемента Уе (см. рис. 1) интегралы (18) определяются аналитически. Используя

(4), (6) в (18) при ау, Ру> 1 найдем

C'

llijk

r\i r\i ( Vaxx Vaxx\ (\rbxx \rbxx \ ( r7cxx r7cxx\

aiaj (X2 - Xl ) (Y2 - Yl ) (Z2 - Zl )

axx b xx Cxx

С

1 вг 1з (Х%У* - Х^* ) (У^ - У^ ) (^Суг - ^уг )

46г] к

Ь

(19)

уг

■'У г

где

— аг + аз + ак — 1, Ьхх — рг + рз + в к + 15

схх — 1г + 1з + 7к + 1 • • • 5 ауг — аг + аз + ак, Ьуг — рг + Рз + в к5 суг — 7г + 7з + 7 к + 1

(Х\,У\, Z\), (Х25У25 И2) — координаты узлов А\, А2 элемента Уе в глобальной системе координат ХУZ (см. рис. 1).

При аз, вз, — 0 вычисления по формулам (14)—(17), (19) упрощаются. Например, если

аз — 0, то в силу (6) имеем д^з/дх — 0, и, следовательно, в силу (18) получим С 11гук — 0 и т. д.

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

Ква — Кав а — 1,..., 3п; в — а, •••, 3п

Пусть срединная поверхность недеформируемого тонкого прямоугольного анизотропного неоднородного конечного элемента Бе высокого порядка, находящегося в плоском напряженном состоянии, расположена в плоскости ХОУ. На рис. 3 показан прямоугольный элемент Бе второго порядка в общей системе координат ХУ и в местной системе

координат £ц.

Рис. 3.

Напряжения ах, ау, тху, деформации ех, еу, 7ху и перемещения и, V элемента Бе связаны соотношениями Гука и Коши [6]:

М — [ОД, (20)

а

ди

дь

ди дь

Єх дх’ Єу ду’ ^ху ду + дх’

(21)

где [Д — матрица упругости; {а}, (е) — векторы напряжений и деформаций, Сгз — модули упругости элемента Бе:

М

<3 х і ' С11 С12 С13

°у ? , Р] = С12 С22 С23

Тху ) С13 С23 Сзз

{є} =

'''/ху

Пусть функции С^ на Бе в системе координат ХУ имеют вид

Сз = ^2 ї'кфк (г%к — действительные числа),

к=1

где

(22)

(23)

фк — хакувк (ак,вк — целые числа).

Отметим, что гладкие функции Сгз, заданные для пластины, всегда можно на узловой сетке данного разбиения представить в виде (26), используя, например, линейные или квадратичные интерполяционные полиномы [2].

Потенциальная энергия деформаций Бе элемента с учетом (22), (23) имеет вид

We

Єх I 1 СііЄх + Сі2Єу + Сі^) + Єу(2С22Єу + Ох~,ху ) +

+ 2ІхуС33Іху( ЗБ.

(24)

Аппроксимирующие функции перемещений и,ь элемента Бе по МКЭ определяются соотношением

и V

[Я М,

(25)

где [Я] — матрица функций формы Бе размерности (2 х 2п), {£} — вектор узловых параметров МКЭ элемента Бе размерности 2п. За счет перестановки параметров вектора {$} матрицу [Я] представим в форме

[Я ]

N. ..яп о... о

0 ... 0 N1. ..Яп

(26)

где Ыа — функции формы элемента Бе(а =1, ... , п), построенные по алгоритмам МКЭ в общей системе координат ХУ

Е

3 = 1

Гаі уз,

где

Уз = х ^

гаі — действительные, аз— целые числа.

Используя представления (21) — (27), из условия

(27)

(28)

дWe

д6г

є

х

Є

у

т

п

определяем верхнюю треугольную часть матрицы жесткости [К] для элемента Бе, коэффициенты которой вычисляются по формулам

Кав — J {С11 Аав + С13(Дав + Дв«) + С33Рав} dS5 Ка+п,в+п — J {С22Рав + С23(Дз;в + Дв«) + С33Аав} (29)

а — 1, • • •, п, в — а, •••, п,

Ка,в+п — J {С12Дав + С13Аав + С23Рав + С33Два} dS■

а, в — 1, • • • , п,

где

. — дЫа дХв р — дЫа дЫв

Аав с\ с\ , рав о о ,

дх дх ду ду

а — 1, • • •, п, в — а, •••, п,

Д дХа дХв а л

Дав — ----^ а, в — 1, • • • , п•

дх ду

Дифференцируя (27) с учетом (28), получим

дХа А з дХа А 3

Ух — — аз ха?-1ув? ^ уу — — вз ха? ув?Ч (32)

^ — а- х"? -V? уз — ^

дх — а х У ^ Уу — ду

Подставляя (30) в (29) с учетом (31), (32), выражения (29) представим в форме

где обозначено

Кав — Сапав + С^ав + С^ + С^,

Ка+п, в+п — С22ав + С23ав + С23ва + Са3ав,

а — 1, • • •, п, в — а, •••, п,

Ка,в+п — С12ав + С13ав + С23ав + С33ав,

а, в = 1, • • • , п,

С11ав — У С11АавdS■ • • • ^ — ] С33Дав^

Яе Яе

п п т

гга

'а, г —С

СПав — Е Е Е ^ Гв? ^А;1 С°11гзк ^

г=1 з=1 к=1

n n m

С33ав = r«irej r^dCЗ3гjk, (34)

г=1 j=1 k=1

Спг3к — У УхУхФкdS■ •••, С33гзк — \ УхУ^уФкdS• (35)

Яе Яе

Итак, согласно (33), (34) коэффициенты матрицы жесткости [К] в явном виде зависят от группы интегралов (35). Для нахождения интегралов (35) достаточно вычислить интегралы типа

J угхузхфк J уу узуфк ^ У угхууфк

Яе Яе Яе

здесь I — 1, • • • , п, ] — 1, • • • , п^

Для прямоугольного Se (рис. 3) эти интегралы определяются аналитически. Подставляя (23), (32) в (35), имеем

Г\1 Г\1 ( ЛТаХХ \гахх \ (\гЬхх \ГЪхх\

Са — агаз (Х2 - Х1 ) (У2 - У1 )

С гзк а Ь ^

ахх Ьхх

-3 агвз(Х2аху - Хаху) (У>Ьху - УЬ1ху)

С гзк — а Ь , ( )

аху Ьху

где ахх аг + аз + ак 1 Ьхх вг + вз + вк + 1 • • • , аху аг + аз + акч Ьху вг + вз + вк;

(Х1,У1), (Х2,У2) — координаты узлов А1, А2 элемента Se (см. рис. 2).

Коэффициенты нижней треугольной части матрицы [К] определяются из условия ее симметрии.

При аг,вг — 0 имеем угх — 0, у у — 0 и, следовательно, в силу (35) получаем С^к — 0,

Сззк — 0 и т.д. Аналогично, как и для трехмерных элементов, определяем интегралы (35) с применением параллельных вычислений для высокоточных треугольных и четырехугольных элементов, используя при этом местные системы координат или Ь-координаты. В данном параллельном алгоритме вычисления интегралов максимальное число параллельных ветвей для однородного элемента Se равно Х — 3п(п + 1)/2.

Рассмотрим экономичную процедуру построения матрицы жесткости конечного элемента Уе высокого порядка формы прямоугольного параллелепипеда для стационарных краевых задач, которые описываются квазигармоническим уравнением [2] (записанным в декартовой системе координат ХУ Z) с переменными коэффициентами

s(K-ш) + i{K'yyI) + l(K-I) + Q = О в VeR <37)

с граничными условиями

у = у0 на S1 (у0 задана) (38)

и (или) на S2:

ду ду ду , ,

Kxx -qX1x + Куу ~Qyly + Kzz ~dz}z + q + h — У^ = 0,

где S = Si + S2 — полная граница области V; q,h = const; lx,ly,lz — направляющие косинусы вектора нормали к границе S2; функции Q = Q(x,y,z), Kxx = Kxx(x,y,z), Kyy = Kyy(x,y,z), Kzz = Kzz(x,y,z) удовлетворяют условиям эллиптичности [Т].

Пусть исходная область V представлена конечными элементами Уе высокого порядка (см. рис. 1) и пусть функции Кхх,Куу,КХХ в области конечных элементов V, являются полиномами вида

Кхх = £

к=1

г1Хик, Куу = Е ик,

к=1

= Е

гк ик,

к=1

где

ик = хакувкг1к;

(40)

(41)

здесь гхх, гк, гк — действительные числа, ак, в к, 7к — целые числа.

Аппроксимирующую функцию у элемента Уе определяем по МКЭ в глобальной системе координат ХУ Z, которая имеет вид [1, 2]

<*1

у = [N1... Яп]

(42)

где

Яа 'У ]га,фз (га, — действительные числа),

з=1

фз = ха6увіг16 (аз, вз,^з — целые числа).

Из условия

дWe

дбі

We

0, і = 1, ... , п, где

Уе

ду

2

у

2

Кхх1 дх) + КуЛ ду) + К*А дг' + 2^У

у

2

ЗБ+

+

ду + -к(у - у^)2

ЗБ;

(43)

(44)

(45)

здесь Se — поверхность элемента Ve.

Используя (42), (43) в (45) с учетом (40), получаем формулы для вычисления коэффициентов Кав(а, в — 1, • • • , п) верхней треугольной части матрицы жесткости [К] элемен-

та Уе

где

Кав = Аав + Вав + Сав + Рав,

а = 1, ..., п, в = а, ..., п,

п п т

Аав = ЕЕ Ег«, г в, гхх^цк,

і=1 з=1 к=1

(46)

п п т

Вав = Е Е Е г*і гв, гкУ Віз ^

і=1 з=1 к=1 п п т

Сав = ЕЕ Ег«.гв,гїСф,

і=1 з=1 к=1

т

т

т

п

1

2

Paf3 h^Y,r*ir>>jP ij •

(47)

i=1 j=1

A,

ijk

dj)i d^j x x

ukdS,

Bijk

Ve

Ve

дфі дф y y

ukdS,

C.

ijk

дфі дф_ dz dz

uk dS, P.

ij

Ve

ФіФ- dS.

(48)

Итак, в силу (46), (47) коэффициенты Кар матрицы жесткости элемента Уе выражаются в явном виде через группу интегралов (48). Отметим, что Аг]к = А]гк, В^к = В]гк, Сгзк = Сзгк, Рг] = Р ]%. Аналогично, как и для трехмерных элементов задачи упругости,

определяются интегралы (4s) с применением параллельных вычислений. Для трехмерной задачи (37) при Kxx = Kyy = Kzz = const в параллельном алгоритме вычисления интегралов (48) максимальное число параллельных ветвей равно N = 2n(n + І).

Для конечного элемента Ye формы прямоугольной призмы интегралы (48) определяются аналитически. Используя (41), (44), в (48) при ai,Pi,Yi ^ І имеем

A,

ijk

aiaj (Xax _ Xlx) (Y“x _ Ya) (Z2“x _ Z^)

ax bx Cx

S

Cjk = YiYj(X? _ Xf) (Y»‘ _ Yb) (Z"z _ Z?), (49)

ijk az bz Cz ’

где

ax = ai + aj + ak _ І, bx = Pi + Pj + вk + І,

Cx = Yi + Yj + Y k + І, • • •,

az = ai + aj + ak + І, bz = Pi + Pj + вk _ І,

Cz = Yi + Yj + Yk _ h

(X20 _ X10) (Y20 _ Y*0) (Z2" _ Z10); (5О)

ij ao bo Co

здесь ао = аг + а.^ + 1, Ьо = вг + в] + 1, со = 7г + Yj + 1, ), (Х2,^2,^2)

координаты узлов А1, А2 элемента Уе в глобальной системе координат ХУZ (см. рис. 1).

При а], в]У{] = 0 вычисления по формулам (46), (47), (49) упрощаются. Например, при а] = 0 в силу (44), (48) имеем Аг]к = 0, при вг = 0 имеем Вг]к = 0 и т. д. Нижняя треугольная часть матрицы [К] элемента Уе определяется из условия ее симметрии.

Расчеты, проведенные на моноЭВМ, показывают, что при численном интегрировании время построения матрицы жесткости элемента Уе второго порядка формы куба для задачи Дирихле (Кхх = Куу = Кгг = 1) по формулам (46)-(48) примерно в 1.8 раз меньше, чем по станрдартной процедуре МКЭ, и в 2.5 раз меньше для элемента Уе третьего порядка.

Результаты данной работы являются обобщением результатов [8, 9].

Список литературы

[1] Зенкевич О. Метод конечных элементов в технике. Мир, М., 1975.

[2] Сегерлинд Л. Применение метода конечных элементов. Мир, М., 1977.

[3] Серебряков В. А. Модель и язык для параллельных вычислений при решении научных задач. Журн. вычисл. матем. и матем. физ., 33, №7, 1993, 1083-1103.

[4] Воеводин В. В. Математические модели и методы в параллельных процессах. Наука, М., 1986.

[5] Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. Мир, М., 1991.

[6] Демидов С. П. Теория упругости. Высшая школа, М., 1979.

[7] Ректорис К. Вариационные методы в математической физике и технике. Мир, М., 1985.

[8] Матвеев А. Д., НЕмировский Ю. В. Энергетический метод определения матриц жесткости двумерных и трехмерных высокоточных элементов. Механика деформируемого твердого тела. ТГУ, Томск, 1988, 95-106.

[9] Матвеев А. Д. Аналитический метод построения жесткостей однородных элементов высокого порядка двух-, трехмерных задач упругости. ВЦ СО РАН, Красноярск, 1993.

Поступила в редакцию 12 мая 1996 г.