G4SmartVoxelHeader Class Reference

#include <G4SmartVoxelHeader.hh>

Public Member Functions
	G4SmartVoxelHeader (G4LogicalVolume *pVolume, G4int pSlice=0)
	~G4SmartVoxelHeader ()
G4int	GetMaxEquivalentSliceNo () const
void	SetMaxEquivalentSliceNo (G4int pMax)
G4int	GetMinEquivalentSliceNo () const
void	SetMinEquivalentSliceNo (G4int pMin)
EAxis	GetAxis () const
EAxis	GetParamAxis () const
G4double	GetMaxExtent () const
G4double	GetMinExtent () const
std::size_t	GetNoSlices () const
G4SmartVoxelProxy *	GetSlice (std::size_t n) const
G4bool	AllSlicesEqual () const
G4bool	operator== (const G4SmartVoxelHeader &pHead) const
	G4SmartVoxelHeader (G4LogicalVolume *pVolume, const G4VoxelLimits &pLimits, const G4VolumeNosVector *pCandidates, G4int pSlice=0)

Protected Member Functions
void	BuildVoxels (G4LogicalVolume *pVolume)
void	BuildReplicaVoxels (G4LogicalVolume *pVolume)
void	BuildConsumedNodes (G4int nReplicas)
void	BuildVoxelsWithinLimits (G4LogicalVolume *pVolume, G4VoxelLimits pLimits, const G4VolumeNosVector *pCandidates)
void	BuildEquivalentSliceNos ()
void	CollectEquivalentNodes ()
void	CollectEquivalentHeaders ()
G4ProxyVector *	BuildNodes (G4LogicalVolume *pVolume, G4VoxelLimits pLimits, const G4VolumeNosVector *pCandidates, EAxis pAxis)
G4double	CalculateQuality (G4ProxyVector *pSlice)
void	RefineNodes (G4LogicalVolume *pVolume, G4VoxelLimits pLimits)

Protected Attributes
G4int	fminEquivalent
G4int	fmaxEquivalent
EAxis	faxis
EAxis	fparamAxis
G4double	fmaxExtent
G4double	fminExtent
G4ProxyVector	fslices

Friends
std::ostream &	operator<< (std::ostream &s, const G4SmartVoxelHeader &h)

Detailed Description

Definition at line 74 of file G4SmartVoxelHeader.hh.

Constructor & Destructor Documentation

G4SmartVoxelHeader() [1/2]

G4SmartVoxelHeader::G4SmartVoxelHeader	(	G4LogicalVolume *	pVolume,
		G4int	pSlice = 0 )

Definition at line 59 of file G4SmartVoxelHeader.cc.

  : fminEquivalent(pSlice),
    fmaxEquivalent(pSlice),
    fparamAxis(kUndefined)
{
  std::size_t nDaughters = pVolume->GetNoDaughters();
 
  // Determine whether daughter is replicated
  //
  if ((nDaughters!=1) || (!pVolume->GetDaughter(0)->IsReplicated()))
  {
    // Daughter not replicated => conventional voxel Build
    // where each daughters extents are computed
    //
    BuildVoxels(pVolume);
  }
  else
  {
    // Single replicated daughter
    //
    BuildReplicaVoxels(pVolume);
  }
}

Referenced by CollectEquivalentHeaders(), GetSlice(), operator<<, operator==(), RefineNodes(), and ~G4SmartVoxelHeader().

~G4SmartVoxelHeader()

G4SmartVoxelHeader::~G4SmartVoxelHeader

(

)

Definition at line 119 of file G4SmartVoxelHeader.cc.

{
  // Manually destroy underlying nodes/headers
  // Delete collected headers and nodes once only
  //
  std::size_t node, proxy, maxNode=fslices.size();
  G4SmartVoxelProxy* lastProxy = nullptr;
  G4SmartVoxelNode *dyingNode, *lastNode=nullptr;
  G4SmartVoxelHeader *dyingHeader, *lastHeader=nullptr;
 
  for (node=0; node<maxNode; ++node)
  {
    if (fslices[node]->IsHeader())
    {
      dyingHeader = fslices[node]->GetHeader();
      if (lastHeader != dyingHeader)
      {
        lastHeader = dyingHeader;
        lastNode = nullptr;
        delete dyingHeader;
      }
    }
    else
    {
      dyingNode = fslices[node]->GetNode();
      if (dyingNode != lastNode)
      {
        lastNode = dyingNode;
        lastHeader = nullptr;
        delete dyingNode;
      }
    }
  }
  // Delete proxies
  //
  for (proxy=0; proxy<maxNode; ++proxy)
  {
    if (fslices[proxy] != lastProxy)
    {
      lastProxy = fslices[proxy];
      delete lastProxy;
    }
  }
  // Don't need to clear slices
  // fslices.clear();
}

G4SmartVoxelHeader() [2/2]

G4SmartVoxelHeader::G4SmartVoxelHeader	(	G4LogicalVolume *	pVolume,
		const G4VoxelLimits &	pLimits,
		const G4VolumeNosVector *	pCandidates,
		G4int	pSlice = 0 )

Definition at line 92 of file G4SmartVoxelHeader.cc.

  : fminEquivalent(pSlice),
    fmaxEquivalent(pSlice),
    fparamAxis(kUndefined)
{
#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "**** G4SmartVoxelHeader::G4SmartVoxelHeader" << G4endl
         << "     Limits " << pLimits << G4endl
         << "     Candidate #s = " ;
  for (auto i=0; i<pCandidates->size(); ++i)
  {
    G4cout << (*pCandidates)[i] << " ";
  }
  G4cout << G4endl;
#endif   
 
  BuildVoxelsWithinLimits(pVolume,pLimits,pCandidates);
}

Member Function Documentation

AllSlicesEqual()

G4bool G4SmartVoxelHeader::AllSlicesEqual

