C++ NodeSet::have_any_mid_node方法代码示例

void NodeSetTest::test_have_any()
{
  NodeSet set;
  set.clear();
  
  CPPUNIT_ASSERT( !set.have_any_corner_node() );
  set.set_corner_node(0);
  CPPUNIT_ASSERT( set.have_any_corner_node() );
  set.clear_corner_node(0);
  CPPUNIT_ASSERT( !set.have_any_corner_node() );
  set.set_corner_node( NodeSet::NUM_CORNER_BITS-1 );
  CPPUNIT_ASSERT( set.have_any_corner_node() );
  CPPUNIT_ASSERT( !set.have_any_mid_node() );

  CPPUNIT_ASSERT( !set.have_any_mid_edge_node() );
  set.set_mid_edge_node(0);
  CPPUNIT_ASSERT( set.have_any_mid_edge_node() );
  CPPUNIT_ASSERT( set.have_any_mid_node() );
  set.clear_mid_edge_node(0);
  CPPUNIT_ASSERT( !set.have_any_mid_edge_node() );
  CPPUNIT_ASSERT( !set.have_any_mid_node() );
  set.set_mid_edge_node( NodeSet::NUM_EDGE_BITS-1 );
  CPPUNIT_ASSERT( set.have_any_mid_edge_node() );
  CPPUNIT_ASSERT( set.have_any_mid_node() );
  set.clear();

  CPPUNIT_ASSERT( !set.have_any_mid_face_node() );
  set.set_mid_face_node(0);
  CPPUNIT_ASSERT( set.have_any_mid_face_node() );
  CPPUNIT_ASSERT( set.have_any_mid_node() );
  set.clear_mid_face_node(0);
  CPPUNIT_ASSERT( !set.have_any_mid_face_node() );
  CPPUNIT_ASSERT( !set.have_any_mid_node() );
  set.set_mid_face_node( NodeSet::NUM_FACE_BITS-1 );
  CPPUNIT_ASSERT( set.have_any_mid_face_node() );
  CPPUNIT_ASSERT( set.have_any_mid_node() );
  set.clear();

  CPPUNIT_ASSERT( !set.have_any_mid_region_node() );
  set.set_mid_region_node(0);
  CPPUNIT_ASSERT( set.have_any_mid_region_node() );
  CPPUNIT_ASSERT( set.have_any_mid_node() );
  set.clear_mid_region_node(0);
  CPPUNIT_ASSERT( !set.have_any_mid_region_node() );
  CPPUNIT_ASSERT( !set.have_any_mid_node() );
  set.set_mid_region_node( NodeSet::NUM_REGION_BITS-1 );
  CPPUNIT_ASSERT( set.have_any_mid_region_node() );
  CPPUNIT_ASSERT( set.have_any_mid_node() );
  set.clear();
}

void TriLagrangeShape::derivatives( Sample loc,
                                    NodeSet nodeset,
                                    size_t* vertex_indices_out,
                                    MsqVector<2>* d_coeff_d_xi_out,
                                    size_t& num_vtx,
                                    MsqError& err ) const
{
  if (!nodeset.have_any_mid_node()) {
    num_vtx = 3;
    get_linear_derivatives( vertex_indices_out, d_coeff_d_xi_out );
    return;
  }

  if (nodeset.have_any_mid_face_node()) {
    MSQ_SETERR(err)("TriLagrangeShape does not support mid-element nodes",
                    MsqError::UNSUPPORTED_ELEMENT);
    return;
  }
  
  switch (loc.dimension) {
    case 0:
      derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    case 1:
      derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    case 2:
      derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    default:
      MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
  }
}

void LinearPyramid::coefficients( Sample loc,
                                  NodeSet nodeset,
                                  double* coeff_out,
                                  size_t* indices_out,
                                  size_t& num_coeff,
                                  MsqError& err ) const
{
  if (nodeset.have_any_mid_node()) {
    MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
    return;
  }
  
  switch (loc.dimension) {
    case 0:
      coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
      break;
    case 1:
      coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
      break;
    case 2:
      coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
      break;
    case 3:
      coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
      break;
    default:
      MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
  }
}

