gmcoarsetopo.cpp

来自「算断裂的」· C++ 代码 · 共 924 行 · 第 1/3 页

              string otherval = b.face_lookup_prop_val(parent2, pvloop.prop());
              if (pvloop.val() != otherval) {
                // ctlog << " not deleting; parents differ in property " << pvloop.prop()
                //  << " values = " << pvloop.val() << "," << otherval << '\n';
                ok_to_delete(fspec) = false;
                break;
              }
            }
            
            if (!ok_to_delete(fspec)) continue;
            for (Brep::Loop_over_propvals_of_face pvloop2(b, parent2);
            pvloop2.notdone();
            ++pvloop2) {
              string otherval = b.face_lookup_prop_val(parent1, pvloop2.prop());
              if (pvloop2.val() != otherval) {
                //ctlog << " not deleting; parents differ in property " << pvloop2.prop()
                //  << " values = " << pvloop2.val() << "," << otherval << '\n';
                ok_to_delete(fspec) = false;
                break;
              }
            }
            if (!ok_to_delete(fspec)) continue;

            // Check that the two parents have the same parents.

            parent1_parentlabels.resize(0);

            for (Brep::Ancestor_Lookup::Loop_over_ancestors ancloop1(anc_lookup, parent1);
            ancloop1.notdone();
            ++ancloop1) {
              Brep::Face_Spec afspec = ancloop1.ancestor_fspec();
              if (afspec.fdim() != fdim + 2) continue;
              if (ok_to_delete(afspec)) continue;
              for (int j = 0; j < ancloop1.occurrence_count(); ++j) 
                parent1_parentlabels.push_back(uf[fdim+1].get_label(afspec.faceind()));
            }

            parent2_parentlabels.resize(0);
            
            for (Brep::Ancestor_Lookup::Loop_over_ancestors ancloop2(anc_lookup, parent2);
            ancloop2.notdone();
            ++ancloop2) {
              Brep::Face_Spec afspec = ancloop2.ancestor_fspec();
              if (afspec.fdim() != fdim + 2) continue;
              if (ok_to_delete(afspec)) continue;
              for (int j = 0; j < ancloop2.occurrence_count(); ++j) 
                parent2_parentlabels.push_back(uf[fdim+1].get_label(afspec.faceind()));
            }
