Creates user defined polyhedron. This function allows to the user to define arbitrary polyhedron. The faces of the polyhedron should be either triangles or planar quadrilateral. Nodes of a face are defined by indexes pointing to the elements in the xyz array. Numeration of the elements in the array starts from 1 (like in fortran). The indexes can be positive or negative. Negative sign means that the corresponding edge is invisible. The normal of the face should be directed to exterior of the polyhedron.
826{
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849 const G4int numSide =
855 double3* xyz;
856 typedef G4int int4[4];
857 int4* faces_vec;
859 {
860
861
862
863 std::vector<G4bool> chopped(
numCorner,
false);
864 std::vector<G4int*> triQuads;
867 while (remaining >= 3)
868 {
869
870
872 G4int iStepper = iStarter;
873 do
874 {
875 if (A < 0) {
A = iStepper; }
876 else if (B < 0) {
B = iStepper; }
877 else if (C < 0) {
C = iStepper; }
878 do
879 {
880 if (++iStepper >=
numCorner) { iStepper = 0; }
881 }
882 while (chopped[iStepper]);
883 }
884 while (C < 0 && iStepper != iStarter);
885
886
887
888
894 {
899 triQuads.push_back(tq);
901 --remaining;
902 }
903 else
904 {
905 do
906 {
907 if (++iStarter >=
numCorner) { iStarter = 0; }
908 }
909 while (chopped[iStarter]);
910 }
911 }
912
913
915 nFaces = numSide *
numCorner + 2 * triQuads.size();
916 faces_vec = new int4[nFaces];
920 for (
G4int iEnd = 0; iEnd < 2; ++iEnd)
921 {
922 for (size_t i = 0; i < triQuads.size(); ++i)
923 {
924
925
927 if (iEnd == 0)
928 {
929 a = triQuads[i][0];
930 b = triQuads[i][1];
931 c = triQuads[i][2];
932 }
933 else
934 {
935 a = triQuads[i][0] + addition;
936 b = triQuads[i][2] + addition;
937 c = triQuads[i][1] + addition;
938 }
939 G4int ab = std::abs(b - a);
940 G4int bc = std::abs(c - b);
941 G4int ca = std::abs(a - c);
942 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
943 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
944 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
945 faces_vec[iface][3] = 0;
946 ++iface;
947 }
948 }
949
950
951
952 xyz = new double3[nNodes];
956 for (
G4int iSide = 0; iSide < numSide; ++iSide)
957 {
959 {
960 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
961 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
963 if (iSide == 0)
964 {
966 {
967 faces_vec[iface][0] = ixyz + 1;
968 faces_vec[iface][1] = -(ixyz +
numCorner + 1);
969 faces_vec[iface][2] = ixyz +
numCorner + 2;
970 faces_vec[iface][3] = ixyz + 2;
971 }
972 else
973 {
974 faces_vec[iface][0] = ixyz + 1;
975 faces_vec[iface][1] = -(ixyz +
numCorner + 1);
976 faces_vec[iface][2] = ixyz + 2;
977 faces_vec[iface][3] = ixyz -
numCorner + 2;
978 }
979 }
980 else if (iSide == numSide - 1)
981 {
983 {
984 faces_vec[iface][0] = ixyz + 1;
985 faces_vec[iface][1] = ixyz +
numCorner + 1;
986 faces_vec[iface][2] = ixyz +
numCorner + 2;
987 faces_vec[iface][3] = -(ixyz + 2);
988 }
989 else
990 {
991 faces_vec[iface][0] = ixyz + 1;
992 faces_vec[iface][1] = ixyz +
numCorner + 1;
993 faces_vec[iface][2] = ixyz + 2;
994 faces_vec[iface][3] = -(ixyz -
numCorner + 2);
995 }
996 }
997 else
998 {
1000 {
1001 faces_vec[iface][0] = ixyz + 1;
1002 faces_vec[iface][1] = -(ixyz +
numCorner + 1);
1003 faces_vec[iface][2] = ixyz +
numCorner + 2;
1004 faces_vec[iface][3] = -(ixyz + 2);
1005 }
1006 else
1007 {
1008 faces_vec[iface][0] = ixyz + 1;
1009 faces_vec[iface][1] = -(ixyz +
numCorner + 1);
1010 faces_vec[iface][2] = ixyz + 2;
1011 faces_vec[iface][3] = -(ixyz -
numCorner + 2);
1012 }
1013 }
1014 ++iface;
1015 ++ixyz;
1016 }
1017 phi += dPhi;
1018 }
1019
1020
1021
1023 {
1024 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1025 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1027 ++ixyz;
1028 }
1029 }
1030 else
1031 {
1034 xyz = new double3[nNodes];
1035 faces_vec = new int4[nFaces];
1038 G4int ixyz = 0, iface = 0;
1039 for (
G4int iSide = 0; iSide < numSide; ++iSide)
1040 {
1042 {
1043 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1044 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1046
1047 if (iSide < numSide - 1)
1048 {
1050 {
1051 faces_vec[iface][0] = ixyz + 1;
1052 faces_vec[iface][1] = -(ixyz +
numCorner + 1);
1053 faces_vec[iface][2] = ixyz +
numCorner + 2;
1054 faces_vec[iface][3] = -(ixyz + 2);
1055 }
1056 else
1057 {
1058 faces_vec[iface][0] = ixyz + 1;
1059 faces_vec[iface][1] = -(ixyz +
numCorner + 1);
1060 faces_vec[iface][2] = ixyz + 2;
1061 faces_vec[iface][3] = -(ixyz -
numCorner + 2);
1062 }
1063 }
1064 else
1065 {
1067 {
1068 faces_vec[iface][0] = ixyz + 1;
1069 faces_vec[iface][1] = -(ixyz +
numCorner - nFaces + 1);
1070 faces_vec[iface][2] = ixyz +
numCorner - nFaces + 2;
1071 faces_vec[iface][3] = -(ixyz + 2);
1072 }
1073 else
1074 {
1075 faces_vec[iface][0] = ixyz + 1;
1076 faces_vec[iface][1] = -(ixyz - nFaces +
numCorner + 1);
1077 faces_vec[iface][2] = ixyz - nFaces + 2;
1078 faces_vec[iface][3] = -(ixyz -
numCorner + 2);
1079 }
1080 }
1081 ++ixyz;
1082 ++iface;
1083 }
1084 phi += dPhi;
1085 }
1086 }
1089 delete [] faces_vec;
1090 delete [] xyz;
1091 if (prob)
1092 {
1093 std::ostringstream message;
1094 message <<
"Problem creating G4Polyhedron for: " <<
GetName();
1095 G4Exception(
"G4GenericPolycone::CreatePolyhedron()",
"GeomSolids1002",
1097 delete polyhedron;
1098 return nullptr;
1099 }
1100 else
1101 {
1102 return polyhedron;
1103 }
1104}
double B(double temperature)
double A(double temperature)
static G4int GetNumberOfRotationSteps()
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])