Initial Implementation of virtual lattice, with contributions by Li Ruihan

include/openmc/cell.h

@@ -60,6 +60,7 @@ class Universe {
 public:
   int32_t id_;            //!< Unique ID
   vector<int32_t> cells_; //!< Cells within this universe
+  int filled_with_triso_base_ = -1;
 
   //! \brief Write universe information to an HDF5 group.
   //! \param group_id An HDF5 group id.
@@ -212,6 +213,15 @@ public:
   int32_t fill_;            //!< Universe # filling this cell
   int32_t n_instances_ {0}; //!< Number of instances of this cell
   GeometryType geom_type_;  //!< Geometric representation type (CSG, DAGMC)
+  bool virtual_lattice_;    //!< If the cell is the base of a virtual triso lattice
+  bool triso_particle_;
+
+
+  // Specification of the viryual lattice
+  vector<double> vl_lower_left_;
+  vector<double> vl_pitch_;
+  vector<int32_t> vl_shape_;
+  vector<vector<int32_t>> vl_triso_distribution_;
 
   //! \brief Index corresponding to this cell in distribcell arrays
   int distribcell_index_ {C_NONE};

include/openmc/surface.h

@@ -88,6 +88,9 @@ public:
   std::shared_ptr<BoundaryCondition> bc_ {nullptr}; //!< Boundary condition
   GeometryType geom_type_;   //!< Geometry type indicator (CSG or DAGMC)
   bool surf_source_ {false}; //!< Activate source banking for the surface?
+  int triso_base_index_;
+  int triso_particle_index_ = -1;
+  bool is_triso_surface_ = false;
 
   explicit Surface(pugi::xml_node surf_node);
   Surface();
@@ -126,6 +129,11 @@ public:
   //!   exactly on the surface.
   virtual double distance(Position r, Direction u, bool coincident) const = 0;
 
+  virtual bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const { return {}; };
+  virtual void connect_to_triso_base(int triso_index, std::string key) {};
+  virtual vector<double> get_center() const { return {}; };
+  virtual double get_radius() const { return {}; };
+
   //! Compute the local outward normal direction of the surface.
   //! \param r A 3D Cartesian coordinate.
   //! \return Normal direction
@@ -164,6 +172,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
   BoundingBox bounding_box(bool pos_side) const;
 
   double x0_;
@@ -182,6 +191,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
   BoundingBox bounding_box(bool pos_side) const;
 
   double y0_;
@@ -200,6 +210,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
   BoundingBox bounding_box(bool pos_side) const;
 
   double z0_;
@@ -218,6 +229,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double A_, B_, C_, D_;
 };
@@ -237,6 +249,7 @@ public:
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
   BoundingBox bounding_box(bool pos_side) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double y0_, z0_, radius_;
 };
@@ -256,6 +269,7 @@ public:
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
   BoundingBox bounding_box(bool pos_side) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double x0_, z0_, radius_;
 };
@@ -275,6 +289,7 @@ public:
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
   BoundingBox bounding_box(bool pos_side) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double x0_, y0_, radius_;
 };
@@ -291,11 +306,16 @@ public:
   explicit SurfaceSphere(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
+  vector<double> get_center() const;
+  double get_radius() const;
+  void connect_to_triso_base(int triso_index, std::string key);
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
   BoundingBox bounding_box(bool pos_side) const;
 
   double x0_, y0_, z0_, radius_;
+  //int triso_base_index_ = -1;
 };
 
 //==============================================================================
@@ -312,6 +332,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double x0_, y0_, z0_, radius_sq_;
 };
@@ -330,6 +351,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double x0_, y0_, z0_, radius_sq_;
 };
@@ -348,6 +370,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double x0_, y0_, z0_, radius_sq_;
 };
@@ -366,6 +389,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
   double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;
@@ -384,6 +408,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double x0_, y0_, z0_, A_, B_, C_;
 };
@@ -401,6 +426,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
 
   double x0_, y0_, z0_, A_, B_, C_;
 };
