ptloc.cpp

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

              ((*pDoc->pmeshelem)[i]).n[j] = Encode(ei, 1); 
              ((*pDoc->pmeshelem)[ei]).n[1] = Encode(i, j);
            } else {
              throw; // Shouldn't be here.
            }
            done = TRUE;
            break;
          } else if(((*pDoc->pmeshelem)[ei]).p[1] == desti) {
            if(((*pDoc->pmeshelem)[ei]).p[2] == orgi) {
              ((*pDoc->pmeshelem)[i]).n[j] = Encode(ei, 0);
              ((*pDoc->pmeshelem)[ei]).n[0] = Encode(i, j);
            } else if(((*pDoc->pmeshelem)[ei]).p[0] == orgi) {   
              ((*pDoc->pmeshelem)[i]).n[j] = Encode(ei, 2); 
              ((*pDoc->pmeshelem)[ei]).n[2] = Encode(i, j);
            } else {
              throw; // Shouldn't be here.
            }
            done = TRUE;
            break;
          } else if(((*pDoc->pmeshelem)[ei]).p[2] == desti) {
            if(((*pDoc->pmeshelem)[ei]).p[0] == orgi) {
              ((*pDoc->pmeshelem)[i]).n[j] = Encode(ei, 1);
              ((*pDoc->pmeshelem)[ei]).n[1] = Encode(i, j);
            } else if(((*pDoc->pmeshelem)[ei]).p[1] == orgi) {  
              ((*pDoc->pmeshelem)[i]).n[j] = Encode(ei, 0); 
              ((*pDoc->pmeshelem)[ei]).n[0] = Encode(i, j);
            } else {
              throw; // Shouldn't be here.
            }
            done = TRUE;
            break;
          } 
        }
        if(!done) {
          // This edge must be a Boundary Edge.
          ((*pDoc->pmeshelem)[i]).n[j] |= BOUNDFLAG;
        }
      } // Finish One Edge        
    } // End of One Element Loop
  } // End of Main Loop
  
  // Build links for Outer Boundary. These links used for fast point 
  //   location functions.
  for( i = 0; i < pDoc->pmeshelem->GetSize(); i ++ ) {
    for( j = 0; j < 3; j++ ) {
      if(((*pDoc->pmeshelem)[i]).n[j] & BOUNDFLAG) {
      	if((((*pDoc->pmeshelem)[i]).n[j] & ~BOUNDFLAG) == INVALIDELEINDEX) {
          TriEdge startedge(i, j);
          pDoc->bdrylinkhead[pDoc->numberofbdrylink] = startedge;
          BuildOuterBoundaryLink(startedge);
          pDoc->numberofbdrylink++;
        }  
      }
    } // End of One Element Loop
  } // End of Main Loop
}

