ptloc.cpp

一个2D电磁场FEM计算的VC++源程序

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// Ptloc.cpp    Implement fast point location routines.                    //
//                                                                         //
// by Si hang Jun, 2001.                                                   //
// [e-mail]                                                        //
//                                                                         //
// This file with it's couple file ptloc.h were specially written for      //
//   the FEMM program(http://groups.yahoo.com/group/femm/) femmview to     //
//   provide a fast point location function to improve the Line Integral   //
//   calculation speed when the mesh size is big.                          //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include <afx.h>
#include <afxtempl.h>
#include "problem.h"      // declared CMeshNode, CMeshElem 

#include "femmview.h"
#include "xyplot.h"
#include "ptloc.h"
#include "femmviewDoc.h"
#include <math.h>

// A large value to indicate a impossible element index.
#define INVALIDELEINDEX    0x7fffffff
// We use the highest bit of a integer as boundary mark.
#define BOUNDFLAG          0x80000000

// some macros for convenience

#define TriEdge_Div4  >> 2
#define TriEdge_Mul4  << 2
#define TriEdge_Mod4  & 03

#define Eleindex(E)   ( (E) TriEdge_Div4 )
#define Orient(E)     ( (E) TriEdge_Mod4 )
#define Encode(T, V)  ( ((T) TriEdge_Mul4) + (V) )

// NOTE: Above 3 macros act like functions!

///////////////////////////////////////////////////////////////////////////////
// Local variables declaration                                               //
// BEGIN                                                                     //
///////////////////////////////////////////////////////////////////////////////

// These two variables should be init before use.

static CFemmviewDoc *pDoc;
// static CArray< CMeshNode, CMeshNode&> *pmeshnode;
// static CArray< CElement, CElement&>   *pmeshelem;

// Fast lookup arrays to speed some of the mesh manipulation primitives.

static int plus1mod3[3] = {1, 2, 0};
static int minus1mod3[3] = {2, 0, 1};

// Variables for fast locate point.
/*
static CArray< CMeshNode, CMeshNode&> *pDoc->pmeshnode;
static CArray< CElement, CElement&>   *pDoc->pmeshelem;
static TriEdge pDoc->recenttri;
static int pDoc->samples;
static unsigned long pDoc->randomseed; 
TriEdge pDoc->bdrylinkhead[512];
int numberofbdrylink;
*/

// Keep a recently visited triangle's index.  Improves point location
//   if proximate points are inserted sequentially.
// static TriEdge recenttri;

// Number of random pDoc->samples for point location.
// static int samples;

// Current random number seed.   
// static unsigned long randomseed;  

// These two variables are used to handle the case where there exist 
//   mutil-isolated-block in the model.  

// Variables for keeping each isolated boundary link's startedge. Accually
//   it should be replaced with a queue, but I don't want introduce too many 
//   data types here for keeping its simplicity. Change the '512' to a larger 
//   number if there exist more than 512 isolated blocks.
//TriEdge bdrylinkhead[512];

// Number of isolated blocks.
// int numberofbdrylink;

///////////////////////////////////////////////////////////////////////////////
// END                                                                       //
// Local variables declaration                                               //
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
// class TriEdge implementation                                              //
// BEGIN                                                                     //
///////////////////////////////////////////////////////////////////////////////

TriEdge::TriEdge() : ele(INVALIDELEINDEX), ori(0)
{
}

TriEdge::TriEdge(int _ele, int _ori) : ele(_ele), ori(_ori)
{
}

TriEdge::TriEdge(TriEdge& src) : ele(src.ele), ori(src.ori)
{
}

TriEdge& TriEdge::operator=(const TriEdge& src)
{
  ele = src.ele;
  ori = src.ori;
  return *this;
}

BOOL TriEdge::operator==(TriEdge& src)
{
  return (ele == src.ele) && (ori == src.ori);
}

BOOL TriEdge::operator!=(TriEdge& src)
{
  return (ele != src.ele) || (ori != src.ori);
}

TriEdge TriEdge::sym()
{
  int E = ((*pDoc->pmeshelem)[ele]).n[ori];
  return TriEdge(Eleindex(E), Orient(E));
}

TriEdge& TriEdge::symself()
{
  int E = ((*pDoc->pmeshelem)[ele]).n[ori];
  ele = Eleindex(E);
  ori = Orient(E);
  return *this;
}

TriEdge TriEdge::lnext()
{
  return TriEdge(ele, plus1mod3[ori]);
}

