G4PolyPhiFace Class Reference

#include <G4PolyPhiFace.hh>

Inheritance diagram for G4PolyPhiFace:

Public Member Functions
	G4PolyPhiFace (const G4ReduciblePolygon *rz, G4double phi, G4double deltaPhi, G4double phiOther)
	~G4PolyPhiFace () override
	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) override
G4double	Distance (const G4ThreeVector &p, G4bool outgoing) override
EInside	Inside (const G4ThreeVector &p, G4double tolerance, G4double *bestDistance) override
G4ThreeVector	Normal (const G4ThreeVector &p, G4double *bestDistance) override
G4double	Extent (const G4ThreeVector axis) override
void	CalculateExtent (const EAxis axis, const G4VoxelLimits &voxelLimit, const G4AffineTransform &tranform, G4SolidExtentList &extentList) override
G4VCSGface *	Clone () override
G4double	SurfaceArea () override
G4double	SurfaceTriangle (const G4ThreeVector &p1, const G4ThreeVector &p2, const G4ThreeVector &p3, G4ThreeVector *p4)
G4ThreeVector	GetPointOnFace () override
	G4PolyPhiFace (__void__ &)
void	Diagnose (G4VSolid *solid)
Public Member Functions inherited from G4VCSGface
	G4VCSGface ()=default
virtual	~G4VCSGface ()=default

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 (const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c)
G4bool	Left (const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c)
G4bool	LeftOn (const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c)
G4bool	Collinear (const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c)
G4bool	IntersectProp (const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c, const G4TwoVector &d)
G4bool	Between (const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c)
G4bool	Intersect (const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c, const 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 = 0
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
  //
  const std::size_t maxEdges = numEdges>0 ? numEdges : 1;
  corners = new G4PolyPhiFaceVertex[maxEdges];
  //
  // 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+maxEdges-1;}
    else
      { corn->prev = helper; }
 
    // Add pointer on next corner
    //
    if( corn < corners+maxEdges-1 )
      { corn->next = corn+1;}
    else
      { corn->next = corners; }
 
    helper = corn;
  } while( ++corn, iterRZ.Next() );
 
  //
  // Allocate edges
  //
  edges = new G4PolyPhiFaceEdge[maxEdges];
 
  //
  // Fill them
  //
  G4double rFact = std::cos(0.5*deltaPhi);
  G4double rFactNormalize = 1.0/std::sqrt(1.0+rFact*rFact);
 
  G4PolyPhiFaceVertex* prev = corners+maxEdges-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+maxEdges );
 
  //
  // Go back and fill in corner "normals"
  //
  G4PolyPhiFaceEdge* prevEdge = edges+maxEdges-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+maxEdges );
  
  //
  // 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 ( )

override

Definition at line 285 of file G4PolyPhiFace.cc.

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

G4PolyPhiFace() [2/3]

G4PolyPhiFace::G4PolyPhiFace

(

const G4PolyPhiFace &

source

)

Definition at line 295 of file G4PolyPhiFace.cc.

{
  CopyStuff( source );
}

G4PolyPhiFace() [3/3]

G4PolyPhiFace::G4PolyPhiFace

(

__void__ &

)

Definition at line 276 of file G4PolyPhiFace.cc.

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

Member Function Documentation

Area2()

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

protected

Definition at line 916 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 ( const G4TwoVector & a, const G4TwoVector & b, const G4TwoVector & c )

protected

Definition at line 969 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 )

overridevirtual

Implements G4VCSGface.

Definition at line 601 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 ( )

inlineoverridevirtual

Implements G4VCSGface.

Collinear()

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

protected

Definition at line 944 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;
 
  const std::size_t maxEdges = (numEdges > 0) ? numEdges : 1;
 
  //
  // Corner dynamic array
  //
  corners = new G4PolyPhiFaceVertex[maxEdges];
  G4PolyPhiFaceVertex *corn = corners,
                      *sourceCorn = source.corners;
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    *corn = *sourceCorn;
  } while( ++sourceCorn, ++corn < corners+maxEdges );
  
  //
  // Edge dynamic array
  //
  edges = new G4PolyPhiFaceEdge[maxEdges];
 
  G4PolyPhiFaceVertex* prev = corners+maxEdges-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+maxEdges );
}

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

