Skip to content

Commit 6d06b51

Browse files
prashjhaprashjha
committed
add ellipse geometry. It requires six parameters, rx, ry, x, y, z, theta
1 parent 68eb2cf commit 6d06b51

File tree

5 files changed

+376
-3
lines changed

5 files changed

+376
-3
lines changed

‎src/model/dem/demModel.cpp‎

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1463,7 +1463,24 @@ double model::DEMModel::createGeometryAtSite(const double &particle_radius,
14631463

14641464
scale = params[0] / rep_geom_params[0];
14651465
}
1466-
} else {
1466+
} else if (rep_geom_p->d_name == "ellipse") {
1467+
1468+
// case - objects requiring six parameters
1469+
size_t num_params = 6;
1470+
1471+
if (params.size() < num_params)
1472+
params.resize(num_params);
1473+
1474+
params[0] = particle_radius;
1475+
params[1] = particle_radius * rep_geom_params[1] / rep_geom_params[0];
1476+
for (int dof = 0; dof < 3; dof++)
1477+
params[dof + 2] = site[dof];
1478+
1479+
params[5] = particle_orient;
1480+
1481+
scale = params[0] / rep_geom_params[0];
1482+
}
1483+
else {
14671484
std::cerr << fmt::format("Error: PeriDEM supports following type "
14681485
"of geometries for particles = {}\n",
14691486
util::io::printStr(util::geometry::acceptable_geometries));

‎src/util/geomObjects.cpp‎

Lines changed: 137 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2045,4 +2045,140 @@ namespace util {
20452045
return oss.str();
20462046
}
20472047

2048-
}// ComplexGeomObject
2048+
}// ComplexGeomObject
2049+
2050+
//
2051+
// Ellipse
2052+
//
2053+
namespace util {
2054+
double util::geometry::Ellipse::volume() const {
2055+
return M_PI * d_a * d_b; // Area of ellipse = π * a * b
2056+
}
2057+
2058+
util::Point util::geometry::Ellipse::center() const {
2059+
return d_x;
2060+
}
2061+
2062+
std::pair<util::Point, util::Point> util::geometry::Ellipse::box() const {
2063+
return box(0.);
2064+
}
2065+
2066+
std::pair<util::Point, util::Point> util::geometry::Ellipse::box(const double &tol) const {
2067+
// For a rotated ellipse, we need to find the extreme points
2068+
double cos_t = std::cos(d_theta);
2069+
double sin_t = std::sin(d_theta);
2070+
2071+
// Maximum extents in x and y directions
2072+
double xmax = std::sqrt(std::pow(d_a * cos_t, 2) + std::pow(d_b * sin_t, 2));
2073+
double ymax = std::sqrt(std::pow(d_a * sin_t, 2) + std::pow(d_b * cos_t, 2));
2074+
2075+
return {d_x - util::Point(xmax + tol, ymax + tol, 0.),
2076+
d_x + util::Point(xmax + tol, ymax + tol, 0.)};
2077+
}
2078+
2079+
double util::geometry::Ellipse::inscribedRadius() const {
2080+
return std::min(d_a, d_b); // Minimum of semi-major and semi-minor axes
2081+
}
2082+
2083+
double util::geometry::Ellipse::boundingRadius() const {
2084+
return d_r; // Maximum of semi-major and semi-minor axes (set in constructor)
2085+
}
2086+
2087+
bool util::geometry::Ellipse::isInside(const util::Point &x) const {
2088+
// Transform point to ellipse coordinate system
2089+
util::Point dx = x - d_x;
2090+
double cos_t = std::cos(d_theta);
2091+
double sin_t = std::sin(d_theta);
2092+
2093+
// Rotate point to align with ellipse axes
2094+
double x_rot = dx.d_x * cos_t + dx.d_y * sin_t;
2095+
double y_rot = -dx.d_x * sin_t + dx.d_y * cos_t;
2096+
2097+
// Check standard ellipse equation (x²/a² + y²/b² ≤ 1)
2098+
return util::isLess(std::pow(x_rot/d_a, 2) + std::pow(y_rot/d_b, 2), 1.0);
2099+
}
2100+
2101+
bool util::geometry::Ellipse::isOutside(const util::Point &x) const {
2102+
return !isInside(x);
2103+
}
2104+
2105+
bool util::geometry::Ellipse::isNear(const util::Point &x, const double &tol) const {
2106+
// Transform point to ellipse coordinate system
2107+
util::Point dx = x - d_x;
2108+
double cos_t = std::cos(d_theta);
2109+
double sin_t = std::sin(d_theta);
2110+
2111+
// Rotate point to align with ellipse axes
2112+
double x_rot = dx.d_x * cos_t + dx.d_y * sin_t;
2113+
double y_rot = -dx.d_x * sin_t + dx.d_y * cos_t;
2114+
2115+
// Calculate normalized distance from ellipse boundary
2116+
double dist = std::pow(x_rot/d_a, 2) + std::pow(y_rot/d_b, 2);
2117+
return util::isLess(std::abs(dist - 1.0), tol/std::min(d_a, d_b));
2118+
}
2119+
2120+
bool util::geometry::Ellipse::isNearBoundary(const util::Point &x,
2121+
const double &tol,
2122+
const bool &within) const {
2123+
return isNear(x, tol) && (within ? isInside(x) : isOutside(x));
2124+
}
2125+
2126+
bool util::geometry::Ellipse::doesIntersect(const util::Point &x) const {
2127+
return isNearBoundary(x, 1.0E-8, false);
2128+
}
2129+
2130+
bool util::geometry::Ellipse::isInside(
2131+
const std::pair<util::Point, util::Point> &box) const {
2132+
// Check if all corners of the box are inside the ellipse
2133+
for (auto p: util::getCornerPoints(2, box))
2134+
if (!isInside(p))
2135+
return false;
2136+
return true;
2137+
}
2138+
2139+
bool util::geometry::Ellipse::isOutside(
2140+
const std::pair<util::Point, util::Point> &box) const {
2141+
// Check if all corners of the box are outside the ellipse
2142+
for (auto p: util::getCornerPoints(2, box))
2143+
if (!isOutside(p))
2144+
return false;
2145+
return true;
2146+
}
2147+
2148+
bool util::geometry::Ellipse::isNear(
2149+
const std::pair<util::Point, util::Point> &box,
2150+
const double &tol) const {
2151+
return util::areBoxesNear(this->box(), box, tol, 2);
2152+
}
2153+
2154+
bool util::geometry::Ellipse::doesIntersect(
2155+
const std::pair<util::Point, util::Point> &box) const {
2156+
// Check if any corner point is inside
2157+
for (auto p: util::getCornerPoints(2, box))
2158+
if (isInside(p))
2159+
return true;
2160+
return false;
2161+
}
2162+
2163+
std::string util::geometry::Ellipse::printStr(int nt, int lvl) const {
2164+
auto tabS = util::io::getTabS(nt);
2165+
std::ostringstream oss;
2166+
2167+
oss << tabS << "------- Ellipse --------" << std::endl << std::endl;
2168+
oss << tabS << "Name = " << d_name << std::endl;
2169+
oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
2170+
oss << tabS << "Semi-major axis = " << d_a << std::endl;
2171+
oss << tabS << "Semi-minor axis = " << d_b << std::endl;
2172+
oss << tabS << "Rotation angle = " << d_theta << " radians" << std::endl;
2173+
oss << std::endl;
2174+
2175+
if (lvl > 0)
2176+
oss << tabS << "Bounding box: "
2177+
<< util::io::printBoxStr(box(0.), nt + 1);
2178+
2179+
if (lvl == 0)
2180+
oss << std::endl;
2181+
2182+
return oss.str();
2183+
}
2184+
} // Ellipse

