Skip to content
Commits on Source (47)
Expand all Hide whitespace changes
Inline Side-by-side
source/geometry/management/include/G4AffineTransform.hh
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4AffineTransform.hh 90701 2015-06-08 09:43:16Z gcosmo $
// $Id: G4AffineTransform.hh 109776 2018-05-09 08:00:11Z gcosmo $
//
//
// class G4AffineTransform
......@@ -53,11 +53,15 @@
// History:
// Paul R C Kent 6 Aug 1996 - initial version
//
// 19.09.96 E.Chernyaev:
// 19.09.1996 E.Tcherniaev:
// - direct access to the protected members of the G4RotationMatrix class
// replaced by access via public access functions
// replaced by access via public access functions
// - conversion of the rotation matrix to angle & axis used to get
// a possibility to remove "friend" from the G4RotationMatrix class
// 06.05.2018 E.Tcherniaev:
// - optimized InverseProduct
// - added methods for inverse transformation: InverseTrasformPoint,
// InverseTransformAxis, InverseNetRotation, InverseNetTranslation
// --------------------------------------------------------------------
#ifndef G4AFFINETRANSFORM_HH
#define G4AFFINETRANSFORM_HH
......@@ -65,99 +69,123 @@
#include "G4Types.hh"
#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"
#include "G4Transform3D.hh"
class G4AffineTransform
{
public:
G4AffineTransform();
inline G4AffineTransform();
public: // with description
G4AffineTransform(const G4ThreeVector &tlate);
inline G4AffineTransform(const G4ThreeVector& tlate);
// Translation only: under t'form translate point at origin by tlate
G4AffineTransform(const G4RotationMatrix &rot);
inline G4AffineTransform(const G4RotationMatrix& rot);
// Rotation only: under t'form rotate by rot
G4AffineTransform(const G4RotationMatrix &rot,
const G4ThreeVector &tlate);
inline G4AffineTransform(const G4RotationMatrix& rot,
const G4ThreeVector& tlate);
// Under t'form: rotate by rot then translate by tlate
G4AffineTransform(const G4RotationMatrix *rot,
const G4ThreeVector &tlate);
inline G4AffineTransform(const G4RotationMatrix* rot,
const G4ThreeVector& tlate);
// Optionally rotate by *rot then translate by tlate - rot may be null
G4AffineTransform operator * (const G4AffineTransform &tf) const;
inline G4AffineTransform(const G4AffineTransform& rhs);
// Copy constructor
inline G4AffineTransform& operator=(const G4AffineTransform& rhs);
// Assignment operator
inline ~G4AffineTransform();
// Destructor
inline G4AffineTransform operator * (const G4AffineTransform& tf) const;
// Compound Transforms:
// tf2=tf2*tf1 equivalent to tf2*=tf1
// Returns compound transformation of self*tf
G4AffineTransform& operator *= (const G4AffineTransform &tf);
inline G4AffineTransform& operator *= (const G4AffineTransform& tf);
// (Modifying) Multiplies self by tf; Returns self reference
// ie. A=AB for a*=b
G4AffineTransform& Product(const G4AffineTransform &tf1,
const G4AffineTransform &tf2);
inline G4AffineTransform& Product(const G4AffineTransform& tf1,
const G4AffineTransform& tf2);
// 'Products' for avoiding (potential) temporaries:
// c.Product(a,b) equivalent to c=a*b
// c.InverseProduct(a*b,b ) equivalent to c=a
// (Modifying) Sets self=tf1*tf2; Returns self reference
G4AffineTransform& InverseProduct(const G4AffineTransform &tf1,
const G4AffineTransform &tf2);
inline G4AffineTransform& InverseProduct(const G4AffineTransform& tf1,
const G4AffineTransform& tf2);
// (Modifying) Sets self=tf1*(tf2^-1); Returns self reference
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const;
inline G4ThreeVector TransformPoint(const G4ThreeVector& vec) const;
// Transform the specified point: returns vec*rot+tlate
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const;
// Transform the specified axis: returns
inline G4ThreeVector InverseTransformPoint(const G4ThreeVector& vec) const;
// Transform the specified point using inverse transformation
inline G4ThreeVector TransformAxis(const G4ThreeVector& axis) const;
// Transform the specified axis: returns vec*rot
void ApplyPointTransform(G4ThreeVector &vec) const;
inline G4ThreeVector InverseTransformAxis(const G4ThreeVector& axis) const;
// Transform the specified axis using inverse transfromation
inline void ApplyPointTransform(G4ThreeVector& vec) const;
// Transform the specified point (in place): sets vec=vec*rot+tlate
void ApplyAxisTransform(G4ThreeVector &axis) const;
inline void ApplyAxisTransform(G4ThreeVector& axis) const;
// Transform the specified axis (in place): sets axis=axis*rot;
G4AffineTransform Inverse() const;
inline G4AffineTransform Inverse() const;
// Return inverse of current transform
G4AffineTransform& Invert();
inline G4AffineTransform& Invert();
// (Modifying) Sets self=inverse of self; Returns self reference
G4AffineTransform& operator +=(const G4ThreeVector &tlate);
G4AffineTransform& operator -=(const G4ThreeVector &tlate);
inline G4AffineTransform& operator +=(const G4ThreeVector& tlate);
inline G4AffineTransform& operator -=(const G4ThreeVector& tlate);
// (Modifying) Adjust net translation by given vector;
// Returns self reference
G4bool operator == (const G4AffineTransform &tf) const;
G4bool operator != (const G4AffineTransform &tf) const;
inline G4bool operator == (const G4AffineTransform& tf) const;
inline G4bool operator != (const G4AffineTransform& tf) const;
G4double operator [] (const G4int n) const;
inline G4double operator [] (const G4int n) const;
G4bool IsRotated() const;
inline G4bool IsRotated() const;
// True if transform includes rotation
G4bool IsTranslated() const;
inline G4bool IsTranslated() const;
// True if transform includes translation
G4RotationMatrix NetRotation() const;
inline G4RotationMatrix NetRotation() const;
inline G4RotationMatrix InverseNetRotation() const;
inline G4ThreeVector NetTranslation() const;
inline G4ThreeVector InverseNetTranslation() const;
G4ThreeVector NetTranslation() const;
inline void SetNetRotation(const G4RotationMatrix& rot);
void SetNetRotation(const G4RotationMatrix &rot);
inline void SetNetTranslation(const G4ThreeVector& tlate);
void SetNetTranslation(const G4ThreeVector &tlate);
inline operator G4Transform3D () const;
// Conversion operator (cast) to G4Transform3D
private:
G4AffineTransform(const G4double prxx,const G4double prxy,const G4double prxz,
const G4double pryx,const G4double pryy,const G4double pryz,
const G4double przx,const G4double przy,const G4double przz,
const G4double ptx, const G4double pty, const G4double ptz );
inline G4AffineTransform(
const G4double prxx, const G4double prxy, const G4double prxz,
const G4double pryx, const G4double pryy, const G4double pryz,
const G4double przx, const G4double przy, const G4double przz,
const G4double ptx, const G4double pty, const G4double ptz);
G4double rxx,rxy,rxz;
G4double ryx,ryy,ryz;
......@@ -165,8 +193,7 @@ private:
G4double tx,ty,tz;
};
std::ostream&
operator << (std::ostream &os, const G4AffineTransform &transf);
std::ostream& operator << (std::ostream& os, const G4AffineTransform& transf);
#include "G4AffineTransform.icc"
......
source/geometry/management/include/G4BoundingEnvelope.hh 0 → 100644
View file @ 292425bd
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id:$
//
//
// class G4BoundingEnvelope
//
// Class description:
//
// Helper class to facilitate calculation of the extent of a solid
// within the limits defined by a G4VoxelLimits object.
//
// The function CalculateExtent() of a particular solid can create
// a G4BoundingEnvelope object that bounds the solid and then call
// CalculateExtent() of the G4BoundingEnvelope object.
//
// Calculation of extent uses G4Transform3D, thus takes into account
// scaling and reflection, if any.
// History:
//
// 2016.05.25 E.Tcherniaev - initial version
//
// --------------------------------------------------------------------
#ifndef G4BOUNDINGENVELOPE_HH
#define G4BOUNDINGENVELOPE_HH
#include <vector>
#include "geomdefs.hh"
#include "G4ThreeVector.hh"
#include "G4VoxelLimits.hh"
#include "G4Transform3D.hh"
#include "G4Point3D.hh"
#include "G4Plane3D.hh"
typedef std::vector<G4ThreeVector> G4ThreeVectorList;
typedef std::vector<G4Point3D> G4Polygon3D;
typedef std::pair<G4Point3D,G4Point3D> G4Segment3D;
class G4BoundingEnvelope
{
public:
G4BoundingEnvelope(const G4ThreeVector& pMin,
const G4ThreeVector& pMax);
// Constructor from an axis aligned bounding box (AABB)
G4BoundingEnvelope(const std::vector<const G4ThreeVectorList*>& polygons);
// Constructor from a sequence of convex polygons, the polygons
// should have equal numbers of vertices except first and last
// polygons which may consist of a single vertex
G4BoundingEnvelope(const G4ThreeVector& pMin,
const G4ThreeVector& pMax,
const std::vector<const G4ThreeVectorList*>& polygons);
// Constructor from AABB and a sequence of polygons
~G4BoundingEnvelope();
// Destructor
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis,
const G4VoxelLimits& pVoxelLimits,
const G4Transform3D& pTransform3D,
G4double& pMin, G4double& pMax) const;
// Analyse the position of the bounding box relative to the voxel.
// It returns "true" in the case where the value of the extent can be
// figured out directly from the dimensions of the bounding box, or
// it is clear that the bounding box and the voxel do not intersect.
// The reply "false" means that further calculations are needed.
G4bool CalculateExtent(const EAxis pAxis,
const G4VoxelLimits& pVoxelLimits,
const G4Transform3D& pTransform3D,
G4double& pMin, G4double& pMax) const;
// Calculate extent of the bounding envelope
private:
void CheckBoundingBox();
// Check correctness of the AABB (axis aligned bounding box)
void CheckBoundingPolygons();
// Check correctness of the sequence of convex polygonal bases
G4double FindScaleFactor(const G4Transform3D& pTransform3D) const;
// Find max scale factor of the transformation
void TransformVertices(const G4Transform3D& pTransform3D,
std::vector<G4Polygon3D*>& pBases) const;
// Create list of transformed polygons
void GetPrismAABB(const G4Polygon3D& pBaseA,
const G4Polygon3D& pBaseB,
G4Segment3D& pAABB) const;
// Find bounding box of a prism
void CreateListOfEdges(const G4Polygon3D& baseA,
const G4Polygon3D& baseB,
std::vector<G4Segment3D>& pEdges) const;
// Create list of edges of a prism
void CreateListOfPlanes(const G4Polygon3D& baseA,
const G4Polygon3D& baseB,
std::vector<G4Plane3D>& pPlanes) const;
// Create list of planes bounding a prism
G4bool ClipEdgesByVoxel(const std::vector<G4Segment3D>& pEdges,
const G4VoxelLimits& pLimits,
G4Segment3D& pExtent) const;
// Clip set of edges by G4VoxelLimits
void ClipVoxelByPlanes(G4int pBits,
const G4VoxelLimits& pLimits,
const std::vector<G4Plane3D>& pPlanes,
const G4Segment3D& pAABB,
G4Segment3D& pExtent) const;
// Clip G4VoxelLimits by set of planes bounding a prism
private:
G4ThreeVector fMin, fMax;
// original bounding box
const std::vector<const G4ThreeVectorList*>* fPolygons;
// ref to original sequence of polygonal bases
};
#endif // G4BOUNDINGENVELOPE_HH
source/geometry/management/include/G4GeomTools.hh 0 → 100644
View file @ 292425bd
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id: $
//
//
// --------------------------------------------------------------------
// class G4GeomTools
//
// Class description:
//
// A collection of utilities which can be helpfull for a wide range
// of geometry-related tasks
// History:
//
// 10.10.2016 E.Tcherniaev: Initial version.
// --------------------------------------------------------------------
#ifndef G4GEOMTOOLS_HH
#define G4GEOMTOOLS_HH
#include <vector>
#include "G4TwoVector.hh"
#include "G4ThreeVector.hh"
typedef std::vector<G4TwoVector> G4TwoVectorList;
typedef std::vector<G4ThreeVector> G4ThreeVectorList;
class G4GeomTools
{
public:
// ==================================================================
// 2D Utilities
// ------------------------------------------------------------------
static G4double TriangleArea(G4double Ax, G4double Ay,
G4double Bx, G4double By,
G4double Cx, G4double Cy);
static G4double TriangleArea(const G4TwoVector& A,
const G4TwoVector& B,
const G4TwoVector& C);
// Calcuate area of 2D triangle, return value is positive if
// vertices of the triangle are given in anticlockwise order,
// otherwise it is negative
static G4double QuadArea(const G4TwoVector& A,
const G4TwoVector& B,
const G4TwoVector& C,
const G4TwoVector& D);
// Calcuate area of 2D quadrilateral, return value is positive if
// vertices of the quadrilateral are given in anticlockwise order,
// otherwise it is negative
static G4double PolygonArea(const G4TwoVectorList& polygon);
// Calcuate area of 2D polygon, return value is positive if
// vertices of the polygon are defined in anticlockwise order,
// otherwise it is negative
static G4bool PointInTriangle(G4double Px, G4double Py,
G4double Ax, G4double Ay,
G4double Bx, G4double By,
G4double Cx, G4double Cy);
// Decide if point (Px,Py) is inside triangle (Ax,Ay)(Bx,By)(Cx,Cy)
static G4bool PointInTriangle(const G4TwoVector& P,
const G4TwoVector& A,
const G4TwoVector& B,
const G4TwoVector& C);
// Decide if point P is inside triangle ABC
/*
static G4bool PointInPolygon(const G4TwoVector& P
const G4TwoVectorList& Polygon);
// Decide if point P is inside Polygon
*/
static G4bool IsConvex(const G4TwoVectorList& polygon);
// Decide if 2D polygon is convex, i.e. all internal angles are
// less than pi
static G4bool TriangulatePolygon(const G4TwoVectorList& polygon,
G4TwoVectorList& result);
static G4bool TriangulatePolygon(const G4TwoVectorList& polygon,
std::vector<G4int>& result);
// Simple implementation of "ear clipping" algorithm for
// triangulation of a simple contour/polygon, it places results
// in a std::vector as triplets of vertices. If triangulation
// is sucsessfull then the function returns true, otherwise false
static void RemoveRedundantVertices(G4TwoVectorList& polygon,
std::vector<G4int>& iout,
G4double tolerance = 0);
// Remove collinear and coincident points from 2D polygon.
// Indices of removed points are available in iout.
static G4bool DiskExtent(G4double rmin, G4double rmax,
G4double startPhi, G4double delPhi,
G4TwoVector& pmin, G4TwoVector& pmax);
// Calculate bounding rectangle of a disk sector,
// it returns false if input parameters do not meet the following:
// rmin >= 0
// rmax > rmin + kCarTolerance
// delPhi > 0 + kCarTolerance
static void DiskExtent(G4double rmin, G4double rmax,
G4double sinPhiStart, G4double cosPhiStart,
G4double sinPhiEnd, G4double cosPhiEnd,
G4TwoVector& pmin, G4TwoVector& pmax);
// Calculate bounding rectangle of a disk sector,
// faster version without check of parameters
static G4double EllipsePerimeter(G4double a,
G4double b);
// Compute the circumference (perimeter) of an ellipse
static G4double EllipticConeLateralArea(G4double a,
G4double b,
G4double h);
// Compute the lateral surface area of an elliptic cone
// ==================================================================
// 3D Utilities
// ------------------------------------------------------------------
static G4ThreeVector TriangleAreaNormal(const G4ThreeVector& A,
const G4ThreeVector& B,
const G4ThreeVector& C);
// Find the normal to the plane of 3D triangle ABC,
// length of the normal is equal to the area of the triangle
static G4ThreeVector QuadAreaNormal(const G4ThreeVector& A,
const G4ThreeVector& B,
const G4ThreeVector& C,
const G4ThreeVector& D);
// Find normal to the plane of 3D quadrilateral ABCD,
// length of the normal is equal to the area of the quadrilateral
static G4ThreeVector PolygonAreaNormal(const G4ThreeVectorList& polygon);
// Find normal to the plane of 3D polygon
// length of the normal is equal to the area of the polygon
/*
static G4bool IsPlanar(const G4ThreeVector& A,
const G4ThreeVector& B,
const G4ThreeVector& C,
const G4ThreeVector& D);
// Decide if 3D quadrilateral ABCD is planar
static G4bool IsPlanar(const G4ThreeVectorList& polygon
const G4ThreeVector& normal);
// Decide if 3D polygon is planar
*/
static G4double DistancePointSegment(const G4ThreeVector& P,
const G4ThreeVector& A,
const G4ThreeVector& B);
// Calculate distance between point P and line segment AB in 3D
static G4ThreeVector ClosestPointOnSegment(const G4ThreeVector& P,
const G4ThreeVector& A,
const G4ThreeVector& B);
// Find point on 3D line segment AB closest to point P
static G4ThreeVector ClosestPointOnTriangle(const G4ThreeVector& P,
const G4ThreeVector& A,
const G4ThreeVector& B,
const G4ThreeVector& C);
// Find point on 3D triangle ABC closest to point P
static G4bool SphereExtent(G4double rmin, G4double rmax,
G4double startTheta, G4double delTheta,
G4double startPhi, G4double delPhi,
G4ThreeVector& pmin, G4ThreeVector& pmax);
// Calculate bounding box of a spherical sector,
// it returns false if input parameters do not meet the following:
// rmin >= 0
// rmax > rmin + kCarTolerance
// startTheta >= 0 && <= pi;
// delTheta > 0 + kCarTolerance
// delPhi > 0 + kCarTolerance
private:
static G4bool CheckSnip(const G4TwoVectorList& contour,
G4int a, G4int b, G4int c,
G4int n, const G4int* V);
// Helper function for TriangulatePolygon()
static G4double comp_ellint_2(G4double e);
// Complete Elliptic Integral of the Second Kind
};
#endif // G4GEOMTOOLS_HH
source/geometry/management/sources.cmake
View file @ 292425bd
......@@ -37,12 +37,14 @@ GEANT4_DEFINE_MODULE(NAME G4geometrymng
G4AffineTransform.icc
G4BlockingList.hh
G4BlockingList.icc
G4BoundingEnvelope.hh
G4ErrorCylSurfaceTarget.hh
G4ErrorPlaneSurfaceTarget.hh
G4ErrorSurfaceTarget.hh
G4ErrorTanPlaneTarget.hh
G4ErrorTarget.hh
G4GeomSplitter.hh
G4GeomTools.hh
G4GeometryManager.hh
G4IdentityTrajectoryFilter.hh
G4LogicalSurface.hh
......@@ -85,11 +87,13 @@ GEANT4_DEFINE_MODULE(NAME G4geometrymng
voxeldefs.hh
SOURCES
G4BlockingList.cc
G4BoundingEnvelope.cc
G4ErrorCylSurfaceTarget.cc
G4ErrorPlaneSurfaceTarget.cc
G4ErrorSurfaceTarget.cc
G4ErrorTanPlaneTarget.cc
G4ErrorTarget.cc
G4GeomTools.cc
G4GeometryManager.cc
G4IdentityTrajectoryFilter.cc
G4LogicalSurface.cc
......
source/geometry/management/src/G4BoundingEnvelope.cc 0 → 100644
View file @ 292425bd
This diff is collapsed.
source/geometry/management/src/G4GeomTools.cc 0 → 100644
View file @ 292425bd
This diff is collapsed.
source/geometry/navigation/include/G4Navigator.hh
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4Navigator.hh 81577 2014-06-03 10:13:36Z gcosmo $
// $Id: G4Navigator.hh 96456 2016-04-15 10:14:01Z gcosmo $
//
//
// class G4Navigator
......@@ -358,8 +358,8 @@ class G4Navigator
protected: // without description
G4double kCarTolerance;
// Geometrical tolerance for surface thickness of shapes.
G4double kCarTolerance, fMinStep, fSqTol;
// Cached tolerances.
//
// BEGIN State information
......
source/geometry/navigation/include/G4Navigator.icc
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4Navigator.icc 73164 2013-08-21 08:55:28Z gcosmo $
// $Id: G4Navigator.icc 109960 2018-05-11 12:54:15Z evc $
//
//
// class G4Navigator Inline implementation
......@@ -52,8 +52,7 @@ G4ThreeVector G4Navigator::GetCurrentLocalCoordinate() const
inline
G4ThreeVector G4Navigator::ComputeLocalAxis(const G4ThreeVector& pVec) const
{
return (fHistory.GetTopTransform().IsRotated())
? fHistory.GetTopTransform().TransformAxis(pVec) : pVec ;
return fHistory.GetTopTransform().TransformAxis(pVec);
}
// ********************************************************************
......@@ -193,9 +192,7 @@ const G4AffineTransform& G4Navigator::GetGlobalToLocalTransform() const
inline
const G4AffineTransform G4Navigator::GetLocalToGlobalTransform() const
{
G4AffineTransform tempTransform;
tempTransform = fHistory.GetTopTransform().Inverse();
return tempTransform;
return fHistory.GetTopTransform().Inverse();
}
// ********************************************************************
......@@ -207,8 +204,7 @@ const G4AffineTransform G4Navigator::GetLocalToGlobalTransform() const
inline
G4ThreeVector G4Navigator::NetTranslation() const
{
G4AffineTransform tf(fHistory.GetTopTransform().Inverse());
return tf.NetTranslation();
return fHistory.GetTopTransform().InverseNetTranslation();
}
// ********************************************************************
......@@ -220,8 +216,7 @@ G4ThreeVector G4Navigator::NetTranslation() const
inline
G4RotationMatrix G4Navigator::NetRotation() const
{
G4AffineTransform tf(fHistory.GetTopTransform().Inverse());
return tf.NetRotation();
return fHistory.GetTopTransform().InverseNetRotation();
}
// ********************************************************************
......@@ -233,10 +228,10 @@ G4RotationMatrix G4Navigator::NetRotation() const
inline
G4GRSVolume* G4Navigator::CreateGRSVolume() const
{
G4AffineTransform tf(fHistory.GetTopTransform().Inverse());
const G4AffineTransform& tf = fHistory.GetTopTransform();
return new G4GRSVolume(fHistory.GetTopVolume(),
tf.NetRotation(),
tf.NetTranslation());
tf.InverseNetRotation(),
tf.InverseNetTranslation());
}
// ********************************************************************
......@@ -248,10 +243,10 @@ G4GRSVolume* G4Navigator::CreateGRSVolume() const
inline
G4GRSSolid* G4Navigator::CreateGRSSolid() const
{
G4AffineTransform tf(fHistory.GetTopTransform().Inverse());
const G4AffineTransform& tf = fHistory.GetTopTransform();
return new G4GRSSolid(fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid(),
tf.NetRotation(),
tf.NetTranslation());
tf.InverseNetRotation(),
tf.InverseNetTranslation());
}
// ********************************************************************
......
source/geometry/navigation/src/G4Navigator.cc
View file @ 292425bd
......@@ -45,6 +45,12 @@
#include "G4VoxelSafety.hh"
// Constant determining how precise normals should be (how close to unit
// vectors). If exceeded, warnings will be issued.
// Can be CLHEP::perMillion (its old default) for geometry checking.
//
static const G4double kToleranceNormalCheck = CLHEP::perThousand;
// ********************************************************************
// Constructor
// ********************************************************************
......@@ -67,6 +73,9 @@ G4Navigator::G4Navigator()
fAbandonThreshold_NoZeroSteps = 25;
kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
fMinStep = 0.05*kCarTolerance;
fSqTol = kCarTolerance*kCarTolerance;
fregularNav.SetNormalNavigation( &fnormalNav );
fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity );
......@@ -283,7 +292,7 @@ G4Navigator::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint,
// o containing volume found
//
while (notKnownContained)
while (notKnownContained) // Loop checking, 07.10.2016, J.Apostolakis
{
if ( fHistory.GetTopVolumeType()!=kReplica )
{
......@@ -535,7 +544,7 @@ G4Navigator::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint,
}
#endif
}
} while (noResult);
} while (noResult); // Loop checking, 07.10.2016, J.Apostolakis
fLastLocatedPointLocal = localPoint;
......@@ -665,7 +674,7 @@ void G4Navigator::SetSavedState()
fSaveState.sEntering = fEntering;
fSaveState.spBlockedPhysicalVolume = fBlockedPhysicalVolume;
fSaveState.sBlockedReplicaNo = fBlockedReplicaNo,
fSaveState.sBlockedReplicaNo = fBlockedReplicaNo;
fSaveState.sLastStepWasZero = fLastStepWasZero;
......@@ -695,7 +704,7 @@ void G4Navigator::RestoreSavedState()
fEntering = fSaveState.sEntering;
fBlockedPhysicalVolume = fSaveState.spBlockedPhysicalVolume;
fBlockedReplicaNo = fSaveState.sBlockedReplicaNo,
fBlockedReplicaNo = fSaveState.sBlockedReplicaNo;
fLastStepWasZero = fSaveState.sLastStepWasZero;
......@@ -761,7 +770,6 @@ G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
fCalculatedExitNormal = false;
// Reset for new step
static const G4double minStep = 0.05*kCarTolerance;
static G4ThreadLocal G4int sNavCScalls=0;
sNavCScalls++;
......@@ -798,7 +806,7 @@ G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
if ( moveLenSq >= kCarTolerance*kCarTolerance )
if ( moveLenSq >= fSqTol )
{
#ifdef G4VERBOSE
ComputeStepLog(pGlobalpoint, moveLenSq);
......@@ -972,7 +980,7 @@ G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
// checked
//
fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
fLastStepWasZero = (Step<minStep);
fLastStepWasZero = (Step<fMinStep);
if (fPushed) { fPushed = fLastStepWasZero; }
// Handle large number of consecutive zero steps
......@@ -1001,12 +1009,15 @@ G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
if ((!fPushed) && (fWarnPush))
{
std::ostringstream message;
message.precision(16);
message << "Track stuck or not moving." << G4endl
<< " Track stuck, not moving for "
<< fNumberZeroSteps << " steps" << G4endl
<< " in volume -" << motherPhysical->GetName()
<< "- at point " << pGlobalpoint << G4endl
<< " direction: " << pDirection << "." << G4endl
<< "- at point " << pGlobalpoint
<< " (local point " << newLocalPoint << ")" << G4endl
<< " direction: " << pDirection
<< " (local direction: " << localDirection << ")." << G4endl
<< " Potential geometry or navigation problem !"
<< G4endl
<< " Trying pushing it of " << Step << " mm ...";
......@@ -1029,7 +1040,10 @@ G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
<< "- at point " << pGlobalpoint << G4endl
<< " direction: " << pDirection << ".";
#ifdef G4VERBOSE
motherPhysical->CheckOverlaps(5000, false);
if ( fWarnPush )
{
motherPhysical->CheckOverlaps(5000, false);
}
#endif
G4Exception("G4Navigator::ComputeStep()", "GeomNav0003",
EventMustBeAborted, message);
......@@ -1135,10 +1149,8 @@ G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
G4int depth= fHistory.GetDepth();
if( depth > 0 )
{
G4AffineTransform GrandMotherToGlobalTransf =
fHistory.GetTransform(depth-1).Inverse();
fExitNormalGlobalFrame =
GrandMotherToGlobalTransf.TransformAxis( fGrandMotherExitNormal );
fHistory.GetTransform(depth-1).InverseTransformAxis( fGrandMotherExitNormal );
}
else
{
......@@ -1324,6 +1336,7 @@ G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
G4ThreeVector ExitNormal(0.,0.,0.);
G4VSolid *currentSolid=0;
G4LogicalVolume *candidateLogical;
if ( fLastTriedStepComputation )
{
// use fLastLocatedPointLocal and next candidate volume
......@@ -1347,7 +1360,7 @@ G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
GetMotherToDaughterTransform( fBlockedPhysicalVolume,
fBlockedReplicaNo,
VolumeType(fBlockedPhysicalVolume) );
G4ThreeVector daughterPointOwnLocal=
G4ThreeVector daughterPointOwnLocal=
MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal );
// OK if it is a parameterised volume
......@@ -1379,7 +1392,9 @@ G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
// Entering the solid ==> opposite
//
ExitNormal = -nextSolidExitNormal;
// First flip ( ExitNormal = -nextSolidExitNormal; )
// and then rotate the the normal to the frame of the mother (current volume)
ExitNormal = MotherToDaughterTransform.InverseTransformAxis( -nextSolidExitNormal );
fCalculatedExitNormal= true;
}
else
......@@ -1439,7 +1454,7 @@ G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
G4VSolid* daughterSolid =fHistory.GetTopVolume()->GetLogicalVolume()
->GetSolid();
ExitNormal= -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal));
if( std::fabs(ExitNormal.mag2()-1.0 ) > CLHEP::perMillion )
if( std::fabs(ExitNormal.mag2()-1.0 ) > kToleranceNormalCheck )
{
G4ExceptionDescription desc;
desc << " Parameters of solid: " << *daughterSolid
......@@ -1569,18 +1584,14 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
{
G4bool validNormal;
G4ThreeVector localNormal, globalNormal;
G4bool usingStored;
usingStored=
fCalculatedExitNormal &&
( ( fLastTriedStepComputation && fExiting) // Just calculated it
|| // No locate in between
(!fLastTriedStepComputation
&& (IntersectPointGlobal-fStepEndPoint).mag2()
< (10.0*kCarTolerance*kCarTolerance)
) // Calculated it 'just' before & then called locate
G4bool usingStored = fCalculatedExitNormal && (
( fLastTriedStepComputation && fExiting ) // Just calculated it
|| // No locate in between
( !fLastTriedStepComputation
&& (IntersectPointGlobal-fStepEndPoint).mag2() < 10.0*fSqTol ) );
// Calculated it 'just' before & then called locate
// but it did not move position
);
if( usingStored )
{
......@@ -1589,7 +1600,7 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
//
globalNormal = fExitNormalGlobalFrame;
G4double normMag2 = globalNormal.mag2();
if( std::fabs ( normMag2 - 1.0 ) < perMillion ) // Value is good
if( std::fabs ( normMag2 - 1.0 ) < perThousand ) // was perMillion ) // Value is good
{
*pNormalCalculated = true; // ComputeStep always computes it if Exiting
// (fExiting==true)
......@@ -1603,7 +1614,18 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
<< " - but |normal| = " << std::sqrt(normMag2)
<< " - and |normal|^2 = " << normMag2 << G4endl
<< " which differs from 1.0 by " << normMag2 - 1.0 << G4endl
<< " n = " << fExitNormalGlobalFrame << G4endl;
<< " n = " << fExitNormalGlobalFrame << G4endl
<< " Global point: " << IntersectPointGlobal << G4endl
<< " Volume: " << fHistory.GetTopVolume()->GetName() << G4endl;
#ifdef G4VERBOSE
G4LogicalVolume* candLog = fHistory.GetTopVolume()->GetLogicalVolume();
if ( candLog )
{
message << " Solid: " << candLog->GetSolid()->GetName()
<< ", Type: " << candLog->GetSolid()->GetEntityType() << G4endl
<< *candLog->GetSolid() << G4endl;
}
#endif
message << "============================================================"
<< G4endl;
G4int oldVerbose = fVerbose;
......@@ -1623,9 +1645,7 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,
&validNormal);
*pNormalCalculated = fCalculatedExitNormal;
G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
globalNormal = localToGlobal.TransformAxis( localNormal );
globalNormal = fHistory.GetTopTransform().InverseTransformAxis( localNormal );
}
}
else
......@@ -1636,7 +1656,7 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
#ifdef G4DEBUG_NAVIGATION
usingStored= false;
if( (!validNormal) && !fCalculatedExitNormal)
if( (!validNormal) && !fCalculatedExitNormal )
{
G4ExceptionDescription edN;
edN << " Calculated = " << fCalculatedExitNormal << G4endl;
......@@ -1654,26 +1674,35 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
#endif
G4double localMag2= localNormal.mag2();
if( validNormal && (std::fabs(localMag2-1.0)) > CLHEP::perMillion )
if( validNormal && (std::fabs(localMag2-1.0)) > kToleranceNormalCheck )
{
G4ExceptionDescription edN;
edN.precision(10);
edN << "G4Navigator::GetGlobalExitNormal: "
<< " Using Local Normal - from call to GetLocalExitNormalAndCheck. "
<< G4endl
<< " Local Exit Normal : " << " || = " << std::sqrt(localMag2)
<< " vec = " << localNormal << G4endl
<< " Global Exit Normal : " << " || = " << globalNormal.mag()
<< " vec = " << globalNormal << G4endl;
edN << " Calculated It = " << fCalculatedExitNormal << G4endl;
<< " vec = " << globalNormal << G4endl
<< " Global point: " << IntersectPointGlobal << G4endl;
edN << " Calculated It = " << fCalculatedExitNormal << G4endl
<< " Volume: " << fHistory.GetTopVolume()->GetName() << G4endl;
#ifdef G4VERBOSE
G4LogicalVolume* candLog = fHistory.GetTopVolume()->GetLogicalVolume();
if ( candLog )
{
edN << " Solid: " << candLog->GetSolid()->GetName()
<< ", Type: " << candLog->GetSolid()->GetEntityType() << G4endl
<< *candLog->GetSolid();
}
#endif
G4Exception("G4Navigator::GetGlobalExitNormal()",
"GeomNav0003",JustWarning, edN,
"Value obtained from new local *solid* is incorrect.");
localNormal = localNormal.unit(); // Should we correct it ??
}
G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
globalNormal = localToGlobal.TransformAxis( localNormal );
globalNormal = fHistory.GetTopTransform().InverseTransformAxis( localNormal );
}
#ifdef G4DEBUG_NAVIGATION
......@@ -1683,12 +1712,11 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
localNormal= GetLocalExitNormalAndCheck(IntersectPointGlobal, &validNormal);
G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
globalNormAgn = localToGlobal.TransformAxis( localNormal );
globalNormAgn = fHistory.GetTopTransform().InverseTransformAxis( localNormal );
// Check the value computed against fExitNormalGlobalFrame
G4ThreeVector diffNorm = globalNormAgn - fExitNormalGlobalFrame;
if( diffNorm.mag2() > perMillion*CLHEP::perMillion)
if( diffNorm.mag2() > kToleranceNormalCheck )
{
G4ExceptionDescription edDfn;
edDfn << "Found difference in normals in case of exiting mother "
......@@ -1702,7 +1730,11 @@ G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
}
}
#endif
// Synchronise stored global exit normal as possibly re-computed here
//
fExitNormalGlobalFrame = globalNormal;
return globalNormal;
}
......@@ -1910,7 +1942,7 @@ G4bool G4Navigator::RecheckDistanceToCurrentBoundary(
if( fEnteredDaughter )
{
if( motherLogical->CharacteriseDaughters() ==kReplica ) { return false; }
if( motherLogical->CharacteriseDaughters() == kReplica ) { return false; }
// Track arrived at boundary of a daughter volume at
// the last call of ComputeStep().
......@@ -2095,12 +2127,12 @@ void G4Navigator::PrintState() const
<< std::setw( 9) << fExiting << " "
<< std::setw( 9) << fEntering << " ";
if ( fBlockedPhysicalVolume==0 )
G4cout << std::setw(15) << "None";
{ G4cout << std::setw(15) << "None"; }
else
G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
G4cout << std::setw( 9) << fBlockedReplicaNo << " "
<< std::setw( 8) << fLastStepWasZero << " "
<< G4endl;
{ G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName(); }
G4cout << std::setw( 9) << fBlockedReplicaNo << " "
<< std::setw( 8) << fLastStepWasZero << " "
<< G4endl;
}
if( fVerbose > 2 )
{
......@@ -2122,11 +2154,11 @@ void G4Navigator::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
// The following checks only make sense if the move is larger
// than the tolerance.
static const G4double fAccuracyForWarning = kCarTolerance,
fAccuracyForException = 1000*kCarTolerance;
const G4double fAccuracyForWarning = kCarTolerance,
fAccuracyForException = 1000*kCarTolerance;
G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
TransformPoint(fLastLocatedPointLocal);
G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().
InverseTransformPoint(fLastLocatedPointLocal);
G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
......
source/geometry/solids/Boolean/include/G4DisplacedSolid.hh
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4DisplacedSolid.hh 83572 2014-09-01 15:23:27Z gcosmo $
// $Id: G4DisplacedSolid.hh 104311 2017-05-24 12:59:46Z gcosmo $
//
//
// class G4DisplacedSolid
......@@ -163,6 +163,9 @@ class G4DisplacedSolid : public G4VSolid
G4AffineTransform* fDirectTransform ;
mutable G4bool fRebuildPolyhedron;
mutable G4Polyhedron* fpPolyhedron;
G4ThreeVector fBoxCenter; // center of the bounding box
G4ThreeVector fBoxSize; // half-dimensions of the bounding box
} ;
#endif
source/geometry/solids/Boolean/include/G4SubtractionSolid.hh
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4SubtractionSolid.hh 66356 2012-12-18 09:02:32Z gcosmo $
// $Id: G4SubtractionSolid.hh 104311 2017-05-24 12:59:46Z gcosmo $
//
//
// class G4SubtractionSolid
......
source/geometry/solids/Boolean/include/G4UnionSolid.hh
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4UnionSolid.hh 66356 2012-12-18 09:02:32Z gcosmo $
// $Id: G4UnionSolid.hh 109443 2018-04-24 06:28:58Z evc $
//
//
// class G4UnionSolid
......@@ -115,6 +115,10 @@ class G4UnionSolid : public G4BooleanSolid
void DescribeYourselfTo ( G4VGraphicsScene& scene ) const ;
G4Polyhedron* CreatePolyhedron () const ;
private:
G4ThreeVector fBoxCenter; // center of the bounding box
G4ThreeVector fBoxSize; // half-dimensions of the bounding box
};
#endif
source/geometry/solids/Boolean/src/G4DisplacedSolid.cc
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4DisplacedSolid.cc 84211 2014-10-10 14:47:30Z gcosmo $
// $Id: G4DisplacedSolid.cc 108785 2018-03-07 08:37:43Z gcosmo $
//
// Implementation for G4DisplacedSolid class for boolean
// operations between other solids
......@@ -34,6 +34,7 @@
// 28.10.98 V.Grichine: created
// 14.11.99 V.Grichine: modifications in CalculateExtent(...) method
// 22.11.00 V.Grichine: new set methods for matrix/vectors
// 28.02.18 E.Tcherniaev: improved contruction from G4DisplacedSolid
//
// --------------------------------------------------------------------
......@@ -57,10 +58,32 @@ G4DisplacedSolid::G4DisplacedSolid( const G4String& pName,
const G4ThreeVector& transVector )
: G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0)
{
fPtrSolid = pSolid ;
fPtrTransform = new G4AffineTransform(rotMatrix,transVector) ;
fPtrTransform->Invert() ;
fDirectTransform = new G4AffineTransform(rotMatrix,transVector) ;
if (pSolid->GetEntityType() == "G4DisplacedSolid")
{
fPtrSolid = ((G4DisplacedSolid*)pSolid)->GetConstituentMovedSolid();
G4AffineTransform t1 = ((G4DisplacedSolid*)pSolid)->GetDirectTransform();
G4AffineTransform t2 = G4AffineTransform(rotMatrix,transVector);
fDirectTransform = new G4AffineTransform(t1*t2);
}
else
{
fPtrSolid = pSolid;
fDirectTransform = new G4AffineTransform(rotMatrix,transVector);
}
fPtrTransform = new G4AffineTransform(fDirectTransform->Inverse());
G4VoxelLimits voxelLimits;
G4AffineTransform affineTransform;
G4double vmin, vmax;
CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setX((vmax + vmin)/2.);
fBoxSize.setX((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setY((vmax + vmin)/2.);
fBoxSize.setY((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setZ((vmax + vmin)/2.);
fBoxSize.setZ((vmax - vmin + kCarTolerance)/2.);
}
/////////////////////////////////////////////////////////////////////////////////
......@@ -72,13 +95,34 @@ G4DisplacedSolid::G4DisplacedSolid( const G4String& pName,
const G4Transform3D& transform )
: G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0)
{
fPtrSolid = pSolid ;
fDirectTransform = new G4AffineTransform(transform.getRotation().inverse(),
transform.getTranslation()) ;
fPtrTransform = new G4AffineTransform(transform.getRotation().inverse(),
transform.getTranslation()) ;
fPtrTransform->Invert() ;
if (pSolid->GetEntityType() == "G4DisplacedSolid")
{
fPtrSolid = ((G4DisplacedSolid*)pSolid)->GetConstituentMovedSolid();
G4AffineTransform t1 = ((G4DisplacedSolid*)pSolid)->GetDirectTransform();
G4AffineTransform t2 = G4AffineTransform(transform.getRotation().inverse(),
transform.getTranslation());
fDirectTransform = new G4AffineTransform(t1*t2);
}
else
{
fPtrSolid = pSolid;
fDirectTransform = new G4AffineTransform(transform.getRotation().inverse(),
transform.getTranslation()) ;
}
fPtrTransform = new G4AffineTransform(fDirectTransform->Inverse());
G4VoxelLimits voxelLimits;
G4AffineTransform affineTransform;
G4double vmin, vmax;
CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setX((vmax + vmin)/2.);
fBoxSize.setX((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setY((vmax + vmin)/2.);
fBoxSize.setY((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setZ((vmax + vmin)/2.);
fBoxSize.setZ((vmax - vmin + kCarTolerance)/2.);
}
///////////////////////////////////////////////////////////////////
......@@ -91,9 +135,32 @@ G4DisplacedSolid::G4DisplacedSolid( const G4String& pName,
const G4AffineTransform directTransform )
: G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0)
{
fPtrSolid = pSolid ;
fDirectTransform = new G4AffineTransform( directTransform );
fPtrTransform = new G4AffineTransform( directTransform.Inverse() ) ;
if (pSolid->GetEntityType() == "G4DisplacedSolid")
{
fPtrSolid = ((G4DisplacedSolid*)pSolid)->GetConstituentMovedSolid();
G4AffineTransform t1 = ((G4DisplacedSolid*)pSolid)->GetDirectTransform();
G4AffineTransform t2 = G4AffineTransform(directTransform);
fDirectTransform = new G4AffineTransform(t1*t2);
}
else
{
fPtrSolid = pSolid;
fDirectTransform = new G4AffineTransform(directTransform);
}
fPtrTransform = new G4AffineTransform(fDirectTransform->Inverse());
G4VoxelLimits voxelLimits;
G4AffineTransform affineTransform;
G4double vmin, vmax;
CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setX((vmax + vmin)/2.);
fBoxSize.setX((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setY((vmax + vmin)/2.);
fBoxSize.setY((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setZ((vmax + vmin)/2.);
fBoxSize.setZ((vmax - vmin + kCarTolerance)/2.);
}
///////////////////////////////////////////////////////////////////
......@@ -127,6 +194,9 @@ G4DisplacedSolid::G4DisplacedSolid(const G4DisplacedSolid& rhs)
{
fPtrTransform = new G4AffineTransform(*(rhs.fPtrTransform));
fDirectTransform = new G4AffineTransform(*(rhs.fDirectTransform));
fBoxCenter = rhs.fBoxCenter;
fBoxSize = rhs.fBoxSize;
}
///////////////////////////////////////////////////////////////
......@@ -152,6 +222,8 @@ G4DisplacedSolid& G4DisplacedSolid::operator = (const G4DisplacedSolid& rhs)
fRebuildPolyhedron = false;
delete fpPolyhedron; fpPolyhedron= 0;
fBoxCenter = rhs.fBoxCenter;
fBoxSize = rhs.fBoxSize;
return *this;
}
......@@ -261,9 +333,9 @@ void G4DisplacedSolid::SetObjectTranslation(const G4ThreeVector& vector)
fRebuildPolyhedron = true;
}
///////////////////////////////////////////////////////////////
//
//////////////////////////////////////////////////////////////////////////
//
// Calculate extent under transform and specified limit
G4bool
G4DisplacedSolid::CalculateExtent( const EAxis pAxis,
......@@ -283,6 +355,12 @@ G4DisplacedSolid::CalculateExtent( const EAxis pAxis,
EInside G4DisplacedSolid::Inside(const G4ThreeVector& p) const
{
// Check if point is outside of bounding box
//
if (std::abs(p.z()-fBoxCenter.z()) > fBoxSize.z()) return kOutside;
//if (std::abs(p.y()-fBoxCenter.y()) > fBoxSize.y()) return kOutside;
//if (std::abs(p.x()-fBoxCenter.x()) > fBoxSize.x()) return kOutside;
G4ThreeVector newPoint = fPtrTransform->TransformPoint(p) ;
return fPtrSolid->Inside(newPoint) ;
}
......@@ -447,8 +525,18 @@ G4Polyhedron*
G4DisplacedSolid::CreatePolyhedron () const
{
G4Polyhedron* polyhedron = fPtrSolid->CreatePolyhedron();
polyhedron
if (polyhedron)
{
polyhedron
->Transform(G4Transform3D(GetObjectRotation(),GetObjectTranslation()));
}
else
{
DumpInfo();
G4Exception("G4DisplacedSolid::CreatePolyhedron()",
"GeomSolids2002", JustWarning,
"No G4Polyhedron for displaced solid");
}
return polyhedron;
}
......
source/geometry/solids/Boolean/src/G4SubtractionSolid.cc
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4SubtractionSolid.cc 95297 2016-02-04 09:30:14Z gcosmo $
// $Id: G4SubtractionSolid.cc 104311 2017-05-24 12:59:46Z gcosmo $
//
// Implementation of methods for the class G4IntersectionSolid
//
......@@ -134,9 +134,9 @@ G4SubtractionSolid::operator = (const G4SubtractionSolid& rhs)
return *this;
}
///////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//
// CalculateExtent
// Calculate extent under transform and specified limit
G4bool
G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
......@@ -146,7 +146,7 @@ G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
G4double& pMax ) const
{
// Since we cannot be sure how much the second solid subtracts
// from the first, we must use the first solid's extent!
// from the first, we must use the first solid's extent!
return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
pTransform, pMin, pMax );
......@@ -159,30 +159,20 @@ G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
{
EInside positionA = fPtrSolidA->Inside(p);
if (positionA == kOutside) return kOutside;
if (positionA == kOutside) { return positionA; } // outside A
EInside positionB = fPtrSolidB->Inside(p);
if(positionA == kInside && positionB == kOutside)
{
return kInside ;
}
else
{
if(( positionA == kInside && positionB == kSurface) ||
( positionB == kOutside && positionA == kSurface) ||
( positionA == kSurface && positionB == kSurface &&
( fPtrSolidA->SurfaceNormal(p) -
fPtrSolidB->SurfaceNormal(p) ).mag2() >
1000.0*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
{
return kSurface;
}
else
{
return kOutside;
}
}
if (positionB == kOutside) { return positionA; }
if (positionB == kInside) { return kOutside; }
if (positionA == kInside) { return kSurface; } // surface B
// Point is on both surfaces
//
static const G4double rtol = 1000*kCarTolerance;
return ((fPtrSolidA->SurfaceNormal(p) -
fPtrSolidB->SurfaceNormal(p)).mag2() > rtol) ? kSurface : kOutside;
}
//////////////////////////////////////////////////////////////
......@@ -193,8 +183,11 @@ G4ThreeVector
G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
{
G4ThreeVector normal;
EInside insideThis= Inside(p);
if( insideThis == kOutside )
EInside InsideA = fPtrSolidA->Inside(p);
EInside InsideB = fPtrSolidB->Inside(p);
if( InsideA == kOutside )
{
#ifdef G4BOOLDEBUG
G4cout << "WARNING - Invalid call [1] in "
......@@ -206,46 +199,41 @@ G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
<< " Point p is outside !" << G4endl;
G4cerr << " p = " << p << G4endl;
#endif
normal = fPtrSolidA->SurfaceNormal(p) ;
}
else
{
EInside InsideA = fPtrSolidA->Inside(p);
EInside InsideB = fPtrSolidB->Inside(p);
if( InsideA == kSurface &&
InsideB != kInside )
else if( InsideA == kSurface &&
InsideB != kInside )
{
normal = fPtrSolidA->SurfaceNormal(p) ;
}
else if( InsideA == kInside &&
InsideB != kOutside )
{
normal = -fPtrSolidB->SurfaceNormal(p) ;
}
else
{
if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
{
normal = fPtrSolidA->SurfaceNormal(p) ;
}
else if( InsideA == kInside &&
InsideB != kOutside )
else
{
normal = -fPtrSolidB->SurfaceNormal(p) ;
}
else
{
if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
{
normal = fPtrSolidA->SurfaceNormal(p) ;
}
else
{
normal = -fPtrSolidB->SurfaceNormal(p) ;
}
#ifdef G4BOOLDEBUG
if(insideThis == kInside)
{
G4cout << "WARNING - Invalid call [2] in "
if(Inside(p) == kInside)
{
G4cout << "WARNING - Invalid call [2] in "
<< "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
<< " Point p is inside !" << G4endl;
G4cout << " p = " << p << G4endl;
G4cerr << "WARNING - Invalid call [2] in "
G4cout << " p = " << p << G4endl;
G4cerr << "WARNING - Invalid call [2] in "
<< "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
<< " Point p is inside !" << G4endl;
G4cerr << " p = " << p << G4endl;
}
#endif
G4cerr << " p = " << p << G4endl;
}
#endif
}
return normal;
}
......@@ -258,7 +246,7 @@ G4double
G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
const G4ThreeVector& v ) const
{
G4double dist = 0.0,disTmp = 0.0 ;
G4double dist = 0.0, dist2 = 0.0, disTmp = 0.0;
#ifdef G4BOOLDEBUG
if( Inside(p) == kInside )
......@@ -296,8 +284,10 @@ G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
if( Inside(p+dist*v) == kOutside )
{
disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
dist += disTmp ;
disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
dist2 = dist+disTmp;
if (dist == dist2) { return dist; } // no progress
dist = dist2 ;
count1++;
if( count1 > 1000 ) // Infinite loop detected
{
......@@ -351,7 +341,9 @@ G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
{
return kInfinity ;
}
dist += disTmp ;
dist2 = dist+disTmp;
if (dist == dist2) { return dist; } // no progress
dist = dist2 ;
count2++;
if( count2 > 1000 ) // Infinite loop detected
{
......
source/geometry/solids/Boolean/src/G4UnionSolid.cc
View file @ 292425bd
......@@ -24,9 +24,9 @@
// ********************************************************************
//
//
// $Id: G4UnionSolid.cc 95297 2016-02-04 09:30:14Z gcosmo $
// $Id: G4UnionSolid.cc 109443 2018-04-24 06:28:58Z evc $
//
// Implementation of methods for the class G4IntersectionSolid
// Implementation of methods for the class G4UnionSolid
//
// History:
//
......@@ -34,6 +34,8 @@
// 28.11.98 V.Grichine: fix while loops in DistToIn/Out
// 27.07.99 V.Grichine: modifications in DistToOut(p,v,...), while -> do-while
// 16.03.01 V.Grichine: modifications in CalculateExtent()
// 17.03.17 E.Tcherniaev: revision of SurfaceNormal()
// 23.04.18 E.Tcherniaev: added extended BBox, yearly return in Inside()
//
// --------------------------------------------------------------------
......@@ -60,6 +62,18 @@ G4UnionSolid:: G4UnionSolid( const G4String& pName,
G4VSolid* pSolidB )
: G4BooleanSolid(pName,pSolidA,pSolidB)
{
G4VoxelLimits voxelLimits;
G4AffineTransform affineTransform;
G4double vmin, vmax;
CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setX((vmax + vmin)/2.);
fBoxSize.setX((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setY((vmax + vmin)/2.);
fBoxSize.setY((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setZ((vmax + vmin)/2.);
fBoxSize.setZ((vmax - vmin + kCarTolerance)/2.);
}
/////////////////////////////////////////////////////////////////////
......@@ -74,6 +88,18 @@ G4UnionSolid::G4UnionSolid( const G4String& pName,
: G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
{
G4VoxelLimits voxelLimits;
G4AffineTransform affineTransform;
G4double vmin, vmax;
CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setX((vmax + vmin)/2.);
fBoxSize.setX((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setY((vmax + vmin)/2.);
fBoxSize.setY((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setZ((vmax + vmin)/2.);
fBoxSize.setZ((vmax - vmin + kCarTolerance)/2.);
}
///////////////////////////////////////////////////////////
......@@ -86,6 +112,18 @@ G4UnionSolid::G4UnionSolid( const G4String& pName,
const G4Transform3D& transform )
: G4BooleanSolid(pName,pSolidA,pSolidB,transform)
{
G4VoxelLimits voxelLimits;
G4AffineTransform affineTransform;
G4double vmin, vmax;
CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setX((vmax + vmin)/2.);
fBoxSize.setX((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setY((vmax + vmin)/2.);
fBoxSize.setY((vmax - vmin + kCarTolerance)/2.);
CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
fBoxCenter.setZ((vmax + vmin)/2.);
fBoxSize.setZ((vmax - vmin + kCarTolerance)/2.);
}
//////////////////////////////////////////////////////////////////
......@@ -113,6 +151,8 @@ G4UnionSolid::~G4UnionSolid()
G4UnionSolid::G4UnionSolid(const G4UnionSolid& rhs)
: G4BooleanSolid (rhs)
{
fBoxCenter = rhs.fBoxCenter;
fBoxSize = rhs.fBoxSize;
}
///////////////////////////////////////////////////////////////
......@@ -129,12 +169,14 @@ G4UnionSolid& G4UnionSolid::operator = (const G4UnionSolid& rhs)
//
G4BooleanSolid::operator=(rhs);
fBoxCenter = rhs.fBoxCenter;
fBoxSize = rhs.fBoxSize;
return *this;
}
///////////////////////////////////////////////////////////////
//
//////////////////////////////////////////////////////////////////////////
//
// Calculate extent under transform and specified limit
G4bool
G4UnionSolid::CalculateExtent( const EAxis pAxis,
......@@ -149,15 +191,18 @@ G4UnionSolid::CalculateExtent( const EAxis pAxis,
touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
pTransform, minA, maxA);
touchesB= fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
pTransform, minB, maxB);
touchesB = fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
pTransform, minB, maxB);
if( touchesA || touchesB )
{
pMin = std::min( minA, minB );
pMax = std::max( maxA, maxB );
out = true ;
}
else out = false ;
else
{
out = false ;
}
return out ; // It exists in this slice if either one does.
}
......@@ -170,78 +215,65 @@ G4UnionSolid::CalculateExtent( const EAxis pAxis,
EInside G4UnionSolid::Inside( const G4ThreeVector& p ) const
{
EInside positionA = fPtrSolidA->Inside(p);
if (positionA == kInside) { return kInside; }
// Check if point is outside of bounding box
//
if (std::abs(p.z()-fBoxCenter.z()) > fBoxSize.z()) return kOutside;
EInside positionA = fPtrSolidA->Inside(p);
if (positionA == kInside) { return positionA; } // inside A
EInside positionB = fPtrSolidB->Inside(p);
if (positionA == kOutside) { return positionB; }
if( positionB == kInside ||
( positionA == kSurface && positionB == kSurface &&
( fPtrSolidA->SurfaceNormal(p) +
fPtrSolidB->SurfaceNormal(p) ).mag2() <
1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
{
return kInside;
}
else
{
if( ( positionB == kSurface ) || ( positionA == kSurface ) )
{ return kSurface; }
else
{ return kOutside; }
}
if (positionB == kInside) { return positionB; } // inside B
if (positionB == kOutside) { return positionA; } // surface A
// Point is on both surfaces
//
static const G4double rtol
= 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
return ((fPtrSolidA->SurfaceNormal(p) +
fPtrSolidB->SurfaceNormal(p)).mag2() < rtol) ? kInside : kSurface;
}
//////////////////////////////////////////////////////////////
//
//
// Get surface normal
G4ThreeVector
G4UnionSolid::SurfaceNormal( const G4ThreeVector& p ) const
{
G4ThreeVector normal;
EInside positionA = fPtrSolidA->Inside(p);
EInside positionB = fPtrSolidB->Inside(p);
#ifdef G4BOOLDEBUG
if( Inside(p) == kOutside )
{
G4cout << "WARNING - Invalid call in "
<< "G4UnionSolid::SurfaceNormal(p)" << G4endl
<< " Point p is outside !" << G4endl;
G4cout << " p = " << p << G4endl;
G4cerr << "WARNING - Invalid call in "
<< "G4UnionSolid::SurfaceNormal(p)" << G4endl
<< " Point p is outside !" << G4endl;
G4cerr << " p = " << p << G4endl;
}
#endif
if (positionA == kSurface &&
positionB == kOutside) return fPtrSolidA->SurfaceNormal(p);
if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside)
{
normal= fPtrSolidA->SurfaceNormal(p) ;
}
else if(fPtrSolidB->Inside(p) == kSurface &&
fPtrSolidA->Inside(p) != kInside)
if (positionA == kOutside &&
positionB == kSurface) return fPtrSolidB->SurfaceNormal(p);
if (positionA == kSurface &&
positionB == kSurface)
{
if (Inside(p) == kSurface)
{
normal= fPtrSolidB->SurfaceNormal(p) ;
G4ThreeVector normalA = fPtrSolidA->SurfaceNormal(p);
G4ThreeVector normalB = fPtrSolidB->SurfaceNormal(p);
return (normalA + normalB).unit();
}
else
{
normal= fPtrSolidA->SurfaceNormal(p) ;
}
#ifdef G4BOOLDEBUG
if(Inside(p)==kInside)
{
G4cout << "WARNING - Invalid call in "
<< "G4UnionSolid::SurfaceNormal(p)" << G4endl
<< " Point p is inside !" << G4endl;
G4cout << " p = " << p << G4endl;
G4cerr << "WARNING - Invalid call in "
<< "G4UnionSolid::SurfaceNormal(p)" << G4endl
<< " Point p is inside !" << G4endl;
G4cerr << " p = " << p << G4endl;
}
G4String surf[3] = { "OUTSIDE", "SURFACE", "INSIDE" };
std::ostringstream message;
G4int oldprc = message.precision(16);
message << "Invalid call of SurfaceNormal(p) for union solid: "
<< GetName() << " !"
<< "\nPoint p" << p << " is " << surf[Inside(p)] << " !!!";
message.precision(oldprc);
G4Exception("G4UnionSolid::SurfaceNormal()", "GeomMgt0001",
JustWarning, message);
#endif
}
return normal;
return fPtrSolidA->SurfaceNormal(p);
}
/////////////////////////////////////////////////////////////
......
source/geometry/solids/CSG/include/G4Box.hh
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4Box.hh 79491 2014-03-05 15:24:29Z gcosmo $
// $Id: G4Box.hh 104311 2017-05-24 12:59:46Z gcosmo $
//
// --------------------------------------------------------------------
// GEANT 4 class header file
......@@ -75,10 +75,12 @@ class G4Box : public G4CSGSolid
const G4int n,
const G4VPhysicalVolume* pRep);
void BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const;
G4bool CalculateExtent(const EAxis pAxis,
const G4VoxelLimits& pVoxelLimit,
const G4AffineTransform& pTransform,
G4double& pmin, G4double& pmax) const;
G4double& pMin, G4double& pMax) const;
// Accessors and modifiers
......@@ -128,18 +130,6 @@ class G4Box : public G4CSGSolid
G4Box& operator=(const G4Box& rhs);
// Copy constructor and assignment operator.
protected: // with description
G4ThreeVectorList*
CreateRotatedVertices(const G4AffineTransform& pTransform) const;
// Create the List of transformed vertices in the format required
// for G4VSolid:: ClipCrossSection and ClipBetweenSections.
protected: // without description
enum ESide {kUndefined,kPX,kMX,kPY,kMY,kPZ,kMZ};
// Codes for faces (kPX= +x face, kMY= -y face, etc...)
private:
G4ThreeVector ApproxSurfaceNormal( const G4ThreeVector& p) const;
......
source/geometry/solids/CSG/include/G4Box.icc
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4Box.icc 66356 2012-12-18 09:02:32Z gcosmo $
// $Id: G4Box.icc 66241 2012-12-13 18:34:42Z gunter $
//
// --------------------------------------------------------------------
// GEANT 4 inline definitions file
......
source/geometry/solids/CSG/include/G4Cons.hh
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4Cons.hh 79491 2014-03-05 15:24:29Z gcosmo $
// $Id: G4Cons.hh 104311 2017-05-24 12:59:46Z gcosmo $
//
//
// --------------------------------------------------------------------
......@@ -105,6 +105,10 @@ class G4Cons : public G4CSGSolid
inline G4double GetZHalfLength() const;
inline G4double GetStartPhiAngle() const;
inline G4double GetDeltaPhiAngle() const;
inline G4double GetSinStartPhi() const;
inline G4double GetCosStartPhi() const;
inline G4double GetSinEndPhi() const;
inline G4double GetCosEndPhi() const;
// Modifiers
......@@ -125,10 +129,12 @@ class G4Cons : public G4CSGSolid
const G4int n,
const G4VPhysicalVolume* pRep );
void BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const;
G4bool CalculateExtent( const EAxis pAxis,
const G4VoxelLimits& pVoxelLimit,
const G4AffineTransform& pTransform,
G4double& pmin, G4double& pmax ) const;
G4double& pMin, G4double& pMax ) const;
EInside Inside( const G4ThreeVector& p ) const;
......@@ -181,9 +187,6 @@ class G4Cons : public G4CSGSolid
private:
G4ThreeVectorList*
CreateRotatedVertices(const G4AffineTransform& pTransform) const;
inline void Initialize();
//
// Reset relevant values to zero
......
source/geometry/solids/CSG/include/G4Cons.icc
View file @ 292425bd
......@@ -24,7 +24,7 @@
// ********************************************************************
//
//
// $Id: G4Cons.icc 83572 2014-09-01 15:23:27Z gcosmo $
// $Id: G4Cons.icc 100457 2016-10-23 11:22:46Z evc $
//
// --------------------------------------------------------------------
// GEANT 4 inline definitions file
......@@ -76,6 +76,30 @@ G4double G4Cons::GetDeltaPhiAngle() const
return fDPhi;
}
inline
G4double G4Cons::GetSinStartPhi() const
{
return sinSPhi;
}
inline
G4double G4Cons::GetCosStartPhi() const
{
return cosSPhi;
}
inline
G4double G4Cons::GetSinEndPhi() const
{
return sinEPhi;
}
inline
G4double G4Cons::GetCosEndPhi() const
{
return cosEPhi;
}
inline
void G4Cons::Initialize()
{
......