latman.c

基于格子Boltzmann方法开源可视化软件的源代码 具有很高的实用价值。对学习LBM方法及其软件开发非常游泳

//##############################################################################
//
// latman.c
//
//  - Lattice Manager.
//
//  - Routines for managing a lattice:
//
//    - construct
//    - init
//    - destruct
//

// void process_matrix( struct lattice_struct *lattice, int **matrix)
//##############################################################################
//
// P R O C E S S   M A T R I X 
//
//  - Process the soil matrix.
//
//  - Convert the full matrix representation into sparse lattice form.
//
//  - This routine is useful when the domain is read from a BMP file
//    representing the full lattice.  If more general mechanisms
//    are developed (e.g., reading parameterized fracture information),
//    then determining the sparse lattice data structure will require a
//    corresponding mechanism (and storing the full lattice, even for
//    pre-/post-processing, may become undesirable if not impossible).
//
void process_matrix( struct lattice_struct *lattice, int **matrix)
{
  // Variable declarations.
  int i,  j;
  int in, jn;
  int ip, jp;
  int ei, ej;
  int NumNodes;
  int n;

  // Ending indices.
  ei = lattice->param.LX-1;
  ej = lattice->param.LY-1;

  // Mark solid nodes that have fluid nodes as neighbors and
  // count the total number of nodes requiring storage.
  NumNodes = 0;
  if( lattice->periodic_y)
  {
    for( j=0; j<=ej; j++)
    {
      jp = ( j<ej) ? ( j+1) : ( 0 );
      jn = ( j>0 ) ? ( j-1) : ( ej);

      for( i=0; i<=ei; i++)
      {
        ip = ( i<ei) ? ( i+1) : ( 0 );
        in = ( i>0 ) ? ( i-1) : ( ei);

        if( matrix[j ][i ] == 0) 
        { 
          NumNodes++;

          if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumNodes++;}
          if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumNodes++;}
          if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumNodes++;}
          if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumNodes++;}
          if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumNodes++;}
          if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumNodes++;}
          if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumNodes++;}
          if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumNodes++;}

        }

      } /* for( i=0; i<=ei; i++) */
    } /* for( j=0; j<=ej; j++) */
  } /* if( lattice->periodic_y) */

  else /* !lattice->periodic_y */
  {
    j  = 0;
    jp = 1;
    for( i=0; i<=ei; i++)
    {
      ip = ( i<ei) ? ( i+1) : ( 0 );
      in = ( i>0 ) ? ( i-1) : ( ei);

      if( matrix[j ][i ] == 0) 
      { 
        NumNodes++;

        if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumNodes++;}
        if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumNodes++;}
        if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumNodes++;}
        //if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumNodes++;}
        if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumNodes++;}
        if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumNodes++;}
        //if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumNodes++;}
        //if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumNodes++;}

      }

    } /* for( i=0; i<=ei; i++) */


    for( j=1; j<ej; j++)
    {
      jp = j+1;
      jn = j-1;

      for( i=0; i<=ei; i++)
      {
        ip = ( i<ei) ? ( i+1) : ( 0 );
        in = ( i>0 ) ? ( i-1) : ( ei);

        if( matrix[j ][i ] == 0) 
        { 
          NumNodes++;

          if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumNodes++;}
          if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumNodes++;}
          if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumNodes++;}
          if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumNodes++;}
          if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumNodes++;}
          if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumNodes++;}
          if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumNodes++;}
          if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumNodes++;}

        }

      } /* for( i=0; i<=ei; i++) */
    } /* for( j=0; j<=ej; j++) */

    j  = ej;
    jn = ej-1;
    for( i=0; i<=ei; i++)
    {
      ip = ( i<ei) ? ( i+1) : ( 0 );
      in = ( i>0 ) ? ( i-1) : ( ei);

      if( matrix[j ][i ] == 0) 
      { 
        NumNodes++;

        if( matrix[j ][ip] == 1) { matrix[j ][ip] = 2; NumNodes++;}
        //if( matrix[jp][i ] == 1) { matrix[jp][i ] = 2; NumNodes++;}
        if( matrix[j ][in] == 1) { matrix[j ][in] = 2; NumNodes++;}
        if( matrix[jn][i ] == 1) { matrix[jn][i ] = 2; NumNodes++;}
        //if( matrix[jp][ip] == 1) { matrix[jp][ip] = 2; NumNodes++;}
        //if( matrix[jp][in] == 1) { matrix[jp][in] = 2; NumNodes++;}
        if( matrix[jn][in] == 1) { matrix[jn][in] = 2; NumNodes++;}
        if( matrix[jn][ip] == 1) { matrix[jn][ip] = 2; NumNodes++;}

      }

    } /* for( i=0; i<=ei; i++) */

  } /* if( lattice->periodic_y) else */

