qpatchtab.cpp

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

// ------------------------------------------------------------------
// qpatchtab.cpp
//
// This file contains routines for the patch table in mesh generation.
// ------------------------------------------------------------------
// 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 "qpatchtab.h"
#include "qancestor.h"
#include "qmatvec.h"
#include "qerr.h"


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

#ifdef DEBUGGING
  ostream* debug_str;
#endif

}


// Static member function to apply three-d hh transform
// to a vector, and then shift it over.  Used to project
// from 3D to 2D.

void QMG::MG::PatchTable::apply_3d_hh_and_shift(Real beta, 
                                                const Point& hhtransform,
                                                Point& vec) {
  Real ip = 0.0;
  {
    for (int j = 0; j < 3; ++j) {
      ip += hhtransform[j] * vec[j];
      }
  }
  {
    for (int j = 1; j < 3; ++j) {
      vec[j - 1] = vec[j] - beta * ip * hhtransform[j];
    }
  }
  vec[2] = 0.0;
}



// member function get_outward_tan_in_parent computes a tangent
// to the specified parent patch of the given patch, that
// points outward from the specified point on the patch.

QMG::Point 
QMG::MG::PatchTable::get_outward_tan_in_parent(const Point& paramcoord,
                                               Index patchind, 
                                               Loop_over_parents& loop_in_progress) const {

  Index parent_index = loop_in_progress.parent_index();
  Point liftcoords = lift_parameters(paramcoord, patchind, loop_in_progress.listind());
  Point paramdirec;
  if (di_ == 2) {
    if (loop_in_progress.lift_index(0) == 0)
      paramdirec[0] = 1.0;
    else
      paramdirec[0] = -1.0;
  } 
  else { // di_ == 3

    int lift0 = loop_in_progress.lift_index(0);
    int lift1 = loop_in_progress.lift_index(1);
    if (patches_[parent_index].ptype == BEZIER_TRIANGLE) {
      if ((lift0 == 2 || lift0 == 3) && (lift1 == 0)) {
        paramdirec[0] = -0.5;
        paramdirec[1] = 1.0;
      }
      else if (lift0 == 0 && (lift1 == 2 || lift1 == 3)) {
        paramdirec[0] = 1.0;
        paramdirec[1] = -0.5;
      }
      else if (lift0 + lift1 == 5) {
        paramdirec[0] = -0.5;
        paramdirec[1] = -0.5;
      }
      else {
        throw_error("Unknown lift index for triangular patch");
      }
    }
    else { //patchtype == BEZIER_QUAD
      if (lift0 == 0) {
        paramdirec[0] = 1.0;
      }
      else if (lift0 == 1) {
        paramdirec[0] = -1.0;
      }
      else {
        paramdirec[0] = 0.0;
      }
      if (lift1 == 0) {
        paramdirec[1] = 1.0;
      }
      else if (lift1 == 1) {
        paramdirec[1] = -1.0;
      }
      else {
        paramdirec[1] = 0.0;
      }
    }
  }
  return patchmath(parent_index).direc_deriv(liftcoords, paramdirec);
}


// member function vtxnbd_patchcone_bdry_tans_
// takes a 2D patch containing a 0D patch, and
// computes the tangents to the two boundary
// edges at the point.

#ifndef BUG_IN_USING
using QMG::pair;
#endif

using QMG::Triple;


Triple<QMG::Point,QMG::Point,QMG::Point>
QMG::MG::PatchTable::
vtxnbd_patchcone_bdry_tans_(Index patchind,
                            Index grandparent_index) const {
  Point bdry_tan_1, bdry_tan_2, gpliftcoord, bdry_tan;
  Point paramdirec, parent_paramcoord;
  bool first_found = false;
  bool second_found = false;
  for (Loop_over_children loop3(*this, grandparent_index);
       loop3.notdone();
       ++loop3) {
    Index uncle_index = loop3.child_index();
    for (Loop_over_children loop4(*this, uncle_index);
         loop4.notdone();
         ++loop4) {
      Index grandchild_index = loop4.child_index();
      if (grandchild_index == patchind) {
        
        int liftind = loop4.lift_index(0);
        if (liftind == 0) 
          paramdirec[0] = 1.0;
        else
          paramdirec[0] = -1.0;
        parent_paramcoord[0] = static_cast<double>(liftind);
        bdry_tan = patchmath(uncle_index).direc_deriv(parent_paramcoord, paramdirec);
        if (!first_found) {
          bdry_tan_1 = bdry_tan;
          first_found = true;
        }
        else if (!second_found) {
          bdry_tan_2 = bdry_tan;
          second_found = true;
          break;
        }
        gpliftcoord = lift_parameters_from_childloop(parent_paramcoord,
                                                     grandparent_index,
                                                     loop3.listind());
        
      }
    }
    if (first_found && second_found) break;
  }
  if (!first_found || !second_found)
    throw_error("Parent/child patch inconsistency");
  bdry_tan_1.normalize();
  bdry_tan_2.normalize();
  return Triple<Point, Point,Point>(bdry_tan_1, bdry_tan_2, gpliftcoord);
}


