G4PolyPhiFace Class Reference

#include <G4PolyPhiFace.hh>

Inheritance diagram for G4PolyPhiFace:

Public Member Functions
	G4PolyPhiFace (const G4ReduciblePolygon *rz, G4double phi, G4double deltaPhi, G4double phiOther)
virtual	~G4PolyPhiFace ()
	G4PolyPhiFace (const G4PolyPhiFace &source)
G4PolyPhiFace &	operator= (const G4PolyPhiFace &source)
G4bool	Intersect (const G4ThreeVector &p, const G4ThreeVector &v, G4bool outgoing, G4double surfTolerance, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal, G4bool &allBehind)
G4double	Distance (const G4ThreeVector &p, G4bool outgoing)
EInside	Inside (const G4ThreeVector &p, G4double tolerance, G4double *bestDistance)
G4ThreeVector	Normal (const G4ThreeVector &p, G4double *bestDistance)
G4double	Extent (const G4ThreeVector axis)
void	CalculateExtent (const EAxis axis, const G4VoxelLimits &voxelLimit, const G4AffineTransform &tranform, G4SolidExtentList &extentList)
G4VCSGface *	Clone ()
G4double	SurfaceArea ()
G4double	SurfaceTriangle (G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4ThreeVector *p4)
G4ThreeVector	GetPointOnFace ()
	G4PolyPhiFace (__void__ &)
void	Diagnose (G4VSolid *solid)
Public Member Functions inherited from G4VCSGface
	G4VCSGface ()
virtual	~G4VCSGface ()
virtual G4bool	Intersect (const G4ThreeVector &p, const G4ThreeVector &v, G4bool outgoing, G4double surfTolerance, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal, G4bool &allBehind)=0
virtual G4double	Distance (const G4ThreeVector &p, G4bool outgoing)=0
virtual EInside	Inside (const G4ThreeVector &p, G4double tolerance, G4double *bestDistance)=0
virtual G4ThreeVector	Normal (const G4ThreeVector &p, G4double *bestDistance)=0
virtual G4double	Extent (const G4ThreeVector axis)=0
virtual void	CalculateExtent (const EAxis axis, const G4VoxelLimits &voxelLimit, const G4AffineTransform &tranform, G4SolidExtentList &extentList)=0
virtual G4VCSGface *	Clone ()=0
virtual G4double	SurfaceArea ()=0
virtual G4ThreeVector	GetPointOnFace ()=0