///////////////////////////////////////////////////////////////////////////////
// END                                                                       //
// Ptloc Preprocess Functions                                                //
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  randomnation()   Generate a random number between 0 and `choices' - 1.   //
//                                                                           //
//  This is a simple linear congruential random number generator.  Hence, it //
//  is a bad random number generator, but good enough for most randomized    //
//  geometric algorithms.                                                    //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

unsigned long Randomnation(unsigned int choices)
{
  pDoc->randomseed = (pDoc->randomseed * 1366l + 150889l) % 714025l;
  return pDoc->randomseed / (714025l / choices + 1);
}

///////////////////////////////////////////////////////////////////////////////
// Do Ptloc Functions                                                        //                        
// BEGIN                                                                     // 
/////////////////////////////////////////////////////////////////////////////// 

// Find if 'searchtri' is belong to a specific boundary link, 'head' is
//   the first edge of the boundary link.

BOOL SearchBdryEdge(TriEdge* head, TriEdge* searchtri)
{
  TriEdge traverseedge;
	
  // Keep where we start search to startedge.
  traverseedge = *head;
  traverseedge.nextboundaryself();
  while(traverseedge != *head) {
    if (traverseedge == *searchtri) return TRUE;
    traverseedge.nextboundaryself();
  }
  return FALSE; 	
}	

BOOL WalkThroughBoundary(CMeshNode& searchpoint, TriEdge* searchtri)
{
  TriEdge startedge;
  CMeshNode startp;
  CMeshNode *torg, *tdest;
  int orgorient, destorient;

  // Keep where we start search to startedge.
  startedge = *searchtri;
  // First we need form a ray that shot from searchtri to searchpoint.
  //   set startp as this ray's origin. Now, for simplicity, I just
  //   choose startedge's midpoint as startp. 
  torg = startedge.org();
  tdest = startedge.dest();
  startp.x = 0.5 * (torg->x + tdest->x);
  startp.y = 0.5 * (torg->y + tdest->y);
  
  // The WalkThroughBoundary Algorithm:
  //   Take startedge's next edge in boundary link as new searchtri. 
  //   Loop Until we encounter the startedge.
  //     Do current edge and ray intersection check.
  //     If they intersected then
  //       return TRUE. 
  //     else
  //       Take current edge's next edge in boundary link as new searchtri. 
  //       goto Loop.
  //   return FALSE.
  searchtri->nextboundaryself(); 
  while(*searchtri != startedge) {
    torg = searchtri->org();
    tdest = searchtri->dest();
	//startp.x = 0.5 * (torg->x + tdest->x);
    //startp.y = 0.5 * (torg->y + tdest->y);
    orgorient = Orient2D(startp, searchpoint, *torg);
    destorient = Orient2D(startp, searchpoint, *tdest);
    if(orgorient * destorient == 0) {
      if((orgorient == 0) && (destorient != 0)) {
        // startp, searchpoint, *torg collinear,
        if((FuzzyComp(torg->x, searchpoint.x) == 0) &&
           (FuzzyComp(torg->y, searchpoint.y) == 0)) {
          // This mean torg is just the searchpoint, we can return TRUE now.
          //   LocateResult should be eOnVertex.
          return TRUE; 	
        } else {   	
          // Check if `torg' is between `startp' and `searchpoint'.     		
          if (((startp.x < torg->x) == (torg->x < searchpoint.x)) && 
              ((startp.y < torg->y) == (torg->y < searchpoint.y))) {
            return TRUE;
          }
        }  
      } else if((orgorient != 0) && (destorient == 0)) {
        // startp, searchpoint, *tdest collinear,
        if((FuzzyComp(tdest->x, searchpoint.x) == 0) &&
           (FuzzyComp(tdest->y, searchpoint.y) == 0)) {
          // This mean tdest is just the searchpoint, we can return TRUE now.
          //   LocateResult should be eOnVertex.
          return TRUE; 	
        } else {   	
          //   Check if `tdest' is between `startp' and `searchpoint'.     		
          if (((startp.x < tdest->x) == (tdest->x < searchpoint.x)) && 
              ((startp.y < tdest->y) == (tdest->y < searchpoint.y))) {
            return TRUE;
          }
        }
      } else {
        // ray parallel with this edge, I don't care this edge.
      }	
    } else if(orgorient * destorient < 0) {
      orgorient = Orient2D(*torg, *tdest, startp);
      destorient = Orient2D(*torg, *tdest, searchpoint);
      if(orgorient * destorient == 0) {
        if(destorient == 0) {
          // This condition is searchpoint just locate on edge torg->tdest.
          //   LocateResult should be eOnEdge.
          return TRUE;
        } else if(orgorient == 0) {
        	// This condition should not happen. Because only the startedge
        	//   will hold this condition. Loop is already break before get 
        	//   here. 
          throw;
        }
      } else if(orgorient * destorient < 0) {
        return TRUE;
      }
    }	
    searchtri->nextboundaryself(); 
  }
  // We can't walk through this edge. return anyway.
  return FALSE;  
}

enum LocateResult PreciseLocate(CMeshNode& searchpoint, TriEdge* searchtri)
{
  TriEdge backtracetri;
  CMeshNode *forg, *fdest, *fapex;
  int orgorient, destorient;
  int moveleft;
  
  if ((searchtri->ele<0) || (searchtri->ele==INVALIDELEINDEX))
  	  return eOutSide; 
  
  // Where are we?
  forg = searchtri->org(); 
  fdest = searchtri->dest(); 
  fapex = searchtri->apex(); 

  while (1) {
    // Check whether the apex is the point we seek.
    if ((FuzzyComp(fapex->x, searchpoint.x) == 0) && 
        (FuzzyComp(fapex->y, searchpoint.y) == 0)) {
      searchtri->lprevself(); 
      pDoc->recenttri = *searchtri;
      return eOnVertex;
    }
    // Does the point lie on the other side of the line defined by the
    //   triangle edge opposite the triangle's destination?
    destorient = Orient2D(*forg, *fapex, searchpoint);
    // Does the point lie on the other side of the line defined by the
    //   triangle edge opposite the triangle's origin?
    orgorient = Orient2D(*fapex, *fdest, searchpoint);
    if (destorient > 0) {
      if (orgorient > 0) {
        // Move left if the inner product of (fapex - searchpoint) and
        //   (fdest - forg) is positive.  This is equivalent to drawing
        //   a line perpendicular to the line (forg, fdest) passing
        //   through `fapex', and determining which side of this line
        //   `searchpoint' falls on.
        moveleft = (fapex->x - searchpoint.x) * (fdest->x - forg->x) +
                   (fapex->y - searchpoint.y) * (fdest->y - forg->y) > 0.0;
      } else {
        moveleft = 1;
      }
    } else {
      if (orgorient > 0) {
        moveleft = 0;
      } else {
        // The point we seek must be on the boundary of or inside this
        //   triangle.
        if (destorient == 0) {
          searchtri->lprevself(); 
          pDoc->recenttri = *searchtri;
          return eOnEdge;
        }
        if (orgorient == 0) {
          searchtri->lnextself(); 
          pDoc->recenttri = *searchtri;
          return eOnEdge; 
        }
        pDoc->recenttri = *searchtri;
        return eInTriangle;
      }
    }

    // Move to another triangle.  Leave a trace `backtracktri' in case
    //   floating-point roundoff or some such bogey causes us to walk
    //   off a boundary of the triangulation.  We can just bounce off
    //   the boundary as if it were an elastic band.
    if (moveleft) {
      backtracetri = searchtri->lprev(); 
      fdest = fapex;
    } else {
      backtracetri = searchtri->lnext();
      forg = fapex;
    }
    *searchtri = backtracetri.sym(); 
    // Note: Here searchtri may be invalid. Let's check it.
    
    // Check for walking off the edge.
    if (!backtracetri.symexist()) {
      // Turn around.
      *searchtri = backtracetri; 
      // try to walk out the boundary to reach the other side boundary
      if(WalkThroughBoundary(searchpoint, searchtri)) {
        fapex = searchtri->apex(); 
      } else {
        pDoc->recenttri = *searchtri;	
        return eOutSide; 
      }
    } else {
      fapex = searchtri->apex(); 
    }
  }
}

