close

Вход

Забыли?

вход по аккаунту

?

Пояснительная записка (14)

код для вставкиСкачать
Министерство образования и науки Российской Федерации
Новосибирский государственный технический университет
Кафедра Прикладной математики
Пояснительная записка к курсовому проекту по курсу
Численных методов
Факультет: ПМИ
Группа: ПМ-12
Студент: Бибаев В.И.
Преподаватели: Персова М.Г., Задорожный А.Г.
Вариант: 89
Новосибирск
2014
Постановка задачи
МКЭ для трехмерной краевой задачи для уравнения эллиптического типа в декартовой системе координат. Базисные функции линейные на призмах с треугольным основанием (основание призмы в координатах (x, z)). Краевые условия всех типов. Коэффициент диффузии λ- кусочно - постоянная функция. Правая часть задана в виде производной некоторой скалярной функции с известными значениями в узлах. Матрицу СЛАУ генерировать в разреженном строчном формате. Для решения СЛАУ использовать МСГ или ЛОС с неполной факторизацией.
Эллиптическая задача для функции u определяется дифференциальным уравнением
-div(λ grad u)+ γu=f
В некоторой расчётной области Ω и границей S c заданными на ней краевыми условиями:
u|_(S_1 )=u_g,
λ ∂u/∂n |_(S_2 )=θ,
λ ∂u/∂n |_(S_3 )+β(u|_(S_3 )-u_β )=0
Расчётная область Ω - параллелепипед с параллельными координатным плоскостям гранями. Граница расчётной области S разбита на три непересекающиеся части, на которых заданы краевые условия различных типов. S=S_1+S_2+S_3
Теоретическая часть
Будем называть пространством H^m множество функций v, который вместе со своими производными до m-го порядка включительно интегрируемы с квадратов на Ω.
Исходное уравнение можно переписать в виде:
R(u)=0,
Где R(u)=-div(λ grad u)+ γu-f
Невязка для исходного уравнения
Потребуем, чтобы невязка дифференциального уравнения была ортогональна (в смысле скалярного произведения пространства L_2 (Ω)≡H^0) некоторому пространству пробных функций v, пространству Φ, т.е:
∫_Ω^ ▒〖(-div(λ grad u)+ γu-f)v dΩ=0〗 ∀v∈Φ
Далее воспользуемся формулой Грина:
∫_Ω^ ▒〖λ grad u grad v dΩ=-∫_Ω^ ▒〖div(λ grad u)v dΩ〗+∫_S^ ▒〖λ ∂u/∂n v dS〗〗
С помощью этой формулы преобразуем слагаемое:
∫_Ω^ ▒〖div(λ grad u)v dΩ〗
Получим:
Интеграл по границе, учитывая краевые условия можно преобразовать:
В качестве Φ выберем H_0^1 - пространство пробных функций v_0∈H^1, которые на границе S_1 удовлетворяют однородным первым краевым условиям. При этом будем считать, что u∈H_g^1, где H_g^1 -множество функций, имеющих суммируемые с квадратом первые производные и удовлетворяющие только первым краевым условиям на границе S_1. С учётом того, что v_0 |_(S_1 )=0:
Далее, заменим пространства H_g^1 и H_0^1 аппроксимирующими их конечномерными подпространствами〖 V〗_g^1 и V_0^1. Причём функции из этого пространства являются элементами одного и того же конечномерного пространства V^h, которое является линейным, натянутым
на базисные функции ψ_i, i=1,...,n
Базисные функции пространства V^h финитные кусочно - полиномиальные. И решение исходной задачи будем искать как линейную комбинацию базисных функций с некоторыми весами q_i.
Получим аппроксимацию последнего выражения в конечномерных подпространствах.
Для этого заменим функцию u∈H_g^1 аппроксимирующей её функцией u^h∈V_g^h, а функцию v_0∈H_0^1 функцией v_0^h∈V_0^h.
Поскольку любая функция v_0^h∈V_0^h может быть представлена в виде линейной комбинации
Подставив это выражение в предыдущее, и сократив на соответствующие веса, получим систему уравнений:
Поскольку функция u^h∈V_g^h может быть представлена в виде линейной комбинации базисных функций этого пространства
Причём n-n_0 компонент вектора весов q=(q_1,q_2,...,q_n ) могут быть фиксированы и определены из условия:
С учётом этого получим измененную СЛАУ для компонент вектора q, которые не были определены из первых краевых условий:
Таким образом получена СЛАУ вида Аq = b
Где компоненты матрицы А и вектора b определяются соотношениями:
Где δ_ij - символ Кронекера.
Теперь определим базисные функции и вид конечных элементов, а так же их нумерацию.
На рисунке изображена расчётная область с конечными элементами в виде прямых призм с треугольным основанием, так же отмечены глобальные номера узлов, эти же номера и соответствуют глобальным базисным функциям.
Нумерация конечных элементов следующая:
0: 0 1 4 9 10 13
1: 0 3 4 9 12 13
2: 1 2 5 10 11 14
3: 1 4 5 10 13 14
4: 3 4 7 12 13 16
5:3 6 7 12 15 16
И. т.д.
Теперь определим локальные номера узлов и базисных функций, а так же вид самих базисных функций
Для этого рассмотрим произвольный конечный элемент:
Определим на прямой треугольной призме Ω_m с вершинами
(x ̂_1, y ̂_1,z ̂_1 ), (x ̂_2, y ̂_1,z ̂_2 ), (x ̂_3, y ̂_1,z ̂_3 ), (x ̂_1, y ̂_2,z ̂_1 ), (x ̂_2, y ̂_2,z ̂_2 ), (x ̂_3, y ̂_2,z ̂_3 )
Шесть линейных базисных функция, равный единице в одной вершине призмы и нулю во всех остальных.
Для этого зададим на двумерном треугольнике, являющимся основанием призмы двумерные линейные базисные функции φ ̂_1 (x,z),φ ̂_2 (x,z), φ ̂_3 (x,z) , которые являются L - координатами оснований. Кроме этого зададим на отрезке [y ̂_1, y ̂_2 ] одномерные линейные базисные функции:
Y_1 (y)= (y ̂_2-y)/h_y , Y_2=(y- y ̂_1)/h_y где h_y=y ̂_2-y ̂_1
Тогда функции ψ_i (x,y,z) можно записать в виде
ψ_i (x,y,z)=φ ̂_(μ(i)) (x,z) Y_(v(i)) (y)
Где целочисленные функции μ(i), v(i) определяются формулами:
μ(i)=i mod 3, v(i)=[i/3]
Функции φ ̂_i (x,z) определяются следующим образом:
φ ̂_i (x,z)=α_0^i+α_1^i x+α_2^i z
А все коэффициенты определяются из следующих соотношений:
detD=(x ̂_2-x ̂_1 )(z ̂_3-z ̂_1 )-(x ̂_3-x ̂_2 )(z ̂_2-z ̂_1 )
[■(α_0^1&α_1^1&α_2^1@α_0^2&α_1^2&α_2^2@α_0^3&α_1^3&α_2^3 )]=1/detD [■(x ̂_2 z ̂_3-x ̂_3 z ̂_2&z ̂_2-z ̂_3&x ̂_3-x ̂_2@x ̂_3 z ̂_1-x ̂_1 z ̂_3&z ̂_3-z ̂_1&x ̂_1-x ̂_3@x ̂_1 z ̂_2-x ̂_2 z ̂_1&z ̂_1-z ̂_2&x ̂_2-x ̂_1 )]
Площадь треугольника в основании определяется так:
mesΩ_m^Δ=1/2 |detD|
Интеграл по треугольнику в основании удобно вычислять по формуле:
∫_(Ω_m^Δ)^ ▒〖(L_1 )^(v_1 ) (L_2 )^(v_2 ) (L_3 )^(v_3 ) dxdz〗=(v_1 !v_2 !v_3 !)/(v_1+v_2+v_3+2)! |detD|
Найдём локальные матрицы жесткости и массы для конечного элемента.
Будем считать, что параметры λ,γ на конечном элементе заменяются на свои осреднённые значения λ ̅, γ ̅, Тогда, учитывая полученные базисные функции:
G ̂_ij=∫_(Ω_m)^ ▒〖λ ̅((∂ψ ̂_i)/∂x (∂ψ ̂_j)/∂x+(∂ψ ̂_i)/∂y (∂ψ ̂_j)/∂y+(∂ψ ̂_i)/∂z (∂ψ ̂_j)/∂z)dΩ〗= =λ ̅(∫_(Ω_m^Δ)^ ▒((∂φ ̂_μ(i) )/∂x (∂φ ̂_μ(j) )/∂x+(∂φ ̂_μ(i) )/∂z (∂φ ̂_μ(j) )/∂z) dxdz∫_(y ̂_1)^(y ̂_2)▒〖Y_v(i) (y) Y_v(j) (y)dy〗+∫_(Ω_m^Δ)^ ▒〖φ ̂_μ(i) φ ̂_μ(i) dxdz∫_(y ̂_1)^(y ̂_2)▒〖(dY_v(i) (y))/dy (dY_v(j) (y))/dy〗〗 dy)
Обозначим матрицы жёсткости и массы соответствующих двумерных и одномерных линейных элементов следующим образом:
G_ij^xz=∫_(Ω_m^Δ)^ ▒((∂φ ̂_i)/∂x (∂φ ̂_j)/∂x+(∂φ ̂_i)/∂z (∂φ ̂_j)/∂z) dΩ=(α_1^i α_1^j+α_2^i α_2^j )mesΩ_m^Δ
M^xz=(mesΩ_m^Δ)/12 [■(2&1&1@1&2&1@1&1&2)]
G^y=1/h_y [■(1&-1@-1&1)]
M^y=h_y/6 [■(2&1@1&2)]
Таким образом, используя эти вспомогательные массивы, можно вычислить локальные матрицы жесткости и массы:
G ̂_ij=λ ̅(G_(μ(i)μ(j))^xz M_(v(i)v(j))^y+M_(μ(i)μ(j))^xz G_(v(i)v(j))^y ),
M ̂_ij=γ ̅M_(μ(i)μ(j))^xz M_v(i)v(j)^y.
В данном варианте правая часть задаётся в виде производной известной функции, значения которой известны во всех узлах сетки.
f=∂g/∂y
g=∑_(j=1)^N▒〖q_i^g ψ ̂_i 〗
Где N - количество узлов сетки.
При таких условиях изменится способ вычисления локального вектора правой части на конечном элементе (потребуется ещё одна вспомогательная матрица).
b ̂_i^Ω=∫_Ω^ ▒〖fψ ̂_i 〗 dΩ=∑_(j=1)^(N^l)▒〖q ̂_j^g ∫_(Ω_l)^ ▒〖(∂ψ ̂_j)/∂x ψ ̂_i 〗〗 dΩ
Где N^l - количество узлов на локальном элементе, а q ̂_i^g - известные веса разложения функции g по локальным базисным функциям конечного элемента〖 Ω〗_l.
Очевидно, что локальный вектор b ̂^Ω можно представить как произведение некоторой матрицы D ̂ вектора весов q ̂^g, где компоненты матрицы D ̂ вычисляются как:
D ̂_ij=∫_(Ω_l)^ ▒〖ψ ̂_i (∂ψ ̂_j)/∂y〗 dΩ=∫_(Ω_m^Δ)^ ▒〖φ ̂_μ(i) φ ̂_μ(i) dxdz∫_(y ̂_1)^(y ̂_2)▒〖Y_v(i) (y) (dY_v(j) (y))/dy〗〗 dy
Заметим, что значение первого интеграла это элемент уже известной матрицы〖 M〗^xz. А значение второго интеграла можно определить из следующей матрицы
M ̃^y=1/2 [■(-1&1@-1&1)]
Получим выражение для элементов матрицы D ̂:
D ̂_ij=M_(μ(i)μ(j))^xz M ̃_v(i)v(j)^y.
Краевые условия
Учёт краевых условий происходит уже после сборки глобальной матрицы. Сначала учитываются естественные условия (вторые и третьи), затем главные. Рассмотрим сначала учёт естественных краевых условий.
Третьи краевые условия:
λ ∂u/∂n |_(S_3 )+β(u|_(S_3 )-u_β )=0
Где β- константа
u_β- некоторая скалярная функция, заданная на границе S_3 своим интерполянтом.
Как было сказано ранее, от этого условия есть вклад в глобальную матрицу и в глобальный вектор правой части.
Добавка в глобальную матрицу:
A_ij=A_ij+∫_(S_3)^ ▒〖βψ_i ψ_j dS〗
Добавка в вектор правой части:
b_i=b_i+∫_(S_3)^ ▒〖βu_β ψ_j dS〗
Границу можно разбить на конечные элементы меньшей размерности (грани). Затем вычислить локальную матрицу, и добавить её к глобальной матрице.
Элемент локальной матрицы имеет вид для k - го конечного элемента с граничным условием:
A ̂_ij=∫_(S_3^k)^ ▒〖βψ ̂_i ψ ̂_j dS〗
Заметим, что граничные конечные элементы для данной задачи бывают двух типов:
Прямоугольник и треугольник. Рассмотрим два эти случая по отдельности, определим локальную нумерацию узлов и базисных функций.
Определим четыре линейные локальные базисные функции, равные единице в одном узле, и нулю в остальных.
T_1 (t)=(t_(p+1)-t)/h_t , T_2=(t-t_p)/h_t Y_1 (t)=(y_(s+1)-y)/h_y , Y_2=(y-y_s)/h_y Базисные функции на конечном элементе:
ψ ̂_1 (t,y)=T_1 (t)Y_1 (y)ψ ̂_2 (t,y)=T_2 (t)Y_1 (y)
ψ ̂_3 (t,y)=T_1 (t)Y_2 (y)ψ ̂_4 (t,y)=T_2 (t)Y_2 (y)
Вычислим локальную матрицу на этом элементе:
A ̂_11=∫_(S_3^k)^ ▒〖βψ ̂_1 ψ ̂_1 dS〗=β∫_(t_p)^(t_(p+1))▒〖T_1^2 dt∫_(y_s)^(y_(s+1))▒〖Y_1^2 dy〗〗=β (h_t h_y)/9
A ̂_12=∫_(S_3^k)^ ▒〖βψ ̂_1 ψ ̂_2 dS〗=β∫_(t_p)^(t_(p+1))▒〖T_1^ T_2 dt∫_(y_s)^(y_(s+1))▒〖Y_1^2 dy〗〗=β (h_t h_y)/18
И т.д могут быть вычислены все 16 элементов этой матрицы: выпишем её:
A ̂=β (h_t h_y)/36 [■(■(4&2@2&4)&■(2&1@1&2)@■(2&1@1&2)&■(4&2@2&4))]
Теперь найдём вклад этого конечного элемента в глобальный вектор правой части
Заметим, что u_β=∑_(i=1)^k▒〖u_β^i ψ ̂_i 〗
k - количество узлов на конечном элементе.
Т.е. вектор (b ) ̂ можно получить как произведение матрицы A ̂ на вектор весов u_β.
b ̂=A ̂u_β
u_β=(u_β^1,u_β^2,...,u_β^k)
После этого локальную матрицу и вектор правой можно добавлять в глобальные.
Теперь рассмотрим случай, когда грань является треугольником:
На этом элементе можно ввести L - координаты Зададим на двумерном треугольнике, двумерные линейные базисные функции φ ̂_1 (x,z),φ ̂_2 (x,z), φ ̂_3 (x,z). Теперь локальная матрица будет 3х3. Элементами этой матрицы являются:
A ̂_ij=∫_(S_3^k)^ ▒〖βψ ̂_i ψ ̂_j dS〗=∫_(Ω_m^Δ)^ ▒〖〖βφ ̂〗_μ(i) φ ̂_μ(i) dxdz〗=βM_ij^xz
Для локального вектора правой части:
b ̂=A ̂u_β
Вторые краевые условия учитываются аналогично, только вместо значений в узлах сетки функции u_β будут находится значения функции θ, так же β≔1.
Учёт главных краевых условий. Происходит уже после сборки глобальной матрицы и глобального вектора. Значения в узлах с первым краевым условием фиксируются:
На диагонали матрицы ставится единица, все элементы этой строки обнуляются, а в вектор правой части на соответствующую позицию ставится значение функции в этом узле. С помощью гауссова исключения ненулевых элементов из матрицы сохраняется свойство её симметричности.
Описание разработанных программ
Определим следующие структуры данных для хранения информации о расчётной области и сетки:
Массив вещественных чисел двойной точности net[N][3] - содержит декартовы координаты узлов сетки.
Целочисленный массив finiteElems[M][6] - содержит глобальные номера узлов принадлежащих к каждому конечному элементу в порядке локальной нумерации.
Два массива вещественных чисел двойной точности lymbda[M] и gama[M] для задания параметров λ,γ на каждом конечном элементе.
Для решения СЛАУ используется метод ЛОС с неполной LU(sq) факторизацией для симметричной матрицы.
Структуры данных для решения СЛАУ:
di - массив вещественных чисел размерности N - для хранения главной диагонали матрицы СЛАУ
gg - массив вещественных чисел, размерности M - для хранения элементов портрета матрицы. Размерность массива gg заранее неизвестна, но известно ограничение сверху - N(N-1)/2 элементов
ig - целочисленный индексный массив, содержит указатели на начала строк в массиве gg. Размерность N + 1
jg - целочисленный индексный массив, содержит номер столбца(строки) для вне диагональных элементов треугольника. Размерность M.
F - вектор правой части, содержит вещественные числа. Размерность N.
x - вектор решения СЛАУ.
Начальное приближение для ЛОС выбирается нулевым.
Программа состоит из следующих модулей:
(отдельная часть) Python-скрипт(netCreator.py), который по заданным границам области и шагам генерирует файлы, содержащие координаты узлов, список конечных элементов, в файл с первыми краевыми условиями на всей границе.
Модуль ввода/вывода(io): содержит функции ввода сетки, списка конечных элементов, коэффициентов дифференциального уравнения, параметров решения СЛАУ.
Модуль генерации портрета матрицы СЛАУ (portrait): содержит функции, реализующие списки и генерацию портрета матрицы.
Модуль создания СЛАУ(createSlau): вычисления локальных матриц и локальных векторов правой части, сборки глобальной матрицы и глобального вектора правой части, и процедуры учёта краевых условий.
Модуль решения СЛАУ(solveSlau): содержит функции для получения факторизации исходной матрицы и другие процедуры, необходимые для получения решения методом ЛОС с LU(sq) факторизацией и симметричной матрицей.
Модуль функций (functions): содержит функцию правой части в виде производной некоторой функции, так же содержит функцию верного решения задачи, для вычисления погрешности вычисления решения.
Описание тестирования программ
Протестируем программу на простых функциях, на которых должны получить точное решение. Будем проводить тесты на сетки с одним внутренним узлом, и первыми краевыми условиями. Результаты для сеток с меньшим шагом в отчёт включать не будем.
u=(x+z)y
f=(x+z)〖y/2〗^2
λ=1, γ=1
Краевые условия всюду первые
Тест проверяет правильность построения портрета, построения локальной матрицы массы (матрица жесткости на этом тесте нулевая), и учёт первых краевых условий.
Результат:
Получили решение без погрешности. Логичный результат т.к. базисные функции линейные, и дают второй порядок сходимости. При дроблении сетки так же строится точное решение.
u=x^2+y^2+z^2
f=(x+z)〖y/2〗^2
λ=1, γ=1
Краевые условия всюду первые
Тест проверяет правильность построение локальной матрицы жесткости, и порядок сходимости решения на полиноме второй степени.
Как и ожидалось, получили точное решение на полиномах второй степени. Проведенные исследования и выводы
Исследуем программу на неравномерной сетке, для этого изменим координаты узлов так, чтобы сетка задавалась следующим образом:
Точки на оси Ох: {0;0.5;2}
На оси Oz: {0;0.9;2}
На оси Oy: {0;1.3;2}
В качестве задачи выберем полином второй степени с первыми краевыми условиями
u=(x+z)y
f=(x+z)〖y/2〗^2
λ=1, γ=1
Проведем исследование на не полиномиальном решении для определения порядка аппроксимации. В качестве полиномиального решения выберем экспоненту. Для простоты ограничимся первыми краевыми условиями, параметры гама и лямбда константы, равные единице. Будем исследовать поведение решения на четырёх вложенных сетках. Погрешность аппроксимации определим по формуле:
√(∑_(i=1)^N▒(q_i^h-q_i^* )^2 )/√(∑_(i=1)^N▒(q_i^* )^2 )
Где q_i^h- компоненты вектора решения. q_i^*- точные значения функции в узлах. N-количество узлов в сетке.
u=e^(x+y+z)
f=-2e^(x+y+z)
λ=1, γ=1
Из этих данных можно сделать вывод, что порядок аппроксимации второй.
Тексты основных модулей программы
Заголовочные файлы:
io.h
#ifndef IO_H
#define IO_H
#include <conio.h>
#include <stdio.h>
#include <math.h>
#include "functions.h"
#include "portrait.h"
#include "createSlau.h"
#include "solveSlau.h"
extern int N, M;
extern int *finiteElems;
extern double *net, *lyambda, *gama;
void inputParameterSlau(int &maxiter, double &eps);
void inputNet();
void inputFinElems();
void outputVector(const char * filename, int size, double *vect);
#endif
portrait.h
#ifndef PORTRAIT_H
#define PORTRAIT_H
#include "io.h"
extern int *finiteElems;
struct list{
int elem;
list *next;
list(){
elem = -1;
next = NULL;
}
};
void generatePortrait(int *ig, int *jg);
#endif
functions.h
#ifndef FUNCTION_H
#define FUNCTION_H
#include "io.h"
double f(double x, double y, double z);
double solution(double x, double y, double z);
#endif
createSlau.h
#ifndef CREATESLAU_H
#define CREATESLAU_H
#include "io.h"
#include "functions.h"
extern int N, M;
extern int *finiteElems;
extern double *net, *lyambda, *gama;
extern int indexOfUnknows(int elem, int localNumber);
void fillMatrix(double *di, double *gg, int *ig, int *jg, double *F);
void firstCondition(double *di,double *gg, int *ig, int *jg, double *F);
void secondCondition(double *F);
void thirdCondition(double *di,double *gg, int *ig, int *jg, double *F);
#endif
solveSlau.h
#ifndef SOLVESLAU_H
#define SOLVESLAU_H
#include "io.h"
extern int N,M,K;
void getSolution(double* di, double* gg, double* F, double* x,int* ig,int* jg, int maxiter, double eps);
#endif
Реализация основных модулей:
main.cpp
#include "io.h"
int N, M, K;
int *finiteElems;
double *net, *lyambda, *gama;
void main(){
int maxiter;
int *ig, *jg;
double eps;
double *di, *gg, *F, *q;
inputParameterSlau(maxiter, eps);
inputNet();
inputFinElems();
F = new double[N];
q = new double[N];
di = new double[N];
ig = new int [N + 1];
jg = new int[N * (N - 1) / 2];
generatePortrait(ig, jg);
K = ig[N];
gg = new double[K];
for(int i = 0; i < K; i++)
gg[i] = 0;
for(int i = 0; i < N; i++)
di[i] = F[i] = q[i] = 0;
fillMatrix(di, gg, ig, jg, F);
secondCondition(F);
thirdCondition(di, gg, ig, jg, F);
firstCondition(di, gg, ig, jg, F);
getSolution(di, gg, F, q, ig, jg, maxiter, eps);
for(int i = 0; i < N; i++)
printf("%.15lf\n", abs(q[i] - solution(net[3 * i], net[3 * i + 1], net[3 * i + 2])));
outputVector("out.txt", N, q);
getch();
delete di, gg, ig, jg, F;
}
io.cpp
#include "io.h"
/*Ввод списка узлов сетки из файла "net.txt"
* Сначала указыватся целое число N - количество узлов сетки
* Затем в файле находятся N строк содержащих тройки вещественных чисел
* x,y,z координаты узлов
* Функция возвращается сслыку на массив, содержащий координаты узлов в пространстве
*/
void inputNet(){
FILE *f;
f = fopen("net.txt", "r");
fscanf(f, "%d", &N);
net = new double[N * 3];
for (int i = 0; i < 3 * N; i += 3)
fscanf(f, "%lf%lf%lf", &net[i], &net[i + 1], &net[i + 2]);
fclose(f);
}
/*Ввод списка конечных элементов из файла "finiteElements.txt"
* Сначала в файле указывается число M конечных элементов
* Затем M строк из 8 чисел, 6 первых целые - номера вершин призмы
* Два оставшихся вещественные, значения коэффиниентов Лямбда и Гама * На данном конечном элементе
* В массив l помещаются значения лямбда для каждого КЭ
* В массив g помещаются значения гама для каждого КЭ
*/
void inputFinElems(){
FILE *f;
double a, b;
f = fopen("finiteElements.txt", "r");
fscanf(f, "%d", &M);
finiteElems = new int[M * 6];
lyambda = new double[M];
gama = new double[M];
for(int i = 0; i < 6 * M; i += 6)
fscanf(f, "%d%d%d%d%d%d%lf%lf", &finiteElems[i], &finiteElems[i + 1], &finiteElems[i + 2], &finiteElems[i + 3],
&finiteElems[i + 4], &finiteElems[i + 5], &lyambda[i / 6], &gama[i / 6]);
fclose(f);
}
/*Ввод параметров решения СЛАУ из файла "SlauParams.txt"
* Параметры - максимальное число итераций, точность решения
*/
void inputParameterSlau(int &maxiter, double &eps){
FILE *f;
f = fopen("SlauParams.txt", "r");
fscanf(f, "%d%lf", &maxiter, &eps);
fclose(f);
}
/*Процедура вывод вектора вещественныз чисел размерности size
* В файл, с названием filename
*/
void outputVector(const char * filename, int size, double *vect){
int i;
double x, y, z, c, max = 0;
FILE *f;
f = fopen(filename, "w");
for(i = 0; i < size; i++){
x = net[3 * i];
y = net[3 * i + 1];
z = net[3 * i + 2];
fprintf(f, "%.15lf\t", vect[i] );
fprintf(f, "%.15lf\t", solution(x,y,z) );
c = abs(vect[i] - solution(x,y,z));
if (c > max) max = c;
fprintf(f, "%.15lf\n", c );
}
fprintf(f, "%.15lf", max);
fclose(f);
}
portrait.cpp
#include "portrait.h"
/*
* Операции со списками. Реализация списка: связный, упорядоченный по возрастанию, с фиктивным звеном
* Элементы этого списка уникальны
*/
/* Добавить элемент x с учётом порядка в список с указателем на голову *head*/
void append (int x, list *head){
list *added, *tmp;
added = new list;
added->elem = x;
for(tmp = head; tmp->next != NULL && tmp->next->elem <= x; tmp = tmp->next);
if(tmp->elem != x){
added->next = tmp-> next;
tmp->next = added;
}
}
/* * Извлечь первый элемент из списка
* При этом этот элемент из списка удаляется
* Результат - Значение функции
*/
int pop(list *head){
int rez;
list *popped;
popped = head->next;
if(popped == NULL)
return -1;
head->next = popped->next;
rez = popped->elem;
delete popped;
return rez;
}
/*
* Определить длину списка с указателем на голову head
* Результат - значение функции
*/
int len(list *head){
int i;
if(head -> next == NULL)
return 0;
for(i = 0; head -> next != NULL; head = head -> next, i++);
return i;
}
/*
* Напечатать содержимое списка
* Процедура для отладки
*/
void printlist(list *head){
list *tmp;
for(tmp = head->next; tmp != NULL; tmp = tmp ->next)
printf("%d ", tmp->elem);
}
/*
* Число узлов в конечном элементе с номером i
* Для данного типа сетки это число константа равная 6
*/
int numberOfUnknows(int i){
return 6;
}
/*
* Получить глобальный номер некоторой базисной * функции с локальным номером localNumber
* для конечного элемента с номером elem
* Результат - значение функции
*/
int indexOfUnknows(int elem, int localNumber){
return finiteElems[elem * 6 + localNumber];
}
/*
* Процедура построения портрета матрицы
* по заданной сетке и массиву конечных элементов
* Результат - заполненые массивы ig, jg
*/
void generatePortrait(int *ig, int *jg){\
int i, ielem, k, j, index1, index2;
list *l;
l = new list[N];
for(ielem = 0; ielem < M; ielem++){
for(i = 0; i < numberOfUnknows(ielem); i++){
k = indexOfUnknows(ielem, i);
for(j = i + 1; j < numberOfUnknows(ielem); j++){
index1 = k;
index2 = indexOfUnknows(ielem, j);
if(index2 > index1){
index1 = index2;
index2 = k;
}
append(index2, &l[index1]);
}
}
}
ig[0] = 0;
for(i = 0; i < N; i++){
printf(" element # %d contains %d: ", i, len(&l[i]));
printlist(&l[i]);
printf("\n");
}
for(i = 0; i < N; i++){
ig[i + 1] = ig[i];
k = pop(&l[i]);
while(k != -1){
jg[ig[i+1]] = k;
ig[i + 1]++;
k = pop(&l[i]);
}
}
}
functions.cpp
#include "functions.h"
double f(double x, double y, double z){
return - 2 * exp(x + y + z);//return (x+z) * y*y / 2.0;
}
double solution(double x, double y, double z){
return exp(x + y + z);
}
createSlau.cpp
#include "createSlau.h"
/*
Функции для построения элементов локальных матриц массы и жесткости
*/
int mu(int i){
return (i % 3);
}
int v(int i){
return (i / 3);
}
/*
* Процедура определия координат (x, y, z) для вершин конечного элемента
* Номера вершин конечного элемента (в порядке локальной нумерации) находятся в массиве L
* k - количество узлов в конечном элементе
*/
void getVertices(double &x1, double &x2, double &x3, double &y1, double &y2, double &z1, double &z2, double &z3, int *L, int k){
int i;
i = L[0];
x1 = net[3 * i];
y1 = net[3 * i + 1];
z1 = net[3 * i + 2];
i = L[1];
x2 = net[3 * i];
z2 = net[3 * i + 2];
i = L[2];
x3 = net[3 * i];
z3 = net[3 * i + 2];
i = L[k - 1];
y2 = net[3 * i + 1];
}
/*
* Получить коэффициенты для L - координат * для верхнего и нижнего основания призмы(по оси y)
* дополнительно вычисляется прощадь оснований призмы
*/
double getAlphaAndSquare(double *alpha, double x1, double x2, double x3, double z1, double z2, double z3){
double det;
det = (x2 - x1) * (z3 - z1) - (x3 - x1) * (z2 - z1);
alpha[0] = (z2 - z3) / det;
alpha[1] = (x3 - x2) / det;
alpha[2] = (z3 - z1) / det;
alpha[3] = (x1 - x3) / det;
alpha[4] = (z1 - z2) / det;
alpha[5] = (x2 - x1) / det;
return abs(det) / 2.0;
}
/*
* Процедура вычисления значений функции в узлах конечного элемента
* Результат поместить в массив p
* Размерности массивов L и p равно k
*/
void getFunctionValues(int k, int* L, double *p){
int i, j;
for(i = 0; i < k; i++){
j = 3 * L[i];
p[i] = f(net[j], net[j + 1], net[j + 2]);
}
}
/*Вспомогательные функции для вычисления локальной матрицы КЭ */
double Gxz(int i, int j, double *alpha){
return alpha[i * 2] * alpha[j * 2] + alpha[i * 2 + 1] * alpha[j * 2 + 1];
}
double My(int i, int j){
if(i == j) return 1.0 / 3;
return 1.0 / 6;
}
double My1(int i, int j){
if(j == 0) return -1.0 / 2;
return 1.0 / 2;
}
double Mxz(int i, int j){
if(i == j) return 1.0 / 6;
return 1.0 / 12;
}
double Gy(int i, int j){
if(i == j)
return 1;
return -1;
}
/*
* Процедура добавления локальной матрицы А, в матрицу в разреженном формате
* k - число узлов в конечном элементе
* A - локальная матрица
* L - глобальные номера узлов КЭ
*/
void addLocalMatrix(double *di, double *gg, int *ig, int *jg, int k, int *L, double *A){
int i, j, ibegin, iend, ind;
for(i = 0; i < k; i++)
di[L[i]] += A[i * k + i];
for(i = 0; i < k; i++){
ibegin = ig[L[i]];
for(j = 0; j < i; j++){
iend = ig[L[i] + 1] - 1;
while(jg[ibegin] != L[j]){
ind = (ibegin + iend) / 2;
if(jg[ind] < L[j])
ibegin = ind + 1;
else
iend = ind;
}
gg[ibegin++] += A[i * k + j];
}
}
}
/*
* Процедура добавления локального вектора в глоабльный
* k - число узлов в конечном элементе
* L - глобальные номера узлов КЭ
* localVector - локальный вектор правой части
* F - глобальный вектор правой части
*/
void addLocalVector(int k, int *L,double *localVector, double *F){
for(int i = 0; i < k; i++)
F[L[i]] += localVector[i];
}
/*Процедура построения локальной матрицы и локального вектораправой части
* для конечного элемента с номером elemNumber
* Результат в массиве A, локальный вектор в массиве b
* elemNumber - номер конечного элемента
* k - число узлов в конечном элементе
* A - Локальная матрица
* b - локальный вектор правой части
* L - глобальные номера узлов КЭ
* alpha - L - координаты оснований
*/
void createLocalMatrixVector(int elemNumber, int k, double *A, double *b, int *L, double *alpha, double *G, double *Mass, double *p){
int i, j;
double x1, x2, x3, y1, y2, z1, z2, z3;
double mes, hy, sum, l, g;
for(i = 0; i < k; i++)
L[i] = indexOfUnknows(elemNumber, i);
getVertices(x1, x2, x3, y1, y2, z1, z2, z3, L, k);
mes = getAlphaAndSquare(alpha, x1, x2, x3, z1, z2, z3);
hy = y2 - y1;
for (i = 0; i < k; i++)
for(j = 0; j < k; j++){
G[i * k + j] = mes * ( hy * Gxz(mu(i), mu(j), alpha) * My(v(i), v(j)) + Mxz(mu(i), mu(j)) * Gy(v(i), v(j)) / hy);
Mass[i * k + j] = mes * hy * Mxz(mu(i), mu(j)) * My(v(i), v(j));
}
l = lyambda[elemNumber];
g = gama[elemNumber];
for(i = 0; i < k * k; i++)
A[i] = l * G[i] + g * Mass[i];
for (i = 0; i < k; i++)
for(j = 0; j < k; j++)
Mass[i * k + j] = mes * Mxz(mu(i), mu(j)) * My1(v(i), v(j));
getFunctionValues(k, L, p);
for(i = 0; i < k; i++){
sum = 0;
for(j = 0; j < k; j++)
sum += p[j] * Mass[i * k + j];
b[i] = sum;
}
}
/*
* Процедура заполнения матрицы СЛАУ и вектора правой части
* массивы di, gg, ig, jg - массивы для разреженного хранения СЛАУ
* F - вектор правой части
* k - число узлов в конечном элементе
* L - глобальные номера узлов КЭ
*/
void fillMatrix(double *di, double *gg, int *ig, int *jg, double *F){
double *A, *b, *alpha, *G, *Mass, *p;
int i, k;
int *L;
k = 6;
alpha = new double[3 * 2];
A = new double[k * k];
L = new int[k];
b = new double[k];
G = new double[k * k];
Mass = new double[k * k];
p = new double[k];
for(i = 0; i < M; i++){
createLocalMatrixVector(i, k, A, b, L, alpha, G, Mass, p);
addLocalMatrix(di, gg, ig, jg, k, L, A);
addLocalVector(k, L, b, F);
}
delete A, b, alpha, L, G, Mass, p;
}
/*
* Процедура заполнения i - ой строки матрицы в разреженном формате нулями
* При этом сохраняется симметричность матрицы
* массивы di, gg, ig, jg - массивы для разреженного хранения СЛАУ
*/
void setNull(int i, double *gg, int *ig, int *jg, double *F){
int k, j, l;
for(k = ig[i]; k < ig[i + 1]; k++){
l = jg[k];
F[l] -= F[i] * gg[k];
gg[k] = 0;
}
for(k = i + 1; k < N; k++){
for(j = ig[k]; j < ig[k + 1] && jg[j] <= i; j++);
if(jg[--j] == i){
F[k] -= F[i] * gg[j];
gg[j] = 0;
}
}
}
/*
* Получить шаги по всем направлениям для прямоугольного конечного элемента
* L - глобальные номера узлов конечного элемента
*/
void getSteps(int *L, double &hx, double &hy, double &hz){
hx = abs(net[3 * L[0]] - net[3 * L[1]]);
hy = abs(net[3 * L[0] + 1] - net[3 * L[2] + 1]);
hz = abs(net[3 * L[0] + 2] - net[3 * L[1] + 2]);
}
double C(int i, int j){
if(i == j)
return 4.0;
if(i + j == 4)
return 1.0;
return 2.0;
}
/*
* Процедура получения вершин треугольника, глобальные номера узлов которого
* находятся в массиве L
*/
void getTriangleVertices(int *L, double &x1, double &x2, double &x3, double &z1, double &z2, double &z3){
x1 = net[3 * L[0]];
x2 = net[3 * L[1]];
x3 = net[3 * L[2]];
z1 = net[3 * L[0] + 2];
z2 = net[3 * L[1] + 2];
z3 = net[3 * L[2] + 2];
}
/*
* Процедура учёта первых краевых условий
* На диагональ ставится единица, внедиагональные элементы i - ой строки зануляются
* в i - й элемент вектора правой части ставится соответствующее знаение краевого условия
* Симметричность матрицы сохраняется
*/
void firstCondition(double *di,double *gg, int *ig, int *jg, double *F){
FILE *f;
int q, i, vertex;
double value;
f = fopen("first_condition.txt", "r");
fscanf(f, "%d", &q);
for(i = 0; i < q; i++){
fscanf(f, "%d%lf", &vertex, &value);
di[vertex] = 1;
F[vertex] = value;
setNull(vertex, gg, ig, jg, F);
}
}
/*
* Процедура учёта вторых краевых условий
* Добавляются вклады в глобальный вектор правой части от конечных элементов меньшей размерности
*/
void secondCondition(double *F){
FILE *f;
int i, j, k, q;
int *L;
double det, hx, hy, hz, ht, x1, x2, x3, z1, z2, z3;
double *teta, *b;
teta = new double[4];
b = new double[4];
L = new int[4];
f = fopen("second_condition.txt", "r");
fscanf(f, "%d", &q);
for(i = 0; i < q; i++){
fscanf(f, "%d%d%d%d%lf%lf%lf%lf", &L[0], &L[1], &L[2], &L[3], &teta[0], &teta[1], &teta[2], &teta[3]);
getSteps(L, hx, hy, hz);
ht = hx ? hx:hz;
for(j = 0; j < 4; j++){
b[j] = 0;
for(k = 0; k < 4; k++)
b[j] += hy * ht * teta[k] * C(j, k) / 36.0;
}
addLocalVector(4, L, b, F);
}
fscanf(f, "%d", &q);
for(i = 0; i < q; i++){
fscanf(f, "%d%d%d%lf%lf%lf", &L[0], &L[1], &L[2], &teta[0], &teta[1], &teta[2]);
getTriangleVertices(L, x1, x2, x3, z1, z2, z3);
det = abs((x2 - x1) * (z3 - z1) - (x3 - x1) * (z2 - z1));
for(j = 0; j < 3; j++){
b[j] = 0;
for(k = 0; k < 3; k++)
b[j] += teta[k] * ((k == j) ? 2.0 : 1.0) * det / 24.0;
}
addLocalVector(3, L, b, F);
}
}
/*
* Процедура учёта вторых краевых условий
* Добавляются вклады в глобальный вектор правой части и в глобальную матрицу
* от конечных элементов меньшей размерности
*/
void thirdCondition(double *di,double *gg, int *ig, int *jg, double *F){
FILE *f;
int i, j, k, q;
int *L;
double det, beta, hx, hy, hz, ht, x1, x2, x3, z1, z2, z3;
double *teta, *b, *A;
A = new double[16];
teta = new double[4];
b = new double[4];
L = new int[4];
f = fopen("third_condition.txt", "r");
fscanf(f, "%d", &q);
for(i = 0; i < q; i++){
fscanf(f, "%lf%d%d%d%d%lf%lf%lf%lf", &beta, &L[0], &L[1], &L[2], &L[3], &teta[0], &teta[1], &teta[2], &teta[3]);
getSteps(L, hx, hy, hz);
ht = hx ? hx:hz;
for(j = 0; j < 4; j++)
for(k = 0; k < 4; k++)
A[j * 4 + k] = beta * hx * ht * C(j, k) / 36.0;
addLocalMatrix(di, gg, ig, jg, 4, L, A);
for(j = 0; j < 4; j++){
b[j] = 0;
for(k = 0; k < 4; k++)
b[j] += A[4 * j + k] * teta[k];
}
addLocalVector(4, L, b, F);
}
fscanf(f, "%d", &q);
for(i = 0; i < q; i++){
fscanf(f, "%lf%d%d%d%lf%lf%lf", &beta, &L[0], &L[1], &L[2], &teta[0], &teta[1], &teta[2]);
getTriangleVertices(L, x1, x2, x3, z1, z2, z3);
det = abs((x2 - x1) * (z3 - z1) - (x3 - x1) * (z2 - z1));
for(j = 0; j < 3; j++)
for(k = 0; k < 3; k++)
A[j * 3 + k] = beta * ((k == j) ? 2.0 : 1.0) * det / 24.0;
addLocalMatrix(di, gg, ig, jg, 3, L, A);
for(j = 0; j < 3; j++){
b[j] = 0;
for(k = 0; k < 3; k++)
b[j] += A[3 * j + k] * teta[k];
}
addLocalVector(3, L, b, F);
}
}
solveSlau.cpp
#include "solveSlau.h"
/*
* Скалярное произведение двух векторов размерности N
*/
double scalarMultiple(double *a, double *b){
int i;
double rez = 0;
for(i = 0; i < N; i++)
rez += a[i] * b[i];
return rez;
}
/*
* Евклидова норма вектора размерности size
*/
double getNorm(double *vector, int size){
double rez = 0;
for (int i = 0; i < size; i++)
rez += vector[i] * vector[i];
return sqrt(rez);
}
/*
* Процедура умножения матрицы в разреженном строчном формате на вектор b.
* Результат помещает в вектор с
* Размерность матрицы и векторов полагается равный N
*/
void multiple(double *di, double *gg, int *ig, int *jg, double *b, double *c){
int i,k, m;
for (i = 0; i < N; i++){
c[i] = b[i] * di[i];
m = ig[i];
for(k = ig[i + 1] - 1;k >= m; k--){
c[i] += gg[k] * b[jg[k]];
c[jg[k]] += gg[k] * b[i];
}
}
}
/*
* Процедура получения LU(sq) разложения для симметричной матрицы
* G - внедиагональные элементы полученный матрицы
* D - диагональные элементы полученной матрицы
*/
void getLU(double *di, double *gg,int *ig,int *jg, double *G, double *D){
double sum;
int i, j, begin, end, first, last, ii, jj, k;
for(i = 0; i < N; i++){
D[i] = di[i];
begin = ig[i];
end = ig[i + 1];
for(j = begin; j < end; j++){
sum = 0;
k = jg[j];
first = ig[k];
last = ig[k + 1];
for(ii = begin, jj = first; ii < j && jj < last;)
if(jg[ii] == jg[jj])
sum += G[ii++] * G[jj++];
else
jg[ii] < jg[jj] ? ii++: jj++;
G[j] = (gg[j] - sum) / D[k];
D[i] -= G[j] * G[j];
}
D[i] = sqrt(D[i]);
}
}
/* * Процедура решения СЛАУ с верхнетреугольной матрицей в разреженном формате
* Вектор правой части - b
* Результат в x
*/
void solveUpperSystem(double *D, double *U, int *ig, int *jg, double *x, double *b){
int i, j, k;
for (i = 0; i < N; i++)
x[i] = b[i];
for(i = N - 1; i >= 0; i--){
x[i] /= D[i];
k = ig[i + 1];
for(j = ig[i]; j < k; j++)
x[jg[j]] -= x[i] * U[j];
}
}
/* * Процедура решения СЛАУ с нижнетреугольной матрицей в разреженном формате
* Вектор правой части - b
* Результат в x
*/
void solveLowerSystem(double *D, double *L, int *ig, int *jg, double *x, double *b){
for (int i = 0; i < N; i++){
x[i] = b[i];
for (int j = ig[i]; j < ig[i + 1]; j++)
x[i] -= x[jg[j]] * L[j];
x[i] /= D[i];
}
}
/*
* Процедура реализующая ЛОС для симметричной матрицы в разреженном формате хранения
* maxiter - максимально допустимое число итераций
* eps - требуемая относительная невязка
* Результат в x
*/
void getSolution(double *di, double *gg, double *F, double *x,int *ig,int *jg, int maxiter, double eps){
int i, j;
double error, alpha, beta, pp, normF;
double *G, *D, *r, *z, *p, *tmp1, *tmp2;
error = 0;
normF = getNorm(F, N);
G = new double[K];
D = new double[N];
r = new double[N];
z = new double[N];
p = new double[N];
tmp1 = new double[N];
tmp2 = new double[N];
getLU(di, gg, ig, jg, G, D);
multiple(di, gg, ig, jg, x, tmp1);
for (i = 0; i < N; i++)
tmp1[i] = F[i] - tmp1[i];
solveLowerSystem(D, G, ig, jg, r, tmp1);
solveUpperSystem(D, G, ig, jg, z, r);
multiple(di, gg, ig, jg, z, tmp1);
solveLowerSystem(D, G, ig, jg, p, tmp1);
for(i = 0; i < N; i++)
error += r[i] * r[i];
error = sqrt(error) / normF;
for(i = 1; (error > eps) && (i < maxiter); i++){
error = 0;
pp = scalarMultiple(p, p);
alpha = scalarMultiple(p, r) / pp;
for(j = 0; j < N; j++){
x[j] += alpha * z[j];
r[j] -= alpha * p[j];
error += r[j] * r[j];
}
solveUpperSystem(D, G, ig, jg, tmp1, r);
multiple(di, gg, ig, jg, tmp1, tmp2);
solveLowerSystem(D, G, ig, jg, tmp1, tmp2);
beta = - scalarMultiple(tmp1, p) / pp;
solveUpperSystem(D, G, ig, jg, tmp2, r);
for(j = 0; j < N; j++){
z[j] = tmp2[j] + beta * z[j];
p[j] = tmp1[j] + beta * p[j];
}
error = sqrt(error) / normF;
printf("iteration #%derror = %.15lf\n", i, error);
}
delete G, D, r, z, p, tmp1, tmp2;
}
Python-скрипт для генерации сетки по заданным границами длине шага
#Создание равномерной сетки из треугольных призм для некоторой геометрической области
#Границы области и шаги по всем направлениям находятся в файле field.txt
#Выходные файлы:
#net.txt - упорядоченный список узлов сетки
#finiteElements.txt - упорядоченный списов конечных элементов
#first_condition.txt - первые краевые условия во всех граничных точках
import math
stepDevider = 1
def f(x, y, z):
return math.exp(x + y + z)
file = open('field.txt', 'r')
buf = file.readline()
buf = buf.split()
Xmin = float(buf[0])
Xmax = float(buf[1])
Xh = float(buf[2]) / stepDevider
buf = file.readline()
buf = buf.split()
Ymin = float(buf[0])
Ymax = float(buf[1])
Yh = float(buf[2]) / stepDevider
buf = file.readline()
buf = buf.split()
Zmin = float(buf[0])
Zmax = float(buf[1])
Zh = float(buf[2]) / stepDevider
file.close()
fileNet = open('net.txt', 'w')
Xnet = []
Ynet = []
Znet = []
N = int((Xmax - Xmin) / Xh) + 1
for i in range(N):
Xnet.append(Xmin + i * Xh)
N = int((Ymax - Ymin) / Yh) + 1
for i in range(N):
Ynet.append(Ymin + i * Yh)
N = int((Zmax - Zmin) / Zh) + 1
for i in range(N):
Znet.append(Zmin + i * Zh)
if Xnet[-1] < Xmax:
Xnet.append(Xmax)
if Ynet[-1] < Ymax:
Ynet.append(Xmax)
if Znet[-1] < Ymax:
Znet.append(Zmax)
Xstep = len(Xnet)
Ystep = len(Ynet)
Zstep = len(Znet)
fileNet.write(str(Xstep * Ystep * Zstep))
for y in Ynet:
for z in Znet:
for x in Xnet:
fileNet.write('\n' + str(x) + ' ' + str(y) + ' ' + str(z))
fileNet.close()
fileFinElements = open('finiteElements.txt', 'w')
fileFinElements.write(str(2 * ((Xstep - 1) * (Ystep - 1) * (Zstep - 1))))
delta = Xstep * Zstep
lower = []
upper = []
basic1 = [0, 1, 1 + Xstep]
basic2 = [0, Xstep, 1 + Xstep]
lower.append(basic1)
lower.append(basic2)
for i in range(Xstep - 2):
basic1 = [basic1[0] + 1, basic1[1] + 1, basic1[2] + 1]
basic2 = [basic2[0] + 1, basic2[1] + 1, basic2[2] + 1]
lower.append(basic1)
lower.append(basic2)
basic1 = []
for x in lower:
basic1.append(x)
for i in range(Zstep - 2):
basic2 = []
for x in basic1:
basic2.append([x[0] + Xstep, x[1] + Xstep, x[2] + Xstep])
for x in basic2:
lower.append(x)
basic1 = basic2
for x in lower:
upper.append([x[0] + delta, x[1] + delta, x[2] + delta])
for i in range(Ystep - 1):
for down, up in zip(upper, lower):
fileFinElements.write('\n' + str(up[0]) + ' ' + str(up[1]) + ' ' + str(up[2]) + ' ' + str(down[0]) + ' '
+ str(down[1]) + ' ' + str(down[2]) + ' 1 1')
lower = upper
upper = []
for x in lower:
upper.append([x[0] + delta, x[1] + delta, x[2] + delta])
fileFinElements.close()
fileCondition = open('first_condition.txt', 'w')
fileNet = open('net.txt', 'r')
N = int(fileNet.readline())
conditions = {}
for i in range(N):
buf = fileNet.readline().split()
x = float(buf[0])
y = float(buf[1])
z = float(buf[2])
if x == Xmin or x == Xmax or y == Ymin or y == Ymax or z == Zmin or z == Zmax:
conditions.update({i: f(x,y,z)})
fileCondition.write(str(len(conditions)))
for key, value in conditions.items():
fileCondition.write('\n' + str(key) + ' ' + str(value))
fileCondition.close()
Документ
Категория
Рефераты
Просмотров
38
Размер файла
407 Кб
Теги
пояснительная, записка
1/--страниц
Пожаловаться на содержимое документа