QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Typedefs for particle and T-matrix indices. 

Former-commit-id: 797052ad7a82a9f7a7186bcaab0bfea7bb962933
This commit is contained in:
Marek Nečada 2019-02-27 11:53:33 +02:00
parent 88b3498cb3
commit a1d14ca1b8
3 changed files with 30 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
qpms/qpms_types.h
View File

@ -6,7 +6,7 @@
#include <complex.h>
#include <stdbool.h>
#include <stddef.h>
//#include <stdint.h>
#include <stdint.h>
#ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
@ -302,5 +302,11 @@ typedef struct qpms_vswf_set_spec_t {
qpms_normalisation_t norm; ///< Normalisation convention. To be set manually if needed.
} qpms_vswf_set_spec_t;
/// T-matrix index used in qpms_scatsys_t and related structures.
typedef int32_t qpms_ss_tmi_t;
/// Particle index used in qpms_scatsys_t and related structures.
typedef int32_t qpms_ss_pi_t;
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
#endif // QPMS_TYPES

47
qpms/scatsystem.c
View File

@ -266,7 +266,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
{
double minx = +INFINITY, miny = +INFINITY, minz = +INFINITY;
double maxx = -INFINITY, maxy = -INFINITY, maxz = -INFINITY;
for (ssize_t i = 0; i < orig->p_count; ++i) {
for (qpms_ss_pi_t i = 0; i < orig->p_count; ++i) {
minx = MIN(minx,orig->p[i].pos.x);
miny = MIN(miny,orig->p[i].pos.y);
minz = MIN(minz,orig->p[i].pos.z);
@ -278,8 +278,8 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
}
// Second, check that there are no duplicit positions in the input system.
for (ssize_t i = 0; i < orig->p_count; ++i)
for (ssize_t j = 0; j < i; ++j)
for (qpms_ss_pi_t i = 0; i < orig->p_count; ++i)
for (qpms_ss_pi_t j = 0; j < i; ++j)
assert(!cart3_isclose(orig->p[i].pos, orig->p[j].pos, 0, QPMS_SCATSYS_LEN_RTOL * lenscale));
// Allocate T-matrix and particle arrays
@ -294,11 +294,11 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
// Copy T-matrices; checking for duplicitie
ptrdiff_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities
qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities
ss->tm_count = 0;
for (ptrdiff_t i = 0; i < orig->tm_count; ++i) {
for (qpms_ss_tmi_t i = 0; i < orig->tm_count; ++i) {
bool found_dupl = false;
ptrdiff_t j;
qpms_ss_tmi_t j;
for (j = 0; j < ss->tm_count; ++j)
if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) {
break;
@ -311,23 +311,23 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
}
// Copy particles, remapping the t-matrix indices
for (ptrdiff_t i = 0; i < orig->p_count; ++(i)) {
for (qpms_ss_pi_t i = 0; i < orig->p_count; ++(i)) {
ss->p[i] = orig->p[i];
ss->p[i].tmatrix_id = tm_dupl_remap[ss->p[i].tmatrix_id];
}
ss->p_count = orig->p_count;
// allocate t-matrix symmetry map
ss->tm_sym_map = malloc(sizeof(ptrdiff_t) * sym->order * sym->order * ss->tm_count);
ss->tm_sym_map = malloc(sizeof(qpms_ss_tmi_t) * sym->order * sym->order * ss->tm_count);
// Extend the T-matrices list by the symmetry operations
for (ptrdiff_t tmi = 0; tmi < ss->tm_count; ++tmi)
for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi)
for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){
const size_t d = ss->tm[tmi]->spec->n;
complex double M[d][d]; // transformation matrix
qpms_irot3_uvswfi_dense(M[0], ss->tm[tmi]->spec, sym->rep3d[gmi]);
qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ss->tm[tmi], M[0]);
ptrdiff_t tmj;
qpms_ss_tmi_t tmj;
for (tmj = 0; tmj < ss->tm_count; ++tmj)
if (qpms_tmatrix_isclose(transformed, ss->tm[tmj], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL))
break;
@ -340,19 +340,19 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix
}
// Possibly free some space using the new ss->tm_count instead of (old) ss->tm_count*sym->order
ss->tm_sym_map = realloc(ss->tm_sym_map, sizeof(ptrdiff_t) * sym->order * ss->tm_count);
ss->tm_sym_map = realloc(ss->tm_sym_map, sizeof(qpms_ss_tmi_t) * sym->order * ss->tm_count);
// tm could be realloc'd as well, but those are just pointers, not likely many.
// allocate particle symmetry map
ss->p_sym_map = malloc(sizeof(ptrdiff_t) * sym->order * sym->order * ss->p_count);
ss->p_sym_map = malloc(sizeof(qpms_ss_pi_t) * sym->order * sym->order * ss->p_count);
// Extend the particle list by the symmetry operations, check that particles mapped by symmetry ops on themselves
// have the correct symmetry
// TODO this could be sped up to O(npart * log(npart)); let's see later if it is needed.
for (ptrdiff_t pi = 0; pi < ss->p_count; ++pi)
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi)
for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi) {
cart3_t newpoint = qpms_irot3_apply_cart3(sym->rep3d[gmi], ss->p[pi].pos);
ptrdiff_t new_tmi = ss->tm_sym_map[gmi + ss->p[pi].tmatrix_id * sym->order]; // transformed t-matrix index
ptrdiff_t pj;
qpms_ss_tmi_t new_tmi = ss->tm_sym_map[gmi + ss->p[pi].tmatrix_id * sym->order]; // transformed t-matrix index
qpms_ss_pi_t pj;
for (pj = 0; pj < ss->p_count; ++pj)
if (cart3_isclose(newpoint, ss->p[pj].pos, 0, ss->lenscale * QPMS_SCATSYS_LEN_RTOL)) {
if (new_tmi != ss->p[pj].tmatrix_id)
@ -369,28 +369,13 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
ss->p_sym_map[gmi + pi * sym->order] = pi;
}
// Possibly free some space using the new ss->p_count instead of (old) ss->p_count*sym->order
ss->p_sym_map = realloc(ss->p_sym_map, sizeof(ptrdiff_t) * sym->order * ss->p_count);
ss->p_sym_map = realloc(ss->p_sym_map, sizeof(qpms_ss_pi_t) * sym->order * ss->p_count);
ss->p = realloc(ss->p, sizeof(qpms_particle_tid_t) * ss->p_count); ss->p_capacity = ss->p_count;
ss->sym = sym;
return ss;
}