void LinearHexahedron::derivatives( Sample loc,
                                    NodeSet nodeset,
                                    size_t* vertex_indices_out,
                                    MsqVector<3>* d_coeff_d_xi_out,
                                    size_t& num_vtx,
                                    MsqError& err ) const
{
  if (nodeset.have_any_mid_node()) {
    MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
    return;
  }
  
  switch (loc.dimension) {
    case 0:
      derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    case 1:
      derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    case 2:
      derivatives_at_mid_face( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    case 3:
      derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    default:
      MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
  }
}

void LinearTetrahedron::derivatives( Sample ,
                                     NodeSet nodeset,
                                     size_t* vertices,
                                     MsqVector<3>* coeff_derivs,
                                     size_t& num_vtx,
                                     MsqError& err ) const
{
  if (nodeset.have_any_mid_node()) {
    MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
    return;
  }
  else {
    num_vtx = 4;
    vertices[0] = 0;
    vertices[1] = 1;
    vertices[2] = 2;
    vertices[3] = 3;
    
    coeff_derivs[0][0] = -1.0;
    coeff_derivs[0][1] = -1.0;
    coeff_derivs[0][2] = -1.0;

    coeff_derivs[1][0] =  1.0;
    coeff_derivs[1][1] =  0.0;
    coeff_derivs[1][2] =  0.0;

    coeff_derivs[2][0] =  0.0;
    coeff_derivs[2][1] =  1.0;
    coeff_derivs[2][2] =  0.0;

    coeff_derivs[3][0] =  0.0;
    coeff_derivs[3][1] =  0.0;
    coeff_derivs[3][2] =  1.0;
  }
}

NodeSet TriLagrangeShape::sample_points( NodeSet ns ) const
{
  if (ns.have_any_mid_node()) {
    ns.set_all_corner_nodes(TRIANGLE);
  }
  else {
    ns.clear();
    ns.set_mid_face_node(0);
  }
  return ns;
}

void LinearTetrahedron::coefficients( Sample location,
                                      NodeSet nodeset,
                                      double* coeff_out,
                                      size_t* indices_out,
                                      size_t& num_coeff,
                                      MsqError& err ) const
{
  if (nodeset.have_any_mid_node()) {
    MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
    return;
  }
  
  switch (location.dimension) {
    case 0:
      num_coeff = 1;
      indices_out[0] = location.number;
      coeff_out[0] = 1.0;
      break;
    case 1:
      num_coeff = 2;
      coeff_out[0] = 0.5;
      coeff_out[1] = 0.5;
      if (location.number < 3) {
        indices_out[0] = location.number;
        indices_out[1] = (location.number+1)%3;
      }
      else {
        indices_out[0] = location.number - 3;
        indices_out[1] = 3;
      }
      break;
    case 2:
      num_coeff = 3;
      indices_out[0] = faces[location.number][0];
      indices_out[1] = faces[location.number][1];
      indices_out[2] = faces[location.number][2];
      coeff_out[0] = coeff_out[1] = coeff_out[2] = coeff_out[3] = MSQ_ONE_THIRD;
      break;
    case 3:
      num_coeff = 4;
      indices_out[0] = 0;
      indices_out[1] = 1;
      indices_out[2] = 2;
      indices_out[3] = 3;
      coeff_out[0] = coeff_out[1] = coeff_out[2] = coeff_out[3] = 0.25;
      break;
    default:
      MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
  }
}

void LinearPyramid::derivatives( Sample loc,
                                 NodeSet nodeset,
                                 size_t* vertex_indices_out,
                                 MsqVector<3>* d_coeff_d_xi_out,
                                 size_t& num_vtx,
                                 MsqError& err ) const
{
  if (nodeset.have_any_mid_node()) {
    MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
    return;
  }
  
  switch (loc.dimension) {
    case 0:
      if (loc.number == 4) {
        set_equal_derivatives( 0.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      }
      else {
        set_corner_derivatives( loc.number, 1.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      }
      break;
    case 1:
      if (loc.number < 4) {
        set_edge_derivatives( loc.number, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      }
      else {
        set_corner_derivatives( loc.number-4, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      }    
      break;
    case 2:
      if (loc.number == 4) {
        set_equal_derivatives( 0.5, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      }
      else {
        set_edge_derivatives( loc.number, 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      }
      break;
    case 3:
      set_equal_derivatives( 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
      break;
    default:
      MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
  }
}

static void derivatives_at_mid_elem( NodeSet nodeset,
                                     size_t* vertices,
                                     MsqVector<2>* derivs,
                                     size_t& num_vtx )
{
    // fast linear case
    // This is provided as an optimization for linear elements.
    // If this block of code were removed, the general-case code
    // below should produce the same result.
  if (!nodeset.have_any_mid_node()) {
    num_vtx = 4;
    vertices[0] = 0; derivs[0][0] = -0.5; derivs[0][1] = -0.5;
    vertices[1] = 1; derivs[1][0] =  0.5; derivs[1][1] = -0.5;
    vertices[2] = 2; derivs[2][0] =  0.5; derivs[2][1] =  0.5;
    vertices[3] = 3; derivs[3][0] = -0.5; derivs[3][1] =  0.5;
    return;
  }
  
  num_vtx = 0;
  
    // N_0
  if (!nodeset.both_edge_nodes(0,3)) {  // if eiter adjacent mid-edge node is missing
    vertices[num_vtx] = 0;
    derivs[num_vtx][0] = nodeset.mid_edge_node(3) ? 0.0 : -0.5;
    derivs[num_vtx][1] = nodeset.mid_edge_node(0) ? 0.0 : -0.5;
    ++num_vtx;
  }
  
    // N_1
  if (!nodeset.both_edge_nodes(0,1)) {  // if eiter adjacent mid-edge node is missing
    vertices[num_vtx] = 1;
    derivs[num_vtx][0] = nodeset.mid_edge_node(1) ? 0.0 :  0.5;
    derivs[num_vtx][1] = nodeset.mid_edge_node(0) ? 0.0 : -0.5;
    ++num_vtx;
  }
  
    // N_2
  if (!nodeset.both_edge_nodes(1,2)) {  // if eiter adjacent mid-edge node is missing
    vertices[num_vtx] = 2;
    derivs[num_vtx][0] = nodeset.mid_edge_node(1) ? 0.0 :  0.5;
    derivs[num_vtx][1] = nodeset.mid_edge_node(2) ? 0.0 :  0.5;
    ++num_vtx;
  }
  
    // N_3
  if (!nodeset.both_edge_nodes(2,3)) {  // if eiter adjacent mid-edge node is missing
    vertices[num_vtx] = 3;
    derivs[num_vtx][0] = nodeset.mid_edge_node(3) ? 0.0 : -0.5;
    derivs[num_vtx][1] = nodeset.mid_edge_node(2) ? 0.0 :  0.5;
    ++num_vtx;
  }
  
    // N_4
  if (nodeset.mid_edge_node(0)) {
    vertices[num_vtx] = 4;
    derivs[num_vtx][0] =  0.0;
    derivs[num_vtx][1] = -1.0;
    ++num_vtx;
  }
  
    // N_5
  if (nodeset.mid_edge_node(1)) {
    vertices[num_vtx] = 5;
    derivs[num_vtx][0] =  1.0;
    derivs[num_vtx][1] =  0.0;
    ++num_vtx;
  }
  
    // N_6
  if (nodeset.mid_edge_node(2)) {
    vertices[num_vtx] = 6;
    derivs[num_vtx][0] =  0.0;
    derivs[num_vtx][1] =  1.0;
    ++num_vtx;
  }
  
    // N_7
  if (nodeset.mid_edge_node(3)) {
    vertices[num_vtx] = 7;
    derivs[num_vtx][0] = -1.0;
    derivs[num_vtx][1] =  0.0;
    ++num_vtx;
  }
  
    // N_8 (mid-quad node) never contributes to Jacobian at element center!!!
}

static void derivatives_at_mid_edge( unsigned edge, 
                                     NodeSet nodeset,
                                     size_t* vertices,
                                     MsqVector<2>* derivs,
                                     size_t& num_vtx )
{
  static const double values[][9] = { {-0.5, -0.5, 0.5,  0.5, -1.0,  2.0,  1.0,  2.0,  4.0},
                                      {-0.5,  0.5, 0.5, -0.5, -2.0,  1.0, -2.0, -1.0, -4.0},
                                      {-0.5, -0.5, 0.5,  0.5, -1.0, -2.0,  1.0, -2.0, -4.0},
                                      {-0.5,  0.5, 0.5, -0.5,  2.0,  1.0,  2.0, -1.0,  4.0} };
  static const double edge_values[][2] = { {-1,  1},
                                           {-1,  1},
                                           { 1, -1},
                                           { 1, -1} }; 
  const unsigned prev_corner = edge;           // index of start vertex of edge
  const unsigned next_corner = (edge+1)%4;     // index of end vertex of edge
  const unsigned is_eta_edge = edge % 2;       // edge is xi = +/- 0
  const unsigned is_xi_edge  = 1 - is_eta_edge;// edge is eta = +/- 0
  //const unsigned mid_edge_node = edge + 4;     // mid-edge node index
  const unsigned prev_opposite = (prev_corner+3)%4; // index of corner adjacent to prev_corner
  const unsigned next_opposite = (next_corner+1)%4; // index of corner adjacent to next_corner;
 
    // First do derivatives along edge (e.g. wrt xi if edge is eta = +/-1)
  num_vtx = 2;
  vertices[0] = prev_corner;
  vertices[1] = next_corner;
  derivs[0][is_eta_edge] = edge_values[edge][0];
  derivs[0][is_xi_edge]  = 0.0;
  derivs[1][is_eta_edge] = edge_values[edge][1];
  derivs[1][is_xi_edge]  = 0.0;
    // That's it for the edge-direction derivatives.  No other vertices contribute.
    
    // Next handle the linear element case.  Handle this as a special case first,
    // so the generalized solution doesn't impact performance for linear elements
    // too much.
  if (!nodeset.have_any_mid_node()) {
    num_vtx = 4;
    vertices[2] = prev_opposite;
    vertices[3] = next_opposite;
    derivs[0][is_xi_edge] = values[edge][prev_corner];
    derivs[1][is_xi_edge] = values[edge][next_corner];
    derivs[2][is_xi_edge] = values[edge][prev_opposite];
    derivs[2][is_eta_edge] = 0.0;
    derivs[3][is_xi_edge] = values[edge][next_opposite];
    derivs[3][is_eta_edge] = 0.0;
    return;
  }
  
    // Initial (linear) contribution for each corner
  double v[4] = { values[edge][0], 
                  values[edge][1], 
                  values[edge][2], 
                  values[edge][3] };

    // If mid-face node is present
  double v8 = 0.0;
  if (nodeset.mid_face_node(0)) {
    v8 = values[edge][8];
    vertices[num_vtx] = 8;
    derivs[num_vtx][is_eta_edge] = 0.0;
    derivs[num_vtx][is_xi_edge] = v8;
    v[0] -= 0.25 * v8;
    v[1] -= 0.25 * v8;
    v[2] -= 0.25 * v8;
    v[3] -= 0.25 * v8;
    ++num_vtx;
  }

    // If mid-edge nodes are present
  for (unsigned i = 0; i < 4; ++i) {
    if (nodeset.mid_edge_node(i)) {
      const double value = values[edge][i+4] - 0.5 * v8;
      if (fabs(value) > 0.125) {
        v[ i     ] -= 0.5 * value;
        v[(i+1)%4] -= 0.5 * value;
        vertices[num_vtx] = i+4;
        derivs[num_vtx][is_eta_edge] = 0.0;
        derivs[num_vtx][is_xi_edge] = value;
        ++num_vtx;
      }
    }
  }

    // update values for adjacent corners
  derivs[0][is_xi_edge] = v[prev_corner];
  derivs[1][is_xi_edge] = v[next_corner];
    // do other two corners
  if (fabs(v[prev_opposite]) > 0.125) {
    vertices[num_vtx] = prev_opposite;
    derivs[num_vtx][is_eta_edge] = 0.0;
    derivs[num_vtx][is_xi_edge] = v[prev_opposite];
    ++num_vtx;
  }
  if (fabs(v[next_opposite]) > 0.125) {
    vertices[num_vtx] = next_opposite;
    derivs[num_vtx][is_eta_edge] = 0.0;
    derivs[num_vtx][is_xi_edge] = v[next_opposite];
    ++num_vtx;
  }
}