(

)

const

Definition at line 1257 of file G4SmartVoxelHeader.cc.

{
  std::size_t noSlices = fslices.size();
  G4SmartVoxelProxy* refProxy;
 
  if (noSlices>1)
  {
    refProxy=fslices[0];
    for (std::size_t i=1; i<noSlices; ++i)
    {
      if (refProxy!=fslices[i])
      {
        return false;
      }
    }
  }
  return true;
}

BuildConsumedNodes()

void G4SmartVoxelHeader::BuildConsumedNodes ( G4int nReplicas )

protected

Definition at line 387 of file G4SmartVoxelHeader.cc.

{
  G4int nNode, nVol;
  G4SmartVoxelNode* pNode;
  G4SmartVoxelProxy* pProxyNode;
 
  // Create and fill nodes in temporary G4NodeVector (on stack)
  //
  G4NodeVector nodeList;
  nodeList.reserve(nReplicas);
  for (nNode=0; nNode<nReplicas; ++nNode)
  {
    pNode = new G4SmartVoxelNode(nNode);
    if (pNode == nullptr)
    {
      G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "GeomMgt0003",
                  FatalException, "Node allocation error.");
    }
    nodeList.push_back(pNode);
  }
  for (nVol=0; nVol<nReplicas; ++nVol)
  {
    nodeList[nVol]->Insert(nVol);   // Insert replication of number
  }                                 // identical to voxel number
 
  // Create & fill proxy List `in place' by modifying instance data fslices
  //
  fslices.clear();
  for (nNode=0; nNode<nReplicas; ++nNode)
  {
    pProxyNode = new G4SmartVoxelProxy(nodeList[nNode]);
    if (pProxyNode == nullptr)
    {
      G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "GeomMgt0003",
                  FatalException, "Proxy node allocation error.");
    }
    fslices.push_back(pProxyNode);
  }
}

Referenced by BuildReplicaVoxels().

BuildEquivalentSliceNos()

void G4SmartVoxelHeader::BuildEquivalentSliceNos ( )

protected

Definition at line 559 of file G4SmartVoxelHeader.cc.

{
  std::size_t sliceNo, minNo, maxNo, equivNo;
  std::size_t maxNode = fslices.size();
  G4SmartVoxelNode *startNode, *sampleNode;
  for (sliceNo=0; sliceNo<maxNode; ++sliceNo)
  {
    minNo = sliceNo;
 
    // Get first node (see preconditions - will throw exception if a header)
    //
    startNode = fslices[minNo]->GetNode();
 
    // Find max equivalent
    //
    for (equivNo=minNo+1; equivNo<maxNode; ++equivNo)
    {
      sampleNode = fslices[equivNo]->GetNode();
      if (!((*startNode) == (*sampleNode))) { break; }
    }
    maxNo = equivNo-1;
    if (maxNo != minNo)
    {
      // Set min and max nos
      //
      for (equivNo=minNo; equivNo<=maxNo; ++equivNo)
      {
        sampleNode = fslices[equivNo]->GetNode();
        sampleNode->SetMinEquivalentSliceNo((G4int)minNo);
        sampleNode->SetMaxEquivalentSliceNo((G4int)maxNo);
      }
      // Advance outer loop to end of equivalent group
      //
      sliceNo = maxNo;
    }
  }
}

Referenced by BuildReplicaVoxels(), and BuildVoxelsWithinLimits().

BuildNodes()

G4ProxyVector * G4SmartVoxelHeader::BuildNodes ( G4LogicalVolume * pVolume, G4VoxelLimits pLimits, const G4VolumeNosVector * pCandidates, EAxis pAxis )

protected

Definition at line 726 of file G4SmartVoxelHeader.cc.

