Современные численные методы решения граничных задач. Основы метода конечных элементов, виды конечных элементов и функций формы

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ

УЧРЕЖДЕНИЕ ОБРАЗОВАНИЯ

ГОМЕЛЬСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ ИМЕНИ П. О. СУХОГО

Факультет автоматизированных и информационных систем

Кафедра «Информационные технологии»

дисциплина «Компьютерные системы конечно-элементных расчётов»

ОТЧЕТ ПО ЛАБОРАТОРНОЙ РАБОТЕ № 2

«Современные численные методы решения граничных задач»

Выполнила: студентка гр. ИТ-31

Принял: преподаватель

Гомель 2012

Лабораторная работа №2

Современные численные методы решения граничных задач

Цель: изучить основы метода конечных элементов, виды конечных элементов и функций формы.

Задание.

Согласно варианта задания (таблица 2) необходимо построить математическую модель указанной физической системы и провести её исследование следующими численными методами: методом конечных разностей, методом граничных элементов и методом конечных элементов. Для решения данной задачи необходимо разработать соответствующее программное обеспечение, которое должно удовлетворять следующим требованиям:

- обеспечить ввод исходных данных с помощью GUI;

- отобразить в виде графиков (двумерного и трёхмерного) результаты решения;

- все результаты решения сохранять как в тестовые файлы (для претендующих на оценки 4-5), так и в файлы специальных форматов (для всех остальных).

Текст программы МКЭ:

using System;

using System.Collections.Generic;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Linq;

using System.Text;

using System.Windows.Forms;

using System.IO;

using System.Threading;

namespace lab2_kskr

