| 843 | } |
| 844 | |
| 845 | void |
| 846 | CFINSForcing::projectTensor(const int data_idx, |
| 847 | const Pointer<Variable<NDIM>> /*var*/, |
| 848 | const double /*data_time*/, |
| 849 | const bool initial_time, |
| 850 | const bool extended_box) |
| 851 | { |
| 852 | Pointer<PatchHierarchy<NDIM>> hierarchy = d_adv_diff_integrator->getPatchHierarchy(); |
| 853 | for (int ln = 0; ln <= hierarchy->getFinestLevelNumber(); ln++) |
| 854 | { |
| 855 | Pointer<PatchLevel<NDIM>> level = hierarchy->getPatchLevel(ln); |
| 856 | for (PatchLevel<NDIM>::Iterator p(level); p; p++) |
| 857 | { |
| 858 | Pointer<Patch<NDIM>> patch = level->getPatch(p()); |
| 859 | if (initial_time) return; |
| 860 | Pointer<CellData<NDIM, double>> data = patch->getPatchData(data_idx); |
| 861 | const Box<NDIM>& box = extended_box ? data->getGhostBox() : patch->getBox(); |
| 862 | for (CellIterator<NDIM> it(box); it; it++) |
| 863 | { |
| 864 | CellIndex<NDIM> i = *it; |
| 865 | MatrixNd tens; |
| 866 | Eigen::SelfAdjointEigenSolver<MatrixNd> eigs; |
| 867 | for (int k = 0; k < NDIM * (NDIM + 1) / 2; ++k) |
| 868 | { |
| 869 | const std::pair<int, int>& idx = voigt_to_tensor_idx(k); |
| 870 | tens(idx.first, idx.second) = tens(idx.second, idx.first) = (*data)(i, k); |
| 871 | } |
| 872 | eigs.computeDirect(tens); |
| 873 | MatrixNd eig_vals(MatrixNd::Zero()); |
| 874 | for (int d = 0; d < NDIM; ++d) |
| 875 | { |
| 876 | eig_vals(d, d) = std::max(eigs.eigenvalues()(d), 0.0); |
| 877 | } |
| 878 | MatrixNd eig_vecs = eigs.eigenvectors(); |
| 879 | tens = eig_vecs * eig_vals * eig_vecs.transpose(); |
| 880 | for (int k = 0; k < NDIM * (NDIM + 1) / 2; ++k) |
| 881 | { |
| 882 | const std::pair<int, int>& idx = voigt_to_tensor_idx(k); |
| 883 | (*data)(i, k) = tens(idx.first, idx.second); |
| 884 | } |
| 885 | } |
| 886 | } |
| 887 | } |
| 888 | return; |
| 889 | } // projectTensor |
| 890 | |
| 891 | void |
| 892 | CFINSForcing::applyGradientDetector(Pointer<BasePatchHierarchy<NDIM>> hierarchy, |
no test coverage detected