Page MenuHome
Differential D5594 Diff 18116 extern/quadriflow/src/hierarchy.cpp

Changeset View

Standalone View

View Options

extern/quadriflow/src/hierarchy.cpp

  • This file was added.
#include "hierarchy.hpp"
#include <fstream>
#include <algorithm>
#include <unordered_map>
#include "config.hpp"
#include "field-math.hpp"
#include <queue>
#include "localsat.hpp"
#include "pcg32/pcg32.h"
#ifdef WITH_TBB
# include "tbb/tbb.h"
# include "pss/parallel_stable_sort.h"
#endif
namespace qflow {
Hierarchy::Hierarchy() {
mAdj.resize(MAX_DEPTH + 1);
mV.resize(MAX_DEPTH + 1);
mN.resize(MAX_DEPTH + 1);
mA.resize(MAX_DEPTH + 1);
mPhases.resize(MAX_DEPTH + 1);
mToLower.resize(MAX_DEPTH);
mToUpper.resize(MAX_DEPTH);
rng_seed = 0;
mCQ.reserve(MAX_DEPTH + 1);
mCQw.reserve(MAX_DEPTH + 1);
mCO.reserve(MAX_DEPTH + 1);
mCOw.reserve(MAX_DEPTH + 1);
}
#undef max
void Hierarchy::Initialize(double scale, int with_scale) {
this->with_scale = with_scale;
generate_graph_coloring_deterministic(mAdj[0], mV[0].cols(), mPhases[0]);
for (int i = 0; i < MAX_DEPTH; ++i) {
DownsampleGraph(mAdj[i], mV[i], mN[i], mA[i], mV[i + 1], mN[i + 1], mA[i + 1], mToUpper[i],
mToLower[i], mAdj[i + 1]);
generate_graph_coloring_deterministic(mAdj[i + 1], mV[i + 1].cols(), mPhases[i + 1]);
if (mV[i + 1].cols() == 1) {
mAdj.resize(i + 2);
mV.resize(i + 2);
mN.resize(i + 2);
mA.resize(i + 2);
mToUpper.resize(i + 1);
mToLower.resize(i + 1);
break;
}
}
mQ.resize(mV.size());
mO.resize(mV.size());
mS.resize(mV.size());
mK.resize(mV.size());
mCO.resize(mV.size());
mCOw.resize(mV.size());
mCQ.resize(mV.size());
mCQw.resize(mV.size());
//Set random seed
srand(rng_seed);
mScale = scale;
for (int i = 0; i < mV.size(); ++i) {
mQ[i].resize(mN[i].rows(), mN[i].cols());
mO[i].resize(mN[i].rows(), mN[i].cols());
mS[i].resize(2, mN[i].cols());
mK[i].resize(2, mN[i].cols());
for (int j = 0; j < mN[i].cols(); ++j) {
Vector3d s, t;
coordinate_system(mN[i].col(j), s, t);
//rand() is not thread safe!
double angle = ((double)rand()) / RAND_MAX * 2 * M_PI;
double x = ((double)rand()) / RAND_MAX * 2 - 1.f;
double y = ((double)rand()) / RAND_MAX * 2 - 1.f;
mQ[i].col(j) = s * std::cos(angle) + t * std::sin(angle);
mO[i].col(j) = mV[i].col(j) + (s * x + t * y) * scale;
if (with_scale) {
mS[i].col(j) = Vector2d(1.0f, 1.0f);
mK[i].col(j) = Vector2d(0.0, 0.0);
}
}
}
#ifdef WITH_CUDA
printf("copy to device...\n");
CopyToDevice();
printf("copy to device finish...\n");
#endif
}
#ifdef WITH_TBB
void Hierarchy::generate_graph_coloring_deterministic(const AdjacentMatrix& adj, int size,
std::vector<std::vector<int>>& phases) {
struct ColorData {
uint8_t nColors;
uint32_t nNodes[256];
ColorData() : nColors(0) {}
};
const uint8_t INVALID_COLOR = 0xFF;
phases.clear();
/* Generate a permutation */
std::vector<uint32_t> perm(size);
std::vector<tbb::spin_mutex> mutex(size);
for (uint32_t i = 0; i < size; ++i) perm[i] = i;
tbb::parallel_for(tbb::blocked_range<uint32_t>(0u, size, GRAIN_SIZE),
[&](const tbb::blocked_range<uint32_t>& range) {
pcg32 rng;
rng.advance(range.begin());
for (uint32_t i = range.begin(); i != range.end(); ++i) {
uint32_t j = i, k = rng.nextUInt(size - i) + i;
if (j == k) continue;
if (j > k) std::swap(j, k);
tbb::spin_mutex::scoped_lock l0(mutex[j]);
tbb::spin_mutex::scoped_lock l1(mutex[k]);
std::swap(perm[j], perm[k]);
}
});
std::vector<uint8_t> color(size, INVALID_COLOR);
ColorData colorData = tbb::parallel_reduce(
tbb::blocked_range<uint32_t>(0u, size, GRAIN_SIZE), ColorData(),
[&](const tbb::blocked_range<uint32_t>& range, ColorData colorData) -> ColorData {
std::vector<uint32_t> neighborhood;
bool possible_colors[256];
for (uint32_t pidx = range.begin(); pidx != range.end(); ++pidx) {
uint32_t i = perm[pidx];
neighborhood.clear();
neighborhood.push_back(i);
// for (const Link *link = adj[i]; link != adj[i + 1]; ++link)
for (auto& link : adj[i]) neighborhood.push_back(link.id);
std::sort(neighborhood.begin(), neighborhood.end());
for (uint32_t j : neighborhood) mutex[j].lock();
std::fill(possible_colors, possible_colors + colorData.nColors, true);
// for (const Link *link = adj[i]; link != adj[i + 1]; ++link) {
for (auto& link : adj[i]) {
uint8_t c = color[link.id];
if (c != INVALID_COLOR) {
while (c >= colorData.nColors) {
possible_colors[colorData.nColors] = true;
colorData.nNodes[colorData.nColors] = 0;
colorData.nColors++;
}
possible_colors[c] = false;
}
}
uint8_t chosen_color = INVALID_COLOR;
for (uint8_t j = 0; j < colorData.nColors; ++j) {
if (possible_colors[j]) {
chosen_color = j;
break;
}
}
if (chosen_color == INVALID_COLOR) {
if (colorData.nColors == INVALID_COLOR - 1)
throw std::runtime_error(
"Ran out of colors during graph coloring! "
"The input mesh is very likely corrupt.");
colorData.nNodes[colorData.nColors] = 1;
color[i] = colorData.nColors++;
} else {
colorData.nNodes[chosen_color]++;
color[i] = chosen_color;
}
for (uint32_t j : neighborhood) mutex[j].unlock();
}
return colorData;
},
[](ColorData c1, ColorData c2) -> ColorData {
ColorData result;
result.nColors = std::max(c1.nColors, c2.nColors);
memset(result.nNodes, 0, sizeof(uint32_t) * result.nColors);
for (uint8_t i = 0; i < c1.nColors; ++i) result.nNodes[i] += c1.nNodes[i];
for (uint8_t i = 0; i < c2.nColors; ++i) result.nNodes[i] += c2.nNodes[i];
return result;
});
phases.resize(colorData.nColors);
for (int i = 0; i < colorData.nColors; ++i) phases[i].reserve(colorData.nNodes[i]);
for (uint32_t i = 0; i < size; ++i) phases[color[i]].push_back(i);
}
#else
void Hierarchy::generate_graph_coloring_deterministic(const AdjacentMatrix& adj, int size,
std::vector<std::vector<int>>& phases) {
phases.clear();
std::vector<uint32_t> perm(size);
for (uint32_t i = 0; i < size; ++i) perm[i] = i;
pcg32 rng;
rng.shuffle(perm.begin(), perm.end());
std::vector<int> color(size, -1);
std::vector<uint8_t> possible_colors;
std::vector<int> size_per_color;
int ncolors = 0;
for (uint32_t i = 0; i < size; ++i) {
uint32_t ip = perm[i];
std::fill(possible_colors.begin(), possible_colors.end(), 1);
for (auto& link : adj[ip]) {
int c = color[link.id];
if (c >= 0) possible_colors[c] = 0;
}
int chosen_color = -1;
for (uint32_t j = 0; j < possible_colors.size(); ++j) {
if (possible_colors[j]) {
chosen_color = j;
break;
}
}
if (chosen_color < 0) {
chosen_color = ncolors++;
possible_colors.resize(ncolors);
size_per_color.push_back(0);
}
color[ip] = chosen_color;
size_per_color[chosen_color]++;
}
phases.resize(ncolors);
for (int i = 0; i < ncolors; ++i) phases[i].reserve(size_per_color[i]);
for (uint32_t i = 0; i < size; ++i) phases[color[i]].push_back(i);
}
#endif
void Hierarchy::DownsampleGraph(const AdjacentMatrix adj, const MatrixXd& V, const MatrixXd& N,
const VectorXd& A, MatrixXd& V_p, MatrixXd& N_p, VectorXd& A_p,
MatrixXi& to_upper, VectorXi& to_lower, AdjacentMatrix& adj_p) {
struct Entry {
int i, j;
double order;
inline Entry() { i = j = -1; };
inline Entry(int i, int j, double order) : i(i), j(j), order(order) {}
inline bool operator<(const Entry& e) const { return order > e.order; }
inline bool operator==(const Entry& e) const { return order == e.order; }
};
int nLinks = 0;
for (auto& adj_i : adj) nLinks += adj_i.size();
std::vector<Entry> entries(nLinks);
std::vector<int> bases(adj.size());
for (int i = 1; i < bases.size(); ++i) {
bases[i] = bases[i - 1] + adj[i - 1].size();
}
#ifdef WITH_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < V.cols(); ++i) {
int base = bases[i];
auto& ad = adj[i];
auto entry_it = entries.begin() + base;
for (auto it = ad.begin(); it != ad.end(); ++it, ++entry_it) {
int k = it->id;
double dp = N.col(i).dot(N.col(k));
double ratio = A[i] > A[k] ? (A[i] / A[k]) : (A[k] / A[i]);
*entry_it = Entry(i, k, dp * ratio);
}
}
#ifdef WITH_TBB
pss::parallel_stable_sort(entries.begin(), entries.end(), std::less<Entry>());
#else
std::stable_sort(entries.begin(), entries.end(), std::less<Entry>());
#endif
std::vector<bool> mergeFlag(V.cols(), false);
int nCollapsed = 0;
for (int i = 0; i < nLinks; ++i) {
const Entry& e = entries[i];
if (mergeFlag[e.i] || mergeFlag[e.j]) continue;
mergeFlag[e.i] = mergeFlag[e.j] = true;
entries[nCollapsed++] = entries[i];
}
int vertexCount = V.cols() - nCollapsed;
// Allocate memory for coarsened graph
V_p.resize(3, vertexCount);
N_p.resize(3, vertexCount);
A_p.resize(vertexCount);
to_upper.resize(2, vertexCount);
to_lower.resize(V.cols());
#ifdef WITH_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < nCollapsed; ++i) {
const Entry& e = entries[i];
const double area1 = A[e.i], area2 = A[e.j], surfaceArea = area1 + area2;
if (surfaceArea > RCPOVERFLOW)
V_p.col(i) = (V.col(e.i) * area1 + V.col(e.j) * area2) / surfaceArea;
else
V_p.col(i) = (V.col(e.i) + V.col(e.j)) * 0.5f;
Vector3d normal = N.col(e.i) * area1 + N.col(e.j) * area2;
double norm = normal.norm();
N_p.col(i) = norm > RCPOVERFLOW ? Vector3d(normal / norm) : Vector3d::UnitX();
A_p[i] = surfaceArea;
to_upper.col(i) << e.i, e.j;
to_lower[e.i] = i;
to_lower[e.j] = i;
}
int offset = nCollapsed;
for (int i = 0; i < V.cols(); ++i) {
if (!mergeFlag[i]) {
int idx = offset++;
V_p.col(idx) = V.col(i);
N_p.col(idx) = N.col(i);
A_p[idx] = A[i];
to_upper.col(idx) << i, -1;
to_lower[i] = idx;
}
}
adj_p.resize(V_p.cols());
std::vector<int> capacity(V_p.cols());
std::vector<std::vector<Link>> scratches(V_p.cols());
#ifdef WITH_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < V_p.cols(); ++i) {
int t = 0;
for (int j = 0; j < 2; ++j) {
int upper = to_upper(j, i);
if (upper == -1) continue;
t += adj[upper].size();
}
scratches[i].reserve(t);
adj_p[i].reserve(t);
}
#ifdef WITH_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < V_p.cols(); ++i) {
auto& scratch = scratches[i];
for (int j = 0; j < 2; ++j) {
int upper = to_upper(j, i);
if (upper == -1) continue;
auto& ad = adj[upper];
for (auto& link : ad) scratch.push_back(Link(to_lower[link.id], link.weight));
}
std::sort(scratch.begin(), scratch.end());
int id = -1;
auto& ad = adj_p[i];
for (auto& link : scratch) {
if (link.id != i) {
if (id != link.id) {
ad.push_back(link);
id = link.id;
} else {
ad.back().weight += link.weight;
}
}
}
}
}
void Hierarchy::SaveToFile(FILE* fp) {
Save(fp, mScale);
Save(fp, mF);
Save(fp, mE2E);
Save(fp, mAdj);
Save(fp, mV);
Save(fp, mN);
Save(fp, mA);
Save(fp, mToLower);
Save(fp, mToUpper);
Save(fp, mQ);
Save(fp, mO);
Save(fp, mS);
Save(fp, mK);
Save(fp, this->mPhases);
}
void Hierarchy::LoadFromFile(FILE* fp) {
Read(fp, mScale);
Read(fp, mF);
Read(fp, mE2E);
Read(fp, mAdj);
Read(fp, mV);
Read(fp, mN);
Read(fp, mA);
Read(fp, mToLower);
Read(fp, mToUpper);
Read(fp, mQ);
Read(fp, mO);
Read(fp, mS);
Read(fp, mK);
Read(fp, this->mPhases);
}
void Hierarchy::UpdateGraphValue(std::vector<Vector3i>& FQ, std::vector<Vector3i>& F2E,
std::vector<Vector2i>& edge_diff) {
FQ = std::move(mFQ[0]);
F2E = std::move(mF2E[0]);
edge_diff = std::move(mEdgeDiff[0]);
}
void Hierarchy::DownsampleEdgeGraph(std::vector<Vector3i>& FQ, std::vector<Vector3i>& F2E,
std::vector<Vector2i>& edge_diff,
std::vector<int>& allow_changes, int level) {
std::vector<Vector2i> E2F(edge_diff.size(), Vector2i(-1, -1));
for (int i = 0; i < F2E.size(); ++i) {
for (int j = 0; j < 3; ++j) {
int e = F2E[i][j];
if (E2F[e][0] == -1)
E2F[e][0] = i;
else
E2F[e][1] = i;
}
}
int levels = (level == -1) ? 100 : level;
mFQ.resize(levels);
mF2E.resize(levels);
mE2F.resize(levels);
mEdgeDiff.resize(levels);
mAllowChanges.resize(levels);
mSing.resize(levels);
mToUpperEdges.resize(levels - 1);
mToUpperOrients.resize(levels - 1);
for (int i = 0; i < FQ.size(); ++i) {
Vector2i diff(0, 0);
for (int j = 0; j < 3; ++j) {
diff += rshift90(edge_diff[F2E[i][j]], FQ[i][j]);
}
if (diff != Vector2i::Zero()) {
mSing[0].push_back(i);
}
}
mAllowChanges[0] = allow_changes;
mFQ[0] = std::move(FQ);
mF2E[0] = std::move(F2E);
mE2F[0] = std::move(E2F);
mEdgeDiff[0] = std::move(edge_diff);
for (int l = 0; l < levels - 1; ++l) {
auto& FQ = mFQ[l];
auto& E2F = mE2F[l];
auto& F2E = mF2E[l];
auto& Allow = mAllowChanges[l];
auto& EdgeDiff = mEdgeDiff[l];
auto& Sing = mSing[l];
std::vector<int> fixed_faces(F2E.size(), 0);
for (auto& s : Sing) {
fixed_faces[s] = 1;
}
auto& toUpper = mToUpperEdges[l];
auto& toUpperOrients = mToUpperOrients[l];
toUpper.resize(E2F.size(), -1);
toUpperOrients.resize(E2F.size(), 0);
auto& nFQ = mFQ[l + 1];
auto& nE2F = mE2F[l + 1];
auto& nF2E = mF2E[l + 1];
auto& nAllow = mAllowChanges[l + 1];
auto& nEdgeDiff = mEdgeDiff[l + 1];
auto& nSing = mSing[l + 1];
for (int i = 0; i < E2F.size(); ++i) {
if (EdgeDiff[i] != Vector2i::Zero()) continue;
if ((E2F[i][0] >= 0 && fixed_faces[E2F[i][0]]) ||
(E2F[i][1] >= 0 && fixed_faces[E2F[i][1]])) {
continue;
}
for (int j = 0; j < 2; ++j) {
int f = E2F[i][j];
if (f < 0)
continue;
for (int k = 0; k < 3; ++k) {
int neighbor_e = F2E[f][k];
for (int m = 0; m < 2; ++m) {
int neighbor_f = E2F[neighbor_e][m];
if (neighbor_f < 0)
continue;
if (fixed_faces[neighbor_f] == 0) fixed_faces[neighbor_f] = 1;
}
}
}
if (E2F[i][0] >= 0)
fixed_faces[E2F[i][0]] = 2;
if (E2F[i][1] >= 0)
fixed_faces[E2F[i][1]] = 2;
toUpper[i] = -2;
}
for (int i = 0; i < E2F.size(); ++i) {
if (toUpper[i] == -2) continue;
if ((E2F[i][0] < 0 || fixed_faces[E2F[i][0]] == 2) && (E2F[i][1] < 0 || fixed_faces[E2F[i][1]] == 2)) {
toUpper[i] = -3;
continue;
}
}
int numE = 0;
for (int i = 0; i < toUpper.size(); ++i) {
if (toUpper[i] == -1) {
if ((E2F[i][0] < 0 || fixed_faces[E2F[i][0]] < 2) && (E2F[i][1] < 0 || fixed_faces[E2F[i][1]] < 2)) {
nE2F.push_back(E2F[i]);
toUpperOrients[i] = 0;
toUpper[i] = numE++;
continue;
}
int f0 = (E2F[i][1] < 0 || fixed_faces[E2F[i][0]] < 2) ? E2F[i][0] : E2F[i][1];
int e = i;
int f = f0;
std::vector<std::pair<int, int>> paths;
paths.push_back(std::make_pair(i, 0));
while (true) {
if (E2F[e][0] == f)
f = E2F[e][1];
else if (E2F[e][1] == f)
f = E2F[e][0];
if (f < 0 || fixed_faces[f] < 2) {
for (int j = 0; j < paths.size(); ++j) {
auto& p = paths[j];
toUpper[p.first] = numE;
int orient = p.second;
if (j > 0) orient = (orient + toUpperOrients[paths[j - 1].first]) % 4;
toUpperOrients[p.first] = orient;
}
nE2F.push_back(Vector2i(f0, f));
numE += 1;
break;
}
int ind0 = -1, ind1 = -1;
int e0 = e;
for (int j = 0; j < 3; ++j) {
if (F2E[f][j] == e) {
ind0 = j;
break;
}
}
for (int j = 0; j < 3; ++j) {
int e1 = F2E[f][j];
if (e1 != e && toUpper[e1] != -2) {
e = e1;
ind1 = j;
break;
}
}
if (ind1 != -1) {
paths.push_back(std::make_pair(e, (FQ[f][ind1] - FQ[f][ind0] + 6) % 4));
} else {
if (EdgeDiff[e] != Vector2i::Zero()) {
printf("Unsatisfied !!!...\n");
printf("%d %d %d: %d %d\n", F2E[f][0], F2E[f][1], F2E[f][2], e0, e);
exit(0);
}
for (auto& p : paths) {
toUpper[p.first] = numE;
toUpperOrients[p.first] = 0;
}
numE += 1;
nE2F.push_back(Vector2i(f0, f0));
break;
}
}
}
}
nEdgeDiff.resize(numE);
nAllow.resize(numE * 2, 1);
for (int i = 0; i < toUpper.size(); ++i) {
if (toUpper[i] >= 0 && toUpperOrients[i] == 0) {
nEdgeDiff[toUpper[i]] = EdgeDiff[i];
}
if (toUpper[i] >= 0) {
int dimension = toUpperOrients[i] % 2;
if (Allow[i * 2 + dimension] == 0)
nAllow[toUpper[i] * 2] = 0;
else if (Allow[i * 2 + dimension] == 2)
nAllow[toUpper[i] * 2] = 2;
if (Allow[i * 2 + 1 - dimension] == 0)
nAllow[toUpper[i] * 2 + 1] = 0;
else if (Allow[i * 2 + 1 - dimension] == 2)
nAllow[toUpper[i] * 2 + 1] = 2;
}
}
std::vector<int> upperface(F2E.size(), -1);
for (int i = 0; i < F2E.size(); ++i) {
Vector3i eid;
for (int j = 0; j < 3; ++j) {
eid[j] = toUpper[F2E[i][j]];
}
if (eid[0] >= 0 && eid[1] >= 0 && eid[2] >= 0) {
Vector3i eid_orient;
for (int j = 0; j < 3; ++j) {
eid_orient[j] = (FQ[i][j] + 4 - toUpperOrients[F2E[i][j]]) % 4;
}
upperface[i] = nF2E.size();
nF2E.push_back(eid);
nFQ.push_back(eid_orient);
}
}
for (int i = 0; i < nE2F.size(); ++i) {
for (int j = 0; j < 2; ++j) {
if (nE2F[i][j] >= 0)
nE2F[i][j] = upperface[nE2F[i][j]];
}
}
for (auto& s : Sing) {
if (upperface[s] >= 0) nSing.push_back(upperface[s]);
}
mToUpperFaces.push_back(std::move(upperface));
if (nEdgeDiff.size() == EdgeDiff.size()) {
levels = l + 1;
break;
}
}
mFQ.resize(levels);
mF2E.resize(levels);
mAllowChanges.resize(levels);
mE2F.resize(levels);
mEdgeDiff.resize(levels);
mSing.resize(levels);
mToUpperEdges.resize(levels - 1);
mToUpperOrients.resize(levels - 1);
}
int Hierarchy::FixFlipSat(int depth, int threshold) {
if (system("which minisat > /dev/null 2>&1")) {
printf("minisat not found, \"-sat\" will not be used!\n");
return 0;
}
if (system("which timeout > /dev/null 2>&1")) {
printf("timeout not found, \"-sat\" will not be used!\n");
return 0;
}
auto& F2E = mF2E[depth];
auto& E2F = mE2F[depth];
auto& FQ = mFQ[depth];
auto& EdgeDiff = mEdgeDiff[depth];
auto& AllowChanges = mAllowChanges[depth];
// build E2E
std::vector<int> E2E(F2E.size() * 3, -1);
for (int i = 0; i < E2F.size(); ++i) {
int f1 = E2F[i][0];
int f2 = E2F[i][1];
int t1 = 0;
int t2 = 2;
if (f1 != -1) while (F2E[f1][t1] != i) t1 += 1;
if (f2 != -1) while (F2E[f2][t2] != i) t2 -= 1;
t1 += f1 * 3;
t2 += f2 * 3;
if (f1 != -1) E2E[t1] = (f2 == -1) ? -1 : t2;
if (f2 != -1) E2E[t2] = (f1 == -1) ? -1 : t1;
}
auto IntegerArea = [&](int f) {
Vector2i diff1 = rshift90(EdgeDiff[F2E[f][0]], FQ[f][0]);
Vector2i diff2 = rshift90(EdgeDiff[F2E[f][1]], FQ[f][1]);
return diff1[0] * diff2[1] - diff1[1] * diff2[0];
};
std::deque<std::pair<int, int>> Q;
std::vector<bool> mark_dedges(F2E.size() * 3, false);
for (int f = 0; f < F2E.size(); ++f) {
if (IntegerArea(f) < 0) {
for (int j = 0; j < 3; ++j) {
if (mark_dedges[f * 3 + j]) continue;
Q.push_back(std::make_pair(f * 3 + j, 0));
mark_dedges[f * 3 + j] = true;
}
}
}
int mark_count = 0;
while (!Q.empty()) {
int e0 = Q.front().first;
int depth = Q.front().second;
Q.pop_front();
mark_count++;
int e = e0, e1;
do {
e1 = E2E[e];
if (e1 == -1) break;
int length = EdgeDiff[F2E[e1 / 3][e1 % 3]].array().abs().sum();
if (length == 0 && !mark_dedges[e1]) {
mark_dedges[e1] = true;
Q.push_front(std::make_pair(e1, depth));
}
e = (e1 / 3) * 3 + (e1 + 1) % 3;
mark_dedges[e] = true;
} while (e != e0);
if (e1 == -1) {
do {
e1 = E2E[e];
if (e1 == -1) break;
int length = EdgeDiff[F2E[e1 / 3][e1 % 3]].array().abs().sum();
if (length == 0 && !mark_dedges[e1]) {
mark_dedges[e1] = true;
Q.push_front(std::make_pair(e1, depth));
}
e = (e1 / 3) * 3 + (e1 + 2) % 3;
mark_dedges[e] = true;
} while (e != e0);
}
do {
e1 = E2E[e];
if (e1 == -1) break;
int length = EdgeDiff[F2E[e1 / 3][e1 % 3]].array().abs().sum();
if (length > 0 && depth + length <= threshold && !mark_dedges[e1]) {
mark_dedges[e1] = true;
Q.push_back(std::make_pair(e1, depth + length));
}
e = e1 / 3 * 3 + (e1 + 1) % 3;
mark_dedges[e] = true;
} while (e != e0);
if (e1 == -1) {
do {
e1 = E2E[e];
if (e1 == -1) break;
int length = EdgeDiff[F2E[e1 / 3][e1 % 3]].array().abs().sum();
if (length > 0 && depth + length <= threshold && !mark_dedges[e1]) {
mark_dedges[e1] = true;
Q.push_back(std::make_pair(e1, depth + length));
}
e = e1 / 3 * 3 + (e1 + 2) % 3;
mark_dedges[e] = true;
} while (e != e0);
}
}
lprintf("[FlipH] Depth %2d: marked = %d\n", depth, mark_count);
std::vector<bool> flexible(EdgeDiff.size(), false);
for (int i = 0; i < F2E.size(); ++i) {
for (int j = 0; j < 3; ++j) {
int dedge = i * 3 + j;
int edgeid = F2E[i][j];
if (mark_dedges[dedge]) {
flexible[edgeid] = true;
}
}
}
for (int i = 0; i < flexible.size(); ++i) {
if (E2F[i][0] == E2F[i][1]) flexible[i] = false;
if (AllowChanges[i] == 0) flexible[i] = false;
}
// Reindexing and solve
int num_group = 0;
std::vector<int> groups(EdgeDiff.size(), -1);
std::vector<int> indices(EdgeDiff.size(), -1);
for (int i = 0; i < EdgeDiff.size(); ++i) {
if (groups[i] == -1 && flexible[i]) {
// group it
std::queue<int> q;
q.push(i);
groups[i] = num_group;
while (!q.empty()) {
int e = q.front();
q.pop();
int f[] = {E2F[e][0], E2F[e][1]};
for (int j = 0; j < 2; ++j) {
if (f[j] == -1) continue;
for (int k = 0; k < 3; ++k) {
int e1 = F2E[f[j]][k];
if (flexible[e1] && groups[e1] == -1) {
groups[e1] = num_group;
q.push(e1);
}
}
}
}
num_group += 1;
}
}
std::vector<int> num_edges(num_group);
std::vector<int> num_flips(num_group);
std::vector<std::vector<int>> values(num_group);
std::vector<std::vector<Vector3i>> variable_eq(num_group);
std::vector<std::vector<Vector3i>> constant_eq(num_group);
std::vector<std::vector<Vector4i>> variable_ge(num_group);
std::vector<std::vector<Vector2i>> constant_ge(num_group);
for (int i = 0; i < groups.size(); ++i) {
if (groups[i] != -1) {
indices[i] = num_edges[groups[i]]++;
values[groups[i]].push_back(EdgeDiff[i][0]);
values[groups[i]].push_back(EdgeDiff[i][1]);
}
}
std::vector<int> num_edges_flexible = num_edges;
std::map<std::pair<int, int>, int> fixed_variables;
for (int i = 0; i < F2E.size(); ++i) {
Vector2i var[3];
Vector2i cst[3];
int gind = 0;
while (gind < 3 && groups[F2E[i][gind]] == -1) gind += 1;
if (gind == 3) continue;
int group = groups[F2E[i][gind]];
int ind[3] = {-1, -1, -1};
for (int j = 0; j < 3; ++j) {
int g = groups[F2E[i][j]];
if (g != group) {
if (g == -1) {
auto key = std::make_pair(F2E[i][j], group);
auto it = fixed_variables.find(key);
if (it == fixed_variables.end()) {
ind[j] = num_edges[group];
values[group].push_back(EdgeDiff[F2E[i][j]][0]);
values[group].push_back(EdgeDiff[F2E[i][j]][1]);
fixed_variables[key] = num_edges[group]++;
} else {
ind[j] = it->second;
}
}
} else {
ind[j] = indices[F2E[i][j]];
}
}
for (int j = 0; j < 3; ++j) assert(ind[j] != -1);
for (int j = 0; j < 3; ++j) {
var[j] = rshift90(Vector2i(ind[j] * 2 + 1, ind[j] * 2 + 2), FQ[i][j]);
cst[j] = var[j].array().sign();
var[j] = var[j].array().abs() - 1;
}
num_flips[group] += IntegerArea(i) < 0;
variable_eq[group].push_back(Vector3i(var[0][0], var[1][0], var[2][0]));
constant_eq[group].push_back(Vector3i(cst[0][0], cst[1][0], cst[2][0]));
variable_eq[group].push_back(Vector3i(var[0][1], var[1][1], var[2][1]));
constant_eq[group].push_back(Vector3i(cst[0][1], cst[1][1], cst[2][1]));
variable_ge[group].push_back(Vector4i(var[0][0], var[1][1], var[0][1], var[1][0]));
constant_ge[group].push_back(Vector2i(cst[0][0] * cst[1][1], cst[0][1] * cst[1][0]));
}
int flip_before = 0, flip_after = 0;
for (int i = 0; i < F2E.size(); ++i) {
int area = IntegerArea(i);
if (area < 0) flip_before++;
}
for (int i = 0; i < num_group; ++i) {
std::vector<bool> flexible(values[i].size(), true);
for (int j = num_edges_flexible[i] * 2; j < flexible.size(); ++j) {
flexible[j] = false;
}
SolveSatProblem(values[i].size(), values[i], flexible, variable_eq[i], constant_eq[i],
variable_ge[i], constant_ge[i]);
}
for (int i = 0; i < EdgeDiff.size(); ++i) {
int group = groups[i];
if (group == -1) continue;
EdgeDiff[i][0] = values[group][2 * indices[i] + 0];
EdgeDiff[i][1] = values[group][2 * indices[i] + 1];
}
for (int i = 0; i < F2E.size(); ++i) {
Vector2i diff(0, 0);
for (int j = 0; j < 3; ++j) {
diff += rshift90(EdgeDiff[F2E[i][j]], FQ[i][j]);
}
assert(diff == Vector2i::Zero());
int area = IntegerArea(i);
if (area < 0) flip_after++;
}
lprintf("[FlipH] FlipArea, Before: %d After %d\n", flip_before, flip_after);
return flip_after;
}
void Hierarchy::PushDownwardFlip(int depth) {
auto& EdgeDiff = mEdgeDiff[depth];
auto& nEdgeDiff = mEdgeDiff[depth - 1];
auto& toUpper = mToUpperEdges[depth - 1];
auto& toUpperOrients = mToUpperOrients[depth - 1];
auto& toUpperFaces = mToUpperFaces[depth - 1];
for (int i = 0; i < toUpper.size(); ++i) {
if (toUpper[i] >= 0) {
int orient = (4 - toUpperOrients[i]) % 4;
nEdgeDiff[i] = rshift90(EdgeDiff[toUpper[i]], orient);
} else {
nEdgeDiff[i] = Vector2i(0, 0);
}
}
auto& nF2E = mF2E[depth - 1];
auto& nFQ = mFQ[depth - 1];
for (int i = 0; i < nF2E.size(); ++i) {
Vector2i diff(0, 0);
for (int j = 0; j < 3; ++j) {
diff += rshift90(nEdgeDiff[nF2E[i][j]], nFQ[i][j]);
}
if (diff != Vector2i::Zero()) {
printf("Fail!!!!!!! %d\n", i);
for (int j = 0; j < 3; ++j) {
Vector2i d = rshift90(nEdgeDiff[nF2E[i][j]], nFQ[i][j]);
printf("<%d %d %d>\n", nF2E[i][j], nFQ[i][j], toUpperOrients[nF2E[i][j]]);
printf("%d %d\n", d[0], d[1]);
printf("%d -> %d\n", nF2E[i][j], toUpper[nF2E[i][j]]);
}
printf("%d -> %d\n", i, toUpperFaces[i]);
exit(1);
}
}
}
void Hierarchy::FixFlip() {
int l = mF2E.size() - 1;
auto& F2E = mF2E[l];
auto& E2F = mE2F[l];
auto& FQ = mFQ[l];
auto& EdgeDiff = mEdgeDiff[l];
auto& AllowChange = mAllowChanges[l];
// build E2E
std::vector<int> E2E(F2E.size() * 3, -1);
for (int i = 0; i < E2F.size(); ++i) {
int v1 = E2F[i][0];
int v2 = E2F[i][1];
int t1 = 0;
int t2 = 2;
if (v1 != -1)
while (F2E[v1][t1] != i) t1 += 1;
if (v2 != -1)
while (F2E[v2][t2] != i) t2 -= 1;
t1 += v1 * 3;
t2 += v2 * 3;
if (v1 != -1)
E2E[t1] = (v2 == -1) ? -1 : t2;
if (v2 != -1)
E2E[t2] = (v1 == -1) ? -1 : t1;
}
auto Area = [&](int f) {
Vector2i diff1 = rshift90(EdgeDiff[F2E[f][0]], FQ[f][0]);
Vector2i diff2 = rshift90(EdgeDiff[F2E[f][1]], FQ[f][1]);
return diff1[0] * diff2[1] - diff1[1] * diff2[0];
};
std::vector<int> valences(F2E.size() * 3, -10000); // comment this line
auto CheckShrink = [&](int deid, int allowed_edge_length) {
// Check if we want shrink direct edge deid so that all edge length is smaller than
// allowed_edge_length
if (deid == -1) {
return false;
}
std::vector<int> corresponding_faces;
std::vector<int> corresponding_edges;
std::vector<Vector2i> corresponding_diff;
int deid0 = deid;
while (deid != -1) {
deid = deid / 3 * 3 + (deid + 2) % 3;
if (E2E[deid] == -1)
break;
deid = E2E[deid];
if (deid == deid0)
break;
}
Vector2i diff = EdgeDiff[F2E[deid / 3][deid % 3]];
do {
corresponding_diff.push_back(diff);
corresponding_edges.push_back(deid);
corresponding_faces.push_back(deid / 3);
// transform to the next face
deid = E2E[deid];
if (deid == -1) {
return false;
}
// transform for the target incremental diff
diff = -rshift90(diff, FQ[deid / 3][deid % 3]);
deid = deid / 3 * 3 + (deid + 1) % 3;
// transform to local
diff = rshift90(diff, (4 - FQ[deid / 3][deid % 3]) % 4);
} while (deid != corresponding_edges.front());
// check diff
if (deid != -1 && diff != corresponding_diff.front()) {
return false;
}
std::unordered_map<int, Vector2i> new_values;
for (int i = 0; i < corresponding_diff.size(); ++i) {
int deid = corresponding_edges[i];
int eid = F2E[deid / 3][deid % 3];
new_values[eid] = EdgeDiff[eid];
}
for (int i = 0; i < corresponding_diff.size(); ++i) {
int deid = corresponding_edges[i];
int eid = F2E[deid / 3][deid % 3];
for (int j = 0; j < 2; ++j) {
if (corresponding_diff[i][j] != 0 && AllowChange[eid * 2 + j] == 0) return false;
}
auto& res = new_values[eid];
res -= corresponding_diff[i];
int edge_thres = allowed_edge_length;
if (abs(res[0]) > edge_thres || abs(res[1]) > edge_thres) {
return false;
}
if ((abs(res[0]) > 1 && abs(res[1]) != 0) || (abs(res[1]) > 1 && abs(res[0]) != 0))
return false;
}
int prev_area = 0, current_area = 0;
for (int f = 0; f < corresponding_faces.size(); ++f) {
int area = Area(corresponding_faces[f]);
if (area < 0) prev_area += 1;
}
for (auto& p : new_values) {
std::swap(EdgeDiff[p.first], p.second);
}
for (int f = 0; f < corresponding_faces.size(); ++f) {
int area = Area(corresponding_faces[f]);
if (area < 0) {
current_area += 1;
}
}
if (current_area < prev_area) {
return true;
}
for (auto& p : new_values) {
std::swap(EdgeDiff[p.first], p.second);
}
return false;
};
std::queue<int> flipped;
for (int i = 0; i < F2E.size(); ++i) {
int area = Area(i);
if (area < 0) {
flipped.push(i);
}
}
bool update = false;
int max_len = 1;
while (!update && max_len <= 2) {
while (!flipped.empty()) {
int f = flipped.front();
if (Area(f) >= 0) {
flipped.pop();
continue;
}
for (int i = 0; i < 3; ++i) {
if (CheckShrink(f * 3 + i, max_len) || CheckShrink(E2E[f * 3 + i], max_len)) {
update = true;
break;
}
}
flipped.pop();
}
max_len += 1;
}
if (update) {
Hierarchy flip_hierarchy;
flip_hierarchy.DownsampleEdgeGraph(mFQ.back(), mF2E.back(), mEdgeDiff.back(),
mAllowChanges.back(), -1);
flip_hierarchy.FixFlip();
flip_hierarchy.UpdateGraphValue(mFQ.back(), mF2E.back(), mEdgeDiff.back());
}
PropagateEdge();
}
void Hierarchy::PropagateEdge() {
for (int level = mToUpperEdges.size(); level > 0; --level) {
auto& EdgeDiff = mEdgeDiff[level];
auto& nEdgeDiff = mEdgeDiff[level - 1];
auto& FQ = mFQ[level];
auto& nFQ = mFQ[level - 1];
auto& F2E = mF2E[level - 1];
auto& toUpper = mToUpperEdges[level - 1];
auto& toUpperFace = mToUpperFaces[level - 1];
auto& toUpperOrients = mToUpperOrients[level - 1];
for (int i = 0; i < toUpper.size(); ++i) {
if (toUpper[i] >= 0) {
int orient = (4 - toUpperOrients[i]) % 4;
nEdgeDiff[i] = rshift90(EdgeDiff[toUpper[i]], orient);
} else {
nEdgeDiff[i] = Vector2i(0, 0);
}
}
for (int i = 0; i < toUpperFace.size(); ++i) {
if (toUpperFace[i] == -1) continue;
Vector3i eid_orient = FQ[toUpperFace[i]];
for (int j = 0; j < 3; ++j) {
nFQ[i][j] = (eid_orient[j] + toUpperOrients[F2E[i][j]]) % 4;
}
}
}
}
void Hierarchy::clearConstraints() {
int levels = mV.size();
if (levels == 0) return;
for (int i = 0; i < levels; ++i) {
int size = mV[i].cols();
mCQ[i].resize(3, size);
mCO[i].resize(3, size);
mCQw[i].resize(size);
mCOw[i].resize(size);
mCQw[i].setZero();
mCOw[i].setZero();
}
}
void Hierarchy::propagateConstraints() {
int levels = mV.size();
if (levels == 0) return;
for (int l = 0; l < levels - 1; ++l) {
auto& N = mN[l];
auto& N_next = mN[l + 1];
auto& V = mV[l];
auto& V_next = mV[l + 1];
auto& CQ = mCQ[l];
auto& CQ_next = mCQ[l + 1];
auto& CQw = mCQw[l];
auto& CQw_next = mCQw[l + 1];
auto& CO = mCO[l];
auto& CO_next = mCO[l + 1];
auto& COw = mCOw[l];
auto& COw_next = mCOw[l + 1];
auto& toUpper = mToUpper[l];
MatrixXd& S = mS[l];
for (uint32_t i = 0; i != mV[l + 1].cols(); ++i) {
Vector2i upper = toUpper.col(i);
Vector3d cq = Vector3d::Zero(), co = Vector3d::Zero();
float cqw = 0.0f, cow = 0.0f;
bool has_cq0 = CQw[upper[0]] != 0;
bool has_cq1 = upper[1] != -1 && CQw[upper[1]] != 0;
bool has_co0 = COw[upper[0]] != 0;
bool has_co1 = upper[1] != -1 && COw[upper[1]] != 0;
if (has_cq0 && !has_cq1) {
cq = CQ.col(upper[0]);
cqw = CQw[upper[0]];
} else if (has_cq1 && !has_cq0) {
cq = CQ.col(upper[1]);
cqw = CQw[upper[1]];
} else if (has_cq1 && has_cq0) {
Vector3d q_i = CQ.col(upper[0]);
Vector3d n_i = CQ.col(upper[0]);
Vector3d q_j = CQ.col(upper[1]);
Vector3d n_j = CQ.col(upper[1]);
auto result = compat_orientation_extrinsic_4(q_i, n_i, q_j, n_j);
cq = result.first * CQw[upper[0]] + result.second * CQw[upper[1]];
cqw = (CQw[upper[0]] + CQw[upper[1]]);
}
if (cq != Vector3d::Zero()) {
Vector3d n = N_next.col(i);
cq -= n.dot(cq) * n;
if (cq.squaredNorm() > RCPOVERFLOW) cq.normalize();
}
if (has_co0 && !has_co1) {
co = CO.col(upper[0]);
cow = COw[upper[0]];
} else if (has_co1 && !has_co0) {
co = CO.col(upper[1]);
cow = COw[upper[1]];
} else if (has_co1 && has_co0) {
double scale_x = mScale;
double scale_y = mScale;
if (with_scale) {
// FIXME
// scale_x *= S(0, i);
// scale_y *= S(1, i);
}
double inv_scale_x = 1.0f / scale_x;
double inv_scale_y = 1.0f / scale_y;
double scale_x_1 = mScale;
double scale_y_1 = mScale;
if (with_scale) {
// FIXME
// scale_x_1 *= S(0, j);
// scale_y_1 *= S(1, j);
}
double inv_scale_x_1 = 1.0f / scale_x_1;
double inv_scale_y_1 = 1.0f / scale_y_1;
auto result = compat_position_extrinsic_4(
V.col(upper[0]), N.col(upper[0]), CQ.col(upper[0]), CO.col(upper[0]),
V.col(upper[1]), N.col(upper[1]), CQ.col(upper[1]), CO.col(upper[1]), scale_x,
scale_y, inv_scale_x, inv_scale_y, scale_x_1, scale_y_1, inv_scale_x_1,
inv_scale_y_1);
cow = COw[upper[0]] + COw[upper[1]];
co = (result.first * COw[upper[0]] + result.second * COw[upper[1]]) / cow;
}
if (co != Vector3d::Zero()) {
Vector3d n = N_next.col(i), v = V_next.col(i);
co -= n.dot(cq - v) * n;
}
#if 0
cqw *= 0.5f;
cow *= 0.5f;
#else
if (cqw > 0) cqw = 1;
if (cow > 0) cow = 1;
#endif
CQw_next[i] = cqw;
COw_next[i] = cow;
CQ_next.col(i) = cq;
CO_next.col(i) = co;
}
}
}
#ifdef WITH_CUDA
#include <cuda_runtime.h>
void Hierarchy::CopyToDevice() {
if (cudaAdj.empty()) {
cudaAdj.resize(mAdj.size());
cudaAdjOffset.resize(mAdj.size());
for (int i = 0; i < mAdj.size(); ++i) {
std::vector<int> offset(mAdj[i].size() + 1, 0);
for (int j = 0; j < mAdj[i].size(); ++j) {
offset[j + 1] = offset[j] + mAdj[i][j].size();
}
cudaMalloc(&cudaAdjOffset[i], sizeof(int) * (mAdj[i].size() + 1));
cudaMemcpy(cudaAdjOffset[i], offset.data(), sizeof(int) * (mAdj[i].size() + 1),
cudaMemcpyHostToDevice);
// cudaAdjOffset[i] = (int*)malloc(sizeof(int) * (mAdj[i].size() + 1));
// memcpy(cudaAdjOffset[i], offset.data(), sizeof(int) * (mAdj[i].size() +
// 1));
cudaMalloc(&cudaAdj[i], sizeof(Link) * offset.back());
// cudaAdj[i] = (Link*)malloc(sizeof(Link) * offset.back());
std::vector<Link> plainlink(offset.back());
for (int j = 0; j < mAdj[i].size(); ++j) {
memcpy(plainlink.data() + offset[j], mAdj[i][j].data(),
mAdj[i][j].size() * sizeof(Link));
}
cudaMemcpy(cudaAdj[i], plainlink.data(), plainlink.size() * sizeof(Link),
cudaMemcpyHostToDevice);
}
}
if (cudaN.empty()) {
cudaN.resize(mN.size());
for (int i = 0; i < mN.size(); ++i) {
cudaMalloc(&cudaN[i], sizeof(glm::dvec3) * mN[i].cols());
// cudaN[i] = (glm::dvec3*)malloc(sizeof(glm::dvec3) * mN[i].cols());
}
}
for (int i = 0; i < mN.size(); ++i) {
cudaMemcpy(cudaN[i], mN[i].data(), sizeof(glm::dvec3) * mN[i].cols(),
cudaMemcpyHostToDevice);
// memcpy(cudaN[i], mN[i].data(), sizeof(glm::dvec3) * mN[i].cols());
}
if (cudaV.empty()) {
cudaV.resize(mV.size());
for (int i = 0; i < mV.size(); ++i) {
cudaMalloc(&cudaV[i], sizeof(glm::dvec3) * mV[i].cols());
// cudaV[i] = (glm::dvec3*)malloc(sizeof(glm::dvec3) * mV[i].cols());
}
}
for (int i = 0; i < mV.size(); ++i) {
cudaMemcpy(cudaV[i], mV[i].data(), sizeof(glm::dvec3) * mV[i].cols(),
cudaMemcpyHostToDevice);
// memcpy(cudaV[i], mV[i].data(), sizeof(glm::dvec3) * mV[i].cols());
}
if (cudaQ.empty()) {
cudaQ.resize(mQ.size());
for (int i = 0; i < mQ.size(); ++i) {
cudaMalloc(&cudaQ[i], sizeof(glm::dvec3) * mQ[i].cols());
// cudaQ[i] = (glm::dvec3*)malloc(sizeof(glm::dvec3) * mQ[i].cols());
}
}
for (int i = 0; i < mQ.size(); ++i) {
cudaMemcpy(cudaQ[i], mQ[i].data(), sizeof(glm::dvec3) * mQ[i].cols(),
cudaMemcpyHostToDevice);
// memcpy(cudaQ[i], mQ[i].data(), sizeof(glm::dvec3) * mQ[i].cols());
}
if (cudaO.empty()) {
cudaO.resize(mO.size());
for (int i = 0; i < mO.size(); ++i) {
cudaMalloc(&cudaO[i], sizeof(glm::dvec3) * mO[i].cols());
// cudaO[i] = (glm::dvec3*)malloc(sizeof(glm::dvec3) * mO[i].cols());
}
}
for (int i = 0; i < mO.size(); ++i) {
cudaMemcpy(cudaO[i], mO[i].data(), sizeof(glm::dvec3) * mO[i].cols(),
cudaMemcpyHostToDevice);
// memcpy(cudaO[i], mO[i].data(), sizeof(glm::dvec3) * mO[i].cols());
}
if (cudaPhases.empty()) {
cudaPhases.resize(mPhases.size());
for (int i = 0; i < mPhases.size(); ++i) {
cudaPhases[i].resize(mPhases[i].size());
for (int j = 0; j < mPhases[i].size(); ++j) {
cudaMalloc(&cudaPhases[i][j], sizeof(int) * mPhases[i][j].size());
// cudaPhases[i][j] = (int*)malloc(sizeof(int) *
// mPhases[i][j].size());
}
}
}
for (int i = 0; i < mPhases.size(); ++i) {
for (int j = 0; j < mPhases[i].size(); ++j) {
cudaMemcpy(cudaPhases[i][j], mPhases[i][j].data(), sizeof(int) * mPhases[i][j].size(),
cudaMemcpyHostToDevice);
// memcpy(cudaPhases[i][j], mPhases[i][j].data(), sizeof(int) *
// mPhases[i][j].size());
}
}
if (cudaToUpper.empty()) {
cudaToUpper.resize(mToUpper.size());
for (int i = 0; i < mToUpper.size(); ++i) {
cudaMalloc(&cudaToUpper[i], mToUpper[i].cols() * sizeof(glm::ivec2));
// cudaToUpper[i] = (glm::ivec2*)malloc(mToUpper[i].cols() *
// sizeof(glm::ivec2));
}
}
for (int i = 0; i < mToUpper.size(); ++i) {
cudaMemcpy(cudaToUpper[i], mToUpper[i].data(), sizeof(glm::ivec2) * mToUpper[i].cols(),
cudaMemcpyHostToDevice);
// memcpy(cudaToUpper[i], mToUpper[i].data(), sizeof(glm::ivec2) *
// mToUpper[i].cols());
}
cudaDeviceSynchronize();
}
void Hierarchy::CopyToHost() {}
#endif
} // namespace qflow