TriEdge& TriEdge::lnextself()
{
  ori = plus1mod3[ori];
  return *this;
}

TriEdge TriEdge::lprev()
{
  return TriEdge(ele, minus1mod3[ori]);
}

TriEdge& TriEdge::lprevself()
{
  ori = minus1mod3[ori];
  return *this;
}

TriEdge TriEdge::onext()
{
  return lprev().symself();
}

TriEdge& TriEdge::onextself()
{
  lprevself().symself();
  return *this;
}

TriEdge TriEdge::oprev()
{
  return sym().lnextself();
}

TriEdge& TriEdge::oprevself()
{
  symself().lnextself();
  return *this;
}

TriEdge TriEdge::dnext()
{
  return sym().lprevself();
}

TriEdge& TriEdge::dnextself()
{
  symself().lprevself();
  return *this;
}

TriEdge TriEdge::dprev()
{
  return lnext().symself();
}

TriEdge& TriEdge::dprevself()
{
  lnextself().symself();
  return *this;
}

TriEdge TriEdge::rnext()
{
  return sym().lnextself().symself();
}

TriEdge& TriEdge::rnextself()
{
  symself().lnextself().symself();
  return *this;
}

TriEdge TriEdge::rprev()
{
  return sym().lprevself().symself();
}

TriEdge& TriEdge::rprevself()
{
  symself().lprevself().symself();
  return *this;
}

CMeshNode* TriEdge::org()
{
  int N = ((*pDoc->pmeshelem)[ele]).p[plus1mod3[ori]]; 
  return &((*pDoc->pmeshnode)[N]);
}

CMeshNode* TriEdge::dest()
{
  int N = ((*pDoc->pmeshelem)[ele]).p[minus1mod3[ori]]; 
  return &((*pDoc->pmeshnode)[N]);
}

CMeshNode* TriEdge::apex()
{
  int N = ((*pDoc->pmeshelem)[ele]).p[ori];  
  return &((*pDoc->pmeshnode)[N]);
}

BOOL TriEdge::symexist()
{
  int E = ((*pDoc->pmeshelem)[ele]).n[ori];
  return !(E & BOUNDFLAG);
}

void TriEdge::bondboundary(TriEdge& tedge)
{
  // Here tedge must also be a boundary edge.
  if(tedge.symexist()) throw;
  // This edge's neigh must be Outer Space.
  int E = ((*pDoc->pmeshelem)[ele]).n[ori];
  if(!(E & BOUNDFLAG)) throw;
  // Now bound tedge as this edge's neigh.
  // Note: We also need keep the BOUNDFLAG in its neigh.
  E = Encode(tedge.ele, tedge.ori);
  E |= BOUNDFLAG;
  ((*pDoc->pmeshelem)[ele]).n[ori] = E;
}

TriEdge TriEdge::nextboundary()
{
  // This edge's neigh must be Outer Space.
  int E = ((*pDoc->pmeshelem)[ele]).n[ori];
  if(!(E & BOUNDFLAG)) throw;
  // Get the true TriEdge handle.
  // Note: We must get off the BOUNDFLAG from its neigh first.
  E = (E & ~BOUNDFLAG);
  return TriEdge(Eleindex(E), Orient(E));
}

TriEdge& TriEdge::nextboundaryself()
{
  // This edge's neigh must be Outer Space.
  int E = ((*pDoc->pmeshelem)[ele]).n[ori];
  if(!(E & BOUNDFLAG)) throw;
  // Get the true TriEdge handle.
  // Note: We must get off the BOUNDFLAG first.
  E = (E & ~BOUNDFLAG);
  ele = Eleindex(E);
  ori = Orient(E);
  return *this;
}
///////////////////////////////////////////////////////////////////////////////
// END                                                                       //
// class TriEdge implementation                                              //
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
// Geometry Predicates functions                                             //
// BEGIN                                                                     //
///////////////////////////////////////////////////////////////////////////////

// Compare 2 floating point values dx and dy.
//   return +1 if dx > dy
//   return -1 if dx < dy
//   return  0 if dx == dy       
int FuzzyComp(const double dx, const double dy)
{
  double dSum  = fabs(dx) + fabs(dy);
  double dDiff = dx - dy;
  const double dTol = 5.e-9;
  const double dFloor = 1.e-6;

  dSum = (dSum > dFloor) ? dSum : dFloor;
  if (dDiff > dTol * dSum) return (1);
  else if (dDiff < - dTol * dSum) return (-1);
  else return (0);
}

