C++ PatchData::element_by_index方法代码示例

void CachingTargetCalculator::notify_sub_patch( PatchData& ,
                                                CachedTargetData& data,
                                                PatchData& subpatch, 
                                                const size_t* ,
                                                const size_t* element_map,
                                                MsqError& err )
{
    // If no cached data for this patch, just return
  if (data.has_data())
    return;
  
    // Create a new cached data object on the subpatch
  CachedTargetData& sub_data = get_data( subpatch );
  sub_data.clear();
  
    // populate the element offset list, and count the total
    // number of cached target matrices.
  sub_data.elementOffsets.resize( subpatch.num_elements() );
  size_t count_2D = 0, count_3D = 0;
  for (size_t i = 0; i < subpatch.num_elements(); ++i) {
    EntityTopology type = subpatch.element_by_index(i).get_element_type();
    size_t& count = (TopologyInfo::dimension( type ) == 2) ? count_2D : count_3D;
    sub_data.elementOffsets[i] = count;
    NodeSet samples = subpatch.get_samples( i );
    count += samples.num_nodes();
  }
  
  const bool orient = have_surface_orient();
  sub_data.targets3D.resize( count_3D );
  if (orient) 
    sub_data.targetsSurface.resize( count_2D );
  else
    sub_data.targets2D.resize( count_2D );

  for (size_t i = 0; i < subpatch.num_elements(); ++i) {
    EntityTopology type = subpatch.element_by_index(i).get_element_type();
    size_t off = sub_data.elementOffsets[i];
    size_t old_off = data.elementOffsets[element_map[i]];
    NodeSet samples = subpatch.get_samples( i );
    size_t count = samples.num_nodes();
   
    if (TopologyInfo::dimension( type ) == 3) 
      memcpy( &sub_data.targets3D[off], &data.targets3D[old_off], count*sizeof(MsqMatrix<3,3>) );
    else if (orient)
      memcpy( &sub_data.targetsSurface[off], &data.targetsSurface[old_off], count*sizeof(MsqMatrix<3,2>) );
    else
      memcpy( &sub_data.targets2D[off], &data.targets2D[old_off], count*sizeof(MsqMatrix<2,2>) );
  }
}

bool FauxMetric::evaluate_with_indices( PatchData& pd, size_t h, double& v, 
                              std::vector<size_t>& indices, MsqError& err )
{
  evaluate( pd, h, v, err );
  indices.resize(3);
  size_t e = ElemSampleQM::elem( h );
  Sample s = ElemSampleQM::sample( h );
  size_t* verts = pd.element_by_index(e).get_vertex_index_array();
  size_t n = pd.element_by_index(e).vertex_count();
  CPPUNIT_ASSERT_EQUAL( (unsigned short)0, s.dimension );
  indices[0] = verts[s.number];
  indices[1] = verts[(s.number+1)%n];
  indices[2] = verts[(s.number+n-1)%n];
  return true;
}

int NonGradient::getPatchDimension(const PatchData &pd, MsqError &err)
{
  const size_t numElt = pd.num_elements();
  unsigned edimMax = 0;  // in case numElt == 0
  for( size_t elementId = 0; elementId < numElt; elementId++)
  {
    const MsqMeshEntity& element = pd.element_by_index( elementId );
    EntityTopology type = element.get_element_type();
    unsigned edim = TopologyInfo::dimension( type );
    if( elementId == 0 ) 
    {
      edimMax = edim;
    } 
    else
    { 
      if(edimMax != edim)
      {
        MSQ_SETERR(err)("A patch has elements of different dimensions", MsqError::INVALID_MESH);
        std::cout << "A patch has elements of different dimensions" << edimMax << "  "  << edim << std::endl;
        if(edimMax < edim)
        {
          edimMax = edim;
        }
      }
    }
  }
  return( edimMax );
}

