From 08b8a1a315e0b7b376cc7dfd430a1d5d5a1824f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sat, 7 Sep 2019 12:39:17 +0300 Subject: [PATCH] WIP generate 3D lattices Former-commit-id: 16c6f7779db68a3169dbc77a6adb9ffdd6ad374a --- qpms/latticegens.c | 167 +++++++++++++++++++++++++++++++++++++++++++++ qpms/lattices.h | 10 ++- 2 files changed, 176 insertions(+), 1 deletion(-) diff --git a/qpms/latticegens.c b/qpms/latticegens.c index 6eed602..295053a 100644 --- a/qpms/latticegens.c +++ b/qpms/latticegens.c @@ -585,3 +585,170 @@ const PGenClassInfo PGen_xyWeb = { }; +//==== PGen_xyzWeb ==== +// 3D lattice generator in the "spiderweb" style, generated in the "perimetre" order, +// not strictly ordered (or limited) by distance from origin. +// The minR and maxR here refer to the TODO WWHAT + +extern const PGenClassInfo PGen_xyzWeb; // forward declaration needed by constructor (may be placed in header file instead) + +// Internal state structure +typedef struct PGen_xyzWeb_StateData { + long i[6]; // real size given by nsteps + unsigned short phase; // 0 to 5 + long layer; + long last_layer; // generation stops when layer > last_layer + double layer_min_height; // this * layer is what minR and maxR are compared to + double minR, maxR; + bool inc_minR, inc_maxR; + cart3_t b[3]; // lattice vectors (TODO needed?) + cart3_t steps[6]; // generated by l3d_shortestBase356 + int nsteps; // real size of steps (3, 5, or 6), depending on lattice type + cart3_t offset; // offset of the zeroth lattice point from origin (TODO will be normalised to the WS cell) + // TODO type rectangular vs. triangular + LatticeFlags lf; +} PGen_xyzWeb_StateData; + +// Constructor +PGen PGen_xyzWeb_new(const cart3_t b[3], double rtol, cart2_t offset, double minR, bool inc_minR, double maxR, bool inc_maxR) { + PGen_xyzWeb_StateData *s = calloc(1, sizeof(PGen_xyzWeb_StateData)); + s->minR = minR; s->maxR = maxR; + s->inc_minR = inc_minR; + s->inc_maxR = inc_maxR; + l3d_reduceBasis(b, &(s->b)); + s->nsteps = l3d_shortestBase356(s->b, s->steps); + s->offset = offset; // TODO shorten into the WS cell ? + s->lf = l3d_detectRightAnles(s->b, rtol); + s->layer_min_height = l3d_inSphereRadius(s->b); + s->layer = ceil(s->minR/s->layer_min_height); + if(!inc_minR && (s->layer * s->layer_min_height) <= minR) + ++(s->layer); + s->i = s->layer; s->j = 0; s->phase = 0; // init indices + s->last_layer = floor(s->maxR/s->layer_min_height); + if(!inc_maxR && (s->last_layer * s->layer_min_height) >= maxR) + --(s->last_layer); + PGen g = {&PGen_xyzWeb, (void *) s}; + return g; +} + +// Destructor +void PGen_xyzWeb_destructor(PGen *g) { + free(g->stateData); + g->stateData = NULL; +} + +// Extractor (3D cartesian, native) +PGenCart3ReturnData PGen_xyzWeb_next_cart3(PGen *g) { + if (g->stateData == NULL) // already destroyed + return PGenCart2DoneVal; + else { + PGen_xyzWeb_StateData * const s = (PGen_xyzWeb_StateData *) g->stateData; + assert(s->layer >= 0); + if (s->layer <= s->last_layer) { + const cart2_t thePoint = cart2_add(s->offset, + cart2_add(cart2_scale(s->i, s->b1), cart2_scale(s->j, s->b2))); + if(s->layer == 0) { // origin is unique, proceed with next layer + ++s->layer; + s->phase = 0; + s->i = s->layer; + s->j = 0; + } + else if(s->lt & (SQUARE | RECTANGULAR)) { + // rectangular or square lattice, four perimeters + switch(s->phase) { + case 0: // initial i = l, j = 0 + --s->i; + ++s->j; + if(s->i <= 0) ++s->phase; + break; + case 1: // initial i = 0, j = l + --s->i; + --s->j; + if(s->j <= 0) ++s->phase; + break; + case 2: // initial i = -l, j = 0 + ++s->i; + --s->j; + if(s->i >= 0) ++s->phase; + break; + case 3: // initial i = 0, j = -l + ++s->i; + ++s->j; + if(s->j >= 0) ++s->phase; + break; + default: + abort(); + } + if(s->phase == 4) { // phase overflow, start new layer + ++s->layer; + s->phase = 0; + s->i = s->layer; + s->j = 0; + } + } else { // non-rectangular lattice, six perimeters + switch(s->phase) { + case 0: + --s->i; + ++s->j; + if(s->i <= 0) ++s->phase; + break; + case 1: + --s->i; + if(s->i + s->j <= 0) ++s->phase; + break; + case 2: + --s->j; + if(s->j <= 0) ++s->phase; + break; + case 3: + ++s->i; + --s->j; + if(s->i >= 0) ++s->phase; + break; + case 4: + ++s->i; + if(s->i + s->j >= 0) ++s->phase; + break; + case 5: + ++s->j; + if(s->j >= 0) ++s->phase; + break; + default: + abort(); + } + if(s->phase == 6) { // phase overflow, start next layer + ++s->layer; + s->phase = 0; + s->i = s->layer; + s->j = 0; + } + } + PGenCart2ReturnData retval = {(PGEN_NOTDONE | PGEN_AT_XY | PGEN_COORDS_CART2), thePoint}; + return retval; + } else { + PGen_destroy(g); + return PGenCart2DoneVal; + } + } +} + +// Class metadata structure; TODO maybe this can rather be done by macro. +const PGenClassInfo PGen_xyzWeb = { + "PGen_xyzWeb", + 3, + PGEN_COORDS_CART3, + NULL,//PGen_xyzWeb_next, + NULL,//PGen_xyzWeb_next_z, + NULL,//PGen_xyzWeb_next_pol, + PGen_next_sph_from_cart3, //NULL,//PGen_xyzWeb_next_sph, + NULL,//PGen_xyzWeb_next_cart2, + PGen_xyzWeb_next_cart3, // native + NULL,//PGen_xyzWeb_fetch, + NULL,//PGen_xyzWeb_fetch_z, + NULL,//PGen_xyzWeb_fetch_pol, + NULL,//PGen_xyzWeb_fetch_sph, + NULL,//PGen_xyzWeb_fetch_cart2, + NULL,//PGen_xyzWeb_fetch_cart3, + PGen_xyzWeb_destructor +}; + diff --git a/qpms/lattices.h b/qpms/lattices.h index bbf515a..4ed41e7 100644 --- a/qpms/lattices.h +++ b/qpms/lattices.h @@ -620,7 +620,11 @@ PGen PGen_1D_new_minMaxR(double period, ///< Distance between points. extern const PGenClassInfo PGen_xyWeb; -PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double rtol, cart2_t offset, +PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double orthogonality_rtol, cart2_t offset, + double minR, bool inc_minR, double maxR, bool inc_maxR); + +extern const PGenClassInfo PGen_xyzWeb; +PGen PGen_xyzWeb_new(const cart3_t basis[3], double orthogonality_rtol, cart3_t offset, double minR, bool inc_minR, double maxR, bool inc_maxR); @@ -691,6 +695,8 @@ int l2d_shortestBase46(cart2_t i1, cart2_t i2, cart2_t *o1, cart2_t *o2, cart2_ // variant int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2, cart2_t *oarr, double rtol); +int l3d_shortestBase356(const cart3_t basis[3], cart3_t shortest356[6], double rtol); + // Determines whether angle between inputs is obtuse bool l2d_is_obtuse_r(cart2_t i1, cart2_t i2, double rtol); bool l2d_is_obtuse(cart2_t i1, cart2_t i2); @@ -729,6 +735,8 @@ double l2d_unitcell_area(cart2_t b1, cart2_t b2); // returns the radius of inscribed circle of a hexagon (or rectangle/square if applicable) created by the shortest base triple double l2d_hexWebInCircleRadius(cart2_t b1, cart2_t b2); +double l3d_inSphereRadius(const cart3_t b[3]); + /* * THE MORE OR LESS OK PART * ========================