qseparation.cpp

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

// ------------------------------------------------------------------
// qseparation.cpp
//
// This file contains routines for separation on boxes in the
// mesh generator.
// ------------------------------------------------------------------
// Copyright (c) 1999 by Cornell University.  All rights reserved.
// 
// See the accompanying file 'Copyright' for authorship information,
// the terms of the license governing this software, and disclaimers
// concerning this software.
// ------------------------------------------------------------------
// This file is part of the QMG software.  
// Version 2.0 of QMG, release date September 3, 1999
// ------------------------------------------------------------------



#ifdef _MSC_VER
#if _MSC_VER == 1200
#include "qmatvec.h"
#endif
#endif

#include "qseparation.h"


namespace {
  using namespace QMG;
  using namespace QMG::MG;

  const int NUM_TEST_DIRECTIONS_2D = 4;
  const Real test_direction_init[] = {0,1,0, 1,0,0, 1,1,0, 1,-1,0,
                                      0,0,1,  0,1,1, 1,0,1, 1,1,1,
                                      0,-1,1, -1,0,1, 1,-1,1, -1,1,1,
                                      -1,-1,1};
}

// This routine tests if a box is crowded.  Returns a bool (is it crowded?) and the face
// of the apex if not.  If it is crowded, returns two interfering faces.

using QMG::Triple;


Triple<QMG::Brep::Face_Spec, QMG::Brep::Face_Spec, bool> 
QMG::MG::Separation_Functor::is_crowded_(const ActiveBox& b,
                                         int phase) {

  Brep::Face_Spec apex1(-1, -1);
  Brep::Face_Spec apex2(-1, -1);
  {
    for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
      PatchTable::Index patchind = b.patch(relnum).patchind();
      Brep::Face_Spec owner = patchtable_.owner(patchind);
      // If there is a brep face of dimension lower than phase, then
      // the box is crowded.
      if (owner.fdim() < phase)
        return Triple<Brep::Face_Spec, Brep::Face_Spec, bool>(owner, owner, true);
      if (apex1.fdim() == -1 || apex1.fdim() > owner.fdim())
        apex1 = owner;
      if (apex1.fdim() == owner.fdim() && apex1 != owner)
        apex2 = owner;
    }
  }
  // If there is no apex, or if the apex dimension is higher than the phase,
  // then the box is not crowded.
  if (apex1.fdim() == -1 || apex1.fdim() > phase)
    return Triple<Brep::Face_Spec, Brep::Face_Spec, bool>(apex1, apex1, false);
  // If two separate apexes, box is crowded.
  if (apex1.fdim() == apex2.fdim())
    return Triple<Brep::Face_Spec, Brep::Face_Spec, bool>(apex1, apex2, true);
  {
    for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
      PatchTable::Index patchind = b.patch(relnum).patchind();
      Brep::Face_Spec owner = patchtable_.owner(patchind);
      bool owner_is_ancestor = false;
      for (Brep::Ancestor_Lookup::Loop_over_ancestors ancloop(lca_, apex1);
      ancloop.notdone();
      ++ancloop) {
        if (ancloop.ancestor_fspec() == owner) {
          owner_is_ancestor = true;
          break;
        }
      }
      if (!owner_is_ancestor)
        return Triple<Brep::Face_Spec, Brep::Face_Spec, bool>(apex1, owner, true);
    }
  }
  return Triple<Brep::Face_Spec, Brep::Face_Spec, bool>(apex1, apex1, false);
}

  
  
  // ------------------------------------------------------------------
  // Tests if a completely interior box is OK for the size function.
  // Evaluates the size function at the box vertices.

QMG::MG::Separation_Functor::Size_Return 
QMG::MG::Separation_Functor::
ok_for_size_function_completely_interior_(const ActiveBox& b,
                                          const Brep::Face_Spec& apex) {
  
  Real boxwidth = b.real_width();
  Size_Return ret;
  for (ActiveBox::Loop_over_vertices it(b);
  it.notdone();
  ++it) {
    const Point& incpoint = it.realpoint();
    pair<Real,bool> rval = size_control_.eval_interior(incpoint, apex);
    Real rqsize = rval.first;
    if (rqsize < boxwidth) {
      ret.val = rqsize;
      ret.fspec = apex;
      ret.ok = false;
      return ret;
    }
    if (rval.second) {
      ret.ok = true;
      return ret;
    }
  }
  ret.ok = true;
  return ret;
}
  
  // ------------------------------------------------------------------
  // Checks if a box is OK for the size function.  Evaluates the size
  // function at all the incidences in the box.

