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

Last change on this file since 1350 was 1316, checked in by garnier, 14 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.