cell_tet10.c

一个用来实现偏微分方程中网格的计算库

    case VTK:
      {
        conn.resize(10);
        conn[0] = this->node(0);
        conn[1] = this->node(1);
        conn[2] = this->node(2);
        conn[3] = this->node(3);
        conn[4] = this->node(4);
        conn[5] = this->node(5);
        conn[6] = this->node(6);
        conn[7] = this->node(7);
        conn[8] = this->node(8);
        conn[9] = this->node(9);
        return;
        /*
           conn.resize(4);
           switch (sc)
           {            
        // Linear sub-tet 0
        case 0:

        conn[0] = this->node(0);
        conn[1] = this->node(4);
        conn[2] = this->node(6);
        conn[3] = this->node(7);

        return;

        // Linear sub-tet 1
        case 1:

        conn[0] = this->node(4);
        conn[1] = this->node(1);
        conn[2] = this->node(5);
        conn[3] = this->node(8);

        return;

        // Linear sub-tet 2
        case 2:

        conn[0] = this->node(5);
        conn[1] = this->node(2);
        conn[2] = this->node(6);
        conn[3] = this->node(9);

        return;

        // Linear sub-tet 3
        case 3:

        conn[0] = this->node(7);
        conn[1] = this->node(8);
        conn[2] = this->node(9);
        conn[3] = this->node(3);

        return;

        // Linear sub-tet 4
        case 4:

        conn[0] = this->node(4);
        conn[1] = this->node(8);
        conn[2] = this->node(6);
        conn[3] = this->node(7);

        return;

        // Linear sub-tet 5
        case 5:

        conn[0] = this->node(4);
        conn[1] = this->node(5);
        conn[2] = this->node(6);
        conn[3] = this->node(8);

        return;

        // Linear sub-tet 6
        case 6:

        conn[0] = this->node(5);
        conn[1] = this->node(9);
        conn[2] = this->node(6);
        conn[3] = this->node(8);

        return;

        // Linear sub-tet 7
        case 7:

        conn[0] = this->node(7);
        conn[1] = this->node(6);
        conn[2] = this->node(9);
        conn[3] = this->node(8);

        return;


        default:

        libmesh_error();
      }
      */
      }

    default:
      libmesh_error();
  }

  libmesh_error();
}



const unsigned short int Tet10::_second_order_vertex_child_number[10] =
{
  99,99,99,99, // Vertices
  0,1,0,0,1,2  // Edges
};



const unsigned short int Tet10::_second_order_vertex_child_index[10] =
{
  99,99,99,99, // Vertices
  1,2,2,3,3,3  // Edges
};



std::pair<unsigned short int, unsigned short int>
Tet10::second_order_child_vertex (const unsigned int n) const
{
  libmesh_assert (n >= this->n_vertices());
  libmesh_assert (n < this->n_nodes());
  return std::pair<unsigned short int, unsigned short int>
    (_second_order_vertex_child_number[n],
     _second_order_vertex_child_index[n]);
}



unsigned short int Tet10::second_order_adjacent_vertex (const unsigned int n,
							const unsigned int v) const
{ 
  libmesh_assert (n >= this->n_vertices());
  libmesh_assert (n <  this->n_nodes());
  libmesh_assert (v < 2);
  return _second_order_adjacent_vertices[n-this->n_vertices()][v]; 
}



const unsigned short int Tet10::_second_order_adjacent_vertices[6][2] = 
{
  {0, 1}, // vertices adjacent to node 4 
  {1, 2}, // vertices adjacent to node 5 
  {0, 2}, // vertices adjacent to node 6 
  {0, 3}, // vertices adjacent to node 7 
  {1, 3}, // vertices adjacent to node 8 
  {2, 3}  // vertices adjacent to node 9 
};





#ifdef ENABLE_AMR