Protected Member Functions
G4bool	InsideEdgesExact (G4double r, G4double z, G4double normSign, const G4ThreeVector &p, const G4ThreeVector &v)
G4bool	InsideEdges (G4double r, G4double z)
G4bool	InsideEdges (G4double r, G4double z, G4double *distRZ2, G4PolyPhiFaceVertex **base3Dnorm=nullptr, G4ThreeVector **head3Dnorm=nullptr)
G4double	ExactZOrder (G4double z, G4double qx, G4double qy, G4double qz, const G4ThreeVector &v, G4double normSign, const G4PolyPhiFaceVertex *vert) const
void	CopyStuff (const G4PolyPhiFace &source)
G4double	Area2 (G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4bool	Left (G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4bool	LeftOn (G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4bool	Collinear (G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4bool	IntersectProp (G4TwoVector a, G4TwoVector b, G4TwoVector c, G4TwoVector d)
G4bool	Between (G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4bool	Intersect (G4TwoVector a, G4TwoVector b, G4TwoVector c, G4TwoVector d)
G4bool	Diagonalie (G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b)
G4bool	InCone (G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b)
G4bool	Diagonal (G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b)
void	EarInit ()
void	Triangulate ()

Protected Attributes
G4int	numEdges
G4PolyPhiFaceEdge *	edges = nullptr
G4PolyPhiFaceVertex *	corners = nullptr
G4ThreeVector	normal
G4ThreeVector	radial
G4ThreeVector	surface
G4ThreeVector	surface_point
G4double	rMin
G4double	rMax
G4double	zMin
G4double	zMax
G4bool	allBehind = false
G4double	kCarTolerance
G4double	fSurfaceArea = 0.0
G4PolyPhiFaceVertex *	triangles = nullptr

Detailed Description

Definition at line 75 of file G4PolyPhiFace.hh.

Constructor & Destructor Documentation

G4PolyPhiFace() [1/3]

G4PolyPhiFace::G4PolyPhiFace	(	const G4ReduciblePolygon *	rz,
		G4double	phi,
		G4double	deltaPhi,
		G4double	phiOther
	)

Definition at line 51 of file G4PolyPhiFace.cc.

{
  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 
  numEdges = rz->NumVertices();
  
  rMin = rz->Amin();
  rMax = rz->Amax();
  zMin = rz->Bmin();
  zMax = rz->Bmax();
 
  //
  // Is this the "starting" phi edge of the two?
  //
  G4bool start = (phiOther > phi);
  
  //
  // Build radial vector
  //
  radial = G4ThreeVector( std::cos(phi), std::sin(phi), 0.0 );
  
  //
  // Build normal
  //
  G4double zSign = start ? 1 : -1;
  normal = G4ThreeVector( zSign*radial.y(), -zSign*radial.x(), 0 );
  
  //
  // Is allBehind?
  //
  allBehind = (zSign*(std::cos(phiOther)*radial.y() - std::sin(phiOther)*radial.x()) < 0);
  
  //
  // Adjacent edges
  //
  G4double midPhi = phi + (start ? +0.5 : -0.5)*deltaPhi;
  G4double cosMid = std::cos(midPhi), 
                 sinMid = std::sin(midPhi);
  //
  // Allocate corners
  //
  corners = new G4PolyPhiFaceVertex[numEdges];
  //
  // Fill them
  //
  G4ReduciblePolygonIterator iterRZ(rz);
  
  G4PolyPhiFaceVertex* corn   = corners;
  G4PolyPhiFaceVertex* helper = corners;
 
  iterRZ.Begin();
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    corn->r = iterRZ.GetA();
    corn->z = iterRZ.GetB();
    corn->x = corn->r*radial.x();
    corn->y = corn->r*radial.y();
 
    // Add pointer on prev corner
    //
    if( corn == corners )
      { corn->prev = corners+numEdges-1;}
    else
      { corn->prev = helper; }
 
    // Add pointer on next corner
    //
    if( corn < corners+numEdges-1 )
      { corn->next = corn+1;}
    else
      { corn->next = corners; }
 
    helper = corn;
  } while( ++corn, iterRZ.Next() );
 
  //
  // Allocate edges
  //
  edges = new G4PolyPhiFaceEdge[numEdges];
 
  //
  // Fill them
  //
  G4double rFact = std::cos(0.5*deltaPhi);
  G4double rFactNormalize = 1.0/std::sqrt(1.0+rFact*rFact);
 
  G4PolyPhiFaceVertex* prev = corners+numEdges-1,
                     * here = corners;
  G4PolyPhiFaceEdge*   edge = edges;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    G4ThreeVector sideNorm;
    
    edge->v0 = prev;
    edge->v1 = here;
 
    G4double dr = here->r - prev->r,
             dz = here->z - prev->z;
                         
    edge->length = std::sqrt( dr*dr + dz*dz );
 
    edge->tr = dr/edge->length;
    edge->tz = dz/edge->length;
    
    if ((here->r < DBL_MIN) && (prev->r < DBL_MIN))
    {
      //
      // Sigh! Always exceptions!
      // This edge runs at r==0, so its adjoing surface is not a
      // PolyconeSide or PolyhedraSide, but the opposite PolyPhiFace.
      //
      G4double zSignOther = start ? -1 : 1;
      sideNorm = G4ThreeVector(  zSignOther*std::sin(phiOther), 
                                -zSignOther*std::cos(phiOther), 0 );
    }
    else
    {
      sideNorm = G4ThreeVector( edge->tz*cosMid,
                                edge->tz*sinMid,
                               -edge->tr*rFact );
      sideNorm *= rFactNormalize;
    }
    sideNorm += normal;
    
    edge->norm3D = sideNorm.unit();
  } while( edge++, prev=here, ++here < corners+numEdges );
 
  //
  // Go back and fill in corner "normals"
  //
  G4PolyPhiFaceEdge* prevEdge = edges+numEdges-1;
  edge = edges;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    //
    // Calculate vertex 2D normals (on the phi surface)
    //
    G4double rPart = prevEdge->tr + edge->tr;
    G4double zPart = prevEdge->tz + edge->tz;
    G4double norm = std::sqrt( rPart*rPart + zPart*zPart );
    G4double rNorm = +zPart/norm;
    G4double zNorm = -rPart/norm;
    
    edge->v0->rNorm = rNorm;
    edge->v0->zNorm = zNorm;
    
    //
    // Calculate the 3D normals.
    //
    // Find the vector perpendicular to the z axis
    // that defines the plane that contains the vertex normal
    //
    G4ThreeVector xyVector;
    
    if (edge->v0->r < DBL_MIN)
    {
      //
      // This is a vertex at r==0, which is a special
      // case. The normal we will construct lays in the
      // plane at the center of the phi opening.
      //
      // We also know that rNorm < 0
      //
      G4double zSignOther = start ? -1 : 1;
      G4ThreeVector normalOther(  zSignOther*std::sin(phiOther), 
                                 -zSignOther*std::cos(phiOther), 0 );
                
      xyVector = - normal - normalOther;
    }
    else
    {
      //
      // This is a vertex at r > 0. The plane
      // is the average of the normal and the
      // normal of the adjacent phi face
      //
      xyVector = G4ThreeVector( cosMid, sinMid, 0 );
      if (rNorm < 0)
        xyVector -= normal;
      else
        xyVector += normal;
    }
    
    //
    // Combine it with the r/z direction from the face
    //
    edge->v0->norm3D = rNorm*xyVector.unit() + G4ThreeVector( 0, 0, zNorm );
  } while(  prevEdge=edge, ++edge < edges+numEdges );
  
  //
  // Build point on surface
  //
  G4double rAve = 0.5*(rMax-rMin),
           zAve = 0.5*(zMax-zMin);
  surface = G4ThreeVector( rAve*radial.x(), rAve*radial.y(), zAve );
}

~G4PolyPhiFace()

G4PolyPhiFace::~G4PolyPhiFace ( )

virtual

Definition at line 284 of file G4PolyPhiFace.cc.

{
  delete [] edges;
  delete [] corners;
}

G4PolyPhiFace() [2/3]

G4PolyPhiFace::G4PolyPhiFace

(

const G4PolyPhiFace &

source

)

Definition at line 294 of file G4PolyPhiFace.cc.

  : G4VCSGface()
{
  CopyStuff( source );
}

G4PolyPhiFace() [3/3]

G4PolyPhiFace::G4PolyPhiFace

(

__void__ &

)

Definition at line 275 of file G4PolyPhiFace.cc.

  : numEdges(0), rMin(0.), rMax(0.), zMin(0.), zMax(0.), kCarTolerance(0.)
{
}

Member Function Documentation

Area2()

G4double G4PolyPhiFace::Area2 ( G4TwoVector  a, G4TwoVector  b, G4TwoVector  c  )

protected

Definition at line 914 of file G4PolyPhiFace.cc.

{
  return ((b.x()-a.x())*(c.y()-a.y())-
          (c.x()-a.x())*(b.y()-a.y()));
}

Referenced by Collinear(), Left(), and LeftOn().

Between()

G4bool G4PolyPhiFace::Between ( G4TwoVector  a, G4TwoVector  b, G4TwoVector  c  )

protected

Definition at line 967 of file G4PolyPhiFace.cc.

{
  if( !Collinear(a,b,c) ) { return false; }
 
  if(a.x()!=b.x())
  {
    return ((a.x()<=c.x())&&(c.x()<=b.x()))||
           ((a.x()>=c.x())&&(c.x()>=b.x()));
  }
  else
  {
    return ((a.y()<=c.y())&&(c.y()<=b.y()))||
           ((a.y()>=c.y())&&(c.y()>=b.y()));
  }
}

Referenced by Intersect().

CalculateExtent()

void G4PolyPhiFace::CalculateExtent ( const EAxis  axis, const G4VoxelLimits &  voxelLimit, const G4AffineTransform &  tranform, G4SolidExtentList &  extentList  )

virtual

Implements G4VCSGface.

Definition at line 599 of file G4PolyPhiFace.cc.

{
  //
  // Construct a (sometimes big) clippable polygon, 
  //
  // Perform the necessary transformations while doing so
  //
  G4ClippablePolygon polygon;
  
  G4PolyPhiFaceVertex* corner = corners;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    G4ThreeVector point( 0, 0, corner->z );
    point += radial*corner->r;
    
    polygon.AddVertexInOrder( transform.TransformPoint( point ) );
  } while( ++corner < corners + numEdges );
  
  //
  // Clip away
  //
  if (polygon.PartialClip( voxelLimit, axis ))
  {
    //
    // Add it to the list
    //
    polygon.SetNormal( transform.TransformAxis(normal) );
    extentList.AddSurface( polygon );
  }
}

Clone()

G4VCSGface * G4PolyPhiFace::Clone ( )

inlinevirtual

Implements G4VCSGface.

Collinear()

G4bool G4PolyPhiFace::Collinear ( G4TwoVector  a, G4TwoVector  b, G4TwoVector  c  )

protected

Definition at line 942 of file G4PolyPhiFace.cc.

{
  return Area2(a,b,c)==0;
}

Referenced by Between(), and IntersectProp().

CopyStuff()

void G4PolyPhiFace::CopyStuff ( const G4PolyPhiFace &  source )

protected

Definition at line 320 of file G4PolyPhiFace.cc.

{
  //
  // The simple stuff
  //
  numEdges  = source.numEdges;
  normal    = source.normal;
  radial    = source.radial;
  surface   = source.surface;
  rMin    = source.rMin;
  rMax    = source.rMax;
  zMin    = source.zMin;
  zMax    = source.zMax;
  allBehind  = source.allBehind;
  triangles  = nullptr;
 
  kCarTolerance = source.kCarTolerance;
  fSurfaceArea = source.fSurfaceArea;
 
  //
  // Corner dynamic array
  //
  corners = new G4PolyPhiFaceVertex[numEdges];
  G4PolyPhiFaceVertex *corn = corners,
                      *sourceCorn = source.corners;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    *corn = *sourceCorn;
  } while( ++sourceCorn, ++corn < corners+numEdges );
  
  //
  // Edge dynamic array
  //
  edges = new G4PolyPhiFaceEdge[numEdges];
 
  G4PolyPhiFaceVertex* prev = corners+numEdges-1,
                     * here = corners;
  G4PolyPhiFaceEdge*   edge = edges,
                   *   sourceEdge = source.edges;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    *edge = *sourceEdge;
    edge->v0 = prev;
    edge->v1 = here;
  } while( ++sourceEdge, ++edge, prev=here, ++here < corners+numEdges );
}

