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

void TriLagrangeShape::coefficients( Sample loc,
                                     NodeSet nodeset,
                                     double* coeff_out,
                                     size_t* indices_out,
                                     size_t& num_coeff,
                                     MsqError& err ) const
{
  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:
      num_coeff = 1;
      indices_out[0] = loc.number;
      coeff_out[0] = 1.0;
      break;
    case 1:
      if (nodeset.mid_edge_node(loc.number)) { // if mid-edge node is present
        num_coeff = 1;
        indices_out[0] = 3+loc.number;
        coeff_out[0] = 1.0;
      }
      else { // no mid node on edge
        num_coeff = 2;
        indices_out[0] = loc.number;
        indices_out[1] = (loc.number+1)%3;
        coeff_out[0] = 0.5;
        coeff_out[1] = 0.5;
      }
      break;
    case 2:
      num_coeff = 3;
      indices_out[0] = 0;
      indices_out[1] = 1;
      indices_out[2] = 2;
      coeff_out[0] = 1.0/3.0;
      coeff_out[1] = 1.0/3.0;
      coeff_out[2] = 1.0/3.0;
      for (int i = 0; i < 3; ++i) { // for each mid-edge node
        if (nodeset.mid_edge_node(i)) { // if node is present
          indices_out[num_coeff] = i+3;
            // coeff for mid-edge node
          coeff_out[num_coeff] = 4.0/9.0;
            // adjust coeff for adj corner nodes
          coeff_out[i]       -= 2.0/9.0;
          coeff_out[(i+1)%3] -= 2.0/9.0;
            // update count
          ++num_coeff;
        }
      }
      break;
    default:
      MSQ_SETERR(err)(MsqError::UNSUPPORTED_ELEMENT,
                  "Request for dimension %d mapping function value"
                  "for a triangular element", loc.dimension);
  }
}

static void compare_coefficients( const double* coeffs,
                                  const size_t* indices,
                                  const double* expected_coeffs,
                                  size_t num_coeff,
                                  unsigned loc, NodeSet bits )
{
    // find the location in the returned list for each node
  size_t revidx[6];
  double test_vals[6];
  for (size_t i = 0; i < 6; ++i) {
    revidx[i] = std::find( indices, indices+num_coeff, i ) - indices;
    test_vals[i] = (revidx[i] == num_coeff) ? 0.0 : coeffs[revidx[i]];
  }

    // Check that index list doesn't contain any nodes not actually
    // present in the element.
  CPPUNIT_ASSERT( bits.mid_edge_node(0) || (revidx[3] == num_coeff) );
  CPPUNIT_ASSERT( bits.mid_edge_node(1) || (revidx[4] == num_coeff) );
  CPPUNIT_ASSERT( bits.mid_edge_node(2) || (revidx[5] == num_coeff) );
    
    // compare expected and actual coefficient values
  ASSERT_VALUES_EQUAL( expected_coeffs[0], test_vals[0], loc, bits );
  ASSERT_VALUES_EQUAL( expected_coeffs[1], test_vals[1], loc, bits );
  ASSERT_VALUES_EQUAL( expected_coeffs[2], test_vals[2], loc, bits );
  ASSERT_VALUES_EQUAL( expected_coeffs[3], test_vals[3], loc, bits );
  ASSERT_VALUES_EQUAL( expected_coeffs[4], test_vals[4], loc, bits );
  ASSERT_VALUES_EQUAL( expected_coeffs[5], test_vals[5], loc, bits );
}

