83 std::ostringstream message;
84 message <<
"TwistedTrapBoxSide is not used as a the side of a box: "
87 G4Exception(
"G4TwistBoxSide::G4TwistBoxSide()",
"GeomSolids0002",
120 fRot.rotateZ( AngleSide ) ;
135 :
G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
136 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
137 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.),
138 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.),
212 G4bool IsParallel = false ;
213 G4bool IsConverged = false ;
255 G4bool tmpisvalid = false ;
257 std::vector<Intersection> xbuf ;
271 if ( (std::fabs(p.z()) <= L) &&
fPhiTwist ) {
276 + (
fTAlph*v.x() + v.y())*std::sin(phi)));
282 * (v.y()*std::cos(phi) - v.x()*std::sin(phi))) / q;
291 xbuf.push_back(xbuftmp) ;
302 areacode[0], isvalid[0],
303 0, validate, &gp, &gv);
330 G4cout <<
"coef = " << c[0] <<
" "
344 for (
G4int i = 0 ; i<num ; i++ ) {
347 G4cout <<
"Solution " << i <<
" : " << srd[i] <<
G4endl ;
349 phi = std::fmod(srd[i] , pihalf) ;
351 u = (2*phiyz + 4*
fDz*phi*v.y() - 2*
fdeltaY*phi*v.z()
363 xbuf.push_back(xbuftmp) ;
366 G4cout <<
"solution " << i <<
" = " << phi <<
" , " << u <<
G4endl ;
385 for (
size_t k = 0 ; k<xbuf.size() ; k++ ) {
388 G4cout <<
"Solution " << k <<
" : "
389 <<
"reconstructed phiR = " << xbuf[k].phi
390 <<
", uR = " << xbuf[k].u <<
G4endl ;
396 IsConverged = false ;
398 for (
G4int i = 1 ; i<maxint ; i++ ) {
401 surfacenormal =
NormAng(phi,u) ;
403 deltaX = ( tmpxx - xxonsurface ).mag() ;
404 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
405 if ( theta < 0.001 ) {
414 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
421 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
424 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
430 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
433 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
445 if (tmpdist >= 0) tmpisvalid =
true;
450 if (tmpdist >= 0) tmpisvalid =
true;
455 "Feature NOT implemented !");
467 xbuf[k].distance = tmpdist ;
468 xbuf[k].areacode = tmpareacode ;
469 xbuf[k].isvalid = tmpisvalid ;
478 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
484 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
489 G4int nxxtmp = xbuf.size() ;
491 if ( nxxtmp<2 || IsParallel ) {
509 xbuf.push_back(xbuftmp) ;
525 xbuf.push_back(xbuftmp) ;
527 for (
size_t k = nxxtmp ; k<xbuf.size() ; k++ ) {
530 G4cout <<
"Solution " << k <<
" : "
531 <<
"reconstructed phiR = " << xbuf[k].phi
532 <<
", uR = " << xbuf[k].u <<
G4endl ;
538 IsConverged = false ;
540 for (
G4int i = 1 ; i<maxint ; i++ ) {
543 surfacenormal =
NormAng(phi,u) ;
545 deltaX = ( tmpxx - xxonsurface ).mag() ;
546 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
547 if ( theta < 0.001 ) {
555 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
562 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
565 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
571 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
574 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
586 if (tmpdist >= 0) tmpisvalid =
true;
591 if (tmpdist >= 0) tmpisvalid =
true;
594 G4Exception(
"G4TwistedBoxSide::DistanceToSurface()",
596 "Feature NOT implemented !");
608 xbuf[k].distance = tmpdist ;
609 xbuf[k].areacode = tmpareacode ;
610 xbuf[k].isvalid = tmpisvalid ;
623 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
626 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
632 for (
size_t i = 0 ; i<xbuf.size() ; i++ ) {
634 distance[i] = xbuf[i].distance;
636 areacode[i] = xbuf[i].areacode ;
637 isvalid[i] = xbuf[i].isvalid ;
640 isvalid[i], nxx, validate, &gp, &gv);
643 G4cout <<
"element Nr. " << i
644 <<
", local Intersection = " << xbuf[i].xx
645 <<
", distance = " << xbuf[i].distance
646 <<
", u = " << xbuf[i].u
647 <<
", phi = " << xbuf[i].phi
648 <<
", isvalid = " << xbuf[i].isvalid
657 G4cout << nxx <<
" possible physical solutions found" <<
G4endl ;
658 for (
G4int k= 0 ; k< nxx ; k++ ) {
659 G4cout <<
"global intersection Point found: " << gxx[k] <<
G4endl ;
715 for (
G4int i = 1 ; i<maxint ; i++ ) {
718 surfacenormal =
NormAng(phiR,uR) ;
720 deltaX = ( xx - xxonsurface ).mag() ;
723 G4cout <<
"i = " << i <<
", distance = " << distance[0] <<
", " << deltaX <<
G4endl ;
730 if ( deltaX <= ctol ) { break ; }
739 if ( phiR > halfphi ) phiR = halfphi ;
740 if ( phiR < -halfphi ) phiR = -halfphi ;
741 if ( uR > uMax ) uR = uMax ;
742 if ( uR < -uMax ) uR = -uMax ;
745 distance[0] = ( p - xx ).mag() ;
746 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
751 G4cout <<
"refined solution " << phiR <<
" , " << uR <<
" , " <<
G4endl ;
760 G4cout <<
"intersection Point found: " << gxx[0] <<
G4endl ;
792 G4cout <<
"GetAreaCode: yprime = " << yprime <<
G4endl ;
793 G4cout <<
"Intervall is " << fYAxisMin <<
" to " << fYAxisMax <<
G4endl ;
808 if (yprime < fYAxisMin + ctol) {
810 if (yprime <= fYAxisMin - ctol) isoutside =
true;
812 }
else if (yprime > fYAxisMax - ctol) {
814 if (yprime >= fYAxisMax + ctol) isoutside =
true;
819 if (xx.z() <
fAxisMin[zaxis] + ctol) {
824 if (xx.z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
826 }
else if (xx.z() >
fAxisMax[zaxis] - ctol) {
831 if (xx.z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
839 areacode = tmpareacode;
848 if (yprime < fYAxisMin ) {
850 }
else if (yprime > fYAxisMax) {
861 }
else if (xx.z() >
fAxisMax[zaxis]) {
875 "Feature NOT implemented !");
928 "Method NOT implemented !");
946 direction = direction.unit();
952 direction = direction.unit();
958 direction = direction.unit();
964 direction = direction.unit();
972 "Feature NOT implemented !");
1033 for ( i = 0 ; i<
n ; i++ ) {
1035 z = -
fDz+i*(2.*
fDz)/(n-1) ;
1039 for ( j = 0 ; j<k ; j++ ) {
1041 nnode =
GetNode(i,j,k,n,iside) ;
1042 u = -b/2 +j*b/(k-1) ;
1045 xyz[nnode][0] = p.x() ;
1046 xyz[nnode][1] = p.y() ;
1047 xyz[nnode][2] = p.z() ;
1049 if ( i<n-1 && j<k-1 ) {
1051 nface =
GetFace(i,j,k,n,iside) ;
static const G4int sAxisZ
G4int GetAreacode(G4int i) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4SurfCurNormal fCurrentNormal
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sOutside
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
G4TwistBoxSide(const G4String &name, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph, G4double AngleSide)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
static double normal(HepRandomEngine *eptr)
static const G4int sC0Min1Min
G4GLOB_DLL std::ostream G4cout
G4bool IsOutside(G4int areacode) const
G4double GetDistance(G4int i) const
static const G4int sAxis1
static const G4int sC0Max1Max
static const G4int sBoundary
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
G4bool IsValid(G4int i) const
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
static const G4int sAxis0
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisMin
virtual G4double GetBoundaryMax(G4double phi)
static const G4int sAxisY
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
static const G4int sInside
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
static const G4int sCorner
static const G4double factor
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sAxisMax
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
const G4double x[NPOINTSGL]
G4ThreeVector GetXX(G4int i) const
virtual G4String GetName() const
virtual void SetBoundaries()
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
virtual ~G4TwistBoxSide()
G4ThreeVector ProjectPoint(const G4ThreeVector &p, G4bool isglobal=false)
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
G4double GetValueB(G4double phi)
G4bool DistanceSort(const Intersection &a, const Intersection &b)
static const G4int sC0Max1Min
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
CurrentStatus fCurStatWithV
G4ThreeVector NormAng(G4double phi, G4double u)
virtual void SetCorners()
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)