• 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 
 • 







.

, [1].
(1950 .).
, , .
;
.
1963 . , ,
-.


[2,3].
,

(),
()

, A
, , n2 - . ,
, [5].
1 ,


, , , ,
, ,
- , .
-
[1,2,3].


. , , ,
,
. . ,
:
1.
. .
2.
, .
3.
, .
.
4 .
,
.
, ,
.

.

. ,
,

-. ,
, , .
1.1

[4,5].
n-
n- .
,
:
, d
.
j bj,
n d1,...,dn.
. n n
, ,
, :
x1=d1/d; x2=d2/d;....; xn-1=dn-1/d; xn=dn/d;
.

,
n . ,
(3) rmin{m,n}, r
, m-r .
, r
r .
, m - r
(3)
r (
),
. r
n
, r- ,
r , r 0, . .
. ,
, , n - r
- .
(4)
xr+1,..., xn. , r
r .
Mr
, .
,
.
.

n

, - Kane
, , .
(5) ,
.
, (5) , . . r ()< n,
.
(5) r ,
. r
, n - r
. , ,
r x1,..., r n - r
.
(5) .
, .
n - r
.
.
n n

(6) .
apq ,
, mi=-aiq/apq
ip ( - , ,
).
i-
, mi;
.
,
q- , apq,
.
p-
, .
, ..
,
, - (
). xi
, .

n .
, [6,7].
.

.
Ax=b (7)
- n, x,b -
.

,
..
=,

,
ij bij :
,
(7) :
CBx=b. (9)
Bx B - x
-, y:
Bx=y. (10)
(9)
:
Cy=b. (11)
ij ,
(7)
.
(11),
:
yi
bij.
yi
(12),
(10).
bij
(8), , ,
:
,
, :
, ,
,
.
.

Ax = b
m x m.
, ,
,
,

, .

, . , ,
, ..



.
- (2)
, ()- .. .
.

, -
.
1.2
(
).


.
().

.
( ).
n n :
=b, (14)
, aii
0 (i = 2, ...,
n), xi x2 -
. . , (14):
; , i == 1, 2, ...,n; j == 1,2,..., n.
(15)
(16)
. . (k+1)-

x(0),...,x(k)
, (15),
, .. [4,6].
.

. (k+1)-
xi (i>1) (k+1)-
xi-1.
,
:
(17)

(17).
.
, .. .
, , k- , (k+1)-
:
k=0,1,...,n
.
(1), ,

,
[4],
:
(18)


. , ,
, ..
.
, .. .

, (19)
- .

, :
(20)

k .

. ,
, (18)
. ,

(1),
(1), -
, .
, (19) (20) [10,11].
2


, , ,
.

. ,
. .
,
,
. , - , . 1,
: 1 2 3 4 5 6 7 1, 2, 5, 6, 7 1, 2, 3, 6 2, 3, 4, 6 3, 4, 5, 6, 7 1, 4, 5, 7 1, 2, 3, 4, 6, 7 1, 4, 5, 6, 7
,
,
,
. . 2
1 [9].

, -
, 1.

-
.
,
.
,

.

, - .
:
, (*)
m (m=1,2,3).

, (*).
3


[12] (. 4).


- .

,
.

FORL [5].
( ) . 5.

- 880 2016 .
,
2,7 .
315 .

Pentium 166 32 .
1.

1. () 280 2.2101 -2.4608 1.3756 -5.2501 1.7406 -2.3489 150 2.2137 -2.4669 1.3904 -5.2572 1.7433 -2.3883
1 ,
,
,
.
.

()
.
, ( )
.

,
, (),
, ,
,
.



.
.
1.
., .
// .: , 1980
2.
., // .:
., 1975
3.
., .
// .: , 1977
4.
.., .., ..
// .: , 1987
5.

.., .. // .:, 1984
6.

.. // .: , 1975
7.
..
// : , 1980
8.

.., ..
// i i 1997.
4.
9.
F.G. Gustavson, Some basic techniques for
solving sparse matrix algorithms, // editer by D.J. Rose and R.A.Willoughby,
Plenum Press, New York, 1972
10.
.,
. , // , ,
1984
11.
D.J. Rose, A graph theoretic study of the
numerical solution of sparse positive definite system of linear equations //
New York, Academic Press, 1972
12.

