u4/tests/U4TreeTessellatedTest.cc

0001 /**
0002  * Verifies U4Tree handling of repeated, nested tessellated solids.
0003  *
0004  * The test builds a Boolean solid containing a tessellated tetrahedron at
0005  * runtime and places its logical volume twice in a simple world. It checks
0006  * that U4Tree selects the enclosing Boolean solid for triangulation, keeps
0007  * both placements out of analytic repeat factors, and records both as global
0008  * triangulated nodes. This is a CPU-side geometry-conversion test; it does not
0009  * run GPU intersection.
0010  */
0011
0012 #include <algorithm>
0013 #include <cstdlib>
0014 #include <iostream>
0015 #include <vector>
0016
0017 #include "G4Box.hh"
0018 #include "G4LogicalVolume.hh"
0019 #include "G4PVPlacement.hh"
0020 #include "G4TessellatedSolid.hh"
0021 #include "G4TriangularFacet.hh"
0022 #include "G4UnionSolid.hh"
0023
0024 #include "OPTICKS_LOG.hh"
0025 #include "U4Material.hh"
0026 #include "U4Tree.h"
0027
0028 #include "NPFold.h"
0029 #include "stree.h"
0030
0031 namespace
0032 {
0033 bool addFacet(G4TessellatedSolid& solid, const G4ThreeVector& a, const G4ThreeVector& b, const G4ThreeVector& c)
0034 {
0035     return solid.AddFacet(new G4TriangularFacet(a, b, c, ABSOLUTE));
0036 }
0037
0038 G4TessellatedSolid* makeTetrahedron()
0039 {
0040     G4ThreeVector p0(0., 0., 0.);
0041     G4ThreeVector p1(4., 0., 0.);
0042     G4ThreeVector p2(0., 4., 0.);
0043     G4ThreeVector p3(0., 0., 4.);
0044
0045     G4TessellatedSolid* solid = new G4TessellatedSolid("NestedTetrahedron");
0046     bool                facetsAdded =
0047         addFacet(*solid, p0, p2, p1) &&
0048         addFacet(*solid, p0, p1, p3) &&
0049         addFacet(*solid, p0, p3, p2) &&
0050         addFacet(*solid, p1, p2, p3);
0051
0052     if (!facetsAdded)
0053         return nullptr;
0054     solid->SetSolidClosed(true);
0055     return solid;
0056 }
0057
0058 int findLvid(const stree& st, const char* name)
0059 {
0060     std::vector<std::string>::const_iterator it = std::find(st.soname.begin(), st.soname.end(), name);
0061     return it == st.soname.end() ? -1 : int(std::distance(st.soname.begin(), it));
0062 }
0063 } // namespace
0064
0065 int main(int argc, char** argv)
0066 {
0067     OPTICKS_LOG(argc, argv);
0068
0069     G4TessellatedSolid* tessellated = makeTetrahedron();
0070     if (tessellated == nullptr)
0071     {
0072         std::cerr << "failed to construct tessellated test geometry" << std::endl;
0073         return EXIT_FAILURE;
0074     }
0075
0076     G4Material* material = U4Material::Get(U4Material::VACUUM);
0077     if (material == nullptr)
0078     {
0079         std::cerr << "failed to construct vacuum material" << std::endl;
0080         return EXIT_FAILURE;
0081     }
0082
0083     G4Box*           component = new G4Box("BooleanBox", 3., 3., 3.);
0084     G4UnionSolid*    nested = new G4UnionSolid("NestedTessellatedUnion", component, tessellated);
0085     G4LogicalVolume* nestedLV = new G4LogicalVolume(nested, material, "NestedTessellatedLV");
0086
0087     G4Box*           worldSolid = new G4Box("WorldBox", 100., 100., 100.);
0088     G4LogicalVolume* worldLV = new G4LogicalVolume(worldSolid, material, "WorldLV");
0089
0090     new G4PVPlacement(nullptr, G4ThreeVector(-20., 0., 0.), nestedLV, "NestedPV0", worldLV, false, 0);
0091     new G4PVPlacement(nullptr, G4ThreeVector(20., 0., 0.), nestedLV, "NestedPV1", worldLV, false, 1);
0092     G4VPhysicalVolume* world = new G4PVPlacement(
0093         nullptr, G4ThreeVector(), worldLV, "WorldPV", nullptr, false, 0);
0094
0095     stree st;
0096     st.FREQ_CUT = 2;
0097     U4Tree* tree = U4Tree::Create(&st, world);
0098     if (tree == nullptr)
0099     {
0100         std::cerr << "U4Tree creation failed" << std::endl;
0101         return EXIT_FAILURE;
0102     }
0103
0104     int              lvid = findLvid(st, "NestedTessellatedUnion");
0105     std::vector<int> nodes;
0106     if (lvid > -1)
0107         st.find_lvid_nodes(nodes, lvid, 'N');
0108
0109     bool nodesAreGlobal = nodes.size() == 2;
0110     for (unsigned i = 0; i < nodes.size(); ++i)
0111         nodesAreGlobal = nodesAreGlobal && st.nds[nodes[i]].repeat_index == 0;
0112
0113     int triangulatedPlacements = 0;
0114     for (unsigned i = 0; i < st.tri.size(); ++i)
0115         if (st.tri[i].lvid == lvid)
0116             triangulatedPlacements += 1;
0117
0118     const NPFold* mesh = lvid > -1 ? st.mesh->get_subfold(st.soname[lvid].c_str()) : nullptr;
0119     bool          passed =
0120         lvid > -1 &&
0121         st.is_force_triangulate(lvid) &&
0122         st.get_num_factor() == 0 &&
0123         nodesAreGlobal &&
0124         triangulatedPlacements == 2 &&
0125         mesh != nullptr &&
0126         mesh->get_meta<int>("lvid", -1) == lvid;
0127
0128     if (!passed)
0129     {
0130         std::cerr
0131             << "nested tessellated solid was not retained as global triangulated geometry"
0132             << " lvid " << lvid
0133             << " factors " << st.get_num_factor()
0134             << " nodes " << nodes.size()
0135             << " triangulated placements " << triangulatedPlacements
0136             << std::endl;
0137         return EXIT_FAILURE;
0138     }
0139
0140     std::cout << "repeated nested tessellated solid remained global and triangulated" << std::endl;
0141     return EXIT_SUCCESS;
0142 }