2933{
2934
2936 if (nnodes == 0)
2937 {
2938 std::cerr
2939 << "HepPolyhedronTetMesh: Empty tetrahedron mesh" << std::endl;
2940 return;
2941 }
2942 G4int ntet = nnodes/4;
2943 if (nnodes != ntet*4)
2944 {
2945 std::cerr << "HepPolyhedronTetMesh: Number of nodes = " << nnodes
2946 << " in tetrahedron mesh is NOT multiple of 4"
2947 << std::endl;
2948 return;
2949 }
2950
2951
2952
2953
2954 std::vector<G4int> iheads(nnodes, -1);
2955 std::vector<std::pair<G4int,G4int>> ipairs(nnodes,std::pair(-1,-1));
2956 for (
G4int i = 0; i < nnodes; ++i)
2957 {
2958
2960 auto key = std::hash<G4double>()(point.
x());
2961 key ^= std::hash<G4double>()(point.
y());
2962 key ^= std::hash<G4double>()(point.
z());
2963 key %= nnodes;
2964
2965 if (iheads[key] < 0)
2966 {
2967 iheads[key] = i;
2968 ipairs[i].first = i;
2969 continue;
2970 }
2971
2972 for (
G4int icur = iheads[key], iprev = 0;;)
2973 {
2974 G4int icheck = ipairs[icur].first;
2975 if (tetrahedra[icheck] == point)
2976 {
2977 ipairs[i].first = icheck;
2978 break;
2979 }
2980 iprev = icur;
2981 icur = ipairs[icur].second;
2982
2983 if (icur < 0)
2984 {
2985 ipairs[i].first = i;
2986 ipairs[iprev].second = i;
2987 break;
2988 }
2989 }
2990 }
2991
2992
2993 struct facet
2994 {
2996 facet() : i1(0), i2(0), i3(0) {};
2998 };
2999 G4int nfacets = nnodes;
3000 std::vector<facet> ifacets(nfacets);
3001 for (
G4int i = 0; i < nfacets; i += 4)
3002 {
3003 G4int i0 = ipairs[i + 0].first;
3004 G4int i1 = ipairs[i + 1].first;
3005 G4int i2 = ipairs[i + 2].first;
3006 G4int i3 = ipairs[i + 3].first;
3007 if (i0 > i1) std::swap(i0, i1);
3008 if (i0 > i2) std::swap(i0, i2);
3009 if (i0 > i3) std::swap(i0, i3);
3010 if (i1 > i2) std::swap(i1, i2);
3011 if (i1 > i3) std::swap(i1, i3);
3016 if (volume > 0.) std::swap(i2, i3);
3017 ifacets[i + 0] = facet(i0, i1, i2);
3018 ifacets[i + 1] = facet(i0, i2, i3);
3019 ifacets[i + 2] = facet(i0, i3, i1);
3020 ifacets[i + 3] = facet(i1, i3, i2);
3021 }
3022
3023
3024 std::fill(iheads.begin(), iheads.end(), -1);
3025 std::fill(ipairs.begin(), ipairs.end(), std::pair(-1,-1));
3026 for (
G4int i = 0; i < nfacets; ++i)
3027 {
3028
3029 G4int key = ifacets[i].i1;
3030 if (iheads[key] < 0)
3031 {
3032 iheads[key] = i;
3033 ipairs[i].first = i;
3034 continue;
3035 }
3036
3037 G4int i2 = ifacets[i].i2, i3 = ifacets[i].i3;
3038 for (
G4int icur = iheads[key], iprev = -1;;)
3039 {
3040 G4int icheck = ipairs[icur].first;
3041 if (ifacets[icheck].i2 == i3 && ifacets[icheck].i3 == i2)
3042 {
3043 if (iprev < 0)
3044 {
3045 iheads[key] = ipairs[icur].second;
3046 }
3047 else
3048 {
3049 ipairs[iprev].second = ipairs[icur].second;
3050 }
3051 ipairs[icur].first = -1;
3052 ipairs[icur].second = -1;
3053 break;
3054 }
3055 iprev = icur;
3056 icur = ipairs[icur].second;
3057
3058 if (icur < 0)
3059 {
3060 ipairs[i].first = i;
3061 ipairs[iprev].second = i;
3062 break;
3063 }
3064 }
3065 }
3066
3067
3068 std::fill(iheads.begin(), iheads.end(), -1);
3069 G4int nver = 0, nfac = 0;
3070 for (
G4int i = 0; i < nfacets; ++i)
3071 {
3072 if (ipairs[i].first < 0) continue;
3073 G4int i1 = ifacets[i].i1;
3074 G4int i2 = ifacets[i].i2;
3075 G4int i3 = ifacets[i].i3;
3076 if (iheads[i1] < 0) iheads[i1] = nver++;
3077 if (iheads[i2] < 0) iheads[i2] = nver++;
3078 if (iheads[i3] < 0) iheads[i3] = nver++;
3079 nfac++;
3080 }
3081
3082
3084 for (
G4int i = 0; i < nnodes; ++i)
3085 {
3086 G4int k = iheads[i];
3087 if (k >= 0)
SetVertex(k + 1, tetrahedra[i]);
3088 }
3089 for (
G4int i = 0, k = 0; i < nfacets; ++i)
3090 {
3091 if (ipairs[i].first < 0) continue;
3092 G4int i1 = iheads[ifacets[i].i1] + 1;
3093 G4int i2 = iheads[ifacets[i].i2] + 1;
3094 G4int i3 = iheads[ifacets[i].i3] + 1;
3096 }
3098}
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
void SetVertex(G4int index, const G4Point3D &v)
void SetFacet(G4int index, G4int iv1, G4int iv2, G4int iv3, G4int iv4=0)
void AllocateMemory(G4int Nvert, G4int Nface)