Changeset 1315 for trunk/source/geometry/solids

Timestamp:

Jun 18, 2010, 11:42:07 AM (15 years ago)

Author:

garnier

Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:

trunk/source/geometry/solids

Files:

• Boolean/History (modified) (2 diffs)
• Boolean/include/G4BooleanSolid.hh (modified) (3 diffs)
• Boolean/src/G4BooleanSolid.cc (modified) (3 diffs)
• Boolean/src/G4IntersectionSolid.cc (modified) (3 diffs)
• Boolean/src/G4SubtractionSolid.cc (modified) (3 diffs)
• Boolean/src/G4UnionSolid.cc (modified) (3 diffs)
• CSG/History (modified) (2 diffs)
• CSG/include/G4Box.hh (modified) (3 diffs)
• CSG/src/G4Box.cc (modified) (41 diffs)
• CSG/src/G4Orb.cc (modified) (9 diffs)
• specific/History (modified) (2 diffs)
• specific/include/G4ExtrudedSolid.hh (modified) (2 diffs)
• specific/include/G4PolyconeSide.hh (modified) (3 diffs)
• specific/include/G4PolyhedraSide.hh (modified) (3 diffs)
• specific/src/G4ExtrudedSolid.cc (modified) (9 diffs)
• specific/src/G4PolyconeSide.cc (modified) (6 diffs)
• specific/src/G4PolyhedraSide.cc (modified) (7 diffs)
• specific/src/G4SolidExtentList.cc (modified) (3 diffs)
• specific/src/G4TessellatedSolid.cc (modified) (3 diffs)
• specific/src/G4TriangularFacet.cc (modified) (5 diffs)

trunk/source/geometry/solids/Boolean/History

@@ -1 +1 @@
 
-$Id: History,v 1.62 2007/10/24 08:19:11 gcosmo Exp $
+$Id: History,v 1.64 2010/05/11 15:04:01 gcosmo Exp $
 -------------------------------------------------------------------
 
@@ -20 +20 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+ May 11, 2010  J.Allison                  geom-bool-V09-03-00
+  - Introduced recursive algorithm in CreatePolyhedron(): it uses
+    HepPolyhedronProcessor from 'graphics_reps', which tries small shifts
+    in an attempt to avoid numerical problems in the calculation of the
+    polyhedron in BooleanProcessor. Recursion allows HepPolyhedronProcessor
+    to try all permutations, also for Booleans of Booleans.
+    Helps in reducing the number of cases of "Error in Boolean processor" for
+    visualization, but still some stubborn cases are left.
 
  Oct 24, 2007  V.Grichine                 geom-bool-V09-00-00

trunk/source/geometry/solids/Boolean/include/G4BooleanSolid.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4BooleanSolid.hh,v 1.15 2006/10/19 15:34:49 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4BooleanSolid.hh,v 1.17 2010/05/11 15:03:45 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 //
@@ -50 +50 @@
 #include "G4Transform3D.hh"
 
-class G4VSolid;
+class HepPolyhedronProcessor;
 
 class G4BooleanSolid : public G4VSolid
@@ -111 +111 @@
     G4VSolid* fPtrSolidB;
 
+    G4Polyhedron* StackPolyhedron(HepPolyhedronProcessor&,
+                                  const G4VSolid*) const;
+      // Stack polyhedra for processing. Return top polyhedron.
+
   private:
 

trunk/source/geometry/solids/Boolean/src/G4BooleanSolid.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4BooleanSolid.cc,v 1.21 2006/10/19 15:34:49 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4BooleanSolid.cc,v 1.23 2010/05/11 15:03:45 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Implementation for the abstract base class for solids created by boolean 
@@ -40 +40 @@
 #include "G4VSolid.hh"
 #include "G4Polyhedron.hh"
+#include "HepPolyhedronProcessor.h"
 #include "Randomize.hh"
 
@@ -231 +232 @@
   return fpPolyhedron;
 }
+
+//////////////////////////////////////////////////////////////////////////
+//
+// Stacks polyhedra for processing. Returns top polyhedron.
+
+G4Polyhedron*
+G4BooleanSolid::StackPolyhedron(HepPolyhedronProcessor& processor,
+                                const G4VSolid* solid) const
+{
+  HepPolyhedronProcessor::Operation operation;
+  const G4String& type = solid->GetEntityType();
+  if (type == "G4UnionSolid")
+    { operation = HepPolyhedronProcessor::UNION; }
+  else if (type == "G4IntersectionSolid")
+    { operation = HepPolyhedronProcessor::INTERSECTION; }
+  else if (type == "G4SubtractionSolid")
+    { operation = HepPolyhedronProcessor::SUBTRACTION; }
+  else
+  {
+    std::ostringstream oss;
+    oss << "Solid - " << solid->GetName()
+        << " - Unrecognised composite solid"
+        << "\n  Returning NULL !";
+    G4Exception("StackPolyhedron()", "InvalidSetup",
+                JustWarning, oss.str().c_str());
+    return 0;
+  }
+
+  G4Polyhedron* top = 0;
+  const G4VSolid* solidA = solid->GetConstituentSolid(0);
+  const G4VSolid* solidB = solid->GetConstituentSolid(1);
+
+  if (solidA->GetConstituentSolid(0))
+  {
+    top = StackPolyhedron(processor, solidA);
+  }
+  else
+  {
+    top = solidA->GetPolyhedron();
+  }
+  G4Polyhedron* operand = solidB->GetPolyhedron();
+  processor.push_back (operation, *operand);
+
+  return top;
+}

trunk/source/geometry/solids/Boolean/src/G4IntersectionSolid.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4IntersectionSolid.cc,v 1.30 2006/11/08 09:37:41 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4IntersectionSolid.cc,v 1.32 2010/05/11 15:03:45 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Implementation of methods for the class G4IntersectionSolid
@@ -50 +50 @@
 #include "G4VGraphicsScene.hh"
 #include "G4Polyhedron.hh"
+#include "HepPolyhedronProcessor.h"
 #include "G4NURBS.hh"
 // #include "G4NURBSbox.hh"
@@ -502 +503 @@
 G4IntersectionSolid::CreatePolyhedron () const 
 {
-  G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
-  G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
-  if (pA && pB)
-  {
-    G4Polyhedron* resultant = new G4Polyhedron (pA->intersect(*pB));
-    return resultant;
-  }
-  else
-  {
-    std::ostringstream oss;
-    oss << "Solid - " << GetName()
-        << " - one of the Boolean components has no" << G4endl
-        << " corresponding polyhedron. Returning NULL !";
-    G4Exception("G4IntersectionSolid::CreatePolyhedron()", "InvalidSetup",
-                JustWarning, oss.str().c_str());
-    return 0;
-  }
+  HepPolyhedronProcessor processor;
+  // Stack components and components of components recursively
+  // See G4BooleanSolid::StackPolyhedron
+  G4Polyhedron* top = StackPolyhedron(processor, this);
+  G4Polyhedron* result = new G4Polyhedron(*top);
+  if (processor.execute(*result)) { return result; }
+  else { return 0; }
 }
 

trunk/source/geometry/solids/Boolean/src/G4SubtractionSolid.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4SubtractionSolid.cc,v 1.31 2007/10/23 14:42:31 grichine Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4SubtractionSolid.cc,v 1.33 2010/05/11 15:03:45 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Implementation of methods for the class G4IntersectionSolid
@@ -48 +48 @@
 #include "G4VGraphicsScene.hh"
 #include "G4Polyhedron.hh"
