Skip to content

Commit

Permalink
misc (#533)
Browse files Browse the repository at this point in the history
* use better policies

* removed redundant copy

* added condensed_matter test

* disable normal check for CondensedMatter16

* disable overlap checking as well
  • Loading branch information
pca006132 authored Aug 23, 2023
1 parent 111350f commit cc83e1b
Show file tree
Hide file tree
Showing 8 changed files with 278 additions and 101 deletions.
2 changes: 1 addition & 1 deletion samples/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ project(samples)

add_library(samples src/bracelet.cpp src/knot.cpp src/menger_sponge.cpp
src/rounded_frame.cpp src/scallop.cpp src/tet_puzzle.cpp
src/gyroid_module.cpp)
src/gyroid_module.cpp src/condensed_matter.cpp)

target_include_directories(samples
PUBLIC ${PROJECT_SOURCE_DIR}/include
Expand Down
2 changes: 2 additions & 0 deletions samples/include/samples.h
Original file line number Diff line number Diff line change
Expand Up @@ -51,5 +51,7 @@ Manifold TetPuzzle(float edgeLength, float gap, int nDivisions);
Manifold Scallop();

Manifold GyroidModule(float size = 20, int n = 20);

Manifold CondensedMatter(int fn = 16);
/** @} */ // end of Samples
} // namespace manifold
149 changes: 149 additions & 0 deletions samples/src/condensed_matter.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,149 @@
// Copyright 2023 The Manifold Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// From
// https://gist.github.com/ochafik/70a6b15e982b7ccd5a79ff9afd99dbcf#file-condensed-matter-scad

#include "samples.h"

