| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // class G4SmartVoxelHeader implementation
27 //
28 // Define G4GEOMETRY_VOXELDEBUG for debugging information on G4cout
29 //
30 // 29.04.02 Use 3D voxelisation for non consuming replication - G.C.
31 // 18.04.01 Migrated to STL vector - G.C.
32 // 12.02.99 Introduction of new quality/smartless: max for (slices/cand) - S.G.
33 // 11.02.99 Voxels at lower levels are now built for collapsed slices - S.G.
34 // 21.07.95 Full implementation, supporting non divided physical volumes - P.K.
35 // 14.07.95 Initial version - stubb definitions only - P.K.
36 // --------------------------------------------------------------------
37
38 #include "G4SmartVoxelHeader.hh"
39
40 #include "G4ios.hh"
41
42 #include "G4LogicalVolume.hh"
43 #include "G4VPhysicalVolume.hh"
44 #include "G4VoxelLimits.hh"
45
46 #include "voxeldefs.hh"
47 #include "G4AffineTransform.hh"
48 #include "G4VSolid.hh"
49 #include "G4VPVParameterisation.hh"
50
51 // ***************************************************************************
52 // Constructor for topmost header, to begin voxel construction at a
53 // given logical volume.
54 // Constructs target List of volumes, calls "Build and refine" constructor.
55 // Assumes all daughters represent single volumes (ie. no divisions
56 // or parametric)
57 // ***************************************************************************
58 //
59 G4SmartVoxelHeader::G4SmartVoxelHeader(G4LogicalVolume* pVolume,
60 G4int pSlice)
61 : fminEquivalent(pSlice),
62 fmaxEquivalent(pSlice),
63 fparamAxis(kUndefined)
64 {
65 std::size_t nDaughters = pVolume->GetNoDaughters();
66
67 // Determine whether daughter is replicated
68 //
69 if ((nDaughters!=1) || (!pVolume->GetDaughter(0)->IsReplicated()))
70 {
71 // Daughter not replicated => conventional voxel Build
72 // where each daughters extents are computed
73 //
74 BuildVoxels(pVolume);
75 }
76 else
77 {
78 // Single replicated daughter
79 //
80 BuildReplicaVoxels(pVolume);
81 }
82 }
83
84 // ***************************************************************************
85 // Protected constructor:
86 // builds and refines voxels between specified limits, considering only
87 // the physical volumes numbered `pCandidates'. `pSlice' is used to set max
88 // and min equivalent slice nos for the header - they apply to the level
89 // of the header, not its nodes.
90 // ***************************************************************************
91 //
92 G4SmartVoxelHeader::G4SmartVoxelHeader(G4LogicalVolume* pVolume,
93 const G4VoxelLimits& pLimits,
94 const G4VolumeNosVector* pCandidates,
95 G4int pSlice)
96 : fminEquivalent(pSlice),
97 fmaxEquivalent(pSlice),
98 fparamAxis(kUndefined)
99 {
100 #ifdef G4GEOMETRY_VOXELDEBUG
101 G4cout << "**** G4SmartVoxelHeader::G4SmartVoxelHeader" << G4endl
102 << " Limits " << pLimits << G4endl
103 << " Candidate #s = " ;
104 for (auto i=0; i<pCandidates->size(); ++i)
105 {
106 G4cout << (*pCandidates)[i] << " ";
107 }
108 G4cout << G4endl;
109 #endif
110
111 BuildVoxelsWithinLimits(pVolume,pLimits,pCandidates);
112 }
113
114 // ***************************************************************************
115 // Destructor:
116 // deletes all proxies and underlying objects.
117 // ***************************************************************************
118 //
119 G4SmartVoxelHeader::~G4SmartVoxelHeader()
120 {
121 // Manually destroy underlying nodes/headers
122 // Delete collected headers and nodes once only
123 //
124 std::size_t node, proxy, maxNode=fslices.size();
125 G4SmartVoxelProxy* lastProxy = nullptr;
126 G4SmartVoxelNode *dyingNode, *lastNode=nullptr;
127 G4SmartVoxelHeader *dyingHeader, *lastHeader=nullptr;
128
129 for (node=0; node<maxNode; ++node)
130 {
131 if (fslices[node]->IsHeader())
132 {
133 dyingHeader = fslices[node]->GetHeader();
134 if (lastHeader != dyingHeader)
135 {
136 lastHeader = dyingHeader;
137 lastNode = nullptr;
138 delete dyingHeader;
139 }
140 }
141 else
142 {
143 dyingNode = fslices[node]->GetNode();
144 if (dyingNode != lastNode)
145 {
146 lastNode = dyingNode;
147 lastHeader = nullptr;
148 delete dyingNode;
149 }
150 }
151 }
152 // Delete proxies
153 //
154 for (proxy=0; proxy<maxNode; ++proxy)
155 {
156 if (fslices[proxy] != lastProxy)
157 {
158 lastProxy = fslices[proxy];
159 delete lastProxy;
160 }
161 }
162 // Don't need to clear slices
163 // fslices.clear();
164 }
165
166 // ***************************************************************************
167 // Equality operator: returns true if contents are equivalent.
168 // Implies a deep search through contained nodes/header.
169 // Compares headers' axes,sizes,extents. Returns false if different.
170 // For each contained proxy, determines whether node/header, compares and
171 // returns if different. Compares and returns if proxied nodes/headers
172 // are different.
173 // ***************************************************************************
174 //
175 G4bool G4SmartVoxelHeader::operator == (const G4SmartVoxelHeader& pHead) const
176 {
177 if ( (GetAxis() == pHead.GetAxis())
178 && (GetNoSlices() == pHead.GetNoSlices())
179 && (GetMinExtent() == pHead.GetMinExtent())
180 && (GetMaxExtent() == pHead.GetMaxExtent()) )
181 {
182 std::size_t node, maxNode;
183 G4SmartVoxelProxy *leftProxy, *rightProxy;
184 G4SmartVoxelHeader *leftHeader, *rightHeader;
185 G4SmartVoxelNode *leftNode, *rightNode;
186
187 maxNode = GetNoSlices();
188 for (node=0; node<maxNode; ++node)
189 {
190 leftProxy = GetSlice(node);
191 rightProxy = pHead.GetSlice(node);
192 if (leftProxy->IsHeader())
193 {
194 if (rightProxy->IsNode())
195 {
196 return false;
197 }
198 else
199 {
200 leftHeader = leftProxy->GetHeader();
201 rightHeader = rightProxy->GetHeader();
202 if (!(*leftHeader == *rightHeader))
203 {
204 return false;
205 }
206 }
207 }
208 else
209 {
210 if (rightProxy->IsHeader())
211 {
212 return false;
213 }
214 else
215 {
216 leftNode = leftProxy->GetNode();
217 rightNode = rightProxy->GetNode();
218 if (!(*leftNode == *rightNode))
219 {
220 return false;
221 }
222 }
223 }
224 }
225 return true;
226 }
227 else
228 {
229 return false;
230 }
231 }
232
233 // ***************************************************************************
234 // Builds voxels for daughters specified volume, in NON-REPLICATED case
235 // o Create List of target volume nos (all daughters; 0->noDaughters-1)
236 // o BuildWithinLimits does Build & also determines mother dimensions.
237 // ***************************************************************************
238 //
239 void G4SmartVoxelHeader::BuildVoxels(G4LogicalVolume* pVolume)
240 {
241 G4VoxelLimits limits; // Create `unlimited' limits object
242 std::size_t nDaughters = pVolume->GetNoDaughters();
243
244 G4VolumeNosVector targetList;
245 targetList.reserve(nDaughters);
246 for (std::size_t i=0; i<nDaughters; ++i)
247 {
248 targetList.push_back((G4int)i);
249 }
250 BuildVoxelsWithinLimits(pVolume, limits, &targetList);
251 }
252
253 // ***************************************************************************
254 // Builds voxels for specified volume containing a single replicated volume.
255 // If axis is not specified (i.e. "kUndefined"), 3D voxelisation is applied,
256 // and the best axis is determined according to heuristics as for placements.
257 // ***************************************************************************
258 //
259 void G4SmartVoxelHeader::BuildReplicaVoxels(G4LogicalVolume* pVolume)
260 {
261 G4VPhysicalVolume* pDaughter = nullptr;
262
263 // Replication data
264 //
265 EAxis axis;
266 G4int nReplicas;
267 G4double width,offset;
268 G4bool consuming;
269
270 // Consistency check: pVolume should contain single replicated volume
271 //
272 if ( (pVolume->GetNoDaughters()==1)
273 && (pVolume->GetDaughter(0)->IsReplicated()) )
274 {
275 // Obtain replication data
276 //
277 pDaughter = pVolume->GetDaughter(0);
278 pDaughter->GetReplicationData(axis,nReplicas,width,offset,consuming);
279 fparamAxis = axis;
280 if ( !consuming )
281 {
282 G4VoxelLimits limits; // Create `unlimited' limits object
283 G4VolumeNosVector targetList;
284 targetList.reserve(nReplicas);
285 for (auto i=0; i<nReplicas; ++i)
286 {
287 targetList.push_back(i);
288 }
289 if (axis != kUndefined)
290 {
291 // Apply voxelisation along the specified axis only
292
293 G4ProxyVector* pSlices=BuildNodes(pVolume,limits,&targetList,axis);
294 faxis = axis;
295 fslices = *pSlices;
296 delete pSlices;
297
298 // Calculate and set min and max extents given our axis
299 //
300 const G4AffineTransform origin;
301 pVolume->GetSolid()->CalculateExtent(faxis, limits, origin,
302 fminExtent, fmaxExtent);
303 // Calculate equivalent nos
304 //
305 BuildEquivalentSliceNos();
306 CollectEquivalentNodes(); // Collect common nodes
307 }
308 else
309 {
310 // Build voxels similarly as for normal placements considering
311 // all three cartesian axes.
312
313 BuildVoxelsWithinLimits(pVolume, limits, &targetList);
314 }
315 }
316 else
317 {
318 // Replication is consuming -> Build voxels directly
319 //
320 // o Cartesian axes - range is -width*nREplicas/2 to +width*nREplicas/2
321 // nReplicas replications result
322 // o Radial axis (rho) = range is 0 to width*nReplicas
323 // nReplicas replications result
324 // o Phi axi - range is offset to offset+width*nReplicas radians
325 //
326 // Equivalent slices no computation & collection not required - all
327 // slices are different
328 //
329 switch (axis)
330 {
331 case kXAxis:
332 case kYAxis:
333 case kZAxis:
334 fminExtent = -width*nReplicas*0.5;
335 fmaxExtent = width*nReplicas*0.5;
336 break;
337 case kRho:
338 fminExtent = offset;
339 fmaxExtent = width*nReplicas+offset;
340 break;
341 case kPhi:
342 fminExtent = offset;
343 fmaxExtent = offset+width*nReplicas;
344 break;
345 default:
346 G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels()",
347 "GeomMgt0002", FatalException, "Illegal axis.");
348 break;
349 }
350 faxis = axis; // Set axis
351 BuildConsumedNodes(nReplicas);
352 if ( (axis==kXAxis) || (axis==kYAxis) || (axis==kZAxis) )
353 {
354 // Sanity check on extent
355 //
356 G4double emin = kInfinity, emax = -kInfinity;
357 G4VoxelLimits limits;
358 G4AffineTransform origin;
359 pVolume->GetSolid()->CalculateExtent(axis, limits, origin, emin, emax);
360 if ( (std::fabs((emin-fminExtent)/fminExtent) +
361 std::fabs((emax-fmaxExtent)/fmaxExtent)) > 0.05)
362 {
363 std::ostringstream message;
364 message << "Sanity check: wrong solid extent." << G4endl
365 << " Replicated geometry, logical volume: "
366 << pVolume->GetName();
367 G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels",
368 "GeomMgt0002", FatalException, message);
369 }
370 }
371 }
372 }
373 else
374 {
375 G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels", "GeomMgt0002",
376 FatalException, "Only one replicated daughter is allowed !");
377 }
378 }
379
380 // ***************************************************************************
381 // Builds `consumed nodes': nReplicas nodes each containing one replication,
382 // numbered in sequence 0->nReplicas-1
383 // o Modifies fslices `in place'
384 // o faxis,fminExtent,fmaxExtent NOT modified.
385 // ***************************************************************************
386 //
387 void G4SmartVoxelHeader::BuildConsumedNodes(G4int nReplicas)
388 {
389 G4int nNode, nVol;
390 G4SmartVoxelNode* pNode;
391 G4SmartVoxelProxy* pProxyNode;
392
393 // Create and fill nodes in temporary G4NodeVector (on stack)
394 //
395 G4NodeVector nodeList;
396 nodeList.reserve(nReplicas);
397 for (nNode=0; nNode<nReplicas; ++nNode)
398 {
399 pNode = new G4SmartVoxelNode(nNode);
400 if (pNode == nullptr)
401 {
402 G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "GeomMgt0003",
403 FatalException, "Node allocation error.");
404 }
405 nodeList.push_back(pNode);
406 }
407 for (nVol=0; nVol<nReplicas; ++nVol)
408 {
409 nodeList[nVol]->Insert(nVol); // Insert replication of number
410 } // identical to voxel number
411
412 // Create & fill proxy List `in place' by modifying instance data fslices
413 //
414 fslices.clear();
415 for (nNode=0; nNode<nReplicas; ++nNode)
416 {
417 pProxyNode = new G4SmartVoxelProxy(nodeList[nNode]);
418 if (pProxyNode == nullptr)
419 {
420 G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "GeomMgt0003",
421 FatalException, "Proxy node allocation error.");
422 }
423 fslices.push_back(pProxyNode);
424 }
425 }
426
427 // ***************************************************************************
428 // Builds and refines voxels between specified limits, considering only
429 // the physical volumes numbered `pCandidates'.
430 // o Chooses axis
431 // o Determines min and max extents (of mother solid) within limits.
432 // ***************************************************************************
433 //
434 void
435 G4SmartVoxelHeader::BuildVoxelsWithinLimits(G4LogicalVolume* pVolume,
436 G4VoxelLimits pLimits,
437 const G4VolumeNosVector* pCandidates)
438 {
439 // Choose best axis for slicing by:
440 // 1. Trying all unlimited cartesian axes
441 // 2. Select axis which gives greatest no slices
442
443 G4ProxyVector *pGoodSlices=nullptr, *pTestSlices;
444 G4double goodSliceScore=kInfinity, testSliceScore;
445 EAxis goodSliceAxis = kXAxis;
446 std::size_t node, maxNode;
447 G4VoxelLimits noLimits;
448
449 // Try all non-limited cartesian axes
450 //
451 for ( EAxis testAxis : { kXAxis, kYAxis, kZAxis } )
452 {
453 if (!pLimits.IsLimited(testAxis))
454 {
455 pTestSlices = BuildNodes(pVolume,pLimits,pCandidates,testAxis);
456 testSliceScore = CalculateQuality(pTestSlices);
457 if ( (pGoodSlices == nullptr) || (testSliceScore<goodSliceScore) )
458 {
459 goodSliceAxis = testAxis;
460 goodSliceScore = testSliceScore;
461 std::swap( pGoodSlices, pTestSlices);
462 }
463 if (pTestSlices != nullptr)
464 {
465 // Destroy pTestSlices and all its contents
466 //
467 maxNode=pTestSlices->size();
468 for (node=0; node<maxNode; ++node)
469 {
470 delete (*pTestSlices)[node]->GetNode();
471 }
472 G4SmartVoxelProxy* tmpProx;
473 while (!pTestSlices->empty()) // Loop checking, 06.08.2015, G.Cosmo
474 {
475 tmpProx = pTestSlices->back();
476 pTestSlices->pop_back();
477 for (auto i=pTestSlices->cbegin(); i!=pTestSlices->cend(); )
478 {
479 if (*i==tmpProx)
480 {
481 i = pTestSlices->erase(i);
482 }
483 else
484 {
485 ++i;
486 }
487 }
488 delete tmpProx;
489 }
490 delete pTestSlices;
491 }
492 }
493 }
494 // Check for error case.. when limits already 3d,
495 // so cannot select a new axis
496 //
497 if (pGoodSlices == nullptr)
498 {
499 G4Exception("G4SmartVoxelHeader::BuildVoxelsWithinLimits()",
500 "GeomMgt0002", FatalException,
501 "Cannot select more than 3 axis for optimisation.");
502 return;
503 }
504
505 //
506 // We have selected pGoodSlices, with a score testSliceScore
507 //
508
509 // Store chosen axis, slice ptr
510 //
511 fslices =* pGoodSlices; // Set slice information, copy ptrs in collection
512 delete pGoodSlices; // Destroy slices vector, but not contained
513 // proxies or nodes
514 faxis = goodSliceAxis;
515
516 #ifdef G4GEOMETRY_VOXELDEBUG
517 G4cout << G4endl << " Volume = " << pVolume->GetName()
518 << G4endl << " Selected axis = " << faxis << G4endl;
519 for (auto islice=0; islice<fslices.size(); ++islice)
520 {
521 G4cout << " Node #" << islice << " = {";
522 for (auto j=0; j<fslices[islice]->GetNode()->GetNoContained(); ++j)
523 {
524 G4cout << " " << fslices[islice]->GetNode()->GetVolume(j);
525 }
526 G4cout << " }" << G4endl;
527 }
528 G4cout << G4endl;
529 #endif
530
531 // Calculate and set min and max extents given our axis
532 //
533 G4VSolid* outerSolid = pVolume->GetSolid();
534 const G4AffineTransform origin;
535 if(!outerSolid->CalculateExtent(faxis,pLimits,origin,fminExtent,fmaxExtent))
536 {
537 outerSolid->CalculateExtent(faxis,noLimits,origin,fminExtent,fmaxExtent);
538 }
539
540 // Calculate equivalent nos
541 //
542 BuildEquivalentSliceNos();
543 CollectEquivalentNodes(); // Collect common nodes
544 RefineNodes(pVolume, pLimits); // Refine nodes creating headers
545
546 // No common headers can exist because collapsed by construction
547 }
548
549 // ***************************************************************************
550 // Calculates and stores the minimum and maximum equivalent neighbour
551 // values for all slices at our level.
552 //
553 // Precondition: all slices are nodes.
554 // For each potential start of a group of equivalent nodes:
555 // o searches forwards in fslices to find group end
556 // o loops from start to end setting start and end slices.
557 // ***************************************************************************
558 //
559 void G4SmartVoxelHeader::BuildEquivalentSliceNos()
560 {
561 std::size_t sliceNo, minNo, maxNo, equivNo;
562 std::size_t maxNode = fslices.size();
563 G4SmartVoxelNode *startNode, *sampleNode;
564 for (sliceNo=0; sliceNo<maxNode; ++sliceNo)
565 {
566 minNo = sliceNo;
567
568 // Get first node (see preconditions - will throw exception if a header)
569 //
570 startNode = fslices[minNo]->GetNode();
571
572 // Find max equivalent
573 //
574 for (equivNo=minNo+1; equivNo<maxNode; ++equivNo)
575 {
576 sampleNode = fslices[equivNo]->GetNode();
577 if (!((*startNode) == (*sampleNode))) { break; }
578 }
579 maxNo = equivNo-1;
580 if (maxNo != minNo)
581 {
582 // Set min and max nos
583 //
584 for (equivNo=minNo; equivNo<=maxNo; ++equivNo)
585 {
586 sampleNode = fslices[equivNo]->GetNode();
587 sampleNode->SetMinEquivalentSliceNo((G4int)minNo);
588 sampleNode->SetMaxEquivalentSliceNo((G4int)maxNo);
589 }
590 // Advance outer loop to end of equivalent group
591 //
592 sliceNo = maxNo;
593 }
594 }
595 }
596
597 // ***************************************************************************
598 // Collects common nodes at our level, deleting all but one to save
599 // memory, and adjusting stored slice pointers appropriately.
600 //
601 // Preconditions:
602 // o the slices have not previously be "collected"
603 // o all of the slices are nodes.
604 // ***************************************************************************
605 //
606 void G4SmartVoxelHeader::CollectEquivalentNodes()
607 {
608 std::size_t sliceNo, maxNo, equivNo;
609 std::size_t maxNode=fslices.size();
610 G4SmartVoxelNode* equivNode;
611 G4SmartVoxelProxy* equivProxy;
612
613 for (sliceNo=0; sliceNo<maxNode; ++sliceNo)
614 {
615 equivProxy=fslices[sliceNo];
616
617 // Assumption (see preconditions): all slices are nodes
618 //
619 equivNode = equivProxy->GetNode();
620 maxNo = equivNode->GetMaxEquivalentSliceNo();
621 if (maxNo != sliceNo)
622 {
623 #ifdef G4GEOMETRY_VOXELDEBUG
624 G4cout << "**** G4SmartVoxelHeader::CollectEquivalentNodes" << G4endl
625 << " Collecting Nodes = "
626 << sliceNo << " - " << maxNo << G4endl;
627 #endif
628 // Do collection between sliceNo and maxNo inclusive
629 //
630 for (equivNo=sliceNo+1; equivNo<=maxNo; ++equivNo)
631 {
632 delete fslices[equivNo]->GetNode();
633 delete fslices[equivNo];
634 fslices[equivNo] = equivProxy;
635 }
636 sliceNo = maxNo;
637 }
638 }
639 }
640
641 // ***************************************************************************
642 // Collects common headers at our level, deleting all but one to save
643 // memory, and adjusting stored slice pointers appropriately.
644 //
645 // Preconditions:
646 // o if a header forms part of a range of equivalent slices
647 // (ie. GetMaxEquivalentSliceNo()>GetMinEquivalentSliceNo()),
648 // it is assumed that all slices in the range are headers.
649 // o this will be true if a constant Expression is used to evaluate
650 // when to refine nodes.
651 // ***************************************************************************
652 //
653 void G4SmartVoxelHeader::CollectEquivalentHeaders()
654 {
655 std::size_t sliceNo, maxNo, equivNo;
656 std::size_t maxNode = fslices.size();
657 G4SmartVoxelHeader *equivHeader, *sampleHeader;
658 G4SmartVoxelProxy *equivProxy;
659
660 for (sliceNo=0; sliceNo<maxNode; ++sliceNo)
661 {
662 equivProxy = fslices[sliceNo];
663 if (equivProxy->IsHeader())
664 {
665 equivHeader = equivProxy->GetHeader();
666 maxNo = equivHeader->GetMaxEquivalentSliceNo();
667 if (maxNo != sliceNo)
668 {
669 // Attempt collection between sliceNo and maxNo inclusive:
670 // look for common headers. All slices between sliceNo and maxNo
671 // are guaranteed to be headers but may not have equal contents
672 //
673 #ifdef G4GEOMETRY_VOXELDEBUG
674 G4cout << "**** G4SmartVoxelHeader::CollectEquivalentHeaders" << G4endl
675 << " Collecting Headers =";
676 #endif
677 for (equivNo=sliceNo+1; equivNo<=maxNo; ++equivNo)
678 {
679 sampleHeader = fslices[equivNo]->GetHeader();
680 if ( (*sampleHeader) == (*equivHeader) )
681 {
682 #ifdef G4GEOMETRY_VOXELDEBUG
683 G4cout << " " << equivNo;
684 #endif
685 // Delete sampleHeader + proxy and replace with equivHeader/Proxy
686 //
687 delete sampleHeader;
688 delete fslices[equivNo];
689 fslices[equivNo] = equivProxy;
690 }
691 else
692 {
693 // Not equal. Set this header to be
694 // the current header for comparisons
695 //
696 equivProxy = fslices[equivNo];
697 equivHeader = equivProxy->GetHeader();
698 }
699
700 }
701 #ifdef G4GEOMETRY_VOXELDEBUG
702 G4cout << G4endl;
703 #endif
704 // Skip past examined slices
705 //
706 sliceNo = maxNo;
707 }
708 }
709 }
710 }
711
712 // ***************************************************************************
713 // Builds the nodes corresponding to slices between the specified limits
714 // and along the specified axis, using candidate volume no.s in the vector
715 // pCandidates. If the `daughters' are replicated volumes (ie. the logical
716 // volume has a single replicated/parameterised volume for a daughter)
717 // the candidate no.s are interpreted as PARAMETERISED volume no.s &
718 // PARAMETERISATIONs are applied to compute transformations & solid
719 // dimensions appropriately. The volume must be parameterised - ie. has a
720 // parameterisation object & non-consuming) - in this case.
721 //
722 // Returns pointer to built node "structure" (guaranteed non NULL) consisting
723 // of G4SmartVoxelNodeProxies referring to G4SmartVoxelNodes.
724 // ***************************************************************************
725 //
726 G4ProxyVector* G4SmartVoxelHeader::BuildNodes(G4LogicalVolume* pVolume,
727 G4VoxelLimits pLimits,
728 const G4VolumeNosVector* pCandidates,
729 EAxis pAxis)
730 {
731 G4double motherMinExtent= kInfinity, motherMaxExtent= -kInfinity,
732 targetMinExtent= kInfinity, targetMaxExtent= -kInfinity;
733 G4VPhysicalVolume* pDaughter = nullptr;
734 G4VPVParameterisation* pParam = nullptr;
735 G4VSolid *targetSolid;
736 G4AffineTransform targetTransform;
737 G4bool replicated;
738 std::size_t nCandidates = pCandidates->size();
739 std::size_t nVol, nNode, targetVolNo;
740 G4VoxelLimits noLimits;
741
742 #ifdef G4GEOMETRY_VOXELDEBUG
743 G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
744 << " Limits = " << pLimits << G4endl
745 << " Axis = " << pAxis << G4endl
746 << " Candidates = " << nCandidates << G4endl;
747 #endif
748
749 // Compute extent of logical volume's solid along this axis
750 // NOTE: results stored locally and not preserved/reused
751 //
752 G4VSolid* outerSolid = pVolume->GetSolid();
753 const G4AffineTransform origin;
754 if( !outerSolid->CalculateExtent(pAxis, pLimits, origin,
755 motherMinExtent, motherMaxExtent) )
756 {
757 outerSolid->CalculateExtent(pAxis, noLimits, origin,
758 motherMinExtent, motherMaxExtent);
759 }
760 G4VolumeExtentVector minExtents(nCandidates,0.);
761 G4VolumeExtentVector maxExtents(nCandidates,0.);
762
763 if ( (pVolume->GetNoDaughters() == 1)
764 && (pVolume->GetDaughter(0)->IsReplicated()) )
765 {
766 // Replication data not required: only parameterisation object
767 // and volume no. List used
768 //
769 pDaughter = pVolume->GetDaughter(0);
770 pParam = pDaughter->GetParameterisation();
771 if (pParam == nullptr)
772 {
773 std::ostringstream message;
774 message << "PANIC! - Missing parameterisation." << G4endl
775 << " Replicated volume with no parameterisation object !";
776 G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
777 FatalException, message);
778 return nullptr;
779 }
780
781 // Setup daughter's transformations
782 //
783 targetTransform = G4AffineTransform(pDaughter->GetRotation(),
784 pDaughter->GetTranslation());
785 replicated = true;
786 }
787 else
788 {
789 replicated = false;
790 }
791
792 // Compute extents
793 //
794 for (nVol=0; nVol<nCandidates; ++nVol)
795 {
796 targetVolNo = (*pCandidates)[nVol];
797 if (!replicated)
798 {
799 pDaughter = pVolume->GetDaughter(targetVolNo);
800
801 // Setup daughter's transformations
802 //
803 targetTransform = G4AffineTransform(pDaughter->GetRotation(),
804 pDaughter->GetTranslation());
805 // Get underlying (and setup) solid
806 //
807 targetSolid = pDaughter->GetLogicalVolume()->GetSolid();
808 }
809 else
810 {
811 // Find solid
812 //
813 targetSolid = pParam->ComputeSolid((G4int)targetVolNo,pDaughter);
814
815 // Setup solid
816 //
817 targetSolid->ComputeDimensions(pParam,(G4int)targetVolNo,pDaughter);
818
819 // Setup transform
820 //
821 pParam->ComputeTransformation((G4int)targetVolNo,pDaughter);
822 targetTransform = G4AffineTransform(pDaughter->GetRotation(),
823 pDaughter->GetTranslation());
824 }
825 // Calculate extents
826 //
827 if(!targetSolid->CalculateExtent(pAxis, pLimits, targetTransform,
828 targetMinExtent, targetMaxExtent))
829 {
830 targetSolid->CalculateExtent(pAxis, noLimits, targetTransform,
831 targetMinExtent,targetMaxExtent);
832 }
833 minExtents[nVol] = targetMinExtent;
834 maxExtents[nVol] = targetMaxExtent;
835
836 #ifdef G4GEOMETRY_VOXELDEBUG
837 G4cout << "---------------------------------------------------" << G4endl
838 << " Volume = " << pDaughter->GetName() << G4endl
839 << " Min Extent = " << targetMinExtent << G4endl
840 << " Max Extent = " << targetMaxExtent << G4endl
841 << "---------------------------------------------------" << G4endl;
842 #endif
843
844 // Check not entirely outside mother when processing toplevel nodes
845 //
846 if ( (!pLimits.IsLimited()) && ((targetMaxExtent<=motherMinExtent)
847 ||(targetMinExtent>=motherMaxExtent)) )
848 {
849 std::ostringstream message;
850 message << "PANIC! - Overlapping daughter with mother volume." << G4endl
851 << " Daughter physical volume "
852 << pDaughter->GetName() << G4endl
853 << " is entirely outside mother logical volume "
854 << pVolume->GetName() << " !!";
855 G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0002",
856 FatalException, message);
857 }
858
859 #ifdef G4GEOMETRY_VOXELDEBUG
860 // Check for straddling volumes when debugging.
861 // If a volume is >kStraddlePercent percent over the mother
862 // boundary, print a warning.
863 //
864 if (!pLimits.IsLimited())
865 {
866 G4double width;
867 G4int kStraddlePercent = 5;
868 width = maxExtents[nVol]-minExtents[nVol];
869 if ( (((motherMinExtent-minExtents[nVol])*100/width) > kStraddlePercent)
870 ||(((maxExtents[nVol]-motherMaxExtent)*100/width) > kStraddlePercent) )
871 {
872 G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
873 << " WARNING : Daughter # " << nVol
874 << " name = " << pDaughter->GetName() << G4endl
875 << " Crosses mother boundary of logical volume, name = "
876 << pVolume->GetName() << G4endl
877 << " by more than " << kStraddlePercent
878 << "%" << G4endl;
879 }
880 }
881 #endif
882 }
883
884 // Extents of all daughters known
885
886 // Calculate minimum slice width, only including volumes inside the limits
887 //
888 G4double minWidth = kInfinity;
889 G4double currentWidth;
890 for (nVol=0; nVol<nCandidates; ++nVol)
891 {
892 // currentWidth should -always- be a positive value. Inaccurate computed extent
893 // from the solid or situations of malformed geometries (overlaps) may lead to
894 // negative values and therefore unpredictable crashes !
895 //
896 currentWidth = std::abs(maxExtents[nVol]-minExtents[nVol]);
897 if ( (currentWidth<minWidth)
898 && (maxExtents[nVol]>=pLimits.GetMinExtent(pAxis))
899 && (minExtents[nVol]<=pLimits.GetMaxExtent(pAxis)) )
900 {
901 minWidth = currentWidth;
902 }
903 }
904
905 // No. of Nodes formula - nearest integer to
906 // mother width/half min daughter width +1
907 //
908 G4double noNodesExactD = ((motherMaxExtent-motherMinExtent)*2.0/minWidth)+1.0;
909
910 // Compare with "smartless quality", i.e. the average number of slices
911 // used per contained volume.
912 //
913 G4double smartlessComputed = noNodesExactD / nCandidates;
914 G4double smartlessUser = pVolume->GetSmartless();
915 G4double smartless = (smartlessComputed <= smartlessUser)
916 ? smartlessComputed : smartlessUser;
917 G4double noNodesSmart = smartless*nCandidates;
918 auto noNodesExactI = G4int(noNodesSmart);
919 G4long noNodes = ((noNodesSmart-noNodesExactI)>=0.5)
920 ? noNodesExactI+1 : noNodesExactI;
921 if( noNodes == 0 ) { noNodes=1; }
922
923 #ifdef G4GEOMETRY_VOXELDEBUG
924 G4cout << " Smartless computed = " << smartlessComputed << G4endl
925 << " Smartless volume = " << smartlessUser
926 << " => # Smartless = " << smartless << G4endl;
927 G4cout << " Min width = " << minWidth
928 << " => # Nodes = " << noNodes << G4endl;
929 #endif
930
931 if (noNodes>kMaxVoxelNodes)
932 {
933 noNodes=kMaxVoxelNodes;
934 #ifdef G4GEOMETRY_VOXELDEBUG
935 G4cout << " Nodes Clipped to = " << kMaxVoxelNodes << G4endl;
936 #endif
937 }
938 G4double nodeWidth = (motherMaxExtent-motherMinExtent)/noNodes;
939
940 // Create G4VoxelNodes. Will Add proxies before setting fslices
941 //
942 auto* nodeList = new G4NodeVector();
943 if (nodeList == nullptr)
944 {
945 G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
946 FatalException, "NodeList allocation error.");
947 return nullptr;
948 }
949 nodeList->reserve(noNodes);
950
951 for (nNode=0; G4long(nNode)<noNodes; ++nNode)
952 {
953 G4SmartVoxelNode *pNode;
954 pNode = new G4SmartVoxelNode((G4int)nNode);
955 if (pNode == nullptr)
956 {
957 G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
958 FatalException, "Node allocation error.");
959 return nullptr;
960 }
961 nodeList->push_back(pNode);
962 }
963
964 // All nodes created (empty)
965
966 // Fill nodes: Step through extent lists
967 //
968 for (nVol=0; nVol<nCandidates; ++nVol)
969 {
970 G4long nodeNo, minContainingNode, maxContainingNode;
971 minContainingNode = (minExtents[nVol]-motherMinExtent)/nodeWidth;
972 maxContainingNode = (maxExtents[nVol]-motherMinExtent)/nodeWidth;
973
974 // Only add nodes that are inside the limits of the axis
975 //
976 if ( (maxContainingNode>=0) && (minContainingNode<noNodes) )
977 {
978 // If max extent is on max boundary => maxContainingNode=noNodes:
979 // should be one less as nodeList has noNodes entries
980 //
981 if (maxContainingNode>=noNodes)
982 {
983 maxContainingNode = noNodes-1;
984 }
985 //
986 // Protection against protruding volumes
987 //
988 if (minContainingNode<0)
989 {
990 minContainingNode = 0;
991 }
992 for (nodeNo=minContainingNode; nodeNo<=maxContainingNode; ++nodeNo)
993 {
994 (*nodeList)[nodeNo]->Insert((*pCandidates)[nVol]);
995 }
996 }
997 }
998
999 // All nodes filled
1000
1001 // Create proxy List : caller has deletion responsibility
1002 // (but we must delete nodeList *itself* - not the contents)
1003 //
1004 auto* proxyList = new G4ProxyVector();
1005 if (proxyList == nullptr)
1006 {
1007 G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
1008 FatalException, "Proxy list allocation error.");
1009 return nullptr;
1010 }
1011 proxyList->reserve(noNodes);
1012
1013 //
1014 // Fill proxy List
1015 //
1016 for (nNode=0; G4long(nNode)<noNodes; ++nNode)
1017 {
1018 // Get rid of possible excess capacity in the internal node vector
1019 //
1020 ((*nodeList)[nNode])->Shrink();
1021 auto* pProxyNode = new G4SmartVoxelProxy((*nodeList)[nNode]);
1022 if (pProxyNode == nullptr)
1023 {
1024 G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
1025 FatalException, "Proxy node allocation failed.");
1026 return nullptr;
1027 }
1028 proxyList->push_back(pProxyNode);
1029 }
1030 delete nodeList;
1031 return proxyList;
1032 }
1033
1034 // ***************************************************************************
1035 // Calculate a "quality value" for the specified vector of voxels.
1036 // The value returned should be >0 and such that the smaller the number
1037 // the higher the quality of the slice.
1038 //
1039 // Preconditions: pSlice must consist of G4SmartVoxelNodeProxies only
1040 // Process:
1041 // o Examine each node in turn, summing:
1042 // no. of non-empty nodes
1043 // no. of volumes in each node
1044 // o Calculate Quality=sigma(volumes in nod)/(no. of non-empty nodes)
1045 // if all nodes empty, return kInfinity
1046 // o Call G4Exception on finding a G4SmartVoxelHeaderProxy
1047 // ***************************************************************************
1048 //
1049 G4double G4SmartVoxelHeader::CalculateQuality(G4ProxyVector *pSlice)
1050 {
1051 G4double quality;
1052 std::size_t nNodes = pSlice->size();
1053 std::size_t noContained, maxContained=0, sumContained=0, sumNonEmptyNodes=0;
1054 G4SmartVoxelNode *node;
1055
1056 for (std::size_t i=0; i<nNodes; ++i)
1057 {
1058 if ((*pSlice)[i]->IsNode())
1059 {
1060 // Definitely a node. Add info to running totals
1061 //
1062 node = (*pSlice)[i]->GetNode();
1063 noContained = node->GetNoContained();
1064 if (noContained != 0)
1065 {
1066 ++sumNonEmptyNodes;
1067 sumContained += noContained;
1068 //
1069 // Calc maxContained for statistics
1070 //
1071 if (noContained>maxContained)
1072 {
1073 maxContained = noContained;
1074 }
1075 }
1076 }
1077 else
1078 {
1079 G4Exception("G4SmartVoxelHeader::CalculateQuality()", "GeomMgt0001",
1080 FatalException, "Not applicable to replicated volumes.");
1081 }
1082 }
1083
1084 // Calculate quality with protection against no non-empty nodes
1085 //
1086 if (sumNonEmptyNodes != 0)
1087 {
1088 quality = sumContained/sumNonEmptyNodes;
1089 }
1090 else
1091 {
1092 quality = kInfinity;
1093 }
1094
1095 #ifdef G4GEOMETRY_VOXELDEBUG
1096 G4cout << "**** G4SmartVoxelHeader::CalculateQuality" << G4endl
1097 << " Quality = " << quality << G4endl
1098 << " Nodes = " << nNodes
1099 << " of which " << sumNonEmptyNodes << " non empty" << G4endl
1100 << " Max Contained = " << maxContained << G4endl;
1101 #endif
1102
1103 return quality;
1104 }
1105
1106 // ***************************************************************************
1107 // Examined each contained node, refines (creates a replacement additional
1108 // dimension of voxels) when there is more than one voxel in the slice.
1109 // Does not refine further if already limited in two dimensions (=> this
1110 // is the third level of limits)
1111 //
1112 // Preconditions: slices (nodes) have been built.
1113 // ***************************************************************************
1114 //
1115 void G4SmartVoxelHeader::RefineNodes(G4LogicalVolume* pVolume,
1116 G4VoxelLimits pLimits)
1117 {
1118 std::size_t refinedDepth=0, minVolumes;
1119 std::size_t maxNode = fslices.size();
1120
1121 if (pLimits.IsXLimited())
1122 {
1123 ++refinedDepth;
1124 }
1125 if (pLimits.IsYLimited())
1126 {
1127 ++refinedDepth;
1128 }
1129 if (pLimits.IsZLimited())
1130 {
1131 ++refinedDepth;
1132 }
1133
1134 // Calculate minimum number of volumes necessary to refine
1135 //
1136 switch (refinedDepth)
1137 {
1138 case 0:
1139 minVolumes=kMinVoxelVolumesLevel2;
1140 break;
1141 case 1:
1142 minVolumes=kMinVoxelVolumesLevel3;
1143 break;
1144 default:
1145 minVolumes=10000; // catch refinedDepth=3 and errors
1146 break;
1147 }
1148
1149 if (refinedDepth<2)
1150 {
1151 std::size_t targetNo, noContainedDaughters, minNo, maxNo, replaceNo, i;
1152 G4double sliceWidth = (fmaxExtent-fminExtent)/maxNode;
1153 G4VoxelLimits newLimits;
1154 G4SmartVoxelNode* targetNode;
1155 G4SmartVoxelProxy* targetNodeProxy;
1156 G4SmartVoxelHeader* replaceHeader;
1157 G4SmartVoxelProxy* replaceHeaderProxy;
1158 G4VolumeNosVector* targetList;
1159 G4SmartVoxelProxy* lastProxy;
1160
1161 for (targetNo=0; targetNo<maxNode; ++targetNo)
1162 {
1163 // Assume all slices are nodes (see preconditions)
1164 //
1165 targetNodeProxy = fslices[targetNo];
1166 targetNode = targetNodeProxy->GetNode();
1167
1168 if (targetNode->GetNoContained() >= minVolumes)
1169 {
1170 noContainedDaughters = targetNode->GetNoContained();
1171 targetList = new G4VolumeNosVector();
1172 if (targetList == nullptr)
1173 {
1174 G4Exception("G4SmartVoxelHeader::RefineNodes()",
1175 "GeomMgt0003", FatalException,
1176 "Target volume node list allocation error.");
1177 return;
1178 }
1179 targetList->reserve(noContainedDaughters);
1180 for (i=0; i<noContainedDaughters; ++i)
1181 {
1182 targetList->push_back(targetNode->GetVolume((G4int)i));
1183 }
1184 minNo = targetNode->GetMinEquivalentSliceNo();
1185 maxNo = targetNode->GetMaxEquivalentSliceNo();
1186
1187 #ifdef G4GEOMETRY_VOXELDEBUG
1188 G4cout << "**** G4SmartVoxelHeader::RefineNodes" << G4endl
1189 << " Refining nodes " << minNo
1190 << " - " << maxNo << " inclusive" << G4endl;
1191 #endif
1192 if (minNo > maxNo) // Delete node and list to be replaced
1193 { // and avoid further action ...
1194 delete targetNode;
1195 delete targetList;
1196 return;
1197 }
1198
1199 // Delete node proxies at start of collected sets of nodes/headers
1200 //
1201 lastProxy=nullptr;
1202 for (replaceNo=minNo; replaceNo<=maxNo; ++replaceNo)
1203 {
1204 if (lastProxy != fslices[replaceNo])
1205 {
1206 lastProxy=fslices[replaceNo];
1207 delete lastProxy;
1208 }
1209 }
1210 // Delete node to be replaced
1211 //
1212 delete targetNode;
1213
1214 // Create new headers + proxies and replace in fslices
1215 //
1216 newLimits = pLimits;
1217 newLimits.AddLimit(faxis,fminExtent+sliceWidth*minNo,
1218 fminExtent+sliceWidth*(maxNo+1));
1219 replaceHeader = new G4SmartVoxelHeader(pVolume,newLimits,
1220 targetList,(G4int)replaceNo);
1221 if (replaceHeader == nullptr)
1222 {
1223 G4Exception("G4SmartVoxelHeader::RefineNodes()", "GeomMgt0003",
1224 FatalException, "Refined VoxelHeader allocation error.");
1225 return;
1226 }
1227 replaceHeader->SetMinEquivalentSliceNo((G4int)minNo);
1228 replaceHeader->SetMaxEquivalentSliceNo((G4int)maxNo);
1229 replaceHeaderProxy = new G4SmartVoxelProxy(replaceHeader);
1230 if (replaceHeaderProxy == nullptr)
1231 {
1232 G4Exception("G4SmartVoxelHeader::RefineNodes()", "GeomMgt0003",
1233 FatalException, "Refined VoxelProxy allocation error.");
1234 return;
1235 }
1236 for (replaceNo=minNo; replaceNo<=maxNo; ++replaceNo)
1237 {
1238 fslices[replaceNo] = replaceHeaderProxy;
1239 }
1240 // Finished replacing current `equivalent' group
1241 //
1242 delete targetList;
1243 targetNo=maxNo;
1244 }
1245 }
1246 }
1247 }
1248
1249 // ***************************************************************************
1250 // Returns true if all slices have equal contents.
1251 // Preconditions: all equal slices have been collected.
1252 // Procedure:
1253 // o checks all slice proxy pointers are equal
1254 // o returns true if only one slice or all slice proxies pointers equal.
1255 // ***************************************************************************
1256 //
1257 G4bool G4SmartVoxelHeader::AllSlicesEqual() const
1258 {
1259 std::size_t noSlices = fslices.size();
1260 G4SmartVoxelProxy* refProxy;
1261
1262 if (noSlices>1)
1263 {
1264 refProxy=fslices[0];
1265 for (std::size_t i=1; i<noSlices; ++i)
1266 {
1267 if (refProxy!=fslices[i])
1268 {
1269 return false;
1270 }
1271 }
1272 }
1273 return true;
1274 }
1275
1276 // ***************************************************************************
1277 // Streaming operator for debugging.
1278 // ***************************************************************************
1279 //
1280 std::ostream& operator << (std::ostream& os, const G4SmartVoxelHeader& h)
1281 {
1282 os << "Axis = " << G4int(h.faxis) << G4endl;
1283 G4SmartVoxelProxy *collectNode=nullptr, *collectHead=nullptr;
1284 std::size_t collectNodeNo = 0;
1285 std::size_t collectHeadNo = 0;
1286 std::size_t i, j;
1287 G4bool haveHeaders = false;
1288
1289 for (i=0; i<h.fslices.size(); ++i)
1290 {
1291 os << "Slice #" << i << " = ";
1292 if (h.fslices[i]->IsNode())
1293 {
1294 if (h.fslices[i]!=collectNode)
1295 {
1296 os << "{";
1297 for (std::size_t k=0; k<h.fslices[i]->GetNode()->GetNoContained(); ++k)
1298 {
1299 os << " " << h.fslices[i]->GetNode()->GetVolume((G4int)k);
1300 }
1301 os << " }" << G4endl;
1302 collectNode = h.fslices[i];
1303 collectNodeNo = i;
1304 }
1305 else
1306 {
1307 os << "As slice #" << collectNodeNo << G4endl;
1308 }
1309 }
1310 else
1311 {
1312 haveHeaders=true;
1313 if (h.fslices[i] != collectHead)
1314 {
1315 os << "Header" << G4endl;
1316 collectHead = h.fslices[i];
1317 collectHeadNo = i;
1318 }
1319 else
1320 {
1321 os << "As slice #" << collectHeadNo << G4endl;
1322 }
1323 }
1324 }
1325
1326 if (haveHeaders)
1327 {
1328 collectHead=nullptr;
1329 for (j=0; j<h.fslices.size(); ++j)
1330 {
1331 if (h.fslices[j]->IsHeader())
1332 {
1333 os << "Header at Slice #" << j << " = ";
1334 if (h.fslices[j] != collectHead)
1335 {
1336 os << G4endl
1337 << (*(h.fslices[j]->GetHeader()));
1338 collectHead = h.fslices[j];
1339 collectHeadNo = j;
1340 }
1341 else
1342 {
1343 os << "As slice #" << collectHeadNo << G4endl;
1344 }
1345 }
1346 }
1347 }
1348 return os;
1349 }
1350