eic-opticks/sysrap/SBnd.h

0001 #pragma once
0002 /**
0003 SBnd.h : Used to fish material properties out of the boundary buffer
0004 =======================================================================
0005
0006 NB: this only works with the standard set of 8 properties, it
0007 does not work with scintillator properties as those are not
0008 stored in the boundary buffer.
0009
0010 Principal user QBnd.hh
0011
0012 * SSim has several functions that perhaps could be relocated here
0013 * Note that SLOG logging machinery doesnt work with header only imps
0014
0015 **/
0016
0017 #include <cassert>
0018 #include <csignal>
0019 #include <string>
0020 #include <vector>
0021 #include <array>
0022 #include <sstream>
0023 #include <set>
0024
0025 #include "NP.hh"
0026 #include "sstr.h"
0027 #include "sdigest.h"
0028 #include "sproplist.h"
0029 #include "sidxname.h"
0030
0031
0032
0033 struct SBnd
0034 {
0035     const NP* bnd ;
0036
0037     //static constexpr const unsigned MISSING = ~0u ;
0038     static constexpr const int MISSING = -1 ;
0039
0040     static constexpr const char* EXAMPLE = "Pyrex/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/Vacuum" ;
0041
0042     const std::vector<std::string>& bnames ;
0043
0044     SBnd(const NP* bnd_);
0045
0046     std::string getItemDigest( int i, int j, int w=8 ) const ;
0047     std::string descBoundary() const ;
0048     int getNumBoundary() const ;
0049     const char* getBoundarySpec(int idx) const ;
0050     void        getBoundarySpec(std::vector<std::string>& names, const int* idx , int num_idx ) const ;
0051
0052     int    getBoundaryIndex(const char* spec) const ;
0053
0054     static constexpr const char* getBoundaryIndices_NOTES = R"(
0055 SBnd::getBoundaryIndices_NOTES
0056 --------------------------------
0057
0058 Failure to look up a boundary index indicates the spec string
0059 specifies a boundary that is not present within the current geometry.
0060 For example check the spec is present in the bnd_names.txt of the bnd array
0061 be used::
0062
0063     $HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/standard/bnd.npy
0064     $HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/standard/bnd_names.txt
0065
0066
0067 )" ;
0068
0069     bool        hasAllBoundaryIndices( const char* bnd_sequence, char delim=',' ) const ;
0070     void        getBoundaryIndices( std::vector<int>& bnd_idx, const char* bnd_sequence, char delim=',' ) const ;
0071     std::string descBoundaryIndices( const std::vector<int>& bnd_idx ) const ;
0072
0073     int    getBoundaryLine(const char* spec, int j) const ;
0074     static int GetMaterialLine( const char* material, const std::vector<std::string>& specs );
0075     int    getMaterialLine( const char* material ) const ;
0076
0077     static std::string DescDigest(const NP* bnd, int w=16) ;
0078
0079     std::string desc() const ;
0080
0081     void getMaterialNames( std::vector<std::string>& names ) const ;
0082     static std::string DescNames( std::vector<std::string>& names ) ;
0083
0084
0085
0086     bool findName( int& i, int& j, const char* qname ) const ;
0087     static bool FindName( int& i, int& j, const char* qname, const std::vector<std::string>& names ) ;
0088
0089     NP* getPropertyGroup(const char* qname, int k=-1) const ;
0090
0091     template<typename T>
0092     void getProperty(std::vector<T>& out, const char* qname, const char* propname ) const ;
0093
0094
0095     static void FillMaterialLine(
0096          std::vector<int>& mtline,
0097          const std::vector<int>& mtindex,
0098          const std::vector<std::string>& mtname,
0099          const std::vector<std::string>& specs
0100     );
0101
0102
0103     static void GetSpecsFromString( std::vector<std::string>& specs , const char* specs_, char delim='\n' );
0104
0105
0106
0107     NP* bd_from_optical(const NP* op ) const ;
0108     NP* mat_from_bd(const NP* bd) const ;
0109     NP* reconstruct_sur() const ;
0110
0111 };
0112
0113
0114 inline SBnd::SBnd(const NP* bnd_)
0115     :
0116     bnd(bnd_),
0117     bnames(bnd->names)
0118 {
0119     int num_bnames = bnames.size() ;
0120     if( num_bnames == 0 ) std::cerr << "SBnd::SBnd no names from bnd " << ( bnd ? bnd->sstr() : "-" ) << std::endl ;
0121     //assert(num_bnames > 0 );
0122 }
0123
0124 inline std::string SBnd::getItemDigest( int i, int j, int w ) const
0125 {
0126     return sdigest::Item(bnd, i, j, w );   // formerly SSim::GetItemDigest
0127 }
0128 inline std::string SBnd::descBoundary() const
0129 {
0130     std::stringstream ss ;
0131     for(int i=0 ; i < int(bnames.size()) ; i++)
0132        ss << std::setw(2) << i << " " << bnames[i] << std::endl ;
0133     std::string s = ss.str();
0134     return s ;
0135 }
0136
0137
0138 inline int SBnd::getNumBoundary() const
0139 {
0140     return bnames.size();
0141 }
0142 inline const char* SBnd::getBoundarySpec(int idx) const
0143 {
0144     assert( idx < int(bnames.size()) );
0145     const std::string& s = bnames[idx];
0146     return s.c_str();
0147 }
0148 inline void SBnd::getBoundarySpec(std::vector<std::string>& names, const int* idx , int num_idx ) const
0149 {
0150     for(int i=0 ; i < num_idx ; i++)
0151     {
0152         int index = idx[i] ;
0153         const char* spec = getBoundarySpec(index);   // 0-based
0154         names.push_back(spec);
0155     }
0156 }
0157
0158
0159 inline int SBnd::getBoundaryIndex(const char* spec) const
0160 {
0161     int idx = MISSING ;
0162     for(int i=0 ; i < int(bnames.size()) ; i++)
0163     {
0164         if(spec && strcmp(bnames[i].c_str(), spec) == 0)
0165         {
0166             idx = i ;
0167             break ;
0168         }
0169     }
0170     return idx ;
0171 }
0172
0173
0174 inline bool SBnd::hasAllBoundaryIndices( const char* bnd_sequence, char delim ) const
0175 {
0176     std::vector<std::string> bnd ;
0177     sstr::SplitTrim(bnd_sequence,delim, bnd );
0178     int missing = 0 ;
0179
0180     for(int i=0 ; i < int(bnd.size()) ; i++)
0181     {
0182         const char* spec = bnd[i].c_str();
0183         if(strlen(spec) == 0) continue ;
0184         int bidx = getBoundaryIndex(spec);
0185         if( bidx == MISSING ) missing += 1 ;
0186     }
0187     return missing == 0 ;
0188 }
0189
0190 inline void SBnd::getBoundaryIndices( std::vector<int>& bnd_idx, const char* bnd_sequence, char delim ) const
0191 {
0192     assert( bnd_idx.size() == 0 );
0193
0194     std::vector<std::string> bnd ;
0195     sstr::SplitTrim(bnd_sequence,delim, bnd );
0196
0197     for(int i=0 ; i < int(bnd.size()) ; i++)
0198     {
0199         const char* spec = bnd[i].c_str();
0200         if(strlen(spec) == 0) continue ;
0201         int bidx = getBoundaryIndex(spec);
0202         if( bidx == MISSING ) std::cerr
0203             << " SBnd::getBoundaryIndices "
0204             << " i " << i << " invalid spec [" << spec << "]"
0205             << std::endl
0206             << getBoundaryIndices_NOTES
0207             << std::endl
0208             ;
0209         assert( bidx != MISSING );
0210         bnd_idx.push_back(bidx) ;
0211     }
0212 }
0213
0214 inline std::string SBnd::descBoundaryIndices( const std::vector<int>& bnd_idx ) const
0215 {
0216     int num_bnd_idx = bnd_idx.size() ;
0217     std::stringstream ss ;
0218     for(int i=0 ; i < num_bnd_idx ; i++)
0219     {
0220         int bidx = bnd_idx[i] ;
0221         const char* spec = getBoundarySpec(bidx);
0222         ss
0223             << " i " << std::setw(3) << i
0224             << " bidx " << std::setw(3) << bidx
0225             << " spec " << spec
0226             << std::endl
0227             ;
0228     }
0229     std::string s = ss.str();
0230     return s ;
0231 }
0232
0233 /**
0234 SBnd::getBoundaryLine
0235 -----------------------
0236
0237 The boundary spec allows to obtain the boundary index,
0238 the boundary line returned is : 4*boundary_index + j
0239
0240 **/
0241
0242 inline int SBnd::getBoundaryLine(const char* spec, int j) const
0243 {
0244     int num_bnames = bnames.size() ;
0245     int idx = getBoundaryIndex(spec);
0246     bool is_missing = idx == MISSING ;
0247     bool is_valid = !is_missing && idx < num_bnames ;
0248
0249     if(!is_valid)
0250     {
0251         std::cerr
0252             << " not is_valid "
0253             << " spec " << spec
0254             << " idx " << idx
0255             << " is_missing " << is_missing
0256             << " bnames.size " << bnames.size()
0257             << std::endl
0258             ;
0259     }
0260
0261     assert( is_valid );
0262     int line = 4*idx + j ;
0263     return line ;
0264 }
0265
0266 /**
0267 SBnd::GetMaterialLine
0268 ----------------------
0269
0270 The spec strings are assumed to be "/" delimited : omat/osur/isur/imat
0271 The first omat or imat line matching the *material* argument is returned.
0272
0273 **/
0274
0275 inline int SBnd::GetMaterialLine( const char* material, const std::vector<std::string>& specs ) // static
0276 {
0277     int line = MISSING ;
0278     for(int i=0 ; i < int(specs.size()) ; i++)
0279     {
0280         std::vector<std::string> elem ;
0281         sstr::Split(specs[i].c_str(), '/', elem );
0282         const char* omat = elem[0].c_str();
0283         const char* imat = elem[3].c_str();
0284
0285         if(strcmp( material, omat ) == 0 )
0286         {
0287             line = i*4 + 0 ;
0288             break ;
0289         }
0290         if(strcmp( material, imat ) == 0 )
0291         {
0292             line = i*4 + 3 ;
0293             break ;
0294         }
0295     }
0296     return line ;
0297 }
0298
0299 /**
0300 SBnd::getMaterialLine
0301 -----------------------
0302
0303 Searches the bname spec for the *material* name in omat or imat slots,
0304 returning the first found.
0305
0306 **/
0307
0308 inline int SBnd::getMaterialLine( const char* material ) const
0309 {
0310     return GetMaterialLine(material, bnames);
0311 }
0312
0313
0314
0315
0316 /**
0317 SBnd::DescDigest
0318 --------------------
0319
0320 bnd with shape (44, 4, 2, 761, 4, )::
0321
0322    ni : boundaries
0323    nj : 0:omat/1:osur/2:isur/3:imat
0324    nk : 0 or 1 property group
0325    nl : wavelengths
0326    nm : payload
0327
0328 ::
0329
0330     2022-04-20 14:53:14.544 INFO  [4031964] [test_DescDigest@133]
0331     5acc01c3 79cfae67 79cfae67 5acc01c3  Galactic///Galactic
0332     5acc01c3 79cfae67 79cfae67 8b22bf98  Galactic///Rock
0333     8b22bf98 79cfae67 79cfae67 5acc01c3  Rock///Galactic
0334     8b22bf98 79cfae67 0a5eab3f c2759ba7  Rock//Implicit_RINDEX_NoRINDEX_pDomeAir_pDomeRock/Air
0335     8b22bf98 79cfae67 79cfae67 8b22bf98  Rock///Rock
0336     8b22bf98 79cfae67 0a5eab3f c2759ba7  Rock//Implicit_RINDEX_NoRINDEX_pExpHall_pExpRockBox/Air
0337     c2759ba7 79cfae67 79cfae67 8b22bf98  Air///Steel
0338
0339 **/
0340
0341
0342 inline std::string SBnd::DescDigest(const NP* bnd, int w )  // static
0343 {
0344     int ni = bnd->shape[0] ;
0345     int nj = bnd->shape[1] ;
0346
0347     const std::vector<std::string>& names = bnd->names ;
0348     assert( int(names.size()) == ni );
0349
0350     std::stringstream ss ;
0351     ss << "SBnd::DescDigest" << std::endl ;
0352     for(int i=0 ; i < ni ; i++)
0353     {
0354         ss << std::setw(3) << i << " " ;
0355         for(int j=0 ; j < nj ; j++)
0356         {
0357             std::string dig = sdigest::Item(bnd, i, j ) ;    // formerly SDigestNP::Item
0358             std::string sdig = dig.substr(0, w);
0359             ss << std::setw(w) << sdig << " " ;
0360         }
0361         ss << " " << names[i] << std::endl ;
0362     }
0363     std::string s = ss.str();
0364     return s ;
0365 }
0366
0367
0368 inline std::string SBnd::desc() const
0369 {
0370     return DescDigest(bnd,8) ;
0371 }
0372
0373 /**
0374 SBnd::getMaterialNames
0375 -----------------------
0376
0377 HMM: name order not the original one
0378
0379 **/
0380
0381 inline void SBnd::getMaterialNames( std::vector<std::string>& names ) const
0382 {
0383     for(int i=0 ; i < int(bnames.size()) ; i++)
0384     {
0385         std::vector<std::string> elem ;
0386         sstr::Split(bnames[i].c_str(), '/', elem );
0387         const char* omat = elem[0].c_str();
0388         const char* imat = elem[3].c_str();
0389
0390         if(std::find(names.begin(), names.end(), omat) == names.end()) names.push_back(omat);
0391         if(std::find(names.begin(), names.end(), imat) == names.end()) names.push_back(imat);
0392     }
0393 }
0394 inline std::string SBnd::DescNames( std::vector<std::string>& names )
0395 {
0396     std::stringstream ss ;
0397     for(int i=0 ; i < int(names.size()) ; i++) ss << names[i] << std::endl ;
0398     std::string s = ss.str();
0399     return s ;
0400 }
0401
0402
0403
0404 /**
0405 SBnd::findName
0406 ----------------
0407
0408 Returns the first (i,j)=(bidx,species) with element name
0409 matching query name *qname*.
0410
0411 bidx
0412     0-based boundary index
0413
0414 species
0415     0,1,2,3 for omat/osur/isur/imat
0416
0417 **/
0418
0419 inline bool SBnd::findName( int& i, int& j, const char* qname ) const
0420 {
0421     return FindName(i, j, qname, bnames ) ;
0422 }
0423
0424 inline bool SBnd::FindName( int& i, int& j, const char* qname, const std::vector<std::string>& names )
0425 {
0426     i = -1 ;
0427     j = -1 ;
0428     for(int b=0 ; b < int(names.size()) ; b++)
0429     {
0430         std::vector<std::string> elem ;
0431         sstr::Split(names[b].c_str(), '/', elem );
0432
0433         for(int s=0 ; s < 4 ; s++)
0434         {
0435             const char* name = elem[s].c_str();
0436             if(strcmp(name, qname) == 0 )
0437             {
0438                 i = b ;
0439                 j = s ;
0440                 return true ;
0441             }
0442         }
0443     }
0444     return false ;
0445 }
0446
0447 /**
0448 SBnd::getPropertyGroup
0449 ------------------------
0450
0451 Returns an array of material or surface properties selected by *qname* eg "Water".
0452 For example with a source bnd array of shape  (44, 4, 2, 761, 4, )
0453 the expected spawned property group  array shape depends on the value of k:
0454
0455 k=-1
0456      (2, 761, 4,)   eight property values across wavelength domain
0457 k=0
0458      (761, 4)       four property values across wavelength domain
0459 k=1
0460      (761, 4)
0461
0462 **/
0463
0464 inline NP* SBnd::getPropertyGroup(const char* qname, int k) const
0465 {
0466     int i, j ;
0467     bool found = findName(i, j, qname);
0468     if(!found) std::raise(SIGINT) ;
0469     assert(found);
0470     return bnd->spawn_item(i,j,k);
0471 }
0472
0473 /**
0474 SBnd::getProperty
0475 ------------------
0476
0477 ::
0478
0479     // <f8(44, 4, 2, 761, 4, )
0480     int species = 0 ;
0481     int group = 1 ;
0482     int wavelength = -1 ;
0483     int prop = 0 ;
0484
0485 **/
0486
0487 template<typename T>
0488 inline void SBnd::getProperty(std::vector<T>& out, const char* qname, const char* propname ) const
0489 {
0490     assert( sizeof(T) == bnd->ebyte );
0491
0492     int boundary, species ;
0493     bool found_qname = findName(boundary, species, qname);
0494     assert(found_qname);
0495     if(!found_qname) std::raise(SIGINT);
0496
0497     //const SBndProp* bp = FindMaterialProp(propname);
0498
0499     const sproplist* pm = sproplist::Material();
0500     const sprop* bp = pm->findProp(propname) ;
0501     assert(bp);
0502     if(!bp) std::raise(SIGINT);
0503
0504     int group = bp->group ;
0505     int prop = bp->prop  ;
0506     int wavelength = -1 ;   // slice dimension
0507
0508     bnd->slice(out, boundary, species, group, wavelength, prop );
0509 }
0510
0511
0512 /**
0513 SBnd::FillMaterialLine
0514 -----------------------
0515
0516 Used from stree::import_bnd
0517
0518 Uses the "specs" boundary name list to convert
0519 all the stree::mtname into st->mtline
0520
0521 These mtline are used to lookup material properties
0522 from the boundary texture array.
0523
0524 **/
0525
0526 inline void SBnd::FillMaterialLine(
0527      std::vector<int>& mtline,
0528      const std::vector<int>& mtindex,
0529      const std::vector<std::string>& mtname,
0530      const std::vector<std::string>& specs )
0531 {
0532     assert( mtindex.size() == mtname.size() );
0533     int num_mt = mtindex.size() ;
0534     mtline.clear();
0535
0536     for(int i=0 ; i < num_mt ; i++)
0537     {
0538         const char* mt = mtname[i].c_str() ;
0539         int mt_line = GetMaterialLine(mt, specs) ;  // unsigned ~0u "MISSING" becomes int -1
0540         mtline.push_back(mt_line);
0541     }
0542 }
0543
0544
0545 inline void SBnd::GetSpecsFromString( std::vector<std::string>& specs , const char* specs_, char delim )
0546 {
0547     std::stringstream ss;
0548     ss.str(specs_)  ;
0549     std::string s;
0550     while (std::getline(ss, s, delim)) if(!sstr::Blank(s.c_str())) specs.push_back(s) ;
0551     std::cout << " specs_ [" << specs_ << "] specs.size " << specs.size()  ;
0552 }
0553
0554
0555
0556 /**
0557 SBnd::bd_from_optical
0558 ----------------------
0559
0560 Hmm, but the bd is new. For reconstruction of old_mat from the old_bnd
0561 I need to use an old_bd. Can recreate that by folding first column of old_optical
0562 and subtracting 1.::
0563
0564     bd = np.array( t.old_optical[:,0].reshape(-1,4), dtype=np.int32 ) - 1
0565
0566 **/
0567
0568 inline NP* SBnd::bd_from_optical(const NP* op ) const
0569 {
0570
0571     /*
0572     //old optical was unsigned:( num_bd*4, 4 )
0573     assert( op && op->uifc == 'u' && op->shape.size() == 2 && op->shape[1] == 4 ) ;
0574     int num_op = op->shape[0] ;
0575     assert( num_op % 4 == 0 );
0576     int num_bd = num_op / 4 ;
0577     */
0578
0579     // new optical is int:(num_bd, 4, 4)
0580     assert( op && op->uifc == 'i' && op->shape.size() == 3 && op->shape[1] == 4 && op->shape[2] == 4 ) ;
0581     int num_bd = op->shape[0] ;
0582
0583
0584     const int* op_v = op->cvalues<int>() ;
0585     int ni = num_bd ;
0586     int nj = 4 ;
0587
0588     NP* bd = NP::Make<int>(ni, nj) ;
0589     int* bd_v = bd->values<int>() ;
0590
0591     for(int i=0 ; i < ni ; i++)
0592     for(int j=0 ; j < nj ; j++)
0593     {
0594         int src_index = (i*nj + j)*4 ;
0595         int dst_index = i*nj + j ;
0596         bd_v[dst_index] = int(op_v[src_index]) - 1 ;
0597     }
0598     bd->set_names(bnd->names) ;
0599     return bd ;
0600 }
0601
0602
0603 /**
0604 SBnd::mat_from_bd
0605 -----------------------
0606
0607 Creating mat from bd obtained from optical 1st column and bnd array
0608 is of course the wrong order. Typically the converse is done
0609 the bnd is created by interleaving together the mat and sur arrays.
0610 However in the old GGeo/X4 workflow the mat and sur arrays were not
0611 persisted with the bnd being created directly.
0612
0613 Hence to facilitate development of mat,sur and bnd in new workflow
0614 it is useful to reconstruct what the old workflow mat and sur would
0615 actually have been by pulling apart the old bnd and putting it together
0616 to form what the missing old mat and sur would have been.
0617
0618 When this is applied to an old bnd the old mat it provides can
0619 be compared with new mat from stree::create_mat which uses sstandard::mat
0620
0621 No gaps for materials, all mt indices in the bd int4 vector form contiguous range::
0622
0623     In [19]: np.unique( np.hstack( (t.bd[:,0], t.bd[:,3]) ) )
0624     Out[19]: array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19], dtype=int32)
0625
0626 **/
0627
0628 inline NP* SBnd::mat_from_bd(const NP* bd) const
0629 {
0630     // 0. check bd array is as expected
0631
0632     assert( bd && bd->uifc == 'i' && bd->shape.size() == 2 && bd->shape[1] == 4 ) ;
0633     assert( int(bd->names.size()) == bd->shape[0] );
0634     const int* bd_v = bd->cvalues<int>();
0635
0636     // 1. check consistency between bd (num_bd, 4) int pointers
0637     //    and bnd data eg (53, 4, 2, 761, 4, )
0638
0639     const std::vector<NP::INT>& bnd_sh = bnd->shape ;
0640     assert( bnd && bnd->uifc == 'f' && bnd->ebyte == 8 );
0641     assert( int(bnd_sh.size()) == 5 ) ;
0642     assert( bnd_sh[0] == bd->shape[0] );
0643     const double* bnd_v = bnd->cvalues<double>() ;
0644
0645
0646     // 2. first bd pass to find the number of material indices
0647
0648     NP::INT ni = bd->shape[0] ;
0649     NP::INT nj = 4 ;
0650
0651     std::set<sidxname, sidxname_ordering> mm ;
0652     for(NP::INT i=0 ; i < ni ; i++)
0653     {
0654         int omat = bd_v[i*nj+0] ;
0655         int imat = bd_v[i*nj+3] ;
0656
0657         const char* bdn = bd->names[i].c_str() ;
0658         std::vector<std::string> elem ;
0659         sstr::Split(bdn, '/', elem );
0660         const char* omat_ = elem[0].c_str();
0661         const char* imat_ = elem[3].c_str();
0662
0663         sidxname om(omat,omat_) ;
0664         sidxname im(imat,imat_) ;
0665
0666         mm.insert(om);
0667         mm.insert(im);
0668     }
0669
0670     std::vector<sidxname> vmm(mm.begin(),mm.end()) ;
0671     int num_mt = vmm.size() ;
0672
0673     // 3. assert that bd (omat,imat) contains contiguous set of material indices
0674
0675     for(int i=0 ; i < num_mt ; i++) assert( i == vmm[i].idx ) ;
0676     for(int i=0 ; i < num_mt ; i++) std::cout << vmm[i].desc() << std::endl;
0677
0678     std::vector<std::string> names ;
0679     for(int i=0 ; i < num_mt ; i++)
0680     {
0681         names.push_back( vmm[i].name ) ;  // HMM null termination ?
0682     }
0683
0684     std::cout
0685         << "SBnd::mat_from_bd"
0686         << " bd->names.size " << bd->names.size()
0687         << " bnd->names.size " << bnd->names.size()
0688         << " bnd->shape[0] " << bnd->shape[0]
0689         << " bnd->sstr() " << bnd->sstr()
0690         << " ni " << ni
0691         << std::endl
0692         ;
0693
0694
0695     // 4. prep mat array
0696
0697     int np = bnd_sh[2]*bnd_sh[3]*bnd_sh[4] ; // num payload values for mat (or sur)
0698
0699     NP* mat = NP::Make<double>(num_mt,bnd_sh[2], bnd_sh[3], bnd_sh[4] ) ;
0700     mat->set_names(names) ;
0701
0702     double* mat_v = mat->values<double>();
0703
0704     // 4. second bd pass to populate reconstructed mat array
0705     //    (note that each mat may be written multiple times
0706     //     but that doesnt matter as all the same)
0707
0708     for(int i=0 ; i < ni ; i++)
0709     {
0710         int omat = bd_v[i*nj+0] ;
0711         assert( omat > -1 && omat < num_mt );
0712         for(int p=0 ; p < np ; p++) mat_v[omat*np+p] = bnd_v[(i*4+0)*np+p] ;
0713
0714         int imat = bd_v[i*nj+3] ;
0715         assert( imat > -1 && imat < num_mt );
0716         for(int p=0 ; p < np ; p++) mat_v[imat*np+p] = bnd_v[(i*4+3)*np+p] ;
0717     }
0718
0719     return mat ;
0720 }
0721
0722
0723 /**
0724 Quite a few gaps for surfaces, not all surfaces are referenced from the bd int4 vector::
0725
0726     In [20]: np.unique( np.hstack( (t.bd[:,1], t.bd[:,2]) ) )
0727     Out[20]: array([-1,  0,  1,  2,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 33, 34, 35, 36, 37, 38, 39], dtype=int32)
0728
0729 **/
0730
0731 inline NP* SBnd::reconstruct_sur() const
0732 {
0733     return nullptr ;
0734 }
0735