+#include "HepPolyhedronProcessor.h"
 #include "G4NURBS.hh"
 // #include "G4NURBSbox.hh"
@@ -463 +464 @@
 G4SubtractionSolid::CreatePolyhedron () const 
 {
-  G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
-  G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
-  if (pA && pB)
-  {
-    G4Polyhedron* resultant = new G4Polyhedron (pA->subtract(*pB));
-    return resultant;
-  }
-  else
-  {
-    std::ostringstream oss;
-    oss << "Solid - " << GetName()
-        << " - one of the Boolean components has no" << G4endl
-        << " corresponding polyhedron. Returning NULL !";
-    G4Exception("G4SubtractionSolid::CreatePolyhedron()", "InvalidSetup",
-                JustWarning, oss.str().c_str());
-    return 0;
-  }
+  HepPolyhedronProcessor processor;
+  // Stack components and components of components recursively
+  // See G4BooleanSolid::StackPolyhedron
+  G4Polyhedron* top = StackPolyhedron(processor, this);
+  G4Polyhedron* result = new G4Polyhedron(*top);
+  if (processor.execute(*result)) { return result; }
+  else { return 0; }
 }
 

trunk/source/geometry/solids/Boolean/src/G4UnionSolid.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4UnionSolid.cc,v 1.35 2007/10/23 14:42:31 grichine Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4UnionSolid.cc,v 1.37 2010/05/11 15:03:45 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Implementation of methods for the class G4IntersectionSolid
@@ -49 +49 @@
 #include "G4VGraphicsScene.hh"
 #include "G4Polyhedron.hh"
+#include "HepPolyhedronProcessor.h"
 #include "G4NURBS.hh"
 // #include "G4NURBSbox.hh"
@@ -454 +455 @@
 G4UnionSolid::CreatePolyhedron () const 
 {
-  G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
-  G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
-  if (pA && pB) {
-    G4Polyhedron* resultant = new G4Polyhedron (pA->add(*pB));
-    return resultant;
-  } else {
-    std::ostringstream oss;
-    oss << GetName() <<
-      ": one of the Boolean components has no corresponding polyhedron.";
-    G4Exception("G4UnionSolid::CreatePolyhedron",
-                "", JustWarning, oss.str().c_str());
-    return 0;
-  }
+  HepPolyhedronProcessor processor;
+  // Stack components and components of components recursively
+  // See G4BooleanSolid::StackPolyhedron
+  G4Polyhedron* top = StackPolyhedron(processor, this);
+  G4Polyhedron* result = new G4Polyhedron(*top);
+  if (processor.execute(*result)) { return result; }
+  else { return 0; }
 }
 

trunk/source/geometry/solids/CSG/History

@@ -1 @@
-$Id: History,v 1.119 2009/11/30 10:20:53 gcosmo Exp $
+$Id: History,v 1.121 2010/05/25 09:19:02 gcosmo Exp $
 -------------------------------------------------------------------
 
@@ -17 +17 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+May 25, 2010  G.Cosmo                    geom-csg-V09-03-01
+- G4Box: more minor fixes in coding style, following code review discussion.
+
+May 18, 2010  G.Cosmo                    geom-csg-V09-03-00
+- G4Box: first implementation of speed improvements and corrections from joint
+  code review of G4Box class: introduced half-tolerance constants and better
+  logic in some conditional statements.
+- Some code cleanup.
+
+Dec 04, 2009  V.Grichine
+- G4Orb: modified DistanceToIn(p,v) to be more consistent with Inside(p).
 
 Nov 30, 2009  V.Grichine                 geom-csg-V09-02-08

trunk/source/geometry/solids/CSG/include/G4Box.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4Box.hh,v 1.17 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4Box.hh,v 1.18 2010/05/18 10:07:52 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // --------------------------------------------------------------------
@@ -101 +101 @@
     std::ostream& StreamInfo(std::ostream& os) const;
 
-  // Functions for visualization
+  // Utilities for visualization
 
     void          DescribeYourselfTo (G4VGraphicsScene& scene) const;
@@ -125 +125 @@
 
     enum ESide {kUndefined,kPX,kMX,kPY,kMY,kPZ,kMZ};
-      // Codes for faces (kPX=plus x face,kMY= minus y face etc)
+      // Codes for faces (kPX= +x face, kMY= -y face, etc...)
 
   private:

trunk/source/geometry/solids/CSG/src/G4Box.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4Box.cc,v 1.49 2010/05/25 10:14:41 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -32 +32 @@
 // Implementation for G4Box class
 //
-//  24.06.98 - V. Grichine: insideEdge in DistanceToIn(p,v)
+//  24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v)
 //  20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
 //  07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0