void MsqMeshEntityTest::test_all_nodes( EntityTopology type, unsigned num_nodes )
{
  const unsigned num_vtx = 27;
  double coords[3*num_vtx] = {0.0};
  size_t conn[num_vtx];
  for (size_t i = 0; i < num_vtx; ++i)
    conn[i] = i;
  bool fixed[num_vtx] = {false};
  CPPUNIT_ASSERT(num_nodes <= num_vtx);
  
  MsqError err;
  PatchData pd;
  size_t n = num_nodes;
  pd.fill( num_nodes, coords, 1, &type, &n, conn, fixed, err );
  ASSERT_NO_ERROR(err);
  
  MsqMeshEntity& elem = pd.element_by_index(0);
  NodeSet all = elem.all_nodes(err);
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( num_nodes, all.num_nodes() );
  CPPUNIT_ASSERT( all.have_any_corner_node() );
  bool mid_edge, mid_face, mid_reg;
  TopologyInfo::higher_order( type, num_nodes, mid_edge, mid_face, mid_reg, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( mid_edge, !!all.have_any_mid_edge_node() );
  CPPUNIT_ASSERT_EQUAL( mid_face, !!all.have_any_mid_face_node() );
  CPPUNIT_ASSERT_EQUAL( mid_reg,  !!all.have_any_mid_region_node() );
  
}

void MappingFunction3D::jacobian( const PatchData& pd,
                                  size_t element_number,
                                  NodeSet nodeset,
                                  Sample location,
                                  size_t* vertex_patch_indices_out,
                                  MsqVector<3>* d_coeff_d_xi_out,
                                  size_t& num_vtx_out,
                                  MsqMatrix<3,3>& jacobian_out,
                                  MsqError& err ) const
{
  const MsqMeshEntity& elem = pd.element_by_index( element_number );
  const size_t* conn = elem.get_vertex_index_array();
  
  derivatives( location, nodeset, vertex_patch_indices_out,
               d_coeff_d_xi_out, num_vtx_out, err ); MSQ_ERRRTN(err);
 
  convert_connectivity_indices( elem.node_count(), vertex_patch_indices_out, 
                                num_vtx_out, err );  MSQ_ERRRTN(err);
 
  jacobian_out.zero();
  size_t w = 0;
  for (size_t r = 0; r < num_vtx_out; ++r) {
    size_t i = conn[vertex_patch_indices_out[r]];
    MsqMatrix<3,1> coords( pd.vertex_by_index( i ).to_array() );
    jacobian_out += coords * transpose(d_coeff_d_xi_out[r]);
    if (i < pd.num_free_vertices()) {
      vertex_patch_indices_out[w] = i;
      d_coeff_d_xi_out[w] = d_coeff_d_xi_out[r];
      ++w;
    }
  }
  num_vtx_out = w;
}

bool AffineMapMetric::evaluate_with_indices( PatchData& pd,
                                             size_t handle,
                                             double& value,
                                             std::vector<size_t>& indices,
                                             MsqError& err )
{
  Sample   s = ElemSampleQM::sample( handle );
  size_t   e = ElemSampleQM::  elem( handle );
  MsqMeshEntity& elem = pd.element_by_index( e );
  EntityTopology type = elem.get_element_type();
  const size_t* conn = elem.get_vertex_index_array();
  
    // this metric only supports sampling at corners
  if (s.dimension != 0) {
    if (type != TRIANGLE && type != TETRAHEDRON) {
      MSQ_SETERR(err)("Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
      return false;
    }
    s.dimension = 0;
    s.number = 0;
  }

  unsigned n;
  const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
  indices.clear();
  if (conn[s.number] < pd.num_free_vertices())
    indices.push_back(conn[s.number]);
  for (unsigned i = 0; i < n; ++i)
    if (conn[adj[i]] < pd.num_free_vertices())
      indices.push_back(conn[adj[i]]);
  
  return evaluate( pd, handle, value, err );
}

bool SizeMetric::evaluate( PatchData& pd, 
                           size_t handle, 
                           double& value, 
                           MsqError& err )
{
  value = pd.element_by_index(handle).compute_unsigned_area( pd, err );
  return true;
}

bool FauxMetric::evaluate( PatchData& pd, size_t h, double& v, MsqError&  )
{
  size_t e = ElemSampleQM::elem( h );
  Sample s = ElemSampleQM::sample( h );
  size_t* verts = pd.element_by_index(e).get_vertex_index_array();
  CPPUNIT_ASSERT_EQUAL( (unsigned short)0, s.dimension );
  v = (double)(verts[s.number]);
  return true;
}

//UNTANGLED
void TerminationCriterionTest::test_untangled_mesh()
{
  MsqPrintError err(std::cout);
  PatchData pd;
  create_twelve_hex_patch( pd, err );
  ASSERT_NO_ERROR(err);
  
    // get two opposite vertices in first hexahedral element
  int vtx1 = pd.element_by_index(0).get_vertex_index_array()[0];
  int vtx2 = pd.element_by_index(0).get_vertex_index_array()[7];
  Vector3D saved_coords = pd.vertex_by_index(vtx2);
  Vector3D opposite_coords = pd.vertex_by_index(vtx1);

    // invert the element
  pd.move_vertex( 2*(opposite_coords-saved_coords), vtx2, err );
  ASSERT_NO_ERROR(err);
  int inverted, samples;
  pd.element_by_index(0).check_element_orientation( pd, inverted, samples, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(inverted > 0);
  
    // check initial termination criterion
  TerminationCriterion tc;
  tc.add_untangled_mesh();
  tc.reset_inner( pd, ofEval, err );
  ASSERT_NO_ERROR(err);
  tc.reset_patch( pd, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(!tc.terminate());
  
    // fix the element
  pd.set_vertex_coordinates( saved_coords, vtx2, err );
  ASSERT_NO_ERROR(err);
  pd.element_by_index(0).check_element_orientation( pd, inverted, samples, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL(0,inverted);
  
    // check that TC recognized untangled mesh
  tc.accumulate_inner( pd, 0.0, 0, err );
  ASSERT_NO_ERROR(err);
  tc.accumulate_patch( pd, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(tc.terminate());
}

bool TQualityMetric::evaluate_internal( PatchData& pd,
                                        size_t handle,
                                        double& value,
                                        size_t* indices,
                                        size_t& num_indices,
                                        MsqError& err )
{
  const Sample s = ElemSampleQM::sample( handle );
  const size_t e = ElemSampleQM::  elem( handle );
  MsqMeshEntity& elem = pd.element_by_index( e );
  EntityTopology type = elem.get_element_type();
  unsigned edim = TopologyInfo::dimension( type );
  const NodeSet bits = pd.non_slave_node_set( e );
  
  bool rval;
  if (edim == 3) { // 3x3 or 3x2 targets ?
    const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
    if (!mf) {
      MSQ_SETERR(err)( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
      return false;
    }

    MsqMatrix<3,3> A, W;
    mf->jacobian( pd, e, bits, s, indices, mDerivs3D, num_indices, A, err );
    MSQ_ERRZERO(err);
    targetCalc->get_3D_target( pd, e, s, W, err ); MSQ_ERRZERO(err);
    const MsqMatrix<3,3> Winv = inverse(W);
    const MsqMatrix<3,3> T = A*Winv;
    rval = targetMetric->evaluate( T, value, err ); MSQ_ERRZERO(err);
#ifdef PRINT_INFO
    print_info<3>( e, s, A, W, A * inverse(W) );
#endif
  }
  else if (edim == 2) {
    MsqMatrix<2,2> W, A;
    MsqMatrix<3,2> S_a_transpose_Theta;
    rval = evaluate_surface_common( pd, s, e, bits, indices, num_indices,
                                 mDerivs2D, W, A, S_a_transpose_Theta, err ); 
    if (MSQ_CHKERR(err) || !rval)
      return false;
    const MsqMatrix<2,2> Winv = inverse(W);
    const MsqMatrix<2,2> T = A * Winv;
    rval = targetMetric->evaluate( T, value, err );
    MSQ_ERRZERO(err);

#ifdef PRINT_INFO
    print_info<2>( e, s, J, Wp, A * inverse(W) );
#endif
  }
  else {
    assert(false);
    return false;
  }
  
  return rval;
}

size_t TerminationCriterion::count_inverted( PatchData& pd, MsqError& err )
{
  size_t num_elem = pd.num_elements();
  size_t count=0;
  int inverted, samples;
  for (size_t i = 0; i < num_elem; i++) {
    pd.element_by_index(i).check_element_orientation(pd, inverted, samples, err);
    if (inverted)
      ++count;
  }
  return count;
}

size_t SmartLaplacianSmoother::num_inverted( PatchData& pd, MsqError& err )
{
  size_t result = 0;
  int inverted, junk;
  for (size_t i = 0; i < pd.num_elements(); ++i) {
    pd.element_by_index(i).check_element_orientation( pd, inverted, junk, err ); 
    MSQ_ERRZERO(err);
    if (inverted)
      ++result;
  }
  return result;
}

void MsqMeshEntityTest::test_unsigned_area_common( EntityTopology type,
                                                   const double* coords,
                                                   double expected )
{
  MsqError err;
  PatchData pd;

  pd.fill( TopologyInfo::corners( type ), coords, 1, type, conn, fixed, err );
  ASSERT_NO_ERROR(err);
  
  double a = pd.element_by_index(0).compute_unsigned_area( pd, err );
  ASSERT_NO_ERROR(err);
  
  CPPUNIT_ASSERT_DOUBLES_EQUAL( expected, a, 1e-8 );
}

void MsqHessianTest::accumulate_entries( const PatchData &pd, 
                                         size_t elem_index,
                                         const Matrix3D* mat3d_array, 
                                         MsqError &err )
{
  const MsqMeshEntity& elem = pd.element_by_index( elem_index );
  const size_t nv = elem.vertex_count();
  const size_t* v = elem.get_vertex_index_array();
  for (size_t r = 0; r < nv; ++r) {
    for (size_t c = r; c < nv; ++c) {
      add( v[r], v[c], *mat3d_array, err );
      CPPUNIT_ASSERT(!err);
      ++mat3d_array;
    }
  }
}

bool EdgeQM::evaluate_with_indices( PatchData& pd,
                                    size_t handle,
                                    double& value,
                                    std::vector<size_t>& indices,
                                    MsqError& err )
{
  const MsqMeshEntity& element = pd.element_by_index( elem(handle) );
  EntityTopology type = element.get_element_type(); 
  const unsigned* verts = TopologyInfo::edge_vertices( type, edge(handle) );
  const size_t* conn = element.get_vertex_index_array();
  indices.clear();
  if (conn[verts[0]] < pd.num_free_vertices())
    indices.push_back( conn[verts[0]] );
  if (conn[verts[1]] < pd.num_free_vertices())
    indices.push_back( conn[verts[1]] );
  return evaluate( pd, handle, value, err );
}