Referenced by G4PolyPhiFace(), and operator=().

Diagnose()

void G4PolyPhiFace::Diagnose

(

G4VSolid *

solid

)

Definition at line 258 of file G4PolyPhiFace.cc.

{
  G4PolyPhiFaceVertex* corner = corners;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    G4ThreeVector test(corner->x, corner->y, corner->z);
    test -= 1E-6*corner->norm3D;
    
    if (owner->Inside(test) != kInside) 
      G4Exception( "G4PolyPhiFace::Diagnose()", "GeomSolids0002",
                   FatalException, "Bad vertex normal found." );
  } while( ++corner < corners+numEdges );
}

Diagonal()

G4bool G4PolyPhiFace::Diagonal ( G4PolyPhiFaceVertex *  a, G4PolyPhiFaceVertex *  b  )

protected

Definition at line 1063 of file G4PolyPhiFace.cc.

{ 
  return InCone(a,b) && InCone(b,a) && Diagonalie(a,b);
}

Referenced by EarInit(), and Triangulate().

Diagonalie()

G4bool G4PolyPhiFace::Diagonalie ( G4PolyPhiFaceVertex *  a, G4PolyPhiFaceVertex *  b  )

protected

Definition at line 1004 of file G4PolyPhiFace.cc.

{
  G4PolyPhiFaceVertex   *corner = triangles;
  G4PolyPhiFaceVertex   *corner_next=triangles;
 
  // For each Edge (corner,corner_next) 
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    corner_next=corner->next;
 
    // Skip edges incident to a of b
    //
    if( (corner!=a)&&(corner_next!=a)
      &&(corner!=b)&&(corner_next!=b) )
    {
       G4TwoVector rz1,rz2,rz3,rz4;
       rz1 = G4TwoVector(a->r,a->z);
       rz2 = G4TwoVector(b->r,b->z);
       rz3 = G4TwoVector(corner->r,corner->z);
       rz4 = G4TwoVector(corner_next->r,corner_next->z);
       if( Intersect(rz1,rz2,rz3,rz4) ) { return false; }
    }
    corner=corner->next;   
   
  } while( corner != triangles );
 
  return true;
}

