WIP generate 3D lattices

Former-commit-id: 16c6f7779db68a3169dbc77a6adb9ffdd6ad374a
This commit is contained in:
Marek Nečada 2019-09-07 12:39:17 +03:00
parent 658c564a0c
commit 08b8a1a315
2 changed files with 176 additions and 1 deletions

View File

@ -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
};

View File

@ -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
* ========================