| 25 | } |
| 26 | |
| 27 | void EBToPVD::EBToPolygon(const Real* problo, const Real* dx, |
| 28 | const Box & bx, Array4<EBCellFlag const> const& flag, |
| 29 | Array4<Real const> const& bcent, |
| 30 | Array4<Real const> const& apx, Array4<Real const> const& apy, Array4<Real const> const& apz) |
| 31 | { |
| 32 | const auto lo = lbound(bx); |
| 33 | const auto hi = ubound(bx); |
| 34 | |
| 35 | for(int k = lo.z; k <= hi.z; ++k) { |
| 36 | for(int j = lo.y; j <= hi.y; ++j) { |
| 37 | for(int i = lo.x; i <= hi.x; ++i) { |
| 38 | // NOTE: do not skip fully enclosed cells (is_covered_cell), as this seems |
| 39 | // to skip thin walls in the domain: |
| 40 | // if(.not.is_regular_cell(flag(i,j,k)) .and. & |
| 41 | // .not.is_covered_cell(flag(i,j,k))) then |
| 42 | |
| 43 | // Instead only look for EBs |
| 44 | // if( .not.is_regular_cell(flag(i,j,k))) then |
| 45 | |
| 46 | // If covered cells are accounted for in this loop, a FPE arises |
| 47 | // since apnorm is zero. |
| 48 | |
| 49 | if(flag(i,j,k).isSingleValued()) { |
| 50 | Real axm = apx(i ,j ,k ); |
| 51 | Real axp = apx(i+1,j ,k ); |
| 52 | Real aym = apy(i ,j ,k ); |
| 53 | Real ayp = apy(i ,j+1,k ); |
| 54 | Real azm = apz(i ,j ,k ); |
| 55 | Real azp = apz(i ,j ,k+1); |
| 56 | |
| 57 | Real adx = (axm-axp) * dx[1] * dx[2]; |
| 58 | Real ady = (aym-ayp) * dx[0] * dx[2]; |
| 59 | Real adz = (azm-azp) * dx[0] * dx[1]; |
| 60 | |
| 61 | Real apnorm = std::sqrt(adx*adx + ady*ady + adz*adz); |
| 62 | Real apnorminv = Real(1.0)/apnorm; |
| 63 | |
| 64 | std::array<Real,3> normal, centroid; |
| 65 | std::array<std::array<Real,3>,8> vertex; |
| 66 | |
| 67 | normal[0] = adx * apnorminv; |
| 68 | normal[1] = ady * apnorminv; |
| 69 | normal[2] = adz * apnorminv; |
| 70 | |
| 71 | // convert bcent to global coordinate system centered at plo |
| 72 | centroid[0] = problo[0] + bcent(i,j,k,0)*dx[0] + (static_cast<Real>(i) + Real(0.5))*dx[0]; |
| 73 | centroid[1] = problo[1] + bcent(i,j,k,1)*dx[1] + (static_cast<Real>(j) + Real(0.5))*dx[1]; |
| 74 | centroid[2] = problo[2] + bcent(i,j,k,2)*dx[2] + (static_cast<Real>(k) + Real(0.5))*dx[2]; |
| 75 | |
| 76 | // vertices of bounding cell (i,j,k) |
| 77 | vertex[0] = {problo[0] + static_cast<Real>(i )*dx[0], problo[1] + static_cast<Real>(j )*dx[1], problo[2] + static_cast<Real>(k )*dx[2]}; |
| 78 | vertex[1] = {problo[0] + static_cast<Real>(i+1)*dx[0], problo[1] + static_cast<Real>(j )*dx[1], problo[2] + static_cast<Real>(k )*dx[2]}; |
| 79 | vertex[2] = {problo[0] + static_cast<Real>(i )*dx[0], problo[1] + static_cast<Real>(j+1)*dx[1], problo[2] + static_cast<Real>(k )*dx[2]}; |
| 80 | vertex[3] = {problo[0] + static_cast<Real>(i+1)*dx[0], problo[1] + static_cast<Real>(j+1)*dx[1], problo[2] + static_cast<Real>(k )*dx[2]}; |
| 81 | vertex[4] = {problo[0] + static_cast<Real>(i )*dx[0], problo[1] + static_cast<Real>(j )*dx[1], problo[2] + static_cast<Real>(k+1)*dx[2]}; |
| 82 | vertex[5] = {problo[0] + static_cast<Real>(i+1)*dx[0], problo[1] + static_cast<Real>(j )*dx[1], problo[2] + static_cast<Real>(k+1)*dx[2]}; |
| 83 | vertex[6] = {problo[0] + static_cast<Real>(i )*dx[0], problo[1] + static_cast<Real>(j+1)*dx[1], problo[2] + static_cast<Real>(k+1)*dx[2]}; |
| 84 | vertex[7] = {problo[0] + static_cast<Real>(i+1)*dx[0], problo[1] + static_cast<Real>(j+1)*dx[1], problo[2] + static_cast<Real>(k+1)*dx[2]}; |
no test coverage detected