@@ -67 +67 @@
   if ( (pX > 2*kCarTolerance)
     && (pY > 2*kCarTolerance)
-    && (pZ > 2*kCarTolerance) )
+    && (pZ > 2*kCarTolerance) )  // limit to thickness of surfaces
   {
     fDx = pX ;
@@ -105 +105 @@
 void G4Box::SetXHalfLength(G4double dx)
 {
-  if(dx > 2*kCarTolerance)
+  if(dx > 2*kCarTolerance)  // limit to thickness of surfaces
+  {
     fDx = dx;
+  }
   else
   {
@@ -122 +124 @@
 void G4Box::SetYHalfLength(G4double dy) 
 {
-  if(dy > 2*kCarTolerance)
+  if(dy > 2*kCarTolerance)  // limit to thickness of surfaces
+  {
     fDy = dy;
+  }
   else
   {
@@ -139 +143 @@
 void G4Box::SetZHalfLength(G4double dz) 
 {
-  if(dz > 2*kCarTolerance)
+  if(dz > 2*kCarTolerance)  // limit to thickness of surfaces
+  {
     fDz = dz;
+  }
   else
   {
@@ -153 +159 @@
   fpPolyhedron = 0;
 } 
-    
-
 
 ////////////////////////////////////////////////////////////////////////
@@ -193 +197 @@
     if (pVoxelLimit.IsXLimited())
     {
-      if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance || 
-           xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance    ) return false ;
+      if ((xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) || 
+          (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance)) { return false ; }
       else
       {
-        if (xMin < pVoxelLimit.GetMinXExtent())
-        {
-          xMin = pVoxelLimit.GetMinXExtent() ;
-        }
-        if (xMax > pVoxelLimit.GetMaxXExtent())
-        {
-          xMax = pVoxelLimit.GetMaxXExtent() ;
-        }
+        xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
+        xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
       }
     }
@@ -213 +211 @@
     if (pVoxelLimit.IsYLimited())
     {
-      if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
-           yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance   ) return false ;
+      if ((yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) ||
+          (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance)) { return false ; }
       else
       {
-        if (yMin < pVoxelLimit.GetMinYExtent())
-        {
-          yMin = pVoxelLimit.GetMinYExtent() ;
-        }
-        if (yMax > pVoxelLimit.GetMaxYExtent())
-        {
-          yMax = pVoxelLimit.GetMaxYExtent() ;
-        }
+        yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
+        yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
       }
     }
@@ -233 +225 @@
     if (pVoxelLimit.IsZLimited())
     {
-      if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
-           zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance   ) return false ;
+      if ((zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) ||
+          (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance)) { return false ; }
       else
       {
-        if (zMin < pVoxelLimit.GetMinZExtent())
-        {
-          zMin = pVoxelLimit.GetMinZExtent() ;
-        }
-        if (zMax > pVoxelLimit.GetMaxZExtent())
-        {
-          zMax = pVoxelLimit.GetMaxZExtent() ;
-        }
+        zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
+        zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
       }
     }
@@ -286 +272 @@
     if (pVoxelLimit.IsLimited(pAxis) == false) 
     {  
-      if ( pMin != kInfinity || pMax != -kInfinity ) 
+      if ( (pMin != kInfinity) || (pMax != -kInfinity) ) 
       {
         existsAfterClip = true ;
@@ -292 +278 @@
         // Add 2*tolerance to avoid precision troubles
 
-        pMin           -= kCarTolerance;
-        pMax           += kCarTolerance;
+        pMin -= kCarTolerance;
+        pMax += kCarTolerance;
       }
     }      
@@ -303 +289 @@
        ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
 
-      if ( pMin != kInfinity || pMax != -kInfinity )
+      if ( (pMin != kInfinity) || (pMax != -kInfinity) )
       {
         existsAfterClip = true ;
@@ -356 +342 @@
 EInside G4Box::Inside(const G4ThreeVector& p) const
 {
+  static const G4double delta=0.5*kCarTolerance;
   EInside in = kOutside ;
-
-  if ( std::fabs(p.x()) <= fDx - kCarTolerance*0.5 )
-  {
-    if (std::fabs(p.y()) <= fDy - kCarTolerance*0.5 )
-    {
-      if      (std::fabs(p.z()) <= fDz - kCarTolerance*0.5 ) in = kInside ;
-      else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
-    }
-    else if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
-    {
-      if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
-    }
-  }
-  else if (std::fabs(p.x()) <= fDx + kCarTolerance*0.5 )
-  {
-    if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
-    {
-      if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5) in = kSurface ;
+  G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
+
+  if ( q.x() <= (fDx - delta) )
+  {
+    if (q.y() <= (fDy - delta) )
+    {
+      if      ( q.z() <= (fDz - delta) ) { in = kInside ;  }
+      else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
+    }
+    else if ( q.y() <= (fDy + delta) )
+    {
+      if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
+    }
+  }
+  else if ( q.x() <= (fDx + delta) )
+  {
+    if ( q.y() <= (fDy + delta) )
+    {
+      if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
     }
   }
@@ -389 +377 @@
 {
   G4double distx, disty, distz ;
-  G4ThreeVector norm ;
+  G4ThreeVector norm(0.,0.,0.);
 
   // Calculate distances as if in 1st octant
@@ -400 +388 @@
   // normals
 
-  const G4double delta    = 0.5*kCarTolerance;
+  static const G4double delta    = 0.5*kCarTolerance;
   const G4ThreeVector nX  = G4ThreeVector( 1.0, 0,0  );
   const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0  );
@@ -415 +403 @@
   {
     noSurfaces ++; 
-    if ( p.x() >= 0.){        // on +X surface
-      normX= nX ;    // G4ThreeVector( 1.0, 0., 0. );
-    }else{
-      normX= nmX;    // G4ThreeVector(-1.0, 0., 0. ); 
-    }
+    if ( p.x() >= 0. )  { normX= nX ; }       // on +X surface : (1,0,0)
+    else                { normX= nmX; }       //                 (-1,0,0)
     sumnorm= normX; 
   }
@@ -426 +411 @@
   {
     noSurfaces ++; 
-    if ( p.y() >= 0.){        // on +Y surface
-      normY= nY;
-    }else{
-      normY = nmY; 
-    }
+    if ( p.y() >= 0. )  { normY= nY;  }       // on +Y surface
+    else                { normY= nmY; }
     sumnorm += normY; 
   }
@@ -437 +419 @@
   {
     noSurfaces ++; 
-    if ( p.z() >= 0.){        // on +Z surface
-      normZ= nZ;
-    }else{
-      normZ = nmZ; 
-    }
-    sumnorm += normZ; 
-  }
-
-  // sumnorm= normX + normY + normZ; 
-  const G4double invSqrt2 = 1.0 / std::sqrt( 2.0); 
-  const G4double invSqrt3 = 1.0 / std::sqrt( 3.0); 
-
-  norm= G4ThreeVector( 0., 0., 0.); 
+    if ( p.z() >= 0. )  { normZ= nZ;  }       // on +Z surface
+    else                { normZ= nmZ; }
+    sumnorm += normZ;
+  }
+
+  static const G4double invSqrt2 = 1.0 / std::sqrt(2.0); 
+  static const G4double invSqrt3 = 1.0 / std::sqrt(3.0); 
+
   if( noSurfaces > 0 )
   { 
-    if( noSurfaces == 1 ){ 
+    if( noSurfaces == 1 )
+    { 
       norm= sumnorm; 
-    }else{
+    }
+    else
+    {
       // norm = sumnorm . unit(); 
-      if( noSurfaces == 2 ) { 
+      if( noSurfaces == 2 )
+      { 
         // 2 surfaces -> on edge 
         norm = invSqrt2 * sumnorm; 
-      } else { 
+      }
+      else
+      { 
         // 3 surfaces (on corner)
         norm = invSqrt3 * sumnorm; 
       }
     }
-  }else{
+  }
+  else
+  {
 #ifdef G4CSGDEBUG
      G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning, 
@@ -483 +468 @@
 {
   G4double distx, disty, distz ;
-  G4ThreeVector norm ;
+  G4ThreeVector norm(0.,0.,0.);
 
   // Calculate distances as if in 1st octant
@@ -495 +480 @@
     if ( distx <= distz )     // Closest to X
     {
-      if ( p.x() < 0 ) norm = G4ThreeVector(-1.0,0,0) ;
-      else             norm = G4ThreeVector( 1.0,0,0) ;
+      if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
+      else             { norm = G4ThreeVector( 1.0,0,0) ; }
     }
     else                      // Closest to Z
     {
-      if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
-      else             norm = G4ThreeVector(0,0, 1.0) ;
+      if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
+      else             { norm = G4ThreeVector(0,0, 1.0) ; }
     }
   }
@@ -508 +493 @@
     if ( disty <= distz )      // Closest to Y
     {
-      if ( p.y() < 0 ) norm = G4ThreeVector(0,-1.0,0) ;
-      else             norm = G4ThreeVector(0, 1.0,0) ;
+      if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
+      else             { norm = G4ThreeVector(0, 1.0,0) ; }
     }
     else                       // Closest to Z
     {
-      if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
-      else             norm = G4ThreeVector(0,0, 1.0) ;
+      if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
+      else             { norm = G4ThreeVector(0,0, 1.0) ; }
     }
   }
@@ -541 +526 @@
 // shape.
 
-G4double G4Box::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const
+G4double G4Box::DistanceToIn(const G4ThreeVector& p,
+                             const G4ThreeVector& v) const
 {
   G4double safx, safy, safz ;
@@ -549 +535 @@
   G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
 
+  static const G4double delta = 0.5*kCarTolerance;
+
   safx = std::fabs(p.x()) - fDx ;     // minimum distance to x surface of shape
   safy = std::fabs(p.y()) - fDy ;
@@ -558 +546 @@
   // travel is in a direction away from the shape.
 
-  if (    ((p.x()*v.x() >= 0.0) && safx > -kCarTolerance*0.5) 
-       || ((p.y()*v.y() >= 0.0) && safy > -kCarTolerance*0.5)
-       || ((p.z()*v.z() >= 0.0) && safz > -kCarTolerance*0.5)   ) 
+  if (    ((p.x()*v.x() >= 0.0) && (safx > -delta)) 
+       || ((p.y()*v.y() >= 0.0) && (safy > -delta))
+       || ((p.z()*v.z() >= 0.0) && (safz > -delta))   ) 
   {
     return kInfinity ;  // travel away or parallel within tolerance
@@ -568 +556 @@
   // X Planes
 
-  if ( v.x())
+  if ( v.x() )  // != 0
   {
     stmp = 1.0/std::fabs(v.x()) ;
@@ -579 +567 @@
     else
     {
-      if (v.x() > 0)  sOut = (fDx - p.x())*stmp ;
-      if (v.x() < 0)  sOut = (fDx + p.x())*stmp ;
+      if (v.x() < 0)  { sOut = (fDx + p.x())*stmp ; }
+      else            { sOut = (fDx - p.x())*stmp ; }
     }
   }
@@ -586 +574 @@
   // Y Planes
 
-  if ( v.y()) 
+  if ( v.y() )  // != 0
   {
     stmp = 1.0/std::fabs(v.y()) ;
@@ -595 +583 @@
       smaxy = (fDy+std::fabs(p.y()))*stmp ;
 
-      if (sminy > smin) smin=sminy ;
-      if (smaxy < smax) smax=smaxy ;
-
-      if (smin >= smax-kCarTolerance*0.5)
+      if (sminy > smin) { smin=sminy ; }
+      if (smaxy < smax) { smax=smaxy ; }
+
+      if (smin >= (smax-delta))
       {
         return kInfinity ;  // touch XY corner
@@ -605 +593 @@
     else
     {
-      if (v.y() > 0)  sOuty = (fDy - p.y())*stmp ;
-      if (v.y() < 0)  sOuty = (fDy + p.y())*stmp ;
-      if( sOuty < sOut ) sOut = sOuty ;
+      if (v.y() < 0)  { sOuty = (fDy + p.y())*stmp ; }
+      else            { sOuty = (fDy - p.y())*stmp ; }
+      if( sOuty < sOut ) { sOut = sOuty ; }
     }     
   }
@@ -613 +601 @@
   // Z planes
 
-  if ( v.z() )
+  if ( v.z() )  // != 0
   {
     stmp = 1.0/std::fabs(v.z()) ;
 
-    if ( safz >= 0.0)
+    if ( safz >= 0.0 )
     {
       sminz = safz*stmp ;
       smaxz = (fDz+std::fabs(p.z()))*stmp ;
 
-      if (sminz > smin) smin = sminz ;
-      if (smaxz < smax) smax = smaxz ;
-
-      if (smin >= smax-kCarTolerance*0.5)
+      if (sminz > smin) { smin = sminz ; }
+      if (smaxz < smax) { smax = smaxz ; }
+
+      if (smin >= (smax-delta))
       { 
         return kInfinity ;    // touch ZX or ZY corners
@@ -632 +620 @@
     else
     {
-      if (v.z() > 0)  sOutz = (fDz - p.z())*stmp ;
-      if (v.z() < 0)  sOutz = (fDz + p.z())*stmp ;
-      if( sOutz < sOut ) sOut = sOutz ;
-    }
-  }
-
-  if ( sOut <= smin + 0.5*kCarTolerance) // travel over edge
+      if (v.z() < 0)  { sOutz = (fDz + p.z())*stmp ; }
+      else            { sOutz = (fDz - p.z())*stmp ; }
+      if( sOutz < sOut ) { sOut = sOutz ; }
+    }
+  }
+
+  if (sOut <= (smin + delta)) // travel over edge
   {
     return kInfinity ;
   }
-  if (smin < 0.5*kCarTolerance)  smin = 0.0 ;
+  if (smin < delta)  { smin = 0.0 ; }
 
   return smin ;
@@ -662 +650 @@
   safez = std::fabs(p.z()) - fDz ;
 
-  if (safex > safe) safe = safex ;
-  if (safey > safe) safe = safey ;
-  if (safez > safe) safe = safez ;
+  if (safex > safe) { safe = safex ; }
+  if (safey > safe) { safe = safey ; }
+  if (safez > safe) { safe = safez ; }
 
   return safe ;
@@ -671 +659 @@
 /////////////////////////////////////////////////////////////////////////
 //
-// Calcluate distance to surface of box from inside
+// Calculate distance to surface of box from inside
 // by calculating distances to box's x/y/z planes.
 // Smallest distance is exact distance to exiting.
@@ -682 +670 @@
 {
   ESide side = kUndefined ;
-  G4double pdist,stmp,snxt;
-
-  if (calcNorm) *validNorm = true ; // All normals are valid
-
-  if (v.x() > 0)   // X planes  
+  G4double pdist,stmp,snxt=kInfinity;
+
+  static const G4double delta = 0.5*kCarTolerance;
+
+  if (calcNorm) { *validNorm = true ; }  // All normals are valid
+
+  if (v.x() > 0)   // X planes
   {
     pdist = fDx - p.x() ;
 
-    if (pdist > kCarTolerance*0.5)
+    if (pdist > delta)
     {
       snxt = pdist/v.x() ;
@@ -697 +687 @@
     else
     {
-      if (calcNorm) *n    = G4ThreeVector(1,0,0) ;
-      return         snxt = 0 ;
-    }
-  }
-  else if (v.x() < 0) 
+      if (calcNorm) { *n   = G4ThreeVector(1,0,0) ; }
+      return        snxt = 0 ;
+    }
+  }
+  else if (v.x() < 0)
   {
     pdist = fDx + p.x() ;
 
-    if (pdist > kCarTolerance*0.5)
+    if (pdist > delta)
     {
       snxt = -pdist/v.x() ;
@@ -712 +702 @@
     else
     {
-      if (calcNorm) *n   = G4ThreeVector(-1,0,0) ;
+      if (calcNorm) { *n   = G4ThreeVector(-1,0,0) ; }
       return        snxt = 0 ;
     }
   }
-  else snxt = kInfinity ;
-
-  if ( v.y() > 0 )   // Y planes  
-  {
-    pdist=fDy-p.y();
-
-    if (pdist>kCarTolerance*0.5)
-    {
-      stmp=pdist/v.y();
-
-      if (stmp<snxt)
-      {
-        snxt=stmp;
-        side=kPY;
-      }
-    }
-    else
-    {
-      if (calcNorm) *n    = G4ThreeVector(0,1,0) ;
-      return         snxt = 0 ;
-    }
-  }
-  else if ( v.y() < 0 ) 
-  {
-    pdist = fDy + p.y() ;
-
-    if (pdist > kCarTolerance*0.5)
-    {
-      stmp=-pdist/v.y();
-
-      if (stmp<snxt)
-      {
-        snxt=stmp;
-        side=kMY;
-      }
-    }
-    else
-    {
-      if (calcNorm) *n    = G4ThreeVector(0,-1,0) ;
-      return         snxt = 0 ;
-    }
-  }
-  if (v.z()>0)        // Z planes 
-  {
-    pdist=fDz-p.z();
-
-    if (pdist > kCarTolerance*0.5)
-    {
-      stmp=pdist/v.z();
+
+  if (v.y() > 0)   // Y planes
+  {
+    pdist = fDy-p.y();
+
+    if (pdist > delta)
+    {
+      stmp = pdist/v.y();
 
       if (stmp < snxt)
       {
-        snxt=stmp;
-        side=kPZ;
+        snxt = stmp;
+        side = kPY;
       }
     }
     else
     {
-      if (calcNorm) *n    = G4ThreeVector(0,0,1) ;
-      return         snxt = 0 ;
-    }
-  }
-  else if (v.z()<0) 
-  {
-    pdist = fDz + p.z() ;
-
-    if (pdist > kCarTolerance*0.5)
-    {
-      stmp=-pdist/v.z();
-
-      if (stmp < snxt)
+      if (calcNorm) { *n   = G4ThreeVector(0,1,0) ; }
+      return        snxt = 0 ;
+    }
+  }
+  else if (v.y() < 0)
+  {
+    pdist = fDy + p.y() ;
+
+    if (pdist > delta)
+    {
+      stmp = -pdist/v.y();
+
+      if ( stmp < snxt )
       {
-        snxt=stmp;
-        side=kMZ;
+        snxt = stmp;
+        side = kMY;
       }
     }
     else
     {
-      if (calcNorm) *n    = G4ThreeVector(0,0,-1) ;
-      return         snxt = 0 ;
-    }
-  }
+      if (calcNorm) { *n   = G4ThreeVector(0,-1,0) ; }
+      return        snxt = 0 ;
+    }
+  }
+
+  if (v.z() > 0)        // Z planes
+  {
+    pdist = fDz-p.z();
+
+    if ( pdist > delta )
+    {
+      stmp = pdist/v.z();
+
+      if ( stmp < snxt )
+      {
+        snxt = stmp;
+        side = kPZ;
+      }
+    }
+    else
+    {
+      if (calcNorm) { *n   = G4ThreeVector(0,0,1) ; } 
+      return        snxt = 0 ;
+    }
+  }
+  else if (v.z() < 0)
+  {
+    pdist = fDz + p.z();
+
+    if ( pdist > delta )
+    {
+      stmp = -pdist/v.z();
+
+      if ( stmp < snxt )
+      {
+        snxt = stmp;
+        side = kMZ;
+      }
+    }
+    else
+    {
+      if (calcNorm) { *n   = G4ThreeVector(0,0,-1) ; }
+      return        snxt = 0 ;
+    }
+  }
+
   if (calcNorm)
   {      
@@ -875 +866 @@
   // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
 
-  if (safx2 < safx1) safe = safx2 ;
-  else               safe = safx1 ;
-  if (safy1 < safe)  safe = safy1 ;
-  if (safy2 < safe)  safe = safy2 ;
-  if (safz1 < safe)  safe = safz1 ;
-  if (safz2 < safe)  safe = safz2 ;
-
-  if (safe < 0) safe = 0 ;
+  if (safx2 < safx1) { safe = safx2; }
+  else               { safe = safx1; }
+  if (safy1 < safe)  { safe = safy1; }
+  if (safy2 < safe)  { safe = safy2; }
+  if (safz1 < safe)  { safe = safz1; }
+  if (safz2 < safe)  { safe = safz2; }
+
+  if (safe < 0) { safe = 0 ; }
   return safe ;  
 }
@@ -979 +970 @@
     py = -fDy +2*fDy*G4UniformRand();
 
-    if(G4UniformRand() > 0.5) pz =  fDz;
-    else                      pz = -fDz;
+    if(G4UniformRand() > 0.5) { pz =  fDz; }
+    else                      { pz = -fDz; }
   }
   else if ( ( select - Sxy ) < Sxz ) 
@@ -987 +978 @@
     pz = -fDz +2*fDz*G4UniformRand();
 
-    if(G4UniformRand() > 0.5) py =  fDy;
-    else                      py = -fDy;
+    if(G4UniformRand() > 0.5) { py =  fDy; }
+    else                      { py = -fDy; }
   }
   else  
@@ -995 +986 @@
     pz = -fDz +2*fDz*G4UniformRand();
 
-    if(G4UniformRand() > 0.5) px =  fDx;
-    else                      px = -fDx;
+    if(G4UniformRand() > 0.5) { px =  fDx; }
+    else                      { px = -fDx; }
   } 
   return G4ThreeVector(px,py,pz);

trunk/source/geometry/solids/CSG/src/G4Orb.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Orb.cc,v 1.30 2009/11/30 10:20:38 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4Orb.cc,v 1.31 2009/12/04 15:39:56 grichine Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // class G4Orb
@@ -296 +296 @@
 
 
-  rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z() ;
+  rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z();
 
   G4double rad = std::sqrt(rad2);
@@ -304 +304 @@
   // sets `in'
   
-  tolRMax = fRmax - fRmaxTolerance*0.5 ;
+  tolRMax = fRmax - fRmaxTolerance*0.5;
     
-  if ( rad <= tolRMax )  { in = kInside ; }
+  if ( rad <= tolRMax )  { in = kInside; }
   else
   {
-    tolRMax = fRmax + fRmaxTolerance*0.5 ;       
-    if ( rad <= tolRMax )  { in = kSurface ; }
-    else                   { in = kOutside ; }
+    tolRMax = fRmax + fRmaxTolerance*0.5;       
+    if ( rad <= tolRMax )  { in = kSurface; }
+    else                   { in = kOutside; }
   }
   return in;
@@ -357 +357 @@
                               const G4ThreeVector& v  ) const
 {
-  G4double snxt = kInfinity ;      // snxt = default return value
-
-  G4double rad2, pDotV3d, tolORMax2, tolIRMax2 ;
-  G4double c, d2, s = kInfinity ;
+  G4double snxt = kInfinity;      // snxt = default return value
+
+  G4double rad, pDotV3d; // , tolORMax2, tolIRMax2;
+  G4double c, d2, s = kInfinity;
 
   const G4double dRmax = 100.*fRmax;
@@ -366 +366 @@
   // General Precalcs
 
-  rad2    = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
-  pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
+  rad    = std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
+  pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z();
 
   // Radial Precalcs
 
-  tolORMax2 = (fRmax+fRmaxTolerance*0.5)*(fRmax+fRmaxTolerance*0.5) ;
-  tolIRMax2 = (fRmax-fRmaxTolerance*0.5)*(fRmax-fRmaxTolerance*0.5) ;
+  // tolORMax2 = (fRmax+fRmaxTolerance*0.5)*(fRmax+fRmaxTolerance*0.5);
+  // tolIRMax2 = (fRmax-fRmaxTolerance*0.5)*(fRmax-fRmaxTolerance*0.5);
 
   // Outer spherical shell intersection
@@ -388 +388 @@
   // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
 
-
-  G4double rad = std::sqrt(rad2);
   c = (rad - fRmax)*(rad + fRmax);
 
-  if ( c > fRmaxTolerance*fRmax )
-  {
-    // If outside tolerant boundary of outer G4Orb
-    // [ should be std::sqrt(rad2) - fRmax > fRmaxTolerance*0.5 ]
-
-    d2 = pDotV3d*pDotV3d - c ;
-
-    if ( d2 >= 0 )
-    {
-      s = -pDotV3d - std::sqrt(d2) ;
-      if ( s >= 0 )
-      {
-        if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on
-        {              // 64 bits systems. Split long distances and recompute
-          G4double fTerm = s-std::fmod(s,dRmax);
-          s = fTerm + DistanceToIn(p+fTerm*v,v);
-        } 
-        return snxt = s;
-      }
-    }
-    else    // No intersection with G4Orb
-    {
-      return snxt = kInfinity;
-    }
-  }
-  else
-  {
-    if ( c > -fRmaxTolerance*fRmax )  // on surface  
-    {
-      d2 = pDotV3d*pDotV3d - c ;             
-      if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
+  if( rad > fRmax-fRmaxTolerance*0.5 ) // not inside in terms of Inside(p)
+  {
+    if ( c > fRmaxTolerance*fRmax )
+    {
+      // If outside tolerant boundary of outer G4Orb in terms of c
+      // [ should be std::sqrt(rad2) - fRmax > fRmaxTolerance*0.5 ]
+
+      d2 = pDotV3d*pDotV3d - c;
+
+      if ( d2 >= 0 )
+      {
+        s = -pDotV3d - std::sqrt(d2);
+        if ( s >= 0 )
+        {
+          if ( s > dRmax ) // Avoid rounding errors due to precision issues seen on
+          {                // 64 bits systems. Split long distances and recompute
+            G4double fTerm = s - std::fmod(s,dRmax);
+            s = fTerm + DistanceToIn(p+fTerm*v,v);
+          } 
+          return snxt = s;
+        }
+      }
+      else    // No intersection with G4Orb
       {
         return snxt = kInfinity;
       }
-      else
-      {
-        return snxt = 0.;
-      }
-    }
+    }
+    else // not outside in terms of c
+    {
+      if ( c > -fRmaxTolerance*fRmax )  // on surface  
+      {
+        d2 = pDotV3d*pDotV3d - c;             
+        if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
+        {
+          return snxt = kInfinity;
+        }
+        else
+        {
+          return snxt = 0.;
+        }
+      }
+    }
+  }
 #ifdef G4CSGDEBUG
-    else // inside ???
-    {
+  else // inside ???
+  {
       G4Exception("G4Orb::DistanceToIn(p,v)", "Notification",
                   JustWarning, "Point p is inside !?");
-    }
+  }
 #endif
-  }
+
   return snxt;
 }
@@ -520 +522 @@
         if(calcNorm)
         {
-          *validNorm = true ;
-          *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
+          *validNorm = true;
+          *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax);
         }
         return snxt = 0;
@@ -528 +530 @@
       {
         snxt = -pDotV3d + std::sqrt(d2);    // second root since inside Rmax
-        side = kRMax ; 
+        side = kRMax; 
       }
     }
@@ -596 +598 @@
   if( Inside(p) == kOutside )
   {
-     G4cout.precision(16) ;
-     G4cout << G4endl ;
+     G4cout.precision(16);
+     G4cout << G4endl;
      DumpInfo();
-     G4cout << "Position:"  << G4endl << G4endl ;
-     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
-     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
-     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
+     G4cout << "Position:"  << G4endl << G4endl;
+     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
+     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
+     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
      G4Exception("G4Orb::DistanceToOut(p)", "Notification", JustWarning, 
                  "Point p is outside !?" );

trunk/source/geometry/solids/specific/History

@@ -1 @@
-$Id: History,v 1.160 2009/11/11 12:25:46 gcosmo Exp $
+$Id: History,v 1.170 2010/06/11 09:42:59 gcosmo Exp $
 -------------------------------------------------------------------
 
@@ -17 +17 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+11-Jun-2010  T.Nikitina (geom-specific-V09-03-07)
+- G4GenericTrap: fixed cases of zero dToIn/sToOut when vertices are collapsed
+  to triangle, line or point. Added new methods handling those specific cases.
+- Added unit test for G4GenericTrap.
+
+10-Jun-2010  I.Hrivnacova
+- G4GenericTrap: fixed parameter names in CalculateExtent() for the test
+  case construction through tessellated facets.
+
+09-Jun-2010  G.Cosmo
+- G4GenericTrap: moved internal methods to private section and reordered
+  in source file. Added missing implementation for IsTwisted() method.
+
+03-Jun-2010  T.Nikitina, G.Cosmo (geom-specific-V09-03-06)
+- G4GenericTrap:
+  o Fixed initialization of fSurfaceArea and fCubicVolume and
+    calculation of surface area.
+  o Fixed error in Inside(p) function, and corrected use std::fabs() instead
+    of std::abs() for floating point values.
+  o Added missing initialisation of fpPolyhedron pointer.
+  o More corrected signatures for use of non-const references for vectors
+    passed as arguments to functions.
+
+02-Jun-2010  G.Cosmo (geom-specific-V09-03-05)
+- G4GenericTrap: use const reference for vector of vertices passed as argument
+  in constructor and accessor.
+
+27-May-2010  T.Nikitina (geom-specific-V09-03-04)
+- First implementation of G4GenericTrap shape, a new solid representing an
+  arbitrary trapezoid with up to 8 vertices standing on two parallel planes
+  perpendicular to the Z axis.
+
+28-Apr-2010  P.R.Truscott (geom-specific-V09-03-03)
+- Fix in G4TriangularFacet and G4TessellatedSolid to correct treatment of
+  optical photon transport related to internal reflection at surface.
+  Addresses problem report #1103.
+
+15-Apr-2010  I.Hrivnacova (geom-specific-V09-03-02)
+- G4ExtrudedSolid: eliminated requirement for clockwise ordering of polygon
+  vertices. Added a check for vertices ordering; if vertices are defined
+  anti-clockwise their ordering is reverted.
+  Fix in polygon facet triangularization for consequent concave vertices.
+
+24-Feb-2010  G.Cosmo (geom-specific-V09-03-01)
+- Adopt caching of Phi in G4PolyconeSide and G4PolyhedraSide to avoid
+  unnecessary consecutive computations on the same point.
+
+10-Feb-2010  G.Cosmo (geom-specific-V09-03-00)
+- Use kInfinity for initialising minimum and maximum allowed extent for
+  G4SolidExtentList of faceted solids.
 
 11-Nov-2009  G.Cosmo (geom-specific-V09-02-08)

trunk/source/geometry/solids/specific/include/G4ExtrudedSolid.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4ExtrudedSolid.hh,v 1.7 2008/02/27 12:32:48 ivana Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4ExtrudedSolid.hh,v 1.8 2010/04/15 10:23:34 ivana Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -45 +45 @@
 // const G4String& pName             - solid name
 // std::vector<G4TwoVector> polygon  - the vertices of the outlined polygon
-//                                     defined in clock-wise order     
+//                                     defined in clockwise or anti-clockwise order     
 // std::vector<ZSection>             - the z-sections defined by
 //                                     z position, offset and scale

trunk/source/geometry/solids/specific/include/G4PolyconeSide.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4PolyconeSide.hh,v 1.12 2008/05/15 11:41:58 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4PolyconeSide.hh,v 1.13 2010/02/24 11:18:25 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -114 +114 @@
   protected:
 
+    G4double DistanceAway( const G4ThreeVector &p, G4bool opposite,
+                                 G4double &distOutside2, G4double *rzNorm=0 );
+      
+    G4bool PointOnCone( const G4ThreeVector &hit, G4double normSign,
+                        const G4ThreeVector &p,
+                        const G4ThreeVector &v, G4ThreeVector &normal );
+
+    void CopyStuff( const G4PolyconeSide &source );
+  
+    static void FindLineIntersect( G4double x1, G4double y1,
+                                   G4double tx1, G4double ty1,
+                                   G4double x2, G4double y2,
+                                 G4double tx2, G4double ty2,
+                                 G4double &x, G4double &y );
+
+    G4double GetPhi( const G4ThreeVector& p );
+
+  protected:
+
     G4double r[2], z[2]; // r, z parameters, in specified order
     G4double startPhi,   // Start phi (0 to 2pi), if phiIsOpen
@@ -136 +155 @@
     G4ThreeVector *corners; // The coordinates of the corners (if phiIsOpen)
 
-    G4double DistanceAway( const G4ThreeVector &p, G4bool opposite,
-                                 G4double &distOutside2, G4double *rzNorm=0 );
-      
-    G4bool PointOnCone( const G4ThreeVector &hit, G4double normSign,
-                        const G4ThreeVector &p,
-                        const G4ThreeVector &v, G4ThreeVector &normal );
-
-    void CopyStuff( const G4PolyconeSide &source );
-  
-    static void FindLineIntersect( G4double x1, G4double y1,
-                                   G4double tx1, G4double ty1,
-                                   G4double x2, G4double y2,
-                                 G4double tx2, G4double ty2,
-                                 G4double &x, G4double &y );
   private:
 
+    std::pair<G4ThreeVector, G4double> fPhi;  // Cached value for phi
     G4double kCarTolerance; // Geometrical surface thickness
     G4double fSurfaceArea;  // Used for surface calculation 

trunk/source/geometry/solids/specific/include/G4PolyhedraSide.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4PolyhedraSide.hh,v 1.11 2008/05/15 11:41:59 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4PolyhedraSide.hh,v 1.12 2010/02/24 11:18:25 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -166 +166 @@
   
     G4int PhiSegment( G4double phi );
-  
+
+    G4double GetPhi( const G4ThreeVector& p );
+
     G4double DistanceToOneSide( const G4ThreeVector &p,
                                 const G4PolyhedraSideVec &vec,
@@ -197 +199 @@
   private:
 
+    std::pair<G4ThreeVector, G4double> fPhi;  // Cached value for phi
     G4double kCarTolerance;  // Geometrical surface thickness
     G4double fSurfaceArea;   // Surface Area 

trunk/source/geometry/solids/specific/src/G4ExtrudedSolid.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4ExtrudedSolid.cc,v 1.18 2008/10/30 11:47:45 ivana Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4ExtrudedSolid.cc,v 1.19 2010/04/15 10:23:34 ivana Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 //
@@ -40 +40 @@
 #include <algorithm>
 #include <cmath>
+#include <iomanip>
 
 #include "G4ExtrudedSolid.hh"
@@ -94 +95 @@
         "Z-sections with the same z position are not supported.");
     }
-  }      
-
+  }  
+  
+  // Check if polygon vertices are defined clockwise
+  // (the area is positive if polygon vertices are defined anti-clockwise)
+  //
+  G4double area = 0.;
+  for ( G4int i=0; i<fNv; ++i ) {
+    G4int j = i+1;
+    if ( j == fNv ) j = 0;
+    area += 0.5 * ( polygon[i].x()*polygon[j].y() - polygon[j].x()*polygon[i].y());
+  }
+ 
   // Copy polygon
   //
-  for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
+  if  ( area < 0. ) {   
+    // Polygon vertices are defined clockwise, we just copy the polygon       
+    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
+  }
+  else {
+    // Polygon vertices are defined anti-clockwise, we revert them
+    //G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
+    //            JustWarning, 
+    //            "Polygon vertices defined anti-clockwise, reverting polygon");      
+    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
+  }
+    
   
   // Copy z-sections
@@ -148 +170 @@
   }
      
+  // Check if polygon vertices are defined clockwise
+  // (the area is positive if polygon vertices are defined anti-clockwise)
+  
+  G4double area = 0.;
+  for ( G4int i=0; i<fNv; ++i ) {
+    G4int j = i+1;
+    if ( j == fNv ) j = 0;
+    area += 0.5 * ( polygon[i].x()*polygon[j].y() - polygon[j].x()*polygon[i].y());
+  }
+ 
   // Copy polygon
   //
-  for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
+  if  ( area < 0. ) {   
+    // Polygon vertices are defined clockwise, we just copy the polygon       
+    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
+  }
+  else {
+    // Polygon vertices are defined anti-clockwise, we revert them
+    //G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
+    //            JustWarning, 
+    //            "Polygon vertices defined anti-clockwise, reverting polygon");      
+    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
+  }
   
   // Copy z-sections
@@ -453 +495 @@
     //
     G4double angle = GetAngle(c2->first, c3->first, c1->first);
-
-    if ( angle > pi )
+    //G4cout << "angle " << angle  << G4endl;
+
+    G4int counter = 0;
+    while ( angle > pi )
     {
       // G4cout << "Skipping concave vertex " << c2->second << G4endl;
@@ -467 +511 @@
       // G4cout << "Looking at triangle : "
       //        << c1->second << "  " << c2->second
-      //        << "  " << c3->second << G4endl;  
-
+      //        << "  " << c3->second << G4endl; 
+      
+      angle = GetAngle(c2->first, c3->first, c1->first); 
+      //G4cout << "angle " << angle  << G4endl;
+      
+      counter++;
+      
+      if ( counter > fNv) {
+        G4Exception("G4ExtrudedSolid::AddGeneralPolygonFacets", "InvalidSetup" ,
+                    FatalException, "Triangularisation has failed.");
+        break;
+      }  
     }
 
@@ -734 +788 @@
 }
 
-
 //_____________________________________________________________________________
 
@@ -751 +804 @@
   for ( G4int i=0; i<fNv; ++i )
   {
-    os << "   vx = " << fPolygon[i].x()/mm << " mm" 
+    os << std::setw(5) << "#" << i 
+       << "   vx = " << fPolygon[i].x()/mm << " mm" 
        << "   vy = " << fPolygon[i].y()/mm << " mm" << G4endl;
   }
@@ -764 +818 @@
   }     
 