Referenced by Diagonal().

Distance()

G4double G4PolyPhiFace::Distance ( const G4ThreeVector &  p, G4bool  outgoing  )

virtual

Implements G4VCSGface.

Definition at line 424 of file G4PolyPhiFace.cc.

{
  G4double normSign = outgoing ? +1 : -1;
  //
  // Correct normal? 
  //
  G4ThreeVector ps = p - surface;
  G4double distPhi = -normSign*normal.dot(ps);
  
  if (distPhi < -0.5*kCarTolerance) 
    return kInfinity;
  else if (distPhi < 0)
    distPhi = 0.0;
  
  //
  // Calculate projected point in r,z
  //
  G4double r = radial.dot(p);
  
  //
  // Are we inside the face?
  //
  G4double distRZ2;
  
  if (InsideEdges( r, p.z(), &distRZ2, 0 ))
  {
    //
    // Yup, answer is just distPhi
    //
    return distPhi;
  }
  else
  {
    //
    // Nope. Penalize by distance out
    //
    return std::sqrt( distPhi*distPhi + distRZ2 );
  }
}  

EarInit()

void G4PolyPhiFace::EarInit ( )

protected

Definition at line 1071 of file G4PolyPhiFace.cc.

{
  G4PolyPhiFaceVertex* corner = triangles;
  G4PolyPhiFaceVertex* c_prev, *c_next;
  
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
     // We need to determine three consecutive vertices
     //
     c_next=corner->next;
     c_prev=corner->prev; 
 
     // Calculation of ears
     //
     corner->ear=Diagonal(c_prev,c_next);   
     corner=corner->next;
 
  } while( corner!=triangles );
}

Referenced by Triangulate().

ExactZOrder()

G4double G4PolyPhiFace::ExactZOrder ( G4double  z, G4double  qx, G4double  qy, G4double  qz, const G4ThreeVector &  v, G4double  normSign, const G4PolyPhiFaceVertex *  vert  ) const

inlineprotected

Referenced by InsideEdgesExact().

Extent()

G4double G4PolyPhiFace::Extent ( const G4ThreeVector  axis )

virtual

Implements G4VCSGface.

Definition at line 579 of file G4PolyPhiFace.cc.

{
  G4double max = -kInfinity;
  
  G4PolyPhiFaceVertex* corner = corners;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    G4double here = axis.x()*corner->r*radial.x()
            + axis.y()*corner->r*radial.y()
            + axis.z()*corner->z;
    if (here > max) max = here;
  } while( ++corner < corners + numEdges );
  
  return max;
}  

GetPointOnFace()

G4ThreeVector G4PolyPhiFace::GetPointOnFace ( )

virtual

Implements G4VCSGface.

Definition at line 902 of file G4PolyPhiFace.cc.

{
  Triangulate();
  return surface_point;
}

InCone()

G4bool G4PolyPhiFace::InCone ( G4PolyPhiFaceVertex *  a, G4PolyPhiFaceVertex *  b  )

protected

Definition at line 1037 of file G4PolyPhiFace.cc.

{
  // a0,a and a1 are consecutive vertices
  //
  G4PolyPhiFaceVertex *a0,*a1;
  a1=a->next;
  a0=a->prev;
 
  G4TwoVector arz,arz0,arz1,brz;
  arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector(a0->r,a0->z);
  arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVector(b->r,b->z);
    
  
  if(LeftOn(arz,arz1,arz0))  // If a is convex vertex
  {
    return Left(arz,brz,arz0)&&Left(brz,arz,arz1);
  }
  else                       // Else a is reflex
  {
    return !( LeftOn(arz,brz,arz1)&&LeftOn(brz,arz,arz0));
  }
}