// member function vtxnbd_select_patchcone_randpoint
// takes a 2D patch containing a 0D patch, and
// selects a random point in the cone defined by the 
// the patch (extending the boundary tangents
// from the point out to infinity).
// Only used for 3D.

QMG::Point
QMG::MG::PatchTable::
vtxnbd_select_patchcone_randpoint(Index patchind,
                                  Index ancestor_index) const {

  Point endpoint2;
  if (gdim(ancestor_index) == 2) {
    Real wt1 = 0.2 + 0.6 * random_real();
    Real wt2 = 1.0 - wt1;
    Triple<Point, Point, Point> rval =
      vtxnbd_patchcone_bdry_tans_(patchind, ancestor_index);
    for (int j = 0; j < 3; ++j)
      endpoint2[j] = wt1 * rval.first[j] + wt2 * rval.second[j];
  }
  else {
    bool found = false;
    Point parent_paramcoord, paramdirec;
    for (Loop_over_children loop3(*this, ancestor_index);
         loop3.notdone();
         ++loop3) {
      if (loop3.child_index() == patchind) {
        int liftind = loop3.lift_index(0);
        if (liftind == 0) 
          paramdirec[0] = 1.0;
        else
          paramdirec[0] = -1.0;
        parent_paramcoord[0] = static_cast<double>(liftind);
        endpoint2 = 
          patchmath(ancestor_index).direc_deriv(parent_paramcoord, paramdirec);
        found = true;
        break;
      }
    }
    if (!found)
      throw_error("Illegal arguments to vtxnbd_select_patchcone_randpoint");
  }
  return endpoint2;
}


// member function vtxnbd_patchcone_meets_seg
// takes a 2D patch containing a 0D patch, and
// computes the intersection of the
// the patch with a segment.  It returns a code and the distance.

pair<QMG::MG::PatchTable::ConeIntersectCode, QMG::Real>
QMG::MG::PatchTable::
vtxnbd_patchcone_meets_seg(Index patchind,
                           Index grandparent_index,
                           const Point& endpoint1,
                           const Point& endpoint2,
                           const Point& segtan,
                           double degen_cutoff) const {

  // Get the two boundary tangents.

  Triple<Point, Point, Point> rval =
    vtxnbd_patchcone_bdry_tans_(patchind, grandparent_index);
  Point gpnormal = 
    patchmath(grandparent_index).normal(rval.third);
  Real denom = Point::inner_product(segtan, gpnormal);
  Real numer = Point::inner_product(endpoint1,gpnormal);
  if (fabs(numer) < degen_cutoff || fabs(denom) < degen_cutoff) {
    return pair<ConeIntersectCode, Real>(DEGEN, 0.0);
  }
  Real newdist = numer / denom;
  if (newdist < 0)
    return pair<ConeIntersectCode, Real>(MISSES_PATCH, 0.0);
  
  // Check if seg goes thru patch here
  
  Point ipoint;
  {
    for (int j = 0; j < 3; ++j) {
      ipoint[j] = (1.0 - newdist) * endpoint1[j] + newdist * endpoint2[j];
    }
  }
  
  Point cross1 = Point::cross_product(rval.first, ipoint);
  Point cross2 = Point::cross_product(ipoint, rval.second);
  Point cross3 = Point::cross_product(rval.first, rval.second);
  Real ip1 = Point::inner_product(cross1, cross3);
  Real ip2 = Point::inner_product(cross2, cross3);
  if (fabs(ip1) < degen_cutoff) {
    Real ip3 = Point::inner_product(ipoint, rval.second);
    if (ip3 > 0)
      return pair<ConeIntersectCode, Real>(HITS_EDGE, newdist);
    else 
      return pair<ConeIntersectCode, Real>(MISSES_PATCH, 0.0);
  }
  if (ip1 < 0) {
    return pair<ConeIntersectCode, Real>(MISSES_PATCH, 0.0);
  }
  if (fabs(ip2) < degen_cutoff) {
    Real ip3 = Point::inner_product(ipoint, rval.first);
    if (ip3 > 0)
      return pair<ConeIntersectCode, Real>(HITS_EDGE, newdist);
    else 
      return pair<ConeIntersectCode, Real>(MISSES_PATCH, 0.0);
  }
  if (ip2 < 0) {
    return pair<ConeIntersectCode, Real>(MISSES_PATCH, 0.0);
  }
  if (denom > 0)
    return pair<ConeIntersectCode, Real>(HITS_PATCH_FRONT, newdist);
  return pair<ConeIntersectCode, Real>(HITS_PATCH_BACK, newdist);
}