#if 0 // Dump the matrix contents to the screen.
  for( j=0; j<=ej; j++)
  {
    for( i=0; i<=ei; i++)
    {
      printf(" %d", matrix[j][i]);
    }
    printf("\n");
  }
#endif

  // Set NumNodes in the lattice.
  lattice->NumNodes = NumNodes;

#if VERBOSITY_LEVEL > 0
  printf("process_matrix() -- NumNodes = %d\n", NumNodes);
#endif /* VERBOSITY_LEVEL > 0 */

  // Allocate memory for NumNodes nodes.
  lattice->node = 
    ( struct node_struct*)malloc( NumNodes*sizeof( struct node_struct));
  assert( lattice->node!=NULL);

  // Allocate memory for NumNodes boundary conditions.
  lattice->bc= 
    ( struct bc_struct*)malloc( NumNodes*sizeof( struct bc_struct));
  assert( lattice->bc!=NULL);

  // Set coordinates and index of lattice nodes.
  // Use matrix entries to store pointers to associated nodes to
  // facilitate finding neighbors on a following traversal.
  n = 0;
  for( j=0; j<=ej; j++)
  {
    for( i=0; i<=ei; i++)
    {
      if( matrix[j][i] != 1)
      {
        // Coordinates.
        lattice->node[n].i = i;
        lattice->node[n].j = j;

        // Index.
        lattice->node[n].n = n;

        // Distinguish solids.
        if( matrix[j][i] == 2) { lattice->bc[n].bc_type |= BC_SOLID_NODE;}
        else { lattice->bc[n].bc_type = 0;}

        // Now use matrix entry to point to the corresponding lattice node.
        // This will facilitate assigning neighbors below.
        matrix[j][i] = (int)&(lattice->node[n]);

//printf("(int)&(lattice->node[%d]) = %d\n", n, (int)&(lattice->node[n]));
//printf("((node_ptr)matrix[%d][%d])->n = %d\n", 
//    j, i, ((node_ptr)matrix[j ][i ])->n);

        n++;

      } /* if( matrix[j][i] != 1) */
    } /* for( i=0; i<=ei; i++) */
  } /* for( j=0; j<=ej; j++) */

#if 0 // Dump the matrix contents to the screen.
  for( j=0; j<=ej; j++)
  {
    for( i=0; i<=ei; i++)
    {
      printf(" %d", matrix[j][i]);
    }
    printf("\n");
  }
