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));
}
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;
}
}
}
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; 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; 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; 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;
}
}
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;
if (argc<5 || argc>6) 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);
return;
}
if (argc==6)
{
op=argv[5];
if (strcmp(op, "/8") && strcmp(op, "/6"))
{
printf("Unknown options %s\n\n", op);
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;
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")) 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")) 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);
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;
}
}
}
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);
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]--;
return true;
}