double.cpp

来自「算断裂的」· C++ 代码 · 共 1,820 行 · 第 1/5 页

  bool FspecCopyNum::operator<(const FspecCopyNum& o) const {
    if (fspec < o.fspec) return true;
    if (fspec > o.fspec) return false;
    return copynum < o.copynum;
  }




  // Routine to determine which side of a surface has a certain
  // simplex.

  int determine_side(const Brep& b,
                     const Brep::Face_Spec& owner,
                     const SimpComplex& sc,
                     Brep::Face_Spec regionspec,
                     int seqno,
                     int which_facet,
                     const map<PatchInBrepAddress,bool>& patch_is_flipped,
                     double scfac,
                     double tol) {
    int di = b.embedded_dim();
    Point startpoint, basepoint, endpoint, tang, coef;
    PatchMath::Workspace workspace;
    vector<PatchMath::IntersectRecord> outputstack;

    {
      SimpComplex::VertexGlobalIndex vnum = 
        sc.node_of_meshface_on_face(regionspec, seqno, which_facet);
      SimpComplex::VertexOrdinalIndex vnumo = sc.global_to_ordinal(vnum);
      for (int j = 0; j < di; ++j)
        startpoint[j] = sc.real_coord_o(vnumo, j);
    }
    int side;
    bool degen = true;
    int trycount = 0;
    while (degen &&  ++trycount < 400) {
      degen = false;

      // Select a point in the base facet of the simplex
      // via interpolation with randomized coeffs.

      Real t = 1.0;
      {
        for (int j = 0; j < di - 1; ++j) {
          coef[j] = 1.0/((double)(2*di)) + 1.0/((double)(2*di)) * random_real();
          t -= coef[j];
        }
      }
      coef[di - 1] = t;
      {
        for (int i = 0; i < di; ++i)
          basepoint[i] = 0.0;
      }
      {
        int k = 0;
        for (int j = 0; j < di + 1; ++j) {
          if (j == which_facet) continue;
          SimpComplex::VertexGlobalIndex vnum = 
            sc.node_of_meshface_on_face(regionspec, seqno, j);
          SimpComplex::VertexOrdinalIndex vnumo = sc.global_to_ordinal(vnum);
          for (int i = 0; i < di; ++i)
            basepoint[i] += coef[k] * sc.real_coord_o(vnumo,i);
          ++k;
        }
      }
      // endpoint extrapolates from startpoint to basepoint by
      // a factor of 5.

      {
        for (int j = 0; j < di; ++j) {
          tang[j] = startpoint[j] - basepoint[j];
          endpoint[j] = startpoint[j] - 5.0 * tang[j];
        }
      }

      // Loop over patches of the face; for each patch find
      // the point where the segment from startpoint to
      // endpoint crosses.  Save the first hit.

      Real first_hit_dist = BIG_REAL;

      for (Brep::Loop_over_patches_of_face ptloop(b, owner);
           ptloop.notdone();
           ++ptloop) {
        outputstack.resize(0, PatchMath::IntersectRecord());
        ptloop.patchmath().intersect_seg(outputstack, workspace,
                                         startpoint, endpoint, scfac, tol, true);
        for (int j = 0; j < outputstack.size(); ++j) {
          if (outputstack[j].degen) {
            degen = true;
            break;
          }
          if (outputstack[j].dist < first_hit_dist) {
            first_hit_dist = outputstack[j].dist;
            Real ip = Point::inner_product(ptloop.patchmath()
                                           .normal(outputstack[j].paramcoord), tang);
            typedef map<PatchInBrepAddress,bool> MPIBB;
            MPIBB::const_iterator it = 
              patch_is_flipped.find(PatchInBrepAddress(owner, ptloop.index()));
            if (it == patch_is_flipped.end()) 
              throw_error("record missing from patch_is_flipped");
            if (it -> second)
              ip = -ip;
            if (ip < 0)
              side = 0;
            else
              side = 1;
          }
        }
        if (degen) break;
      }
      if (first_hit_dist > 10)
        degen = true;
    }
    if (degen)
      throw_error("Unrecoverable degeneracy in 'determine_side'");
    return side;
  }

  
}


// This routine doubles the internal boundary faces of a brep.
// The faces to double are given by the second argument.

