diff --git a/algoim/organizer/loader.hpp b/algoim/organizer/loader.hpp index be67ade..1934613 100644 --- a/algoim/organizer/loader.hpp +++ b/algoim/organizer/loader.hpp @@ -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& 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; @@ -198,30 +179,30 @@ public: algoim::organizer::BlobTree tree; algoim::organizer::Blob blob0; - blob0.isPrimitive = 1; - blob0.nodeOp = 0; - blob0.inOut = 0; + blob0.isPrimitive = 1; + blob0.nodeOp = 0; + blob0.inOut = 0; blob0.oneChildInOut = 0; - blob0.isLeft = 1; - blob0.ancestor = 2; + blob0.isLeft = 1; + blob0.ancestor = 2; tree.structure.push_back(blob0); algoim::organizer::Blob blob1; - blob1.isPrimitive = 1; - blob1.nodeOp = 0; - blob1.inOut = 0; + blob1.isPrimitive = 1; + blob1.nodeOp = 0; + blob1.inOut = 0; blob1.oneChildInOut = 0; - blob1.isLeft = 0; - blob1.ancestor = 0; + blob1.isLeft = 0; + blob1.ancestor = 0; tree.structure.push_back(blob1); algoim::organizer::Blob blob2; - blob2.isPrimitive = 0; - blob2.nodeOp = 3; // no set - blob2.inOut = 0; + blob2.isPrimitive = 0; + blob2.nodeOp = 3; // no set + blob2.inOut = 0; blob2.oneChildInOut = 0; - blob2.isLeft = 0; - blob2.ancestor = 0; + blob2.isLeft = 0; + blob2.ancestor = 0; tree.structure.push_back(blob2); tree.primitiveNodeIdx.push_back(0); @@ -239,38 +220,28 @@ public: */ algoim::organizer::VisiblePrimitiveRep unionNode(const algoim::organizer::VisiblePrimitiveRep& rep1, const algoim::organizer::VisiblePrimitiveRep& rep2, - const bool modify = true) + const bool modify = true) { - auto tree = createEmptyBlobTree(); + auto tree = createEmptyBlobTree(); tree.structure[2].nodeOp = 0; const std::vector reps = {rep1, rep2}; - std::vector minimalReps; + std::vector minimalReps; algoim::organizer::mergeSubtree2Leaf(tree, minimalReps, reps); algoim::organizer::VisiblePrimitiveRep result; result.subBlobTree = tree; - result.aabb = rep1.aabb; + 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()); } @@ -287,38 +258,28 @@ public: */ algoim::organizer::VisiblePrimitiveRep intersectNode(const algoim::organizer::VisiblePrimitiveRep& rep1, const algoim::organizer::VisiblePrimitiveRep& rep2, - const bool modify = true) + const bool modify = true) { - auto tree = createEmptyBlobTree(); + auto tree = createEmptyBlobTree(); tree.structure[2].nodeOp = 1; const std::vector reps = {rep1, rep2}; - std::vector minimalReps; + std::vector minimalReps; algoim::organizer::mergeSubtree2Leaf(tree, minimalReps, reps); algoim::organizer::VisiblePrimitiveRep result; result.subBlobTree = tree; - result.aabb = rep1.aabb; + 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()); } @@ -335,37 +296,27 @@ public: */ algoim::organizer::VisiblePrimitiveRep differentNode(const algoim::organizer::VisiblePrimitiveRep& rep1, const algoim::organizer::VisiblePrimitiveRep& rep2, - const bool modify = true) + const bool modify = true) { - auto tree = createEmptyBlobTree(); + auto tree = createEmptyBlobTree(); tree.structure[2].nodeOp = 2; const std::vector reps = {rep1, rep2}; - std::vector minimalReps; + std::vector minimalReps; algoim::organizer::mergeSubtree2Leaf(tree, minimalReps, reps); algoim::organizer::VisiblePrimitiveRep result; result.subBlobTree = tree; - result.aabb = rep1.aabb; + 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()); } @@ -375,21 +326,16 @@ public: void unionNode(const BodyTag body1, const BodyTag body2) { - auto& rep1 = this->m_allVisible[body1]; - auto& rep2 = this->m_allVisible[body2]; + auto& rep1 = this->m_allVisible[body1]; + 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; @@ -397,17 +343,14 @@ public: void intersectNode(const BodyTag body1, const BodyTag body2) { - auto& rep1 = this->m_allVisible[body1]; - auto& rep2 = this->m_allVisible[body2]; + auto& rep1 = this->m_allVisible[body1]; + 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; @@ -415,16 +358,13 @@ public: void differentNode(const BodyTag body1, const BodyTag body2) { - auto& rep1 = this->m_allVisible[body1]; - auto& rep2 = this->m_allVisible[body2]; + auto& rep1 = this->m_allVisible[body1]; + 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); } @@ -434,23 +374,21 @@ public: void offset(const BodyTag body, const Direction3D& directrion, const double length) { algoim::uvector scale = 1; - algoim::uvector bias = -directrion.getUVector3Data(); + algoim::uvector bias = -directrion.getUVector3Data(); 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) { - auto& rep = this->m_allVisible[body]; - auto halfPlane = this->createHalfPlane(basePoint, -normal); - auto result = this->intersectNode(rep, halfPlane, false); + auto& rep = this->m_allVisible[body]; + auto halfPlane = this->createHalfPlane(basePoint, -normal); + auto result = this->intersectNode(rep, halfPlane, false); this->m_allVisible[body] = result; } @@ -461,29 +399,22 @@ public: * @return The polygonal column */ algoim::organizer::VisiblePrimitiveRep createPolygonalColumnWithoutTopBottom(const std::vector& points, - const Vector3D& extusion) + const Vector3D& extusion) { int pointNumber = points.size(); std::vector vertices; - std::vector indices; - std::vector indexInclusiveScan; + std::vector indices; + std::vector 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); @@ -493,17 +424,14 @@ public: indexInclusiveScan.push_back(index); } - algoim::organizer::MeshDesc polygonalColumn(vertices, indices, indexInclusiveScan); + algoim::organizer::MeshDesc polygonalColumn(vertices, indices, indexInclusiveScan); algoim::organizer::VisiblePrimitiveRep result; result.tensors.resize(pointNumber, algoim::tensor3(nullptr, 3)); std::vector*> temp; 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; } @@ -517,12 +445,12 @@ public: * @return The cylinder column */ algoim::organizer::VisiblePrimitiveRep createCylinderWithoutTopBottom(const Point3D& origion, - const double radius, - const double length, - const int alignAxis) + const double radius, + const double length, + const int alignAxis) { algoim::uvector3 ext = 3; - ext(alignAxis) = 1; + ext(alignAxis) = 1; algoim::organizer::VisiblePrimitiveRep result; result.tensors.resize(1, algoim::tensor3(nullptr, ext)); @@ -530,7 +458,7 @@ public: algoim::algoimSparkAllocHeapVector(resultTemp, result.tensors); this->m_allPointer.push_back(resultTemp[0]); - algoim::organizer::CylinderDesc cylinderDesc(origion.getUVector3Data(), radius, length, alignAxis); + algoim::organizer::CylinderDesc cylinderDesc(origion.getUVector3Data(), radius, length, alignAxis); algoim::organizer::VisiblePrimitiveRep cylinder; cylinder.tensors.resize(3, algoim::tensor3(nullptr, 3)); cylinder.tensors[0].ext_ = ext; @@ -538,7 +466,7 @@ public: algoim::organizer::makeCylinder(cylinderDesc, cylinder); result.tensors[0] = cylinder.tensors[0]; - result.aabb = cylinder.aabb; + result.aabb = cylinder.aabb; result.subBlobTree.primitiveNodeIdx.push_back(0); result.subBlobTree.structure.push_back(algoim::organizer::Blob{1, 0, 0, 0, 0, 0}); @@ -571,102 +499,92 @@ public: * @param[in] extusion The Stretch direction and length */ BodyTag addExtrudeWithTwoPoint(const std::vector& points, - const std::vector& bulges, - const Vector3D& extusion) + const std::vector& bulges, + const Vector3D& extusion) { - auto normal = Direction3D(extusion); + auto normal = Direction3D(extusion); auto& point1 = points[0]; auto& point2 = points[1]; - auto bulge1 = bulges[0]; - auto bulge2 = bulges[1]; + auto bulge1 = bulges[0]; + auto bulge2 = bulges[1]; algoim::organizer::VisiblePrimitiveRep rep1, rep2, result; auto halfDistance = point1.getDistance(point2) / 2.0; - auto middlePoint = point1.getMiddlePoint(point2); + auto middlePoint = point1.getMiddlePoint(point2); auto middleToOrigion1 = normal.cross(Direction3D(point2 - point1)); auto middleToOrigion2 = normal.cross(Direction3D(point1 - point2)); /* 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."); } auto getPrimitive = [this, normal, halfDistance, middlePoint, extusion, alignAxis]( const Point3D& point1, const Point3D& point2, const double bulge) { - auto middleToOrigion = normal.cross(Direction3D(point2 - point1)); - double sinHalfTheta = 2 * bulge / (1 + bulge * bulge); - double radius = halfDistance / sinHalfTheta; - double scalar = std::sqrt(radius * radius - halfDistance * halfDistance); + auto middleToOrigion = normal.cross(Direction3D(point2 - point1)); + 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}); - rep2 = getPrimitive(point2, point1, bulge2); + // 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); - rep1 = getPrimitive(point1, point2, bulge1); - rep2 = this->createHalfPlane(point2, -Direction3D{middleToOrigion1}); + // GJJ + rep1 = getPrimitive(point1, point2, bulge1); + 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); } } auto halfPlane1 = createHalfPlane(points[0], -normal); auto halfPlane2 = createHalfPlane(points[0] + extusion, normal); - result = this->intersectNode(result, halfPlane1); - result = this->intersectNode(result, halfPlane2); + result = this->intersectNode(result, halfPlane1); + result = this->intersectNode(result, halfPlane2); this->m_allVisible.push_back(result); return this->m_allVisible.size() - 1; @@ -682,69 +600,51 @@ 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]; + 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]; } /* Compute the origion and radius */ - auto halfDistance = point1.getDistance(point2) / 2.0; - auto middlePoint = point1.getMiddlePoint(point2); + auto halfDistance = point1.getDistance(point2) / 2.0; + auto middlePoint = point1.getMiddlePoint(point2); auto middleToOrigion = normal.cross(Direction3D(point2 - point1)); double sinHalfTheta = 2 * bulge / (1 + bulge * bulge); - double radius = halfDistance / sinHalfTheta; - double scalar = std::sqrt(radius * radius - halfDistance * halfDistance); + double radius = halfDistance / sinHalfTheta; + double scalar = std::sqrt(radius * radius - halfDistance * halfDistance); /* Determine whether to merge or subtract */ /* The operation is merge if flag is true, otherwise it is subtract */ bool flag; /* */ - auto centroidPoint = computePolygonCentroid(points); + 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,18 +667,15 @@ 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); + auto halfPlane = createHalfPlane(middlePoint, middleToOrigion); auto subtraction = this->intersectNode(cylinder, halfPlane); /* base + (cylinder - column) */ base = this->unionNode(base, subtraction); - } - else - { + } else { /* difference */ /* base - cylinder */ base = this->differentNode(base, cylinder); @@ -794,13 +684,10 @@ public: auto halfPlane1 = createHalfPlane(points[0], -normal); 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; - } + base = this->intersectNode(base, halfPlane1); + base = this->intersectNode(base, halfPlane2); + // 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*> 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; @@ -881,44 +752,30 @@ public: { auto normal = Direction3D{offset}; /* Determine which axis is aligned */ - int alignAxis; - if (normal.isParallel(Direction3D(1, 0, 0))) - { + int alignAxis; + 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."); } algoim::uvector3 ext = 3; - ext(alignAxis) = 1; + ext(alignAxis) = 1; algoim::organizer::VisiblePrimitiveRep cylinder; cylinder.tensors.resize(3, algoim::tensor3(nullptr, 3)); cylinder.tensors[0].ext_ = ext; std::vector*> 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; @@ -937,9 +794,8 @@ public: assert(rep.tensors.size() == rep.aabbs.size()); std::vector minimalReps; - algoim::uvector3 range = rep.aabb.max - rep.aabb.min; - for (size_t i = 0; i < rep.tensors.size(); i++) - { + algoim::uvector3 range = rep.aabb.max - rep.aabb.min; + for (size_t i = 0; i < rep.tensors.size(); i++) { auto& tenser = rep.tensors[i]; algoim::organizer::detail::powerTransformation(range, rep.aabb.min, tenser); @@ -949,31 +805,24 @@ public: minimalReps.push_back(minimalRep); } - auto size = rep.aabb.size(); - algoim::organizer::Scene scene{minimalReps, + auto size = rep.aabb.size(); + 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); - } - int cnt = 1; + for (int i = 0; i < minimalReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); } + int cnt = 1; std::vector leaves; algoim::organizer::buildOcTreeV0(scene, rootNode, leaves, 1, cnt, 0, -1); std::cout << "octree built over" << std::endl; // basicTask(scene, leaves[14], q); - int i = 0; - int q = 10; + int i = 0; + int q = 10; double volume = 0; - double area = 0; - for (const auto& leaf : leaves) - { + double area = 0; + 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,11 +834,9 @@ public: return std::make_pair(area, volume); } - void output(const BodyTag& tag) - { - } + void output(const BodyTag& tag) {} private: std::vector m_allVisible; - std::vector*> m_allPointer; + std::vector*> m_allPointer; }; diff --git a/examples/example_loader.cpp b/examples/example_loader.cpp index 2d07889..4f3e9dd 100644 --- a/examples/example_loader.cpp +++ b/examples/example_loader.cpp @@ -5,10 +5,10 @@ void loaderTest1() Loader loader; std::vector points; - std::vector bulges; - Vector3D extusion; - Point3D topPoint, bottomPoint, basePoint, leftBottomPoint; - Direction3D direction; + std::vector bulges; + Vector3D extusion; + Point3D topPoint, bottomPoint, basePoint, leftBottomPoint; + Direction3D direction; double offset, radius1, radius2, length, width, height, areaDifference, volumeDifference; @@ -21,10 +21,10 @@ void loaderTest1() points.push_back(Point3D{32449.999998877582, -3.8392495305561452e-07, 0.0000000000000000}); bulges.push_back(0.99999999999999989); bulges.push_back(0.99999999999999989); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag1 = loader.addExtrude(points, bulges, extusion); - loader.getAreaAndVolume(tag1); + loader.getAreaAndVolume(tag1); /* 体2 */ points.clear(); @@ -37,10 +37,10 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag2 = loader.addExtrude(points, bulges, extusion); - loader.getAreaAndVolume(tag2); + loader.getAreaAndVolume(tag2); /* 体3 */ points.clear(); @@ -53,31 +53,31 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag3 = loader.addExtrude(points, bulges, extusion); - loader.getAreaAndVolume(tag3); + loader.getAreaAndVolume(tag3); /* 体4 */ - topPoint = Point3D{-1.1224183253943920e-06, -3.8391322798592142e-07, 0.0000000000000000}; + topPoint = Point3D{-1.1224183253943920e-06, -3.8391322798592142e-07, 0.0000000000000000}; bottomPoint = Point3D{-1.1224183253943920e-06, -3.8391322798592142e-07, 3300.0000000000000000}; - radius1 = 32450.000000000000; - radius2 = 33539.000000000000; - auto tag4 = loader.addCone(topPoint, bottomPoint, radius1, radius2); + radius1 = 32450.000000000000; + radius2 = 33539.000000000000; + auto tag4 = loader.addCone(topPoint, bottomPoint, radius1, radius2); - loader.getAreaAndVolume(tag4); + loader.getAreaAndVolume(tag4); /* 体3被切割 */ basePoint = Point3D{-1.1224183253943920e-06, -3.8391322798592142e-07, 0.0000000000000000}; direction = Direction3D{0, 0, 1}; loader.split(tag3, basePoint, direction); - loader.getAreaAndVolume(tag3); + loader.getAreaAndVolume(tag3); /* 体3和体4布尔差 */ loader.differentNode(tag3, tag4); - loader.getAreaAndVolume(tag3); + loader.getAreaAndVolume(tag3); /* 体5 */ points.clear(); @@ -90,7 +90,7 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag5 = loader.addExtrude(points, bulges, extusion); /* 体5被切割 */ @@ -112,7 +112,7 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag6 = loader.addExtrude(points, bulges, extusion); /* 体6被切割 */ @@ -137,7 +137,7 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag7 = loader.addExtrude(points, bulges, extusion); /* 体8 */ @@ -151,15 +151,15 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag8 = loader.addExtrude(points, bulges, extusion); /* 体9 */ - topPoint = Point3D{-1.1224183253943920e-06, -3.8392403600706615e-07, 0.0000000000000000}; + topPoint = Point3D{-1.1224183253943920e-06, -3.8392403600706615e-07, 0.0000000000000000}; bottomPoint = Point3D{-1.1224183253943920e-06, -3.8392403600706615e-07, 3300.0000000000000000}; - radius1 = 32450.000000000000; - radius2 = 33539.000000000000; - auto tag9 = loader.addCone(topPoint, bottomPoint, radius1, radius2); + radius1 = 32450.000000000000; + radius2 = 33539.000000000000; + auto tag9 = loader.addCone(topPoint, bottomPoint, radius1, radius2); /* 体8被切割 */ basePoint = Point3D{-1.1224183253943920e-06, -3.8392403600706615e-07, 3300.0000000000000000}; @@ -180,7 +180,7 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag10 = loader.addExtrude(points, bulges, extusion); /* 体10被切割 */ @@ -202,7 +202,7 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(-0.99999999999999989); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 3300.0000000000000}; auto tag11 = loader.addExtrude(points, bulges, extusion); /* 体11被切割 */ @@ -224,7 +224,7 @@ void loaderTest1() /* 体1 Z轴偏移 -3600 */ direction = Direction3D{0, 0, 1}; - offset = -3600; + offset = -3600; loader.offset(tag1, direction, offset); /* 新的体11 */ @@ -234,7 +234,7 @@ void loaderTest1() points.push_back(Point3D{32049.999998877582, -3.8387224776670337e-07, 0.0000000000000000}); bulges.push_back(0.99999999999999989); bulges.push_back(0.99999999999999989); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; auto newTag11 = loader.addExtrude(points, bulges, extusion); /* 体0 */ @@ -244,7 +244,7 @@ void loaderTest1() points.push_back(Point3D{11801.337366082498, 5999.2409883221590, 0.0000000000000000}); bulges.push_back(0.99999999999999989); bulges.push_back(0.99999999999999989); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; auto tag0 = loader.addExtrude(points, bulges, extusion); /* 新的体11和体0布尔差 */ @@ -252,7 +252,7 @@ void loaderTest1() /* 新的体11 Z轴偏移 -3500 */ direction = Direction3D{0, 0, 1}; - offset = -3500; + offset = -3500; loader.offset(newTag11, direction, offset); /* 体1和新的体11判交 */ @@ -265,12 +265,12 @@ void loaderTest1() points.push_back(Point3D{-673.66711640879771, 5999.2409883221790, 0.0000000000000000}); bulges.push_back(0.99999999999999989); bulges.push_back(0.99999999999999989); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; auto tag12 = loader.addExtrude(points, bulges, extusion); /* 体12 Z轴偏移 -3500 */ direction = Direction3D{0, 0, 1}; - offset = -3500; + offset = -3500; loader.offset(tag12, direction, offset); /* 体1和体12判交 */ @@ -286,7 +286,7 @@ void loaderTest1() points.push_back(Point3D{2401.3334598319489, 5999.2409883221790, 0.0000000000000000}); bulges.push_back(0.99999999999999989); bulges.push_back(0.99999999999999989); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; auto tag13 = loader.addExtrude(points, bulges, extusion); /* 体14 */ @@ -296,7 +296,7 @@ void loaderTest1() points.push_back(Point3D{676.33403607269543, 5999.2409883221790, 0.0000000000000000}); bulges.push_back(-0.99999999999999989); bulges.push_back(-0.99999999999999989); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; auto tag14 = loader.addExtrude(points, bulges, extusion); /* 体13和体14布尔差 */ @@ -304,7 +304,7 @@ void loaderTest1() /* 体13 Z轴偏移-3500 */ direction = Direction3D{0, 0, 1}; - offset = -3500; + offset = -3500; loader.offset(tag13, direction, offset); /* 体1和体13判交 */ @@ -321,17 +321,17 @@ void loaderTest1() after = loader.getAreaAndVolume(tag1); /* 获取 体1 布尔前后的 体积差 和 面积差 */ - areaDifference = after.first - before.first; + areaDifference = after.first - before.first; volumeDifference = after.second - before.second; /* 定义空体 扣减体1 */ auto subTag1 = loader.addEmpty(); - auto cycle = [&loader, subTag1](const Point3D& point) { + auto cycle = [&loader, subTag1](const Point3D& point) { std::vector points; - std::vector bulges; - Vector3D extusion; - Direction3D direction; - double offset; + std::vector bulges; + Vector3D extusion; + Direction3D direction; + double offset; /* 循环体1 */ points.clear(); @@ -344,18 +344,18 @@ void loaderTest1() bulges.push_back(0.0000000000000000); bulges.push_back(0.0000000000000000); bulges.push_back(0.0000000000000000); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; auto cycleTag1 = loader.addExtrude(points, bulges, extusion); /* 循环体1 偏移+ 传入的X,Y,Z坐标值 */ loader.offset(cycleTag1, Direction3D{point}, point.length()); /* 循环体2 */ - auto leftBottomPoint = Point3D{point.m_x - 1000, point.m_y - 1000, point.m_z + 0.0000000000005}; - double length = 2000.0000000000036; - double width = 2000.0000000000036; - double height = 600; - auto cycleTag2 = loader.addBox(leftBottomPoint, length, width, height); + auto leftBottomPoint = Point3D{point.m_x - 1000, point.m_y - 1000, point.m_z + 0.0000000000005}; + double length = 2000.0000000000036; + double width = 2000.0000000000036; + double height = 600; + auto cycleTag2 = loader.addBox(leftBottomPoint, length, width, height); /* 循环体3 */ points.clear(); @@ -364,12 +364,12 @@ void loaderTest1() points.push_back(Point3D{32049.999998877582, -3.8387224776670337e-07, 0.0000000000000000}); bulges.push_back(0.99999999999999989); bulges.push_back(0.99999999999999989); - extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; + extusion = Vector3D{0.0000000000000000, 0.0000000000000000, 600.0000000000000}; auto cycleTag3 = loader.addExtrude(points, bulges, extusion); /* 循环体3 Z轴偏移 -2900 */ direction = Direction3D{0, 0, 1}; - offset = -2900; + offset = -2900; loader.offset(cycleTag3, direction, offset); /* 循环体3和体2布尔交 */ @@ -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); @@ -445,7 +442,7 @@ void loaderTest1() after = loader.getAreaAndVolume(tag1); /* 获取 体1 布尔前后的 体积差 和 面积差 */ - areaDifference = after.first - before.first; + areaDifference = after.first - before.first; volumeDifference = after.second - before.second; /* 体15：Ent1.bool */ @@ -459,7 +456,7 @@ void loaderTest1() points.push_back(Point3D{2501.2960250298465, 5999.2409883221790, -3500.0000000000000}); bulges.push_back(0.99999999999999989); bulges.push_back(0.99999999999999989); - extusion = Vector3D{-0.0000000000000000, -0.0000000000000000, -100.00000000000000}; + extusion = Vector3D{-0.0000000000000000, -0.0000000000000000, -100.00000000000000}; auto tag16 = loader.addExtrude(points, bulges, extusion); before = loader.getAreaAndVolume(tag1); @@ -470,7 +467,7 @@ void loaderTest1() after = loader.getAreaAndVolume(tag1); /* 获取 体1 布尔前后的 体积差 和 面积差 */ - areaDifference = after.first - before.first; + areaDifference = after.first - before.first; volumeDifference = after.second - before.second; } diff --git a/gjj/analyticSolution.py b/gjj/analyticSolution.py index 09d85ec..1087007 100644 --- a/gjj/analyticSolution.py +++ b/gjj/analyticSolution.py @@ -1,35 +1,52 @@ import math + class Sphere: - def __init__(self, radius, center): - self.radius = radius - self.center = center + def __init__(self, radius, center): + self.radius = radius + self.center = center # specify class for the input parameters def caseSpheresBOnAFace(rA, rB): - #intersection: - rAsq = rA**2 - rBsq = rB**2 - crownHeightA = rBsq / (2 * rA) - crownHeightB = rB - crownHeightA - crownVolumeA = math.pi * crownHeightA**2 * (3 * rA - crownHeightA) / 3 - crownVolumeB = math.pi * crownHeightB**2 * (3 * rB - crownHeightB) / 3 - volumeInter = crownVolumeA + crownVolumeB - - #union: - volumeA = 4/3 * math.pi * rAsq * rA - volumeB = 4/3 * math.pi * rBsq * rB - volumeUnion = volumeA + volumeB - volumeInter - - #difference: - volumeDiff = volumeA - volumeInter - - return volumeInter, volumeUnion, volumeDiff - - -if __name__=="__main__": - volumeInter, volumeUnion, volumeDiff = caseSpheresBOnAFace(800000, 1000) - print("volumeInter = ", volumeInter) - print("volumeUnion = ", volumeUnion) - print("volumeDiff = ", volumeDiff) \ No newline at end of file + # intersection: + rAsq = rA**2 + rBsq = rB**2 + crownHeightA = rBsq / (2 * rA) + crownHeightB = rB - crownHeightA + crownVolumeA = math.pi * crownHeightA**2 * (3 * rA - crownHeightA) / 3 + crownVolumeB = math.pi * crownHeightB**2 * (3 * rB - crownHeightB) / 3 + volumeInter = crownVolumeA + crownVolumeB + + # union: + volumeA = 4 / 3 * math.pi * rAsq * rA + volumeB = 4 / 3 * math.pi * rBsq * rB + volumeUnion = volumeA + volumeB - volumeInter + + # difference: + volumeDiff = volumeA - volumeInter + + 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) + print(caseExtrude()) diff --git a/gjj/test.py b/gjj/test.py new file mode 100644 index 0000000..e69de29