QMG::MG::Separation_Functor::Size_Return 
QMG::MG::Separation_Functor::ok_for_size_function_(const ActiveBox& b, 
                                                   const Brep::Face_Spec& apex) {
  
  
  Real boxwidth = b.real_width();
  Point incpoint_real;
  Point incpoint_par;
  // *debug_str << "in ok boxwidth = " << boxwidth << '\n';
  
  Size_Return ret;
  for (int relnum = 0; relnum < b.num_patch(); relnum++) {
    Patch_in_ActiveBox patch = b.patch(relnum);
    PatchTable::Index patchind = patch.patchind();
    int numinc = patch.num_incidence();
    for (int ii = 0; ii < numinc; ii++) {
      IncidenceTable::Index inc = patch.incidence(ii);
      incpoint_real = inc_table_.real_coord(inc);
      incpoint_par = inc_table_.param_coord(inc);
      pair<Real,bool> rval = 
        size_control_.eval(incpoint_real, incpoint_par, patchind);
      Real rqsize = rval.first;
      // *debug_str << " rqsize = " << rqsize << '\n';
      if (rqsize < boxwidth) {
        // *debug_str << " return false\n";
        ret.val = rqsize;
        ret.fspec = apex;
        ret.ok = false;
        return ret;
      }
      if (rval.second) {
        // *debug_str << " return true 1\n";
        ret.ok = true;
        return ret;
      }
    }
  }
  // *debug_str << "return true 2\n";
  ret.ok = true;
  return ret;
}




// ------------------------------------------------------------------
// helper function for ok_for_curvature.  For a particular incidence,
// updates the normal_variation and param_variation arrays to 
// keep track of how much the normal and parametric coordinates
// vary over the box.
// Results stored in track_min_max_.

void 
QMG::MG::Separation_Functor::process_incidence0_(const Brep::Face_Spec& owner,
                                                 const Point& normal_or_tan) {
  if (logstr_.verbosity() >= 9)
    logstr_.str() << " normal est owner = " << owner << " nrml = " << normal_or_tan << '\n';
  
  for (int l = 0; l < num_test_directions_actual_; ++l) {
    Real ip = Point::inner_product(test_directions_[l], normal_or_tan);
    if (ip < track_max_min_(owner).lower_bound[l])
      track_max_min_(owner).lower_bound[l] = ip;
    if (ip > track_max_min_(owner).upper_bound[l])
      track_max_min_(owner).upper_bound[l] = ip;
  }
}


// Assuming track_max_min is already initialized (by the previous routine)
// this routine checks whether a new vector is valid for track_max_min_.

bool
QMG::MG::Separation_Functor::recheck_normal_(const Brep::Face_Spec& owner,
                                             const Point& normal_or_tan) {
  Real rq1 = rqvariance_(owner);
  for (int l = 0; l < num_test_directions_actual_; ++l) {
    Real ip = Point::inner_product(test_directions_[l], normal_or_tan);
    Real mtp = (track_max_min_(owner).lower_bound[l] + 
      track_max_min_(owner).upper_bound[l]) / 2.0;
    if (fabs(ip - mtp) > rq1)
      return false;
  }
  return true;
}




// For an incidence in a box, figure out the normal 
// and track it in track_min_max.

void 
QMG::MG::Separation_Functor::
process_incidence1_(PatchTable::Index patchind,
                    const Point& paramcoord,
                    const Point& realcoord) {  
  if (logstr_.verbosity() >= 10) {
    logstr_.str() << " reached process incidence1 pat#" << patchind << "# paramcoord = "
      << paramcoord << " realcoord = " << realcoord << '\n';
  }
  
  Brep::Face_Spec owner = patchtable_.owner(patchind);
  int pdim = patchtable_.gdim(patchind);
  if (owner.fdim() == 0 || pdim != owner.fdim())
    return;


  Point normal_or_tan;
  if (pdim == 2) {
    normal_or_tan = 
      patchtable_.patchmath(patchind).normal(paramcoord);
  }
  else {
    Point param_direc;
    param_direc[0] = 1.0;
    normal_or_tan = patchtable_.patchmath(patchind).
      direc_deriv(paramcoord, param_direc);
    normal_or_tan.normalize();
    if (!patchtable_.is_forward(patchind)) {
      for (int j = 0; j < di_; ++j)
        normal_or_tan[j] = -normal_or_tan[j];
    }
  }
  if (logstr_.verbosity() >= 10) {
    logstr_.str() << "owner = " << owner
      << " normal from direct eval = " << normal_or_tan << '\n';
  }
  process_incidence0_(owner, normal_or_tan);
}



// In ok_for_curvature_ we do a search on patches to see
// if the patches in a box are connected.  This routine
// determines whether a patch may be used for extending the search.