‎src/util/geomObjects.h‎

Lines changed: 170 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2053,6 +2053,176 @@ namespace util {
20532053
double d_r;
20542054
};
20552055

2056+
/*!
2057+
* @brief Defines ellipse
2058+
*/
2059+
class Ellipse : public GeomObject {
2060+
2061+
public:
2062+
/*!
2063+
* @brief Constructor
2064+
*/
2065+
Ellipse()
2066+
: GeomObject("ellipse", ""),
2067+
d_x(util::Point()),
2068+
d_a(0.),
2069+
d_b(0.),
2070+
d_theta(0.),
2071+
d_r(0.) {};
2072+
2073+
/*!
2074+
* @brief Constructor
2075+
*
2076+
* @param a Semi-major axis length
2077+
* @param b Semi-minor axis length
2078+
* @param x Center point
2079+
* @param theta Rotation angle in radians (counterclockwise from x-axis)
2080+
* @param description Description of object (e.g., further classification or any tag)
2081+
*/
2082+
Ellipse(double a, double b, util::Point x = util::Point(),
2083+
double theta = 0., std::string description = "")
2084+
: GeomObject("ellipse", description),
2085+
d_x(x),
2086+
d_a(a),
2087+
d_b(b),
2088+
d_theta(theta),
2089+
d_r(std::max(a, b)) {};
2090+
2091+
/*!
2092+
* @copydoc GeomObject::volume() const
2093+
*/
2094+
double volume() const override;
2095+
2096+
/*!
2097+
* @copydoc GeomObject::center() const
2098+
*/
2099+
util::Point center() const override;
2100+
2101+
/*!
2102+
* @copydoc GeomObject::box() const
2103+
*/
2104+
std::pair<util::Point, util::Point> box() const override;
2105+
2106+
/*!
2107+
* @copydoc GeomObject::box(const double &tol) const
2108+
*/
2109+
std::pair<util::Point, util::Point>
2110+
box(const double &tol) const override;
2111+
2112+
/*!
2113+
* @copydoc GeomObject::inscribedRadius() const
2114+
*/
2115+
double inscribedRadius() const override;
2116+
2117+
/*!
2118+
* @copydoc GeomObject::boundingRadius() const
2119+
*/
2120+
double boundingRadius() const override;
2121+
2122+
/**
2123+
* @name Interaction with point
2124+
*/
2125+
/** @{*/
2126+
2127+
/*!
2128+
* @copydoc GeomObject::isInside(const util::Point &x) const
2129+
*/
2130+
bool isInside(const util::Point &x) const override;
2131+
2132+
/*!
2133+
* @copydoc GeomObject::isOutside(const util::Point &x) const
2134+
*/
2135+
bool isOutside(const util::Point &x) const override;
2136+
2137+
/*!
2138+
* @copydoc GeomObject::isNear(const util::Point &x, const double &tol) const
2139+
*/
2140+
bool isNear(const util::Point &x, const double &tol) const override;
2141+
2142+
/*!
2143+
* @copydoc GeomObject::isNearBoundary(const util::Point &x, const double &tol,
2144+
const bool &within) cons
2145+
*/
2146+
bool isNearBoundary(const util::Point &x, const double &tol,
2147+
const bool &within) const override;
2148+
2149+
/*!
2150+
* @copydoc GeomObject::doesIntersect(const util::Point &x) const
2151+
*/
2152+
bool doesIntersect(const util::Point &x) const override;
2153+
2154+
/** @}*/
2155+
2156+
/**
2157+
* @name Interaction with box
2158+
*/
2159+
/** @{*/
2160+
2161+
/*!
2162+
* @copydoc GeomObject::isInside(
2163+
const std::pair<util::Point, util::Point> &box) const
2164+
*/
2165+
bool
2166+
isInside(
2167+
const std::pair<util::Point, util::Point> &box) const override;
2168+
2169+
/*!
2170+
* @copydoc GeomObject::isOutside(
2171+
const std::pair<util::Point, util::Point> &box) const
2172+
*/
2173+
bool
2174+
isOutside(
2175+
const std::pair<util::Point, util::Point> &box) const override;
2176+
2177+
/*!
2178+
* @copydoc GeomObject::isNear(const std::pair<util::Point, util::Point> &box,
2179+
const double &tol) const
2180+
*/
2181+
bool isNear(const std::pair<util::Point, util::Point> &box,
2182+
const double &tol) const override;
2183+
2184+
/*!
2185+
* @copydoc GeomObject::doesIntersect(
2186+
const std::pair<util::Point, util::Point> &box) const
2187+
*/
2188+
bool doesIntersect(
2189+
const std::pair<util::Point, util::Point> &box) const override;
2190+
2191+
/** @}*/
2192+
2193+
/*!
2194+
* @copydoc GeomObject::printStr(int nt, int lvl) const
2195+
*/
2196+
std::string printStr(int nt, int lvl) const override;
2197+
2198+
/*!
2199+
* @copydoc GeomObject::print(int nt, int lvl) const
2200+
*/
2201+
void print(int nt, int lvl) const override {
2202+
std::cout << printStr(nt, lvl);
2203+
};
2204+
2205+
/*!
2206+
* @copydoc GeomObject::print() const
2207+
*/
2208+
void print() const override { print(0, 0); };
2209+
2210+
/*! @brief Center */
2211+
util::Point d_x;
2212+
2213+
/*! @brief Semi-major axis length */
2214+
double d_a;
2215+
2216+
/*! @brief Semi-minor axis length */
2217+
double d_b;
2218+
2219+
/*! @brief Rotation angle in radians (counterclockwise from x-axis) */
2220+
double d_theta;
2221+
2222+
/*! @brief Bounding radius (max of semi-major and semi-minor axis) */
2223+
double d_r;
2224+
};
2225+
20562226
/*!
20572227
* @brief Defines sphere
20582228
*/

0 commit comments

Comments
 (0)