+/*
+  // Triangles (for debogging)
+  os << G4endl; 
+  os << " Triangles:" << G4endl;
+  os << " Triangle #   vertex1   vertex2   vertex3" << G4endl;
+
+  G4int counter = 0;
+  std::vector< std::vector<G4int> >::const_iterator it;
+  for ( it = fTriangles.begin(); it != fTriangles.end(); it++ ) {
+     std::vector<G4int> triangle = *it;
+     os << std::setw(10) << counter++ 
+        << std::setw(10) << triangle[0] << std::setw(10)  << triangle[1]  << std::setw(10)  << triangle[2] 
+        << G4endl;
+  }          
+*/
   return os;
 }  

trunk/source/geometry/solids/specific/src/G4PolyconeSide.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4PolyconeSide.cc,v 1.22 2009/11/11 12:23:37 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4PolyconeSide.cc,v 1.24 2010/02/24 11:31:49 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -67 +67 @@
   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
   fSurfaceArea = 0.0;
+  fPhi.first = G4ThreeVector(0,0,0);
+  fPhi.second= 0.0;
 
   //
@@ -492 +494 @@
   if (phiIsOpen)
   {
-    G4double phi = axis.phi();
+    G4double phi = GetPhi(axis);
     while( phi < startPhi ) phi += twopi;
     
@@ -853 +855 @@
 }
 
