GJJ
/
HeteQuadrature
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

debug

master
gjj 8 months ago
parent
e41ddafec1
commit
6ffbea5b92
4 changed files with 296 additions and 435 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 387
      algoim/organizer/loader.hpp
  2. 5
      examples/example_loader.cpp
  3. 25
      gjj/analyticSolution.py
  4. 0
      gjj/test.py

387
algoim/organizer/loader.hpp
View File

@ -8,6 +8,7 @@ class Vector3D
{
public:
Vector3D() = default;
Vector3D(const double x, const double y, const double z)
{
this->m_x = x;
@ -15,17 +16,11 @@ public:
this->m_z = z;
}
double length() const
{
return std::sqrt(this->m_x * this->m_x + this->m_y * this->m_y + this->m_z * this->m_z);
}
double length() const { return std::sqrt(this->m_x * this->m_x + this->m_y * this->m_y + this->m_z * this->m_z); }
Vector3D operator/(const double scalar) const
{
if (scalar == 0)
{
throw std::runtime_error("Division by zero error.");
}
if (scalar == 0) { throw std::runtime_error("Division by zero error."); }
return Vector3D(this->m_x / scalar, this->m_y / scalar, this->m_z / scalar);
}
@ -37,15 +32,11 @@ public:
Vector3D cross(const Vector3D& other) const
{
return Vector3D(this->m_y * other.m_z - this->m_z * other.m_y,
this->m_z * other.m_x - this->m_x * other.m_z,
return Vector3D(this->m_y * other.m_z - this->m_z * other.m_y, this->m_z * other.m_x - this->m_x * other.m_z,
this->m_x * other.m_y - this->m_y * other.m_x);
}
double dot(const Vector3D& other) const
{
return this->m_x * other.m_x + this->m_y * other.m_y + this->m_z * other.m_z;
}
double dot(const Vector3D& other) const { return this->m_x * other.m_x + this->m_y * other.m_y + this->m_z * other.m_z; }
algoim::uvector3 getUVector3Data() const
{
@ -63,6 +54,7 @@ class Direction3D : public Vector3D
{
public:
Direction3D() = default;
Direction3D(const double x, const double y, const double z)
{
this->m_x = x;
@ -83,23 +75,16 @@ public:
Vector3D cross(const Direction3D& other) const
{
return Vector3D(this->m_y * other.m_z - this->m_z * other.m_y,
this->m_z * other.m_x - this->m_x * other.m_z,
return Vector3D(this->m_y * other.m_z - this->m_z * other.m_y, this->m_z * other.m_x - this->m_x * other.m_z,
this->m_x * other.m_y - this->m_y * other.m_x);
}
double dot(const Direction3D& other) const
{
return this->m_x * other.m_x + this->m_y * other.m_y + this->m_z * other.m_z;
}
double dot(const Direction3D& other) const { return this->m_x * other.m_x + this->m_y * other.m_y + this->m_z * other.m_z; }
void normalized()
{
double length = this->length();
if (std::abs(length) < 1e-8)
{
throw std::runtime_error("Cannot normalize a zero-length vector.");
}
if (std::abs(length) < 1e-8) { throw std::runtime_error("Cannot normalize a zero-length vector."); }
this->m_x /= length;
this->m_y /= length, this->m_z /= length;
@ -111,16 +96,14 @@ public:
return std::abs(cross.length()) < 1e-8;
}
Direction3D operator-() const
{
return Direction3D(-this->m_x, -this->m_y, -this->m_z);
}
Direction3D operator-() const { return Direction3D(-this->m_x, -this->m_y, -this->m_z); }
};
class Point3D : public Vector3D
{
public:
Point3D() = default;
Point3D(const double x, const double y, const double z)
{
this->m_x = x;
@ -155,9 +138,8 @@ public:
double getDistance(const Point3D& other) const
{
return std::sqrt((other.m_x - this->m_x) * (other.m_x - this->m_x) +
(other.m_y - this->m_y) * (other.m_y - this->m_y) +
(other.m_z - this->m_z) * (other.m_z - this->m_z));
return std::sqrt((other.m_x - this->m_x) * (other.m_x - this->m_x) + (other.m_y - this->m_y) * (other.m_y - this->m_y)
+ (other.m_z - this->m_z) * (other.m_z - this->m_z));
}
Point3D getMiddlePoint(const Point3D& other) const
@ -179,8 +161,7 @@ public:
Point3D computePolygonCentroid(const std::vector<Point3D>& points) const
{
double centroidX = 0, centroidY = 0, centroidZ = 0;
for (const auto& point : points)
{
for (const auto& point : points) {
centroidX += point.m_x;
centroidY += point.m_y;
centroidZ += point.m_z;
@ -253,24 +234,14 @@ public:
result.aabb = rep1.aabb;
result.aabb.extend(rep2.aabb);
if (modify)
{
for (auto& iter : minimalReps)
{
if (modify) {
for (auto& iter : minimalReps) {
result.tensors.push_back(iter.tensor);
result.aabbs.push_back(iter.aabb);
}
}
else
{
for (auto& iter : rep1.tensors)
{
result.tensors.push_back(iter);
}
for (auto& iter : rep2.tensors)
{
result.tensors.push_back(iter);
}
} else {
for (auto& iter : rep1.tensors) { result.tensors.push_back(iter); }
for (auto& iter : rep2.tensors) { result.tensors.push_back(iter); }
result.aabbs = rep1.aabbs;
result.aabbs.insert(result.aabbs.end(), rep2.aabbs.begin(), rep2.aabbs.end());
}
@ -301,24 +272,14 @@ public:
result.aabb = rep1.aabb;
result.aabb.intersect(rep2.aabb);
if (modify)
{
for (auto& iter : minimalReps)
{
if (modify) {
for (auto& iter : minimalReps) {
result.tensors.push_back(iter.tensor);
result.aabbs.push_back(iter.aabb);
}
}
else
{
for (auto& iter : rep1.tensors)
{
result.tensors.push_back(iter);
}
for (auto& iter : rep2.tensors)
{
result.tensors.push_back(iter);
}
} else {
for (auto& iter : rep1.tensors) { result.tensors.push_back(iter); }
for (auto& iter : rep2.tensors) { result.tensors.push_back(iter); }
result.aabbs = rep1.aabbs;
result.aabbs.insert(result.aabbs.end(), rep2.aabbs.begin(), rep2.aabbs.end());
}
@ -348,24 +309,14 @@ public:
result.subBlobTree = tree;
result.aabb = rep1.aabb;
if (modify)
{
for (auto& iter : minimalReps)
{
if (modify) {
for (auto& iter : minimalReps) {
result.tensors.push_back(iter.tensor);
result.aabbs.push_back(iter.aabb);
}
}
else
{
for (auto& iter : rep1.tensors)
{
result.tensors.push_back(iter);
}
for (auto& iter : rep2.tensors)
{
result.tensors.push_back(iter);
}
} else {
for (auto& iter : rep1.tensors) { result.tensors.push_back(iter); }
for (auto& iter : rep2.tensors) { result.tensors.push_back(iter); }
result.aabbs = rep1.aabbs;
result.aabbs.insert(result.aabbs.end(), rep2.aabbs.begin(), rep2.aabbs.end());
}
@ -379,17 +330,12 @@ public:
auto& rep2 = this->m_allVisible[body2];
algoim::organizer::VisiblePrimitiveRep result;
if (rep1.isEmpty())
{
if (rep1.isEmpty()) {
result = rep2;
}
else if (rep2.isEmpty())
{
} else if (rep2.isEmpty()) {
result = rep1;
}
else
{
result = this->differentNode(rep1, rep2, false);
} else {
result = this->unionNode(rep1, rep2, false);
}
this->m_allVisible[body1] = result;
@ -401,13 +347,10 @@ public:
auto& rep2 = this->m_allVisible[body2];
algoim::organizer::VisiblePrimitiveRep result;
if (rep1.isEmpty() || rep2.isEmpty())
{
if (rep1.isEmpty() || rep2.isEmpty()) {
result = algoim::organizer::VisiblePrimitiveRep{};
}
else
{
result = this->differentNode(rep1, rep2, false);
} else {
result = this->intersectNode(rep1, rep2, false);
}
this->m_allVisible[body1] = result;
@ -419,12 +362,9 @@ public:
auto& rep2 = this->m_allVisible[body2];
algoim::organizer::VisiblePrimitiveRep result;
if (rep1.isEmpty() || rep2.isEmpty())
{
if (rep1.isEmpty() || rep2.isEmpty()) {
result = rep1;
}
else
{
} else {
result = this->differentNode(rep1, rep2, false);
}
@ -438,12 +378,10 @@ public:
auto& rep = this->m_allVisible[body];
for (auto& iter : rep.tensors)
{
algoim::organizer::detail::powerTransformation(scale, bias, iter);
}
for (auto& iter : rep.tensors) { algoim::organizer::detail::powerTransformation(scale, bias, iter); }
rep.aabb += directrion.getUVector3Data();
for (auto& aabb : rep.aabbs) { aabb += directrion.getUVector3Data(); }
}
void split(const BodyTag body, const Point3D& basePoint, const Direction3D& normal)
@ -470,20 +408,13 @@ public:
std::vector<int> indexInclusiveScan;
/* All bottom point */
for (int i = 0; i < pointNumber; i++)
{
vertices.push_back(points[i].getUVector3Data());
}
for (int i = 0; i < pointNumber; i++) { vertices.push_back(points[i].getUVector3Data()); }
/* All top point */
for (int i = 0; i < pointNumber; i++)
{
vertices.push_back((points[i] + extusion).getUVector3Data());
}
for (int i = 0; i < pointNumber; i++) { vertices.push_back((points[i] + extusion).getUVector3Data()); }
/* Side face */
int index = 0;
for (int i = 0; i < pointNumber; i++)
{
for (int i = 0; i < pointNumber; i++) {
indices.push_back(i);
indices.push_back((i + 1) % pointNumber);
indices.push_back((i + 1) % pointNumber + pointNumber);
@ -500,10 +431,7 @@ public:
algoim::algoimSparkAllocHeapVector(temp, result.tensors);
algoim::organizer::makeMesh(polygonalColumn, result);
for (auto& pointer : temp)
{
this->m_allPointer.push_back(pointer);
}
for (auto& pointer : temp) { this->m_allPointer.push_back(pointer); }
return result;
}
@ -590,20 +518,13 @@ public:
/* Determine which axis is aligned */
int alignAxis;
if (normal.isParallel(Direction3D(1, 0, 0)))
{
if (normal.isParallel(Direction3D(1, 0, 0))) {
alignAxis = 0;
}
else if (normal.isParallel(Direction3D(0, 1, 0)))
{
} else if (normal.isParallel(Direction3D(0, 1, 0))) {
alignAxis = 1;
}
else if (normal.isParallel(Direction3D(0, 0, 1)))
{
} else if (normal.isParallel(Direction3D(0, 0, 1))) {
alignAxis = 2;
}
else
{
} else {
throw std::runtime_error("Non align axis cylinder.");
}
@ -613,52 +534,49 @@ public:
double sinHalfTheta = 2 * bulge / (1 + bulge * bulge);
double radius = halfDistance / sinHalfTheta;
double scalar = std::sqrt(radius * radius - halfDistance * halfDistance);
// GJJ
if (std::abs(bulge) > 1.0) { scalar *= -1; }
auto origion = middlePoint + middleToOrigion * scalar;
/* Create the cylinder face */
return this->createCylinderWithoutTopBottom(origion, radius, extusion.length(), alignAxis);
};
if (std::abs(bulge1) <= 1e-8)
{
if (std::abs(bulge1) <= 1e-8) {
assert(std::abs(bulge2) > 1e-8);
rep1 = this->createHalfPlane(point1, -Direction3D{middleToOrigion2});
// GJJ
rep1 = this->createHalfPlane(point1, Direction3D{middleToOrigion2});
// rep1 = this->createHalfPlane(point1, -Direction3D{middleToOrigion2});
rep2 = getPrimitive(point2, point1, bulge2);
result = this->intersectNode(rep1, rep2);
}
else if (std::abs(bulge2) <= 1e-8)
{
} else if (std::abs(bulge2) <= 1e-8) {
assert(std::abs(bulge1) > 1e-8);
// GJJ
rep1 = getPrimitive(point1, point2, bulge1);
rep2 = this->createHalfPlane(point2, -Direction3D{middleToOrigion1});
rep2 = this->createHalfPlane(point2, Direction3D{middleToOrigion1});
// rep2 = this->createHalfPlane(point2, -Direction3D{middleToOrigion1});
result = this->intersectNode(rep1, rep2);
}
else
{
} else {
rep1 = getPrimitive(point1, point2, bulge1);
rep2 = getPrimitive(point2, point1, bulge2);
/* if the bulge == 1 and bulge == 1, it is a cylinder */
if (std::abs(bulge1 - 1.0) < 1e-8 && std::abs(bulge2 - 1.0) < 1e-8)
{
if (std::abs(bulge1 - 1.0) < 1e-8 && std::abs(bulge2 - 1.0) < 1e-8) {
result = getPrimitive(point1, point2, bulge1);
}
/* if the bulge == -1 and bulge == -1, it is a cylinder */
else if (std::abs(bulge1 + 1.0) < 1e-8 && std::abs(bulge2 + 1.0) < 1e-8)
{
else if (std::abs(bulge1 + 1.0) < 1e-8 && std::abs(bulge2 + 1.0) < 1e-8) {
result = getPrimitive(point1, point2, bulge1);
}
/* if the bulge1 and bulge2 has the same symbol, it is merge */
else if (bulge1 * bulge2 > 0.0)
{
else if (bulge1 * bulge2 > 0.0) {
result = this->intersectNode(rep1, rep2);
}
else if (bulge1 > 0.0)
{
// GJJ
// } else if (bulge1 > 0.0) {
} else if (std::abs(bulge1) > std::abs(bulge2)) {
result = this->differentNode(rep1, rep2);
}
else
{
} else {
result = this->differentNode(rep2, rep1);
}
}
@ -682,32 +600,22 @@ public:
{
int pointNumber = points.size();
assert(pointNumber >= 2);
if (pointNumber == 2)
{
return addExtrudeWithTwoPoint(points, bulges, extusion);
}
if (pointNumber == 2) { return addExtrudeWithTwoPoint(points, bulges, extusion); }
/* Get base polygonal column */
auto base = createPolygonalColumnWithoutTopBottom(points, extusion);
auto normal = Direction3D(extusion);
for (int i = 0; i < points.size(); i++)
{
for (int i = 0; i < points.size(); i++) {
/* Get point and bulge data */
auto bulge = bulges[i];
if (std::abs(bulge) < 1e-8)
{
continue;
}
if (std::abs(bulge) < 1e-8) { continue; }
auto& point1 = points[i];
Point3D point2;
if (i + 1 == points.size())
{
if (i + 1 == points.size()) {
point2 = points[0];
}
else
{
} else {
point2 = points[i + 1];
}
@ -727,24 +635,16 @@ public:
auto centroidPoint = computePolygonCentroid(points);
auto middleToCentroid = Direction3D(centroidPoint - middlePoint);
/* out */
if (middleToCentroid.dot(middleToOrigion) > 0.0)
{
if (middleToCentroid.dot(middleToOrigion) > 0.0) {
/* |bulge| > 1 */
if (std::abs(bulge) > 1.0 + 1e-8)
{
scalar *= -1;
}
if (std::abs(bulge) > 1.0 + 1e-8) { scalar *= -1; }
flag = true;
}
/* in */
else
{
else {
/* |bulge| > 1 */
if (std::abs(bulge) > 1.0 + 1e-8)
{
scalar *= -1;
}
if (std::abs(bulge) > 1.0 + 1e-8) { scalar *= -1; }
flag = false;
}
@ -753,20 +653,13 @@ public:
/* Determine which axis is aligned */
int alignAxis;
if (normal.isParallel(Direction3D(1, 0, 0)))
{
if (normal.isParallel(Direction3D(1, 0, 0))) {
alignAxis = 0;
}
else if (normal.isParallel(Direction3D(0, 1, 0)))
{
} else if (normal.isParallel(Direction3D(0, 1, 0))) {
alignAxis = 1;
}
else if (normal.isParallel(Direction3D(0, 0, 1)))
{
} else if (normal.isParallel(Direction3D(0, 0, 1))) {
alignAxis = 2;
}
else
{
} else {
throw std::runtime_error("Non align axis cylinder.");
}
@ -774,8 +667,7 @@ public:
auto cylinder = createCylinderWithoutTopBottom(origion, radius, extusion.length(), alignAxis);
/* Perform union and difference operations on the basic prismatic faces */
if (flag)
{
if (flag) {
/* Union */
/* cylinder - half plane */
auto halfPlane = createHalfPlane(middlePoint, middleToOrigion);
@ -783,9 +675,7 @@ public:
/* base + (cylinder - column) */
base = this->unionNode(base, subtraction);
}
else
{
} else {
/* difference */
/* base - cylinder */
base = this->differentNode(base, cylinder);
@ -796,11 +686,8 @@ public:
auto halfPlane2 = createHalfPlane(points[0] + extusion, normal);
base = this->intersectNode(base, halfPlane1);
base = this->intersectNode(base, halfPlane2);
for (size_t i = 0; i < base.aabbs.size(); i++)
{
base.aabbs[i] = base.aabb;
}
// gjj
for (size_t i = 0; i < base.aabbs.size(); i++) { base.aabbs[i] = base.aabb; }
this->m_allVisible.push_back(base);
@ -809,44 +696,34 @@ public:
BodyTag addCone(const Point3D& topPoint, const Point3D& bottomPoint, const double radius1, const double radius2)
{
auto bottomToTop = topPoint - bottomPoint;
auto normal = Direction3D{bottomToTop};
// gjj
// auto bottomToTop = topPoint - bottomPoint;
// auto normal = Direction3D{bottomToTop};
/* Determine which axis is aligned */
int alignAxis;
if (normal.isParallel(Direction3D(1, 0, 0)))
{
alignAxis = 0;
}
else if (normal.isParallel(Direction3D(0, 1, 0)))
{
alignAxis = 1;
}
else if (normal.isParallel(Direction3D(0, 0, 1)))
{
alignAxis = 2;
}
else
{
throw std::runtime_error("Non align axis cone.");
}
auto coneDesc = algoim::organizer::ConeDesc{bottomPoint.getUVector3Data(), radius1, radius2, alignAxis};
// int alignAxis;
// if (normal.isParallel(Direction3D(1, 0, 0))) {
// alignAxis = 0;
// } else if (normal.isParallel(Direction3D(0, 1, 0))) {
// alignAxis = 1;
// } else if (normal.isParallel(Direction3D(0, 0, 1))) {
// alignAxis = 2;
// } else {
// throw std::runtime_error("Non align axis cone.");
// }
// auto coneDesc = algoim::organizer::ConeDesc{bottomPoint.getUVector3Data(), radius1, radius2, alignAxis};
auto coneDesc =
algoim::organizer::ConeDesc{bottomPoint.getUVector3Data(), topPoint.getUVector3Data(), radius2, radius1};
algoim::organizer::VisiblePrimitiveRep cone;
cone.tensors.resize(3, algoim::tensor3(nullptr, 3));
std::vector<algoim::SparkStack<algoim::real>*> temp;
algoim::algoimSparkAllocHeapVector(temp, cone.tensors);
algoim::organizer::makeCone(coneDesc, cone);
for (auto& pointer : temp)
{
this->m_allPointer.push_back(pointer);
}
for (auto& pointer : temp) { this->m_allPointer.push_back(pointer); }
for (int i = 0; i < 3; i++)
{
cone.aabbs.push_back(cone.aabb);
}
for (int i = 0; i < 3; i++) { cone.aabbs.push_back(cone.aabb); }
this->m_allVisible.push_back(cone);
return this->m_allVisible.size() - 1;
@ -863,15 +740,9 @@ public:
algoim::algoimSparkAllocHeapVector(temp, box.tensors);
algoim::organizer::makeMesh(boxDesc, box);
for (auto& pointer : temp)
{
this->m_allPointer.push_back(pointer);
}
for (auto& pointer : temp) { this->m_allPointer.push_back(pointer); }
for (int i = 0; i < 6; i++)
{
box.aabbs.push_back(box.aabb);
}
for (int i = 0; i < 6; i++) { box.aabbs.push_back(box.aabb); }
this->m_allVisible.push_back(box);
return this->m_allVisible.size() - 1;
@ -882,20 +753,13 @@ public:
auto normal = Direction3D{offset};
/* Determine which axis is aligned */
int alignAxis;
if (normal.isParallel(Direction3D(1, 0, 0)))
{
if (normal.isParallel(Direction3D(1, 0, 0))) {
alignAxis = 0;
}
else if (normal.isParallel(Direction3D(0, 1, 0)))
{
} else if (normal.isParallel(Direction3D(0, 1, 0))) {
alignAxis = 1;
}
else if (normal.isParallel(Direction3D(0, 0, 1)))
{
} else if (normal.isParallel(Direction3D(0, 0, 1))) {
alignAxis = 2;
}
else
{
} else {
throw std::runtime_error("Non align axis cylinder.");
}
@ -907,18 +771,11 @@ public:
cylinder.tensors[0].ext_ = ext;
std::vector<algoim::SparkStack<algoim::real>*> temp;
algoim::algoimSparkAllocHeapVector(temp, cylinder.tensors);
algoim::organizer::CylinderDesc cylinderDesc(
bottomOrigion.getUVector3Data(), radius, offset.length(), alignAxis);
algoim::organizer::CylinderDesc cylinderDesc(bottomOrigion.getUVector3Data(), radius, offset.length(), alignAxis);
for (auto& pointer : temp)
{
this->m_allPointer.push_back(pointer);
}
for (auto& pointer : temp) { this->m_allPointer.push_back(pointer); }
for (int i = 0; i < 3; i++)
{
cylinder.aabbs.push_back(cylinder.aabb);
}
for (int i = 0; i < 3; i++) { cylinder.aabbs.push_back(cylinder.aabb); }
this->m_allVisible.push_back(cylinder);
return this->m_allVisible.size() - 1;
@ -938,8 +795,7 @@ public:
std::vector<algoim::organizer::MinimalPrimitiveRep> minimalReps;
algoim::uvector3 range = rep.aabb.max - rep.aabb.min;
for (size_t i = 0; i < rep.tensors.size(); i++)
{
for (size_t i = 0; i < rep.tensors.size(); i++) {
auto& tenser = rep.tensors[i];
algoim::organizer::detail::powerTransformation(range, rep.aabb.min, tenser);
@ -953,10 +809,7 @@ public:
algoim::organizer::Scene scene{minimalReps,
algoim::organizer::AABB(rep.aabb.min - size * 1e-5, rep.aabb.max + size * 1e-5)};
algoim::organizer::OcTreeNode rootNode(0, 1, rep.subBlobTree);
for (int i = 0; i < minimalReps.size(); ++i)
{
rootNode.polyIntersectIndices.emplace_back(i);
}
for (int i = 0; i < minimalReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); }
int cnt = 1;
std::vector<algoim::organizer::OcTreeNode> leaves;
algoim::organizer::buildOcTreeV0(scene, rootNode, leaves, 1, cnt, 0, -1);
@ -967,13 +820,9 @@ public:
int q = 10;
double volume = 0;
double area = 0;
for (const auto& leaf : leaves)
{
for (const auto& leaf : leaves) {
auto basicRes = algoim::organizer::basicTask(scene, leaf, q);
if (std::isinf(basicRes.volume))
{
std::cout << "inf volume when solving leaf: " << i << std::endl;
}
if (std::isinf(basicRes.volume)) { std::cout << "inf volume when solving leaf: " << i << std::endl; }
volume += basicRes.volume;
std::cout << "Solved leaves: " << ++i << "/" << leaves.size()
<< ", primitive cnt: " << leaf.polyIntersectIndices.size() << ", volume: " << volume << std::endl;
@ -985,9 +834,7 @@ public:
return std::make_pair(area, volume);
}
void output(const BodyTag& tag)
{
}
void output(const BodyTag& tag) {}
private:
std::vector<algoim::organizer::VisiblePrimitiveRep> m_allVisible;

5
examples/example_loader.cpp
View File

@ -432,10 +432,7 @@ void loaderTest1()
points.push_back(Point3D{23121.637452736089, 11781.092611367474, -2900.0000000000000});
points.push_back(Point3D{11781.053468361089, -23121.619302695028, -2900.0000000000000});
points.push_back(Point3D{-23121.583250388911, 11781.089681679974, -2900.0000000000000});
for (auto& point : points)
{
cycle(point);
}
for (auto& point : points) { cycle(point); }
before = loader.getAreaAndVolume(tag1);

25
gjj/analyticSolution.py
View File

@ -1,5 +1,6 @@
import math
class Sphere:
def __init__(self, radius, center):
self.radius = radius
@ -28,8 +29,24 @@ def caseSpheresBOnAFace(rA, rB):
return volumeInter, volumeUnion, volumeDiff
def caseExtrude():
x1 = -32450.000001122418
y1 = -3.8391231093737305e07
x2 = 32449.999998877582
y2 = -3.8392495305561452e07
x3 = 33538.999998877582
y3 = 3.8392534654100421e07
x4 = -33539.000001122418
y4 = 3.8391228016184551e07
area = 0.5 * abs(
x1 * y2 + x2 * y3 + x3 * y4 + x4 * y1 - (y1 * x2 + y2 * x3 + y3 * x4 + y4 * x1)
)
return area * 3300
if __name__ == "__main__":
volumeInter, volumeUnion, volumeDiff = caseSpheresBOnAFace(800000, 1000)
print("volumeInter = ", volumeInter)
print("volumeUnion = ", volumeUnion)
print("volumeDiff = ", volumeDiff)
# volumeInter, volumeUnion, volumeDiff = caseSpheresBOnAFace(800000, 1000)
# print("volumeInter = ", volumeInter)
# print("volumeUnion = ", volumeUnion)
# print("volumeDiff = ", volumeDiff)
print(caseExtrude())

0
gjj/test.py
View File

Write Preview
Loading…
Cancel
Save