bool 
QMG::MG::Separation_Functor::acceptable_for_search_(PatchTable::Index patchind,
                                                    const Point& bmidpoint,
                                                    Real bwidth) {
  if (logstr_.verbosity() >= 13) {
    logstr_.str() << "Reached acceptable_for_search_ pat#"
      << patchind << " bmidpoint = " << bmidpoint << " bwidth = " << bwidth << '\n';
  }
  
  if (acceptable_for_search_cache_[patchind] >= 0) {
    if (logstr_.verbosity() >= 13) {
      logstr_.str() << " returning cached value " << acceptable_for_search_cache_[patchind]
        << '\n';
    }
    return (acceptable_for_search_cache_[patchind] > 0);
  }
    
  // If the patch is in the box, it's OK.
  
  if (map_to_relnum_[patchind] >= 0) {
    acceptable_for_search_cache_[patchind] = 1;
    if (logstr_.verbosity() >= 13) {
      logstr_.str() << " is in box, return true\n";
    }
    return true;
  }
  
  int gdim1 = patchtable_.gdim(patchind);
  
  // Can search only on 0 or 1-dim patches.
  
  if (gdim1 > 1) {
    if (logstr_.verbosity() >= 13) {
      logstr_.str() << " is hi dim, return false\n";
    }
    acceptable_for_search_cache_[patchind] = 0;
    return false;
  }
  // 1-dml patches always OK.
  
  if (gdim1 == 1) {    
    if (logstr_.verbosity() >= 13) {
      logstr_.str() << " is dim 1, return false\n";
    }
    acceptable_for_search_cache_[patchind] = 1;
    return true;
  }
  
  // The only case left is the patch is a vertex and outside b itself.
  // Check if we can continue searching from a vertex
  // checking (1) whether it is in a larger box around b and (2)
  // whether its normals satisfy the rqvariance_ requirement.
  
  Point paramcoord;
  
  Point realcoord = 
    patchtable_.patchmath(patchind).get_real_point(paramcoord);
  
  for (int j = 0; j < di_; ++j) {
    if (fabs(realcoord[j] - bmidpoint[j]) > 2 * bwidth) {   
      if (logstr_.verbosity() >= 13) {
        logstr_.str() << " is outside box realcoord = " << realcoord << " return false\n";
      }
      acceptable_for_search_cache_[patchind] = 0;
      return false;
    }
  }
  
  
  Point normal_or_tan, param_direc;
  param_direc[0] = 1;
  
  for (PatchTable::Loop_over_parents ploop1(patchtable_, patchind);
  ploop1.notdone();
  ++ploop1) {
    PatchTable::Index ppatchind = ploop1.parent_index();
    Point param1 = patchtable_.lift_parameters(paramcoord, patchind, ploop1.listind());
    Brep::Face_Spec eowner = patchtable_.owner(ppatchind); 
    if (eowner.fdim() == 1) {
      normal_or_tan = patchtable_.patchmath(ppatchind).
        direc_deriv(param1, param_direc);
      normal_or_tan.normalize();
      if (!patchtable_.is_forward(patchind)) {
        for (int j = 0; j < di_; ++j)
          normal_or_tan[j] = -normal_or_tan[j];
      }
      if (!recheck_normal_(eowner, normal_or_tan)) {
        if (logstr_.verbosity() >= 13) {
          logstr_.str() << " failed to confirm normal on parent pat#" << ppatchind
            << "# owner " << eowner << " return false\n";
        }
        acceptable_for_search_cache_[patchind] = 0;
        return false;
      }
    }
    for (PatchTable::Loop_over_parents ploop2(patchtable_, ppatchind);
    ploop2.notdone();
    ++ploop2) {
      PatchTable::Index gppatchind = ploop2.parent_index();
      Point param2 = patchtable_.lift_parameters(param1, ppatchind, ploop2.listind());
      Brep::Face_Spec fowner = patchtable_.owner(gppatchind);
      normal_or_tan = patchtable_.patchmath(gppatchind).normal(param2);
      if (!recheck_normal_(fowner, normal_or_tan)) {
        if (logstr_.verbosity() >= 13) {
          logstr_.str() << " failed to confirm normal on grandparent pat#" << gppatchind
            << "# owner " << fowner << " return false\n";
        }
        acceptable_for_search_cache_[patchind] = 0;
        return false;
      }
    }
  }
  if (logstr_.verbosity() >= 13) {
    logstr_.str() << "Ran thru tests, all passed\n";
  }
  acceptable_for_search_cache_[patchind] = 1;
  return true;
}


  

// ------------------------------------------------------------------
// function determines whether the curvature of patches in a box
// is excessive.


QMG::MG::Separation_Functor::Size_Return 
QMG::MG::Separation_Functor::
ok_for_curvature_(const ActiveBox& b) {
  Real boxwidth = b.real_width();
  int oldverb = logstr_.verbosity();
  /*

  IntCoord p1 = b.boxaddress().lowerleft();
  if (p1[0] == 0x20000000U && p1[1] == 0x18000000U) {
    logstr_.str() << " reached special box";
    b.full_dump2(logstr_.str(), patchtable_, inc_table_);
    logstr_.change_verbosity(15);
  }
  */



  Point normal;
  if (logstr_.verbosity() >= 8)
    logstr_.str() << "In ok for curvature, b = " << b << '\n';

  // Initialize the rqvariance field for each face.

  {
    for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
      Patch_in_ActiveBox patch = b.patch(relnum);
      PatchTable::Index patchind = patch.patchind();      
      Brep::Face_Spec owner = patchtable_.owner(patchind);
      rqvariance_(owner) = BIG_REAL;
    }