+//
+// GetPhi
+//
+// Calculate Phi for a given 3-vector (point), if not already cached for the
+// same point, in the attempt to avoid consecutive computation of the same
+// quantity
+//
+G4double G4PolyconeSide::GetPhi( const G4ThreeVector& p )
+{
+  G4double val=0.;
+
+  if (fPhi.first != p)
+  {
+    val = p.phi();
+    fPhi.first = p;
+    fPhi.second = val;
+  }
+  else
+  {
+    val = fPhi.second;
+  }
+  return val;
+}
 
 //
@@ -922 +947 @@
     // Finally, check phi
     //
-    G4double phi = p.phi();
+    G4double phi = GetPhi(p);
     while( phi < startPhi ) phi += twopi;
     
@@ -978 +1003 @@
     // PolyPhiFace. See PolyPhiFace::InsideEdgesExact
     //
-    G4double phi = hit.phi();
+    G4double phi = GetPhi(hit);
     while( phi < startPhi-phiTolerant ) phi += twopi;
     

trunk/source/geometry/solids/specific/src/G4PolyhedraSide.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4PolyhedraSide.cc,v 1.15 2008/05/15 11:41:59 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4PolyhedraSide.cc,v 1.17 2010/02/24 11:31:49 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -64 +64 @@
                                         G4bool isAllBehind )
 {
-
   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
   fSurfaceArea=0.;
+  fPhi.first = G4ThreeVector(0,0,0);
+  fPhi.second= 0.0;
+
   //
   // Record values
@@ -561 +563 @@
   // Try the closest phi segment first
   //
-  G4int iPhi = ClosestPhiSegment( p.phi() );
+  G4int iPhi = ClosestPhiSegment( GetPhi(p) );
   
   G4ThreeVector pdotc = p - vecs[iPhi].center;
@@ -594 +596 @@
   // Which phi segment is closest to this point?
   //
-  G4int iPhi = ClosestPhiSegment( p.phi() );
+  G4int iPhi = ClosestPhiSegment( GetPhi(p) );
   
   G4double norm;
@@ -624 +626 @@
   // Which phi segment is closest to this point?
   //
-  G4int iPhi = ClosestPhiSegment( p.phi() );
+  G4int iPhi = ClosestPhiSegment( GetPhi(p) );
 
   //
@@ -656 +658 @@
   // Which phi segment, if any, does the axis belong to
   //
-  iPhi = PhiSegment( axis.phi() );
+  iPhi = PhiSegment( GetPhi(axis) );
   
   if (iPhi < 0)
@@ -971 +973 @@
 
 //
+// GetPhi
+//
+// Calculate Phi for a given 3-vector (point), if not already cached for the
+// same point, in the attempt to avoid consecutive computation of the same
+// quantity
+//
+G4double G4PolyhedraSide::GetPhi( const G4ThreeVector& p )
+{
+  G4double val=0.;
+
+  if (fPhi.first != p)
+  {
+    val = p.phi();
+    fPhi.first = p;
+    fPhi.second = val;
+  }
+  else
+  {
+    val = fPhi.second;
+  }
+  return val;
+}
+
+
+//
 // DistanceToOneSide
 //

trunk/source/geometry/solids/specific/src/G4SolidExtentList.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4SolidExtentList.cc,v 1.5 2007/05/11 13:54:29 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4SolidExtentList.cc,v 1.6 2010/02/10 16:38:37 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -51 +51 @@
   axis = kZAxis;
   limited = false;
-  minLimit = -DBL_MAX;
-  maxLimit = +DBL_MAX;
+  minLimit = -kInfinity;
+  maxLimit = +kInfinity;
 }
 