Referenced by Diagonal().

Inside()

EInside G4PolyPhiFace::Inside ( const G4ThreeVector &  p, G4double  tolerance, G4double *  bestDistance  )

virtual

Implements G4VCSGface.

Definition at line 466 of file G4PolyPhiFace.cc.

{
  //
  // Get distance along phi, which if negative means the point
  // is nominally inside the shape.
  //
  G4ThreeVector ps = p - surface;
  G4double distPhi = normal.dot(ps);
  
  //
  // Calculate projected point in r,z
  //
  G4double r = radial.dot(p);
  
  //
  // Are we inside the face?
  //
  G4double distRZ2;
  G4PolyPhiFaceVertex* base3Dnorm = nullptr;
  G4ThreeVector*       head3Dnorm = nullptr;
  
  if (InsideEdges( r, p.z(), &distRZ2, &base3Dnorm, &head3Dnorm ))
  {
    //
    // Looks like we're inside. Distance is distance in phi.
    //
    *bestDistance = std::fabs(distPhi);
    
    //
    // Use distPhi to decide fate
    //
    if (distPhi < -tolerance) return kInside;
    if (distPhi <  tolerance) return kSurface;
    return kOutside;
  }
  else
  {
    //
    // We're outside the extent of the face,
    // so the distance is penalized by distance from edges in RZ
    //
    *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
    
    //
    // Use edge normal to decide fate
    //
    G4ThreeVector cc( base3Dnorm->r*radial.x(),
                      base3Dnorm->r*radial.y(),
                      base3Dnorm->z );
    cc = p - cc;
    G4double normDist = head3Dnorm->dot(cc);
    if ( distRZ2 > tolerance*tolerance )
    {
      //
      // We're far enough away that kSurface is not possible
      //
      return normDist < 0 ? kInside : kOutside;
    }
    
    if (normDist < -tolerance) return kInside;
    if (normDist <  tolerance) return kSurface;
    return kOutside;
  }
}  

InsideEdges() [1/2]

G4bool G4PolyPhiFace::InsideEdges ( G4double  r, G4double  z  )

protected

Definition at line 778 of file G4PolyPhiFace.cc.

{
  //
  // Quick check of extent
  //
  if ( r < rMin || r > rMax ) return false;
  if ( z < zMin || z > zMax ) return false;
  
  //
  // More thorough check
  //
  G4double notUsed;
  
  return InsideEdges( r, z, &notUsed, 0 );
}

Referenced by Distance(), Inside(), InsideEdges(), and Normal().

InsideEdges() [2/2]

G4bool G4PolyPhiFace::InsideEdges ( G4double  r, G4double  z, G4double *  distRZ2, G4PolyPhiFaceVertex **  base3Dnorm = nullptr, G4ThreeVector **  head3Dnorm = nullptr  )

protected

Definition at line 798 of file G4PolyPhiFace.cc.

{
  G4double bestDistance2 = kInfinity;
  G4bool answer = false;
  
  G4PolyPhiFaceEdge* edge = edges;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    G4PolyPhiFaceVertex* testMe = nullptr;
    G4PolyPhiFaceVertex* v0_vertex = edge->v0;
    G4PolyPhiFaceVertex* v1_vertex = edge->v1;
    //
    // Get distance perpendicular to the edge
    //
    G4double dr = (r-v0_vertex->r), dz = (z-v0_vertex->z);
 
    G4double distOut = dr*edge->tz - dz*edge->tr;
    G4double distance2 = distOut*distOut;
    if (distance2 > bestDistance2) continue;        // No hope!
 
    //
    // Check to see if normal intersects edge within the edge's boundary
    //
    G4double q = dr*edge->tr + dz*edge->tz;
 
    //
    // If it doesn't, penalize distance2 appropriately
    //
    if (q < 0)
    {
      distance2 += q*q;
      testMe = v0_vertex;
    }
    else if (q > edge->length)
    {
      G4double s2 = q-edge->length;
      distance2 += s2*s2;
      testMe = v1_vertex;
    }
    
    //
    // Closest edge so far?
    //
    if (distance2 < bestDistance2)
    {
      bestDistance2 = distance2;
      if (testMe != nullptr)
      {
        G4double distNorm = dr*testMe->rNorm + dz*testMe->zNorm;
        answer = (distNorm <= 0);
        if (base3Dnorm != nullptr)
        {
          *base3Dnorm = testMe;
          *head3Dnorm = &testMe->norm3D;
        }
      }
      else
      {
        answer = (distOut <= 0);                        
        if (base3Dnorm != nullptr)
        {
          *base3Dnorm = v0_vertex;
          *head3Dnorm = &edge->norm3D;
        }
      }
    }
  } while( ++edge < edges + numEdges );
  
  *bestDist2 = bestDistance2;
  return answer;
}

InsideEdgesExact()

G4bool G4PolyPhiFace::InsideEdgesExact ( G4double  r, G4double  z, G4double  normSign, const G4ThreeVector &  p, const G4ThreeVector &  v  )

protected

Definition at line 648 of file G4PolyPhiFace.cc.