.., .., ..,
// .:, 1978
1
, -
.
#include
<math.h>
#include
<stdio.h>
#include
<string.h>
#include
<stdlib.h>
#include
<fstream.h>
#include
"matrix.h"
#define
BASE3D_4 4
#define
BASE3D_8 8
#define BASE3D_10
10
const double Eps
= 1.0E-10;
DWORD CurrentType
= BASE3D_10;
void
PrintHeader(void)
{
printf("Command format: AConvert
-<switch> <file1.in> <file2.in> <file3.out>
[/Options]\n");
printf("Switch: -t10 -
Tetraedr(10)\n");
printf(" -c8 -
Cube(8)\n");
printf(" -s4 -
Shell(4)\n");
printf(" -s8 -
Shell(8)\n\n");
printf("Optins: /8 - convert
Tetraedr(10)->8*Tetraedr(4)\n");
printf(" /6 - convert Cube(8)->6*Tetraedr(4)\n");
}
bool Output(char*
fname,Vector<double>& x,Vector<double>&
y,Vector<double>& z, Matrix<DWORD>& tr, DWORD n,
DWORD NumNewPoints,DWORD
ntr,Matrix<DWORD>& Bounds,DWORD CountBn)
{
char*
Label = "NTRout";
int
type = CurrentType,
type1 = (type==BASE3D_4 ||
type==BASE3D_10) ? 3 : 4;
DWORD
NewSize,
i,
j;
ofstream Out;
if (type==BASE3D_4) type1 = 3;
else if (type==BASE3D_8) type1 = 4;
else if (type==BASE3D_10) type1 = 6;
Out.open(fname,ios::out | ios:: binary);
if (Out.bad()) return true;
Out.write((const char*)Label,6 *
sizeof(char));
if (Out.fail()) return true;
Out.write((const
char*)&type,sizeof(int));
if (Out.fail()) return true;
Out.write((const
char*)&CountBn,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
Out.write((const char*)&(NewSize = n +
NumNewPoints),sizeof(DWORD));
if (Out.fail()) return true;
Out.write((const
char*)&(NumNewPoints),sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
for (DWORD i = 0; i < n; i++)
{
Out.write((const
char*)&x[i],sizeof(double));
Out.write((const
char*)&y[i],sizeof(double));
Out.write((const
char*)&z[i],sizeof(double));
if (Out.fail())
{
Out.close();
return true;
}
}
for (i = 0; i < NumNewPoints; i++)
{
Out.write((const char*)&x[n +
i],sizeof(double));
Out.write((const char*)&y[n +
i],sizeof(double));
if (Out.fail())
{
Out.close();
return true;
}
}
Out.write((const
char*)&(ntr),sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
for (i = 0; i < ntr; i++)
for (j = 0; j < (DWORD)type; j++)
{
DWORD out = tr[i][j];
Out.write((const char*)&out,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
}
for (i = 0; i < CountBn; i++)
for (j = 0; j < (DWORD)type1; j++)
{
DWORD out = Bounds[i][j];
Out.write((const char*)&out,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
}
{
//*********************
// Create Links
printf("Create links...\r");
Vector<DWORD> Link(n),
Size(n);
Matrix<DWORD> Links(n,n);
DWORD Count;
int type = CurrentType;
for (DWORD i = 0; i < n; i++)
{
for (DWORD j = 0; j < ntr; j++)
for (DWORD k = 0; k < (DWORD)type;
k++)
if (tr[j][k] == i)
for (DWORD m = 0; m < (DWORD)type;
m++) Link[tr[j][m]] = 1;
Count = 0;
for (DWORD m = 0; m < n; m++)
if (Link[m]) Count++;
Size[i] = Count;
Count = 0;
for (DWORD m = 0; m < n; m++)
if (Link[m])
Links[i][Count++] = m;
//Set zero
Link.ReSize(n);
}
// Output
//*********************
for (DWORD i = 0; i < n; i++)
{
DWORD Sz = Size[i];
Out.write((const
char*)&Sz,sizeof(DWORD));
for (DWORD j = 0; j < Sz; j++)
Out.write((const char*)&(Links[i][j]),sizeof(DWORD));
}
//*********************
}
printf(" \r");
printf("Points: %ld\n",n);
printf("FE: %ld\n",ntr);
Out.close();
return false;
}
bool Test(DWORD*
a,DWORD* b)
{
bool
result;
int
NumPoints = 3;
if (CurrentType == BASE3D_8) NumPoints = 4;
else if (CurrentType == BASE3D_10) NumPoints
= 6;
for (int i = 0; i < NumPoints; i++)
{
result = false;
for (int j = 0; j < NumPoints; j++)
if (b[j] == a[i])
{
result = true;
break;
}
if (result == false) return false;
}
return true;
}
void
Convert(Vector<double>& X,Vector<double>&
Y,Vector<double>& Z, Matrix<DWORD>& FE,DWORD
NumTr,Matrix<DWORD>& Bounds,DWORD& BnCount)
{
int
cData8[6][5] = {{0,4,5,1,7},
{6,2,3,7,0},
{4,6,7,5,0},
{2,0,1,3,5},
{1,5,7,3,4},
{6,4,0,2,1}},
cData4[4][4] = {{0,1,2,3},
{1,3,2,0},
{3,0,2,1},
{0,3,1,2}},
cData10[4][7] =
{{0,1,2,4,5,6,3},
{0,1,3,4,8,7,2},

{1,3,2,8,9,5,0},

{0,2,3,6,9,7,1}},
cData[6][7],
Data[6],
l,
Num1,
Num2,
m;
DWORD
i,
j,
p[6],
pp[6],
Index;
Matrix<DWORD> BoundList(4 * NumTr,6);
double
cx,
cy,
cz,
x1,
y1,
z1,
x2,
y2,
z2,
x3,
y3,
z3;
Bounds.ReSize(4 * NumTr,6);
switch (CurrentType)
{
case BASE3D_4:
Num1 = 4;
Num2 = 3;
for (l = 0; l < Num1; l++)
for (m = 0; m < Num2+1; m++)
cData[l][m] = cData4[l][m];
break;
case BASE3D_8:
Num1 = 6;
Num2 = 4;
for (l = 0; l < Num1; l++)
for (m = 0; m < Num2 + 1; m++)
cData[l][m] = cData8[l][m];
break;
case BASE3D_10:
Num1 = 4;
Num2 = 6;
for (l = 0; l < Num1; l++)
for (m = 0; m < Num2+1; m++)
cData[l][m] = cData10[l][m];
}
printf("Create bounds...\r");
for (i = 0; i < NumTr - 1; i++)
for (int j = 0; j < Num1; j++)
if (!BoundList[i][j])
{
for (l = 0; l < Num2; l++)
p[l]
= FE[i][cData[j][l]];
for (DWORD k = i + 1; k < NumTr; k++)
for (int m = 0; m < Num1; m++)
if (!BoundList[k][m])
{
for (int l = 0; l < Num2; l++)
pp[l] = FE[k][cData[m][l]];
if (Test(p,pp))
BoundList[i][j] = BoundList[k][m] =
1;
}
}
for (i = 0; i < NumTr; i++)
for (j = 0; j < (DWORD)Num1; j++)
if (BoundList[i][j] == 0)
{
if (CurrentType == BASE3D_4)
{
cx = X[FE[i][cData[j][3]]];
cy = Y[FE[i][cData[j][3]]];
cz = Z[FE[i][cData[j][3]]];
}
else
if (CurrentType == BASE3D_10)
{
cx = X[FE[i][cData[j][6]]];
cy = Y[FE[i][cData[j][6]]];
cz = Z[FE[i][cData[j][6]]];
}
else
{
cx = X[FE[i][cData[j][4]]];
cy = Y[FE[i][cData[j][4]]];
cz = Z[FE[i][cData[j][4]]];
}
x1 = X[FE[i][cData[j][0]]];
y1 = Y[FE[i][cData[j][0]]];
z1 = Z[FE[i][cData[j][0]]];
x2 = X[FE[i][cData[j][1]]];
y2 = Y[FE[i][cData[j][1]]];
z2 = Z[FE[i][cData[j][1]]];
x3 = X[FE[i][cData[j][2]]];
y3 = Y[FE[i][cData[j][2]]];
z3 = Z[FE[i][cData[j][2]]];
for (l = 0; l < Num2; l++)
Data[l] = cData[j][l];
if ( ((cx-x1)*(y2-y1)*(z3-z1) +
(cy-y1)*(z2-z1)*(x3-x1) + (y3-y1)*(cz-z1)*(x2-x1) -
(x3-x1)*(y2-y1)*(cz-z1) - (y3-y1)*(z2-z1)*(cx-x1) -
(cy-y1)*(z3-z1)*(x2-x1)) > 0)
{
if (CurrentType == BASE3D_4)
{
Data[0] = cData[j][0];
Data[1] = cData[j][2];
Data[2] = cData[j][1];
}
else
if (CurrentType == BASE3D_10)
{
Data[0] = cData[j][0];
Data[1] = cData[j][2];
Data[2] = cData[j][1];
Data[3] = cData[j][5];
Data[5] = cData[j][3];
}
else
{
Data[0] = cData[j][0];
Data[1] = cData[j][3];
Data[2] = cData[j][2];
Data[3] = cData[j][1];
}
}
for (l = 0; l < Num2; l++)
Bounds[BnCount][l] = FE[i][Data[l]];
BnCount++;
}
}
void main(int
argc,char** argv)
{
char *input1,
*input2,
*input3,
*op = "",
*sw;
bool CreateFile(char*,char*,char*,char*);
printf("ANSYS->FORL file convertor.
ZSU(c) 1998.\n\n");
if (argc < 5 || argc > 6)
{
PrintHeader();
return;
}
sw
= argv[1];
input1 = argv[2];
input2 = argv[3];
input3 = argv[4];
if (!strcmp(sw,"-t10"))
CurrentType = BASE3D_10;
else
if (!strcmp(sw,"-c8"))
CurrentType = BASE3D_8;
else
{
printf("Unknown switch
%s\n\n",sw);
PrintHeader();
return;
}
if (argc == 6)
{
op = argv[5];
if (strcmp(op,"/8") &&
strcmp(op,"/6"))
{
printf("Unknown options
%s\n\n",op);
PrintHeader();
return;
}
}
if (CreateFile(input1,input2,input3,op))
printf("OK\n");
}
bool
CreateFile(char* fn1,char* fn2,char* fn3,char* Op)
{
FILE *in1,
*in2,
*in3;
Vector<double> X(1000),
Y(1000),
Z(1000);
DWORD NumPoints,
NumFE,
NumBounds = 0,
tmp;
Matrix<DWORD> FE(1000,10),
Bounds;
bool
ReadTetraedrData(char*,char*,FILE*,FILE*,Vector<double>&,Vector<double>&,Vector<double>&,

Matrix<DWORD>&,DWORD&,DWORD&),

ReadCubeData(char*,char*,FILE*,FILE*,Vector<double>&,Vector<double>&,Vector<double>&,

Matrix<DWORD>&,DWORD&,DWORD&);
void Convert824(Matrix<DWORD>&,DWORD&),

Convert1024(Matrix<DWORD>&,DWORD&);
if ((in1 = fopen(fn1,"r")) ==
NULL)
{
printf("Unable open file
%s",fn1);
return false;
}
if ((in2 = fopen(fn2,"r")) ==
NULL)
{
printf("Unable open file
%s",fn2);
return false;
}
if (CurrentType == BASE3D_10)
{
if
(!ReadTetraedrData(fn1,fn2,in1,in2,X,Y,Z,FE,NumPoints,NumFE)) return false;
if (!strcmp(Op,"/8"))
{
// Create 8*Tetraedr(4)
Convert1024(FE,NumFE);
}

Convert(X,Y,Z,FE,NumFE,Bounds,NumBounds);
return
!Output(fn3,X,Y,Z,FE,NumPoints,0,NumFE,Bounds,NumBounds);
}
if (CurrentType == BASE3D_8)
{
if
(!ReadCubeData(fn1,fn2,in1,in2,X,Y,Z,FE,NumPoints,NumFE)) return false;
if (!strcmp(Op,"/6"))
{
// Create 6*Tetraedr(4)
Convert824(FE,NumFE);
}

Convert(X,Y,Z,FE,NumFE,Bounds,NumBounds);
return
!Output(fn3,X,Y,Z,FE,NumPoints,0,NumFE,Bounds,NumBounds);
}
return false;
}
void
Convert824(Matrix<DWORD>& FE,DWORD& NumFE)
{
Matrix<DWORD> nFE(6 * NumFE,4);
DWORD
data[][4] = {
{ 0,2,3,6 },
{ 4,5,1,7 },
{ 0,4,1,3 },
{ 6,7,3,4 },
{ 1,3,7,4 },
{ 0,4,6,3 }
},
Current = 0;
for (DWORD i = 0; i < NumFE; i++)
for (DWORD j = 0; j < 6; j++)
{
for (DWORD k = 0; k < 4; k++)
nFE[Current][k] = FE[i][data[j][k]];
Current++;
}
CurrentType = BASE3D_4;
NumFE = Current;
FE
= nFE;
}
void Convert1024(Matrix<DWORD>&
FE,DWORD& NumFE)
{
Matrix<DWORD> nFE(8 * NumFE,4);
DWORD
data[][4] = {
{ 3,7,8,9 },
{ 0,4,7,6 },
{ 2,5,9,6 },
{ 7,9,8,6 },
{ 4,8,5,1 },
{ 4,5,8,6 },
{ 7,8,4,6 },
{ 8,9,5,6 }
},
Current = 0;
for (DWORD i = 0; i < NumFE; i++)
for (DWORD j = 0; j < 8; j++)
{
for (DWORD k = 0; k < 4; k++)
nFE[Current][k] = FE[i][data[j][k]];
Current++;
}
CurrentType = BASE3D_4;
NumFE = Current;
FE
= nFE;
}
bool
ReadTetraedrData(char* fn1,char* fn2,FILE* in1,FILE*
in2,Vector<double>& X,Vector<double>& Y,Vector<double>&
Z,
Matrix<DWORD>&
FE,DWORD& NumPoints,DWORD& NumFE)
{
double tx,
ty,
tz;
char TextBuffer[1001];
DWORD Current = 0L,
tmp;
while (!feof(in1))
{
if (fgets(TextBuffer,1000,in1) == NULL)
{
if (feof(in1)) break;
printf("Unable read file
%s",fn1);
fclose(in1);
fclose(in2);
return false;
}
if (sscanf(TextBuffer,"%ld %lf %lf
%lf", &NumPoints,&tx,&ty,&tz) != 4) continue;
X[Current] = tx;
Y[Current] = ty;
Z[Current] = tz;
if (++Current == 999)
{
Vector<double> t1 = X,
t2 = Y,
t3 = Z;
X.ReSize(2 * X.Size());
Y.ReSize(2 * X.Size());
Z.ReSize(2 * X.Size());
for (DWORD i = 0; i < Current;
i++)
{
X[i] = t1[i];
Y[i] = t2[i];
Z[i] = t3[i];
}
}
if (Current % 100 == 0)
printf("Line:
%ld\r",Current);
}
fclose(in1);
printf(" \r");
NumPoints = Current;
Current = 0L;
while (!feof(in2))
{
if (fgets(TextBuffer,1000,in2) == NULL)
{
if (feof(in2)) break;
printf("Unable read file
%s",fn2);
fclose(in2);
return false;
}
if (sscanf(TextBuffer,"%d %d %d %d
%d %ld %ld %ld %ld %ld %ld %ld %ld",

&tmp,&tmp,&tmp,&tmp,&tmp,
&FE[Current][0],&FE[Current][1],&FE[Current][2],&FE[Current][3],

&FE[Current][4],&FE[Current][5],&FE[Current][6],&FE[Current][7])
!= 13) continue;
if (fgets(TextBuffer,1000,in2) == NULL)
{
printf("Unable read file
%s",fn2);
fclose(in2);
return false;
}
if (sscanf(TextBuffer,"%ld
%ld",&FE[Current][8],&FE[Current][9]) != 2)
{
printf("Unable read file
%s",fn2);
fclose(in2);
return false;
}
{
if (fabs((tx = 0.5*(X[FE[Current][0] -
1] + X[FE[Current][1] - 1])) - X[FE[Current][4] - 1]) > Eps)
X[FE[Current][4] - 1] = tx;
if (fabs((ty = 0.5*(Y[FE[Current][0] -
1] + Y[FE[Current][1] - 1])) - Y[FE[Current][4] - 1]) > Eps)
Y[FE[Current][4] - 1] = ty;
if (fabs((tz = 0.5*(Z[FE[Current][0] -
1] + Z[FE[Current][1] - 1])) - Z[FE[Current][4] - 1]) > Eps)
Z[FE[Current][4] - 1] = tz;
if (fabs((tx = 0.5*(X[FE[Current][2] -
1] + X[FE[Current][1] - 1])) - X[FE[Current][5] - 1]) > Eps)
X[FE[Current][5] - 1] = tx;
if (fabs((ty = 0.5*(Y[FE[Current][2] -
1] + Y[FE[Current][1] - 1])) - Y[FE[Current][5] - 1]) > Eps)
Y[FE[Current][5] - 1] = ty;
if (fabs((tz = 0.5*(Z[FE[Current][2] -
1] + Z[FE[Current][1] - 1])) - Z[FE[Current][5] - 1]) > Eps)
Z[FE[Current][5] - 1] = tz;
if (fabs((tx = 0.5*(X[FE[Current][0] -
1] + X[FE[Current][2] - 1])) - X[FE[Current][6] - 1]) > Eps)
X[FE[Current][6] - 1] = tx;
if (fabs((ty = 0.5*(Y[FE[Current][0] -
1] + Y[FE[Current][2] - 1])) - Y[FE[Current][6] - 1]) > Eps)
Y[FE[Current][6] - 1] = ty;
if (fabs((tz = 0.5*(Z[FE[Current][0] -
1] + Z[FE[Current][2] - 1])) - Z[FE[Current][6] - 1]) > Eps)
Z[FE[Current][6] - 1] = tz;
if (fabs((tx = 0.5*(X[FE[Current][0] -
1] + X[FE[Current][3] - 1])) - X[FE[Current][7] - 1]) > Eps)
X[FE[Current][7] - 1] = tx;
if (fabs((ty = 0.5*(Y[FE[Current][0] -
1] + Y[FE[Current][3] - 1])) - Y[FE[Current][7] - 1]) > Eps)
Y[FE[Current][7] - 1] = ty;
if (fabs((tz = 0.5*(Z[FE[Current][0] -
1] + Z[FE[Current][3] - 1])) - Z[FE[Current][7] - 1]) > Eps)
Z[FE[Current][7] - 1] = tz;
if (fabs((tx = 0.5*(X[FE[Current][3] -
1] + X[FE[Current][1] - 1])) - X[FE[Current][8] - 1]) > Eps)
X[FE[Current][8] - 1] = tx;
if (fabs((ty = 0.5*(Y[FE[Current][3] -
1] + Y[FE[Current][1] - 1])) - Y[FE[Current][8] - 1]) > Eps)
Y[FE[Current][8] - 1] = ty;
if (fabs((tz = 0.5*(Z[FE[Current][3] -
1] + Z[FE[Current][1] - 1])) - Z[FE[Current][8] - 1]) > Eps)
Z[FE[Current][8] - 1] = tz;
if (fabs((tx = 0.5*(X[FE[Current][3] -
1] + X[FE[Current][2] - 1])) - X[FE[Current][9] - 1]) > Eps)
X[FE[Current][9] - 1] = tx;
if (fabs((ty = 0.5*(Y[FE[Current][3] -
1] + Y[FE[Current][2] - 1])) - Y[FE[Current][9] - 1]) > Eps)
Y[FE[Current][9] - 1] = ty;
if (fabs((tz = 0.5*(Z[FE[Current][3] -
1] + Z[FE[Current][2] - 1])) - Z[FE[Current][9] - 1]) > Eps)
Z[FE[Current][9] - 1] = tz;
}
if (++Current == 999)
{
Matrix<DWORD> t = FE;
FE.ReSize(2 * FE.Size1(),10);
for (DWORD i = 0; i < Current;
i++)
for (DWORD j = 0; j < 10; j++)
FE[i][j] = t[i][j];
}
if (Current % 100 == 0)
printf("Line:
%ld\r",Current);
}
NumFE = Current;
for (DWORD i = 0; i < NumFE; i++)
for (DWORD j = 0; j < 10; j++)
FE[i][j]--;
printf(" \r");
return true;
}
bool
ReadCubeData(char* fn1,char*fn2,FILE* in1,FILE* in2,Vector<double>&
X,Vector<double>& Y,Vector<double>& Z,
Matrix<DWORD>&
FE,DWORD& NumPoints,DWORD& NumFE)
{
double tx,
ty,
tz;
char TextBuffer[1001];
DWORD Current = 0L,
tmp;
while (!feof(in1))
{
if (fgets(TextBuffer,1000,in1) == NULL)
{
if (feof(in1)) break;
printf("Unable read file
%s",fn1);
fclose(in1);
fclose(in2);
return false;
}
if (sscanf(TextBuffer,"%ld %lf %lf
%lf", &NumPoints,&tx,&ty,&tz) != 4) continue;
X[Current] = tx;
Y[Current] = ty;
Z[Current] = tz;
if (++Current == 999)
{
Vector<double> t1 = X,
t2 = Y,
t3 = Z;
X.ReSize(2 * X.Size());
Y.ReSize(2 * X.Size());
Z.ReSize(2 * X.Size());
for (DWORD i = 0; i < Current;
i++)
{
X[i] = t1[i];
Y[i] = t2[i];
Z[i] = t3[i];
}
}
if (Current % 100 == 0)
printf("Line:
%ld\r",Current);
}
fclose(in1);
printf(" \r");
NumPoints = Current;
Current = 0L;
while (!feof(in2))
{
if (fgets(TextBuffer,1000,in2) == NULL)
{
if (feof(in2)) break;
printf("Unable read file
%s",fn2);
fclose(in2);
return false;
}
if (sscanf(TextBuffer,"%d %d %d %d
%d %ld %ld %ld %ld %ld %ld %ld %ld",

&tmp,&tmp,&tmp,&tmp,&tmp,

&FE[Current][0],&FE[Current][1],&FE[Current][3],&FE[Current][2],

&FE[Current][4],&FE[Current][5],&FE[Current][7],&FE[Current][6])
!= 13) continue;
if (++Current == 999)
{
Matrix<DWORD> t = FE;
FE.ReSize(2 * FE.Size1(),10);
for (DWORD i = 0; i < Current;
i++)
for (DWORD j = 0; j < 10; j++)
FE[i][j] = t[i][j];
}
if (Current % 100 == 0)
printf("Line:
%ld\r",Current);}
NumFE = Current;
for (DWORD i = 0; i < NumFE; i++)
for (DWORD j = 0; j < 10; j++)
FE[i][j]--;
printf(" \r");
return true;}
2.

,
.
#include
"matrix.h"
class RVector
{
private:
Vector<double> Buffer;
public:
RVector(void) {}
~RVector() {}
RVector(DWORD Size) {
Buffer.ReSize(Size); }
RVector(RVector& right) { Buffer =
right.Buffer; }
RVector(Vector<double>& right) {
Buffer = right; }
DWORD Size(void) { return Buffer.Size(); }
void
ReSize(DWORD Size) { Buffer.ReSize(Size); }
double& operator [] (DWORD i) { return Buffer[i]; }
RVector& operator = (RVector&
right) { Buffer = right.Buffer; return *this; }
RVector& operator =
(Vector<double>& right) { Buffer = right; return *this; }
void Sub(RVector&);
void Sub(RVector&,double);
void Mul(double);
void Add(RVector&);
friend double
Norm(RVector&,RVector&);
};
class TSMatrix
{
private:
Vector<double> Right;
Vector<double>* Array;
Vector<DWORD>* Links;
uint Dim;
DWORD Size;
public:
TSMatrix(void) { Size = 0; Dim = 0; Array =
NULL; Links = NULL; }
TSMatrix(Vector<DWORD>*,DWORD,uint);
~TSMatrix(void) { if (Array) delete []
Array; }
Vector<double>& GetRight(void) {
return Right; }
DWORD GetSize(void) { return Size; }
uint
GetDim(void) { return Dim; }
Vector<double>& GetVector(DWORD
i) { return Array[i]; }
Vector<DWORD>* GetLinks(void) { return Links; }
void
SetLinks(Vector<DWORD>* l) { Links = l; }
void
Add(Matrix<double>&,Vector<DWORD>&);
void Add(DWORD I, DWORD L, DWORD J, DWORD
K, double v)
{
DWORD Row = I,
Col = L * Links[I].Size() * Dim +
Find(I,J) * Dim + K;
Array[Row][Col] += v;
}
void Add(DWORD I, double v)
{
Right[I] += v;
}
DWORD Find(DWORD,DWORD);
void
Restore(Matrix<double>&);
void
Set(DWORD,DWORD,double,bool);
void
Set(DWORD Index1,DWORD Index2,double value)
{
DWORD I
= Index1 / Dim,
L
= Index1 % Dim,
J
= Index2 / Dim,
K
= Index2 % Dim,
Pos = Find(I,J),
Row = I,
Col;
if (Pos == DWORD(-1)) return;
Col = L * Links[I].Size() * Dim +
Find(I,J) * Dim + K;
Array[Row][Col] = value;
}
bool
Get(DWORD Index1,DWORD Index2,double& value)
{
DWORD I
= Index1 / Dim,
L
= Index1 % Dim,
J
= Index2 / Dim,
K
= Index2 % Dim,
Pos = Find(I,J),
Row = I,
Col;
value = 0;
if (Pos == DWORD(-1)) return false;

Col = L * Links[I].Size() * Dim + Find(I,J) * Dim + K;
value = Array[Row][Col];
return true;
}
void
Mul(RVector&,RVector&);
double
Mul(DWORD,RVector&);
void
write(ofstream&);
void
read(ifstream&);
};
class RMatrix
{
private:
Vector<double> Buffer;
DWORD size;
public:
RMatrix(DWORD sz) { size = sz;
Buffer.ReSize(size*(size + 1)*0.5); }
~RMatrix() {}
DWORD Size(void) { return size; }
double& Get(DWORD i,DWORD j) { return
Buffer[(2*size + 1 - i)*0.5*i + j - i]; }
};
//************************
#include "smatrix.h"
double Norm(RVector& Left,RVector& Right)
{
double
Ret = 0;
for
(DWORD i = 0; i < Left.Size(); i++)
Ret +=
Left[i] * Right[i];
return
Ret;
}
void RVector::Sub(RVector& Right)
{
for
(DWORD i = 0; i < Size(); i++)

(*this)[i] -= Right[i];
}
void RVector::Add(RVector& Right)
{
for
(DWORD i = 0; i < Size(); i++)

(*this)[i] += Right[i];
}
void RVector::Mul(double koff)
{
for
(DWORD i = 0; i < Size(); i++)
(*this)[i]
*= koff;
}
void RVector::Sub(RVector& Right,double koff)
{
for
(DWORD i = 0; i < Size(); i++)

(*this)[i] -= Right[i]*koff;
}
TSMatrix::TSMatrix(Vector<DWORD>* links,
DWORD size, uint dim)
{
Dim = dim;

Links = links;
Size = size;

Right.ReSize(Dim * Size);
Array =
new Vector<double>[Size];
for
(DWORD i = 0; i < Size; i++)

Array[i].ReSize(Links[i].Size() * Dim * Dim);
}
void TSMatrix::Add(Matrix<double>&
FEMatr,Vector<DWORD>& FE)
{
double
Res;

DWORD RRow;
for (DWORD
i = 0L; i < FE.Size(); i++) for (DWORD l
= 0L; l < Dim; l++) for (DWORD j = 0L; j < FE.Size(); j++) for (DWORD k = 0L; k < Dim; k++)
{
Res =
FEMatr[i * Dim + l][j * Dim + k];
if
(Res) Add(FE[i],l,FE[j],k,Res);
}
for
(DWORD i = 0L; i < FE.Size(); i++) for (DWORD l
= 0L; l < Dim; l++) {
RRow =
FE[UINT(i % (FE.Size()))] * Dim + l;
Res =
FEMatr[i * Dim + l][FEMatr.Size1()];
if
(Res) Add(RRow,Res); }
}
DWORD TSMatrix::Find(DWORD I,DWORD J)
{
DWORD i;
for (i =
0; i < Links[I].Size(); i++)
if
(Links[I][i] == J) return i;
return
DWORD(-1);
}
void TSMatrix::Restore(Matrix<double>&
Matr)
{
DWORD i,
j,

NRow,

NPoint,

NLink,

Pos;

Matr.ReSize(Size * Dim,Size * Dim + 1);
for (i =
0; i < Size; i++)
for (j =
0; j < Array[i].Size(); j++)
{

NRow = j / (Array[i].Size() /
Dim); // Number of row

NPoint = (j - NRow *
(Array[i].Size() / Dim)) / Dim; // Number of points

NLink = j % Dim; // Number
of link

Pos = Links[i][NPoint];

Matr[i * Dim + NRow][Pos * Dim + NLink] = Array[i][j];
}
for (i =
0; i < Right.Size(); i++) Matr[i][Matr.Size1()] = Right[i];
}
void
TSMatrix::Set(DWORD Index,DWORD Position,double Value,bool Case)
{

DWORD Row = Index,

Col = Position *
Links[Index].Size() * Dim + Find(Index,Index) * Dim + Position,
i;
double
koff = Array[Row][Col],

val;
if
(!Case)

Right[Dim * Index + Position] = Value;
else
{

Right[Index * Dim + Position] = Value * koff;
for (i
= 0L; i < Size * Dim; i++)
if (i
!= Index * Dim + Position)
{

Set(Index * Dim + Position,i,0);

Set(i,Index * Dim + Position,0);
if (Get(i,Index * Dim + Position,val))

Right[i] -= val * Value;
}
}
}
void
TSMatrix::Mul(RVector& Arr,RVector& Res)
{
DWORD i,
j,

NRow,

NPoint,

NLink,

Pos;

Res.ReSize(Arr.Size());
for (i =
0; i < Size; i++)
for (j =
0; j < Array[i].Size(); j++)
{

NRow = j / (Array[i].Size() /
Dim);

NPoint = (j - NRow *
(Array[i].Size() / Dim)) / Dim;

NLink = j % Dim;

Pos = Links[i][NPoint];
Res[i
* Dim + NRow] += Arr[Pos * Dim + NLink] * Array[i][j];
}
}
double TSMatrix::Mul(DWORD Index,RVector& Arr)
{

DWORD j,

I = Index / Dim,

L = Index % Dim,

Start = L * (Array[I].Size() / Dim),

Stop = Start + (Array[I].Size()
/ Dim),

NRow,

NPoint,

NLink,

Pos;
double
Res = 0;
for (j =
Start; j < Stop; j++)
{

NRow = j / (Array[I].Size() /
Dim);

NPoint = (j - NRow *
(Array[I].Size() / Dim)) / Dim;

NLink = j % Dim;

Pos = Links[I][NPoint];
Res +=
Arr[Pos * Dim + NLink] * Array[I][j];
}
return
Res;
}
void
TSMatrix::write(ofstream& Out)
{
DWORD
ColSize;

Out.write((char*)&(Dim),sizeof(DWORD));

Out.write((char*)&(Size),sizeof(DWORD));
for
(DWORD i = 0; i < Size; i++)
{

ColSize = Array[i].Size();

Out.write((char*)&(ColSize),sizeof(DWORD));
for
(DWORD j = 0; j < ColSize; j++)

Out.write((char*)&(Array[i][j]),sizeof(double));
}
for
(DWORD i = 0; i < Size * Dim; i++)

Out.write((char*)&(Right[i]),sizeof(double));
}
void TSMatrix::read(ifstream& In)
{
DWORD
ColSize;

In.read((char*)&(Dim),sizeof(DWORD));

In.read((char*)&(Size),sizeof(DWORD));
if
(Array) delete [] Array;
Array =
new Vector<double>[Size];

Right.ReSize(Size * Dim);
for
(DWORD i = 0; i < Size; i++)
{

In.read((char*)&(ColSize),sizeof(DWORD));

Array[i].ReSize(ColSize);
for
(DWORD j = 0; j < ColSize; j++)

In.read((char*)&(Array[i][j]),sizeof(double));
}
for (DWORD
i = 0; i < Size * Dim; i++)

In.read((char*)&(Right[i]),sizeof(double));
}

      ©2010