void NodeSetTest::test_set_get_clear_simple()
{
  NodeSet set;
  for (unsigned i = 0; i < NodeSet::NUM_CORNER_BITS; ++i) {
    CPPUNIT_ASSERT( !set.corner_node(i) );
    set.set_corner_node( i );
    CPPUNIT_ASSERT( set.corner_node(i) );
    set.clear_corner_node( i );
    CPPUNIT_ASSERT( !set.corner_node(i) );
  }
  for (unsigned i = 0; i < NodeSet::NUM_EDGE_BITS; ++i) {
    CPPUNIT_ASSERT( !set.mid_edge_node(i) );
    set.set_mid_edge_node( i );
    CPPUNIT_ASSERT( set.mid_edge_node(i) );
    set.clear_mid_edge_node( i );
    CPPUNIT_ASSERT( !set.mid_edge_node(i) );
  }
  for (unsigned i = 0; i < NodeSet::NUM_FACE_BITS; ++i) {
    CPPUNIT_ASSERT( !set.mid_face_node(i) );
    set.set_mid_face_node( i );
    CPPUNIT_ASSERT( set.mid_face_node(i) );
    set.clear_mid_face_node( i );
    CPPUNIT_ASSERT( !set.mid_face_node(i) );
  }
  for (unsigned i = 0; i < NodeSet::NUM_REGION_BITS; ++i) {
    CPPUNIT_ASSERT( !set.mid_region_node(i) );
    set.set_mid_region_node( i );
    CPPUNIT_ASSERT( set.mid_region_node(i) );
    set.clear_mid_region_node( i );
    CPPUNIT_ASSERT( !set.mid_region_node(i) );
  }
}

static void derivatives_at_corner( unsigned corner, 
                                   NodeSet nodeset,
                                   size_t* vertices,
                                   MsqVector<2>* derivs,
                                   size_t& num_vtx )
{
  static const unsigned xi_adj_corners[]  = { 1, 0, 3, 2 };
  static const unsigned xi_adj_edges[]    = { 0, 0, 2, 2 };
  static const unsigned eta_adj_corners[] = { 3, 2, 1, 0 };
  static const unsigned eta_adj_edges[]   = { 3, 1, 1, 3 };
  
  static const double corner_xi[]  = { -1,  1,  1, -1 }; // xi values by corner
  static const double corner_eta[] = { -1, -1,  1,  1 }; // eta values by corner
  static const double other_xi[]   = {  1, -1, -1,  1 }; // xi values for adjacent corner in xi direction
  static const double other_eta[]  = {  1,  1, -1, -1 }; // eta values for adjcent corner in eta direction
  static const double mid_xi[]     = {  2, -2, -2,  2 }; // xi values for mid-node in xi direction
  static const double mid_eta[]    = {  2,  2, -2, -2 }; // eta values for mid-node in eta direction
  
  num_vtx = 3;
  vertices[0] = corner;
  vertices[1] = xi_adj_corners[corner];
  vertices[2] = eta_adj_corners[corner];

  derivs[0][0] = corner_xi [corner];
  derivs[0][1] = corner_eta[corner];
  derivs[1][0] = other_xi  [corner];
  derivs[1][1] = 0.0;
  derivs[2][0] = 0.0;
  derivs[2][1] = other_eta [corner];

  if (nodeset.mid_edge_node(xi_adj_edges[corner])) {
    vertices[num_vtx] = 4 + xi_adj_edges[corner];
    derivs[num_vtx][0] = 2.0*mid_xi[corner];
    derivs[num_vtx][1] = 0.0;
    derivs[0][0] -= mid_xi[corner];
    derivs[1][0] -= mid_xi[corner];
    ++num_vtx;
  }

  if (nodeset.mid_edge_node(eta_adj_edges[corner])) {
    vertices[num_vtx] = 4 + eta_adj_edges[corner];
    derivs[num_vtx][0] = 0.0;
    derivs[num_vtx][1] = 2.0*mid_eta[corner];
    derivs[0][1] -= mid_eta[corner];
    derivs[2][1] -= mid_eta[corner];
    ++num_vtx;
  }
}