{
  //
  // Quick check of extent
  //
  if ( (r < rMin-kCarTolerance)
    || (r > rMax+kCarTolerance) ) return false;
       
  if ( (z < zMin-kCarTolerance)
    || (z > zMax+kCarTolerance) ) return false;
  
  //
  // Exact check: loop over all vertices
  //
  G4double qx = p.x() + v.x(),
           qy = p.y() + v.y(),
           qz = p.z() + v.z();
 
  G4int answer = 0;
  G4PolyPhiFaceVertex *corn = corners, 
                      *prev = corners+numEdges-1;
 
  G4double cornZ, prevZ;
  
  prevZ = ExactZOrder( z, qx, qy, qz, v, normSign, prev );
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    //
    // Get z order of this vertex, and compare to previous vertex
    //
    cornZ = ExactZOrder( z, qx, qy, qz, v, normSign, corn );
    
    if (cornZ < 0)
    {
      if (prevZ < 0) continue;
    }
    else if (cornZ > 0)
    {
      if (prevZ > 0) continue;
    }
    else
    {
      //
      // By chance, we overlap exactly (within precision) with 
      // the current vertex. Continue if the same happened previously
      // (e.g. the previous vertex had the same z value)
      //
      if (prevZ == 0) continue;
      
      //
      // Otherwise, to decide what to do, we need to know what is
      // coming up next. Specifically, we need to find the next vertex
      // with a non-zero z order.
      //
      // One might worry about infinite loops, but the above conditional
      // should prevent it
      //
      G4PolyPhiFaceVertex *next = corn;
      G4double nextZ;
      do    // Loop checking, 13.08.2015, G.Cosmo
      {
        ++next;
        if (next == corners+numEdges) next = corners;
 
        nextZ = ExactZOrder( z, qx, qy, qz, v, normSign, next );
      } while( nextZ == 0 );
      
      //
      // If we won't be changing direction, go to the next vertex
      //
      if (nextZ*prevZ < 0) continue;
    }
  
      
    //
    // We overlap in z with the side of the face that stretches from
    // vertex "prev" to "corn". On which side (left or right) do
    // we lay with respect to this segment?
    //  
    G4ThreeVector qa( qx - prev->x, qy - prev->y, qz - prev->z ),
                  qb( qx - corn->x, qy - corn->y, qz - corn->z );
 
    G4double aboveOrBelow = normSign*qa.cross(qb).dot(v);
    
    if (aboveOrBelow > 0) 
      ++answer;
    else if (aboveOrBelow < 0)
      --answer;
    else
    {
      //
      // A precisely zero answer here means we exactly
      // intersect (within roundoff) the edge of the face.
      // Return true in this case.
      //
      return true;
    }
  } while( prevZ = cornZ, prev=corn, ++corn < corners+numEdges );
  
  return answer!=0;
}

Referenced by Intersect().

Intersect() [1/2]

G4bool G4PolyPhiFace::Intersect ( const G4ThreeVector &  p, const G4ThreeVector &  v, G4bool  outgoing, G4double  surfTolerance, G4double &  distance, G4double &  distFromSurface, G4ThreeVector &  normal, G4bool &  allBehind  )

virtual

Implements G4VCSGface.

Definition at line 369 of file G4PolyPhiFace.cc.

{
  G4double normSign = outgoing ? +1 : -1;
  
  //
  // These don't change
  //
  isAllBehind = allBehind;
  aNormal = normal;
 
  //
  // Correct normal? Here we have straight sides, and can safely ignore
  // intersections where the dot product with the normal is zero.
  //
  G4double dotProd = normSign*normal.dot(v);
  
  if (dotProd <= 0) return false;
 
  //
  // Calculate distance to surface. If the side is too far
  // behind the point, we must reject it.
  //
  G4ThreeVector ps = p - surface;
  distFromSurface = -normSign*ps.dot(normal);
    
  if (distFromSurface < -surfTolerance) return false;
 
  //
  // Calculate precise distance to intersection with the side
  // (along the trajectory, not normal to the surface)
  //
  distance = distFromSurface/dotProd;
 
  //
  // Calculate intersection point in r,z
  //
  G4ThreeVector ip = p + distance*v;
  
  G4double r = radial.dot(ip);
  
  //
  // And is it inside the r/z extent?
  //
  return InsideEdgesExact( r, ip.z(), normSign, p, v );
}

Referenced by Diagonalie().

Intersect() [2/2]

G4bool G4PolyPhiFace::Intersect ( G4TwoVector  a, G4TwoVector  b, G4TwoVector  c, G4TwoVector  d  )

protected

Definition at line 986 of file G4PolyPhiFace.cc.

{
 if( IntersectProp(a,b,c,d) )
   { return true; }
 else if( Between(a,b,c)||
          Between(a,b,d)||
          Between(c,d,a)||
          Between(c,d,b) )
   { return true; }
 else
   { return false; }
}

IntersectProp()

G4bool G4PolyPhiFace::IntersectProp ( G4TwoVector  a, G4TwoVector  b, G4TwoVector  c, G4TwoVector  d  )

protected

Definition at line 952 of file G4PolyPhiFace.cc.