@@ -417,6 +443,7 @@ public:
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
+  bool triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const;
   void to_hdf5_inner(hid_t group_id) const;
 
   double x0_, y0_, z0_, A_, B_, C_;

openmc/cell.py

@@ -112,6 +112,12 @@ class Cell(IDManagerMixin):
         self._num_instances = None
         self._volume = None
         self._atoms = None
+        self._triso_particle = False
+        self.virtual_lattice = False
+        self.lower_left = None
+        self.pitch = None
+        self.shape = None
+
 
     def __contains__(self, point):
         if self.region is None:
@@ -572,6 +578,14 @@ class Cell(IDManagerMixin):
         """
         element = ET.Element("cell")
         element.set("id", str(self.id))
+        if self._triso_particle:
+            element.set("triso_particle", 'true')
+        if self.virtual_lattice:
+            element.set("virtual_lattice", str(self.virtual_lattice))
+            element.set("lower_left", ' '.join(map(str, self.lower_left)))
+            element.set("pitch", ' '.join(map(str, self.pitch)))
+            element.set("shape", ' '.join(map(str, self.shape)))
+
 
         if len(self._name) > 0:
             element.set("name", str(self.name))

openmc/model/triso.py

@@ -55,6 +55,7 @@ class TRISO(openmc.Cell):
 
     def __init__(self, outer_radius, fill, center=(0., 0., 0.)):
         self._surface = openmc.Sphere(r=outer_radius)
+        self._triso_particle = False
         super().__init__(fill=fill, region=-self._surface)
         self.center = np.asarray(center)
 
@@ -800,7 +801,7 @@ class _SphericalShell(_Container):
             q[:] = (q - c)*ll[0]/r + c
 
 
-def create_triso_lattice(trisos, lower_left, pitch, shape, background):
+def create_triso_lattice(trisos, lower_left, pitch, shape, background, virtual=False):
     """Create a lattice containing TRISO particles for optimized tracking.
 
     Parameters
