Смекни!
smekni.com
Математика

Алгоритм компактного хранения и решения СЛАУ высокого порядка (стр. 3 из 5)

Предложеннаяв работе методика компактного хранения матриц коэффициентов СЛАУ ииспользования метода Ланцоша позволили на примере решения контактных задачдобиться существенной экономии процессорного времени и затрат оперативнойпамяти.

ПРИЛОЖЕНИЕ 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]&bsol;n");

  printf("Switch: -t10 -Tetraedr(10)&bsol;n");

  printf("        -c8  - Cube(8)&bsol;n");

  printf("        -s4  -Shell(4)&bsol;n");

  printf("        -s8  -Shell(8)&bsol;n&bsol;n");

  printf("Optins: /8 - convertTetraedr(10)->8*Tetraedr(4)&bsol;n");

  printf("        /6 - convertCube(8)->6*Tetraedr(4)&bsol;n");

}

bool Output(char*fname,Vector<double>& x,Vector<double>&y,Vector<double>& z, Matrix<DWORD>& tr, DWORD n,

            DWORD NumNewPoints,DWORDntr,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((constchar*)&type,sizeof(int));

  if (Out.fail()) return true;

  Out.write((constchar*)&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((constchar*)&(NumNewPoints),sizeof(DWORD));

  if (Out.fail())

   {

     Out.close();

     return true;

   }

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

   {

     Out.write((constchar*)&x[i],sizeof(double));

     Out.write((constchar*)&y[i],sizeof(double));

     Out.write((constchar*)&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((constchar*)&(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((constchar*)&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((constchar*)&out,sizeof(DWORD));

      if (Out.fail())

        {

          Out.close();

          return true;

        }

    }

{

  //*********************

  // Create Links

  printf("Create links...&bsol;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((constchar*)&Sz,sizeof(DWORD));

     for (DWORD j = 0; j < Sz; j++)

      Out.write((constchar*)&(Links[i][j]),sizeof(DWORD));

   }

  //*********************

}

  printf("                           &bsol;r");

  printf("Points: %ld&bsol;n",n);

  printf("FE:     %ld&bsol;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...&bsol;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 fileconvertor. ZSU(c) 1998.&bsol;n&bsol;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&bsol;n&bsol;n",sw);

      PrintHeader();

      return;

    }

   if (argc == 6)

    {

      op = argv[5];

      if (strcmp(op,"/8")&& strcmp(op,"/6"))

       {

         printf("Unknown options%s&bsol;n&bsol;n",op);

         PrintHeader();

         return;

       }

    }

   if (CreateFile(input1,input2,input3,op))

    printf("OK&bsol;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"))

        {

 
назад
читать дальше
1
2
3
4
5
НАПИШИТЕ НАМ
Copyright © Smekni.com. All rigths reserved.