commit 18160a766ef033d06df92467857fc86235dcdbf4
parent 668eedbee627f6d490e4345d0f85573f56f37a05
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 22 Mar 2016 14:29:59 +0100
Update the ellipsoid distribution
The "a" and "b" semi-principal axis are now equals. Add the ellipsoid
distribution with an equivalent sphere. The radius of the equivalent
sphere follows a given distribution and the a/c ratio is fixed.
Diffstat:
4 files changed, 107 insertions(+), 62 deletions(-)
diff --git a/doc/schiff-geometry.5 b/doc/schiff-geometry.5
@@ -28,10 +28,10 @@ in triangular meshes with respect to the \fBslices\fR attribute, i.e. the number
of discrete divisions along 2PI. By default \fBslices\fR is set to 64.
.PP
An ellipsoidal distribution of the soft particle shapes is controlled by the
-distribution of the 3 semi\-principal axis \fBa\fR, \fBb\fR and \fBc\fR of the
+distribution of the 2 semi\-principal axis \fBa\fR and \fBc\fR of the
\fBellipsoid\fR equation:
.IP " " 4
-(x/\fBa\fR)^2 + (y/\fBb\fR)^2 + (z/\fBc\fR)^2 = 1
+(x/\fBa\fR)^2 + (y/\fBa\fR)^2 + (z/\fBc\fR)^2 = 1
.PP
A cylindrical distribution can be defined in two ways. The first manner is to
directly control the distribution of the \fBheight\fR and the \fBradius\fR of
@@ -83,7 +83,6 @@ the example section for illustrations of such alternatives.
.TP
<\fIellipsoid\-geometry\fR>
\fBa:\fR <\fIdistribution\fR>
- \fBb:\fR <\fIdistribution\fR>
\fBc:\fR <\fIdistribution\fR>
.TP
<\fIcylinder\-geometry\fR> ::=
@@ -140,7 +139,6 @@ Soft particles are ellipsoids whose one of its semi-principal axis is
distributed with respect to a lognormal distribution:
\fBellipsoid:
a: 1.0
- b: 2.1
c:
lognormal:
sigma: 0.2
diff --git a/src/schiff_args.c b/src/schiff_args.c
@@ -105,8 +105,10 @@ geometry_release(struct schiff_geometry* geom)
switch(geom->type) {
case SCHIFF_ELLIPSOID:
param_release(&geom->data.ellipsoid.a);
- param_release(&geom->data.ellipsoid.b);
param_release(&geom->data.ellipsoid.c);
+ case SCHIFF_ELLIPSOID_AS_SPHERE:
+ param_release(&geom->data.ellipsoid.radius);
+ break;
case SCHIFF_CYLINDER:
param_release(&geom->data.cylinder.radius);
param_release(&geom->data.cylinder.height);
@@ -532,9 +534,10 @@ parse_yaml_ellipsoid
enum {
PROBA = BIT(0),
LENGTH_a = BIT(1),
- LENGTH_b = BIT(2),
- LENGTH_c = BIT(3),
- SLICES = BIT(4)
+ LENGTH_c = BIT(2),
+ RADIUS = BIT(3),
+ ASPECT_RATIO = BIT(4),
+ SLICES = BIT(5)
};
int mask = 0; /* Register the parsed histogram attributes */
size_t nattrs;
@@ -558,74 +561,97 @@ parse_yaml_ellipsoid
val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value);
ASSERT(key->type == YAML_SCALAR_NODE);
+ #define SETUP_MASK(Flag, Name) { \
+ if(mask & Flag) { \
+ log_err(filename, node, "the "Name" is already defined.\n"); \
+ return RES_BAD_ARG; \
+ } \
+ mask |= Flag; \
+ } (void)0
+
+ /* Probability of the distribution */
if(!strcmp((char*)key->data.scalar.value, "proba")) {
- if(mask & PROBA) {
- log_err(filename, node, "the sphere proba is already defined.\n");
- return RES_BAD_ARG;
- }
- mask |= PROBA;
+ SETUP_MASK(PROBA, "sphere proba");
res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, geom_proba);
/* # slices used to discretized the triangular ellipsoid */
} else if(!strcmp((char*)key->data.scalar.value, "slices")) {
- if(mask & SLICES) {
- log_err(filename, key,
- "the ellipsoid number of slices is already defined.\n");
- return RES_BAD_ARG;
- }
- mask |= SLICES;
+ SETUP_MASK(SLICES, "ellipsoid number of slices");
res = parse_yaml_uint
(filename, val, 4, 32768, &geom->data.ellipsoid.nslices);
- }
- /* Length of the ellipsoid semi-principal axis */
- #define PARSE_SEMI_PRINCIPAL_AXIS(Axis) { \
- if(mask & LENGTH_ ## Axis) { \
- log_err(filename, node, \
- "the ellipsoid \""STR(Axis)"\" parameter is alread defined.\n"); \
- return RES_BAD_ARG; \
- } \
- mask |= LENGTH_ ## Axis; \
- res = parse_yaml_param_distribution \
- (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.Axis); \
- } (void)0
- else if(!strcmp((char*)key->data.scalar.value, "a")) {
- PARSE_SEMI_PRINCIPAL_AXIS(a);
- } else if(!strcmp((char*)key->data.scalar.value, "b")) {
- PARSE_SEMI_PRINCIPAL_AXIS(b);
+ /* equivalent sphere radius */
+ } else if(!strcmp((char*)key->data.scalar.value, "radius")) {
+ SETUP_MASK(RADIUS, "equivalent sphere radius of the ellipsoid");
+ res = parse_yaml_param_distribution
+ (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.radius);
+
+ } else if(!strcmp((char*)key->data.scalar.value, "aspect_ratio")) {
+ SETUP_MASK(ASPECT_RATIO, "a/b aspect ratio of the ellipsoid");
+ res = parse_yaml_double
+ (filename, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.a_over_c);
+
+ /* Length of the ellipsoid "a" semi-principal axis */
+ } else if(!strcmp((char*)key->data.scalar.value, "a")) {
+ SETUP_MASK(LENGTH_a, "ellipsoid \"a\" parameter");
+ res = parse_yaml_param_distribution
+ (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.a);
+
+ /* Length of the ellipsoid "c" semi-principal axis */
} else if(!strcmp((char*)key->data.scalar.value, "c")) {
- PARSE_SEMI_PRINCIPAL_AXIS(c);
- }
- #undef PARSE_SEMI_PRINCIPAL_AXIS
+ SETUP_MASK(LENGTH_c, "ellipsoid \"c\" parameter");
+ res = parse_yaml_param_distribution
+ (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.c);
/* Error */
- else {
- log_err(filename, key, "unkown sphere attribute `%s'.\n",
+ } else {
+ log_err(filename, key, "unknown sphere attribute `%s'.\n",
key->data.scalar.value);
return RES_BAD_ARG;
}
if(res != RES_OK) return res;
}
+ #undef SETUP_MASK
+
+ /* Ensure the completeness of the ellipsoid distribution */
+ if((mask & (LENGTH_a | LENGTH_c))) {
+ if(mask & (RADIUS | ASPECT_RATIO)) {
+ log_err(filename, node,
+ "the ellipsoid semi-principal axis and its equivalent sphere cannot "
+ "be defined simultaneously.\n");
+ return RES_BAD_ARG;
+ } else if(!(mask & LENGTH_a)) {
+ log_err(filename, node,
+ "missing the length of the semi principal axis \"a\".\n");
+ return RES_BAD_ARG;
+ } else if(!(mask & LENGTH_c)) {
+ log_err(filename, node,
+ "missing the length of the semi principal axis \"c\".\n");
+ return RES_BAD_ARG;
+ }
+ geom->type = SCHIFF_ELLIPSOID;
+ } else if (mask & (RADIUS | ASPECT_RATIO)) {
+ if(mask & (LENGTH_a | LENGTH_c)) {
+ log_err(filename, node,
+ "the ellipsoid semi-principal axis and its equivalent sphere cannot "
+ "be defined simultaneously.\n");
+ return RES_BAD_ARG;
+ } else if(!(mask & RADIUS)) {
+ log_err(filename, node, "missing the equivalent sphere radius.\n");
+ return RES_BAD_ARG;
+ } else if(!(mask & ASPECT_RATIO)) {
+ log_err(filename, node, "missing the a/c ratio of the ellipsoid.\n");
+ return RES_BAD_ARG;
+ }
+ geom->type = SCHIFF_ELLIPSOID_AS_SPHERE;
+ }
- /* Ensure the completness of the ellipsoid distribution */
- #define CHECK_SEMI_PRINCIPAL_AXIS(Axis) { \
- if(!(mask & LENGTH_ ## Axis)) { \
- log_err(filename, node, \
- "missing the length of the semi principal axis \""STR(Axis)"\".\n"); \
- return RES_BAD_ARG; \
- } \
- } (void)0
- CHECK_SEMI_PRINCIPAL_AXIS(a);
- CHECK_SEMI_PRINCIPAL_AXIS(b);
- CHECK_SEMI_PRINCIPAL_AXIS(c);
- #undef CHECK_SEMI_PRINCIPAL_AXIS
if(!(mask & PROBA)) { /* Default proba */
*geom_proba = 1.0;
}
if(!(mask & SLICES)) { /* Default number of slices */
geom->data.ellipsoid.nslices = SCHIFF_ELLIPSOID_DEFAULT.nslices;
}
- geom->type = SCHIFF_ELLIPSOID;
return RES_OK;
}
diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c
@@ -262,16 +262,27 @@ geometry_sample_ellipsoid
const struct schiff_ellipsoid* ellipsoid;
struct mesh_context mesh_ctx;
struct s3d_vertex_data attrib;
- double a, b, c, a2, b2, c2;
+ double a, b, c, a2, b2, c2, r_eq;
size_t nverts, nprims;
ASSERT(rng && geom && mesh && shape);
- ASSERT(geom->type == SCHIFF_ELLIPSOID);
ellipsoid = &geom->data.ellipsoid;
+ switch(geom->type) {
+ case SCHIFF_ELLIPSOID:
+ a = eval_param(&ellipsoid->a, rng);
+ c = eval_param(&ellipsoid->c, rng);
+ break;
+ case SCHIFF_ELLIPSOID_AS_SPHERE:
+ r_eq = eval_param(&ellipsoid->radius, rng);
+ /* a_over_c^(1/3) could be precomputed */
+ a = r_eq * pow(ellipsoid->a_over_c, 1.0/3.0);
+ c = a / ellipsoid->a_over_c;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
- a = eval_param(&ellipsoid->a, rng);
- b = eval_param(&ellipsoid->b, rng);
- c = eval_param(&ellipsoid->c, rng);
+ /* Schiff exposes ellipsoid whose a and b semi-principal axis are equals */
+ b = a;
a2 = a*a;
b2 = b*b;
c2 = c*c;
@@ -318,7 +329,7 @@ geometry_sample_cylinder
break;
case SCHIFF_CYLINDER_AS_SPHERE:
r = eval_param(&cylinder->radius, rng);
- mesh_ctx.data.cylinder.radius =
+ mesh_ctx.data.cylinder.radius = /* The denominator can be pre-computed */
(float)(r / pow(3.0 / (2.0*cylinder->aspect_ratio), 1.0/3.0));
mesh_ctx.data.cylinder.height =
(float)(2.f * mesh_ctx.data.cylinder.radius / cylinder->aspect_ratio);
@@ -422,6 +433,7 @@ geometry_sample
switch(geom->type) {
case SCHIFF_ELLIPSOID:
+ case SCHIFF_ELLIPSOID_AS_SPHERE:
res = geometry_sample_ellipsoid(rng, geom, mesh, shape);
break;
case SCHIFF_CYLINDER:
@@ -516,6 +528,7 @@ schiff_geometry_distribution_init
/* Generate the mesh template and setup its distribution context */
switch(geoms[i].type) {
case SCHIFF_ELLIPSOID:
+ case SCHIFF_ELLIPSOID_AS_SPHERE:
res = schiff_mesh_init_sphere_polar(&mem_default_allocator,
&ctx->meshes[i], geoms[i].data.ellipsoid.nslices);
break;
diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h
@@ -38,22 +38,30 @@ struct schiff_param {
enum schiff_geometry_type {
SCHIFF_ELLIPSOID,
+ SCHIFF_ELLIPSOID_AS_SPHERE,/* The ellipsoid volume is controled by a sphere */
SCHIFF_CYLINDER,
- SCHIFF_CYLINDER_AS_SPHERE, /* The cylinder volume is control by a sphere */
+ SCHIFF_CYLINDER_AS_SPHERE,/* The cylinder volume is controled by a sphere */
SCHIFF_SPHERE,
SCHIFF_SUPER_SHAPE,
SCHIFF_NONE
};
+/* (x/a)^2 + (y/a)^2 + (z/c)^2 = 1 */
struct schiff_ellipsoid {
+ /* In use by SCHIFF_ELLIPSOID */
struct schiff_param a;
- struct schiff_param b;
struct schiff_param c;
+
+ /* In use by SCHIFF_ELLIPSOID_AS_SPHERE
+ * a = radius * a_over_c^1/3; c = a / a_over_c */
+ struct schiff_param radius; /* Radius of an equivalent sphere */
+ double a_over_c; /* Fixed aspect ratio */
+
unsigned nslices;
};
#define SCHIFF_ELLIPSOID_DEFAULT__ \
- {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, 64}
+ {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, 0, 64}
static const struct schiff_ellipsoid SCHIFF_ELLIPSOID_DEFAULT =
SCHIFF_ELLIPSOID_DEFAULT__;
