source: trunk/source/geometry/solids/specific/test/testG4Paraboloid.cc

Last change on this file was 1316, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 13.3 KB
Line 
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//
27#include <assert.h>
28
29#include "G4GeometryTolerance.hh"
30#include "G4Paraboloid.hh"
31#include "G4Polyhedron.hh"
32#include "Randomize.hh"
33int main()
34{
35 G4ThreeVector testPoint;
36
37 G4Paraboloid paraboloid1("Paraboloid1", 1., 0., 3.);
38 G4Paraboloid paraboloid2("Paraboloid2", .01, 2., 3.);
39
40 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
41
42 // Test Inside function.
43 assert(paraboloid1.Inside(G4ThreeVector(2.598076211, 1.5, 1.)) == kSurface);
44 assert(paraboloid1.Inside(G4ThreeVector((2 * kCarTolerance + 3.) * std::sqrt(3.) / 2., (2. * kCarTolerance + 3.) / 2., 1.)) == kOutside);
45 assert(paraboloid1.Inside(G4ThreeVector((- 2 * kCarTolerance + 3.) * std::sqrt(3.) / 2., (- 2. * kCarTolerance + 3.) / 2., 1.)) == kSurface);
46
47 assert(paraboloid1.Inside(G4ThreeVector(std::sqrt(4.5) / std::sqrt(2.), std::sqrt(4.5) / std::sqrt(2.), 0.)) == kSurface);
48 assert(paraboloid1.Inside(G4ThreeVector(( 2 * kCarTolerance + std::sqrt(4.5)) / std::sqrt(2.), ( 2. * kCarTolerance + std::sqrt(4.5)) / std::sqrt(2.), 0.)) == kOutside);
49 assert(paraboloid1.Inside(G4ThreeVector((- 2 * kCarTolerance + std::sqrt(4.5)) / std::sqrt(2.), (- 2. * kCarTolerance + std::sqrt(4.5)) / std::sqrt(2.), 0.)) == kInside);
50
51 assert(paraboloid1.Inside(G4ThreeVector(0., 0., -1.)) == kSurface);
52 assert(paraboloid1.Inside(G4ThreeVector(0., 0., 1.)) == kSurface);
53 assert(paraboloid1.Inside(G4ThreeVector(0., 0., - 1. - 2 * kCarTolerance)) == kOutside);
54 assert(paraboloid1.Inside(G4ThreeVector(0., 0., 1. + 2 * kCarTolerance)) == kOutside);
55 assert(paraboloid1.Inside(G4ThreeVector(0., 0., - 1. + 2 * kCarTolerance)) == kInside);
56 assert(paraboloid1.Inside(G4ThreeVector(0., 0., 1. - 2 * kCarTolerance)) == kInside);
57 assert(paraboloid1.Inside(G4ThreeVector((- 2 * kCarTolerance) / std::sqrt(2.), ( 2. * kCarTolerance ) / std::sqrt(2.), -1.)) == kOutside);
58
59 assert(paraboloid2.Inside(G4ThreeVector(2.598076211, 1.5, .01)) == kSurface);
60 assert(paraboloid2.Inside(G4ThreeVector((2 * kCarTolerance + 3.) * std::sqrt(3.) / 2., (2. * kCarTolerance + 3.) / 2., .01)) == kOutside);
61 assert(paraboloid2.Inside(G4ThreeVector((- 2 * kCarTolerance + 3.) * std::sqrt(3.) / 2., (- 2. * kCarTolerance + 3.) / 2., .01)) == kSurface);
62
63 assert(paraboloid2.Inside(G4ThreeVector(std::sqrt(6.5) / std::sqrt(2.), std::sqrt(6.5) / std::sqrt(2.), 0.)) == kSurface);
64 assert(paraboloid2.Inside(G4ThreeVector(( 2. * kCarTolerance + std::sqrt(6.5)) / std::sqrt(2.), ( 2. * kCarTolerance + std::sqrt(6.5)) / std::sqrt(2.), 0.)) == kOutside);
65 assert(paraboloid2.Inside(G4ThreeVector((- 2. * kCarTolerance + std::sqrt(6.5)) / std::sqrt(2.), (- 2. * kCarTolerance + std::sqrt(6.5)) / std::sqrt(2.), 0.)) == kInside);
66
67 assert(paraboloid2.Inside(G4ThreeVector(0., 0., -.01)) == kSurface);
68 assert(paraboloid2.Inside(G4ThreeVector(0., 0., .01)) == kSurface);
69 assert(paraboloid2.Inside(G4ThreeVector(0., 0., - .01 - 2 * kCarTolerance)) == kOutside);
70 assert(paraboloid2.Inside(G4ThreeVector(0., 0., .01 + 2 * kCarTolerance)) == kOutside);
71 assert(paraboloid2.Inside(G4ThreeVector(0., 0., - .01 + 2 * kCarTolerance)) == kInside);
72 assert(paraboloid2.Inside(G4ThreeVector(0., 0., .01 - 2 * kCarTolerance)) == kInside);
73 assert(paraboloid2.Inside(G4ThreeVector((- 2. - 2 * kCarTolerance) / std::sqrt(2.), (2. + 2. * kCarTolerance ) / std::sqrt(2.), -.01)) == kOutside);
74 assert(paraboloid2.Inside(G4ThreeVector((- 2. + 2 * kCarTolerance) / std::sqrt(2.), (2. - 2. * kCarTolerance ) / std::sqrt(2.), -.01)) == kSurface);
75 assert(paraboloid2.Inside(G4ThreeVector((- 2.) / std::sqrt(2.), (2.) / std::sqrt(2.), -.01)) == kSurface);
76
77 // Test DistanceToIn(p, v) function.
78 assert(std::fabs(paraboloid1.DistanceToIn( G4ThreeVector(3.,2.,10.),
79 G4ThreeVector(-0.099014754299999993558678568206233, -0.099014754299999993558678568206233, -0.99014754279999994679428709787317))
80 - 9.0895544461477406628091557649896) < kCarTolerance);
81 assert(std::fabs(paraboloid1.DistanceToIn(G4ThreeVector(1., 0., 1.), G4ThreeVector(0., 0., -1.))) < kCarTolerance);
82 assert(std::fabs(paraboloid1.DistanceToIn(G4ThreeVector(0., 0., -2.), G4ThreeVector(0., 0., 1.)) - 1) < kCarTolerance);
83 assert(std::fabs(paraboloid1.DistanceToIn(G4ThreeVector(-2., 0., -1.), G4ThreeVector(1/std::sqrt(2.), 0., 1/std::sqrt(2.))) - 1/ std::sqrt(2.)) < kCarTolerance);
84 assert(std::fabs(paraboloid1.DistanceToIn(G4ThreeVector(1., 0., 1.), G4ThreeVector(1., 0., 0.)) ) >= kInfinity);
85 assert(std::fabs(paraboloid1.DistanceToIn(G4ThreeVector(1., 0., 1.), G4ThreeVector(1., 0., -0.0000001)) ) <= kCarTolerance);
86
87 assert(std::fabs(paraboloid2.DistanceToIn( G4ThreeVector(3.,2.,10.),
88 G4ThreeVector(-0.099014754299999993558678568206233, -0.099014754299999993558678568206233, -0.99014754279999994679428709787317))
89 - 10.0894054352) < kCarTolerance);
90 assert(std::fabs(paraboloid2.DistanceToIn(G4ThreeVector(1., 0., .01), G4ThreeVector(0., 0., -1.))) < kCarTolerance);
91 assert(std::fabs(paraboloid2.DistanceToIn(G4ThreeVector(0., 0., -1.), G4ThreeVector(0., 0., 1.)) - 0.99) < kCarTolerance);
92 assert(std::fabs(paraboloid2.DistanceToIn(G4ThreeVector(-2.2, 0., -.01), G4ThreeVector(1/std::sqrt(2.), 0., 1/std::sqrt(2.))) - 0.004669633692) < kCarTolerance);
93 assert(std::fabs(paraboloid2.DistanceToIn(G4ThreeVector(1., 0., .01), G4ThreeVector(1., 0., 0.)) ) >= kInfinity);
94 assert(std::fabs(paraboloid2.DistanceToIn(G4ThreeVector(1., 0., .01), G4ThreeVector(1., 0., -0.0000001)) ) <= kCarTolerance);
95
96 // Test DistanceToOut(p, v, ...) function.
97 assert(std::fabs(paraboloid1.DistanceToOut(G4ThreeVector(1., 0., 1.), G4ThreeVector(0., 0., -1.), false, NULL, NULL) - 1.7777777777777776) < kCarTolerance);
98 assert(std::fabs(paraboloid1.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(0., 0., -1.), false, NULL, NULL) - 1) < kCarTolerance);
99 assert(std::fabs(paraboloid1.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.), false, NULL, NULL) - 2.1213203435) < kCarTolerance);
100 assert(std::fabs(paraboloid1.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(1. / std::sqrt(2.), 0., 1. / std::sqrt(2.)), false, NULL, NULL) - 1.4142135623) < kCarTolerance);
101 assert(std::fabs(paraboloid1.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(1. / std::sqrt(2.), 0., -1. / std::sqrt(2.)), false, NULL, NULL) - 1.1912334058) < kCarTolerance);
102 assert(std::fabs(paraboloid1.DistanceToOut(G4ThreeVector(1., 0., 1.), G4ThreeVector(1., 0., 0.), false, NULL, NULL) - 2.) < kCarTolerance);
103 assert(std::fabs(paraboloid1.DistanceToOut(G4ThreeVector(1., 0., 1.), G4ThreeVector(1., 0., 0.00000001), false, NULL, NULL)) < kCarTolerance);
104
105 assert(std::fabs(paraboloid2.DistanceToOut(G4ThreeVector(2.5, 0., .01), G4ThreeVector(0., 0., -1.), false, NULL, NULL) - 0.010999999999999999) < kCarTolerance);
106 assert(std::fabs(paraboloid2.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(0., 0., -1.), false, NULL, NULL) - 0.01) < kCarTolerance);
107 assert(std::fabs(paraboloid2.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.), false, NULL, NULL) - 2.5495097567) < kCarTolerance);
108 assert(std::fabs(paraboloid2.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(1. / std::sqrt(2.), 0., 1. / std::sqrt(2.)), false, NULL, NULL) - 0.01414213562) < kCarTolerance);
109 assert(std::fabs(paraboloid2.DistanceToOut(G4ThreeVector(0., 0., 0.), G4ThreeVector(1. / std::sqrt(2.), 0., -1. / std::sqrt(2.)), false, NULL, NULL) - 0.01414213562) < kCarTolerance);
110 assert(std::fabs(paraboloid2.DistanceToOut(G4ThreeVector(1., 0., -.01), G4ThreeVector(1., 0., 0.), false, NULL, NULL) - 1.) < kCarTolerance);
111 assert(std::fabs(paraboloid2.DistanceToOut(G4ThreeVector(1., 0., .01), G4ThreeVector(1., 0., 0.00000001), false, NULL, NULL)) < kCarTolerance);
112
113 // Test volume function.
114 assert(std::fabs(paraboloid1.GetCubicVolume() - 28.274334) < 1e-6);
115 assert(std::fabs(paraboloid2.GetCubicVolume() - 0.408407) < 1e-6);
116
117 // Test area function.
118 // G4cout<<paraboloid1.GetPolyhedron()->GetSurfaceArea()<<G4endl;
119 assert(std::fabs(paraboloid1.GetSurfaceArea() - 66.758844) < 1e-6);
120 assert(std::fabs(paraboloid2.GetSurfaceArea() - 56.551935) < 1e-6);
121
122 // Test SurfaceNormal(p) function.
123 assert((paraboloid1.SurfaceNormal(G4ThreeVector(1., 2., 1.)) - G4ThreeVector(0., 0., 1.)).r() < kCarTolerance);
124 assert((paraboloid1.SurfaceNormal(G4ThreeVector(1.732050808, 2.449489743, 1.)) - G4ThreeVector(0.40824829,0.57735027,0.70710678)).r() < 1e-5);
125 assert((paraboloid1.SurfaceNormal(G4ThreeVector(1.573213272, 1.573213272, 0.1)) - G4ThreeVector(0.49718308,0.49718308,-0.71106819)).r() < 1e-5);
126 assert((paraboloid1.SurfaceNormal(G4ThreeVector(0., 0., -1.)) - G4ThreeVector(0., 0., -1.)).r() < kCarTolerance);
127
128 assert((paraboloid2.SurfaceNormal(G4ThreeVector(1., 2., 1.)) - G4ThreeVector(0., 0., 1.)).r() < kCarTolerance);
129 assert((paraboloid2.SurfaceNormal(G4ThreeVector(1.732050808, 2.449489743, 0.01)) - G4ThreeVector(0.40824829,0.57735027,0.70710678)).r() < 1e-5);
130 assert((paraboloid2.SurfaceNormal(G4ThreeVector(1.658312395, 2., 0.001)) - G4ThreeVector(0.013263635,0.015996545,-0.99978407)).r() < 1e-5);
131 assert((paraboloid2.SurfaceNormal(G4ThreeVector(0., 0., -.01)) - G4ThreeVector(0., 0., -1.)).r() < kCarTolerance);
132
133
134 // Add Test Distance To Out for User Problem with Substraction Of Paraboloid :: Problem Report N1015
135 G4Paraboloid paraboloid3("Paraboloid1", 72., 0., 240.);
136 // G4Paraboloid paraboloid3("Paraboloid1", 72., 210., 240.);//This test is Ok, Intersection with DZ working correctly
137 G4double distOut1;
138 distOut1=paraboloid3.DistanceToOut(G4ThreeVector(200., 0., 72.), G4ThreeVector(0.0,0.,-1.), false, NULL, NULL);
139 assert((distOut1-44.0)<1e-5);
140 //G4cout<<"Test3::distOut="<<distOut1<<G4endl;
141
142 distOut1=paraboloid3.DistanceToOut(G4ThreeVector(200., 0., 72.), G4ThreeVector(0.000000000001,0.,-1.), false, NULL, NULL);
143 assert((distOut1-44.0)<1e-5);
144 //G4cout<<"Test3::distOut="<<distOut1<<G4endl;
145 // Add Rotation On the Surface around "Problem" Point(200.,0.,72.)
146
147 G4double tolA=0.25e-7;
148 G4ThreeVector In (200.,0.,72.);
149 EInside in;
150 G4ThreeVector Dir,vDir,pSurf;
151 Dir=G4ThreeVector(10*tolA,0.,-1.);
152 vDir=Dir.unit();
153 in=paraboloid3.Inside(In);
154 if(in!=kSurface)G4cout<<"Problem ::Point p"<<In<<" is NOT On Surface"<<G4endl;
155
156 for(G4int i=0; i<100; i++){
157
158 G4double distOut=paraboloid3.DistanceToOut(In,vDir);
159 G4cout.precision(16);
160 pSurf=In+vDir*distOut;
161 in=paraboloid3.Inside(pSurf);
162 if(in!=kSurface) G4cout<<"Problem :: ps is Not on the Surface DistToOut="<<distOut<<" dir="<<vDir<<G4endl;
163 Dir=Dir-G4ThreeVector(tolA,0,0);
164 vDir=Dir.unit();
165
166 }
167 // Rotation from the Lower DZ plane
168 Dir=G4ThreeVector(10*tolA,0.,1.);
169 vDir=Dir.unit();
170 In= G4ThreeVector (0.,0.,-72.);
171 for(G4int i=0; i<100; i++){
172
173 G4double distOut=paraboloid3.DistanceToOut(In,vDir);
174 G4cout.precision(16);
175 pSurf=In+vDir*distOut;
176 in=paraboloid3.Inside(pSurf);
177 if(in!=kSurface) G4cout<<"Problem :: ps is Not on the Surface DistToOut="<<distOut<<" dir="<<vDir<<G4endl;
178 Dir=Dir-G4ThreeVector(tolA,0,0);
179 vDir=Dir.unit();
180 }
181 // Paraboloid with R1!=0 and with Point situated On both Surfaces : Z and Parabolic
182 G4Paraboloid paraboloid5("Paraboloid1", 10., 0., 100.);
183 distOut1=paraboloid5.DistanceToOut(G4ThreeVector(100., 0., 10.), G4ThreeVector(0.0,0.,-1.), false, NULL, NULL);
184 assert((distOut1-0.0)<1e-5);
185
186 Dir=G4ThreeVector(-0.9,0.,-0.1);
187 vDir=Dir.unit();
188 distOut1=paraboloid5.DistanceToOut(G4ThreeVector(100., 0., 10.), vDir, false, NULL, NULL);
189 pSurf=G4ThreeVector(100.0,0.,10.)+vDir*distOut1;
190 in=paraboloid5.Inside(pSurf);
191 if(in!=kSurface)G4cout<<"Problem :: ps is Not on the Surface "<<pSurf<<" vDir="<<vDir<<G4endl;
192
193 return 0;
194}
Note: See TracBrowser for help on using the repository browser.