{
  G4double motherMinExtent= kInfinity, motherMaxExtent= -kInfinity,
           targetMinExtent= kInfinity, targetMaxExtent= -kInfinity;
  G4VPhysicalVolume* pDaughter = nullptr;
  G4VPVParameterisation* pParam = nullptr;
  G4VSolid *targetSolid;
  G4AffineTransform targetTransform;
  G4bool replicated;
  std::size_t nCandidates = pCandidates->size();
  std::size_t nVol, nNode, targetVolNo;
  G4VoxelLimits noLimits;
 
#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
         << "     Limits = " << pLimits << G4endl
         << "       Axis = " << pAxis << G4endl
         << " Candidates = " << nCandidates << G4endl;
#endif
 
  // Compute extent of logical volume's solid along this axis
  // NOTE: results stored locally and not preserved/reused
  //
  G4VSolid* outerSolid = pVolume->GetSolid();
  const G4AffineTransform origin;
  if( !outerSolid->CalculateExtent(pAxis, pLimits, origin,
                                   motherMinExtent, motherMaxExtent) )
  {
    outerSolid->CalculateExtent(pAxis, noLimits, origin,
                                motherMinExtent, motherMaxExtent);
  }
  G4VolumeExtentVector minExtents(nCandidates,0.);
  G4VolumeExtentVector maxExtents(nCandidates,0.);
 
  if ( (pVolume->GetNoDaughters() == 1)
    && (pVolume->GetDaughter(0)->IsReplicated()) )
  {
    // Replication data not required: only parameterisation object 
    // and volume no. List used
    //
    pDaughter = pVolume->GetDaughter(0);
    pParam = pDaughter->GetParameterisation();
    if (pParam == nullptr)
    {
      std::ostringstream message;
      message << "PANIC! - Missing parameterisation." << G4endl
              << "         Replicated volume with no parameterisation object !";
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                  FatalException, message);
      return nullptr;
    }
 
    // Setup daughter's transformations
    //
    targetTransform = G4AffineTransform(pDaughter->GetRotation(),
                                        pDaughter->GetTranslation());
    replicated = true;
  }
    else
  {
    replicated = false;
  }
    
  // Compute extents
  //
  for (nVol=0; nVol<nCandidates; ++nVol)
  {
    targetVolNo = (*pCandidates)[nVol];
    if (!replicated)
    {
      pDaughter = pVolume->GetDaughter(targetVolNo);
 
      // Setup daughter's transformations
      //
      targetTransform = G4AffineTransform(pDaughter->GetRotation(),
                                          pDaughter->GetTranslation());
      // Get underlying (and setup) solid
      //
      targetSolid = pDaughter->GetLogicalVolume()->GetSolid();
    }
    else
    {
      // Find  solid
      //
      targetSolid = pParam->ComputeSolid((G4int)targetVolNo,pDaughter);
 
      // Setup solid
      //
      targetSolid->ComputeDimensions(pParam,(G4int)targetVolNo,pDaughter);
 
      // Setup transform
      //
      pParam->ComputeTransformation((G4int)targetVolNo,pDaughter);
      targetTransform = G4AffineTransform(pDaughter->GetRotation(),
                                          pDaughter->GetTranslation());
    }
    // Calculate extents
    //
    if(!targetSolid->CalculateExtent(pAxis, pLimits, targetTransform,
                                     targetMinExtent, targetMaxExtent))
    {
      targetSolid->CalculateExtent(pAxis, noLimits, targetTransform,
                                   targetMinExtent,targetMaxExtent);
    }
    minExtents[nVol] = targetMinExtent;
    maxExtents[nVol] = targetMaxExtent;
 
#ifdef G4GEOMETRY_VOXELDEBUG
   G4cout << "---------------------------------------------------" << G4endl
          << "     Volume = " << pDaughter->GetName() << G4endl
          << " Min Extent = " << targetMinExtent << G4endl
          << " Max Extent = " << targetMaxExtent << G4endl
          << "---------------------------------------------------" << G4endl;
#endif
 
    // Check not entirely outside mother when processing toplevel nodes
    //
    if ( (!pLimits.IsLimited()) && ((targetMaxExtent<=motherMinExtent)
                                  ||(targetMinExtent>=motherMaxExtent)) )
    {
      std::ostringstream message;
      message << "PANIC! - Overlapping daughter with mother volume." << G4endl
              << "         Daughter physical volume "
              << pDaughter->GetName() << G4endl
              << "         is entirely outside mother logical volume "
              << pVolume->GetName() << " !!";
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0002",
                  FatalException, message);
    }
 
#ifdef G4GEOMETRY_VOXELDEBUG
    // Check for straddling volumes when debugging.
    // If a volume is >kStraddlePercent percent over the mother
    // boundary, print a warning.
    //
    if (!pLimits.IsLimited())
    {
      G4double width;
      G4int kStraddlePercent = 5;
      width = maxExtents[nVol]-minExtents[nVol];
      if ( (((motherMinExtent-minExtents[nVol])*100/width) > kStraddlePercent)
         ||(((maxExtents[nVol]-motherMaxExtent)*100/width) > kStraddlePercent) )
      {
        G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
               << "     WARNING : Daughter # " << nVol
               << " name = " << pDaughter->GetName() << G4endl
               << "     Crosses mother boundary of logical volume, name = " 
               << pVolume->GetName() << G4endl
               << "     by more than " << kStraddlePercent 
               << "%" << G4endl;
      }
    }
#endif
  }
 
  // Extents of all daughters known
 
  // Calculate minimum slice width, only including volumes inside the limits
  //
  G4double minWidth = kInfinity;
  G4double currentWidth;
  for (nVol=0; nVol<nCandidates; ++nVol)
  {
    // currentWidth should -always- be a positive value. Inaccurate computed extent
    // from the solid or situations of malformed geometries (overlaps) may lead to
    // negative values and therefore unpredictable crashes !
    //
    currentWidth = std::abs(maxExtents[nVol]-minExtents[nVol]);
    if ( (currentWidth<minWidth)
      && (maxExtents[nVol]>=pLimits.GetMinExtent(pAxis))
      && (minExtents[nVol]<=pLimits.GetMaxExtent(pAxis)) )
    {
      minWidth = currentWidth;
    }
  }
 
  // No. of Nodes formula - nearest integer to
  // mother width/half min daughter width +1
  //
  G4double noNodesExactD = ((motherMaxExtent-motherMinExtent)*2.0/minWidth)+1.0;
 
  // Compare with "smartless quality", i.e. the average number of slices
  // used per contained volume.
  //
  G4double smartlessComputed = noNodesExactD / nCandidates;
  G4double smartlessUser = pVolume->GetSmartless();
  G4double smartless = (smartlessComputed <= smartlessUser)
                       ? smartlessComputed : smartlessUser;
  G4double noNodesSmart = smartless*nCandidates;
  auto     noNodesExactI = G4int(noNodesSmart);
  G4long   noNodes = ((noNodesSmart-noNodesExactI)>=0.5)
                     ? noNodesExactI+1 : noNodesExactI;
  if( noNodes == 0 ) { noNodes=1; }
 
#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "     Smartless computed = " << smartlessComputed << G4endl
         << "     Smartless volume = " << smartlessUser
         << " => # Smartless = " << smartless << G4endl;
  G4cout << "     Min width = " << minWidth
         << " => # Nodes = " << noNodes << G4endl;