// Used for the point location scheme of Mucke, Saias, and Zhu, to decide
//   how large a random sample of triangles to inspect.
#define SAMPLEFACTOR 11

enum LocateResult Locate(CMeshNode& searchpoint, TriEdge* searchtri, void *pd)
{
  CMeshNode *torg, *tdest;
  double searchdist, dist;
  TriEdge sampletri;
  enum LocateResult locres;
  long samplenum;
  int ahead;
  int i;
  
  pDoc=(CFemmviewDoc *)pd;

  // Record the distance from the suggested starting triangle to the
  //   point we seek.
  if((searchtri->ele == INVALIDELEINDEX) || (searchtri->ele < 0) || 
     (searchtri->ele >= pDoc->pmeshelem->GetSize()) ) {
    searchtri->ele = 0;
  } 
  torg = searchtri->org(); 
  searchdist = (searchpoint.x - torg->x) * (searchpoint.x - torg->x)
             + (searchpoint.y - torg->y) * (searchpoint.y - torg->y);
  // If a recently encountered triangle has been recorded and has not been
  //   deallocated, test it as a good starting point.
  if (pDoc->recenttri.ele != INVALIDELEINDEX) {
    torg = pDoc->recenttri.org(); 
    if ((FuzzyComp(torg->x, searchpoint.x) == 0) && 
        (FuzzyComp(torg->y, searchpoint.y) == 0)) {
      *searchtri = pDoc->recenttri; 
      return eOnVertex;
    }
    dist = (searchpoint.x - torg->x) * (searchpoint.x - torg->x)
         + (searchpoint.y - torg->y) * (searchpoint.y - torg->y);
    if (dist < searchdist) {
      *searchtri = pDoc->recenttri; 
      searchdist = dist;
    }
  } 
  