QMG::MG::Brep_Under_Construction
QMG::MG::double_brep(const Brep& b,
                     const vector<Brep::Face_Spec>& faces_to_double,
                     double tol,
                     // other return variables:
                     Brep::Face_Labeling<Brep::FaceIndex>& oldfspec_to_new_map,
                     Brep::Face_Labeling<int>& face_num_copies,
                     map<PatchInBrepAddress, bool>& patch_is_flipped) {
  double scfac = b.scale_factor();
  AnglePatchIndPair::tol = tol;
  SegHitRecord::tol = tol;
  using std::ofstream;
  using std::endl;
  ostream logfile(0);
  Logstream logstr(logfile, 1);
  OnErrorFlushLog _unnamed_(logstr);
  
  // initialize the random number generator.
  
  random_real(1);
  int di = b.embedded_dim();
  small_initialize_static_data(di);
  Point::InitStaticData(di);
#ifdef DEBUGGING
  Local_to_PatchMath_CPP::dump_everything = false;
  Local_to_PatchMath_CPP::debug_str = &logstr.str();
#endif
  
  logstr.str() << "tol = " << tol << " scfac = " << scfac << '\n';
  PatchTable patchtable(b, logstr, scfac, tol);
  
  // Make a mapping from original brep's patchind index
  // to patchtable's patch index.

  map<PatchInBrepAddress, PatchTable::Index> lookup_patchind; 
  {
    for (PatchTable::Index patchind = 0; 
         patchind < patchtable.numpatch(); 
         ++patchind) {
      Brep::Face_Spec owner = patchtable.owner(patchind);
      if (patchtable.gdim(patchind) != owner.fdim()) continue;
      Brep::PatchIndex patchind1 = patchtable.orig_index(patchind);
      lookup_patchind[PatchInBrepAddress(owner,patchind1)] = patchind;
    }
  }
  
  // Mark faces that must be doubled.

  Brep::Face_Labeling<My_bool> must_be_doubled(b, My_bool(false));
  Brep::Face_Labeling<Brep::Face_Spec> 
    parent_of_doubled_face(b, Brep::Face_Spec());
  Brep::Ancestor_Lookup anc_lookup(b);

  // used only in 3D

  vector<bool> is_ib_ancestor_of_vtx(b.level_size(di - 1));


  {
    int ftds = faces_to_double.size();
    for (int j = 0; j < ftds; ++j) {
      Brep::Face_Spec fspec = faces_to_double[j];
      if (fspec.fdim() != di - 1)
        throw_error("Only faces of dimension di-1 may be doubled");
      must_be_doubled(fspec) = true;
      Brep::Face_Spec parent_fspec(-1,-1);
      for (Brep::Ancestor_Lookup::Loop_over_ancestors ancloop(anc_lookup, fspec);
           ancloop.notdone();
           ++ancloop) {
        Brep::Face_Spec afspec = ancloop.ancestor_fspec();
        if (afspec.fdim() != fspec.fdim() + 1) continue;
        if (ancloop.occurrence_count() % 2 != 0)
        throw_error("Specified face is not an internal boundary of parent topological entity");
        if (parent_fspec.fdim() >= 0)
          throw_error("Specified face occurs as boundary of more than one parent topological entity");
        parent_fspec = afspec;
      }
      if (parent_fspec.fdim() < 0)
        throw_error("Specified face is not a boundary of any parent topological entity");
      parent_of_doubled_face(fspec) = parent_fspec;
    }
  }

  int patchtable_size = patchtable.numpatch();
  vector<bool> patch_is_doubled(patchtable_size);

  // Mark doubled patches.
  {
    for (PatchTable::Index patchind = 0; patchind < patchtable_size; ++patchind) {
      patch_is_doubled[patchind] =  (patchtable.gdim(patchind) == di - 1) &&
        must_be_doubled(patchtable.owner(patchind));
    }
  }

  vector<AnglePatchIndPair> bigsortlist;
  vector<int> pos_in_bigsortlist(patchtable_size);
  vector<int> len_in_bigsortlist(patchtable_size, 0);
          
  // Make a table that holds the patches around
  // each bezier curve (3D) or point (2D) in order of angle.

  {
    vector<AnglePatchIndPair> sortlist;
    for (PatchTable::Index patchind = 0; patchind < patchtable_size; ++patchind) {
      if (patchtable.gdim(patchind) != di - 2) continue;

      // For a patch of dimension di-2, find all its parents
      // and properly list them in the big sort list  data structure.
      
      
      Point hhtransform;
      Point paramcoord;
      Real beta;
      
      // Get the angle that the parent patch makes with respect
      // to this patch.

      if (di == 3) {
        Point direc;
        direc[0] = 1.0;
        paramcoord[0] = 0.5;
        Point tangent = 
          patchtable.patchmath(patchind).direc_deriv(paramcoord, direc);
        beta = Matrix::compute_hh_transform(&hhtransform[0], &tangent[0], 1, 3);
      }

      sortlist.resize(0, AnglePatchIndPair());
      for (PatchTable::Loop_over_parents parloop(patchtable, patchind);
           parloop.notdone();
           ++parloop) {
        PatchTable::Index ppatchind = parloop.parent_index();
        Point parent_tan = 
          patchtable.get_outward_tan_in_parent(paramcoord,
                                               patchind, 
                                               parloop);
        Point liftcoord = patchtable.lift_parameters(paramcoord,
                                                     patchind,
                                                     parloop.listind());
        Point parent_normal = patchtable.patchmath(ppatchind).
          normal(liftcoord);
        if (di == 3) {
          PatchTable::apply_3d_hh_and_shift(beta, hhtransform, parent_tan);
          PatchTable::apply_3d_hh_and_shift(beta, hhtransform, parent_normal);
        }
        Real tanang = atan2(parent_tan[1], parent_tan[0]);
        Real nrmang = atan2(parent_normal[1], parent_normal[0]);
        if (nrmang < tanang)
          nrmang += 2 * PI;
        bool is_forward =  fabs(nrmang - tanang - PI_OVER_2) < PI_OVER_2;

        Brep::Face_Spec front = patchtable.front_side(ppatchind);
        Brep::Face_Spec back = patchtable.back_side(ppatchind);

        if (!patch_is_doubled[ppatchind]) {
          
          sortlist.push_back(AnglePatchIndPair(tanang,
                                               ppatchind, 
                                               0, 
                                               is_forward,
                                               is_forward? front: back,
                                               is_forward? back : front));
        }
        else {
          
          // if this is a doubled patch, put two copies
          // of the patch in the sort list.  Assign epsilon so that
          // the copy with the clockwise normal comes second.
          
          Brep::Face_Spec nullspec(di,-1);
          sortlist.push_back(AnglePatchIndPair(tanang, 
                                               ppatchind, 
                                               0,
                                               is_forward,
                                               is_forward? front : nullspec,
                                               is_forward? nullspec : front));

          sortlist.push_back(AnglePatchIndPair(tanang, 
                                               ppatchind, 
                                               1,
                                               is_forward,
                                               is_forward? nullspec : back,
                                               is_forward? back : nullspec));
        }
      }

      // Sort the parent patches in order of angle.
      
      sort(sortlist.begin(), sortlist.end());
      int sortlistlen = sortlist.size();

      // Find a starting point for the sweep.
      // A starting point has the brep on its clockwise side.

      int sweepstart = -1;
      int both_covered_count = 0;
      {
        Brep::Face_Spec prev_region_above(-1,-1), first_region_below(-1,-1);
        for (int ii = 0; ii < sortlistlen; ++ii) {
          PatchTable::Index ppatchind = sortlist[ii].patchind;
          int wc = sortlist[ii].which_copy;
          bool isfor = sortlist[ii].is_forward;
          Brep::Face_Spec region_above = sortlist[ii].region_above;
          Brep::Face_Spec region_below = sortlist[ii].region_below;
          if (ii > 0 && region_below != prev_region_above)
            throw_error("Front/back inconsistency");
          if (region_below.faceind() < 0 && region_above.faceind() >= 0)
            sweepstart = ii;
          if (region_below.faceind() >= 0 && region_above.faceind() >= 0)
            ++both_covered_count;
          prev_region_above = region_above;
          if (ii == 0)
            first_region_below = region_below;
        }
        if (first_region_below != prev_region_above)
          throw_error("Front back inconsistency 2");
      }
      if (sweepstart < 0) {
        if (both_covered_count == sortlistlen) {
          sweepstart = 0;
        }
        else {
          throw_error("ridge does not touch brep region");
        }
      }
      int base = bigsortlist.size();
      bigsortlist.resize(base + sortlistlen, AnglePatchIndPair());
      pos_in_bigsortlist[patchind] = base;
      len_in_bigsortlist[patchind] = sortlistlen;
      {
        for (int ii = 0; ii < sortlistlen; ++ii) {
          int srcii = (ii + sweepstart) % sortlistlen;
          bigsortlist[base + ii] = sortlist[srcii];
        }
      }
    }
  }