APBS 3.0.0
Loading...
Searching...
No Matches
vcsm.c
Go to the documentation of this file.
1
57#include "vcsm.h"
58
59/* Inlineable methods */
60#if !defined(VINLINE_VCSM)
61
62VPUBLIC Valist* Vcsm_getValist(Vcsm *thee) {
63
64 VASSERT(thee != VNULL);
65 return thee->alist;
66
67}
68
69VPUBLIC int Vcsm_getNumberAtoms(Vcsm *thee, int isimp) {
70
71 VASSERT(thee != VNULL);
72 VASSERT(thee->initFlag);
73 return thee->nsqm[isimp];
74
75}
76
77VPUBLIC 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
88VPUBLIC 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
99VPUBLIC 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
109VPUBLIC 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
119VPUBLIC 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
129VPUBLIC 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
136VPUBLIC 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
147VPUBLIC 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
170VPUBLIC 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
292VPUBLIC 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
300VPUBLIC 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
326VPUBLIC 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}
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition valist.c:115
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition valist.c:105
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition vatom.c:63
VPUBLIC int Vcsm_getNumberAtoms(Vcsm *thee, int isimp)
Get number of atoms associated with a simplex.
Definition vcsm.c:69
VPUBLIC int Vcsm_getAtomIndex(Vcsm *thee, int iatom, int isimp)
Get ID of particular atom in a simplex.
Definition vcsm.c:88
VPUBLIC void Vcsm_dtor2(Vcsm *thee)
FORTRAN stub to destroy Vcsm object.
Definition vcsm.c:300
VPUBLIC Vcsm * Vcsm_ctor(Valist *alist, Gem *gm)
Construct Vcsm object.
Definition vcsm.c:136
VPUBLIC SS * Vcsm_getSimplex(Vcsm *thee, int isimp, int iatom)
Get particular simplex associated with an atom.
Definition vcsm.c:109
VPUBLIC int Vcsm_update(Vcsm *thee, SS **simps, int num)
Update the charge-simplex and simplex-charge maps after refinement.
Definition vcsm.c:326
VPUBLIC Vatom * Vcsm_getAtom(Vcsm *thee, int iatom, int isimp)
Get particular atom associated with a simplex.
Definition vcsm.c:77
VPUBLIC int Vcsm_ctor2(Vcsm *thee, Valist *alist, Gem *gm)
FORTRAN stub to construct Vcsm object.
Definition vcsm.c:147
VPUBLIC void Vcsm_init(Vcsm *thee)
Initialize charge-simplex map with mesh and atom data.
Definition vcsm.c:170
VPUBLIC int Vcsm_getNumberSimplices(Vcsm *thee, int iatom)
Get number of simplices associated with an atom.
Definition vcsm.c:99
VPUBLIC unsigned long int Vcsm_memChk(Vcsm *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition vcsm.c:129
VPUBLIC void Vcsm_dtor(Vcsm **thee)
Destroy Vcsm object.
Definition vcsm.c:292
VPUBLIC Valist * Vcsm_getValist(Vcsm *thee)
Get atom list.
Definition vcsm.c:62
VPUBLIC int Vcsm_getSimplexIndex(Vcsm *thee, int isimp, int iatom)
Get index particular simplex associated with an atom.
Definition vcsm.c:119
Container class for list of atom objects.
Definition valist.h:78
Contains public data members for Vatom class/module.
Definition vatom.h:84
Charge-simplex map class.
Definition vcsm.h:89
int ** sqm
Definition vcsm.h:97
int initFlag
Definition vcsm.h:112
Gem * gm
Definition vcsm.h:94
int nsimp
Definition vcsm.h:105
int * nqsm
Definition vcsm.h:111
int natom
Definition vcsm.h:92
int msimp
Definition vcsm.h:107
int * nsqm
Definition vcsm.h:104
int ** qsm
Definition vcsm.h:109
Vmem * vmem
Definition vcsm.h:114
Valist * alist
Definition vcsm.h:91
Contains declarations for the Vcsm class.