#endif
 
  if (noNodes>kMaxVoxelNodes)
  {
    noNodes=kMaxVoxelNodes;
#ifdef G4GEOMETRY_VOXELDEBUG
    G4cout << "     Nodes Clipped to = " << kMaxVoxelNodes << G4endl;
#endif   
  }
  G4double nodeWidth = (motherMaxExtent-motherMinExtent)/noNodes;
 
  // Create G4VoxelNodes. Will Add proxies before setting fslices
  //
  auto* nodeList = new G4NodeVector();
  if (nodeList == nullptr)
  {
    G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                FatalException, "NodeList allocation error.");
    return nullptr;
  }
  nodeList->reserve(noNodes);
 
  for (nNode=0; G4long(nNode)<noNodes; ++nNode)
  {
    G4SmartVoxelNode *pNode;
    pNode = new G4SmartVoxelNode((G4int)nNode);
    if (pNode == nullptr)
    {
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                  FatalException, "Node allocation error.");
      return nullptr;
    }
    nodeList->push_back(pNode);
  }
 
  // All nodes created (empty)
 
  // Fill nodes: Step through extent lists
  //
  for (nVol=0; nVol<nCandidates; ++nVol)
  {
    G4long nodeNo, minContainingNode, maxContainingNode;
    minContainingNode = (minExtents[nVol]-motherMinExtent)/nodeWidth;
    maxContainingNode = (maxExtents[nVol]-motherMinExtent)/nodeWidth;
 
    // Only add nodes that are inside the limits of the axis
    //
    if ( (maxContainingNode>=0) && (minContainingNode<noNodes) )
    {
      // If max extent is on max boundary => maxContainingNode=noNodes:
      // should be one less as nodeList has noNodes entries
      //
      if (maxContainingNode>=noNodes)
      {
        maxContainingNode = noNodes-1;
      }
      //
      // Protection against protruding volumes
      //
      if (minContainingNode<0)
      {
        minContainingNode = 0;
      }
      for (nodeNo=minContainingNode; nodeNo<=maxContainingNode; ++nodeNo)
      {
        (*nodeList)[nodeNo]->Insert((*pCandidates)[nVol]);
      }
    }
  }
 
  // All nodes filled
 
  // Create proxy List : caller has deletion responsibility
  // (but we must delete nodeList *itself* - not the contents)
  //
  auto* proxyList = new G4ProxyVector();
  if (proxyList == nullptr)
  {
    G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                FatalException, "Proxy list allocation error.");
    return nullptr;
  }
  proxyList->reserve(noNodes);
 
  //
  // Fill proxy List
  //
  for (nNode=0; G4long(nNode)<noNodes; ++nNode)
  {
    // Get rid of possible excess capacity in the internal node vector
    //
    ((*nodeList)[nNode])->Shrink();
    auto* pProxyNode = new G4SmartVoxelProxy((*nodeList)[nNode]);
    if (pProxyNode == nullptr)
    {
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                  FatalException, "Proxy node allocation failed.");
      return nullptr;
    }
    proxyList->push_back(pProxyNode);
  }
  delete nodeList;
  return proxyList;
}

Referenced by BuildReplicaVoxels(), and BuildVoxelsWithinLimits().

BuildReplicaVoxels()

void G4SmartVoxelHeader::BuildReplicaVoxels ( G4LogicalVolume * pVolume )

protected

Definition at line 259 of file G4SmartVoxelHeader.cc.

{
  G4VPhysicalVolume* pDaughter = nullptr;
 
  // Replication data
  //
  EAxis axis;
  G4int nReplicas;
  G4double width,offset;
  G4bool consuming;
 
  // Consistency check: pVolume should contain single replicated volume
  //
  if ( (pVolume->GetNoDaughters()==1)
    && (pVolume->GetDaughter(0)->IsReplicated()) )
  {
    // Obtain replication data
    //
    pDaughter = pVolume->GetDaughter(0);
    pDaughter->GetReplicationData(axis,nReplicas,width,offset,consuming);
    fparamAxis = axis;
    if ( !consuming )
    {
      G4VoxelLimits limits;   // Create `unlimited' limits object
      G4VolumeNosVector targetList;
      targetList.reserve(nReplicas);
      for (auto i=0; i<nReplicas; ++i)
      {
        targetList.push_back(i);
      }
      if (axis != kUndefined)
      {
        // Apply voxelisation along the specified axis only
 
        G4ProxyVector* pSlices=BuildNodes(pVolume,limits,&targetList,axis);
        faxis = axis;
        fslices = *pSlices;
        delete pSlices;
 
        // Calculate and set min and max extents given our axis
        //
        const G4AffineTransform origin;
        pVolume->GetSolid()->CalculateExtent(faxis, limits, origin,
                                             fminExtent, fmaxExtent);
        // Calculate equivalent nos
        //
        BuildEquivalentSliceNos();
        CollectEquivalentNodes();   // Collect common nodes
      }
      else
      {
        // Build voxels similarly as for normal placements considering
        // all three cartesian axes.
 
        BuildVoxelsWithinLimits(pVolume, limits, &targetList);
      }
    }
    else
    {
      // Replication is consuming -> Build voxels directly
      //
      // o Cartesian axes - range is -width*nREplicas/2 to +width*nREplicas/2
      //                    nReplicas replications result
      // o Radial axis (rho) = range is 0 to width*nReplicas
      //                    nReplicas replications result
      // o Phi axi       - range is offset to offset+width*nReplicas radians
      //
      // Equivalent slices no computation & collection not required - all
      // slices are different
      //
      switch (axis)
      {
        case kXAxis:
        case kYAxis:
        case kZAxis:
          fminExtent = -width*nReplicas*0.5;
          fmaxExtent =  width*nReplicas*0.5;
          break;
        case kRho:
          fminExtent = offset;
          fmaxExtent = width*nReplicas+offset;
          break;
        case kPhi:
          fminExtent = offset;
          fmaxExtent = offset+width*nReplicas;
          break;
        default:
          G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels()",
                      "GeomMgt0002", FatalException, "Illegal axis.");
          break;
      }  
      faxis = axis;   // Set axis
      BuildConsumedNodes(nReplicas);
      if ( (axis==kXAxis) || (axis==kYAxis) || (axis==kZAxis) )
      {
        // Sanity check on extent
        //
        G4double emin = kInfinity, emax = -kInfinity;
        G4VoxelLimits limits;
        G4AffineTransform origin;
        pVolume->GetSolid()->CalculateExtent(axis, limits, origin, emin, emax);
        if ( (std::fabs((emin-fminExtent)/fminExtent) +
              std::fabs((emax-fmaxExtent)/fmaxExtent)) > 0.05)
        {
          std::ostringstream message;
          message << "Sanity check: wrong solid extent." << G4endl
                  << "        Replicated geometry, logical volume: "
                  << pVolume->GetName();
          G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels",
                      "GeomMgt0002", FatalException, message);
        }
      }
    }
  }
  else
  {
    G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels", "GeomMgt0002",
                FatalException, "Only one replicated daughter is allowed !");
  }
}