static void derivatives_at_mid_edge( unsigned edge, 
                                     NodeSet nodeset,
                                     size_t* vertices,
                                     MsqVector<2>* derivs,
                                     size_t& num_vtx )
{
    // The mid-edge behavior is rather strange.
    // A corner vertex contributes to the jacobian
    // at the mid-point of the opposite edge unless
    // one, but *not* both of the adjacent mid-edge
    // nodes is present.
    
    // The value for each corner is incremented when 
    // a mid-side node is present.  If the value is
    // 0 when finished, the corner doesn't contribute.
    // Initialize values to 0 for corners adjacent to
    // edge so they are never zero.
  int corner_count[3] = { 1, 1, 1 };
  corner_count[(edge+2)%3] = -1;
  
    // begin with derivatives for linear tri
  double corner_derivs[3][2] = { {-1.0,-1.0 },
                                 { 1.0, 0.0 },
                                 { 0.0, 1.0 } };
  
    // do mid-side nodes first
  num_vtx = 0;
  for (unsigned i = 0; i < 3; ++i) {
    if (nodeset.mid_edge_node(i)) {
      vertices[num_vtx] =  i+3 ;
      derivs[num_vtx][0] =  edr[i][edge] ;
      derivs[num_vtx][1] =  eds[i][edge] ;
      ++num_vtx;

      int a = (i+1)%3;
      corner_derivs[i][0] -= 0.5 * edr[i][edge];
      corner_derivs[i][1] -= 0.5 * eds[i][edge];
      corner_derivs[a][0] -= 0.5 * edr[i][edge];
      corner_derivs[a][1] -= 0.5 * eds[i][edge];
      ++corner_count[i];
      ++corner_count[a];
    }
  }
  
    // now add corner nodes to list
  for (unsigned i = 0; i < 3; ++i) {
    if (corner_count[i]) {
      vertices[num_vtx] = i;
      derivs[num_vtx][0] = corner_derivs[i][0];
      derivs[num_vtx][1] = corner_derivs[i][1];
      ++num_vtx;
    }
  }
}

static void get_coeff( NodeSet nodeset, const double* rs, double* coeffs )
{
  for (int i = 0; i < 6; ++i) 
    coeffs[i] = (*N[i])(rs[0], rs[1]);
  for (int i = 0; i < 3; ++i) 
    if (!nodeset.mid_edge_node(i)) {
      coeffs[i]      += 0.5 * coeffs[i+3];
      coeffs[(i+1)%3] += 0.5 * coeffs[i+3];
      coeffs[i+3] = 0;
    }
}

static void get_derivs( NodeSet nodeset, const double* rs, double* derivs )
{
  for (int i = 0; i < 6; ++i) {
    derivs[2*i  ] = (*dNdr[i])(rs[0], rs[1]);
    derivs[2*i+1] = (*dNds[i])(rs[0], rs[1]);
  }
  for (int i = 0; i < 3; ++i) 
    if (!nodeset.mid_edge_node(i)) {
      derivs[2*i]        += 0.5 * derivs[2*i+6];
      derivs[2*i+1]      += 0.5 * derivs[2*i+7];
      int j = (i+1)%3;
      derivs[2*j]        += 0.5 * derivs[2*i+6];
      derivs[2*j+1]      += 0.5 * derivs[2*i+7];
      derivs[2*i+6] = 0;
      derivs[2*i+7] = 0;
    }
}

static void derivatives_at_mid_elem( NodeSet nodeset,
                                     size_t* vertices,
                                     MsqVector<2>* derivs,
                                     size_t& num_vtx )
{
  get_linear_derivatives( vertices, derivs );
  num_vtx = 3;
  for (unsigned i = 0; i < 3; ++i) {
    if (nodeset.mid_edge_node(i)) {
      vertices[num_vtx] = i+3;
      derivs[num_vtx][0] = fdr[i];
      derivs[num_vtx][1] = fds[i];
      ++num_vtx;
       
      int a = (i+1)%3;
      derivs[i][0] -= 0.5 * fdr[i];
      derivs[i][1] -= 0.5 * fds[i];
      derivs[a][0] -= 0.5 * fdr[i];
      derivs[a][1] -= 0.5 * fds[i];
    }
  }
}