#endif

  // Assign neighbors.
  if( lattice->periodic_y)
  {
    n = 0;
    for( j=0; j<=ej; j++)
    {
      jp = ( j<ej) ? ( j+1) : ( 0 );
      jn = ( j>0 ) ? ( j-1) : ( ej);

      for( i=0; i<=ei; i++)
      {
        ip = ( i<ei) ? ( i+1) : ( 0 );
        in = ( i>0 ) ? ( i-1) : ( ei);
         
        if( matrix[j ][i ] != 1)
        {

          if( matrix[j ][ip] != 1) 
          { lattice->node[n].nn[0] = ((node_ptr)matrix[j ][ip])->n;}
          else
          { lattice->node[n].nn[0] = -1;}

          if( matrix[jp][i ] != 1) 
          { lattice->node[n].nn[1] = ((node_ptr)matrix[jp][i ])->n;}
          else
          { lattice->node[n].nn[1] = -1;}

          if( matrix[j ][in] != 1) 
          { lattice->node[n].nn[2] = ((node_ptr)matrix[j ][in])->n;}
          else
          { lattice->node[n].nn[2] = -1;}

          if( matrix[jn][i ] != 1) 
          { lattice->node[n].nn[3] = ((node_ptr)matrix[jn][i ])->n;}
          else
          { lattice->node[n].nn[3] = -1;}

          if( matrix[jp][ip] != 1) 
          { lattice->node[n].nn[4] = ((node_ptr)matrix[jp][ip])->n;}
          else
          { lattice->node[n].nn[4] = -1;}

          if( matrix[jp][in] != 1) 
          { lattice->node[n].nn[5] = ((node_ptr)matrix[jp][in])->n;}
          else
          { lattice->node[n].nn[5] = -1;}

          if( matrix[jn][in] != 1) 
          { lattice->node[n].nn[6] = ((node_ptr)matrix[jn][in])->n;}
          else
          { lattice->node[n].nn[6] = -1;}

          if( matrix[jn][ip] != 1) 
          { lattice->node[n].nn[7] = ((node_ptr)matrix[jn][ip])->n;}
          else
          { lattice->node[n].nn[7] = -1;}

          n++;

        } /* if( matrix[j][i] != 1) */
      } /* for( i=0; i<=ei; i++) */
    } /* for( j=0; j<=ej; j++) */
  } /* if( lattice->periodic_y) */

  else /* !lattice->periodic_y */
  {
    n = 0;

    j  = 0;
    jp = 1;
    for( i=0; i<=ei; i++)
    {
      ip = ( i<ei) ? ( i+1) : ( 0 );
      in = ( i>0 ) ? ( i-1) : ( ei);
       
      if( matrix[j ][i ] != 1)
      {

        if( matrix[j ][ip] != 1) 
        { lattice->node[n].nn[0] = ((node_ptr)matrix[j ][ip])->n;}
        else
        { lattice->node[n].nn[0] = -1;}

        if( matrix[jp][i ] != 1) 
        { lattice->node[n].nn[1] = ((node_ptr)matrix[jp][i ])->n;}
        else
        { lattice->node[n].nn[1] = -1;}

        if( matrix[j ][in] != 1) 
        { lattice->node[n].nn[2] = ((node_ptr)matrix[j ][in])->n;}
        else
        { lattice->node[n].nn[2] = -1;}

        //if( matrix[jn][i ] != 1) 
        //{ lattice->node[n].nn[3] = ((node_ptr)matrix[jn][i ])->n;}
        //else
        //{ lattice->node[n].nn[3] = -1;}
        lattice->node[n].nn[3] = -1;

        if( matrix[jp][ip] != 1) 
        { lattice->node[n].nn[4] = ((node_ptr)matrix[jp][ip])->n;}
        else
        { lattice->node[n].nn[4] = -1;}

        if( matrix[jp][in] != 1) 
        { lattice->node[n].nn[5] = ((node_ptr)matrix[jp][in])->n;}
        else
        { lattice->node[n].nn[5] = -1;}

        //if( matrix[jn][in] != 1) 
        //{ lattice->node[n].nn[6] = ((node_ptr)matrix[jn][in])->n;}
        //else
        //{ lattice->node[n].nn[6] = -1;}
        lattice->node[n].nn[6] = -1;

        //if( matrix[jn][ip] != 1) 
        //{ lattice->node[n].nn[7] = ((node_ptr)matrix[jn][ip])->n;}
        //else
        //{ lattice->node[n].nn[7] = -1;}
        lattice->node[n].nn[7] = -1;

        n++;

      } /* if( matrix[j][i] != 1) */
    } /* for( i=0; i<=ei; i++) */

    for( j=1; j<ej; j++)
    {
      jp =j+1;
      jn =j-1;

      for( i=0; i<=ei; i++)