BuildVoxels()

void G4SmartVoxelHeader::BuildVoxels ( G4LogicalVolume * pVolume )

protected

Definition at line 239 of file G4SmartVoxelHeader.cc.

{
  G4VoxelLimits limits;   // Create `unlimited' limits object
  std::size_t nDaughters = pVolume->GetNoDaughters();
 
  G4VolumeNosVector targetList;
  targetList.reserve(nDaughters);
  for (std::size_t i=0; i<nDaughters; ++i)
  {
    targetList.push_back((G4int)i);
  }
  BuildVoxelsWithinLimits(pVolume, limits, &targetList);
}

BuildVoxelsWithinLimits()

void G4SmartVoxelHeader::BuildVoxelsWithinLimits ( G4LogicalVolume * pVolume, G4VoxelLimits pLimits, const G4VolumeNosVector * pCandidates )

protected

Definition at line 435 of file G4SmartVoxelHeader.cc.

{
  // Choose best axis for slicing by:
  // 1. Trying all unlimited cartesian axes
  // 2. Select axis which gives greatest no slices
 
  G4ProxyVector *pGoodSlices=nullptr, *pTestSlices;
  G4double goodSliceScore=kInfinity, testSliceScore;
  EAxis goodSliceAxis = kXAxis;
  std::size_t node, maxNode; 
  G4VoxelLimits noLimits;
 
  // Try all non-limited cartesian axes
  //
  for ( EAxis testAxis : { kXAxis, kYAxis, kZAxis } )
  {
    if (!pLimits.IsLimited(testAxis))
    {
      pTestSlices = BuildNodes(pVolume,pLimits,pCandidates,testAxis);
      testSliceScore = CalculateQuality(pTestSlices);
      if ( (pGoodSlices == nullptr) || (testSliceScore<goodSliceScore) )
      {
        goodSliceAxis  = testAxis;
        goodSliceScore = testSliceScore;
        std::swap( pGoodSlices, pTestSlices);
      }
      if (pTestSlices != nullptr)
      {
        // Destroy pTestSlices and all its contents
        //
        maxNode=pTestSlices->size();
        for (node=0; node<maxNode; ++node)
        {
          delete (*pTestSlices)[node]->GetNode();
        }
        G4SmartVoxelProxy* tmpProx;
        while (!pTestSlices->empty())  // Loop checking, 06.08.2015, G.Cosmo
        {
          tmpProx = pTestSlices->back();
          pTestSlices->pop_back();
          for (auto i=pTestSlices->cbegin(); i!=pTestSlices->cend(); )
          {
            if (*i==tmpProx)
            {
              i = pTestSlices->erase(i);
            }
            else
            {
              ++i;
            }
          }
          delete tmpProx; 
        }
        delete pTestSlices;
      }
    }
  }
  // Check for error case.. when limits already 3d,
  // so cannot select a new axis
  //
  if (pGoodSlices == nullptr)
  {
    G4Exception("G4SmartVoxelHeader::BuildVoxelsWithinLimits()",
                "GeomMgt0002", FatalException,
                "Cannot select more than 3 axis for optimisation.");
    return;
  }
 
  // 
  // We have selected pGoodSlices, with a score testSliceScore
  //
 
  // Store chosen axis, slice ptr
  //
  fslices =* pGoodSlices; // Set slice information, copy ptrs in collection
  delete pGoodSlices;     // Destroy slices vector, but not contained
                          // proxies or nodes
  faxis = goodSliceAxis;
 
#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << G4endl << "     Volume = " << pVolume->GetName()
         << G4endl << "     Selected axis = " << faxis << G4endl;
  for (auto islice=0; islice<fslices.size(); ++islice)
  {
    G4cout << "     Node #" << islice << " = {";
    for (auto j=0; j<fslices[islice]->GetNode()->GetNoContained(); ++j)
    {
      G4cout << " " << fslices[islice]->GetNode()->GetVolume(j);
    }
    G4cout << " }" << G4endl;
  }
  G4cout << G4endl;
#endif
 
  // Calculate and set min and max extents given our axis
  //
  G4VSolid* outerSolid = pVolume->GetSolid();
  const G4AffineTransform origin;
  if(!outerSolid->CalculateExtent(faxis,pLimits,origin,fminExtent,fmaxExtent))
  {
    outerSolid->CalculateExtent(faxis,noLimits,origin,fminExtent,fmaxExtent);
  }
 
  // Calculate equivalent nos
  //
  BuildEquivalentSliceNos();
  CollectEquivalentNodes();      // Collect common nodes
  RefineNodes(pVolume, pLimits); // Refine nodes creating headers
 
  // No common headers can exist because collapsed by construction
}

Referenced by BuildReplicaVoxels(), BuildVoxels(), and G4SmartVoxelHeader().

CalculateQuality()

G4double G4SmartVoxelHeader::CalculateQuality ( G4ProxyVector * pSlice )

protected

Definition at line 1049 of file G4SmartVoxelHeader.cc.