const float Tet10::_embedding_matrix[8][10][10] =
{
  // embedding matrix for child 0
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    1.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.}, // 0
    {    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.,    0.}, // 1
    {    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.}, // 2
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.}, // 3
    { 0.375,-0.125,    0.,    0.,  0.75,    0.,    0.,    0.,    0.,    0.}, // 4
    {    0.,-0.125,-0.125,    0.,   0.5,  0.25,   0.5,    0.,    0.,    0.}, // 5
    { 0.375,    0.,-0.125,    0.,    0.,    0.,  0.75,    0.,    0.,    0.}, // 6
    { 0.375,    0.,    0.,-0.125,    0.,    0.,    0.,  0.75,    0.,    0.}, // 7
    {    0.,-0.125,    0.,-0.125,   0.5,    0.,    0.,   0.5,  0.25,    0.}, // 8
    {    0.,    0.,-0.125,-0.125,    0.,    0.,   0.5,   0.5,    0.,  0.25}  // 9
  },

  // embedding matrix for child 1
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.,    0.}, // 0
    {    0.,    1.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.}, // 1
    {    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.}, // 2
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.}, // 3
    {-0.125, 0.375,    0.,    0.,  0.75,    0.,    0.,    0.,    0.,    0.}, // 4
    {    0., 0.375,-0.125,    0.,    0.,  0.75,    0.,    0.,    0.,    0.}, // 5
    {-0.125,    0.,-0.125,    0.,   0.5,   0.5,  0.25,    0.,    0.,    0.}, // 6
    {-0.125,    0.,    0.,-0.125,   0.5,    0.,    0.,  0.25,   0.5,    0.}, // 7
    {    0., 0.375,    0.,-0.125,    0.,    0.,    0.,    0.,  0.75,    0.}, // 8
    {    0.,    0.,-0.125,-0.125,    0.,   0.5,    0.,    0.,   0.5,  0.25}  // 9
  },

  // embedding matrix for child 2
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.}, // 0
    {    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.}, // 1
    {    0.,    0.,    1.,    0.,    0.,    0.,    0.,    0.,    0.,    0.}, // 2
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.}, // 3
    {-0.125,-0.125,    0.,    0.,  0.25,   0.5,   0.5,    0.,    0.,    0.}, // 4
    {    0.,-0.125, 0.375,    0.,    0.,  0.75,    0.,    0.,    0.,    0.}, // 5
    {-0.125,    0., 0.375,    0.,    0.,    0.,  0.75,    0.,    0.,    0.}, // 6
    {-0.125,    0.,    0.,-0.125,    0.,    0.,   0.5,  0.25,    0.,   0.5}, // 7
    {    0.,-0.125,    0.,-0.125,    0.,   0.5,    0.,    0.,  0.25,   0.5}, // 8
    {    0.,    0., 0.375,-0.125,    0.,    0.,    0.,    0.,    0.,  0.75}  // 9
  },

  // embedding matrix for child 3
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.}, // 0
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.}, // 1
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.}, // 2
    {    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.,    0.,    0.}, // 3
    {-0.125,-0.125,    0.,    0.,  0.25,    0.,    0.,   0.5,   0.5,    0.}, // 4
    {    0.,-0.125,-0.125,    0.,    0.,  0.25,    0.,    0.,   0.5,   0.5}, // 5
    {-0.125,    0.,-0.125,    0.,    0.,    0.,  0.25,   0.5,    0.,   0.5}, // 6
    {-0.125,    0.,    0., 0.375,    0.,    0.,    0.,  0.75,    0.,    0.}, // 7
    {    0.,-0.125,    0., 0.375,    0.,    0.,    0.,    0.,  0.75,    0.}, // 8
    {    0.,    0.,-0.125, 0.375,    0.,    0.,    0.,    0.,    0.,  0.75}  // 9
  },

  // embedding matrix for child 4
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.,    0.}, // 0
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.}, // 1
    {    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.}, // 2
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.}, // 3
    {-0.125,    0.,    0.,-0.125,   0.5,    0.,    0.,  0.25,   0.5,    0.}, // 4
    {-0.125,-0.125,-0.125,-0.125,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25}, // 5
    {    0.,-0.125,-0.125,    0.,   0.5,  0.25,   0.5,    0.,    0.,    0.}, // 6
    {    0.,-0.125,    0.,-0.125,   0.5,    0.,    0.,   0.5,  0.25,    0.}, // 7
    {-0.125,-0.125,    0.,    0.,  0.25,    0.,    0.,   0.5,   0.5,    0.}, // 8
    {    0.,    0.,-0.125,-0.125,    0.,    0.,   0.5,   0.5,    0.,  0.25}  // 9
  },

  // embedding matrix for child 5
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.,    0.}, // 0
    {    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.}, // 1
    {    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.}, // 2
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.}, // 3
    {-0.125,    0.,-0.125,    0.,   0.5,   0.5,  0.25,    0.,    0.,    0.}, // 4
    {-0.125,-0.125,    0.,    0.,  0.25,   0.5,   0.5,    0.,    0.,    0.}, // 5
    {    0.,-0.125,-0.125,    0.,   0.5,  0.25,   0.5,    0.,    0.,    0.}, // 6
    {-0.125,    0.,    0.,-0.125,   0.5,    0.,    0.,  0.25,   0.5,    0.}, // 7
    {    0.,    0.,-0.125,-0.125,    0.,   0.5,    0.,    0.,   0.5,  0.25}, // 8
    {-0.125,-0.125,-0.125,-0.125,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25}  // 9
  },

  // embedding matrix for child 6
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.}, // 0
    {    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.,    0.}, // 1
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.}, // 2
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.}, // 3
    {-0.125,-0.125,    0.,    0.,  0.25,   0.5,   0.5,    0.,    0.,    0.}, // 4
    {    0.,-0.125,    0.,-0.125,    0.,   0.5,    0.,    0.,  0.25,   0.5}, // 5
    {-0.125,    0.,    0.,-0.125,    0.,    0.,   0.5,  0.25,    0.,   0.5}, // 6
    {-0.125,-0.125,-0.125,-0.125,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25}, // 7
    {    0.,    0.,-0.125,-0.125,    0.,   0.5,    0.,    0.,   0.5,  0.25}, // 8
    {    0.,-0.125,-0.125,    0.,    0.,  0.25,    0.,    0.,   0.5,   0.5}  // 9
  },

  // embedding matrix for child 7
  {
    //    0      1      2      3      4      5      6      7      8      9
    {    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.,    0.}, // 0
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.}, // 1
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.}, // 2
    {    0.,    0.,    0.,    0.,    0.,    0.,    0.,    1.,    0.,    0.}, // 3
    {-0.125,-0.125,-0.125,-0.125,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25}, // 4
    {    0.,-0.125,-0.125,    0.,    0.,  0.25,    0.,    0.,   0.5,   0.5}, // 5
    {-0.125,    0.,    0.,-0.125,    0.,    0.,   0.5,  0.25,    0.,   0.5}, // 6
    {    0.,    0.,-0.125,-0.125,    0.,    0.,   0.5,   0.5,    0.,  0.25}, // 7
    {-0.125,-0.125,    0.,    0.,  0.25,    0.,    0.,   0.5,   0.5,    0.}, // 8
    {-0.125,    0.,-0.125,    0.,    0.,    0.,  0.25,   0.5,    0.,   0.5}  // 9
  }
};






