From e84483b2e8eb413a04b6db5889095a962db07350 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 26 Aug 2018 19:36:45 +0000 Subject: [PATCH] A stupid shift function for points2d_rordered_t Former-commit-id: bbba5408081965d8550d78fed5b7569c5cff95b8 --- qpms/lattices.h | 24 +++++- qpms/lattices2d.c | 103 ++++++++++++++++++++++++++ tests/lattice/2d_triangular_lattice.c | 5 ++ 3 files changed, 128 insertions(+), 4 deletions(-) diff --git a/qpms/lattices.h b/qpms/lattices.h index 987655a..d7012d9 100644 --- a/qpms/lattices.h +++ b/qpms/lattices.h @@ -105,7 +105,7 @@ int l2d_cellCornersWS_arr(cart2_t i1, cart2_t i2, cart2_t *oarr, double rtol); // Reciprocal bases void l2d_reciprocalBasis1(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2); -void l2D_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2); +void l2d_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2); /* @@ -133,10 +133,26 @@ typedef struct { */ } points2d_rordered_t; -// returns a copy but scaled by a factor -points2d_rordered_t *points2d_rordered_scale(const points2d_rordered_t *orig, double factor); -void points2d_rordered_free(points2d_rordered_t *); // use only for result of points2d_rordered_scale +// sorts arbitrary points and creates points2d_rordered_t +points2d_rordered_t *points2d_rordered_frompoints(const point2d *orig_base, + size_t nmemb, double rtol, double atol); + +// returns a copy but shifted by a constant (actually in a stupid way, but whatever) +points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig, + point2d shift, double rtol, double atol); + +// returns a copy but scaled by a factor +points2d_rordered_t *points2d_rordered_scale(const points2d_rordered_t *orig, + double factor); + + +/* The destructor: use only for results of + * - points2D_rordered_frompoints, + * - points2d_rordered_shift, + * - points2d_rordered_scale. + */ +void points2d_rordered_free(points2d_rordered_t *); static inline point2d points2d_rordered_get_point(const points2d_rordered_t *ps, int r_order, int i) { assert(i >= 0); diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c index f161a48..657843e 100644 --- a/qpms/lattices2d.c +++ b/qpms/lattices2d.c @@ -1,6 +1,7 @@ #include "lattices.h" #include #include +#include typedef struct { int i, j; @@ -8,6 +9,7 @@ typedef struct { static inline int sqi(int x) { return x*x; } +static inline double sqf(double x) { return x*x; } void points2d_rordered_free(points2d_rordered_t *p) { free(p->rs); @@ -80,6 +82,107 @@ points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig, return p; } + +static inline double pr2(const point2d p) { + return sqf(p.x) + sqf(p.y); +} + +static inline double prn(const point2d p) { + return sqrt(pr2(p)); +} + +static int point2d_cmp_by_r2(const void *p1, const void *p2) { + const point2d *z1 = (point2d *) p1, *z2 = (point2d *) p2; + double dif = pr2(*z1) - pr2(*z2); + if(dif > 0) return 1; + else if(dif < 0) return -1; + else return 0; +} + +static points2d_rordered_t *points2d_rordered_frompoints_c(point2d *orig_base, const size_t nmemb, + const double rtol, const double atol, bool copybase) +{ + // TODO should the rtol and atol relate to |r| or r**2? (Currently: |r|) + assert(rtol >= 0); + assert(atol >= 0); + points2d_rordered_t *p = malloc(sizeof(points2d_rordered_t)); + if(nmemb == 0) { + p->nrs = 0; + p->rs = NULL; + p->base = NULL; + p->r_offsets = NULL; + return p; + } + + if (copybase) { + p->base = malloc(nmemb * sizeof(point2d)); + memcpy(p->base, orig_base, nmemb * sizeof(point2d)); + } else + p->base = orig_base; + + qsort(p->base, nmemb, sizeof(point2d), point2d_cmp_by_r2); + + // first pass: determine the number of "different" r's. + size_t rcount = 0; + double rcur = -INFINITY; + double rcurmax = -INFINITY; + for (size_t i = 0; i < nmemb; ++i) + if ((rcur = prn(p->base[i])) > rcurmax) { + ++rcount; + rcurmax = rcur * (1 + rtol) + atol; + } + + p->nrs = rcount; + // TODO check malloc return values + p->rs = malloc(rcount * sizeof(double)); + p->r_offsets = malloc((rcount+1) * sizeof(ptrdiff_t)); + + + // second pass: fill teh rs; + size_t ri = 0; + size_t rcurcount = 0; + rcur = prn(p->base[0]); + rcurmax = rcur * (1 + rtol) + atol; + double rcursum = 0; + p->r_offsets[0] = 0; + for (size_t i = 0; i < nmemb; ++i) { + rcur = prn(p->base[i]); + if (rcur > rcurmax) { + p->rs[ri] = rcursum / rcurcount; // average of the accrued r's within tolerance + ++ri; + p->r_offsets[ri] = i; //r_offsets[ri-1] + rcurcount (is the same) + rcurcount = 0; + rcursum = 0; + rcurmax = rcur * (1 + rtol) + atol; + } + rcursum += rcur; + ++rcurcount; + } + p->r_offsets[rcount] = nmemb; + + return p; +} + +points2d_rordered_t *points2d_rordered_frompoints(const point2d *orig_base, const size_t nmemb, + const double rtol, const double atol) +{ + return points2d_rordered_frompoints_c((point2d *)orig_base, nmemb, rtol, atol, true); +} + +points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig, const point2d shift, + double rtol, double atol) +{ + size_t n = (orig->nrs > 0) ? + orig->r_offsets[orig->nrs] - orig->r_offsets[0] : 0; + + point2d * shifted = malloc(n * sizeof(point2d)); + + for(size_t i = 0; i < n; ++i) + shifted[i] = cart2_add(orig->base[i], shift); + + return points2d_rordered_frompoints_c(shifted, n,rtol, atol, false); +} + /* * EQUILATERAL TRIANGULAR LATTICE */ diff --git a/tests/lattice/2d_triangular_lattice.c b/tests/lattice/2d_triangular_lattice.c index 7d9ef74..bb196d6 100644 --- a/tests/lattice/2d_triangular_lattice.c +++ b/tests/lattice/2d_triangular_lattice.c @@ -1,5 +1,6 @@ #include #include +#include void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) { FILE *f = fopen(filename, "w"); @@ -65,6 +66,10 @@ int main() { dump_points2d_rordered(&(h->ps), "hex_h_empty.out"); honeycomb_lattice_gen_extend_to_steps(h, 7); dump_points2d_rordered(&(h->ps), "hex_h_s7.out"); + point2d shift = {0, h->h}; + p = points2d_rordered_shift(&(h->ps), shift, DBL_EPSILON, h->h * DBL_EPSILON); + dump_points2d_rordered(p, "hex_h_s7_shifted.out"); + points2d_rordered_free(p); honeycomb_lattice_gen_extend_to_steps(h, 120); dump_points2d_rordered(&(h->ps), "hex_h_s120.out"); honeycomb_lattice_gen_free(h);