{
  G4double quality;
  std::size_t nNodes = pSlice->size();
  std::size_t noContained, maxContained=0, sumContained=0, sumNonEmptyNodes=0;
  G4SmartVoxelNode *node;
 
  for (std::size_t i=0; i<nNodes; ++i)
  {
    if ((*pSlice)[i]->IsNode())
    {
      // Definitely a node. Add info to running totals
      //
      node = (*pSlice)[i]->GetNode();
      noContained = node->GetNoContained();
      if (noContained != 0)
      {
        ++sumNonEmptyNodes;
        sumContained += noContained;
        //
        // Calc maxContained for statistics
        //
        if (noContained>maxContained)
        {
          maxContained = noContained;
        }
      }
    }
    else
    {
      G4Exception("G4SmartVoxelHeader::CalculateQuality()", "GeomMgt0001",
                  FatalException, "Not applicable to replicated volumes.");
    }
  }
 
  // Calculate quality with protection against no non-empty nodes
  //
  if (sumNonEmptyNodes != 0)
  {
    quality = sumContained/sumNonEmptyNodes;
  }
  else
  {
    quality = kInfinity;
  }
 
#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "**** G4SmartVoxelHeader::CalculateQuality" << G4endl
         << "     Quality = " << quality << G4endl
         << "     Nodes = " << nNodes 
         << " of which " << sumNonEmptyNodes << " non empty" << G4endl
         << "     Max Contained = " << maxContained << G4endl;
#endif
 
  return quality;
}

Referenced by BuildVoxelsWithinLimits().

CollectEquivalentHeaders()

void G4SmartVoxelHeader::CollectEquivalentHeaders ( )

protected

Definition at line 653 of file G4SmartVoxelHeader.cc.

{
  std::size_t sliceNo, maxNo, equivNo;
  std::size_t maxNode = fslices.size();
  G4SmartVoxelHeader *equivHeader, *sampleHeader;
  G4SmartVoxelProxy *equivProxy;
 
  for (sliceNo=0; sliceNo<maxNode; ++sliceNo)
  {
    equivProxy = fslices[sliceNo];
    if (equivProxy->IsHeader())
    {
      equivHeader = equivProxy->GetHeader();
      maxNo = equivHeader->GetMaxEquivalentSliceNo();
      if (maxNo != sliceNo)
      {
        // Attempt collection between sliceNo and maxNo inclusive:
        // look for common headers. All slices between sliceNo and maxNo
        // are guaranteed to be headers but may not have equal contents
        //
#ifdef G4GEOMETRY_VOXELDEBUG
        G4cout << "**** G4SmartVoxelHeader::CollectEquivalentHeaders" << G4endl
               << "     Collecting Headers =";
#endif
        for (equivNo=sliceNo+1; equivNo<=maxNo; ++equivNo)
        {
          sampleHeader = fslices[equivNo]->GetHeader();
          if ( (*sampleHeader) == (*equivHeader) )
          {
#ifdef G4GEOMETRY_VOXELDEBUG
            G4cout << " " << equivNo;
#endif
            // Delete sampleHeader + proxy and replace with equivHeader/Proxy
            //
            delete sampleHeader;
            delete fslices[equivNo];
            fslices[equivNo] = equivProxy;
          }
          else
          {
            // Not equal. Set this header to be
            // the current header for comparisons
            //
            equivProxy  = fslices[equivNo];
            equivHeader = equivProxy->GetHeader();
          }
 
        }
#ifdef G4GEOMETRY_VOXELDEBUG
        G4cout << G4endl;
#endif
        // Skip past examined slices
        //
        sliceNo = maxNo;
      }
    }
  }
}

CollectEquivalentNodes()

void G4SmartVoxelHeader::CollectEquivalentNodes ( )

protected

Definition at line 606 of file G4SmartVoxelHeader.cc.

{
  std::size_t sliceNo, maxNo, equivNo;
  std::size_t maxNode=fslices.size();
  G4SmartVoxelNode* equivNode;
  G4SmartVoxelProxy* equivProxy;
 
  for (sliceNo=0; sliceNo<maxNode; ++sliceNo)
  {
    equivProxy=fslices[sliceNo];
 
    // Assumption (see preconditions): all slices are nodes
    //
    equivNode = equivProxy->GetNode();
    maxNo = equivNode->GetMaxEquivalentSliceNo();
    if (maxNo != sliceNo)
    {
#ifdef G4GEOMETRY_VOXELDEBUG
      G4cout << "**** G4SmartVoxelHeader::CollectEquivalentNodes" << G4endl
             << "     Collecting Nodes = " 
             << sliceNo << " - " << maxNo << G4endl;
#endif
      // Do collection between sliceNo and maxNo inclusive
      //
      for (equivNo=sliceNo+1; equivNo<=maxNo; ++equivNo)
      {
        delete fslices[equivNo]->GetNode();
        delete fslices[equivNo];
        fslices[equivNo] = equivProxy;
      }
      sliceNo = maxNo;
    }
  }
}

Referenced by BuildReplicaVoxels(), and BuildVoxelsWithinLimits().

GetAxis()

EAxis G4SmartVoxelHeader::GetAxis

(

)

const

Referenced by G4VoxelNavigation::LocateNextVoxel(), operator==(), and G4VoxelSafety::SafetyForVoxelHeader().

GetMaxEquivalentSliceNo()

G4int G4SmartVoxelHeader::GetMaxEquivalentSliceNo

(

)

const

Referenced by CollectEquivalentHeaders(), and G4VoxelSafety::SafetyForVoxelHeader().

GetMaxExtent()

G4double G4SmartVoxelHeader::GetMaxExtent

(

)

const

Referenced by operator==(), and G4VoxelSafety::SafetyForVoxelHeader().

GetMinEquivalentSliceNo()

G4int G4SmartVoxelHeader::GetMinEquivalentSliceNo

(

)

const

Referenced by G4VoxelSafety::SafetyForVoxelHeader().

GetMinExtent()

G4double G4SmartVoxelHeader::GetMinExtent

(

)

const