// member function determine_ray_front_face takes a ray
// coming out of a patch and determines which side of a full-dimensional
// patch sees the ray.  This is complicated if the input patch
// is low dimensional, because then we have to determine
// sidedness from a computation involving the
// patch's ancestors.


pair<QMG::MG::PatchTable::Index, QMG::MG::PatchTable::ConeIntersectCode>
QMG::MG::PatchTable::determine_ray_front_face(Index patchind, 
                                              const Point& param_coord,
                                              const Point& ray_direction,
                                              double tol) const {
#ifdef DEBUGGING
  bool trace_on = false;

  if (trace_on) {
    *debug_str << "Reached determine ray_front_face patchind = pat#" << patchind
      << "#\nray_direction = " << ray_direction << " param_coord = ";
    for (int j = 0; j < this -> gdim(patchind); ++j)
      *debug_str << param_coord[j]  << " ";
    *debug_str << "\n";
  }
#endif

  typedef pair<Index, ConeIntersectCode> ReturnvalType;

  Real degen_cutoff = sqrt(tol);

  Point normalized_direction = ray_direction;
  normalized_direction.normalize();
  Index result_patchind;
  ConeIntersectCode result_code;
  bool prev_empty;
  int gdim1 = gdim(patchind);


  // If the patch has dimension di-1, then the side is determined by
  // whether the inner product of the normal to the patch with the ray
  // is positive or negative.

  if (gdim1 == di_ - 1) {
    Point normal = patchmath(patchind).normal(param_coord);
    Real ip = Point::inner_product(normal, ray_direction);


#ifdef DEBUGGING
    if (trace_on) {
      *debug_str << "Case 1, normal = " << normal << " ip = " << ip << "\n";
    }
#endif

    if (fabs(ip) < degen_cutoff) {
#ifdef DEBUGGING
      if (trace_on)
        *debug_str << " returning degen\n";
#endif
      return ReturnvalType(0,DEGEN);
    }
    else if (ip > 0) {
#ifdef DEBUGGING
      if (trace_on)
        *debug_str << " returning " << front_side(patchind) << "\n";
#endif
      return ReturnvalType(patchind, HITS_PATCH_FRONT);
    }
    else {
#ifdef DEBUGGING
      if (trace_on)
        *debug_str << " returning " << back_side(patchind) << "\n";
#endif
      return ReturnvalType(patchind, HITS_PATCH_BACK);
    }
  }

  // If we failed that first test but the patch has no parents, then
  // it is an isolated internal boundary and we already classified the
  // region containing it.

  if (patches_[patchind].num_parents == 0) {
#ifdef DEBUGGING
    if (trace_on) 
      *debug_str << "Case 2 returning " << front_side(patchind) << "\n";
#endif
    return ReturnvalType(patchind, HITS_PATCH_FRONT);
  }

  // If the patch is dimension di-2 then it is either a segment between
  // two or more patches in R^3 or a point between two or more curves
  // in R^2.  In either case, we determine sidedness by finding the closest
  // counterclockwise patch of higher dimension to the ray, and using the