namespace {
using namespace manifold;
using namespace glm;

constexpr float AtomicRadiusN2 = 0.65;
constexpr float BondPairN2 = 1.197;
constexpr float AtomicRadiusSi = 1.1;
constexpr float LatticeCellSizeSi = 5.4309;
constexpr float fccOffset = 0.25;
constexpr float AtomicRadiusC = 0.7;
constexpr float LatticeCellSizeC = 3.65;
constexpr float cellLenA = 2.464;
constexpr float cellLenB = cellLenA;
constexpr float cellLenC = 6.711;
constexpr float cellAngleA = 90;
constexpr float cellAngleB = cellAngleA;
constexpr float cellAngleC = 120;
constexpr float LayerSeperationC = 3.364;

Manifold bond(int fn, vec3 p1 = {0, 0, 0}, vec3 p2 = {1, 1, 1}, float ar1 = 1.0,
float ar2 = 2.0) {
float cyR = std::min(ar1, ar2) / 5.0;
float dist = length(p1 - p2);
vec3 cyC = (p1 + p2) / 2.0f;
float beta = degrees(acos((p1.z - p2.z) / dist));
float gamma = degrees(atan2(p1.y - p2.y, p1.x - p2.x));
vec3 rot = {0.0, beta, gamma};
return Manifold::Cylinder(dist, cyR, -1, fn, true)
.Rotate(rot.x, rot.y, rot.z)
.Translate(cyC);
}

Manifold bondPair(int fn, float d = 0.0, float ar = 1.0) {
float axD = pow(d, 1.0 / 3.0);
vec3 p1 = {+axD, -axD, -axD};
vec3 p2 = {-axD, +axD, +axD};
Manifold sphere = Manifold::Sphere(ar, fn);
return sphere.Translate(p1) + sphere.Translate(p2) + bond(fn, p1, p2, ar, ar);
}

Manifold hexagonalClosePacked(int fn, vec3 dst = {1.0, 1.0, 1.0},
float ar = 1.0) {
std::vector<Manifold> parts;
vec3 p1 = {0, 0, 0};
parts.push_back(Manifold::Sphere(ar, fn));
float baseAg = 30;
vec3 ag = {baseAg, baseAg + 120, baseAg + 240};
vec3 points[] = {{cosd(ag.x) * dst.x, sind(ag.x) * dst.x, 0},
{cosd(ag.y) * dst.y, sind(ag.y) * dst.y, 0},
{cosd(ag.z) * dst.z, sind(ag.z) * dst.z, 0}};
for (vec3 p2 : points) {
parts.push_back(Manifold::Sphere(ar, fn).Translate(p2));
parts.push_back(bond(fn, p1, p2, ar, ar));
}
return Manifold::BatchBoolean(parts, OpType::Add);
}

Manifold fccDiamond(int fn, float ar = 1.0, float unitCell = 2.0,
float fccOffset = 0.25) {
std::vector<Manifold> parts;
float huc = unitCell / 2.0;
float od = fccOffset * unitCell;
vec3 interstitial[] = {
{+od, +od, +od}, {+od, -od, -od}, {-od, +od, -od}, {-od, -od, +od}};
vec3 corners[] = {{+huc, +huc, +huc},
{+huc, -huc, -huc},
{-huc, +huc, -huc},
{-huc, -huc, +huc}};
vec3 fcc[] = {{+huc, 0, 0}, {-huc, 0, 0}, {0, +huc, 0},
{0, -huc, 0}, {0, 0, +huc}, {0, 0, -huc}};
for (auto p : corners) parts.push_back(Manifold::Sphere(ar, fn).Translate(p));

for (auto p : fcc) parts.push_back(Manifold::Sphere(ar, fn).Translate(p));
for (auto p : interstitial)
parts.push_back(Manifold::Sphere(ar, fn).Translate(p));

vec3 bonds[][2] = {{interstitial[0], corners[0]}, {interstitial[0], fcc[0]},
{interstitial[0], fcc[2]}, {interstitial[0], fcc[4]},
{interstitial[1], corners[1]}, {interstitial[1], fcc[0]},
{interstitial[1], fcc[3]}, {interstitial[1], fcc[5]},
{interstitial[2], corners[2]}, {interstitial[2], fcc[1]},
{interstitial[2], fcc[2]}, {interstitial[2], fcc[5]},
{interstitial[3], corners[3]}, {interstitial[3], fcc[1]},
{interstitial[3], fcc[3]}, {interstitial[3], fcc[4]}};
for (auto b : bonds) parts.push_back(bond(fn, b[0], b[1], ar, ar));

return Manifold::BatchBoolean(parts, OpType::Add);
}

Manifold SiCell(int fn, float x = 1.0, float y = 1.0, float z = 1.0) {
return fccDiamond(fn, AtomicRadiusSi, LatticeCellSizeSi, fccOffset)
.Translate({LatticeCellSizeSi * x, LatticeCellSizeSi * y,
LatticeCellSizeSi * z});
}

Manifold SiN2Cell(int fn, float x = 1.0, float y = 1.0, float z = 1.0) {
float n2Offset = LatticeCellSizeSi / 8;
return bondPair(fn, BondPairN2, AtomicRadiusN2)
.Translate({LatticeCellSizeSi * x - n2Offset,
LatticeCellSizeSi * y + n2Offset,
LatticeCellSizeSi * z + n2Offset}) +
SiCell(fn, x, y, z);
}

Manifold GraphiteCell(int fn, vec3 xyz = {1.0, 1.0, 1.0}) {
vec3 loc = {(cellLenA * xyz.x * cosd(30) * 2),
((cellLenB * sind(30)) + cellLenC) * xyz.y, xyz.z};
return hexagonalClosePacked(fn, {cellLenA, cellLenB, cellLenC}, AtomicRadiusC)
.Translate(loc);
}

} // namespace