Referenced by G4VoxelNavigation::ComputeVoxelSafety(), G4VoxelNavigation::LocateNextVoxel(), operator==(), and G4VoxelSafety::SafetyForVoxelHeader().

GetNoSlices()

std::size_t G4SmartVoxelHeader::GetNoSlices

(

)

const

Referenced by G4SmartVoxelStat::CountHeadsAndNodes(), operator==(), and G4VoxelSafety::SafetyForVoxelHeader().

GetParamAxis()

EAxis G4SmartVoxelHeader::GetParamAxis

(

)

const

GetSlice()

G4SmartVoxelProxy * G4SmartVoxelHeader::GetSlice

(

std::size_t

n

)

const

Referenced by G4SmartVoxelStat::CountHeadsAndNodes(), operator==(), and G4VoxelSafety::SafetyForVoxelHeader().

operator==()

G4bool G4SmartVoxelHeader::operator==

(

const G4SmartVoxelHeader &

pHead

)

const

Definition at line 175 of file G4SmartVoxelHeader.cc.

{
  if ( (GetAxis()      == pHead.GetAxis())
    && (GetNoSlices()  == pHead.GetNoSlices())
    && (GetMinExtent() == pHead.GetMinExtent())
    && (GetMaxExtent() == pHead.GetMaxExtent()) )
  {
    std::size_t node, maxNode;
    G4SmartVoxelProxy *leftProxy, *rightProxy;
    G4SmartVoxelHeader *leftHeader, *rightHeader;
    G4SmartVoxelNode *leftNode, *rightNode;
 
    maxNode = GetNoSlices();
    for (node=0; node<maxNode; ++node)
    {
      leftProxy  = GetSlice(node);
      rightProxy = pHead.GetSlice(node);
      if (leftProxy->IsHeader())
      {
        if (rightProxy->IsNode())
        {
          return false;
        }
        else
        {
          leftHeader  = leftProxy->GetHeader();
          rightHeader = rightProxy->GetHeader();
          if (!(*leftHeader == *rightHeader))
          {
            return false;
          }
        }
      }
      else
      {
        if (rightProxy->IsHeader())
        {
          return false;
        }
        else
        {
          leftNode  = leftProxy->GetNode();
          rightNode = rightProxy->GetNode();
          if (!(*leftNode == *rightNode))
          {
            return false;
          }
        }
      }
    }
    return true;
  }
  else
  {
    return false;
  }
}

RefineNodes()

void G4SmartVoxelHeader::RefineNodes ( G4LogicalVolume * pVolume, G4VoxelLimits pLimits )

protected

Definition at line 1115 of file G4SmartVoxelHeader.cc.

{
  std::size_t refinedDepth=0, minVolumes;
  std::size_t maxNode = fslices.size();
 
  if (pLimits.IsXLimited()) 
  {
    ++refinedDepth;
  }
  if (pLimits.IsYLimited()) 
  {
    ++refinedDepth;
  }
  if (pLimits.IsZLimited()) 
  {
    ++refinedDepth;
  }
 
  // Calculate minimum number of volumes necessary to refine
  //
  switch (refinedDepth)
  {
    case 0:
      minVolumes=kMinVoxelVolumesLevel2;
      break;
    case 1:
      minVolumes=kMinVoxelVolumesLevel3;
      break;
    default:
      minVolumes=10000;   // catch refinedDepth=3 and errors
      break;
  }
 
  if (refinedDepth<2)
  {
    std::size_t targetNo, noContainedDaughters, minNo, maxNo, replaceNo, i;
    G4double sliceWidth = (fmaxExtent-fminExtent)/maxNode;
    G4VoxelLimits newLimits;
    G4SmartVoxelNode* targetNode;
    G4SmartVoxelProxy* targetNodeProxy;
    G4SmartVoxelHeader* replaceHeader;
    G4SmartVoxelProxy* replaceHeaderProxy;
    G4VolumeNosVector* targetList;
    G4SmartVoxelProxy* lastProxy;
      
    for (targetNo=0; targetNo<maxNode; ++targetNo)
    {
      // Assume all slices are nodes (see preconditions)
      //
      targetNodeProxy = fslices[targetNo];
      targetNode = targetNodeProxy->GetNode();
 
      if (targetNode->GetNoContained() >= minVolumes)
      {
        noContainedDaughters = targetNode->GetNoContained();
        targetList = new G4VolumeNosVector();
        if (targetList == nullptr)
        {
          G4Exception("G4SmartVoxelHeader::RefineNodes()",
                      "GeomMgt0003", FatalException,
                      "Target volume node list allocation error.");
          return;
        }
        targetList->reserve(noContainedDaughters);
        for (i=0; i<noContainedDaughters; ++i)
        {
          targetList->push_back(targetNode->GetVolume((G4int)i));
        }
        minNo = targetNode->GetMinEquivalentSliceNo();
        maxNo = targetNode->GetMaxEquivalentSliceNo();
 
#ifdef G4GEOMETRY_VOXELDEBUG
        G4cout << "**** G4SmartVoxelHeader::RefineNodes" << G4endl
               << "     Refining nodes " << minNo 
               << " - " << maxNo << " inclusive" << G4endl;
#endif
        if (minNo > maxNo)    // Delete node and list to be replaced
        {                     // and avoid further action ...
          delete targetNode;
          delete targetList;
          return;
        }
 
        // Delete node proxies at start of collected sets of nodes/headers
        //
        lastProxy=nullptr;
        for (replaceNo=minNo; replaceNo<=maxNo; ++replaceNo)
        {
          if (lastProxy != fslices[replaceNo])
          {
            lastProxy=fslices[replaceNo];
            delete lastProxy;
          }
        }
        // Delete node to be replaced
        //
        delete targetNode;
 
        // Create new headers + proxies and replace in fslices
        //
        newLimits = pLimits;
        newLimits.AddLimit(faxis,fminExtent+sliceWidth*minNo,
                           fminExtent+sliceWidth*(maxNo+1));
        replaceHeader = new G4SmartVoxelHeader(pVolume,newLimits,
                                               targetList,(G4int)replaceNo);
        if (replaceHeader == nullptr)
        {
          G4Exception("G4SmartVoxelHeader::RefineNodes()", "GeomMgt0003",
                      FatalException, "Refined VoxelHeader allocation error.");
          return;
        }
        replaceHeader->SetMinEquivalentSliceNo((G4int)minNo);
        replaceHeader->SetMaxEquivalentSliceNo((G4int)maxNo);
        replaceHeaderProxy = new G4SmartVoxelProxy(replaceHeader);
        if (replaceHeaderProxy == nullptr)
        {
          G4Exception("G4SmartVoxelHeader::RefineNodes()", "GeomMgt0003",
                      FatalException, "Refined VoxelProxy allocation error.");
          return;
        }
        for (replaceNo=minNo; replaceNo<=maxNo; ++replaceNo)
        {
          fslices[replaceNo] = replaceHeaderProxy;
        }
        // Finished replacing current `equivalent' group
        //
        delete targetList;
        targetNo=maxNo;
      }
    }
  }
}

