source: trunk/source/geometry/solids/BREPS/src/G4Surface.cc@ 970

Last change on this file since 970 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 12.4 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// $Id: G4Surface.cc,v 1.17 2007/07/16 08:06:55 gcosmo Exp $
28// GEANT4 tag $Name: HEAD $
29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4Surface.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4Surface.hh"
38#include "G4CompositeCurve.hh"
39#include "G4GeometryTolerance.hh"
40
41G4Surface::G4Surface()
42 : FLT_MAXX(kInfinity), FLT_EPSILO(0.0001)
43{
44 AdvancedFace=0;
45 active = 1;
46 distance = 1.0e20;
47 Type = 0;
48 bbox = 0;
49 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
50}
51
52G4Surface::~G4Surface()
53{
54}
55
56G4int G4Surface::operator==( const G4Surface& s )
57{
58 return origin == s.origin;
59}
60
61G4String G4Surface::GetEntityType() const
62{
63 return G4String("Surface");
64}
65
66const char* G4Surface::Name() const
67{
68 return "G4Surface";
69}
70
71G4int G4Surface::MyType() const
72{
73 return Type;
74}
75
76void G4Surface::InitBounded()
77{
78}
79
80G4double G4Surface::GetUHit() const
81{
82 return uhit;
83}
84
85G4double G4Surface::GetVHit() const
86{
87 return vhit;
88}
89
90//void G4Surface::read_surface(fstream& tmp){;}
91
92G4Point3D G4Surface::Evaluation(const G4Ray&)
93{
94 return closest_hit;
95}
96
97G4int G4Surface::Evaluate(const G4Ray&)
98{
99 return 0;
100}
101
102void G4Surface::Reset()
103{
104 Intersected = 0;
105 active = 1;
106 distance = kInfinity;
107}
108
109void G4Surface::SetBoundaries(G4CurveVector* boundaries)
110{
111 surfaceBoundary.Init(*boundaries);
112 InitBounded();
113}
114
115void G4Surface::CalcBBox()
116{
117 // Finds the bounds of the surface iow
118 // calculates the bounds for a bounding box
119 // to the surface. The bounding box is used
120 // for a preliminary check of intersection.
121
122 bbox = new G4BoundingBox3D(surfaceBoundary.BBox().GetBoxMin(),
123 surfaceBoundary.BBox().GetBoxMax());
124 // old implementation
125 // G4Point3d BoundaryMax = OuterBoundary->GetBoundsMax();
126 // G4Point3d BoundaryMin = OuterBoundary->GetBoundsMin();
127 // bbox = new G4BoundingBox( BoundaryMin, BoundaryMax);
128 // return;
129}
130
131G4Vector3D G4Surface::Normal( const G4Vector3D& ) const
132{ // return the Normal unit vector to a Surface at the point p on
133 // (or nearly on) the Surface.
134 // The default is not well defined, so return ( 0, 0, 0 ).
135 return G4Vector3D( 0.0, 0.0, 0.0 );
136}
137
138
139G4int G4Surface::Intersect(const G4Ray&)
140{
141 G4int Result = 0;
142
143 G4Exception("G4Surface::Intersect()", "NotImplemented",
144 FatalException, "Sorry, not yet implemented.");
145
146#ifdef NEW_IMPLEMENTATION
147 // get the intersection
148 // Result = Intersect(rayref);
149
150 // Check that the point is within the polyline
151 // Get Normal at Hitpoint
152 const G4Vector3D& Vec = Normal(closest_hit);
153 G4Ray Normal(closest_hit, Vec);
154
155 // Project points & Hit
156 // OuterBoundary->ProjectBoundaryTo2D(Normal.GetPlane(1),
157 // Normal.GetPlane(2), 0);
158
159
160
161
162 G4Point3D Hit = closest_hit.Project(Normal.GetPlane(1),
163 Normal.GetPlane(2) );
164 // Check point in polygon
165 // Result = OuterBoundary->Inside(Hit, rayref);
166
167#endif
168 return Result;
169
170}
171
172
173G4double G4Surface::ClosestDistanceToPoint(const G4Point3D& Pt)
174{
175 // in fact, a squared distance is returned
176
177 // a bit suspicious, this function
178 // the distance is almost always an overestimate
179 G4double pointDistance= kInfinity;
180 G4double tmpDistance;
181 const G4CurveVector& bounds= surfaceBoundary.GetBounds();
182
183 G4int entr = bounds.size();
184
185 for (G4int i=0; i<entr; i++)
186 {
187 G4Curve* c= bounds[i];
188
189 if (c->GetEntityType() == "G4CompositeCurve")
190 {
191 G4CompositeCurve* cc= (G4CompositeCurve*)c;
192 const G4CurveVector& segments= cc->GetSegments();
193 for (size_t i=0; i<segments.size(); i++)
194 {
195 G4Curve* ccc= segments[i];
196 tmpDistance= (G4Point3D(Pt.x(), Pt.y(), Pt.z())-ccc->GetEnd()).mag2();
197 if (pointDistance > tmpDistance)
198 {
199 pointDistance= tmpDistance;
200 }
201 }
202
203 }
204 else
205 {
206 tmpDistance= (G4Point3D(Pt.x(), Pt.y(), Pt.z())-c->GetEnd()).mag2();
207 if (pointDistance > tmpDistance)
208 {
209 pointDistance= tmpDistance;
210 }
211 }
212 }
213
214 // L. Broglia
215 // Be carreful ! pointdistance is the squared distance
216 return std::sqrt(pointDistance);
217
218 // G4double PointDistance=kInfinity;
219 // G4double TmpDistance=0;
220 // PointDistance = OuterBoundary->ClosestDistanceToPoint(Pt);
221 // TmpDistance =0;
222 // for(G4int a=0;a<NumberOfInnerBoundaries;a++)
223 // {
224 // TmpDistance = InnerBoundary[a]->ClosestDistanceToPoint(Pt);
225 // if(PointDistance > TmpDistance) PointDistance = TmpDistance;
226 // }
227 // return PointDistance;
228
229 //G4double G4Boundary::ClosestDistanceToPoint(const G4ThreeVec& Pt)
230 //{
231 // G4double PointDistance = kInfinity;
232 // G4double TmpDistance = 0;
233 // for(G4int a =0; a < NumberOfPoints;a++)
234 // {
235 // G4Point3d& Pt2 = GetPoint(a);
236 // TmpDistance = Pt2.Distance(Pt);
237 // if(PointDistance > TmpDistance)PointDistance = TmpDistance;
238 // }
239 // return PointDistance;
240 //}
241}
242
243
244std::ostream& operator<<( std::ostream& os, const G4Surface& )
245{
246 // overwrite output operator << to Print out Surface objects
247 // using the PrintOn function defined below
248 // s.PrintOn( os );
249 return os;
250}
251
252
253G4double G4Surface::HowNear( const G4Vector3D& x ) const
254{
255 // Distance from the point x to a Surface.
256 // The default for a Surface is the distance from the point to the origin.
257 G4Vector3D p = G4Vector3D( x - origin );
258 return p.mag();
259}
260
261void G4Surface::Project()
262{
263}
264
265void G4Surface::CalcNormal()
266{
267}
268
269G4int G4Surface::IsConvex() const
270{
271 return -1;
272}
273
274G4int G4Surface::GetConvex() const
275{
276 return 0;
277}
278
279G4int G4Surface::GetNumberOfPoints() const
280{
281 return 0;
282}
283
284const G4Point3D& G4Surface::GetPoint(G4int) const
285{
286 const G4Point3D* tmp= new G4Point3D(0,0,0);
287 return *tmp;
288}
289
290G4Ray* G4Surface::Norm()
291{
292 return (G4Ray*)0;
293}
294
295void G4Surface::Project (G4double& Coord,
296 const G4Point3D& Pt2,
297 const G4Plane& Pl1)
298{
299 Coord = Pt2.x()*Pl1.a + Pt2.y()*Pl1.b + Pt2.z()*Pl1.c - Pl1.d;
300}
301
302/*
303G4double G4Surface::distanceAlongRay( G4int which_way, const G4Ray* ry,
304 G4ThreeVec& p ) const
305{ // Distance along a Ray (straight line with G4ThreeVec) to leave or enter
306 // a Surface. The input variable which_way should be set to +1 to indicate
307 // leaving a Surface, -1 to indicate entering a Surface.
308 // p is the point of intersection of the Ray with the Surface.
309 // This is a default function which just gives the distance
310 // between the origin of the Ray and the origin of the Surface.
311 // Since a generic Surface doesn't have a well-defined Normal, no
312 // further checks are Done.
313
314 // This should always be overwritten for derived classes so Print out
315 // a warning message if this is called.
316 G4cout << "WARNING from Surface::distanceAlongRay\n"
317 << " This function should be overwritten by a derived class.\n"
318 << " Using the Surface base class default.\n";
319 p = GetOrigin();
320 G4ThreeVec d = ry->Position() - p;
321 return d.Magnitude();
322}
323
324G4double G4Surface::distanceAlongHelix( G4int which_way, const Helix* hx,
325 G4ThreeVec& p ) const
326{ // Distance along a Helix to leave or enter a Surface.
327 // The input variable which_way should be set to +1 to indicate
328 // leaving a Surface, -1 to indicate entering a Surface.
329 // p is the point of intersection of the Helix with the Surface.
330 // This is a default function which just gives the distance
331 // between the origin of the Helix and the origin of the Surface.
332 // Since a generic Surface doesn't have a well-defined Normal, no
333 // further checks are Done.
334
335 // This should always be overwritten for derived classes so Print out
336 // a warning message if this is called.
337 G4cout << "WARNING from Surface::distanceAlongHelix\n"
338 << " This function should be overwritten by a derived class.\n"
339 << " Using the Surface base class default.\n";
340 p = GetOrigin();
341 G4ThreeVec d = hx->position() - p;
342 return d.Magnitude();
343}
344
345
346G4ThreeVec G4Surface::Normal() const
347{ // return the Normal unit vector to a Surface
348 // (This is only meaningful for Surfaces for which the Normal does
349 // not depend on location on the Surface).
350 // The default is not well defined, so return ( 0, 0, 0 ).
351 return G4ThreeVec( 0.0, 0.0, 0.0 );
352}
353
354
355G4ThreeVec G4Surface::Normal( const G4ThreeVec& ) const
356{ // return the Normal unit vector to a Surface at the point p on
357 // (or nearly on) the Surface.
358 // The default is not well defined, so return ( 0, 0, 0 ).
359 return G4ThreeVec( 0.0, 0.0, 0.0 );
360}
361
362G4int G4Surface::Inside( const G4ThreeVec& ) const
363{ // return 0 if point p is outside Surface, 1 if Inside
364 // default is not well defined, so return 0
365 return 0;
366}
367
368void G4Surface::move( const G4ThreeVec& p )
369{ // translate origin of Surface by vector p
370 origin += p;
371}
372
373void G4Surface::rotate( G4double alpha, G4double beta,
374 G4double gamma, G4ThreeMat& m, G4int inverse )
375{ // rotate Surface first about global x-axis by angle alpha,
376 // second about global y-axis by angle beta,
377 // and third about global z-axis by angle gamma
378 // by creating and using G4ThreeMat objects
379 // angles are assumed to be given in radians
380 // returns also the overall rotation matrix for use by subclasses
381 // if inverse is non-zero, the order of rotations is reversed
382 // for a generic Surface, only the origin is rotated
383// G4double ax[3][3] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. };
384 G4double ax[3][3];
385 G4double ay[3][3];
386 G4double az[3][3];
387// G4double ay[3][3] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. };
388// G4double az[3][3] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. };
389 ax[0][0] = 1.;
390 ax[1][1] = std::cos( alpha );
391 ax[2][2] = ax[1][1];
392 ax[2][1] = std::sin( alpha );
393 ax[1][2] = -ax[2][1];
394 ay[1][1] = 1.;
395 ay[0][0] = std::cos( beta );
396 ay[2][2] = ay[0][0];
397 ay[0][2] = std::sin( beta );
398 ay[2][0] = -ay[0][2];
399 az[2][2] = 1.;
400 az[0][0] = std::cos( gamma );
401 az[1][1] = az[0][0];
402 az[1][0] = std::sin( gamma );
403 az[0][1] = -az[1][0];
404 G4ThreeMat &Rx = *new G4ThreeMat( ax ); // x-rotation matrix
405 G4ThreeMat &Ry = *new G4ThreeMat( ay ); // y-rotation matrix
406 G4ThreeMat &Rz = *new G4ThreeMat( az ); // z-rotation matrix
407 if ( inverse )
408 m = Rx * ( Ry * Rz );
409 else
410 m = Rz * ( Ry * Rx );
411 origin = m * origin;
412}
413
414void G4Surface::rotate( G4double alpha, G4double beta,
415 G4double gamma, G4int inverse )
416{ // rotate Surface first about global x-axis by angle alpha,
417 // second about global y-axis by angle beta,
418 // and third about global z-axis by angle gamma
419 // by creating and using G4ThreeMat objects
420 // angles are assumed to be given in radians
421 // if inverse is non-zero, the order of rotations is reversed
422 G4ThreeMat m;
423// Just call the above function to do this rotation
424 rotate( alpha, beta, gamma, m, inverse );
425}
426
427*/
Note: See TracBrowser for help on using the repository browser.