Форум программистов [Powered by Invision Power Board]

Код

unit lsm3d;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, Math, StdCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Memo1: TMemo;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

const Points = 5;
      h = 0.1;
      xmax = Round(100/h);
      ymax = Round(100/h);
      zmax = Round(100/h);
      dt = 0.00025;
      b = 1;
      dtau = 0.1*h;
var
  Form1      : TForm1;
  dnosign    : array[1..xmax,1..ymax,1..zmax] of single;
  dsign      : array[1..xmax,1..ymax,1..zmax] of single;
  phi        : array[1..xmax,1..ymax,1..zmax] of single;
  phi1        : array[1..xmax,1..ymax,1..zmax] of single;
  X,Y,Z      : array[1..Points] of single;
  graf        : array[1..xmax,1..ymax,1..zmax] of shortint;
  



implementation

{$R *.dfm}


procedure TForm1.Button1Click(Sender: TObject);
  var i,j,k,q, maxx, maxy, maxz, minx, miny, minz : Integer;
      temp1,temp2 : single;
      a1,a2,a3,a4,a5,a6 : Integer;
      grad_plus, grad_minus, phi_x_plus, phi_x_minus, phi_y_plus, phi_y_minus, phi_z_plus, phi_z_minus : single;
begin


  X[1]  := 70.0; X[2]  := 69.7; X[3]  := 69.1; X[4]  := 67.7; X[5]  := 66.2;
  Y[1]  := 50.1; Y[2]  := 53.2; Y[3]  := 56.2; Y[4]  := 59.1; Y[5]  := 61.8;
  Z[1]  := 30.0; Z[2]  := 30.3; Z[3]  := 31.0; Z[4]  := 32.2; Z[5]  := 33.8;

  // Вычисление расстояния от каждой точки сетки до множества данных
  for i:=1 to xmax  do
    for j:=1 to ymax  do
      for q:=1 to zmax do
        begin
          temp1:=1000;
          for k:=1 to Points do begin
            temp2 := sqrt((i*h - X[k])*(i*h - X[k]) + (j*h - Y[k])*(j*h - Y[k]) + (q*h-Z[k])*(q*h-Z[k]));
            if temp1>temp2 then temp1 := temp2;
          end;
          dnosign[i,j,q] := temp1;
        end;

  // Нахождение минимальных и максимальных элементов в массивах X и Y
  minx:=1; miny:=1; minz := 1; maxx:=1; maxy:=1; maxz:=1;
  for k:=2 to Points do
    begin
      if X[k]<X[minx] then minx := k;
      if Y[k]<Y[miny] then miny := k;
      if Z[k]<Z[minz] then minz := k;
      if X[k]>X[maxx] then maxx := k;
      if Y[k]>Y[maxy] then maxy := k;
      if Z[k]>Z[maxz] then maxz := k;
    end;

  // Задание границ квадрата, ограничивающего исходное множество данных
  a1 := Round(X[minx]/h)-Round(500*h);
  a2 := Round(X[maxx]/h)+Round(500*h);
  a3 := Round(Y[miny]/h)-Round(500*h);
  a4 := Round(Y[maxy]/h)+Round(500*h);
  a5 := Round(Z[minz]/h)-Round(500*h);
  a6 := Round(Z[maxz]/h)+Round(500*h);

  // Вычисление расстояния со знаком
  for i:=1 to xmax  do
    for j:=1 to ymax  do
      for q:=1 to zmax do
        begin
          if((i>a1)and(i<a2)and(j>a3)and(j<a4)and(q>a5)and(q<a6)) then
            dsign[i,j,q] := - dnosign[i,j,q]
          else
            dsign[i,j,q] := dnosign[i,j,q];
        end;

  // Задание начального приближения
  for i:=1 to xmax  do
    for j:=1 to ymax  do
      for q:=1 to zmax do
        phi[i,j,q] := dsign[i,j,q];


  // Решение уравнения LS
  for k:= 1 to 12000 do
    begin
      
    for i:=1 to xmax  do
      for j:=1 to ymax  do
        for q:=1 to zmax do

          begin

          // Вычисление производных
          
          if(i>1)and(i<xmax)and(j>1)and(j<ymax)and(q>1)and(q<zmax) then
          begin
            phi_x_plus  := (phi[i+1,j,q]-phi[i,j,q])/h;
            phi_x_minus := (phi[i,j,q]-phi[i-1,j,q])/h;
            phi_y_plus  := (phi[i,j+1,q]-phi[i,j,q])/h;
            phi_y_minus := (phi[i,j,q]-phi[i,j-1,q])/h;
            phi_z_plus  := (phi[i,j,q+1]-phi[i,j,q])/h;
            phi_z_minus := (phi[i,j,q]-phi[i,j,q-1])/h;
          end;


          grad_plus  := sqrt(max(phi_x_minus,0)*max(phi_x_minus,0) + min(phi_x_plus,0)*min(phi_x_plus,0) + max(phi_y_minus,0)*max(phi_y_minus,0) + min(phi_y_plus,0)*min(phi_y_plus,0) + max(phi_z_minus,0)*max(phi_z_minus,0) + min(phi_z_plus,0)*min(phi_z_plus,0));

          grad_minus := sqrt(max(phi_x_plus,0)*max(phi_x_plus,0) + min(phi_x_minus,0)*min(phi_x_minus,0) + max(phi_y_plus,0)*max(phi_y_plus,0) + min(phi_y_minus,0)*min(phi_y_minus,0) + max(phi_z_plus,0)*max(phi_z_plus,0) + min(phi_z_minus,0)*min(phi_z_minus,0));

          phi1[i,j,q] := phi[i,j,q] - dt*(max(dsign[i,j,q],0)*grad_plus + min(dsign[i,j,q],0)*grad_minus);

        end;//i,j,q

    end;//k

    for i:=1 to xmax  do
      for j:=1 to ymax  do
        for q :=1 to zmax do
          begin
            if(phi1[i,j,q]<=0) then graf[i,j,q]:=-1
            else graf[i,j,q] := 1;
          end;

  end;//procedure

end.