Referenced by BuildVoxelsWithinLimits().

SetMaxEquivalentSliceNo()

void G4SmartVoxelHeader::SetMaxEquivalentSliceNo

(

G4int

pMax

)

Referenced by RefineNodes().

SetMinEquivalentSliceNo()

void G4SmartVoxelHeader::SetMinEquivalentSliceNo

(

G4int

pMin

)

Referenced by RefineNodes().

Friends And Related Symbol Documentation

operator<<

std::ostream & operator<< ( std::ostream & s, const G4SmartVoxelHeader & h )

friend

Definition at line 1280 of file G4SmartVoxelHeader.cc.

{
  os << "Axis = " << G4int(h.faxis) << G4endl;
  G4SmartVoxelProxy *collectNode=nullptr, *collectHead=nullptr;
  std::size_t collectNodeNo = 0;
  std::size_t collectHeadNo = 0;
  std::size_t i, j;
  G4bool haveHeaders = false;
 
  for (i=0; i<h.fslices.size(); ++i)
  {
    os << "Slice #" << i << " = ";
    if (h.fslices[i]->IsNode())
    {
      if (h.fslices[i]!=collectNode)
      {
        os << "{";
        for (std::size_t k=0; k<h.fslices[i]->GetNode()->GetNoContained(); ++k)
        {
          os << " " << h.fslices[i]->GetNode()->GetVolume((G4int)k);
        }
        os << " }" << G4endl;
        collectNode = h.fslices[i];
        collectNodeNo = i;
      }
      else
      {
        os << "As slice #" << collectNodeNo << G4endl;
      }
    }
    else
    {
      haveHeaders=true;
      if (h.fslices[i] != collectHead)
      {
        os << "Header" << G4endl;
        collectHead = h.fslices[i];
        collectHeadNo = i;
      }
      else
      {
        os << "As slice #" << collectHeadNo << G4endl;
      }
    }
  }
 
  if (haveHeaders)
  {
    collectHead=nullptr;
    for (j=0; j<h.fslices.size(); ++j)
    {
      if (h.fslices[j]->IsHeader())
      {
        os << "Header at Slice #" << j << " = ";
        if (h.fslices[j] != collectHead)
        {
          os << G4endl 
             << (*(h.fslices[j]->GetHeader()));
          collectHead = h.fslices[j];
          collectHeadNo = j;
        }
        else
        {
          os << "As slice #" << collectHeadNo << G4endl;
        }
      }
    }
  }
  return os;
}

Member Data Documentation

faxis

EAxis G4SmartVoxelHeader::faxis

protected

Definition at line 187 of file G4SmartVoxelHeader.hh.

Referenced by BuildReplicaVoxels(), BuildVoxelsWithinLimits(), operator<<, and RefineNodes().

fmaxEquivalent

G4int G4SmartVoxelHeader::fmaxEquivalent

protected

Definition at line 184 of file G4SmartVoxelHeader.hh.

Referenced by G4SmartVoxelHeader(), and G4SmartVoxelHeader().

fmaxExtent

G4double G4SmartVoxelHeader::fmaxExtent

protected

Definition at line 190 of file G4SmartVoxelHeader.hh.

Referenced by BuildReplicaVoxels(), BuildVoxelsWithinLimits(), and RefineNodes().

fminEquivalent

G4int G4SmartVoxelHeader::fminEquivalent

protected

Definition at line 183 of file G4SmartVoxelHeader.hh.

Referenced by G4SmartVoxelHeader(), and G4SmartVoxelHeader().

fminExtent

G4double G4SmartVoxelHeader::fminExtent

protected

Definition at line 191 of file G4SmartVoxelHeader.hh.

Referenced by BuildReplicaVoxels(), BuildVoxelsWithinLimits(), and RefineNodes().

fparamAxis

EAxis G4SmartVoxelHeader::fparamAxis

protected

Definition at line 187 of file G4SmartVoxelHeader.hh.

Referenced by BuildReplicaVoxels(), G4SmartVoxelHeader(), and G4SmartVoxelHeader().

fslices

G4ProxyVector G4SmartVoxelHeader::fslices

protected

Definition at line 194 of file G4SmartVoxelHeader.hh.

Referenced by AllSlicesEqual(), BuildConsumedNodes(), BuildEquivalentSliceNos(), BuildReplicaVoxels(), BuildVoxelsWithinLimits(), CollectEquivalentHeaders(), CollectEquivalentNodes(), operator<<, RefineNodes(), and ~G4SmartVoxelHeader().

---

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

• G4SmartVoxelHeader.hh
• G4SmartVoxelHeader.cc