{
  if( Collinear(a,b,c) || Collinear(a,b,d)||
      Collinear(c,d,a) || Collinear(c,d,b) )  { return false; }
 
  G4bool Positive;
  Positive = !(Left(a,b,c))^!(Left(a,b,d));
  return Positive && (!Left(c,d,a)^!Left(c,d,b));
}

Referenced by Intersect().

Left()

G4bool G4PolyPhiFace::Left ( G4TwoVector  a, G4TwoVector  b, G4TwoVector  c  )

protected

Definition at line 924 of file G4PolyPhiFace.cc.

{
  return Area2(a,b,c)>0;
}

Referenced by InCone(), and IntersectProp().

LeftOn()

G4bool G4PolyPhiFace::LeftOn ( G4TwoVector  a, G4TwoVector  b, G4TwoVector  c  )

protected

Definition at line 933 of file G4PolyPhiFace.cc.

{
  return Area2(a,b,c)>=0;
}

Referenced by InCone().

Normal()

G4ThreeVector G4PolyPhiFace::Normal ( const G4ThreeVector &  p, G4double *  bestDistance  )

virtual

Implements G4VCSGface.

Definition at line 538 of file G4PolyPhiFace.cc.

{
  //
  // Get distance along phi, which if negative means the point
  // is nominally inside the shape.
  //
  G4double distPhi = normal.dot(p);
 
  //
  // Calculate projected point in r,z
  //
  G4double r = radial.dot(p);
  
  //
  // Are we inside the face?
  //
  G4double distRZ2;
  
  if (InsideEdges( r, p.z(), &distRZ2, 0 ))
  {
    //
    // Yup, answer is just distPhi
    //
    *bestDistance = std::fabs(distPhi);
  }
  else
  {
    //
    // Nope. Penalize by distance out
    //
    *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
  }
  
  return normal;
}

operator=()

G4PolyPhiFace & G4PolyPhiFace::operator=

(

const G4PolyPhiFace &

source

)

Definition at line 304 of file G4PolyPhiFace.cc.

{
  if (this == &source)  { return *this; }
 
  delete [] edges;
  delete [] corners;
  
  CopyStuff( source );
  
  return *this;
}

SurfaceArea()

G4double G4PolyPhiFace::SurfaceArea ( )

virtual

Implements G4VCSGface.

Definition at line 894 of file G4PolyPhiFace.cc.

{
  if ( fSurfaceArea==0. ) { Triangulate(); }
  return fSurfaceArea;
}

SurfaceTriangle()

G4double G4PolyPhiFace::SurfaceTriangle	(	G4ThreeVector	p1,
		G4ThreeVector	p2,
		G4ThreeVector	p3,
		G4ThreeVector *	p4
	)

Definition at line 876 of file G4PolyPhiFace.cc.

{
  G4ThreeVector v, w;
  
  v = p3 - p1;
  w = p1 - p2;
  G4double lambda1 = G4UniformRand();
  G4double lambda2 = lambda1*G4UniformRand();
 
  *p4=p2 + lambda1*w + lambda2*v;
  return 0.5*(v.cross(w)).mag();
}

Referenced by Triangulate().

Triangulate()

void G4PolyPhiFace::Triangulate ( )

protected

Definition at line 1094 of file G4PolyPhiFace.cc.

{ 
  // The copy of Polycone is made and this copy is reordered in order to 
  // have a list of triangles. This list is used for GetPointOnFace().
 
  G4PolyPhiFaceVertex* tri_help = new G4PolyPhiFaceVertex[numEdges];
  triangles = tri_help;
  G4PolyPhiFaceVertex* triang = triangles;
 
  std::vector<G4double> areas;
  std::vector<G4ThreeVector> points;
  G4double area=0.;
  G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
  v2=triangles;
 
  // Make copy for prev/next for triang=corners
  //
  G4PolyPhiFaceVertex* helper  = corners;
  G4PolyPhiFaceVertex* helper2 = corners;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    triang->r = helper->r;
    triang->z = helper->z;
    triang->x = helper->x;
    triang->y= helper->y;
 
    // add pointer on prev corner
    //
    if( helper==corners )
      { triang->prev=triangles+numEdges-1; }
    else
      { triang->prev=helper2; }
 
    // add pointer on next corner
    //
    if( helper<corners+numEdges-1 )
      { triang->next=triang+1; }
    else
      { triang->next=triangles; }
    helper2=triang;
    helper=helper->next;
    triang=triang->next;
 
  } while( helper!=corners );
 
  EarInit();
 
  G4int n=numEdges;
  G4int i=0;
  G4ThreeVector p1,p2,p3,p4;
  const G4int max_n_loops=numEdges*10000; // protection against infinite loop
 
  // Each step of outer loop removes one ear
  //
  while(n>3)      // Loop checking, 13.08.2015, G.Cosmo
  {               // Inner loop searches for one ear
    v2=triangles; 
    do    // Loop checking, 13.08.2015, G.Cosmo
    {
      if(v2->ear)  // Ear found. Fill variables
      {
        // (v1,v3) is diagonal
        //
        v3=v2->next; v4=v3->next;
        v1=v2->prev; v0=v1->prev;
        
        // Calculate areas and points
 
        p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
        p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
        p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
 
        G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
        points.push_back(p4);
        areas.push_back(result1);
        area=area+result1;
 
        // Update earity of diagonal endpoints
        //
        v1->ear=Diagonal(v0,v3);
        v3->ear=Diagonal(v1,v4);
 
        // Cut off the ear v2 
        // Has to be done for a copy and not for real PolyPhiFace
        //
        v1->next=v3;
        v3->prev=v1;
        triangles=v3; // In case the head was v2
        --n;
 
        break; // out of inner loop
      }        // end if ear found
 
      v2=v2->next;
    
    } while( v2!=triangles );
 
    ++i;
    if(i>=max_n_loops)
    {
      G4Exception( "G4PolyPhiFace::Triangulation()",
                   "GeomSolids0003", FatalException,
                   "Maximum number of steps is reached for triangulation!" );
    }
  }   // end outer while loop
 
  if(v2->next)
  {
     // add last triangle
     //
     v2=v2->next;
     p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
     p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
     p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
     G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
     points.push_back(p4);
     areas.push_back(result1);
     area=area+result1;
  }
 
  // Surface Area is stored
  //
  fSurfaceArea = area;
  
  // Second Step: choose randomly one surface
  //
  G4double chose = area*G4UniformRand();
   
  // Third Step: Get a point on choosen surface
  //
  G4double Achose1, Achose2;
  Achose1=0; Achose2=0.; 
  i=0;
  do     // Loop checking, 13.08.2015, G.Cosmo
  {
    Achose2+=areas[i];
    if(chose>=Achose1 && chose<Achose2)
    {
      G4ThreeVector point;
       point=points[i] ;
       surface_point=point;
       break;     
    }
    i++; Achose1=Achose2;
  } while( i<numEdges-2 );
 
  delete [] tri_help;
  tri_help = nullptr;
}