/*
            {
              ctlog << " parent labels of parent1 = ";
              for (int jj = 0; jj < parent1_parentlabels.size(); ++jj) {
                ctlog << parent1_parentlabels[jj] << ' ';
              }
              ctlog << " of parent2 = ";
              for (int jj2 = 0; jj2 < parent2_parentlabels.size(); ++jj2) {
                ctlog << parent2_parentlabels[jj2] << ' ';
              }
              ctlog << '\n';
            }
            */


            if (parent1_parentlabels.size() != parent2_parentlabels.size()) {
              ok_to_delete(fspec) = false;
              // ctlog << " not deleting; different length parent labels\n";
              continue;
            }

            sort(parent1_parentlabels.begin(), parent1_parentlabels.end());
            sort(parent2_parentlabels.begin(), parent2_parentlabels.end());
            for (int jj = 0; jj < parent1_parentlabels.size(); ++jj) {
              if (parent1_parentlabels[jj] != parent2_parentlabels[jj]) {
                ok_to_delete(fspec) = false;
                // ctlog << " not deleting; different parents\n";
                break;
              }
            }
            if (!ok_to_delete(fspec)) continue;


            // If the number of ancestors is 2, loop over patches of
            // the face to see if the parent patches have sufficiently
            // close normals..


            for (Brep::Loop_over_patches_of_face patchloop(b, fspec);
            patchloop.notdone();
            ++patchloop) {
              Brep::PatchIndex ind = patchloop.index();
              PatchTable::Index patchind = 
                lookup_patchind[PatchInBrepAddress(fspec,ind)];
              int numends = fdim + 1;
              for (int whichend = 0; whichend < numends; ++whichend) {
                // Find the two ancestors of this patch.
                int pt_anc_count = 0;
                Point prev_normal_or_tan;
                Point param;
                param[0] = static_cast<double>(whichend);
                Point paramdirec;
                paramdirec[0] = 1.0;
                for (PatchTable::Loop_over_parents parloop(patchtable, patchind);
                parloop.notdone();
                ++parloop) {
                  PatchTable::Index ppatchind = parloop.parent_index();
                  Brep::PatchIndex thisparentind = patchtable.orig_index(ppatchind);
                  Brep::Face_Spec thisparent = patchtable.owner(ppatchind);
                  if (patchtable.gdim(ppatchind) != fdim + 1 ||
                    thisparent.fdim() != fdim + 1)
                    throw_error("patch loop inconsistency");
                  if (ok_to_delete(thisparent)) continue;
                  if (thisparent != parent1 && thisparent != parent2)
                    throw_error("patch loop problem 2");
                  // Lift the coordinates, and compute the normal.
                  // Do this differently depending on the type of patch.
                  int lift0 = parloop.lift_index(0);
                  int lift1 = (fdim == 1)? parloop.lift_index(1) : 0;
                  Point liftparam = patchtable.lift_parameters(param, 
                    patchind, parloop.listind());
                  Point normal_or_tan;
                  if (fdim == 0) {
                    normal_or_tan = b.patchmath(thisparent, thisparentind).
                      direc_deriv(liftparam, paramdirec);
                    if (normal_or_tan.l2norm() == 0)
                      throw_error("Vanishing tangent on curve");
                    normal_or_tan.normalize();
                  }
                  else {
                    normal_or_tan = b.patchmath(thisparent, thisparentind).                        
                      normal(liftparam);
                  }

                  // Flip sign if necessary.

                  if (lookup_orientation_[lift0*10+lift1] != (pt_anc_count == 0)) {
                    for (int j = 0; j < di; ++j)
                      normal_or_tan[j] = -normal_or_tan[j];
                  }
                  if (pt_anc_count == 0)
                    prev_normal_or_tan = normal_or_tan;
                  else {
                    if (!twovec_ok_for_curve_bound(prev_normal_or_tan, 
                      normal_or_tan, allowable_curvature))
                      ok_to_delete(fspec) = false;
                  }
                  /*
                  ctlog << " normal on patch " << patchtable.orig_index(patchind) 
                    <<  " end " << whichend
                    << " parent patch " << thisparent << ':' << thisparentind
                    << " (count = " << pt_anc_count << ") is"
                    << normal_or_tan << '\n';
                   
                  if (!ok_to_delete(fspec))
                    ctlog << " not ok to delete --curvature\n";
                    */

                  ++pt_anc_count;
                }
                if (pt_anc_count != 2) {
                  throw_error("Patch table problem 3");
                }
                if (!ok_to_delete(fspec)) break;
              }
              if (!ok_to_delete(fspec)) break;
            }
          }
          else { // anc_count == 0

            // ctlog << " no parents\n";
            if (fdim == 0 && di == 3) {
              // In this special case, must loop over patches 
              // containing this node to check if they
              // line up.

              init_maxmin_();
              PatchTable::Index patchind = 
                lookup_patchind[PatchInBrepAddress(fspec,0)];
              int parcount1 = 0;
              int parcount2 = 0;
              for (PatchTable::Loop_over_parents parloop(patchtable, patchind);
              parloop.notdone();
              ++parloop) {
                ++parcount1;
                if (parloop.listind() > 0) continue;
                PatchTable::Index firstcurve = parloop.parent_index();
                PatchTable::Index prev2D = -1;
                PatchTable::Index currcurve = firstcurve;
                int currendp = parloop.lift_index(0);
                already_seen.reset(My_bool(false));

                while (!already_seen[currcurve]) {
/*
                  ctlog << " looping over parent patch " << currcurve
                    << " dim " << patchtable.gdim(currcurve) << " owner "
                    << patchtable.owner(currcurve);
                   
                  if (patchtable.owner(currcurve).fdim() == patchtable.gdim(currcurve))
                    ctlog << ':' << patchtable.orig_index(currcurve);
                  ctlog << '\n';
                  */

                  // Loop over parents of currcurve; find a parent
                  // not equal to prev2D;
                  
                  already_seen[currcurve] = true;
                  int parcount3 = 0;
                  PatchTable::Index nextcurve = -1;
                  PatchTable::Index nextprev2d;
                  int nextendp;
                  for (PatchTable::Loop_over_parents parloop2(patchtable, currcurve);
                  parloop2.notdone();
                  ++parloop2) {
                    ++parcount3;
                    PatchTable::Index gpindex = parloop2.parent_index();
                    if (gpindex == prev2D)
                      continue;
                    for (PatchTable::Loop_over_children chloop(patchtable, gpindex);
                    chloop.notdone();
                    ++chloop) {
                      PatchTable::Index uncleindex = chloop.child_index();
                      if (uncleindex == currcurve) continue;
                      for (PatchTable::Loop_over_children chloop2(patchtable, uncleindex);
                      chloop2.notdone();
                      ++chloop2) {
                        if (chloop2.child_index() == patchind) {
                          nextprev2d = gpindex;
                          nextcurve = uncleindex;
                          nextendp = chloop2.lift_index(0);
                        }
                      }
                    }
                  }
                  if (parcount3 != 2) {
                    ok_to_delete(fspec) = false;
                    // ctlog << " parent patch does not have two parents; not OK to delete\n"; 
                    break;
                  }
                  if (nextcurve < 0)
                    throw_error("Nextcurv not found?");
                  Point paramcoord, direc;
                  paramcoord[0] = static_cast<double>(currendp);
                  direc[0] = (currendp)? -1.0 : 1.0;
                  Point tan1 = patchtable.patchmath(currcurve).
                    direc_deriv(paramcoord, direc);
                  paramcoord[0] = static_cast<double>(nextendp);
                  direc[0] = (nextendp)? -1.0 : 1.0;
                  Point tan2 = patchtable.patchmath(nextcurve).
                    direc_deriv(paramcoord, direc);
                  Point nrml = Point::cross_product(tan1, tan2);
                  if (nrml.l2norm() == 0)
                    throw_error("Zero normal?");
                  nrml.normalize();
                  // ctlog << " crossprod with nextcurve " << nextcurve << " tan yields " << nrml
                  //  << '\n';
                  bool ok1 = update_maxmin_(nrml, allowable_curvature);
                  if (!ok1) {
                    ok_to_delete(fspec) = false;
                    // ctlog << " curvature not met; not OK to delete\n";
                    break;
                  }
                  currcurve = nextcurve;
                  currendp = nextendp;
                  prev2D = nextprev2d;
                  ++parcount2;
                }
                if (currcurve != firstcurve) {
                  // ctlog << " didn't return to first patch currcurve = " << currcurve
                  //  << " firstcurve = " << firstcurve << "; not ok to delete\n";
                  ok_to_delete(fspec) = false;
                }
                if (!ok_to_delete(fspec)) break;
              }
              if (parcount1 != parcount2)
                ok_to_delete(fspec) = false;
            }
          }
          if (!ok_to_delete(fspec)) continue;
          // ctlog << "OK to delete " << fspec;
          if (anc_count == 2) {
            uf[fdim + 1].merge(parent1.faceind(), parent2.faceind());
            // ctlog << " merged " << parent1 << "," << parent2 << " under label "
            //  << uf[fdim+1].get_label(parent1.faceind());
          }
          // ctlog << '\n';
        }
      }
    }



    // uf gives numbers to indicate which faces are merged.
    // Now make a map from old faces to new.
    // We renumber the uf labels to skip over deleted faces.


    Brep::Face_Labeling<Brep::FaceIndex> old_to_new_map(b,0);
    vector<int> new_level_size(di + 1);

    {
      vector<Brep::FaceIndex> renumuf;
      for (int fdim = 0; fdim <= di; ++fdim) {
        int numlabels = uf[fdim].pack_labels();
        renumuf.resize(numlabels);
        for (int j = 0; j < numlabels; ++j)
          renumuf[j] = -1;
        Brep::FaceIndex nextnewfaceind = 0;
        for (Brep::Face_Spec_Loop_Over_Faces_Of_Dim fspec(b, fdim);
        fspec.notdone();