@@ -824,9 +825,16 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background):
 
     """
 
+    if virtual:
+        real_pitch = copy.deepcopy(pitch)
+        real_shape = copy.deepcopy(shape)
+        pitch = [real_pitch[i]*real_shape[i] for i in range(len(real_pitch))]
+        shape = [1 for i in range(len(real_shape))]
+    
     lattice = openmc.RectLattice()
     lattice.lower_left = lower_left
     lattice.pitch = pitch
+    
 
     indices = list(np.broadcast(*np.ogrid[:shape[2], :shape[1], :shape[0]]))
     triso_locations = {idx: [] for idx in indices}
@@ -843,6 +851,8 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background):
                                                 y0=t._surface.y0,
                                                 z0=t._surface.z0)
                 t_copy.region = -t_copy._surface
+                if virtual:
+                    t_copy._triso_particle = True
                 triso_locations[idx].append(t_copy)
             else:
                 warnings.warn('TRISO particle is partially or completely '
@@ -857,6 +867,13 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background):
         else:
             background_cell = openmc.Cell(fill=background)
 
+        if virtual:
+            background_cell.virtual_lattice = True
+            background_cell.pitch = real_pitch
+            background_cell.shape = real_shape
+            background_cell.lower_left = [-pitch[i]/2 for i in range(len(pitch))]
+            
+
         u = openmc.Universe()
         u.add_cell(background_cell)
         for t in triso_list:

src/cell.cpp

@@ -109,6 +109,62 @@ vector<int32_t> tokenize(const std::string region_spec)
   return tokens;
 }
 
+vector<vector<int32_t>> \
+  generate_triso_distribution(\
+    vector<int> lattice_shape,\
+    vector<double> lattice_pitch,\
+    vector<double> lattice_lower_left,\
+    vector<std::int32_t> cell_rpn, int id)
+{
+  vector<vector<int32_t>> triso_distribution(lattice_shape[0]*lattice_shape[1]*lattice_shape[2]);
+  vector<double> mesh_center(3);
+  vector<int> mesh_ind(3);
+  /*
+  for (int i=0; i<lattice_shape[0]; i++) {
+    for (int j=0; j<lattice_shape[1]; j++) {
+      for (int k=0; k<lattice_shape[2]; k++) {
+        mesh_center[0]=(i+0.5)*lattice_pitch[0]+lattice_lower_left[0];
+        mesh_center[1]=(j+0.5)*lattice_pitch[1]+lattice_lower_left[1];
+        mesh_center[2]=(k+0.5)*lattice_pitch[2]+lattice_lower_left[2];
+        for (int32_t token : cell_rpn) {
+          if (token >= OP_UNION) continue;
+          if (model::surfaces[abs(token) - 1]->triso_in_mesh(mesh_center, lattice_pitch)) {
+            triso_distribution[i+j*lattice_shape[0]+k*lattice_shape[0]*lattice_shape[1]].push_back(token);
+            model::surfaces[abs(token) - 1]->connect_to_triso_base(id, "base");
+          }
+        }
+      }
+    }
+  }*/
+  for (int32_t token : cell_rpn) {
+    if (token >= OP_UNION) continue;
+    vector<double> triso_center=model::surfaces[abs(token) - 1]->get_center();
+    for (int i=0; i<3; i++) {
+      mesh_ind[i]=floor((triso_center[i]-lattice_lower_left[i])/lattice_pitch[i]);
+    }
+    for (int i=mesh_ind[0]-1; i<=mesh_ind[0]+1; i++) {
+      for (int j=mesh_ind[1]-1; j<=mesh_ind[1]+1; j++) {
+        for (int k=mesh_ind[2]-1; k<=mesh_ind[2]+1; k++) {
+          if (i < 0 || i >= lattice_shape[0] ||\
+              j < 0 || j >= lattice_shape[1] ||\
+              k < 0 || k >= lattice_shape[2]) continue;
+          mesh_center[0]=(i+0.5)*lattice_pitch[0]+lattice_lower_left[0];
+          mesh_center[1]=(j+0.5)*lattice_pitch[1]+lattice_lower_left[1];
+          mesh_center[2]=(k+0.5)*lattice_pitch[2]+lattice_lower_left[2];
+          if (model::surfaces[abs(token) - 1]->triso_in_mesh(mesh_center, lattice_pitch)) {
+            triso_distribution[i+j*lattice_shape[0]+k*lattice_shape[0]*lattice_shape[1]].push_back(token);
+            model::surfaces[abs(token) - 1]->connect_to_triso_base(id, "base");
+          }
+        }
+      }
+    }
+  }
+
+
+  return triso_distribution;
+}
+
+
 //==============================================================================
 //! Convert infix region specification to Reverse Polish Notation (RPN)
 //!
@@ -211,6 +267,37 @@ void Universe::to_hdf5(hid_t universes_group) const
 
 bool Universe::find_cell(Particle& p) const
 {
+  if (filled_with_triso_base_ != -1) {
+    Cell& c {*model::cells[model::cell_map[filled_with_triso_base_]]};
+    vector<int> lat_ind(3);
+    Position r {p.r_local()};
+    lat_ind[0]=floor((r.x-c.vl_lower_left_[0])/c.vl_pitch_[0]);
+    lat_ind[1]=floor((r.y-c.vl_lower_left_[1])/c.vl_pitch_[1]);
+    lat_ind[2]=floor((r.z-c.vl_lower_left_[2])/c.vl_pitch_[2]);
+    int32_t i_univ = p.coord(p.n_coord() - 1).universe;
+    for (int token : c.vl_triso_distribution_[lat_ind[0]+lat_ind[1]*c.vl_shape_[0]+lat_ind[2]*c.vl_shape_[0]*c.vl_shape_[1]]) {
+      vector<double> triso_center=model::surfaces[abs(token) - 1]->get_center();
+      double triso_radius=model::surfaces[abs(token) - 1]->get_radius();
+      if (model::cells[model::cell_map[model::surfaces[abs(token) - 1]->triso_base_index_]]->universe_!= i_univ) continue;
+      if (abs(token)==abs(p.surface())) {
+        if (p.surface() < 0) {
+          p.coord(p.n_coord() - 1).cell = model::cell_map[model::surfaces[abs(token) - 1]->triso_particle_index_];
+          return true;
+        } else {
+          p.coord(p.n_coord() - 1).cell = model::cell_map[filled_with_triso_base_];
+          return true;
+        }
+      }
+      if (pow(r.x-triso_center[0],2)+pow(r.y-triso_center[1],2)+pow(r.z-triso_center[2],2) < pow(triso_radius,2)) {
+        p.coord(p.n_coord() - 1).cell = model::cell_map[model::surfaces[abs(token) - 1]->triso_particle_index_];
+        return true;
+      }
+    }
+    if (model::cells[model::cell_map[filled_with_triso_base_]]->universe_== i_univ) {
+      p.coord(p.n_coord() - 1).cell = model::cell_map[filled_with_triso_base_];
+      return true;
+    }
+  }
   const auto& cells {
     !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
 
@@ -481,6 +568,29 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
     universe_ = 0;
   }
 
+  // Check if the cell is the base of a virtual triso lattice
+  bool virtual_lattice_present = check_for_node(cell_node, "virtual_lattice");
+  if (virtual_lattice_present) {
+    virtual_lattice_ = get_node_value_bool(cell_node, "virtual_lattice");
+    if (virtual_lattice_) {
+      if (check_for_node(cell_node, "lower_left") && check_for_node(cell_node, "pitch") && check_for_node(cell_node, "shape")) {
+        vl_lower_left_ = get_node_array<double>(cell_node, "lower_left");
+        vl_pitch_ = get_node_array<double>(cell_node, "pitch");
+        vl_shape_ = get_node_array<int>(cell_node, "shape");
+      } else {
+        fatal_error(fmt::format("Lower_left, pitch and shape of the virtual lattice must be specified for cell {}", id_));
+      }
+    }
+  } else {
+    virtual_lattice_ = false;
+  }
+  
+  if (check_for_node(cell_node, "triso_particle")) {
+    triso_particle_ = get_node_value_bool(cell_node, "triso_particle");
+  } else {
+    triso_particle_ = false;
+  }
+
   // Make sure that either material or fill was specified, but not both.
   bool fill_present = check_for_node(cell_node, "fill");
   bool material_present = check_for_node(cell_node, "material");
@@ -603,6 +713,25 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
   }
   rpn_.shrink_to_fit();
 
+  if (virtual_lattice_) {
+    vl_triso_distribution_ = generate_triso_distribution(vl_shape_, vl_pitch_, vl_lower_left_, rpn_, id_);
+    write_message("successfully generated virtual lattice", 5);
+    for (int tri_sur_id : vl_triso_distribution_[floor(vl_shape_[0]/2)+floor(vl_shape_[1]/2)*vl_shape_[0]+\
+    floor(vl_shape_[2]/2)*vl_shape_[0]*vl_shape_[1]])
+    {
+      std::cout<<model::surfaces[abs(tri_sur_id)-1]->id_<<std::endl;
+
+    }
+  }
+
+  if (triso_particle_) {
+    if (rpn_.size() != 1) {
+      fatal_error(fmt::format("Wrong surface definition of triso particle cell {}", id_));
+    } else {
+      model::surfaces[abs(rpn_[0]) - 1]->connect_to_triso_base(id_, "particle");
+    }
+  }
+
   // Read the translation vector.
   if (check_for_node(cell_node, "translation")) {
     if (fill_ == C_NONE) {
@@ -644,26 +773,136 @@ std::pair<double, int32_t> CSGCell::distance(
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
+  double min_dis_vl;
+  int32_t i_surf_vl;
+  if (virtual_lattice_) {
+    /*
+    double tol_dis = p.collision_distance();
+    double u_value = sqrt(pow(u.x,2)+pow(u.y,2)+pow(u.z,2)); //don't know if u has been normalized
+    vector<double> norm_u = {u.x/u_value, u.y/u_value, u.z/u_value};
+    vector<double> r_end = {r.x+tol_dis*norm_u[0], r.y+tol_dis*norm_u[1], r.z+tol_dis*norm_u[2]};
+    double temp_pos_x, temp_pos_y, temp_pos_z;
+    vector<vector<int>> passed_lattice={{floor((r.x-vl_lower_left_[0])/vl_pitch_[0]),\
+                                          floor((r.y-vl_lower_left_[1])/vl_pitch_[1]),\
+                                          floor((r.z-vl_lower_left_[2])/vl_pitch_[2]), 0}};
+    int index_start;
+    int index_end;
+    for (int i = 0; i < 3; i++){
+      if (passed_lattice[0][i] == vl_shape_[i]) {
+        passed_lattice[0][i] = vl_shape_[i]-1;
+      }
+    }
 
-  for (int32_t token : rpn_) {
-    // Ignore this token if it corresponds to an operator rather than a region.
-    if (token >= OP_UNION)
-      continue;
+    if (u.x > 0) {
+      index_start = ceil((r.x-vl_lower_left_[0])/vl_pitch_[0]);
+      index_end = floor((r_end[0]-vl_lower_left_[0])/vl_pitch_[0]);
+      if (index_start <= index_end && index_start < vl_shape_[0]) {
+        for (int i = index_start; i <= index_end; i++) {
+          if (i >= vl_shape_[0]) break;
+          temp_pos_x = i*vl_pitch_[0]+vl_lower_left_[0];
+          temp_pos_y = (temp_pos_x-r.x)*norm_u[1]/norm_u[0]+r.y;
+          temp_pos_z = (temp_pos_x-r.x)*norm_u[2]/norm_u[0]+r.z;
+          passed_lattice.push_back({i, floor((temp_pos_y-vl_lower_left_[1])/vl_pitch_[1]),\
+                                    floor((temp_pos_z-vl_lower_left_[2])/vl_pitch_[2]),\
+                                    (temp_pos_x-r.x)/norm_u[0]})
+        }
+      }
+    } 
+    */
 
-    // Calculate the distance to this surface.
-    // Note the off-by-one indexing
-    bool coincident {std::abs(token) == std::abs(on_surface)};
-    double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
+    double max_dis = p->collision_distance();
+    double tol_dis = 0;
+    vector<double> dis_to_bou(3), dis_to_bou_max(3);
+    double u_value = sqrt(pow(u.x,2)+pow(u.y,2)+pow(u.z,2)); //don't know if u has been normalized
+    vector<double> norm_u = {u.x/u_value, u.y/u_value, u.z/u_value};
+    vector<int> lat_ind(3);
+    vector<double> temp_pos ={r.x,r.y,r.z};
+    int loop_time;
+    for (int i=0; i<3; i++) {
+      lat_ind[i]=floor((temp_pos[i]-vl_lower_left_[i])/vl_pitch_[i]);
+    }
 
-    // Check if this distance is the new minimum.
-    if (d < min_dist) {
-      if (min_dist - d >= FP_PRECISION * min_dist) {
-        min_dist = d;
-        i_surf = -token;
+    dis_to_bou = {INFTY,INFTY,INFTY};
+    for (int i=0; i<3; i++) {
+      if (norm_u[i] > 0) {
+        dis_to_bou[i]=std::abs(((lat_ind[i]+1)*vl_pitch_[i]+vl_lower_left_[i]-temp_pos[i])/norm_u[i]);
+        dis_to_bou_max[i]=vl_pitch_[i]/norm_u[i];
+      } else if (norm_u[i] < 0){
+        dis_to_bou[i]=std::abs((lat_ind[i]*vl_pitch_[i]+vl_lower_left_[i]-temp_pos[i])/norm_u[i]);
+        dis_to_bou_max[i]=-vl_pitch_[i]/norm_u[i];
+      }
+    }
+
+    while(true) {
+      
+      
+
+      if (lat_ind[0] < 0 || lat_ind[0] >= vl_shape_[0] ||\
+          lat_ind[1] < 0 || lat_ind[1] >= vl_shape_[1] ||\
+          lat_ind[2] < 0 || lat_ind[2] >= vl_shape_[2]) break;
+
+      
+      for (int token : vl_triso_distribution_[lat_ind[0]+lat_ind[1]*vl_shape_[0]+lat_ind[2]*vl_shape_[0]*vl_shape_[1]]) {
+        bool coincident {std::abs(token) == std::abs(on_surface)};
+        double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
+        if (d < min_dist) {
+          if (min_dist - d >= FP_PRECISION * min_dist) {
+            min_dist = d;
+            i_surf = -token;
+          }
+        }
+      }
+      
+
+      int mes_bou_crossed=0;
+      if  (dis_to_bou[1] < dis_to_bou[0]) {
+        mes_bou_crossed=1;
+      }
+      if  (dis_to_bou[2] < dis_to_bou[mes_bou_crossed]) {
+        mes_bou_crossed=2;
+      }
+
+      tol_dis = dis_to_bou[mes_bou_crossed];
+
+      if (min_dist < tol_dis) {
+        break;
+      }
+
+      if (norm_u[mes_bou_crossed] > 0) {
+        lat_ind[mes_bou_crossed] += 1;
+      } else {
+        lat_ind[mes_bou_crossed] += -1;
+      }
+
+      dis_to_bou[mes_bou_crossed] += dis_to_bou_max[mes_bou_crossed];
+      
+      if (tol_dis > max_dis) {
+        break;
+      }
+    }
+
+    
+  } else {
+    for (int32_t token : rpn_) {
+      // Ignore this token if it corresponds to an operator rather than a region.
+      if (token >= OP_UNION)
+        continue;
+
+      // Calculate the distance to this surface.
+      // Note the off-by-one indexing
+      bool coincident {std::abs(token) == std::abs(on_surface)};
+      double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
+
+      // Check if this distance is the new minimum.
+      if (d < min_dist) {
+        if (min_dist - d >= FP_PRECISION * min_dist) {
+          min_dist = d;
+          i_surf = -token;
+        }
       }
     }
   }
-
+  
   return {min_dist, i_surf};
 }
 
@@ -1065,6 +1304,9 @@ void read_cells(pugi::xml_node node)
     } else {
       model::universes[it->second]->cells_.push_back(i);
     }
+    if (model::cells[i]->virtual_lattice_) {
+      model::universes[it->second]->filled_with_triso_base_=model::cells[i]->id_;
+    }
   }
   model::universes.shrink_to_fit();
 

src/geometry.cpp

@@ -148,7 +148,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
       const auto& univ {model::universes[i_universe]};
       found = univ->find_cell(p);
     }
-
+    
     if (!found) {
       return found;
     }

src/particle.cpp

@@ -184,7 +184,6 @@ void Particle::event_calculate_xs()
 void Particle::event_advance()
 {
   // Find the distance to the nearest boundary
-  boundary() = distance_to_boundary(*this);
 
   // Sample a distance to collision
   if (type() == ParticleType::electron || type() == ParticleType::positron) {
@@ -195,6 +194,8 @@ void Particle::event_advance()
     collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
   }
 
+  boundary() = distance_to_boundary(*this);
+
   // Select smaller of the two distances
   double distance = std::min(boundary().distance, collision_distance());
 
@@ -396,7 +397,10 @@ void Particle::cross_surface()
   if (settings::verbosity >= 10 || trace()) {
     write_message(1, "    Crossing surface {}", surf->id_);
   }
-
+  /*
+  if (surf->id_ >= 11675 && surf->id_ <=23340) {
+    write_message(1, "    {}", surf->id_);
+  }*/
   if (surf->surf_source_ && simulation::current_batch == settings::n_batches) {
     SourceSite site;
     site.r = r();
@@ -449,8 +453,96 @@ void Particle::cross_surface()
   }
 #endif
 
-  if (neighbor_list_find_cell(*this))
-    return;
+  
+  if (surf->is_triso_surface_) {
+    if (surface() > 0){
+      for (int i = n_coord(); i < model::n_coord_levels; i++) {
+        coord(i).reset();
+      }
+      coord(n_coord() - 1).cell = model::cell_map[model::surfaces[i_surface - 1]->triso_base_index_];
+    } else if (surface() < 0) {
+      for (int i = n_coord(); i < model::n_coord_levels; i++) {
+        coord(i).reset();
+      }
+      if (model::surfaces[i_surface - 1]->triso_particle_index_==-1) {
+        fatal_error(fmt::format("Particle cell of surface {} is not defined", model::surfaces[i_surface - 1]->id_));
+      }
+      coord(n_coord() - 1).cell = model::cell_map[model::surfaces[i_surface - 1]->triso_particle_index_];
+    }
+
+    //find material
+    bool found=true;
+    int i_cell = coord(n_coord() - 1).cell;
+    for (;; ++n_coord()) {
+      if (i_cell == C_NONE) {
+        int i_universe = coord(n_coord() - 1).universe;
+        const auto& univ {model::universes[i_universe]};
+        
+        if (univ->filled_with_triso_base_ != -1) {
+          coord(n_coord() - 1).cell = model::cell_map[univ->filled_with_triso_base_];
+          found=true;
+        } else {
+          found = univ->find_cell(*this);
+        }
+        if (!found) {
+          break;
+        }
+      }
+
+      i_cell = coord(n_coord() - 1).cell;
+
+      Cell& c {*model::cells[i_cell]};
+      if (c.type_ == Fill::MATERIAL) {
+        // Found a material cell which means this is the lowest coord level.
+
+        cell_instance() = 0;
+        // Find the distribcell instance number.
+        if (c.distribcell_index_ >= 0) {
+          cell_instance() = cell_instance_at_level(*this, n_coord() - 1);
+        }
+
+        // Set the material and temperature.
+        material_last() = material();
+        if (c.material_.size() > 1) {
+          material() = c.material_[cell_instance()];
+        } else {
+          material() = c.material_[0];
+        }
+        sqrtkT_last() = sqrtkT();
+        if (c.sqrtkT_.size() > 1) {
+          sqrtkT() = c.sqrtkT_[cell_instance()];
+        } else {
+          sqrtkT() = c.sqrtkT_[0];
+        }
+        return;
+
+      } else if (c.type_ == Fill::UNIVERSE) {
+        //========================================================================
+        //! Found a lower universe, update this coord level then search the next.
+
+        // Set the lower coordinate level universe.
+        auto& coor {coord(n_coord())};
+        coor.universe = c.fill_;
+
+        // Set the position and direction.
+        coor.r = r_local();
+        coor.u = u_local();
+
+        // Apply translation.
+        coor.r -= c.translation_;
+
+        // Apply rotation.
+        if (!c.rotation_.empty()) {
+          coor.rotate(c.rotation_);
+        }
+        i_cell = C_NONE;
+      }
+    }
+  } else {
+      if (neighbor_list_find_cell(*this))
+      return;
+  }
+  
 
   // ==========================================================================
   // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS

src/surface.cpp

@@ -313,10 +313,20 @@ BoundingBox SurfaceXPlane::bounding_box(bool pos_side) const
   }
 }
 
+bool SurfaceXPlane::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 //==============================================================================
 // SurfaceYPlane implementation
 //==============================================================================
 
+bool SurfaceYPlane::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, y0_);
@@ -357,6 +367,11 @@ BoundingBox SurfaceYPlane::bounding_box(bool pos_side) const
 // SurfaceZPlane implementation
 //==============================================================================
 
+bool SurfaceZPlane::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, z0_);
@@ -396,6 +411,10 @@ BoundingBox SurfaceZPlane::bounding_box(bool pos_side) const
 //==============================================================================
 // SurfacePlane implementation
 //==============================================================================
+bool SurfacePlane::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
 
 SurfacePlane::SurfacePlane(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
@@ -515,6 +534,11 @@ Direction axis_aligned_cylinder_normal(
 // SurfaceXCylinder implementation
 //==============================================================================
 
+bool SurfaceXCylinder::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node)
   : CSGSurface(surf_node)
 {
@@ -558,6 +582,11 @@ BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const
 // SurfaceYCylinder implementation
 //==============================================================================
 
+bool SurfaceYCylinder::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node)
   : CSGSurface(surf_node)
 {
@@ -602,6 +631,11 @@ BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const
 // SurfaceZCylinder implementation
 //==============================================================================
 
+bool SurfaceZCylinder::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
   : CSGSurface(surf_node)
 {
@@ -720,6 +754,66 @@ BoundingBox SurfaceSphere::bounding_box(bool pos_side) const
   }
 }
 
+void SurfaceSphere::connect_to_triso_base(int triso_index, std::string key)
+{
+  if (key=="base") {
+    triso_base_index_=triso_index;
+    is_triso_surface_=true;
+  } else if (key=="particle") {
+    triso_particle_index_=triso_index;
+  }
+}
+
+vector<double> SurfaceSphere::get_center() const {
+  return {x0_,y0_,z0_};
+}
+
+double SurfaceSphere::get_radius() const {
+  return radius_;
+}
+
+bool SurfaceSphere::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  double dis_x;
+  double dis_y;
+  double dis_z;
+  double x_min=mesh_center[0]-lattice_pitch[0]/2;
+  double x_max=mesh_center[0]+lattice_pitch[0]/2;
+  double y_min=mesh_center[1]-lattice_pitch[1]/2;
+  double y_max=mesh_center[1]+lattice_pitch[1]/2;
+  double z_min=mesh_center[2]-lattice_pitch[2]/2;
+  double z_max=mesh_center[2]+lattice_pitch[2]/2;
+  if (x0_>=x_min && x0_<=x_max) {
+    dis_x=0;
+  } else if (x0_<x_min) {
+    dis_x=pow(x_min-x0_, 2);
+  } else {
+    dis_x=pow(x_max-x0_, 2);
+  }
+
+  if (y0_>=y_min && y0_<=y_max) {
+    dis_y=0;
+  } else if (y0_<y_min) {
+    dis_y=pow(y_min-y0_, 2);
+  } else {
+    dis_y=pow(y_max-y0_, 2);
+  }
+
+  if (z0_>=z_min && z0_<=z_max) {
+    dis_z=0;
+  } else if (z0_<z_min) {
+    dis_z=pow(z_min-z0_, 2);
+  } else {
+    dis_z=pow(z_max-z0_, 2);
+  }
+
+  if (sqrt(dis_x+dis_y+dis_z) < radius_) {
+    return true;
+  } else {
+    return false;
+  }
+}
+
 //==============================================================================
 // Generic functions for x-, y-, and z-, cones
 //==============================================================================
@@ -812,6 +906,11 @@ Direction axis_aligned_cone_normal(
 // SurfaceXCone implementation
 //==============================================================================
 
+bool SurfaceXCone::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
@@ -844,6 +943,11 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
 // SurfaceYCone implementation
 //==============================================================================
 
+bool SurfaceYCone::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
@@ -876,6 +980,11 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
 // SurfaceZCone implementation
 //==============================================================================
 
+bool SurfaceZCone::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
@@ -908,6 +1017,11 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
 // SurfaceQuadric implementation
 //==============================================================================
 
+bool SurfaceQuadric::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, A_, B_, C_, D_, E_, F_, G_, H_, J_, K_);
@@ -1058,6 +1172,11 @@ double torus_distance(double x1, double x2, double x3, double u1, double u2,
 // SurfaceXTorus implementation
 //==============================================================================
 
+bool SurfaceXTorus::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceXTorus::SurfaceXTorus(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, A_, B_, C_);
@@ -1111,6 +1230,11 @@ Direction SurfaceXTorus::normal(Position r) const
 // SurfaceYTorus implementation
 //==============================================================================
 
+bool SurfaceYTorus::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
+
 SurfaceYTorus::SurfaceYTorus(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, A_, B_, C_);
@@ -1163,6 +1287,10 @@ Direction SurfaceYTorus::normal(Position r) const
 //==============================================================================
 // SurfaceZTorus implementation
 //==============================================================================
+bool SurfaceZTorus::triso_in_mesh(vector<double> mesh_center, vector<double> lattice_pitch) const
+{
+  return false;
+}
 
 SurfaceZTorus::SurfaceZTorus(pugi::xml_node surf_node) : CSGSurface(surf_node)
 {