void QuadLagrangeShape::coefficients( Sample loc,
                                      NodeSet nodeset,
                                      double* coeff_out,
                                      size_t* indices_out,
                                      size_t& num_coeff,
                                      MsqError& err ) const
{
  switch (loc.dimension) {
    case 0:
      num_coeff = 1;
      indices_out[0] = loc.number;
      coeff_out[0] = 1.0;
      break;
    case 1:
      coeff_out[0] = coeff_out[1] = coeff_out[2] =
      coeff_out[3] = coeff_out[4] = coeff_out[5] = 
      coeff_out[6] = coeff_out[7] = coeff_out[8] = 0.0;
      if (nodeset.mid_edge_node(loc.number)) {  
          // if mid-edge node is present
        num_coeff = 1;
        indices_out[0] = loc.number+4;
        coeff_out[0] = 1.0;
      }
      else {
          // If mid-edge node is not present, mapping function value
          // for linear edge is even weight of adjacent vertices.
        num_coeff = 2;
        indices_out[0] = loc.number;
        indices_out[1] = (loc.number+1)%4;
        coeff_out[0] = 0.5;
        coeff_out[1] = 0.5;
      }
      break;
    case 2:
      if (nodeset.mid_face_node(0)) { // if quad center node is present
        num_coeff = 1;
        indices_out[0] = 8;
        coeff_out[0] = 1.0;
      } 
      else {
          // for linear element, (no mid-edge nodes), all corners contribute 1/4.
        num_coeff = 4;
        indices_out[0] = 0;
        indices_out[1] = 1;
        indices_out[2] = 2;
        indices_out[3] = 3;
        coeff_out[0] = 0.25;
        coeff_out[1] = 0.25;
        coeff_out[2] = 0.25;
        coeff_out[3] = 0.25;
          // add in contribution for any mid-edge nodes present
        for (int i = 0; i < 4; ++i) { // for each edge
          if (nodeset.mid_edge_node(i))
          {
            indices_out[num_coeff] = i+4;
            coeff_out[num_coeff] = 0.5;
            coeff_out[ i     ] -= 0.25;
            coeff_out[(i+1)%4] -= 0.25;
            ++num_coeff;
          }
        }
      }
      break;
    default:
      MSQ_SETERR(err)(MsqError::UNSUPPORTED_ELEMENT,
                  "Request for dimension %d mapping function value"
                  "for a quadrilateral element", loc.dimension);
  }
}

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;
  }
}

static void derivatives_at_corner( unsigned corner,
                                   NodeSet nodeset,
                                   size_t* vertices,
                                   MsqVector<2>* derivs,
                                   size_t& num_vtx )
{
  num_vtx = 3;
  get_linear_derivatives( vertices, derivs );
  switch (corner) {
    case 0:
      if (nodeset.mid_edge_node(0)) {
        vertices[num_vtx] = 3;
        derivs[num_vtx][0] = 4.0;
        derivs[num_vtx][1] = 0.0;
        ++num_vtx;
        derivs[0][0] -= 2.0;
        derivs[1][0] -= 2.0;
      }
      if (nodeset.mid_edge_node(2)) {
        vertices[num_vtx] = 5;
        derivs[num_vtx][0] = 0.0;
        derivs[num_vtx][1] = 4.0;
        ++num_vtx;
        derivs[0][1] -= 2.0;
        derivs[2][1] -= 2.0;
      }
      break;
    
    case 1:
      if (nodeset.mid_edge_node(0)) {
        vertices[num_vtx] = 3;
        derivs[num_vtx][0] = -4.0;
        derivs[num_vtx][1] = -4.0;
        ++num_vtx;
        derivs[0][0] += 2.0;
        derivs[0][1] += 2.0;
        derivs[1][0] += 2.0;
        derivs[1][1] += 2.0;
      }
      if (nodeset.mid_edge_node(1)) {
        vertices[num_vtx] = 4;
        derivs[num_vtx][0] = 0.0;
        derivs[num_vtx][1] = 4.0;
        ++num_vtx;
        derivs[1][1] -= 2.0;
        derivs[2][1] -= 2.0;
      }
      break;
    
    case 2:
      if (nodeset.mid_edge_node(1)) {
        vertices[num_vtx] = 4;
        derivs[num_vtx][0] = 4.0;
        derivs[num_vtx][1] = 0.0;
        ++num_vtx;
        derivs[1][0] -= 2.0;
        derivs[2][0] -= 2.0;
      }
      if (nodeset.mid_edge_node(2)) {
        vertices[num_vtx] = 5;
        derivs[num_vtx][0] = -4.0;
        derivs[num_vtx][1] = -4.0;
        ++num_vtx;
        derivs[0][0] += 2.0;
        derivs[0][1] += 2.0;
        derivs[2][0] += 2.0;
        derivs[2][1] += 2.0;
      }
      break;
  }
}