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

AddEdge(P2, P3);

AddEdge(P1, P3);

end;

function TModel.AddVertex(aX, aY: Extended): TRealPoint;

begin

//Если такая точка уже есть

Result := FindPoint(aX, aY);

if Result <> nil then exit; 

Result := TRealPoint.Create(Vertexes.Count, aX, aY);

Vertexes.Add(Result);

{}{}{}

end;

procedure TModel.Clear;

begin

//Очищаем все списки модели

Vertexes.Clear;

Edges.Clear;

Trias.Clear;

end;

constructor TModel.Create;

begin

inherited Create;

Vertexes := TList.Create;

Edges := TList.Create;

Trias := TList.Create;

end;

//----------------------------------------------//ПРОЦЕДУРА РИСОВАНИЯ ОБЪЕКТА

//----------------------------------------------procedure TModel.Draw(C: TCanvas; Scale: Extended);

var

I: Integer;

P: TRealPoint;

L: TRealLine;

T: TTria;

procedure Marker(X, Y: Extended);

begin

C.Ellipse(Round(X*Scale), Round(Y*Scale), Round(X*Scale)+5, Round(Y*Scale)+5);

end;

procedure Circle(X, Y, R: Extended);

begin

C.Ellipse(Round((X-R)*Scale), Round((Y-R)*Scale), Round((X+R)*Scale), Round((Y+R)*Scale));

end;

procedure Line(X1, Y1, X2, Y2: Extended);

begin

C.MoveTo(Round(X1*Scale), Round(Y1*Scale));

C.LineTo(Round(X2*Scale), Round(Y2*Scale));

end;

procedure Triangle(X1, Y1, X2, Y2, X3, Y3: Extended);

begin

C.Polygon([

Point(Round(X1*Scale), Round(Y1*Scale)),

Point(Round(X2*Scale), Round(Y2*Scale)),

Point(Round(X3*Scale), Round(Y3*Scale))

]);

end;

var

KD: Extended;

begin

C.Brush.Color := clWhite;

C.FillRect(Rect(0,0, 1000,1000));

C.Pen.Color := clBlack;

C.Brush.Color := clYellow;

if Kmax = 0 then

KD := 0

else

KD := 5/Kmax;

for I := 0 to Trias.Count - 1 do

begin

T := TTria(Trias[I]);

Triangle(

T.Points[1].X, T.Points[1].Y,

T.Points[2].X, T.Points[2].Y,

T.Points[3].X, T.Points[3].Y);

end;

for I := 0 to Edges.Count - 1 do

begin

L := TRealLine(Edges[I]);

Line(L.FirstPoint.X, L.FirstPoint.Y, L.EndPoint.X, L.EndPoint.Y);

end;

C.Pen.Color := clBlue;

for I := 0 to Trias.Count - 1 do

begin

C.Brush.Color := clRed;

T := TTria(Trias[I]);

Circle((T.X(1)+T.X(2)+T.X(3))/3, (T.Y(1)+T.Y(2)+T.Y(3))/3, T.EpsSum*KD);

end;

end;

//------------------------------КОНЕЦ ПРОЦЕДУРЫ РИСОВАНИЯ

function TModel.FindPoint(aX, aY: Extended): TRealPoint;

var

I: Integer;

begin

for I := 0 to Vertexes.Count - 1 do

with TRealPoint(Vertexes[I])do

if (Abs(X - aX)<1e-3)and(Abs(Y - aY)<1e-3)then

begin

Result := TRealPoint(Vertexes[I]);

exit

end;

Result := nil;   

end;

//Воспомогательные функции - смещение переменной DXi в матрице А

function DXOffset(P: Integer): Integer;

begin

Result := P*2+1;

end;

//Воспомогательные функции - смещение переменной DYi в матрице А

function DYOffset(P: Integer): Integer;

begin

Result := P*2+2;

end;

//ГЛАВНАЯ ПРОЦЕДУРА РЕШЕНИЯ СЛАУ МКЭ

procedure TModel.SolveMKE(E, v, d: Extended; Method: TMethod);

var

SL: TSlau;

VertexID, I, N, VertexIndex, TriaID, PointInTria: Integer;

P: TRealPoint;

T: TTria;

CurEqX, CurEqY: TVector;

EPS,DX,X: TVector;

B, K: TMatrix;

ForceSum, ForceX, ForceY: Extended;

begin

//Определим размер системы (по двы уравнения на каждый узел)

N := 2*Vertexes.Count;

SL := TSlau.Create(N);

//Установим длину воспомогательных векторов 

SetLength(CurEqX,N+1);

SetLength(CurEqY,N+1);

//Цикл по всем узлам объекта

for VertexID := 0 to Vertexes.Count - 1 do

begin

//Присвоим Р - текущий узел

P := TRealPoint(Vertexes[VertexID]);

//Очистим векторы текущих уравнений

for I := 1 to N do

begin

CurEqX[I] := 0;

CurEqY[I] := 0;

end;

//Для каждого треугольника Т, включающего в себя текущий узел

for TriaID := 0 to Trias.Count - 1 do

begin

T := TTria(Trias[TriaID]);

//Найдем положение текущего узла в треугольнике (1 2 или 3)

VertexIndex := T.FindPoint(P);

//Если оно меньше 0 - узел не входит в треугольник, идем дальше

if VertexIndex <= 0 then continue;

//Сформируем матрицу жесткости треугольника Т

T.KMatrix(E, v, d, K);

//Для каждой из вершин треугольника

for PointInTria := 1 to 3 do

begin

//Добавим влияние сил в вершинах треугольника

//На смещение текущего! узла P (в треугольнике он является VertexIndex-м)

//Поэтому умножаем перемещение PointInTria-го узла на К[VertexIndex,PointInTria]