14
qpms/scatsystem.h
View File

@ -193,7 +193,7 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
typedef struct qpms_particle_tid_t {
// Does it make sense to ever use other than cartesian coords for this?
cart3_t pos; ///< Particle position in cartesian coordinates.
ptrdiff_t tmatrix_id; ///< T-matrix index
qpms_ss_tmi_t tmatrix_id; ///< T-matrix index
} qpms_particle_tid_t;
struct qpms_finite_group_t;
@ -201,17 +201,17 @@ struct qpms_finite_group_t;
typedef struct qpms_scatsys_t {
// TODO does bspec belong here?
qpms_tmatrix_t **tm; ///< T-matrices in the system
ptrdiff_t tm_count; ///< Number of all different T-matrices
ptrdiff_t tm_capacity; ///< Capacity of tm[].
qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices
qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[].
qpms_particle_tid_t *p; ///< Particles.
ptrdiff_t p_count; ///< Size of particles array.
ptrdiff_t p_capacity; ///< Capacity of p[].
qpms_ss_pi_t p_count; ///< Size of particles array.
qpms_ss_pi_t p_capacity; ///< Capacity of p[].
//TODO the index types do not need to be so big.
struct qpms_finite_group_t *sym; ///< Symmetry group of the array
ptrdiff_t *p_sym_map; ///< Which particles map onto which by the symmetry ops.
qpms_ss_pi_t *p_sym_map; ///< Which particles map onto which by the symmetry ops.
///< p_sym_map[idi + pi * sym->order] gives the index of pi-th particle under the idi'th sym op.
ptrdiff_t *tm_sym_map; ///< Which t-matrices map onto which by the symmetry ops. Lookup by tm_sum_map[idi + tmi * sym->order].
qpms_ss_tmi_t *tm_sym_map; ///< Which t-matrices map onto which by the symmetry ops. Lookup by tm_sum_map[idi + tmi * sym->order].
// TODO shifted origin of the symmetry group etc.
// TODO some indices for fast operations here.