Diagnose()

void G4PolyPhiFace::Diagnose

(

G4VSolid *

solid

)

Definition at line 259 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 1065 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 1006 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 )

overridevirtual

Implements G4VCSGface.

Definition at line 426 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, nullptr ))
  {
    //
    // 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 1073 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 )

overridevirtual

Implements G4VCSGface.

Definition at line 581 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 ( )

overridevirtual

Implements G4VCSGface.

Definition at line 904 of file G4PolyPhiFace.cc.

{
  Triangulate();
  return surface_point;
}

InCone()

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

protected

Definition at line 1039 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 )

overridevirtual

Implements G4VCSGface.

Definition at line 468 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 780 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, nullptr );
}

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 800 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 650 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 )

overridevirtual

Implements G4VCSGface.

Definition at line 371 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 ( const G4TwoVector & a, const G4TwoVector & b, const G4TwoVector & c, const G4TwoVector & d )

protected

Definition at line 988 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 ( const G4TwoVector & a, const G4TwoVector & b, const G4TwoVector & c, const G4TwoVector & d )

protected

Definition at line 954 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)) != 0);
}

Referenced by Intersect().

Left()

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

protected

Definition at line 926 of file G4PolyPhiFace.cc.

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

Referenced by InCone(), and IntersectProp().

LeftOn()

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

protected

Definition at line 935 of file G4PolyPhiFace.cc.

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

Referenced by InCone().

Normal()

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

overridevirtual

Implements G4VCSGface.

Definition at line 540 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, nullptr ))
  {
    //
    // 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 ( )

overridevirtual

Implements G4VCSGface.

Definition at line 896 of file G4PolyPhiFace.cc.

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

SurfaceTriangle()

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

Definition at line 878 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 1096 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().
 
  const std::size_t maxEdges = (numEdges > 0) ? numEdges : 1;
  auto tri_help = new G4PolyPhiFaceVertex[maxEdges];
  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+maxEdges-1; }
    else
      { triang->prev=helper2; }
 
    // add pointer on next corner
    //
    if( helper<corners+maxEdges-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 != nullptr)
  {
     // 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 218 of file G4PolyPhiFace.hh.

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

corners

G4PolyPhiFaceVertex* G4PolyPhiFace::corners = nullptr

protected

Definition at line 210 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 209 of file G4PolyPhiFace.hh.

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

fSurfaceArea

G4double G4PolyPhiFace::fSurfaceArea = 0.0

protected

Definition at line 221 of file G4PolyPhiFace.hh.

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

kCarTolerance

G4double G4PolyPhiFace::kCarTolerance

protected

Definition at line 220 of file G4PolyPhiFace.hh.

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

normal

G4ThreeVector G4PolyPhiFace::normal

protected

Definition at line 211 of file G4PolyPhiFace.hh.

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

numEdges

G4int G4PolyPhiFace::numEdges = 0

protected

Definition at line 208 of file G4PolyPhiFace.hh.

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

radial

G4ThreeVector G4PolyPhiFace::radial

protected

Definition at line 212 of file G4PolyPhiFace.hh.

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

rMax

G4double G4PolyPhiFace::rMax

protected

Definition at line 216 of file G4PolyPhiFace.hh.

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

rMin

G4double G4PolyPhiFace::rMin

protected

Definition at line 216 of file G4PolyPhiFace.hh.

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

surface

G4ThreeVector G4PolyPhiFace::surface

protected

Definition at line 213 of file G4PolyPhiFace.hh.

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

surface_point

G4ThreeVector G4PolyPhiFace::surface_point

protected

Definition at line 214 of file G4PolyPhiFace.hh.

Referenced by GetPointOnFace(), and Triangulate().

triangles

G4PolyPhiFaceVertex* G4PolyPhiFace::triangles = nullptr

protected

Definition at line 222 of file G4PolyPhiFace.hh.

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

zMax

G4double G4PolyPhiFace::zMax

protected

Definition at line 217 of file G4PolyPhiFace.hh.

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

zMin

G4double G4PolyPhiFace::zMin

protected

Definition at line 217 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