float Tet10::embedding_matrix (const unsigned int i,
			      const unsigned int j,
			      const unsigned int k) const
{
  // Check for uninitialized diagonal selection
  if (this->_diagonal_selection==INVALID_DIAG)
    {
      Real diag_01_23 = (this->point(0)+this->point(1)-this->point(2)-this->point(3)).size_sq();
      Real diag_02_13 = (this->point(0)-this->point(1)+this->point(2)-this->point(3)).size_sq();
      Real diag_03_12 = (this->point(0)-this->point(1)-this->point(2)+this->point(3)).size_sq();

      this->_diagonal_selection=DIAG_02_13;

      if (diag_01_23 < diag_02_13 || diag_03_12 < diag_02_13)
      {
        if (diag_01_23 < diag_03_12)
          this->_diagonal_selection=DIAG_01_23;

        else
          this->_diagonal_selection=DIAG_03_12;
      }
    }

  // Permuted j and k indices
  unsigned int
    jp=j,
    kp=k;

  if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
    {
      // Permute j, k
      if (jp<3) jp=(jp+static_cast<unsigned int>(this->_diagonal_selection))%3;
      else if (jp>3) jp=(jp-1+static_cast<unsigned int>(this->_diagonal_selection))%3+1+3*((jp-1)/3);
      if (kp<3) kp=(kp+static_cast<unsigned int>(this->_diagonal_selection))%3;
      else if (kp>3) kp=(kp-1+static_cast<unsigned int>(this->_diagonal_selection))%3+1+3*((kp-1)/3);
    }

  // Call embedding matrx with permuted indices
  return this->_embedding_matrix[i][jp][kp];
}

#endif // #ifdef ENABLE_AMR