Referenced by GetPointOnFace(), and SurfaceArea().

Member Data Documentation

allBehind

G4bool G4PolyPhiFace::allBehind = false

protected

Definition at line 221 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), G4PolyPhiFace(), and Intersect().

corners

G4PolyPhiFaceVertex* G4PolyPhiFace::corners = nullptr

protected

Definition at line 213 of file G4PolyPhiFace.hh.

Referenced by CalculateExtent(), CopyStuff(), Diagnose(), Extent(), G4PolyPhiFace(), InsideEdgesExact(), operator=(), Triangulate(), and ~G4PolyPhiFace().

edges

G4PolyPhiFaceEdge* G4PolyPhiFace::edges = nullptr

protected

Definition at line 212 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), G4PolyPhiFace(), InsideEdges(), operator=(), and ~G4PolyPhiFace().

fSurfaceArea

G4double G4PolyPhiFace::fSurfaceArea = 0.0

protected

Definition at line 224 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), SurfaceArea(), and Triangulate().

kCarTolerance

G4double G4PolyPhiFace::kCarTolerance

protected

Definition at line 223 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), Distance(), G4PolyPhiFace(), and InsideEdgesExact().

normal

G4ThreeVector G4PolyPhiFace::normal

protected

Definition at line 214 of file G4PolyPhiFace.hh.

Referenced by CalculateExtent(), CopyStuff(), Distance(), G4PolyPhiFace(), Inside(), Intersect(), and Normal().

numEdges

G4int G4PolyPhiFace::numEdges

protected

Definition at line 211 of file G4PolyPhiFace.hh.

Referenced by CalculateExtent(), CopyStuff(), Diagnose(), Extent(), G4PolyPhiFace(), InsideEdges(), InsideEdgesExact(), and Triangulate().

radial

G4ThreeVector G4PolyPhiFace::radial

protected

Definition at line 215 of file G4PolyPhiFace.hh.

Referenced by CalculateExtent(), CopyStuff(), Distance(), Extent(), G4PolyPhiFace(), Inside(), Intersect(), and Normal().

rMax

G4double G4PolyPhiFace::rMax

protected

Definition at line 219 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), G4PolyPhiFace(), InsideEdges(), and InsideEdgesExact().

rMin

G4double G4PolyPhiFace::rMin

protected

Definition at line 219 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), G4PolyPhiFace(), and InsideEdgesExact().

surface

G4ThreeVector G4PolyPhiFace::surface

protected

Definition at line 216 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), Distance(), G4PolyPhiFace(), Inside(), and Intersect().

surface_point

G4ThreeVector G4PolyPhiFace::surface_point

protected

Definition at line 217 of file G4PolyPhiFace.hh.

Referenced by GetPointOnFace(), and Triangulate().

triangles

G4PolyPhiFaceVertex* G4PolyPhiFace::triangles = nullptr

protected

Definition at line 225 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), Diagonalie(), EarInit(), and Triangulate().

zMax

G4double G4PolyPhiFace::zMax

protected

Definition at line 220 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), G4PolyPhiFace(), InsideEdges(), and InsideEdgesExact().

zMin

G4double G4PolyPhiFace::zMin

protected

Definition at line 220 of file G4PolyPhiFace.hh.

Referenced by CopyStuff(), G4PolyPhiFace(), and InsideEdgesExact().

---

The documentation for this class was generated from the following files:

• G4PolyPhiFace.hh
• G4PolyPhiFace.cc