From a3e008db095ade9f52a78c62aec157bac727b82a Mon Sep 17 00:00:00 2001 From: shriwise Date: Thu, 17 May 2018 18:42:23 -0500 Subject: [PATCH] Some hacking to get a simple geometry running. The simulations agree though :tada: --- src/cad.cpp | 8 ++++---- src/cell.cpp | 37 +++++++++++++++++++++++++++++++++---- src/geometry.F90 | 30 +++++++++++++++++++++++++++++- src/tracking.F90 | 25 ++++++++++++++++++++++++- 4 files changed, 90 insertions(+), 10 deletions(-) diff --git a/src/cad.cpp b/src/cad.cpp index be008174e..2271aa856 100644 --- a/src/cad.cpp +++ b/src/cad.cpp @@ -23,12 +23,12 @@ void load_cad_geometry_c() { // set cell ids using global IDs openmc::CADCell* c = new openmc::CADCell(); - c->id_ = DAGMC->id_by_index(3, i); + c->id_ = DAGMC->id_by_index(3, i+1); c->dagmc_ptr = DAGMC; c->universe_ = cad_univ_id; // set to zero for now c->material_.push_back(40); // TEMPORARY openmc::cells_c.push_back(c); - openmc::cell_dict[c->id_] = i; + openmc::cell_dict[c->id_] = c->id; // Populate the Universe vector and dict auto it = openmc::universe_dict.find(cad_univ_id); @@ -51,11 +51,11 @@ void load_cad_geometry_c() { // set cell ids using global IDs openmc::CADSurface* s = new openmc::CADSurface(); - s->id = DAGMC->id_by_index(2, i); + s->id = DAGMC->id_by_index(2, i+1); s->dagmc_ptr = DAGMC; s->bc = openmc::BC_TRANSMIT; openmc::surfaces_c[i] = s; - openmc::surface_dict[s->id] = i; + openmc::surface_dict[s->id] = s->id; } return; diff --git a/src/cell.cpp b/src/cell.cpp index 6060edc7a..c0a28cedb 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -598,14 +598,26 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const CADCell::CADCell() : Cell{} {}; std::pair CADCell::distance(const double xyz[3], const double uvw[3], int32_t on_surface) const { - std::cout << "Cell Distance" << std::endl; - std::pair result(1.0, 0); + + moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id); + moab::EntityHandle hit_surf; + double dist; + dagmc_ptr->ray_fire(vol, xyz, uvw, hit_surf, dist); + + int surf_idx = dagmc_ptr->index_by_handle(hit_surf); + + std::pair result(dist, surf_idx); + return result; } bool CADCell::contains(const double xyz[3], const double uvw[3], int32_t on_surface) const { - std::cout << "Cell Contains" << std::endl; - return true; + moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id); + + int result = 0; + dagmc_ptr->point_in_volume(vol, xyz, result, uvw); + + return bool(result); } void CADCell::to_hdf5(hid_t group_id) const { return; } @@ -749,6 +761,23 @@ extern "C" { int cell_type(Cell* c) {return c->type_;} +#ifdef CAD + int32_t next_cell(CADCell* cur_cell, CADSurface *surf_xed ) { + moab::EntityHandle surf = surf_xed->dagmc_ptr->entity_by_id(2,surf_xed->id); + moab::EntityHandle vol = cur_cell->dagmc_ptr->entity_by_id(3,cur_cell->id); + + moab::EntityHandle new_vol; + cur_cell->dagmc_ptr->next_vol(surf, vol, new_vol); + + return cur_cell->dagmc_ptr->index_by_handle(new_vol); + } + + bool is_implicit_complement(CADCell *c) { + moab::EntityHandle handle = c->dagmc_ptr->entity_by_id(3,c->id); + return c->dagmc_ptr->is_implicit_complement(handle); + } +#endif + int32_t cell_universe(Cell* c) {return c->universe_;} int32_t cell_fill(Cell* c) {return c->fill_;} diff --git a/src/geometry.F90 b/src/geometry.F90 index b5dee81fa..b6212d73a 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -55,10 +55,38 @@ module geometry subroutine neighbor_lists() bind(C) end subroutine neighbor_lists - end interface + + function next_cell_c(current_cell, surface_crossed) & + bind(C, name="next_cell") result(new_cell) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: current_cell + type(C_PTR), intent(in), value :: surface_crossed + integer(C_INT32_T) :: new_cell + end function next_cell_c + + function is_implicit_complement_C(cell) & + bind(C, name="is_implicit_complement") result(res) + import C_PTR, C_BOOL + type(C_PTR), intent(in), value :: cell + logical(C_BOOL) :: res + end function is_implicit_complement_C + end interface contains + function next_cell(c, s) result(new_cell) + type(Cell), intent(in) :: c + type(Surface), intent(in) :: s + integer :: new_cell + new_cell = next_cell_c(c%ptr, s%ptr) + end function next_cell + + function is_implicit_complement(c) result(res) + type(Cell), intent(in) :: c + logical:: res + res = is_implicit_complement_c(c%ptr) + end function is_implicit_complement + !=============================================================================== ! FIND_CELL determines what cell a source particle is in within a particular ! universe. If the base universe is passed, the particle should be found as long diff --git a/src/tracking.F90 b/src/tracking.F90 index d93ac9ab9..67c76f183 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -6,7 +6,7 @@ module tracking use error, only: warning, write_message use geometry_header, only: cells use geometry, only: find_cell, distance_to_boundary, cross_lattice,& - check_cell_overlap + check_cell_overlap, next_cell, is_implicit_complement use material_header, only: materials, Material use message_passing use mgxs_interface @@ -309,11 +309,34 @@ contains real(8) :: norm ! "norm" of surface normal real(8) :: xyz(3) ! Saved global coordinate integer :: i_surface ! index in surfaces + integer :: i_cell ! index of new cell logical :: rotational ! if rotational periodic BC applied logical :: found ! particle found in universe? class(Surface), pointer :: surf class(Surface), pointer :: surf2 ! periodic partner surface + +#ifdef CAD + if (dagmc) then + i_cell = next_cell(cells(p % last_cell(1)), surfaces(ABS(p % surface))) + ! save material and temp + p % last_material = p % material + p % last_sqrtkT = p % sqrtKT + ! set new cell value + p % coord(1) % cell = i_cell + p % cell_instance = 1 + p % material = cells(i_cell) % material(1) + p % sqrtKT = cells(i_cell) % sqrtKT(1) + if (is_implicit_complement(cells(i_cell))) then + p % alive = .false. + end if + + return + + end if +#endif + + i_surface = abs(p % surface) surf => surfaces(i_surface) if (verbosity >= 10 .or. trace) then