eic-opticks/u4/U4Mesh.h

0001 #pragma once
0002 /**
0003 U4Mesh.h : Polygonize and Serialize Solids into Triangles/Quads
0004 ==================================================================
0005
0006 This selects and simplifies some parts of the former x4/X4Mesh.{hh,cc}
0007 and it better integrated with the pyvista PolyData face format
0008 allowing 3D visualization with quads as well as triangles.
0009
0010 +------------+------------------------------------------------------------------------+------------------------------------+
0011 | field      |  content                                                               |  thoughts                          |
0012 +============+========================================================================+====================================+
0013 | (NP*)vtx   |  (num_vtx, 3) array of float coordinates                               | basis                              |
0014 +------------+------------------------------------------------------------------------+------------------------------------+
0015 | (NP*)fpd   |  flat vtx index array : funny pyvista irregular 3/4 vertex face format | pv useful, includes quads, DROP?   |
0016 +------------+------------------------------------------------------------------------+------------------------------------+
0017 | (NP*)face  |  (num_face,4) tri/quad vtx index array using f.face[:,3] == -1 for tri | DROP ?                             |
0018 +------------+------------------------------------------------------------------------+------------------------------------+
0019 | (NP*)tri   |  (num_tri,3) tri vtx index array                                       | standard                           |
0020 +------------+------------------------------------------------------------------------+------------------------------------+
0021 | (NP*)tpd   |  flat vtx index array : funny pyvista using all triangles              | same as tri, pv useful             |
0022 +------------+------------------------------------------------------------------------+------------------------------------+
0023
0024 For use with OpenGL rendering its natural to use "vtx" and "tri".
0025
0026 **/
0027
0028 #include <map>
0029 #include "G4Polyhedron.hh"
0030 #include "ssys.h"
0031 #include "NPX.h"
0032 #include "NPFold.h"
0033
0034
0035 struct U4Mesh
0036 {
0037     static constexpr const char* U4Mesh__NumberOfRotationSteps = "U4Mesh__NumberOfRotationSteps" ;
0038     static constexpr const char* _NumberOfRotationSteps_DUMP = "U4Mesh__NumberOfRotationSteps_DUMP" ;
0039
0040
0041     const G4VSolid* solid ;
0042     const char* entityType ;
0043     const char* solidName ;
0044     int numberOfRotationSteps ;
0045
0046     G4Polyhedron*   poly ;
0047     std::map<int,int> v2fc ;
0048     int             errvtx ;
0049     bool            do_init_vtx_disqualify ;
0050
0051     int             nv_poly, nv, nf ;
0052     NP*             vtx ;
0053     double*         vtx_ ;
0054     NP*             fpd ;  // funny pyvista irregular 3/4 vertex face format
0055
0056     NP*             face ;  // can include tri and quad
0057     int*            face_ ;
0058     int             nt ;
0059     NP*             tri ;  // only triangle indices (more standard gltf suitable layout than other face formats)
0060     int*            tri_ ;
0061     NP*             tpd ;  // funny pyvista face format, but with all 3 vertex
0062
0063     static void Save(const G4VSolid* solid, const char* base ); // base is often "$FOLD"
0064     static void Save(const G4VSolid* solid, const char* base, const char* name);
0065     static NPFold* MakeFold(
0066        const std::vector<const G4VSolid*>& solids,
0067        const std::vector<std::string>& keys
0068       );
0069     static NPFold* Serialize(const G4VSolid* solid) ;
0070     static const char* EType(const G4VSolid* solid);
0071     static const char* SolidName(const G4VSolid* solid);
0072
0073     static std::string FormEKey( const char* prefix, const char* key, const char* val, int idx=-1 );
0074     static int NumberOfRotationSteps(const char* entityType, const char* solidname );
0075     static int NumberOfRotationSteps_solidName_STARTING( const char* _solidName );
0076
0077
0078     static G4Polyhedron* CreatePolyhedron(const G4VSolid* solid, int num);
0079
0080     U4Mesh(const G4VSolid* solid);
0081     void init();
0082
0083     void init_vtx_face_count();
0084     int getVertexFaceCount(int iv) const ;
0085     std::string desc_vtx_face_count() const ;
0086
0087     void init_vtx();
0088     void init_fpd();
0089     NPFold* serialize() const ;
0090     void save(const char* base) const ;
0091
0092     std::string id() const  ;
0093     std::string desc() const  ;
0094
0095     void init_face();
0096     void init_tri();
0097     void init_tpd();
0098
0099 };
0100
0101 inline void U4Mesh::Save(const G4VSolid* solid, const char* base) // static
0102 {
0103     NPFold* fold = U4Mesh::Serialize(solid) ;
0104     fold->save(base);
0105 }
0106
0107 inline void U4Mesh::Save(const G4VSolid* solid, const char* base, const char* name) // static
0108 {
0109     NPFold* fold = U4Mesh::Serialize(solid) ;
0110     fold->set_meta<std::string>("name", name);
0111     fold->save(base, name);
0112 }
0113
0114
0115 /**
0116 U4Mesh::MakeFold
0117 ----------------
0118
0119 **/
0120
0121 inline NPFold* U4Mesh::MakeFold(
0122     const std::vector<const G4VSolid*>& solids,
0123     const std::vector<std::string>& keys
0124    ) // static
0125 {
0126     NPFold* mesh = new NPFold ;
0127     int num_solid = solids.size();
0128     int num_key = keys.size();
0129     assert( num_solid == num_key );
0130
0131     for(int i=0 ; i < num_solid ; i++)
0132     {
0133         int lvid = i ;
0134         const G4VSolid* so = solids[i];
0135         const char* _key = keys[i].c_str();
0136
0137         NPFold* sub = Serialize(so) ;
0138         sub->set_meta<int>("lvid", lvid );
0139
0140         mesh->add_subfold( _key, sub );
0141     }
0142     return mesh ;
0143 }
0144
0145 inline NPFold* U4Mesh::Serialize(const G4VSolid* solid) // static
0146 {
0147     U4Mesh mesh(solid);
0148     return mesh.serialize() ;
0149 }
0150
0151 inline const char* U4Mesh::EType(const G4VSolid* solid)  // static
0152 {
0153     G4GeometryType _etype = solid->GetEntityType();  // G4GeometryType typedef for G4String
0154     return strdup(_etype.c_str()) ;
0155 }
0156
0157
0158 inline const char* U4Mesh::SolidName(const G4VSolid* solid) // static
0159 {
0160     G4String name = solid->GetName();  // forced to duplicate as this returns by value
0161     return strdup(name.c_str());
0162 }
0163
0164 inline std::string U4Mesh::FormEKey( const char* prefix, const char* key, const char* val, int idx  ) // static
0165 {
0166     std::stringstream ss ;
0167     ss << prefix
0168        << "_"
0169        << ( key ? key : "-" )
0170        << "_"
0171        << ( val ? val : "-" )
0172        ;
0173
0174     if( idx > -1 ) ss << "_" << idx ;
0175     std::string str = ss.str();
0176     return str ;
0177 }
0178
0179 /**
0180 U4Mesh::NumberOfRotationSteps
0181 --------------------------------
0182
0183 This static method is invoked by the U4Mesh ctor for every solid
0184 in the geometry with (entityType, solidName) arguments such as
0185 ("G4Torus", "s_EMFcoil_holder_ring26_seg8"). Based on the existence
0186 of envvar keys derived from these inputs values for the number of
0187 rotation steps are obtained or zero if no such keys exist.
0188
0189
0190
0191 Returns envvar configured value if entityType or solidName envvars
0192 are defined, otherwise returns zero.
0193 If both entityType and solidName envvars exist
0194 the solidName value is used as that is more specific.
0195
0196 Example envvar keys::
0197
0198    export U4Mesh__NumberOfRotationSteps_entityType_G4Torus=48
0199    export U4Mesh__NumberOfRotationSteps_solidName_myTorus=48
0200    export U4Mesh__NumberOfRotationSteps_solidName_sTarget=96
0201
0202
0203 Check which solids RotationSteps were customized from the persisted mesh fold::
0204
0205     cd ~/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/mesh
0206     grep Rotation ./NPFold_meta.txt
0207     sSurftube_0V1_0/NPFold_meta.txt:numberOfRotationSteps:480
0208     sSurftube_0V1_1/NPFold_meta.txt:numberOfRotationSteps:480
0209     sSurftube_10V1_0/NPFold_meta.txt:numberOfRotationSteps:480
0210     sSurftube_10V1_1/NPFold_meta.txt:numberOfRotationSteps:480
0211     sSurftube_11V1_0/NPFold_meta.txt:numberOfRotationSteps:480
0212     ...
0213     sTarget/NPFold_meta.txt:numberOfRotationSteps:240
0214     svacSurftube_0V1_0/NPFold_meta.txt:numberOfRotationSteps:480
0215     svacSurftube_0V1_1/NPFold_meta.txt:numberOfRotationSteps:480
0216     svacSurftube_10V1_0/NPFold_meta.txt:numberOfRotationSteps:480
0217     svacSurftube_10V1_1/NPFold_meta.txt:numberOfRotationSteps:480
0218
0219
0220 **/
0221
0222
0223
0224
0225 inline int U4Mesh::NumberOfRotationSteps(const char* _entityType, const char* _solidName )
0226 {
0227     std::string entityType_ekey = FormEKey(U4Mesh__NumberOfRotationSteps, "entityType" , _entityType );
0228     std::string solidName_ekey = FormEKey(U4Mesh__NumberOfRotationSteps, "solidName" , _solidName );
0229     int num_entityType = ssys::getenvint( entityType_ekey.c_str(), 0 );
0230     int num_solidName  = ssys::getenvint( solidName_ekey.c_str(), 0 );
0231     int num_solidName_STARTING = NumberOfRotationSteps_solidName_STARTING( _solidName );
0232
0233     enum { UNRESOLVED, SOLIDNAME, SOLIDNAME_STARTING, ENTITY_TYPE } ;
0234     int resolve = UNRESOLVED ;
0235
0236     int num = 0 ;    // resolution order from most to least specific
0237     if( num_solidName > 0 )
0238     {
0239         num = num_solidName ;
0240         resolve = SOLIDNAME ;
0241     }
0242     else if( num_solidName_STARTING > 0 )
0243     {
0244         num = num_solidName_STARTING ;
0245         resolve = SOLIDNAME_STARTING ;
0246     }
0247     else if( num_entityType > 0 )
0248     {
0249         num = num_entityType ;
0250         resolve = ENTITY_TYPE ;
0251     }
0252
0253     bool DUMP = ssys::getenvbool(_NumberOfRotationSteps_DUMP);
0254
0255     if(DUMP && num > 0) std::cout
0256          << "U4Mesh::NumberOfRotationSteps\n"
0257          << " entityType_ekey[" << entityType_ekey << "]"
0258          << " solidName_ekey[" << solidName_ekey << "]"
0259          << " num_entityType " << num_entityType
0260          << " num_solidName " << num_solidName
0261          << " num_solidName_STARTING " << num_solidName_STARTING
0262          << " num " << num
0263          << " resolve " << resolve
0264          << "\n"
0265          ;
0266
0267     return num  ;
0268 }
0269
0270 /**
0271 U4Mesh::NumberOfRotationSteps_solidName_STARTING
0272 --------------------------------------------------
0273
0274 Returns non-zero value when the solidName matches the value of
0275 several pfx envvars which are paired with corresponding val
0276 envvars that provide the integer value. The first matching
0277 envvar pair provides the value.
0278
0279 Example of the envvar pairs with idx 0 and 1::
0280
0281     export U4Mesh__NumberOfRotationSteps_solidName_STARTING_pfx_0=s_EMFcoil_holder_ring
0282     export U4Mesh__NumberOfRotationSteps_solidName_STARTING_val_0=480
0283
0284     export U4Mesh__NumberOfRotationSteps_solidName_STARTING_pfx_1=s_EMFcoil_holder_ring
0285     export U4Mesh__NumberOfRotationSteps_solidName_STARTING_val_1=16
0286
0287 The value from the first match is used
0288
0289 MAYBE: use sregex.h to generalize this if the need arises
0290 **/
0291
0292 inline int U4Mesh::NumberOfRotationSteps_solidName_STARTING( const char* _solidName )
0293 {
0294     int num = 0 ;
0295     for(int idx=0 ; idx < 3 ; idx++)
0296     {
0297         std::string pfx_ekey = FormEKey(U4Mesh__NumberOfRotationSteps, "solidName", "STARTING_pfx", idx );
0298         std::string val_ekey = FormEKey(U4Mesh__NumberOfRotationSteps, "solidName", "STARTING_val", idx );
0299         bool has_pfx = ssys::hasenv_(pfx_ekey.c_str()) ;
0300         bool has_val = ssys::hasenv_(val_ekey.c_str()) ;
0301         bool has_both = has_pfx && has_val ;
0302         if(!has_both) continue ;
0303         std::string pfx = ssys::getenvvar( pfx_ekey.c_str(), "" );\
0304         const char* q = pfx.c_str();
0305         bool match = sstr::StartsWith(_solidName,q);
0306         if(match)
0307         {
0308             num  = ssys::getenvint( val_ekey.c_str(), 0 );
0309             break ;
0310         }
0311     }
0312     return num ;
0313 }
0314
0315
0316
0317 inline G4Polyhedron* U4Mesh::CreatePolyhedron(const G4VSolid* solid, int numberOfRotationSteps )  // static
0318 {
0319     if(numberOfRotationSteps > 0) G4Polyhedron::SetNumberOfRotationSteps(numberOfRotationSteps);
0320     G4Polyhedron* _poly = solid->CreatePolyhedron();
0321     G4Polyhedron::SetNumberOfRotationSteps(24);
0322     return _poly ;
0323 }
0324
0325
0326 inline U4Mesh::U4Mesh(const G4VSolid* solid_):
0327     solid(solid_),
0328     entityType(EType(solid)),
0329     solidName(SolidName(solid)),
0330     numberOfRotationSteps(NumberOfRotationSteps(entityType,solidName)),
0331     poly(CreatePolyhedron(solid, numberOfRotationSteps )),
0332     errvtx(0),
0333     do_init_vtx_disqualify(false),
0334     nv_poly(poly->GetNoVertices()),
0335     nv(0),
0336     nf(poly->GetNoFacets()),
0337     vtx(nullptr),
0338     vtx_(nullptr),
0339     fpd(nullptr)
0340    ,face(NP::Make<int>(nf,4)),
0341     face_(face->values<int>()),
0342     nt(0),
0343     tri(nullptr),
0344     tri_(nullptr),
0345     tpd(nullptr)
0346 {
0347     init();
0348 }
0349
0350 inline void U4Mesh::init()
0351 {
0352     init_vtx_face_count();
0353
0354     init_vtx() ;
0355     init_fpd() ;
0356
0357     init_face() ;
0358     init_tri() ;
0359     init_tpd() ;
0360
0361 }
0362
0363 /**
0364 U4Mesh::init_vtx_face_count
0365 ----------------------------
0366
0367 For each vtx count how many faces it is referenced from
0368 and keep that in a map. This is used to
0369 exclude vertices that are not referenced
0370 by any facet in the output vtx array.
0371
0372 This is needed as find that the polygonization of G4Orb
0373 via G4PolyhedronSphere includes a vertex (0,0,0)
0374 that is not referenced from any facet which causes
0375 issues for generation of normals, yeilding (nan,nan,nan)
0376 due to the attempt at smooth normalization.
0377
0378 But the disqualify skips necessary verts for polycone ?
0379
0380
0381 **/
0382 inline void U4Mesh::init_vtx_face_count()
0383 {
0384     for(int i=0 ; i < nf ; i++)
0385     {
0386         G4int nedge;
0387         G4int ivertex[4];
0388         G4int iedgeflag[4];
0389         G4int ifacet[4];
0390         poly->GetFacet(i+1, nedge, ivertex, iedgeflag, ifacet);
0391         assert( nedge == 3 || nedge == 4  );
0392
0393         for(int j=0 ; j < nedge ; j++)
0394         {
0395             G4int iv = ivertex[j] - 1 ;
0396             if(v2fc.count(iv) == 0)
0397             {
0398                v2fc[iv] = 1 ;
0399             }
0400             else
0401             {
0402                 v2fc[iv] += 1 ;
0403             }
0404         }
0405     }
0406
0407     nv = do_init_vtx_disqualify ? v2fc.size() : nv_poly ;
0408     //std::cout << desc_vtx_face_count() ;
0409 }
0410
0411
0412 /**
0413 U4Mesh::getVertexFaceCount
0414 ---------------------------
0415
0416 Returns the number of faces that use the provided 0-based vertex index.
0417 Note that the count is obtained from the original triangles and quads,
0418 so will not exactly match the count derived after splitting quads into triangles.
0419
0420 **/
0421
0422 inline int U4Mesh::getVertexFaceCount(int iv) const
0423 {
0424     return v2fc.count(iv) == 1 ? v2fc.at(iv) : 0  ; // map can only have 0 or 1 occurrences of iv key
0425 }
0426
0427 inline std::string U4Mesh::desc_vtx_face_count() const
0428 {
0429     std::stringstream ss ;
0430     ss << "U4Mesh::desc_vtx_face_count" << std::endl ;
0431     typedef std::map<int,int> MII ;
0432     for(MII::const_iterator it = v2fc.begin() ; it != v2fc.end() ; it++)
0433     {
0434         int iv = it->first ;
0435         int fc = it->second ;
0436         int vfc = getVertexFaceCount(iv) ;
0437         assert( fc == vfc );
0438         ss << iv << ":" << fc << std::endl ;
0439     }
0440
0441     ss << "nv_poly:" << nv_poly << "nv:" << nv << std::endl ;
0442     ss << " getVertexFaceCount(0) " << getVertexFaceCount(0) << std::endl ;
0443     ss << " getVertexFaceCount(nv_poly-1) :" << getVertexFaceCount(nv_poly-1) << std::endl ;
0444     ss << " getVertexFaceCount(nv-1)      :"  << getVertexFaceCount(nv-1) << std::endl ;
0445
0446     std::string str = ss.str() ;
0447     return str ;
0448 }
0449
0450 inline void U4Mesh::init_vtx()
0451 {
0452     vtx = NP::Make<double>(nv, 3) ;
0453     vtx_ = vtx->values<double>() ;
0454
0455     int nv_count = 0 ;
0456     for(int iv=0 ; iv < nv_poly ; iv++)
0457     {
0458         int vfc = getVertexFaceCount(iv) ;
0459         G4Point3D point = poly->GetVertex(iv+1) ;
0460         bool collect = do_init_vtx_disqualify ? vfc > 0 : true ;
0461         if( !collect )
0462         {
0463
0464             std::cout
0465                 << "U4Mesh::init_vtx "
0466                 << id()
0467                 << " : DISQUALIFIED VTX NOT INCLUDED IN ANY FACET "
0468                 << " iv " << iv
0469                 << " vfc " << vfc
0470                 << " point [ "
0471                 << point.x()
0472                 << ","
0473                 << point.y()
0474                 << ","
0475                 << point.z()
0476                 << "]"
0477                 << std::endl
0478                 ;
0479         }
0480         else
0481         {
0482             vtx_[3*nv_count+0] = point.x() ;
0483             vtx_[3*nv_count+1] = point.y() ;
0484             vtx_[3*nv_count+2] = point.z() ;
0485             nv_count += 1 ;
0486         }
0487     }
0488     if( do_init_vtx_disqualify )
0489     {
0490         assert( nv_count == nv );
0491     }
0492 }
0493
0494 /**
0495 U4Mesh::init_fpd
0496 --------------------
0497
0498 The format needed by general pyvista PolyData faces
0499 (with tri and quad for example) requires an irregular array
0500 that is not easy to form from the existing regular face array.
0501 But its trivial to form the array in C++, so here it is.
0502 For example the faces of a quad and two triangles::
0503
0504     faces = np.hstack([[4, 0, 1, 2, 3],[3, 0, 1, 4],[3, 1, 2, 4]])
0505
0506 Create pyvista PolyData from the fpd and vtx array with the below::
0507
0508     pd = pv.PolyData(f.vtx, f.fpd)
0509
0510 For use of that in context see u4/tests/U4Mesh_test.py
0511
0512 **/
0513 inline void U4Mesh::init_fpd()
0514 {
0515     std::vector<int> _fpd ;
0516     for(int i=0 ; i < nf ; i++)
0517     {
0518         G4int nedge;
0519         G4int ivertex[4];
0520         G4int iedgeflag[4];
0521         G4int ifacet[4];
0522         poly->GetFacet(i+1, nedge, ivertex, iedgeflag, ifacet);
0523         assert( nedge == 3 || nedge == 4  );
0524
0525         _fpd.push_back(nedge);
0526         for(int j=0 ; j < nedge ; j++)
0527         {
0528             int jvtx = ivertex[j]-1 ;
0529             if( jvtx >= nv ) std::cout
0530                 << "U4Mesh::init_fpd "
0531                 << id()
0532                 << " ERR jvtx >= nv  (j, jvtx, nv, nv_poly ) "
0533                 << "(" << j << "," << jvtx << "," << nv << "," << nv_poly << ")"
0534                 << std::endl
0535                 ;
0536             if( jvtx >= nv ) errvtx += 1 ;
0537             _fpd.push_back(jvtx);
0538         }
0539     }
0540     fpd = NPX::Make<int>(_fpd);
0541 }
0542
0543 inline NPFold* U4Mesh::serialize() const
0544 {
0545     NPFold* fold = new NPFold ;
0546     fold->add("vtx",  vtx );
0547     fold->add("fpd",  fpd );
0548
0549     fold->add("face", face );
0550     fold->add("tri",  tri );
0551     fold->add("tpd",  tpd );
0552
0553     fold->set_meta<std::string>("entityType", entityType);
0554     fold->set_meta<std::string>("solidName",  solidName);
0555     if( numberOfRotationSteps > 0 )
0556     {
0557         fold->set_meta<int>("numberOfRotationSteps",  numberOfRotationSteps );
0558     }
0559     return fold ;
0560 }
0561
0562 inline void U4Mesh::save(const char* base) const
0563 {
0564     NPFold* fold = serialize();
0565     fold->save(base);
0566 }
0567
0568 inline std::string U4Mesh::id() const
0569 {
0570     std::stringstream ss ;
0571     ss << solidName ;
0572     std::string str = ss.str();
0573     return str ;
0574 }
0575
0576 inline std::string U4Mesh::desc() const
0577 {
0578     std::stringstream ss ;
0579     ss << "U4Mesh::desc"
0580        << " solidName " << solidName
0581        << " nv " << nv
0582        << " nv_poly " << nv_poly
0583        << " nf " << nf
0584        << std::endl
0585        ;
0586     std::string str = ss.str();
0587     return str ;
0588 }
0589
0590
0591 /**
0592 U4Mesh::init_face
0593 -------------------
0594
0595 * G4Polyhedron providing 1-based vertex indices,
0596   convert those to more standard 0-based
0597
0598 **/
0599
0600 inline void U4Mesh::init_face()
0601 {
0602     nt = 0 ;
0603     for(int i=0 ; i < nf ; i++)
0604     {
0605         G4int nedge;
0606         G4int ivertex[4];
0607         G4int iedgeflag[4];
0608         G4int ifacet[4];
0609         poly->GetFacet(i+1, nedge, ivertex, iedgeflag, ifacet);
0610         assert( nedge == 3 || nedge == 4  );
0611
0612         face_[i*4+0] = ivertex[0]-1 ;
0613         face_[i*4+1] = ivertex[1]-1 ;
0614         face_[i*4+2] = ivertex[2]-1 ;
0615         face_[i*4+3] = nedge == 4 ? ivertex[3]-1 : -1 ;
0616
0617         switch(nedge)  // count triangles for init_tri
0618         {
0619             case 3: nt += 1 ; break ;
0620             case 4: nt += 2 ; break ;
0621         }
0622     }
0623 }
0624
0625
0626 /**
0627 U4Mesh::init_tri
0628 -----------------
0629
0630 Quads are split into two tri::
0631
0632            w-------z
0633            |       |
0634            |       |
0635            |       |
0636            x-------y
0637
0638            +-------z
0639            |     / |
0640            |   /   |   1:(x->y->z)   0,1,2
0641            | /  (1)|
0642            x-------y
0643
0644            w-------z
0645            | (2) / |
0646            |   /   |   2:(x->z->w)  0,2,3
0647            | /     |
0648            x-------+
0649
0650
0651 HMM: this assumes repeatability of GetFacet ?
0652
0653 * G4Polyhedron provides 1-based vertex indices,
0654   convert those to much more standard 0-based
0655
0656 Note that using tri from pyvista is not so convenient
0657 as have to change format to that needed by pv.PolyData::
0658
0659     def PolyData_FromTRI(f):
0660         tri = np.zeros( (len(f.tri), 4 ), dtype=np.int32 )
0661         tri[:,0] = 3
0662         tri[:,1:] = f.tri
0663         pd = pv.PolyData(f.vtx, tri)
0664         return pd
0665
0666
0667 HMM: JUST RANDOMLY DOING THIS LIABLE TO
0668 GIVE MIXED WINDING ORDER CW/CCW ?
0669
0670
0671 **/
0672
0673 inline void U4Mesh::init_tri()
0674 {
0675     assert( nt > 0 );
0676     tri = NP::Make<int>(nt, 3 );
0677     tri_ = tri->values<int>() ;
0678     int jt = 0 ;
0679
0680     for(int i=0 ; i < nf ; i++)
0681     {
0682         G4int nedge;
0683         G4int ivertex[4];
0684         G4int iedgeflag[4];
0685         G4int ifacet[4];
0686         poly->GetFacet(i+1, nedge, ivertex, iedgeflag, ifacet);
0687         assert( nedge == 3 || nedge == 4  );
0688
0689         if( nedge == 3 )
0690         {
0691             assert( jt < nt );
0692             for(int j=0 ; j < 3 ; j++) tri_[jt*3+j] = ivertex[j]-1 ;
0693             jt += 1 ;
0694         }
0695         else if( nedge == 4 )
0696         {
0697             // adhoc splitting without regard for how the
0698             // vertices are positioned in space ?
0699             // will that messup the CW/CCW winding order
0700
0701             assert( jt < nt );
0702             tri_[jt*3+0] = ivertex[0]-1 ; // x
0703             tri_[jt*3+1] = ivertex[1]-1 ; // y
0704             tri_[jt*3+2] = ivertex[2]-1 ; // z
0705             jt += 1 ;
0706
0707             assert( jt < nt );
0708             tri_[jt*3+0] = ivertex[0]-1 ; // x
0709             tri_[jt*3+1] = ivertex[2]-1 ; // z
0710             tri_[jt*3+2] = ivertex[3]-1 ; // w
0711
0712             jt += 1 ;
0713         }
0714     }
0715     assert( nt == jt );
0716 }
0717
0718 inline void U4Mesh::init_tpd()
0719 {
0720     std::vector<int> _tpd ;
0721     for(int i=0 ; i < nf ; i++)
0722     {
0723         G4int nedge;
0724         G4int ivertex[4];
0725         G4int iedgeflag[4];
0726         G4int ifacet[4];
0727         poly->GetFacet(i+1, nedge, ivertex, iedgeflag, ifacet);
0728         assert( nedge == 3 || nedge == 4  );
0729
0730         if( nedge == 3 )
0731         {
0732             _tpd.push_back(3) ;
0733             for(int j=0 ; j < 3 ; j++) _tpd.push_back(ivertex[j]-1);
0734         }
0735         else if( nedge == 4 )
0736         {
0737             _tpd.push_back(3) ;
0738             _tpd.push_back(ivertex[0]-1) ;  // x
0739             _tpd.push_back(ivertex[1]-1) ;  // y
0740             _tpd.push_back(ivertex[2]-1) ;  // z
0741
0742             _tpd.push_back(3) ;
0743             _tpd.push_back(ivertex[0]-1);  // x
0744             _tpd.push_back(ivertex[2]-1);  // z
0745             _tpd.push_back(ivertex[3]-1);  // w
0746         }
0747     }
0748     tpd = NPX::Make<int>(_tpd);
0749 }
0750
0751