G4QuadrangularFacet.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd,   *
00020 // * subject to DEFCON 705 IPR conditions.                            *
00021 // * By using,  copying,  modifying or  distributing the software (or *
00022 // * any work based  on the software)  you  agree  to acknowledge its *
00023 // * use  in  resulting  scientific  publications,  and indicate your *
00024 // * acceptance of all terms of the Geant4 Software license.          *
00025 // ********************************************************************
00026 //
00027 //
00028 // $Id: G4QuadrangularFacet.cc 67011 2013-01-29 16:17:41Z gcosmo $
00029 //
00030 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00031 //
00032 // CHANGE HISTORY
00033 // --------------
00034 //
00035 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
00036 // 12 October 2012, M Gayer, CERN
00037 //                  New implementation reducing memory requirements by 50%,
00038 //                  and considerable CPU speedup together with the new
00039 //                  implementation of G4TessellatedSolid.
00040 //
00041 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00042 
00043 #include "G4QuadrangularFacet.hh"
00044 #include "geomdefs.hh"
00045 #include "Randomize.hh"
00046  
00047 using namespace std;
00048 
00050 //
00051 // !!!THIS IS A FUDGE!!!  IT'S TWO ADJACENT G4TRIANGULARFACETS
00052 // --- NOT EFFICIENT BUT PRACTICAL.
00053 //
00054 G4QuadrangularFacet::G4QuadrangularFacet (const G4ThreeVector &vt0,
00055                                           const G4ThreeVector &vt1,
00056                                           const G4ThreeVector &vt2,
00057                                           const G4ThreeVector &vt3,
00058                                                 G4FacetVertexType vertexType)
00059 {
00060   G4ThreeVector e1, e2, e3;
00061 
00062   SetVertex(0, vt0);
00063   if (vertexType == ABSOLUTE)
00064   {
00065     SetVertex(1, vt1);
00066     SetVertex(2, vt2);
00067     SetVertex(3, vt3);
00068 
00069     e1 = vt1 - vt0;
00070     e2 = vt2 - vt0;
00071     e3 = vt3 - vt0;
00072   }
00073   else
00074   {
00075     SetVertex(1, vt0 + vt1);
00076     SetVertex(2, vt0 + vt2);
00077     SetVertex(3, vt0 + vt3);
00078 
00079     e1 = vt1;
00080     e2 = vt2;
00081     e3 = vt3;
00082   }
00083   G4double length1 = e1.mag();
00084   G4double length2 = (GetVertex(2)-GetVertex(1)).mag();
00085   G4double length3 = (GetVertex(3)-GetVertex(2)).mag();
00086   G4double length4 = e3.mag();
00087 
00088   G4ThreeVector normal1 = e1.cross(e2).unit();
00089   G4ThreeVector normal2 = e2.cross(e3).unit(); 
00090 
00091   bool isDefined = (length1 > kCarTolerance && length2 > kCarTolerance &&
00092     length3 > kCarTolerance && length4 > kCarTolerance &&
00093     normal1.dot(normal2) >= 0.9999999999);
00094 
00095   if (isDefined)
00096   {
00097     fFacet1 = G4TriangularFacet (GetVertex(0),GetVertex(1),
00098                                  GetVertex(2),ABSOLUTE);
00099     fFacet2 = G4TriangularFacet (GetVertex(0),GetVertex(2),
00100                                  GetVertex(3),ABSOLUTE);
00101 
00102     G4TriangularFacet facet3 (GetVertex(0),GetVertex(1),GetVertex(3),ABSOLUTE);
00103     G4TriangularFacet facet4 (GetVertex(1),GetVertex(2),GetVertex(3),ABSOLUTE);
00104 
00105     G4ThreeVector normal12 = fFacet1.GetSurfaceNormal()
00106                            + fFacet2.GetSurfaceNormal();
00107     G4ThreeVector normal34 = facet3.GetSurfaceNormal()
00108                            + facet4.GetSurfaceNormal();
00109     G4ThreeVector normal = 0.25 * (normal12 + normal34);
00110 
00111     fFacet1.SetSurfaceNormal (normal);
00112     fFacet2.SetSurfaceNormal (normal);
00113 
00114     G4ThreeVector vtmp = 0.5 * (e1 + e2);
00115     fCircumcentre = GetVertex(0) + vtmp;
00116     G4double radiusSqr = vtmp.mag2();
00117     fRadius = std::sqrt(radiusSqr);
00118   }
00119   else
00120   {
00121     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
00122                 "GeomSolids0002", JustWarning,
00123                 "Length of sides of facet are too small or sides not planar.");
00124     G4cout << G4endl;
00125     G4cout << "P0 = " << GetVertex(0) << G4endl;
00126     G4cout << "P1 = " << GetVertex(1) << G4endl;
00127     G4cout << "P2 = " << GetVertex(2) << G4endl;
00128     G4cout << "P3 = " << GetVertex(3) << G4endl;
00129     G4cout << "Side lengths = P0->P1" << length1 << G4endl;    
00130     G4cout << "Side lengths = P1->P2" << length2 << G4endl;    
00131     G4cout << "Side lengths = P2->P3" << length3 << G4endl;    
00132     G4cout << "Side lengths = P3->P0" << length4 << G4endl;    
00133     G4cout << G4endl;
00134     fRadius = 0.0;
00135   }
00136 }
00137 
00139 //
00140 G4QuadrangularFacet::~G4QuadrangularFacet ()
00141 {
00142 }
00143 
00145 //
00146 G4QuadrangularFacet::G4QuadrangularFacet (const G4QuadrangularFacet &rhs)
00147   : G4VFacet(rhs)
00148 {
00149   fFacet1 = rhs.fFacet1;
00150   fFacet2 = rhs.fFacet2;
00151   fRadius = 0.0;
00152 }
00153 
00155 //
00156 G4QuadrangularFacet &
00157 G4QuadrangularFacet::operator=(const G4QuadrangularFacet &rhs)
00158 {
00159   if (this == &rhs)
00160     return *this;
00161 
00162   fFacet1 = rhs.fFacet1;
00163   fFacet2 = rhs.fFacet2;
00164   fRadius = 0.0;
00165 
00166   return *this;
00167 }
00168 
00170 //
00171 G4VFacet *G4QuadrangularFacet::GetClone ()
00172 {
00173   G4QuadrangularFacet *c = new G4QuadrangularFacet (GetVertex(0), GetVertex(1),
00174                                                     GetVertex(2), GetVertex(3),
00175                                                     ABSOLUTE);
00176   return c;
00177 }
00178 
00180 //
00181 G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector &p)
00182 {
00183   G4ThreeVector v1 = fFacet1.Distance(p);
00184   G4ThreeVector v2 = fFacet2.Distance(p);
00185 
00186   if (v1.mag2() < v2.mag2()) return v1;
00187   else return v2;
00188 }
00189 
00191 //
00192 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
00193                                               G4double)
00194 {  
00195   G4double dist = Distance(p).mag();
00196   return dist;
00197 }
00198 
00200 //
00201 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p, G4double,
00202                                         const G4bool outgoing)
00203 {
00204   G4double dist;
00205 
00206   G4ThreeVector v = Distance(p);
00207   G4double dir = v.dot(GetSurfaceNormal());
00208   if ( ((dir > dirTolerance) && (!outgoing))
00209     || ((dir < -dirTolerance) && outgoing))
00210     dist = kInfinity;
00211   else 
00212     dist = v.mag();
00213   return dist;
00214 }
00215 
00217 //
00218 G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
00219 {
00220   G4double ss  = 0;
00221 
00222   for (G4int i = 0; i <= 3; ++i)
00223   {
00224     G4double sp = GetVertex(i).dot(axis);
00225     if (sp > ss) ss = sp;
00226   }
00227   return ss;
00228 }
00229 
00231 //
00232 G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector &p,
00233                                        const G4ThreeVector &v,
00234                                              G4bool outgoing,
00235                                              G4double &distance,
00236                                              G4double &distFromSurface,
00237                                              G4ThreeVector &normal)
00238 {
00239   G4bool intersect =
00240     fFacet1.Intersect(p,v,outgoing,distance,distFromSurface,normal);
00241   if (!intersect) intersect =
00242     fFacet2.Intersect(p,v,outgoing,distance,distFromSurface,normal);
00243   if (!intersect)
00244   {
00245     distance = distFromSurface = kInfinity;
00246     normal.set(0,0,0);
00247   }
00248   return intersect;
00249 }
00250 
00252 //
00253 // Auxiliary method to get a random point on surface
00254 //
00255 G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
00256 {
00257   G4ThreeVector pr = (G4RandFlat::shoot(0.,1.) < 0.5)
00258                    ? fFacet1.GetPointOnFace() : fFacet2.GetPointOnFace();
00259   return pr;
00260 }
00261 
00263 //
00264 // Auxiliary method for returning the surface area
00265 //
00266 G4double G4QuadrangularFacet::GetArea()
00267 {
00268   G4double area = fFacet1.GetArea() + fFacet2.GetArea();
00269   return area;
00270 }
00271 
00273 //
00274 G4String G4QuadrangularFacet::GetEntityType () const
00275 {
00276   return "G4QuadrangularFacet";
00277 }
00278 
00280 //
00281 G4ThreeVector G4QuadrangularFacet::GetSurfaceNormal () const
00282 {
00283   return fFacet1.GetSurfaceNormal();
00284 }

Generated on Mon May 27 17:49:40 2013 for Geant4 by  doxygen 1.4.7