namespace manifold {
Manifold CondensedMatter(int fn) {
std::vector<Manifold> parts;
float siOffset = 3.0 * LatticeCellSizeSi / 8.0;
for (int x = -3; x <= 3; x++)
for (int y = -1; y <= 2; y++)
parts.push_back(
GraphiteCell(fn, {x + (y % 2 == 0 ? 0.0 : 0.5), y,
LayerSeperationC * 0.5 + LatticeCellSizeSi * 1.5})
.Translate({0, -siOffset, 0})
.Rotate(0, 0, 45));

float xyPlane[] = {-2, -1, 0, +1, +2};
for (float x : xyPlane)
for (float y : xyPlane) parts.push_back(SiN2Cell(fn, x, y, 1));
return Manifold::BatchBoolean(parts, OpType::Add);
}
} // namespace manifold
77 changes: 38 additions & 39 deletions src/manifold/src/boolean3.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -91,18 +91,19 @@ struct CopyFaceEdges {
};

SparseIndices Filter11(const Manifold::Impl &inP, const Manifold::Impl &inQ,
const SparseIndices &p1q2, const SparseIndices &p2q1,
ExecutionPolicy policy) {
const SparseIndices &p1q2, const SparseIndices &p2q1) {
SparseIndices p1q1(3 * p1q2.size() + 3 * p2q1.size());
for_each_n(policy, zip(countAt(0), countAt(0)), p1q2.size(),
for_each_n(autoPolicy(p1q2.size()), zip(countAt(0), countAt(0)), p1q2.size(),
CopyFaceEdges<false>({p1q2, p1q1, inQ.halfedge_}));
for_each_n(policy, zip(countAt(p1q2.size()), countAt(0)), p2q1.size(),
CopyFaceEdges<true>({p2q1, p1q1, inP.halfedge_}));
p1q1.Unique(policy);
for_each_n(autoPolicy(p2q1.size()), zip(countAt(p1q2.size()), countAt(0)),
p2q1.size(), CopyFaceEdges<true>({p2q1, p1q1, inP.halfedge_}));
p1q1.Unique();
return p1q1;
}

bool Shadows(float p, float q, float dir) { return p == q ? dir < 0 : p < q; }
inline bool Shadows(float p, float q, float dir) {
return p == q ? dir < 0 : p < q;
}

/**
* Since this function is called from two different places, it is necessary that
Expand All @@ -115,7 +116,7 @@ bool Shadows(float p, float q, float dir) { return p == q ? dir < 0 : p < q; }
* compiled function (they must agree on CPU or GPU). This is now taken care of
* by the shared policy_ member.
*/
thrust::pair<int, glm::vec2> Shadow01(
inline thrust::pair<int, glm::vec2> Shadow01(
const int p0, const int q1, VecView<const glm::vec3> vertPosP,
VecView<const glm::vec3> vertPosQ, VecView<const Halfedge> halfedgeQ,
const float expandP, VecView<const glm::vec3> normalP, const bool reverse) {
Expand Down Expand Up @@ -271,12 +272,12 @@ struct Kernel11 {
std::tuple<Vec<int>, Vec<glm::vec4>> Shadow11(SparseIndices &p1q1,
const Manifold::Impl &inP,
const Manifold::Impl &inQ,
float expandP,
ExecutionPolicy policy) {
float expandP) {
Vec<int> s11(p1q1.size());
Vec<glm::vec4> xyzz11(p1q1.size());

for_each_n(policy, zip(countAt(0), xyzz11.begin(), s11.begin()), p1q1.size(),
for_each_n(autoPolicy(p1q1.size()),
zip(countAt(0), xyzz11.begin(), s11.begin()), p1q1.size(),
Kernel11({inP.vertPos_, inQ.vertPos_, inP.halfedge_, inQ.halfedge_,
expandP, inP.vertNormal_, p1q1}));

Expand All @@ -289,10 +290,10 @@ struct Kernel02 {
VecView<const glm::vec3> vertPosP;
VecView<const Halfedge> halfedgeQ;
VecView<const glm::vec3> vertPosQ;
const bool forward;
const float expandP;
VecView<const glm::vec3> vertNormalP;
const SparseIndices &p0q2;
const bool forward;

void operator()(thrust::tuple<int, int &, float &> inout) {
const int p0 = p0q2.Get(thrust::get<0>(inout), !forward);
Expand Down Expand Up @@ -365,15 +366,15 @@ struct Kernel02 {
std::tuple<Vec<int>, Vec<float>> Shadow02(const Manifold::Impl &inP,
const Manifold::Impl &inQ,
SparseIndices &p0q2, bool forward,
float expandP,
ExecutionPolicy policy) {
float expandP) {
Vec<int> s02(p0q2.size());
Vec<float> z02(p0q2.size());

auto vertNormalP = forward ? inP.vertNormal_ : inQ.vertNormal_;
for_each_n(policy, zip(countAt(0), s02.begin(), z02.begin()), p0q2.size(),
Kernel02({inP.vertPos_, inQ.halfedge_, inQ.vertPos_, forward,
expandP, vertNormalP, p0q2}));
for_each_n(autoPolicy(p0q2.size()), zip(countAt(0), s02.begin(), z02.begin()),
p0q2.size(),
Kernel02({inP.vertPos_, inQ.halfedge_, inQ.vertPos_, expandP,
vertNormalP, p0q2, forward}));

p0q2.KeepFinite(z02, s02);

Expand Down Expand Up @@ -474,12 +475,13 @@ std::tuple<Vec<int>, Vec<glm::vec3>> Intersect12(
const Manifold::Impl &inP, const Manifold::Impl &inQ, const Vec<int> &s02,
const SparseIndices &p0q2, const Vec<int> &s11, const SparseIndices &p1q1,
const Vec<float> &z02, const Vec<glm::vec4> &xyzz11, SparseIndices &p1q2,
bool forward, ExecutionPolicy policy) {
bool forward) {
Vec<int> x12(p1q2.size());
Vec<glm::vec3> v12(p1q2.size());

for_each_n(
policy, zip(countAt(0), x12.begin(), v12.begin()), p1q2.size(),
autoPolicy(p1q2.size()), zip(countAt(0), x12.begin(), v12.begin()),
p1q2.size(),
Kernel12({p0q2.AsVec64(), s02, z02, p1q1.AsVec64(), s11, xyzz11,
inP.halfedge_, inQ.halfedge_, inP.vertPos_, forward, p1q2}));

Expand All @@ -489,10 +491,11 @@ std::tuple<Vec<int>, Vec<glm::vec3>> Intersect12(
};

Vec<int> Winding03(const Manifold::Impl &inP, Vec<int> &vertices, Vec<int> &s02,
bool reverse, ExecutionPolicy policy) {
bool reverse) {
// verts that are not shadowed (not in p0q2) have winding number zero.
Vec<int> w03(inP.NumVert(), 0);
// checking is slow, so just sort and reduce
auto policy = autoPolicy(vertices.size());
stable_sort_by_key(policy, vertices.begin(), vertices.end(), s02.begin());
Vec<int> w03val(w03.size());
Vec<int> w03vert(w03.size());
Expand All @@ -513,10 +516,7 @@ Vec<int> Winding03(const Manifold::Impl &inP, Vec<int> &vertices, Vec<int> &s02,
namespace manifold {
Boolean3::Boolean3(const Manifold::Impl &inP, const Manifold::Impl &inQ,
OpType op)
: inP_(inP),
inQ_(inQ),
expandP_(op == OpType::Add ? 1.0 : -1.0),
policy_(autoPolicy(glm::max(inP.NumEdge(), inQ.NumEdge()))) {
: inP_(inP), inQ_(inQ), expandP_(op == OpType::Add ? 1.0 : -1.0) {
// Symbolic perturbation:
// Union -> expand inP
// Difference, Intersection -> contract inP
Expand All @@ -538,25 +538,24 @@ Boolean3::Boolean3(const Manifold::Impl &inP, const Manifold::Impl &inQ,
p1q2_ = inQ_.EdgeCollisions(inP_);
p2q1_ = inP_.EdgeCollisions(inQ_, true); // inverted

policy_ = autoPolicy(glm::max(p1q2_.size(), p2q1_.size()));
p1q2_.Sort(policy_);
p1q2_.Sort();
PRINT("p1q2 size = " << p1q2_.size());

p2q1_.Sort(policy_);
p2q1_.Sort();
PRINT("p2q1 size = " << p2q1_.size());

// Level 2
// Find vertices that overlap faces in XY-projection
SparseIndices p0q2 = inQ.VertexCollisionsZ(inP.vertPos_);
p0q2.Sort(policy_);
p0q2.Sort();
PRINT("p0q2 size = " << p0q2.size());

SparseIndices p2q0 = inP.VertexCollisionsZ(inQ.vertPos_, true); // inverted
p2q0.Sort(policy_);
p2q0.Sort();
PRINT("p2q0 size = " << p2q0.size());

// Find involved edge pairs from Level 3
SparseIndices p1q1 = Filter11(inP_, inQ_, p1q2_, p2q1_, policy_);
SparseIndices p1q1 = Filter11(inP_, inQ_, p1q2_, p2q1_);
PRINT("p1q1 size = " << p1q1.size());

#ifdef MANIFOLD_DEBUG
Expand All @@ -570,41 +569,41 @@ Boolean3::Boolean3(const Manifold::Impl &inP, const Manifold::Impl &inQ,
// each edge, keeping only those whose intersection exists.
Vec<int> s11;
Vec<glm::vec4> xyzz11;
std::tie(s11, xyzz11) = Shadow11(p1q1, inP, inQ, expandP_, policy_);
std::tie(s11, xyzz11) = Shadow11(p1q1, inP, inQ, expandP_);
PRINT("s11 size = " << s11.size());

// Build up Z-projection of vertices onto triangles, keeping only those that
// fall inside the triangle.
Vec<int> s02;
Vec<float> z02;
std::tie(s02, z02) = Shadow02(inP, inQ, p0q2, true, expandP_, policy_);
std::tie(s02, z02) = Shadow02(inP, inQ, p0q2, true, expandP_);
PRINT("s02 size = " << s02.size());

Vec<int> s20;
Vec<float> z20;
std::tie(s20, z20) = Shadow02(inQ, inP, p2q0, false, expandP_, policy_);
std::tie(s20, z20) = Shadow02(inQ, inP, p2q0, false, expandP_);
PRINT("s20 size = " << s20.size());

// Level 3
// Build up the intersection of the edges and triangles, keeping only those
// that intersect, and record the direction the edge is passing through the
// triangle.
std::tie(x12_, v12_) = Intersect12(inP, inQ, s02, p0q2, s11, p1q1, z02,
xyzz11, p1q2_, true, policy_);
std::tie(x12_, v12_) =
Intersect12(inP, inQ, s02, p0q2, s11, p1q1, z02, xyzz11, p1q2_, true);
PRINT("x12 size = " << x12_.size());

std::tie(x21_, v21_) = Intersect12(inQ, inP, s20, p2q0, s11, p1q1, z20,
xyzz11, p2q1_, false, policy_);
std::tie(x21_, v21_) =
Intersect12(inQ, inP, s20, p2q0, s11, p1q1, z20, xyzz11, p2q1_, false);
PRINT("x21 size = " << x21_.size());

Vec<int> p0 = p0q2.Copy(false);
p0q2.Resize(0);
Vec<int> q0 = p2q0.Copy(true);
p2q0.Resize(0);
// Sum up the winding numbers of all vertices.
w03_ = Winding03(inP, p0, s02, false, policy_);
w03_ = Winding03(inP, p0, s02, false);

w30_ = Winding03(inQ, q0, s20, true, policy_);
w30_ = Winding03(inQ, q0, s20, true);

#ifdef MANIFOLD_DEBUG
intersections.Stop();
Expand Down
1 change: 0 additions & 1 deletion src/manifold/src/boolean3.h
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,5 @@ class Boolean3 {
SparseIndices p1q2_, p2q1_;
Vec<int> x12_, x21_, w03_, w30_;
Vec<glm::vec3> v12_, v21_;
ExecutionPolicy policy_;
};
} // namespace manifold
Loading

0 comments on commit cc83e1b

Please sign in to comment.