// Orientation 2D, input 3 points dA, dB, dC.
//   return +1 if they are occur in counterclockwise order.
//   return -1 if they are occur in clockwise order.
//   return  0 if they are collinear. 
int Orient2D(const double dAx, const double dAy, 
             const double dBx, const double dBy,
             const double dCx, const double dCy)
{             
  // Currently no exact arithmetic
  double dDet = ((dBx - dAx) * (dCy - dAy) - (dCx - dAx) * (dBy - dAy));
  // Scale the determinant by the mean of the magnitudes of vectors:
  double distAB = sqrt((dAx - dBx) * (dAx - dBx) + (dAy - dBy) * (dAy - dBy)),
         distBC = sqrt((dBx - dCx) * (dBx - dCx) + (dBy - dCy) * (dBy - dCy)),
         distCA = sqrt((dCx - dAx) * (dCx - dAx) + (dCy - dAy) * (dCy - dAy));
  double dScale = (distAB + distBC + distCA) / 3; 
  dDet /= (dScale * dScale);
  if (dDet > 1.e-10) return 1;
  else if (dDet < -1.e-10) return -1;
  else return 0;
}

int Orient2D(CMeshNode& dA, CMeshNode& dB, CMeshNode& dC)
{
  return Orient2D(dA.x, dA.y, dB.x, dB.y, dC.x, dC.y);
}

///////////////////////////////////////////////////////////////////////////////
// END                                                                       //
// Geometry Predicates functions                                             //
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
// Ptloc Preprocess Functions                                                //
// BEGIN                                                                     //
///////////////////////////////////////////////////////////////////////////////

void BuildOuterBoundaryLink(TriEdge& startedge)
{
  TriEdge thisedge, nextedge;

  thisedge = startedge;
  // Get startedge's next edge in the same element.
  nextedge = thisedge.lnext();  
  while(TRUE) {
    // Spin nextedge around its org until it encounter an Outer Boundary.
    while(nextedge.symexist()) { 
      nextedge.oprevself();  
    }
    // Add nextedge to link(after thisedge).
    thisedge.bondboundary(nextedge);
    // Check to see if we can leave.
    if(nextedge == startedge) break;
    thisedge = nextedge;
    nextedge = thisedge.lnext();
  }
}

void BuildElemMap(CArray< CMeshNode, CMeshNode&> *pmeshnodes,
                  CArray< CElement, CElement&> *pmeshelems,
                  int *pNumList, int **ppConList, void *pd)
{
  pDoc=(CFemmviewDoc *) pd;
  // Init variables for Point location functions.
  pDoc->recenttri.ele = INVALIDELEINDEX;
  pDoc->samples = 1; 
  pDoc->randomseed = 1;
  pDoc->numberofbdrylink = 0;
  // Init pDoc->pmeshnode, pDoc->pmeshelem.
  pDoc->pmeshnode = pmeshnodes;
  pDoc->pmeshelem = pmeshelems;
  	
  int i, j;

  // Init all elements' neigh to be unfinished.
  for(i = 0; i < pDoc->pmeshelem->GetSize(); i ++) {
    for(j = 0; j < 3; j ++)
	  ((*pDoc->pmeshelem)[i]).n[j] = INVALIDELEINDEX;
  }
  
  int orgi, desti;
  int ei, ni;
  BOOL done;
  // Loop all elements, to find and set there neighs.
  for(i = 0; i < pDoc->pmeshelem->GetSize(); i ++) {
    for(j = 0; j < 3; j ++) {
      if(((*pDoc->pmeshelem)[i]).n[j] == INVALIDELEINDEX) {
        // Get this edge's org and dest node index, 
        orgi = ((*pDoc->pmeshelem)[i]).p[plus1mod3[j]];
        desti = ((*pDoc->pmeshelem)[i]).p[minus1mod3[j]];
        done = FALSE;
        // Find this edge's neigh from the org node's list
        for(ni = 0; ni < pNumList[orgi]; ni ++) {
          // Find a Element around org node contained dest node of this edge. 
          ei = ppConList[orgi][ni];
          if(ei == i) continue; // Skip myself.
          // Check this Element's 3 vert to see if there exist dest node.
          if(((*pDoc->pmeshelem)[ei]).p[0] == desti) {
            if(((*pDoc->pmeshelem)[ei]).p[1] == orgi) {
              ((*pDoc->pmeshelem)[i]).n[j] = Encode(ei, 2);
              ((*pDoc->pmeshelem)[ei]).n[2] = Encode(i, j);
            } else if (((*pDoc->pmeshelem)[ei]).p[2] == orgi) {