APBS  1.5
vcsm.c
Go to the documentation of this file.
1 
57 #include "vcsm.h"
58 
59 /* Inlineable methods */
60 #if !defined(VINLINE_VCSM)
61 
62 VPUBLIC Valist* Vcsm_getValist(Vcsm *thee) {
63 
64  VASSERT(thee != VNULL);
65  return thee->alist;
66 
67 }
68 
69 VPUBLIC int Vcsm_getNumberAtoms(Vcsm *thee, int isimp) {
70 
71  VASSERT(thee != VNULL);
72  VASSERT(thee->initFlag);
73  return thee->nsqm[isimp];
74 
75 }
76 
77 VPUBLIC Vatom* Vcsm_getAtom(Vcsm *thee, int iatom, int isimp) {
78 
79 
80  VASSERT(thee != VNULL);
81  VASSERT(thee->initFlag);
82 
83  VASSERT(iatom < (thee->nsqm)[isimp]);
84  return Valist_getAtom(thee->alist, (thee->sqm)[isimp][iatom]);
85 
86 }
87 
88 VPUBLIC int Vcsm_getAtomIndex(Vcsm *thee, int iatom, int isimp) {
89 
90 
91  VASSERT(thee != VNULL);
92  VASSERT(thee->initFlag);
93 
94  VASSERT(iatom < (thee->nsqm)[isimp]);
95  return (thee->sqm)[isimp][iatom];
96 
97 }
98 
99 VPUBLIC int Vcsm_getNumberSimplices(Vcsm *thee, int iatom) {
100 
101 
102  VASSERT(thee != VNULL);
103  VASSERT(thee->initFlag);
104 
105  return (thee->nqsm)[iatom];
106 
107 }
108 
109 VPUBLIC SS* Vcsm_getSimplex(Vcsm *thee, int isimp, int iatom) {
110 
111 
112  VASSERT(thee != VNULL);
113  VASSERT(thee->initFlag);
114 
115  return Gem_SS(thee->gm, (thee->qsm)[iatom][isimp]);
116 
117 }
118 
119 VPUBLIC int Vcsm_getSimplexIndex(Vcsm *thee, int isimp, int iatom) {
120 
121 
122  VASSERT(thee != VNULL);
123  VASSERT(thee->initFlag);
124 
125  return (thee->qsm)[iatom][isimp];
126 
127 }
128 
129 VPUBLIC unsigned long int Vcsm_memChk(Vcsm *thee) {
130  if (thee == VNULL) return 0;
131  return Vmem_bytes(thee->vmem);
132 }
133 
134 #endif /* if !defined(VINLINE_VCSM) */
135 
136 VPUBLIC Vcsm* Vcsm_ctor(Valist *alist, Gem *gm) {
137 
138  /* Set up the structure */
139  Vcsm *thee = VNULL;
140  thee = (Vcsm*)Vmem_malloc(VNULL, 1, sizeof(Vcsm) );
141  VASSERT( thee != VNULL);
142  VASSERT( Vcsm_ctor2(thee, alist, gm));
143 
144  return thee;
145 }
146 
147 VPUBLIC int Vcsm_ctor2(Vcsm *thee, Valist *alist, Gem *gm) {
148 
149  VASSERT( thee != VNULL );
150 
151  /* Memory management object */
152  thee->vmem = Vmem_ctor("APBS:VCSM");
153 
154  /* Set up the atom list and grid manager */
155  if( alist == VNULL) {
156  Vnm_print(2,"Vcsm_ctor2: got null pointer to Valist object!\n");
157  return 0;
158  }
159  thee->alist = alist;
160  if( gm == VNULL) {
161  Vnm_print(2,"Vcsm_ctor2: got a null pointer to the Gem object!\n");
162  return 0;
163  }
164  thee->gm = gm;
165 
166  thee->initFlag = 0;
167  return 1;
168 }
169 
170 VPUBLIC void Vcsm_init(Vcsm *thee) {
171 
172  /* Counters */
173  int iatom,
174  jatom,
175  isimp,
176  jsimp,
177  gotSimp;
178  /* Atomic information */
179  Vatom *atom;
180  double *position;
181  /* Simplex/Vertex information */
182  SS *simplex;
183  /* Basis function values */
184 
185  if (thee == VNULL) {
186  Vnm_print(2, "Vcsm_init: Error! Got NULL thee!\n");
187  VASSERT(0);
188  }
189  if (thee->gm == VNULL) {
190  Vnm_print(2, "Vcsm_init: Error! Got NULL thee->gm!\n");
191  VASSERT(0);
192  }
193  thee->nsimp = Gem_numSS(thee->gm);
194  if (thee->nsimp <= 0) {
195  Vnm_print(2, "Vcsm_init: Error! Got %d simplices!\n", thee->nsimp);
196  VASSERT(0);
197  }
198  thee->natom = Valist_getNumberAtoms(thee->alist);
199 
200  /* Allocate and initialize space for the first dimensions of the
201  * simplex-charge map, the simplex array, and the counters */
202  thee->sqm = (int**)Vmem_malloc(thee->vmem, thee->nsimp, sizeof(int *));
203  VASSERT(thee->sqm != VNULL);
204  thee->nsqm = (int*)Vmem_malloc(thee->vmem, thee->nsimp, sizeof(int));
205  VASSERT(thee->nsqm != VNULL);
206  for (isimp=0; isimp<thee->nsimp; isimp++) (thee->nsqm)[isimp] = 0;
207 
208  /* Count the number of charges per simplex. */
209  for (iatom=0; iatom<thee->natom; iatom++) {
210  atom = Valist_getAtom(thee->alist, iatom);
211  position = Vatom_getPosition(atom);
212  gotSimp = 0;
213  for (isimp=0; isimp<thee->nsimp; isimp++) {
214  simplex = Gem_SS(thee->gm, isimp);
215  if (Gem_pointInSimplex(thee->gm, simplex, position)) {
216  (thee->nsqm)[isimp]++;
217  gotSimp = 1;
218  }
219  }
220  }
221 
222  /* @todo Combine the following two loops? - PCE */
223  /* Allocate the space for the simplex-charge map */
224  for (isimp=0; isimp<thee->nsimp; isimp++) {
225  if ((thee->nsqm)[isimp] > 0) {
226  thee->sqm[isimp] = (int*)Vmem_malloc(thee->vmem, (thee->nsqm)[isimp],
227  sizeof(int));
228  VASSERT(thee->sqm[isimp] != VNULL);
229  }
230  }
231 
232  /* Finally, set up the map */
233  for (isimp=0; isimp<thee->nsimp; isimp++) {
234  jsimp = 0;
235  simplex = Gem_SS(thee->gm, isimp);
236  for (iatom=0; iatom<thee->natom; iatom++) {
237  atom = Valist_getAtom(thee->alist, iatom);
238  position = Vatom_getPosition(atom);
239  /* Check to see if the atom's in this simplex */
240  if (Gem_pointInSimplex(thee->gm, simplex, position)) {
241  /* Assign the entries in the next vacant spot */
242  (thee->sqm)[isimp][jsimp] = iatom;
243  jsimp++;
244  }
245  }
246  }
247 
248  thee->msimp = thee->nsimp;
249 
250  /* Allocate space for the charge-simplex map */
251  thee->qsm = (int**)Vmem_malloc(thee->vmem, thee->natom, sizeof(int *));
252  VASSERT(thee->qsm != VNULL);
253  thee->nqsm = (int*)Vmem_malloc(thee->vmem, thee->natom, sizeof(int));
254  VASSERT(thee->nqsm != VNULL);
255  for (iatom=0; iatom<thee->natom; iatom++) (thee->nqsm)[iatom] = 0;
256  /* Loop through the list of simplices and count the number of times
257  * each atom appears */
258  for (isimp=0; isimp<thee->nsimp; isimp++) {
259  for (iatom=0; iatom<thee->nsqm[isimp]; iatom++) {
260  jatom = thee->sqm[isimp][iatom];
261  thee->nqsm[jatom]++;
262  }
263  }
264  /* Do a TIME-CONSUMING SANITY CHECK to make sure that each atom was
265  * placed in at simplex */
266  for (iatom=0; iatom<thee->natom; iatom++) {
267  if (thee->nqsm[iatom] == 0) {
268  Vnm_print(2, "Vcsm_init: Atom %d not placed in simplex!\n", iatom);
269  VASSERT(0);
270  }
271  }
272  /* Allocate the appropriate amount of space for each entry in the
273  * charge-simplex map and clear the counter for re-use in assignment */
274  for (iatom=0; iatom<thee->natom; iatom++) {
275  thee->qsm[iatom] = (int*)Vmem_malloc(thee->vmem, (thee->nqsm)[iatom],
276  sizeof(int));
277  VASSERT(thee->qsm[iatom] != VNULL);
278  thee->nqsm[iatom] = 0;
279  }
280  /* Assign the simplices to atoms */
281  for (isimp=0; isimp<thee->nsimp; isimp++) {
282  for (iatom=0; iatom<thee->nsqm[isimp]; iatom++) {
283  jatom = thee->sqm[isimp][iatom];
284  thee->qsm[jatom][thee->nqsm[jatom]] = isimp;
285  thee->nqsm[jatom]++;
286  }
287  }
288 
289  thee->initFlag = 1;
290 }
291 
292 VPUBLIC void Vcsm_dtor(Vcsm **thee) {
293  if ((*thee) != VNULL) {
294  Vcsm_dtor2(*thee);
295  Vmem_free(VNULL, 1, sizeof(Vcsm), (void **)thee);
296  (*thee) = VNULL;
297  }
298 }
299 
300 VPUBLIC void Vcsm_dtor2(Vcsm *thee) {
301  int i;
302 
303  if ((thee != VNULL) && thee->initFlag) {
304 
305  for (i=0; i<thee->msimp; i++) {
306  if (thee->nsqm[i] > 0) Vmem_free(thee->vmem, thee->nsqm[i],
307  sizeof(int), (void **)&(thee->sqm[i]));
308  }
309  for (i=0; i<thee->natom; i++) {
310  if (thee->nqsm[i] > 0) Vmem_free(thee->vmem, thee->nqsm[i],
311  sizeof(int), (void **)&(thee->qsm[i]));
312  }
313  Vmem_free(thee->vmem, thee->msimp, sizeof(int *),
314  (void **)&(thee->sqm));
315  Vmem_free(thee->vmem, thee->msimp, sizeof(int),
316  (void **)&(thee->nsqm));
317  Vmem_free(thee->vmem, thee->natom, sizeof(int *),
318  (void **)&(thee->qsm));
319  Vmem_free(thee->vmem, thee->natom, sizeof(int),
320  (void **)&(thee->nqsm));
321 
322  }
323  Vmem_dtor(&(thee->vmem));
324 }
325 
326 VPUBLIC int Vcsm_update(Vcsm *thee, SS **simps, int num) {
327 
328  /* Counters */
329  int isimp, jsimp, iatom, jatom, atomID, simpID;
330  int nsimps, gotMem;
331  /* Object info */
332  Vatom *atom;
333  SS *simplex;
334  double *position;
335  /* Lists */
336  int *qParent, nqParent;
337  int **sqmNew, *nsqmNew;
338  int *affAtoms, nAffAtoms;
339  int *dnqsm, *nqsmNew, **qsmNew;
340 
341  VASSERT(thee != VNULL);
342  VASSERT(thee->initFlag);
343 
344  /* If we don't have enough memory to accommodate the new entries,
345  * add more by doubling the existing amount */
346  isimp = thee->nsimp + num - 1;
347  gotMem = 0;
348  while (!gotMem) {
349  if (isimp > thee->msimp) {
350  isimp = 2 * isimp;
351  thee->nsqm = (int*)Vmem_realloc(thee->vmem, thee->msimp, sizeof(int),
352  (void **)&(thee->nsqm), isimp);
353  VASSERT(thee->nsqm != VNULL);
354  thee->sqm = (int**)Vmem_realloc(thee->vmem, thee->msimp, sizeof(int *),
355  (void **)&(thee->sqm), isimp);
356  VASSERT(thee->sqm != VNULL);
357  thee->msimp = isimp;
358  } else gotMem = 1;
359  }
360  /* Initialize the nsqm entires we just allocated */
361  for (isimp = thee->nsimp; isimp<thee->nsimp+num-1 ; isimp++) {
362  thee->nsqm[isimp] = 0;
363  }
364 
365  thee->nsimp = thee->nsimp + num - 1;
366 
367  /* There's a simple case to deal with: if simps[0] didn't have a
368  * charge in the first place */
369  isimp = SS_id(simps[0]);
370  if (thee->nsqm[isimp] == 0) {
371  for (isimp=1; isimp<num; isimp++) {
372  thee->nsqm[SS_id(simps[isimp])] = 0;
373  }
374  return 1;
375  }
376 
377  /* The more complicated case has occured; the parent simplex had one or
378  * more charges. First, generate the list of affected charges. */
379  isimp = SS_id(simps[0]);
380  nqParent = thee->nsqm[isimp];
381  qParent = thee->sqm[isimp];
382 
383  sqmNew = (int**)Vmem_malloc(thee->vmem, num, sizeof(int *));
384  VASSERT(sqmNew != VNULL);
385  nsqmNew = (int*)Vmem_malloc(thee->vmem, num, sizeof(int));
386  VASSERT(nsqmNew != VNULL);
387  for (isimp=0; isimp<num; isimp++) nsqmNew[isimp] = 0;
388 
389  /* Loop throught the affected atoms to determine how many atoms each
390  * simplex will get. */
391  for (iatom=0; iatom<nqParent; iatom++) {
392 
393  atomID = qParent[iatom];
394  atom = Valist_getAtom(thee->alist, atomID);
395  position = Vatom_getPosition(atom);
396  nsimps = 0;
397 
398  jsimp = 0;
399 
400  for (isimp=0; isimp<num; isimp++) {
401  simplex = simps[isimp];
402  if (Gem_pointInSimplex(thee->gm, simplex, position)) {
403  nsqmNew[isimp]++;
404  jsimp = 1;
405  }
406  }
407 
408  VASSERT(jsimp != 0);
409  }
410 
411  /* Sanity check that we didn't lose any atoms... */
412  iatom = 0;
413  for (isimp=0; isimp<num; isimp++) iatom += nsqmNew[isimp];
414  if (iatom < nqParent) {
415  Vnm_print(2,"Vcsm_update: Lost %d (of %d) atoms!\n",
416  nqParent - iatom, nqParent);
417  VASSERT(0);
418  }
419 
420  /* Allocate the storage */
421  for (isimp=0; isimp<num; isimp++) {
422  if (nsqmNew[isimp] > 0) {
423  sqmNew[isimp] = (int*)Vmem_malloc(thee->vmem, nsqmNew[isimp],
424  sizeof(int));
425  VASSERT(sqmNew[isimp] != VNULL);
426  }
427  }
428 
429  /* Assign charges to simplices */
430  for (isimp=0; isimp<num; isimp++) {
431 
432  jsimp = 0;
433  simplex = simps[isimp];
434 
435  /* Loop over the atoms associated with the parent simplex */
436  for (iatom=0; iatom<nqParent; iatom++) {
437 
438  atomID = qParent[iatom];
439  atom = Valist_getAtom(thee->alist, atomID);
440  position = Vatom_getPosition(atom);
441  if (Gem_pointInSimplex(thee->gm, simplex, position)) {
442  sqmNew[isimp][jsimp] = atomID;
443  jsimp++;
444  }
445  }
446  }
447 
448  /* Update the QSM map using the old and new SQM lists */
449  /* The affected atoms are those contained in the parent simplex; i.e.
450  * thee->sqm[SS_id(simps[0])] */
451  affAtoms = thee->sqm[SS_id(simps[0])];
452  nAffAtoms = thee->nsqm[SS_id(simps[0])];
453  /* Each of these atoms will go somewhere else; i.e., the entries in
454  * thee->qsm are never destroyed and thee->nqsm never decreases.
455  * However, it is possible that a subdivision could cause an atom to be
456  * shared by two child simplices. Here we record the change, if any,
457  * in the number of simplices associated with each atom. */
458  dnqsm = (int*)Vmem_malloc(thee->vmem, nAffAtoms, sizeof(int));
459  VASSERT(dnqsm != VNULL);
460  nqsmNew = (int*)Vmem_malloc(thee->vmem, nAffAtoms, sizeof(int));
461  VASSERT(nqsmNew != VNULL);
462  qsmNew = (int**)Vmem_malloc(thee->vmem, nAffAtoms, sizeof(int*));
463  VASSERT(qsmNew != VNULL);
464  for (iatom=0; iatom<nAffAtoms; iatom++) {
465  dnqsm[iatom] = -1;
466  atomID = affAtoms[iatom];
467  for (isimp=0; isimp<num; isimp++) {
468  for (jatom=0; jatom<nsqmNew[isimp]; jatom++) {
469  if (sqmNew[isimp][jatom] == atomID) dnqsm[iatom]++;
470  }
471  }
472  VASSERT(dnqsm[iatom] > -1);
473  }
474  /* Setup the new entries in the array */
475  for (iatom=0;iatom<nAffAtoms; iatom++) {
476  atomID = affAtoms[iatom];
477  qsmNew[iatom] = (int*)Vmem_malloc(thee->vmem,
478  (dnqsm[iatom] + thee->nqsm[atomID]),
479  sizeof(int));
480  nqsmNew[iatom] = 0;
481  VASSERT(qsmNew[iatom] != VNULL);
482  }
483  /* Fill the new entries in the array */
484  /* First, do the modified entries */
485  for (isimp=0; isimp<num; isimp++) {
486  simpID = SS_id(simps[isimp]);
487  for (iatom=0; iatom<nsqmNew[isimp]; iatom++) {
488  atomID = sqmNew[isimp][iatom];
489  for (jatom=0; jatom<nAffAtoms; jatom++) {
490  if (atomID == affAtoms[jatom]) break;
491  }
492  if (jatom < nAffAtoms) {
493  qsmNew[jatom][nqsmNew[jatom]] = simpID;
494  nqsmNew[jatom]++;
495  }
496  }
497  }
498  /* Now do the unmodified entries */
499  for (iatom=0; iatom<nAffAtoms; iatom++) {
500  atomID = affAtoms[iatom];
501  for (isimp=0; isimp<thee->nqsm[atomID]; isimp++) {
502  for (jsimp=0; jsimp<num; jsimp++) {
503  simpID = SS_id(simps[jsimp]);
504  if (thee->qsm[atomID][isimp] == simpID) break;
505  }
506  if (jsimp == num) {
507  qsmNew[iatom][nqsmNew[iatom]] = thee->qsm[atomID][isimp];
508  nqsmNew[iatom]++;
509  }
510  }
511  }
512 
513  /* Replace the existing entries in the table. Do the QSM entires
514  * first, since they require affAtoms = thee->sqm[simps[0]] */
515  for (iatom=0; iatom<nAffAtoms; iatom++) {
516  atomID = affAtoms[iatom];
517  Vmem_free(thee->vmem, thee->nqsm[atomID], sizeof(int),
518  (void **)&(thee->qsm[atomID]));
519  thee->qsm[atomID] = qsmNew[iatom];
520  thee->nqsm[atomID] = nqsmNew[iatom];
521  }
522  for (isimp=0; isimp<num; isimp++) {
523  simpID = SS_id(simps[isimp]);
524  if (thee->nsqm[simpID] > 0) Vmem_free(thee->vmem, thee->nsqm[simpID],
525  sizeof(int), (void **)&(thee->sqm[simpID]));
526  thee->sqm[simpID] = sqmNew[isimp];
527  thee->nsqm[simpID] = nsqmNew[isimp];
528  }
529 
530  Vmem_free(thee->vmem, num, sizeof(int *), (void **)&sqmNew);
531  Vmem_free(thee->vmem, num, sizeof(int), (void **)&nsqmNew);
532  Vmem_free(thee->vmem, nAffAtoms, sizeof(int *), (void **)&qsmNew);
533  Vmem_free(thee->vmem, nAffAtoms, sizeof(int), (void **)&nqsmNew);
534  Vmem_free(thee->vmem, nAffAtoms, sizeof(int), (void **)&dnqsm);
535 
536 
537  return 1;
538 
539 
540 }
int * nqsm
Definition: vcsm.h:111
Contains declarations for the Vcsm class.
VPUBLIC int Vcsm_getSimplexIndex(Vcsm *thee, int isimp, int iatom)
Get index particular simplex associated with an atom.
Definition: vcsm.c:119
VPUBLIC Vcsm * Vcsm_ctor(Valist *alist, Gem *gm)
Construct Vcsm object.
Definition: vcsm.c:136
VPUBLIC int Vcsm_update(Vcsm *thee, SS **simps, int num)
Update the charge-simplex and simplex-charge maps after refinement.
Definition: vcsm.c:326
Valist * alist
Definition: vcsm.h:91
VPUBLIC int Vcsm_ctor2(Vcsm *thee, Valist *alist, Gem *gm)
FORTRAN stub to construct Vcsm object.
Definition: vcsm.c:147
VPUBLIC int Vcsm_getNumberSimplices(Vcsm *thee, int iatom)
Get number of simplices associated with an atom.
Definition: vcsm.c:99
int initFlag
Definition: vcsm.h:112
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition: valist.c:105
int ** qsm
Definition: vcsm.h:109
VPUBLIC SS * Vcsm_getSimplex(Vcsm *thee, int isimp, int iatom)
Get particular simplex associated with an atom.
Definition: vcsm.c:109
VPUBLIC Vatom * Vcsm_getAtom(Vcsm *thee, int iatom, int isimp)
Get particular atom associated with a simplex.
Definition: vcsm.c:77
VPUBLIC Valist * Vcsm_getValist(Vcsm *thee)
Get atom list.
Definition: vcsm.c:62
VPUBLIC void Vcsm_dtor(Vcsm **thee)
Destroy Vcsm object.
Definition: vcsm.c:292
Charge-simplex map class.
Definition: vcsm.h:89
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition: valist.c:115
int natom
Definition: vcsm.h:92
int ** sqm
Definition: vcsm.h:97
int msimp
Definition: vcsm.h:107
VPUBLIC int Vcsm_getAtomIndex(Vcsm *thee, int iatom, int isimp)
Get ID of particular atom in a simplex.
Definition: vcsm.c:88
int * nsqm
Definition: vcsm.h:104
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition: vatom.c:63
VPUBLIC unsigned long int Vcsm_memChk(Vcsm *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition: vcsm.c:129
Contains public data members for Vatom class/module.
Definition: vatom.h:84
Container class for list of atom objects.
Definition: valist.h:78
int nsimp
Definition: vcsm.h:105
VPUBLIC int Vcsm_getNumberAtoms(Vcsm *thee, int isimp)
Get number of atoms associated with a simplex.
Definition: vcsm.c:69
Vmem * vmem
Definition: vcsm.h:114
VPUBLIC void Vcsm_init(Vcsm *thee)
Initialize charge-simplex map with mesh and atom data.
Definition: vcsm.c:170
VPUBLIC void Vcsm_dtor2(Vcsm *thee)
FORTRAN stub to destroy Vcsm object.
Definition: vcsm.c:300
Gem * gm
Definition: vcsm.h:94