{

public partial class Form1 : Form

{

public int xr, yr;

public double AB, CD, G, ny, Sa, stkv;

public double dx, dy, t, h;

public double[] x, y, q, xx,F, alpha;

public double xn, xk, yn, yk;

public double[,] A, B, KG, Q, k, k1, k2, ObrA, E, U, V;

public Form1()

{

InitializeComponent();

ny = 0.33f;

h = 0.01f;

}

public void INIT_k(int nomer)

{

for (int i = 0; i < 8; i++)

{

for (int j = 0; j < 8; j++)

{

A[i, j] = 0;

ObrA[i, j] = 0;

if (i < 3)

{

B[i, j] = 0;

Q[i, j] = 0;

if ((i == 2 || j == 2) && j < 3)

E[i, j] = 0;

}

}

}

for (int ii = 0; ii < 4; ii++)

{

U[ii, 0] = 1;

U[ii, 1] = x[ii];

U[ii, 2] = y[ii];

U[ii, 3] = x[ii] * y[ii];

}

if (nomer == xr + 1 || nomer == xr * 2 - 2 || nomer == xr * yr - xr - 2 || nomer == xr * yr - 2 * xr + 1)

{               

G = 2 * (double)Math.Pow(10, 0);

}

else

G = 2 * (double)Math.Pow(10, 11);

E[0, 0] = E[1, 1] = G / (1 - ny * ny);

E[0, 1] = E[1, 0] = G * ny / (1 - ny * ny);

E[2, 2] = 0.5f * G * (1 - ny) / (1 - ny * ny);

int ot = 0, to = 4;

int kk = 0;

for (int i = 0; i < 8; i++)

{

for (int j = ot, u = 0; j < to; j++, u++)

{

A[i, j] = U[kk, u];

}

if ((i + 4) % 2 == 1)

{

ot = 0; to = 4;

kk++;

}

else

{

ot = 4; to = 8;

}

}

M1();

}

public void M1()

{

int n = 8;

double[,] m = new double[8, 8];

for (int i = 0; i < n; i++)

{

for (int u = 0; u < 8; u++)

for (int j = 0; j < 8; j++)

m[u, j] = A[u, j];

for (int j = 0; j < n; j++)

if (j == i)

q[j] = 1;

else

q[j] = 0;

MethodGausa(m, q, xx);

for (int j = 0; j < n; j++)

{

ObrA[j, i] = xx[j];

}

}

}

public void DrawingElement()

{

Graphics gr = panel1.CreateGraphics();

int size_x = panel1.Width;

int size_y = panel1.Height;

Pen p1 = Pens.Red;

float shagx = size_x*1.0f / xr;

float shagy = size_y*1.0f / yr;

gr.FillRectangle(Brushes.Aqua, shagx, shagy, shagx, shagy);

gr.FillRectangle(Brushes.Aqua, shagx, size_y - 2 * shagy, shagx, shagy);

gr.FillRectangle(Brushes.Aqua, size_x - 2 * shagx, shagy, shagx, shagy);

gr.FillRectangle(Brushes.Aqua, size_x - 2 * shagx, size_y - 2 * shagy, shagx, shagy);

for (float yyn = 0; yyn < size_y; yyn += shagy)

gr.DrawLine(p1,0,yyn,size_x,yyn);

for (float xxn = 0; xxn < size_x; xxn += shagx)

gr.DrawLine(p1, xxn, 0, xxn, size_y-1);

gr.DrawLine(p1, 0, size_y-1, size_x, size_y-1);

gr.DrawLine(p1, size_x-1, 0, size_x-1, size_y - 1);

}

public void DrawingElement2(double[] matr)

{

Graphics gr2 = panel2.CreateGraphics();

int size_x = panel2.Width;

int size_y = panel2.Height;

Pen p1 = Pens.Red;

int i = 1;

float sx = 0;

float kfr = 30000;

float shagx = size_x * 1.0f / xr, shagx1 = size_x-1;

float shagy = size_y * 1.0f / yr, shagy1 = size_y - 1;

for (float yyn = size_y+1; yyn >= 0; yyn -= shagy)

{

sx = 0;

shagx1 = shagx;

for (float xxn = 0; xxn < size_x; xxn += shagx)

{

gr2.DrawLine(p1, sx, yyn - (float)matr[i]*kfr, shagx1-1, yyn - (float)matr[i+2]*kfr);

sx = shagx1;

shagx1 += shagx;

i += 2;                   

}

i += 2;               

}

i = 0;

gr2.DrawLine(p1, 0, size_y-1, size_x-1, size_y - 1);

gr2.DrawLine(p1, size_x - 1, 0, size_x - 1, size_y - 1);

for (float xxn = 0; xxn < size_x; xxn += shagx)

gr2.DrawLine(p1, xxn, 0, xxn, size_y - 1);

}

public void DrawingElement3(double[] matr)

{

Graphics gr2 = panel2.CreateGraphics();

int size_x = panel2.Width;

int size_y = panel2.Height;

Pen p1 = Pens.Blue;

int i = 1;

float sx = 0;

float kfr = 29000;

float shagx = size_x * 1.0f / xr, shagx1 = size_x - 1;

float shagy = size_y * 1.0f / yr, shagy1 = size_y - 1;

for (float yyn = size_y + 1; yyn >= 0; yyn -= shagy)

{

sx = 0;

shagx1 = shagx;

for (float xxn = 0; xxn < size_x; xxn += shagx)

{

gr2.DrawLine(p1, sx, yyn - (float)matr[i] * kfr, shagx1 - 1, yyn - (float)matr[i + 2] * kfr);

sx = shagx1;

shagx1 += shagx;

i += 2;

}

kfr -= 1000;

i += 2;

}

i = 0;

gr2.DrawLine(p1, 0, size_y - 1, size_x - 1, size_y - 1);

gr2.DrawLine(p1, size_x - 1, 0, size_x - 1, size_y - 1);

for (float xxn = 0; xxn < size_x; xxn += shagx)

gr2.DrawLine(p1, xxn, 0, xxn, size_y - 1);

}       

private void button2_Click(object sender, EventArgs e)

{

AB = 1; CD = 1; xr = 6; yr = 5; xn = 0; yn = 0; xk = AB + xn; yk = CD + yn;

dx = (xk - xn) / xr; dy = (yk - yn) / yr;

t = 0.1f; stkv = AB / 10;

Sa = dx * dy * 0.01f;

DrawingElement();

alpha = new double[2 * (xr + 1) * (yr + 1)];

k = new double[8, 8]; k1 = new double[8, 8];

F = new double[2 * (xr + 1) * (yr + 1)];

KG = new double[2 * (xr + 1) * (yr + 1), 2 * (xr + 1) * (yr + 1)];

k2 = new double[8, 8];

A = new double[8, 8]; B = new double[3, 8]; x = new double[4]; y = new double[4];

ObrA = new double[8, 8]; Q = new double[3, 8]; E = new double[3, 3];

q = new double[8]; xx = new double[8]; U = new double[4, 4];

AA = new double[2 * (xr + 1) * (yr + 1), 2 * (xr + 1) * (yr + 1)];

BB = new double[2 * (xr + 1) * (yr + 1)];           

for (int i = 0; i < KG.GetLength(0); i++)

{

F[i] = 0;

BB[i] = 0;

for (int j = 0; j < KG.GetLength(1); j++)

{

KG[i, j] = 0;

AA[i, j] = 0;

}

}

button2.Enabled = false;

//  MessageBox.Show("Инициализация выполнена!");

}

public double[,] LK(double[] x1, double[] y1)

{

double[,] temp1 = new double[8, 8];

double[,] temp = new double[8, 8];

Q_Q(x1, y1);

for (int ii = 0; ii < 8; ii++)

for (int i = 0; i < 8; i++)

{

double sum = 0;

for (int j = 0; j < 8; j++)

{

sum += k2[j, i] * ObrA[j, ii];

}

temp1[ii, i] = sum;

}

//FileWrite(temp1);

for (int ii = 0; ii < 8; ii++)

for (int i = 0; i < 8; i++)

{

double sum = 0;

for (int j = 0; j < 8; j++)

{

sum += temp1[ii, j] * ObrA[j, i];

}

temp[i, ii] = sum;

}

// FileWrite(temp);

return temp;

}

public void Q_Q(double[] x1, double[] y1)

{

for (int i = 0; i < 8; i++)

for (int j = 0; j < 8; j++)

k2[i, j] = 0;

double koef = G / (ny * ny - 1);

k2[1, 1] = k2[6, 6] = -koef * Sa;

k2[1, 6] = k2[6, 1] = ny * k2[1, 1];

k2[2, 5] = k2[5, 2] = k2[2, 2] = k2[5, 5] = koef * (ny - 1) * Sa / 2;

k2[2, 3] = k2[3, 2] = k2[3, 5] = k2[5, 3] = (x[1] * x[1] - x[0] * x[0]) * h * koef * dy * (ny - 1) / 4;

k2[1, 3] = k2[3, 1] = -(y[2] * y[2] - y[0] * y[0]) * koef * h * dx / 2;

k2[3, 3] = ((x[1] * x[1] * x[1] - x[0] * x[0] * x[0]) * koef * h * dy * (ny - 1) / 6) - ((y[2] * y[2] * y[2] - y[0] * y[0] * y[0]) * koef * h * dx / 3);

k2[1, 7] = k2[7, 1] = -koef * h * dy * (x[1] * x[1] - x[0] * x[0]) * ny / 2;

k2[2, 7] = k2[7, 2] = (y[2] * y[2] - y[1] * y[1]) * koef * dx * h * (ny - 1) / 4;

k2[3, 6] = k2[6, 3] = k2[1, 3] * ny;

k2[3, 7] = k2[7, 3] = ((x[1] * x[1] - x[0] * x[0]) * (y[2] * y[2] - y[0] * y[0]) * h * koef * (ny - 1) / 8) - ((x[1] * x[1] - x[0] * x[0]) * h * (y[2] * y[2] - y[0] * y[0]) * koef * ny / 4);

k2[5, 7] = k2[7, 5] = k2[2, 7];

k2[6, 7] = k2[7, 6] = k2[1, 7] / ny;

k2[7, 7] = ((y[2] * y[2] * y[2] - y[0] * y[0] * y[0]) * h * koef * dx * (ny - 1) / 6) - ((x[1] * x[1] * x[1] - x[0] * x[0] * x[0]) * koef * h * dy / 3);

}

private void button1_Click(object sender, EventArgs e)

{

double[,] uzlK = new double[(xr + 1) * (yr + 1) * 2, 8];//номера

int[,] uzl = new int[(yr + 1) * (xr + 1) * 2, 8];//координаты вершин

int koord = -1; int L = 0;

for (double yyn = yn; yyn <= (yk-dy); yyn += dy)

{

for (double xxn = xn; xxn <= (xk-dx); xxn += dx)

{

x[0] = x[2] = xxn; x[1] = x[3] = xxn + dx;

y[0] = y[1] = yyn; y[2] = y[3] = yyn + dy;

koord++;

for (int u = 0; u < 4; u++)

{

uzlK[koord, u * 2] = x[u];

uzlK[koord, u * 2 + 1] = y[u];

}

uzl[koord, 0] = 2 * L;

uzl[koord, 1] = 2 * L + 1;

uzl[koord, 2] = 2 * (L + 1);

uzl[koord, 3] = 2 * (L + 1) + 1;

uzl[koord, 4] = 2 * (L + xr + 1);

uzl[koord, 5] = 2 * (L + xr + 1) + 1;

uzl[koord, 6] = 2 * (L + xr + 2);

uzl[koord, 7] = 2 * (L + xr + 2) + 1;

INIT_k(koord);

k1 = LK(x, y);                   

for (int indG = 0; indG < 8; indG++)

{

for (int ind2 = 0; ind2 < 8; ind2++)

{

// MessageBox.Show(uzl[koord, indG].ToString()+"<=>"+uzl[koord, ind2].ToString());

KG[uzl[koord, indG], uzl[koord, ind2]] += k1[indG, ind2];

}

}                   

L++;

}

L++;

}

//  FileWrite(KG);

for (int i = 0; i < KG.GetLength(0); i += 2 * (xr + 1))

{

int uu = i + xr * 2;

for (int j = 0; j < KG.GetLength(0); j++)

{

KG[i, j] = 0;

KG[j, i] = 0;

KG[i + 1, j] = 0;

KG[j, i + 1] = 0;

KG[uu, j] = 0;

KG[j, uu] = 0;

KG[uu + 1, j] = 0;

KG[j, uu + 1] = 0;

}

KG[uu, uu] = 1;

KG[uu + 1, uu + 1] = 1;

KG[i, i] = 1;

KG[i + 1, i + 1] = 1;

}

F[F.Length - xr - 1] = -1000000;

StreamWriter fw = new StreamWriter("F.txt");

for (int i = 0; i < F.Length; i++)

{

fw.WriteLine("{0:f3}", F[i]);

}

fw.Close();

for (int i = 0; i < 2 * xr+2; i += 2)

{

for (int j = 0; j < KG.GetLength(0); j++)

{

KG[i, j] = 0;

KG[j, i] = 0;

KG[i + 1, j] = 0;

KG[j, i + 1] = 0;

}

KG[i, i] = 1;

KG[i + 1, i + 1] = 1;

}

FileWrite(KG);

button1.Enabled = false;

}       

public void FileWrite(double[,] matr)

{

StreamWriter fw = new StreamWriter("KG.txt");

for (int i = 0; i < matr.GetLength(0); i++)

{

for (int j = 0; j < matr.GetLength(1); j++)

fw.Write("{0,9:f2} ", matr[i, j]);

fw.WriteLine();

}

fw.Close();

}

public void FileWrite(double[] matr)

{

StreamWriter fw = new StreamWriter("u.txt");

StreamWriter fw2 = new StreamWriter("v.txt");

StreamWriter fw3 = new StreamWriter("x.txt");

int i = 0;

for (double yyn = yn; yyn <= yk; yyn += dy)

{

for (double xxn = xn; xxn <= xk; xxn += dx)

{

fw.WriteLine("{0,10:f7}{1,10:f7}{2,20:f7}", xxn, yyn, matr[i]);

fw2.WriteLine("{0,10:f7}{1,10:f7}{2,20:f7}", xxn, yyn, matr[i + 1]);

fw3.WriteLine("{0,5:f2}{1,5:f2}{2,30}", xxn, yyn, matr[i]);

fw3.WriteLine("{0,5:f2}{1,5:f2}{2,30}", xxn, yyn, matr[i + 1]);

i+=2;

}

}

fw2.Close();

fw.Close();

fw3.Close();

}

public void MethodGausa(double[,] a, double[] b, double[] xx)

{

int n = 8;

for (int kk = 0; kk < n; kk++)

{

double max = Math.Abs(a[kk, kk]);

int nmax = kk;

for (int i = kk + 1; i < n; i++)

if (Math.Abs(a[i, kk]) > max)

{

max = Math.Abs(a[i, kk]);

nmax = i;

}

for (int j = 0; j < n; j++)

{

double c = a[kk, j];

a[kk, j] = a[nmax, j];

a[nmax, j] = c;

}

double v = b[kk];

b[kk] = b[nmax];

b[nmax] = v;

for (int i = kk + 1; i < n; i++)

{

double m = a[i, kk] / a[kk, kk];

for (int j = kk; j < n; j++)

a[i, j] = a[i, j] - m * a[kk, j];

b[i] = b[i] - m * b[kk];

}

}

if (a[n - 1, n - 1] == 0)

{

MessageBox.Show("Нельзя вычислить обратную матрицу");              

}

else

for (int i = n - 1; i >= 0; i--)

{

double s = 0;

for (int j = i + 1; j < n; j++)

s = s + a[i, j] * xx[j];

xx[i] = (b[i] - s) / a[i, i];

}

}

private void button3_Click(object sender, EventArgs e)

{

double[] resultat;

PramHod(KG.GetLength(0));

resultat = ObratHod(ref KG, ref F);

FileWrite(resultat);

button3.Enabled = false;

DrawingElement2(resultat);

if(xr<7)

DrawingElement3(resultat);

// MessageBox.Show("Вычисления сохранены в файл");

}

public void PramHod(int n)

{

for (int kk = 0; kk < n; kk++)

{

double max = Math.Abs(KG[kk, kk]);

int nmax = kk;

for (int i = kk + 1; i < n; i++)

if (Math.Abs(KG[i, kk]) > max)

{

max = Math.Abs(KG[i, kk]);

nmax = i;

}

for (int j = 0; j < n; j++)

{

double c = KG[kk, j];

KG[kk, j] = KG[nmax, j];

KG[nmax, j] = c;

}

double v = F[kk];

F[kk] = F[nmax];

F[nmax] = v;

for (int i = kk + 1; i < n; i++)

{

double m = KG[i, kk] / KG[kk, kk];

for (int j = kk; j < n; j++)

KG[i, j] = KG[i, j] - m * KG[kk, j];

F[i] = F[i] - m * F[kk];

}

}

for (int kk = 0; kk < n; kk++)

{

for (int i = 0; i < kk; i++)

KG[kk, i] = 0;

}

}

public double[] ObratHod(ref double[,] a, ref double[] b)

{

double[] xx = new double[KG.GetLength(0)];

int i, j;

for (i = 0; i < a.GetLength(0); i++)

xx[i] = 0;

for (i = 0; i < a.GetLength(0); i++)

for (j = 0; j < i; j++)

a[i, j] = 0;

i = 0; j = 0;

double sum, perem;

int n = a.GetLength(0);

int m = a.GetLength(1) + 1;

double[,] aa = new double[n, m];

for (i = 0; i < n; i++)

for (j = 0; j < n; j++)

aa[i, j] = a[i, j];

for (i = 0; i < n; i++)

{

aa[i, m - 1] = b[i];

xx[i] = 0;

}

int t = n;

for (i = n - 1; i >= 0; i--)

{

sum = 0;

for (int p = 0; p < m; )

{

if (aa[i, p] == 0)

p++;

else

{

j = m - 2;

while (j != p)

{

sum += aa[i, j] * xx[j];

j--;

}

perem = (aa[i, m - 1] - sum) / aa[i, p];

xx[t - 1] = perem;

t--;

break;

}

}

}

return xx;

}

//MKR

double[,] AA;

public double[] BB,mkr;

private void button4_Click(object sender, EventArgs e)

{

double[,] A1={

{13,-6,1,-7,2,0},

{-6,12,-6,2,-7,2},

{1,-6,13,0,2,-7},

{-7,2,0,22,-8,1},

{2,-7,2,-8,21,-8},

{0,2,-7,1,-8,22}

};

double[] b1 = { 14,-3,0,15,3,15};

GLeft();

GRight();

Initial();

Centr();

int n = 6;

for (int i = 0; i < 6; i++)

b1[i] *= dx * dy;

BB[BB.Length-xr-1]=-1000000;

for (int kk = 0; kk < n; kk++)

{

double max = Math.Abs(A1[kk, kk]);

int nmax = kk;

for (int i = kk + 1; i < n; i++)

if (Math.Abs(A1[i, kk]) > max)

{

max = Math.Abs(A1[i, kk]);

nmax = i;

}

for (int j = 0; j < n; j++)

{

double c = A1[kk, j];

A1[kk, j] = A1[nmax, j];

A1[nmax, j] = c;

}

double v = b1[kk];

b1[kk] = b1[nmax];

b1[nmax] = v;

for (int i = kk + 1; i < n; i++)

{

double m = A1[i, kk] / A1[kk, kk];

for (int j = kk; j < n; j++)

A1[i, j] = A1[i, j] - m * A1[kk, j];

b1[i] = b1[i] - m * b1[kk];

}

}

for (int kk = 0; kk < n; kk++)

{

for (int i = 0; i < kk; i++)

A1[kk, i] = 0;

}          

mkr = ObratHod(ref A1, ref b1);

FileWrite(mkr);

//DrawingElement2(mkr);

button4.Enabled = false;

}

public void Centr()

{

int num = 0;

double dyy = 1.0f / (dy * dy);

double dxx = 1.0f / (dx * dx);          

for (int i = 1; i < yr; i++)

{

for (int j = 1; j < xr; j++)

{

num = j + i * (xr + 1);

AA[num, num] = -2*(dxx+dyy);

AA[num+xr+1, num+xr+1] = -2 * (dxx + dyy);

if (i == 1)

{

if (j == 1)

{

AA[num, num + 1] = dxx;

BB[num] = dyy * BB[num - xr - 1] + dxx * BB[num - 1];

}

else

{

AA[num, num - 1] = dxx;

if (j == xr - 1)

{

BB[num] = dyy * BB[num - xr - 1] + dxx * BB[num + 1];

}

else

{

AA[num, num + 1] = dxx;

BB[num] = dyy * BB[num - xr - 1];                                    

}

}

}

else

{

AA[num, num - xr - 1] = dyy;

if (j == 1)

{

AA[num, num + 1] = dxx;                          

BB[num] = dxx * BB[num - 1];

}

else

{

AA[num, num - 1] = dxx;

if (j == xr - 1)

{

BB[num] = dxx * BB[num + 1];

}

else

{

AA[num, num + 1] = dxx;

BB[num] = 0;

}

}

}

}

}

for (int j = 1; j < xr; j++)

{

num = j + yr * (xr + 1);

AA[num, num] = dyy;

AA[num - xr - 1, num - xr - 1] = -2 * (dxx + dyy);

AA[num - 2 * xr - 2, num - 2*xr - 2] = dyy;

if (j == 1)

{

AA[num, num + 1] = dxx;

BB[num] = dxx * BB[num - 1];

}

else

{

AA[num, num - 1] = dxx;

if (j == xr - 1)

{

BB[num] = dxx * BB[num + 1];

}

else

{

AA[num, num + 1] = dxx;

BB[num] = 0;

}

}

}

}

public void Initial()

{

for (int i = 0; i < BB.Length; i++)

{

AA[i, i] = 1;

BB[i] = 0.5f;// Int32.Parse(comboBox1.SelectedItem.ToString());

}

}

public void GLeft()

{

for (int i = 0; i < BB.Length; i += xr + 1)

{

AA[i, i] = 1;

BB[i] = 0;// Int32.Parse(listBox1.SelectedItem.ToString());

}

}

public void GRight()

{

for (int i = xr; i < BB.Length; i += xr + 1)

{

AA[i, i] = 1;

BB[i] = 0;// Int32.Parse(domainUpDown1.Text.ToString());

}

}

}

}

Результат выполнения программы:

Выводы: во время выполнения лабораторной работы была построена