@@ -72 +72 @@
   else
   {
-    minLimit = -DBL_MAX;
-    maxLimit = +DBL_MAX;
+    minLimit = -kInfinity;
+    maxLimit = +kInfinity;
   }
 }

trunk/source/geometry/solids/specific/src/G4TessellatedSolid.cc

@@ -25 +25 @@
 // ********************************************************************
 //
-// $Id: G4TessellatedSolid.cc,v 1.19 2009/04/27 08:06:27 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4TessellatedSolid.cc,v 1.20 2010/04/28 16:21:21 flei Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -42 +42 @@
 // CHANGE HISTORY
 // --------------
+//
+// 12 April 2010      P R Truscott, QinetiQ, bug fixes to treat optical
+//                    photon transport, in particular internal reflection
+//                    at surface.
 //
 // 14 November 2007   P R Truscott, QinetiQ & Stan Seibert, U Texas
@@ -670 +674 @@
     if ((*f)->Intersect(p,v,false,dist,distFromSurface,normal))
     {
+//
+//
+// Set minDist to the new distance to current facet if distFromSurface is in
+// positive direction and point is not at surface.  If the point is within
+// 0.5*kCarTolerance of the surface, then force distance to be zero and
+// leave member function immediately (for efficiency), as proposed by & credit
+// to Akira Okumura.
+//
       if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist)
       {
         minDist  = dist;
+      }
+      else if (-0.5*kCarTolerance <= dist && dist <= 0.5*kCarTolerance)
+      {
+        return 0.0;
       }
     }