  // The number of random pDoc->samples taken is proportional to the cube root of
  //   the number of triangles in the mesh.  The next bit of code assumes
  //   that the number of triangles increases monotonically.
  while (SAMPLEFACTOR * pDoc->samples * pDoc->samples * pDoc->samples < pDoc->pmeshelem->GetSize()) {
    pDoc->samples++;
  }
  for(i = 0; i < pDoc->samples; i ++) {
    samplenum = Randomnation(pDoc->pmeshelem->GetSize());
    sampletri.ele = samplenum; 
    torg = sampletri.org(); 
    dist = (searchpoint.x - torg->x) * (searchpoint.x - torg->x)
         + (searchpoint.y - torg->y) * (searchpoint.y - torg->y);
    if (dist < searchdist) {
      *searchtri = sampletri; 
      searchdist = dist;
    }
  }
  
  // Here we need check if the searchtri is a boundary edge, If so, we need 
  //   get a non boundary edge as the searchtri, otherwise the sym() primitive 
  //   will cause exception. 
  if(!searchtri->symexist()) {
    TriEdge nextedge;
    nextedge = searchtri->lnext();
    while(!nextedge.symexist()) {
      nextedge.lnextself();
      if( nextedge == *searchtri) {
      	// It seems there only exist one elem.
      	return eOutSide;
      }	
    }
    *searchtri = nextedge;
  }
  
  // Where are we?
  torg = searchtri->org(); 
  tdest = searchtri->dest(); 
  // Check the starting triangle's vertices.
  if ((FuzzyComp(torg->x, searchpoint.x) == 0) && 
      (FuzzyComp(torg->y, searchpoint.y) == 0)) {
    pDoc->recenttri = *searchtri;
    return eOnVertex;
  }
  if ((FuzzyComp(tdest->x, searchpoint.x) == 0) && 
      (FuzzyComp(tdest->y, searchpoint.y) == 0)) {
    searchtri->lnextself();
    pDoc->recenttri = *searchtri;
    return eOnVertex;
  }
  // Orient `searchtri' to fit the preconditions of calling preciselocate().
  ahead = Orient2D(*torg, *tdest, searchpoint);
  if (ahead < 0) {
    // Turn around so that `searchpoint' is to the left of the
    //   edge specified by `searchtri'.
    searchtri->symself(); 
  } else if (ahead == 0) {
    // Check if `searchpoint' is between `torg' and `tdest'.
    if (((torg->x < searchpoint.x) == (searchpoint.x < tdest->x))
        && ((torg->y < searchpoint.y) == (searchpoint.y < tdest->y))) {
      pDoc->recenttri = *searchtri;
      return eOnEdge;
    }
  }
  locres = PreciseLocate(searchpoint, searchtri);
  if (locres == eOutSide) {
    // is there exist multi-isolated-block in the model?
	if (pDoc->numberofbdrylink > 1) {
		
      for (i = 0; i < pDoc->numberofbdrylink; i++) {
        if (!SearchBdryEdge(&pDoc->bdrylinkhead[i], searchtri)) {
          // Walk to another isolated block in model. 
          *searchtri = pDoc->bdrylinkhead[i];
          torg = searchtri->org();
          tdest = searchtri->dest();
          // Orient `searchtri' to fit the preconditions of calling 
          //   preciselocate().
          ahead = Orient2D(*torg, *tdest, searchpoint);
          if (ahead < 0) {
            // Turn around so that `searchpoint' is to the left of the
            //   edge specified by `searchtri'.
            searchtri->symself(); 
          }
          locres = PreciseLocate(searchpoint, searchtri);
          if (locres != eOutSide) return locres;
        }  	
      }
    }
  }	
  return locres; 
}

/////////////////////////////////////////////////////////////////////////////// 
// END                                                                       //
// Do Ptloc Functions                                                        //                        
///////////////////////////////////////////////////////////////////////////////