| 36 | namespace |
| 37 | { |
| 38 | std::vector<double> |
| 39 | compute_heaviside_integrals(Pointer<HierarchyMathOps> hier_math_ops, int phi_idx, double ncells) |
| 40 | { |
| 41 | const int wgt_cc_idx = hier_math_ops->getCellWeightPatchDescriptorIndex(); |
| 42 | Pointer<PatchHierarchy<NDIM>> patch_hier = hier_math_ops->getPatchHierarchy(); |
| 43 | |
| 44 | const int hier_finest_ln = patch_hier->getFinestLevelNumber(); |
| 45 | double vol_phase1 = 0.0; |
| 46 | double vol_phase2 = 0.0; |
| 47 | double integral_delta = 0.0; |
| 48 | for (int ln = 0; ln <= hier_finest_ln; ++ln) |
| 49 | { |
| 50 | Pointer<PatchLevel<NDIM>> patch_level = patch_hier->getPatchLevel(ln); |
| 51 | for (PatchLevel<NDIM>::Iterator p(patch_level); p; p++) |
| 52 | { |
| 53 | Pointer<Patch<NDIM>> patch = patch_level->getPatch(p()); |
| 54 | const Box<NDIM>& patch_box = patch->getBox(); |
| 55 | |
| 56 | Pointer<CellData<NDIM, double>> phi_data = patch->getPatchData(phi_idx); |
| 57 | Pointer<CellData<NDIM, double>> wgt_data = patch->getPatchData(wgt_cc_idx); |
| 58 | |
| 59 | // Get grid spacing information |
| 60 | Pointer<CartesianPatchGeometry<NDIM>> patch_geom = patch->getPatchGeometry(); |
| 61 | const double* const patch_dx = patch_geom->getDx(); |
| 62 | double cell_size = 1.0; |
| 63 | for (int d = 0; d < NDIM; ++d) cell_size *= patch_dx[d]; |
| 64 | cell_size = std::pow(cell_size, 1.0 / static_cast<double>(NDIM)); |
| 65 | const double alpha = ncells * cell_size; |
| 66 | |
| 67 | for (Box<NDIM>::Iterator it(patch_box); it; it++) |
| 68 | { |
| 69 | CellIndex<NDIM> ci(it()); |
| 70 | |
| 71 | const double phi = (*phi_data)(ci); |
| 72 | const double dv = (*wgt_data)(ci); |
| 73 | |
| 74 | // smoothed delta and Heaviside functions |
| 75 | const double h_phi = IBTK::smooth_heaviside(phi, alpha); |
| 76 | const double h_prime = IBTK::smooth_delta(phi, alpha); |
| 77 | |
| 78 | vol_phase1 += (1.0 - h_phi) * dv; |
| 79 | vol_phase2 += h_phi * dv; |
| 80 | integral_delta += h_prime * dv; |
| 81 | } |
| 82 | } |
| 83 | } |
| 84 | |
| 85 | std::vector<double> integrals{ vol_phase1, vol_phase2, integral_delta }; |
| 86 | IBTK_MPI::sumReduction(&integrals[0], integrals.size()); |
| 87 | |
| 88 | return integrals; |
| 89 | } // compute_heaviside_integrals |
| 90 | |
| 91 | std::vector<double> |
| 92 | compute_heaviside_integrals(Pointer<HierarchyMathOps> hier_math_ops, int phi_idx, int psi_idx, double ncells) |
no test coverage detected