trunk/source/geometry/solids/specific/src/G4TriangularFacet.cc

@@ -25 +25 @@
 // ********************************************************************
 //
-// $Id: G4TriangularFacet.cc,v 1.12 2008/11/13 08:25:07 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4TriangularFacet.cc,v 1.13 2010/04/28 16:21:21 flei Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -56 +56 @@
 //                  plane of the triangle.
 //
+// 12 April 2010    P R Truscott, QinetiQ, bug fixes to treat optical
+//                  photon transport, in particular internal reflection
+//                  at surface.
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -350 +353 @@
 //
 //
-// Do a heck for rounding errors in the distance-squared.
+// Do a check for rounding errors in the distance-squared.  It appears that
+// the conventional methods for calculating sqrDist breaks down when very near
+// to or at the surface (as required by transport).  We'll therefore also use
+// the magnitude-squared of the vector displacement.  (Note that I've also
+// tried to get around this problem by using the existing equations for
+//
+//    sqrDist = function(a,b,c,d,s,t)
+//
+// and use a more accurate addition process which minimises errors and
+// breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
+// doesn't work.
+// Calculation from u = D + s*E[0] + t*E[1] is less efficient, but appears
+// more robust.
 //
   if (sqrDist < 0.0) { sqrDist = 0.0; }
-
-  return D + s*E[0] + t*E[1];
+  G4ThreeVector u = D + s*E[0] + t*E[1];
+  G4double u2     = u.mag2();
+//
+//
+// The following (part of the roundoff error check) is from Oliver Merle's
+// updates.
+//
+  if ( sqrDist > u2 ) sqrDist = u2;
+
+  return u;
 }
 
@@ -542 +565 @@
 // we intersect.
 //
-      distance = -0.5*kCarTolerance;
+//      distance = -0.5*kCarTolerance;
+      distance = 0.0;
       normal   = surfaceNormal;
       return true;
@@ -728 +752 @@
   return area;
 }
+////////////////////////////////////////////////////////////////////////
+//