1 /* Steven Andrews, started 10/22/2001.
2 This is a library of functions for the Smoldyn program.
3 See documentation called SmoldynManual.pdf and SmoldynCodeDoc.pdf, and the Smoldyn
4 website, which is at www.smoldyn.org.
5 Copyright 2003-2016 by Steven Andrews. This work is distributed under the terms
6 of the Gnu Lesser General Public License (LGPL). */
7
8
9 #include <stdio.h>
10 #include <math.h>
11 #include <string.h>
12 #include "random2.h"
13 #include "RnSort.h"
14 #include "string2.h"
15 #include "Zn.h"
16
17 #include "smoldyn.h"
18 #include "smoldynfuncs.h"
19
20 #ifdef OPTION_VCELL
21 #include <algorithm>
22 using std::max;
23 using std::min;
24 #include <fstream>
25 #include <ostream>
26 #include <sstream>
27 #include <iostream>
28
29 extern "C" {
30 #include "zlib.h" }
31 #endif
32
33
34 int compartsupdateparams_original(simptr sim);
35 int compartsupdateparams_volumeSample(simptr sim);
36
37
38 /******************************************************************************/
39 /********************************* Compartments *******************************/
40 /******************************************************************************/
41
42
43 /******************************************************************************/
44 /****************************** Local declarations ****************************/
45 /******************************************************************************/
46
47 // enumerated types
48 enum CmptLogic compartstring2cl(char *string);
49 char *compartcl2string(enum CmptLogic cls,char *string);
50
51 // low level utilities
52
53 // memory management
54 compartptr compartalloc(void);
55 void compartfree(compartptr cmpt);
56 compartssptr compartssalloc(compartssptr cmptss,int maxcmpt);
57
58 // data structure output
59
60 // structure set up
61 int compartupdatebox(simptr sim,compartptr cmpt,boxptr bptr,double volfrac);
62 int compartsupdateparams(simptr sim);
63 int compartsupdatelists(simptr sim);
64
65 /******************************************************************************/
66 /********************************* enumerated types ***************************/
67 /******************************************************************************/
68
69
70 /* compartstring2cl */
compartstring2cl(char * string)71 enum CmptLogic compartstring2cl(char *string) {
72 enum CmptLogic ans;
73
74 if(!strcmp(string,"equal")) ans=CLequal;
75 else if(!strcmp(string,"equalnot")) ans=CLequalnot;
76 else if(!strcmp(string,"and")) ans=CLand;
77 else if(!strcmp(string,"or")) ans=CLor;
78 else if(!strcmp(string,"xor")) ans=CLxor;
79 else if(!strcmp(string,"andnot")) ans=CLandnot;
80 else if(!strcmp(string,"ornot")) ans=CLornot;
81 else ans=CLnone;
82 return ans; }
83
84
85 /* compartcl2string */
compartcl2string(enum CmptLogic cls,char * string)86 char *compartcl2string(enum CmptLogic cls,char *string) {
87 if(cls==CLequal) strcpy(string,"equal");
88 else if(cls==CLequalnot) strcpy(string,"equalnot");
89 else if(cls==CLand) strcpy(string,"and");
90 else if(cls==CLor) strcpy(string,"or");
91 else if(cls==CLxor) strcpy(string,"xor");
92 else if(cls==CLandnot) strcpy(string,"andnot");
93 else if(cls==CLornot) strcpy(string,"ornot");
94 else strcpy(string,"none");
95 return string; }
96
97
98 /******************************************************************************/
99 /****************************** low level utilities ***************************/
100 /******************************************************************************/
101
102 #ifdef OPTION_VCELL
103 //called by posincompart and compartsupdateparams_volumeSample
getCompartmentID(char * cmptName,VolumeSamplesPtr vSamplesPtr)104 unsigned char getCompartmentID(char* cmptName, VolumeSamplesPtr vSamplesPtr)
105 {
106 int c;
107 for(c=0;c<vSamplesPtr->nCmptIDPair;c++) {
108 if(!strcmp(cmptName, vSamplesPtr->compartmentIDPairPtr[c].name))
109 {
110 return vSamplesPtr->compartmentIDPairPtr[c].pixel;
111 }
112 }
113 return 0; //TODO: if can't find the compartment ID, what should we do? ??
114 }
115 #endif
116
117
118 /* posincompart */
posincompart(simptr sim,double * pos,compartptr cmpt,int useoldpos)119 int posincompart(simptr sim,double *pos,compartptr cmpt,int useoldpos) {
120 int s,p,k,incmpt,pcross,cl,incmptl;
121 enum PanelShape ps;
122 surfaceptr srf;
123 double crsspt[DIMMAX];
124 enum CmptLogic sym;
125
126 #ifdef OPTION_VCELL
127 int dim;
128 VolumeSamples* volumeSamplePtr;
129 const int numNeighbors =3;
130 int sampleIdxNeighbors[DIMMAX][numNeighbors];
131
132 incmpt = 0;
133 int sameResultCount = 0;
134
135 //for smoldyn simulations with volumeSample
136 if(sim->volumeSamplesPtr != NULL) {
137 volumeSamplePtr = sim->volumeSamplesPtr;
138 dim = sim->dim;
139 unsigned char cmptID = getCompartmentID(cmpt->cname, sim->volumeSamplesPtr);
140
141 for (int idim = 0; idim < dim; idim ++) {
142 double h = volumeSamplePtr->size[idim]/volumeSamplePtr->num[idim];
143 int idx = (int)floor((pos[idim] - volumeSamplePtr->origin[idim])/h);
144 idx = max(0,idx);//clamp index lower bound to 0
145 idx = min((volumeSamplePtr->num[idim]-1),idx);//clamp index upper bound to number of volume sample at specific dimension
146 sampleIdxNeighbors[idim][0] = std::max<int>(0,(idx-1));
147 sampleIdxNeighbors[idim][1] = idx;
148 sampleIdxNeighbors[idim][2] = std::min<int>((volumeSamplePtr->num[idim]-1),(idx+1)); }
149
150 //int numNeighbors =3;
151 for (int kk = 0; kk < (dim < 3 ? 1 : numNeighbors); kk ++) {
152 for (int jj = 0; jj < (dim < 2 ? 1 : numNeighbors); jj ++) {
153 for (int ii = 0; ii < numNeighbors; ii ++) {
154 int i = sampleIdxNeighbors[0][ii];
155 int j = dim < 2 ? 0 : sampleIdxNeighbors[1][jj];
156 int k = dim < 3 ? 0 : sampleIdxNeighbors[2][kk];
157 int posInSample = i + j*volumeSamplePtr->num[0] + k*volumeSamplePtr->num[0]*volumeSamplePtr->num[1];
158 unsigned char cmptIDFoundInSample = volumeSamplePtr->volsamples[posInSample];
159 if(cmptID == cmptIDFoundInSample){
160 sameResultCount ++; }}}}
161 if (sameResultCount == (int)pow((double)numNeighbors, dim)) {
162 incmpt = 1; }}
163
164 if((incmpt == 0 && sameResultCount > 0) || sim->volumeSamplesPtr == NULL) //for original smoldyn simulations(without volumeSample)
165 #endif
166
167 {
168 incmpt=0;
169 for(k=0;k<cmpt->npts && incmpt==0;k++) {
170 pcross=0;
171 for(s=0;s<cmpt->nsrf && !pcross;s++) {
172 srf=cmpt->surflist[s];
173 for(ps=(enum PanelShape)0;ps<PSMAX&&!pcross;ps=(enum PanelShape)(ps+1))
174 for(p=0;p<srf->npanel[ps]&&!pcross;p++)
175 if(lineXpanel(pos,cmpt->points[k],srf->panels[ps][p],sim->dim,crsspt,NULL,NULL,NULL,NULL,NULL,useoldpos))
176 pcross=1; }
177 if(pcross==0) incmpt=1; }
178
179 for(cl=0;cl<cmpt->ncmptl;cl++) {
180 incmptl=posincompart(sim,pos,cmpt->cmptl[cl],0);
181 sym=cmpt->clsym[cl];
182 if(sym==CLequal) incmpt=incmptl;
183 else if(sym==CLequalnot) incmpt=!incmptl;
184 else if(sym==CLand) incmpt=incmpt&&incmptl;
185 else if(sym==CLor) incmpt=incmpt||incmptl;
186 else if(sym==CLxor) incmpt=(incmpt!=incmptl);
187 else if(sym==CLandnot) incmpt=incmpt&&!incmptl;
188 else if(sym==CLornot) incmpt=incmpt||!incmptl; }
189 }
190 return incmpt; }
191
192
193 /* compartrandpos */
compartrandpos(simptr sim,double * pos,compartptr cmpt)194 int compartrandpos(simptr sim,double *pos,compartptr cmpt) {
195 static int ptmax=10000;
196 int d,dim,i,done,k,bc;
197 boxptr bptr;
198
199 if(cmpt->npts==0&&cmpt->ncmptl==0) return 1;
200 dim=sim->dim;
201
202 done=0;
203 if(cmpt->nbox) {
204 bc=intrandpD(cmpt->nbox,cmpt->cumboxvol);
205 bptr=cmpt->boxlist[bc];
206 for(i=0;i<ptmax&&!done;i++) {
207 boxrandpos(sim,pos,bptr);
208 if(posincompart(sim,pos,cmpt,0)) done=1; }}
209 else {
210 for(i=0;i<ptmax&&!done;i++) {
211 for(d=0;d<dim;d++) pos[d]=unirandCCD(sim->wlist[2*d]->pos,sim->wlist[2*d+1]->pos);
212 if(posincompart(sim,pos,cmpt,0)) done=1; }}
213 if(!done&&cmpt->npts>0) {
214 k=intrand(cmpt->npts);
215 for(d=0;d<sim->dim;d++) pos[d]=cmpt->points[k][d];
216 done=1; }
217 if(!done) return 1;
218 return 0; }
219
220
221 /* fromHex */
fromHex(const char * src)222 unsigned char fromHex(const char* src) {
223 char chs[5];
224 chs[0] = '0';
225 chs[1] = 'x';
226 chs[2] = src[0];
227 chs[3] = src[1];
228 chs[4] = 0;
229 int v;
230 sscanf(chs, "%x", &v);
231 return (unsigned char)v;
232 }
233
234
235 #ifdef OPTION_VCELL
236 using std::istringstream;
237 /* loadHighResVolumeSamples */
loadHighResVolumeSamples(simptr sim,ParseFilePtr * pfpptr,char * line2)238 int loadHighResVolumeSamples(simptr sim,ParseFilePtr *pfpptr,char *line2) { //?? needs to be added to user manual
239 char errstring[STRCHAR];
240
241 char word[STRCHAR];
242 ParseFilePtr pfp = *pfpptr;
243 int pfpcode;
244 *pfpptr=pfp;
245
246 VolumeSamplesPtr volumeSamplesPtr = new VolumeSamples;
247 sim->volumeSamplesPtr = volumeSamplesPtr;
248 volumeSamplesPtr->num[0] = volumeSamplesPtr->num[1] = volumeSamplesPtr->num[2] = 1;
249
250 int done = 0;
251 string pixelLine;
252 while(!done) {
253 pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
254
255 *pfpptr=pfp;
256 CHECKS(pfpcode!=3,"%s",errstring);
257
258 if(pfpcode==0); // already taken care of
259
260 else if(pfpcode==2) { // end reading
261 done = 1;
262 } else if(pfpcode==3) { // error
263 CHECKS(0,"SMOLDYN BUG: parsing error");
264 } else if(!strcmp(word,"end_highResVolumeSamples")) { // end_jms
265 CHECKS(!line2,"unexpected text following end_highResVolumeSamplesFile");
266 break;
267 } else {
268 if (!strcmp(word, "Size")) {
269 istringstream lineInput(line2);
270 lineInput >> volumeSamplesPtr->size[0] >> volumeSamplesPtr->size[1] >> volumeSamplesPtr->size[2];
271 } else if (!strcmp(word, "Origin")) {
272 istringstream lineInput(line2);
273 lineInput >> volumeSamplesPtr->origin[0] >> volumeSamplesPtr->origin[1] >> volumeSamplesPtr->origin[2];
274 } else if (!strcmp(word, "CompartmentHighResPixelMap")) {
275 istringstream lineInput(line2);
276 lineInput >> volumeSamplesPtr->nCmptIDPair;
277 volumeSamplesPtr->compartmentIDPairPtr = new CompartmentIdentifierPair[volumeSamplesPtr->nCmptIDPair];
278 for (int i = 0; i < volumeSamplesPtr->nCmptIDPair; i ++) {
279 pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
280 strcpy(volumeSamplesPtr->compartmentIDPairPtr[i].name, word);
281 istringstream lineInput1(line2);
282 int pixel;
283 lineInput1 >> pixel;
284 volumeSamplesPtr->compartmentIDPairPtr[i].pixel = pixel;
285 }
286 } else if (!strcmp(word, "VolumeSamples")) {
287 istringstream lineInput(line2);
288 lineInput >> volumeSamplesPtr->num[0] >> volumeSamplesPtr->num[1] >> volumeSamplesPtr->num[2];
289 } else {//volume sample data with no token in front
290 pixelLine += word;
291 }
292 }
293 }
294
295 #ifdef HAVE_ZLIB
296 unsigned long inflated_len;
297 const char* compressed_hex;
298 unsigned char* bytes_from_compressed;
299 int compressed_len, returnVal;
300 long numVolumeSamples;
301
302 numVolumeSamples = volumeSamplesPtr->num[0] * volumeSamplesPtr->num[1] * volumeSamplesPtr->num[2];
303 //cout<<pixelLine<<endl;
304 compressed_len = pixelLine.size();
305 if (compressed_len <= 1) {
306 throw "CartesianMesh::readGeometryFile() : invalid compressed volume";
307 }
308
309 compressed_hex = pixelLine.c_str();
310 //volumeSamples compressed, changed from byte to short
311 bytes_from_compressed = new unsigned char[compressed_len+1];
312 memset(bytes_from_compressed, 0, (compressed_len+1) * sizeof(unsigned char));
313 for (int i = 0, j = 0; i < compressed_len; i += 2, j ++) {
314 bytes_from_compressed[j] = fromHex(compressed_hex + i);
315 }
316
317 volumeSamplesPtr->volsamples = new unsigned char[numVolumeSamples];
318 memset(volumeSamplesPtr->volsamples, 0, numVolumeSamples * sizeof(unsigned char));
319
320 inflated_len = numVolumeSamples;
321 returnVal = uncompress(volumeSamplesPtr->volsamples, &inflated_len, bytes_from_compressed, compressed_len);
322 //check if the data can be properly uncompressed, Z_OK should be returned if successfully uncompressed.
323 if(returnVal == Z_MEM_ERROR)
324 {
325 throw "There is not enough memory to uncompress volume samples.";
326 }
327 else if(returnVal == Z_BUF_ERROR)
328 {
329 throw "There is not enough room in the output buffer for volume samples.";
330 }
331 else if(returnVal == Z_DATA_ERROR)
332 {
333 throw "The volume sample input data was corrupted.";
334 }
335 //check if we get the expected number of volume samples after uncompression
336 if ((long)inflated_len != numVolumeSamples) {
337 throw "loadHighResVolumeSamples : unexpected number of volume samples";
338 }
339 //for debug purpose
340 /*for (unsigned long i = 0; i < inflated_len; i ++) {
341 if (volumeSamplesPtr->volsamples[i] == 6) {
342 cout << "volume sample at " << i<< " is " << ((int)volumeSamplesPtr->volsamples[i])<< endl;
343 }
344 else if (volumeSamplesPtr->volsamples[i] == 16) {
345 cout << "volume sample at " << i<< " is " << ((int)volumeSamplesPtr->volsamples[i])<< endl;
346 }
347 }*/
348 return 0;
349
350 #else
351 throw "loadHighResVolumeSamples : function unavailable because compile does not include zlib";
352 #endif
353
354 failure: // failure
355 return 1;
356 }
357 #endif
358
359
360 /******************************************************************************/
361 /******************************* memory management ****************************/
362 /******************************************************************************/
363
364 /* compartalloc */
compartalloc(void)365 compartptr compartalloc(void) {
366 compartptr cmpt;
367
368 cmpt=(compartptr) malloc(sizeof(struct compartstruct));
369 CHECKMEM(cmpt);
370 cmpt->cname=NULL;
371 cmpt->selfindex=-1;
372 cmpt->nsrf=0;
373 cmpt->surflist=NULL;
374 cmpt->npts=0;
375 cmpt->points=NULL;
376 cmpt->ncmptl=0;
377 cmpt->cmptl=NULL;
378 cmpt->clsym=NULL;
379 cmpt->volume=0;
380 cmpt->maxbox=0;
381 cmpt->nbox=0;
382 cmpt->boxlist=NULL;
383 cmpt->boxfrac=NULL;
384 cmpt->cumboxvol=NULL;
385 return cmpt;
386 failure:
387 simLog(NULL,10,"Failed to allocate memory in compartalloc");
388 return NULL; }
389
390
391 /* compartfree */
compartfree(compartptr cmpt)392 void compartfree(compartptr cmpt) {
393 int k;
394
395 if(!cmpt) return;
396 free(cmpt->cumboxvol);
397 free(cmpt->boxfrac);
398 free(cmpt->boxlist);
399 free(cmpt->clsym);
400 free(cmpt->cmptl);
401 if(cmpt->npts&&cmpt->points)
402 for(k=0;k<cmpt->npts;k++) free(cmpt->points[k]);
403 free(cmpt->points);
404 free(cmpt->surflist);
405 free(cmpt);
406 return; }
407
408
409 /* compartssalloc */
compartssalloc(compartssptr cmptss,int maxcmpt)410 compartssptr compartssalloc(compartssptr cmptss,int maxcmpt) {
411 int c;
412 char **newnames;
413 compartptr *newcmptlist;
414
415 CHECKBUG(maxcmpt>0,"maxcmpt, in compartssalloc, needs to be >0");
416
417 newnames=NULL;
418 newcmptlist=NULL;
419
420 if(!cmptss) { // new allocation
421 cmptss=(compartssptr) malloc(sizeof(struct compartsuperstruct));
422 CHECKMEM(cmptss);
423 cmptss->condition=SCinit;
424 cmptss->sim=NULL;
425 cmptss->maxcmpt=0;
426 cmptss->ncmpt=0;
427 cmptss->cnames=NULL;
428 cmptss->cmptlist=NULL; }
429 else { // minor check
430 if(maxcmpt<cmptss->maxcmpt) return cmptss; }
431
432 if(maxcmpt>cmptss->maxcmpt) { // allocate new compartment names and compartments
433 CHECKMEM(newnames=(char**) calloc(maxcmpt,sizeof(char*)));
434 for(c=0;c<maxcmpt;c++) newnames[c]=NULL;
435 for(c=0;c<cmptss->maxcmpt;c++) newnames[c]=cmptss->cnames[c];
436 for(;c<maxcmpt;c++)
437 CHECKMEM(newnames[c]=EmptyString());
438
439 CHECKMEM(newcmptlist=(compartptr*) calloc(maxcmpt,sizeof(compartptr))); // compartment list
440 for(c=0;c<maxcmpt;c++) newcmptlist[c]=NULL;
441 for(c=0;c<cmptss->maxcmpt;c++) newcmptlist[c]=cmptss->cmptlist[c];
442 for(;c<maxcmpt;c++) {
443 CHECKMEM(newcmptlist[c]=compartalloc());
444 newcmptlist[c]->cmptss=cmptss;
445 newcmptlist[c]->cname=newnames[c];
446 newcmptlist[c]->selfindex=c; }}
447
448 cmptss->maxcmpt=maxcmpt;
449 free(cmptss->cnames);
450 cmptss->cnames=newnames;
451 free(cmptss->cmptlist);
452 cmptss->cmptlist=newcmptlist;
453
454 return cmptss;
455
456 failure:
457 compartssfree(cmptss);
458 simLog(NULL,10,"%s","Failed to allocated memory in compartssalloc");
459 return NULL; }
460
461
462 /* compartssfree */
compartssfree(compartssptr cmptss)463 void compartssfree(compartssptr cmptss) {
464 int c;
465
466 if(!cmptss) return;
467 if(cmptss->maxcmpt && cmptss->cmptlist)
468 for(c=0;c<cmptss->maxcmpt;c++) compartfree(cmptss->cmptlist[c]);
469 free(cmptss->cmptlist);
470 if(cmptss->maxcmpt && cmptss->cnames)
471 for(c=0;c<cmptss->maxcmpt;c++) free(cmptss->cnames[c]);
472 free(cmptss->cnames);
473 free(cmptss);
474 return; }
475
476
477 /******************************************************************************/
478 /***************************** data structure output **************************/
479 /******************************************************************************/
480
481 /* compartoutput */
compartoutput(simptr sim)482 void compartoutput(simptr sim) {
483 compartssptr cmptss;
484 compartptr cmpt;
485 int c,dim,s,k,d,cl;
486 char string[STRCHAR];
487
488 cmptss=sim->cmptss;
489 if(!cmptss) return;
490 simLog(sim,2,"COMPARTMENT PARAMETERS\n");
491 dim=sim->dim;
492 simLog(sim,2," Compartments allocated: %i, compartments defined: %i\n",cmptss->maxcmpt,cmptss->ncmpt);
493 for(c=0;c<cmptss->ncmpt;c++) {
494 cmpt=cmptss->cmptlist[c];
495 simLog(sim,2," Compartment: %s\n",cmptss->cnames[c]);
496 simLog(sim,2," %i bounding surfaces:\n",cmpt->nsrf);
497 for(s=0;s<cmpt->nsrf;s++)
498 simLog(sim,2," %s\n",cmpt->surflist[s]->sname);
499 simLog(sim,2," %i interior-defining points:\n",cmpt->npts);
500 for(k=0;k<cmpt->npts;k++) {
501 simLog(sim,2," %i: (",k);
502 for(d=0;d<dim-1;d++)
503 simLog(sim,2,"%g,",cmpt->points[k][d]);
504 simLog(sim,2,"%g)\n",cmpt->points[k][d]); }
505 simLog(sim,2," %i logically combined compartments\n",cmpt->ncmptl);
506 for(cl=0;cl<cmpt->ncmptl;cl++)
507 simLog(sim,2," %s %s\n",compartcl2string(cmpt->clsym[cl],string),cmpt->cmptl[cl]->cname);
508 simLog(sim,2," volume: %g\n",cmpt->volume);
509 simLog(sim,2," %i virtual boxes listed\n",cmpt->nbox); }
510 simLog(sim,2,"\n");
511 return; }
512
513
514 /* writecomparts */
writecomparts(simptr sim,FILE * fptr)515 void writecomparts(simptr sim,FILE *fptr) {
516 compartssptr cmptss;
517 compartptr cmpt;
518 int c,s,k,d,cl;
519 char string[STRCHAR];
520
521 cmptss=sim->cmptss;
522 if(!cmptss) return;
523 fprintf(fptr,"# Compartment parameters\n");
524 fprintf(fptr,"max_compartment %i\n",cmptss->maxcmpt);
525 for(c=0;c<cmptss->ncmpt;c++) {
526 cmpt=cmptss->cmptlist[c];
527 fprintf(fptr,"start_compartment %s\n",cmpt->cname);
528 for(s=0;s<cmpt->nsrf;s++)
529 fprintf(fptr,"surface %s\n",cmpt->surflist[s]->sname);
530 for(k=0;k<cmpt->npts;k++) {
531 fprintf(fptr,"point");
532 for(d=0;d<sim->dim;d++)
533 fprintf(fptr," %g",cmpt->points[k][d]);
534 fprintf(fptr,"\n"); }
535 for(cl=0;cl<cmpt->ncmptl;cl++)
536 fprintf(fptr,"compartment %s %s\n",compartcl2string(cmpt->clsym[cl],string),cmpt->cmptl[cl]->cname);
537 fprintf(fptr,"end_compartment\n\n"); }
538 return; }
539
540
541 /* checkcompartparams */
checkcompartparams(simptr sim,int * warnptr)542 int checkcompartparams(simptr sim,int *warnptr) {
543 int error,warn,c;
544 compartssptr cmptss;
545 compartptr cmpt;
546 char string[STRCHAR];
547
548 error=warn=0;
549 cmptss=sim->cmptss;
550 if(!cmptss) {
551 if(warnptr) *warnptr=warn;
552 return 0; }
553
554 if(cmptss->condition!=SCok) {
555 warn++;
556 simLog(sim,7," WARNING: compartment structure %s\n",simsc2string(cmptss->condition,string)); }
557
558 for(c=0;c<cmptss->ncmpt;c++) {
559 cmpt=cmptss->cmptlist[c];
560 if(cmpt->volume<=0) {warn++;simLog(sim,5," WARNING: compartment %s has 0 volume\n",cmpt->cname);}
561 if(cmpt->nbox==0) {warn++;simLog(sim,5," WARNING: compartment %s overlaps no virtual boxes\n",cmpt->cname);}
562 if(cmpt->nbox>0&&cmpt->cumboxvol[cmpt->nbox-1]!=cmpt->volume) {error++;simLog(sim,10," BUG: compartment %s box volumes do not add to compartment volume\n",cmpt->cname);} }
563 if(warnptr) *warnptr=warn;
564 return error; }
565
566
567 /******************************************************************************/
568 /******************************** structure set up ****************************/
569 /******************************************************************************/
570
571
572 /* compartsetcondition */
compartsetcondition(compartssptr cmptss,enum StructCond cond,int upgrade)573 void compartsetcondition(compartssptr cmptss,enum StructCond cond,int upgrade) {
574 if(!cmptss) return;
575 if(upgrade==0 && cmptss->condition>cond) cmptss->condition=cond;
576 else if(upgrade==1 && cmptss->condition<cond) cmptss->condition=cond;
577 else if(upgrade==2) cmptss->condition=cond;
578 if(cmptss->sim && cmptss->condition<cmptss->sim->condition) {
579 cond=cmptss->condition;
580 simsetcondition(cmptss->sim,cond==SCinit?SClists:cond,0); }
581 return; }
582
583
584 /* compartenablecomparts */
compartenablecomparts(simptr sim,int maxcmpt)585 int compartenablecomparts(simptr sim,int maxcmpt) {
586 compartssptr cmptss;
587
588 if(sim->cmptss) // check for redundant function call
589 if(maxcmpt==-1 || sim->cmptss->maxcmpt==maxcmpt)
590 return 0;
591 cmptss=compartssalloc(sim->cmptss,maxcmpt<0?5:maxcmpt);
592 if(!cmptss) return 1;
593 sim->cmptss=cmptss;
594 cmptss->sim=sim;
595 compartsetcondition(sim->cmptss,SClists,0);
596 return 0; }
597
598
599 /* compartaddcompart */
compartaddcompart(simptr sim,const char * cmptname)600 compartptr compartaddcompart(simptr sim,const char *cmptname) {
601 int er,c;
602 compartssptr cmptss;
603 compartptr cmpt;
604
605 if(!sim->cmptss) {
606 er=compartenablecomparts(sim,-1);
607 if(er) return NULL; }
608 cmptss=sim->cmptss;
609
610 c=stringfind(cmptss->cnames,cmptss->ncmpt,cmptname);
611 if(c<0) {
612 if(cmptss->ncmpt==cmptss->maxcmpt) {
613 er=compartenablecomparts(sim,cmptss->ncmpt*2+1);
614 if(er) return NULL; }
615 c=cmptss->ncmpt++;
616 strncpy(cmptss->cnames[c],cmptname,STRCHAR-1);
617 cmptss->cnames[c][STRCHAR-1]='\0';
618 cmpt=cmptss->cmptlist[c];
619 compartsetcondition(cmptss,SClists,0); }
620 else
621 cmpt=cmptss->cmptlist[c];
622
623 return cmpt; }
624
625
626 /* compartaddsurf */
compartaddsurf(compartptr cmpt,surfaceptr srf)627 int compartaddsurf(compartptr cmpt,surfaceptr srf) {
628 int s;
629 surfaceptr *newsurflist;
630
631 newsurflist=(surfaceptr*) calloc(cmpt->nsrf+1,sizeof(surfaceptr));
632 if(!newsurflist) return 1;
633 for(s=0;s<cmpt->nsrf;s++) {
634 if(cmpt->surflist[s]==srf) {free(newsurflist);return 2;}
635 newsurflist[s]=cmpt->surflist[s]; }
636 newsurflist[s]=srf;
637 cmpt->nsrf++;
638 free(cmpt->surflist);
639 cmpt->surflist=newsurflist;
640 cmpt->nbox=0;
641 cmpt->volume=0;
642 compartsetcondition(cmpt->cmptss,SCparams,0);
643 return 0; }
644
645
646 /* compartaddpoint */
compartaddpoint(compartptr cmpt,int dim,double * point)647 int compartaddpoint(compartptr cmpt,int dim,double *point) {
648 int d,k;
649 double **newpoints;
650
651 CHECKMEM(newpoints=(double**) calloc(cmpt->npts+1,sizeof(double*)));
652 for(k=0;k<cmpt->npts;k++)
653 newpoints[k]=cmpt->points[k];
654 CHECKMEM(newpoints[k]=(double*) calloc(dim,sizeof(double)));
655 for(d=0;d<dim;d++) newpoints[k][d]=point[d];
656 cmpt->npts++;
657 free(cmpt->points);
658 cmpt->points=newpoints;
659 compartsetcondition(cmpt->cmptss,SCparams,0);
660 cmpt->nbox=0;
661 cmpt->volume=0;
662 return 0;
663
664 failure:
665 if(newpoints) free(newpoints);
666 simLog(NULL,10,"Failed to allocate memory in compartaddpoint");
667 return 1; }
668
669
670 /* compartaddcmptl */
compartaddcmptl(compartptr cmpt,compartptr cmptl,enum CmptLogic sym)671 int compartaddcmptl(compartptr cmpt,compartptr cmptl,enum CmptLogic sym) {
672 int cl;
673 compartptr *newcmptl;
674 enum CmptLogic *newclsym;
675
676 if(cmpt==cmptl) return 2;
677 newcmptl=(compartptr*) calloc(cmpt->ncmptl+1,sizeof(compartptr));
678 if(!newcmptl) return 1;
679 newclsym=(enum CmptLogic*) calloc(cmpt->ncmptl+1,sizeof(enum CmptLogic));
680 if(!newclsym) {free(newcmptl);return 1;}
681 for(cl=0;cl<cmpt->ncmptl;cl++) {
682 newcmptl[cl]=cmpt->cmptl[cl];
683 newclsym[cl]=cmpt->clsym[cl]; }
684 newcmptl[cl]=cmptl;
685 newclsym[cl]=sym;
686 cmpt->ncmptl++;
687 free(cmpt->cmptl);
688 free(cmpt->clsym);
689 cmpt->cmptl=newcmptl;
690 cmpt->clsym=newclsym;
691 compartsetcondition(cmpt->cmptss,SCparams,0);
692 cmpt->nbox=0;
693 cmpt->volume=0;
694 return 0; }
695
696
697 /* compartupdatebox */
compartupdatebox(simptr sim,compartptr cmpt,boxptr bptr,double volfrac)698 int compartupdatebox(simptr sim,compartptr cmpt,boxptr bptr,double volfrac) {
699 int ptsmax=100; // number of random points for volume determination
700 int bc,max,ptsin,i,bc2;
701 double pos[DIMMAX],volfrac2,*newboxfrac,*newcumboxvol,boxvol,vol;
702 boxptr *newboxlist;
703
704 newboxlist=NULL;
705 newboxfrac=NULL;
706 newcumboxvol=NULL;
707
708 for(bc=0;bc<cmpt->nbox && cmpt->boxlist[bc]!=bptr;bc++); // check for box already in cmpt
709 if(bc<cmpt->nbox && volfrac==-2) return 0; // box is listed and volume ok, so return
710
711 if(volfrac<=0) { // find actual volume fraction
712 ptsin=0;
713 for(i=0;i<ptsmax;i++) {
714 boxrandpos(sim,pos,bptr);
715 if(posincompart(sim,pos,cmpt,0)) ptsin++; }
716 volfrac2=(double)ptsin/(double)ptsmax; }
717 else if(volfrac>1) volfrac2=1;
718 else volfrac2=volfrac;
719
720 if(volfrac2==0) {
721 if(bc==cmpt->nbox) return 0; // box not listed and 0 volume, so return
722 cmpt->nbox--; // box was listed, so it needs to be removed
723 if(cmpt->nbox==0) {
724 cmpt->volume=0;
725 return 2; } // last box was removed
726 cmpt->boxlist[bc]=cmpt->boxlist[cmpt->nbox];
727 cmpt->boxfrac[bc]=cmpt->boxfrac[cmpt->nbox];
728 boxvol=sim->boxs->boxvol;
729 vol=(bc==0)?0:cmpt->cumboxvol[bc-1];
730 for(bc2=bc;bc2<cmpt->nbox;bc2++) {
731 vol+=boxvol*cmpt->boxfrac[bc2];
732 cmpt->cumboxvol[bc2]=vol; }
733 cmpt->volume=vol;
734 return 2; }
735
736 if(bc<cmpt->nbox) { // box was listed, so just update volume
737 if(cmpt->boxfrac[bc]==volfrac2) return 0; // volume was ok, so return
738 cmpt->boxfrac[bc]=volfrac2; // volume not ok, so update it
739 boxvol=sim->boxs->boxvol;
740 vol=(bc==0)?0:cmpt->cumboxvol[bc-1];
741 for(bc2=bc;bc2<cmpt->nbox;bc2++) {
742 vol+=boxvol*cmpt->boxfrac[bc2];
743 cmpt->cumboxvol[bc2]=vol; }
744 cmpt->volume=vol;
745 return 3; }
746
747 if(cmpt->nbox==cmpt->maxbox) { // expand box list
748 max=cmpt->maxbox>0?cmpt->maxbox*2:1;
749 CHECKMEM(newboxlist=(boxptr*) calloc(max,sizeof(boxptr)));
750 CHECKMEM(newboxfrac=(double*) calloc(max,sizeof(double)));
751 CHECKMEM(newcumboxvol=(double*) calloc(max,sizeof(double)));
752 for(bc2=0;bc2<cmpt->nbox;bc2++) {
753 newboxlist[bc2]=cmpt->boxlist[bc2];
754 newboxfrac[bc2]=cmpt->boxfrac[bc2];
755 newcumboxvol[bc2]=cmpt->cumboxvol[bc2]; }
756 for(;bc2<max;bc2++) {
757 newboxlist[bc2]=NULL;
758 newboxfrac[bc2]=0;
759 newcumboxvol[bc2]=0; }
760 cmpt->maxbox=max;
761 free(cmpt->boxlist);
762 free(cmpt->boxfrac);
763 free(cmpt->cumboxvol);
764 cmpt->boxlist=newboxlist;
765 cmpt->boxfrac=newboxfrac;
766 cmpt->cumboxvol=newcumboxvol; }
767
768 bc=cmpt->nbox++; // put box into cmpt
769 cmpt->boxlist[bc]=bptr;
770 cmpt->boxfrac[bc]=volfrac2;
771 cmpt->volume+=sim->boxs->boxvol*cmpt->boxfrac[bc];
772 cmpt->cumboxvol[bc]=cmpt->volume;
773 return 1;
774
775 failure:
776 free(newboxlist);
777 free(newboxfrac);
778 free(newcumboxvol);
779 simLog(sim,10,"%s","Failed to allocate memory in compartupdatebox");
780 return -1; }
781
782
783 /* compartupdatebox, the volume fraction here is the actual volume fraction for the box inside the compartment*/
784 /* volfrac is actual volume fraction, for logic compartment the volfrac is assigned -2 */
785 #ifdef OPTION_VCELL
compartupdatebox_volumeSample(simptr sim,compartptr cmpt,boxptr bptr,double volfrac)786 int compartupdatebox_volumeSample(simptr sim,compartptr cmpt,boxptr bptr,double volfrac) {
787 /*int ptsmax=100;*/ // number of random points for volume determination
788 int bc,max,bc2;
789 double volfrac2,*newboxfrac,*newcumboxvol,boxvol,vol;
790 boxptr *newboxlist;
791
792 newboxlist=NULL;
793 newboxfrac=NULL;
794 newcumboxvol=NULL;
795
796 for(bc=0;bc<cmpt->nbox && cmpt->boxlist[bc]!=bptr;bc++); // check for box already in cmpt
797 if(bc<cmpt->nbox && volfrac==-2) return 0; // box is listed and volume ok, so return
798
799 if(volfrac>1) volfrac2=1;
800 else volfrac2=volfrac; //the actual volume faction
801
802 if(volfrac2==0) {
803 if(bc==cmpt->nbox) return 0; // box not listed and 0 volume, so return
804 cmpt->nbox--; // box was listed, so it needs to be removed
805 if(cmpt->nbox==0) {
806 cmpt->volume=0;
807 return 2; } // last box was removed
808 cmpt->boxlist[bc]=cmpt->boxlist[cmpt->nbox];
809 cmpt->boxfrac[bc]=cmpt->boxfrac[cmpt->nbox];
810 boxvol=sim->boxs->boxvol;
811 vol=(bc==0)?0:cmpt->cumboxvol[bc-1];
812 for(bc2=bc;bc2<cmpt->nbox;bc2++) {
813 vol+=boxvol*cmpt->boxfrac[bc2];
814 cmpt->cumboxvol[bc2]=vol; }
815 cmpt->volume=vol;
816 return 2; }
817
818 if(bc<cmpt->nbox) { // box was listed, so just update volume
819 if(cmpt->boxfrac[bc]==volfrac2) return 0; // volume was ok, so return
820 cmpt->boxfrac[bc]=volfrac2; // volume not ok, so update it
821 boxvol=sim->boxs->boxvol;
822 vol=(bc==0)?0:cmpt->cumboxvol[bc-1];
823 for(bc2=bc;bc2<cmpt->nbox;bc2++) {
824 vol+=boxvol*cmpt->boxfrac[bc2];
825 cmpt->cumboxvol[bc2]=vol; }
826 cmpt->volume=vol;
827 return 3; }
828
829 if(cmpt->nbox==cmpt->maxbox) { // expand box list
830 max=cmpt->maxbox>0?cmpt->maxbox*2:1;
831 CHECKMEM(newboxlist=(boxptr*)calloc(max,sizeof(boxptr)));
832 CHECKMEM(newboxfrac=(double*)calloc(max,sizeof(double)));
833 CHECKMEM(newcumboxvol=(double*)calloc(max,sizeof(double)));
834 for(bc2=0;bc2<cmpt->nbox;bc2++) {
835 newboxlist[bc2]=cmpt->boxlist[bc2];
836 newboxfrac[bc2]=cmpt->boxfrac[bc2];
837 newcumboxvol[bc2]=cmpt->cumboxvol[bc2]; }
838 for(;bc2<max;bc2++) {
839 newboxlist[bc2]=NULL;
840 newboxfrac[bc2]=0;
841 newcumboxvol[bc2]=0; }
842 cmpt->maxbox=max;
843 free(cmpt->boxlist);
844 free(cmpt->boxfrac);
845 free(cmpt->cumboxvol);
846 cmpt->boxlist=newboxlist;
847 cmpt->boxfrac=newboxfrac;
848 cmpt->cumboxvol=newcumboxvol; }
849
850 bc=cmpt->nbox++; // put box into cmpt
851 cmpt->boxlist[bc]=bptr;
852 cmpt->boxfrac[bc]=volfrac2;
853 cmpt->volume+=sim->boxs->boxvol*cmpt->boxfrac[bc];
854 cmpt->cumboxvol[bc]=cmpt->volume;
855 return 1;
856
857 failure:
858 free(newboxlist);
859 free(newboxfrac);
860 free(newcumboxvol);
861 simLog(sim,10,"%s","Failed to allocate memory in compartupdatebox_volumeSample");
862 return -1; }
863 #endif
864
865
866 /* compartreadstring */
compartreadstring(simptr sim,ParseFilePtr pfp,compartptr cmpt,const char * word,char * line2)867 compartptr compartreadstring(simptr sim,ParseFilePtr pfp,compartptr cmpt,const char *word,char *line2) {
868 char nm[STRCHAR],nm1[STRCHAR];
869 int s,er,cl,dim,itct;
870 double v1[DIMMAX];
871 enum CmptLogic sym;
872 int nvar;
873 char **varnames;
874 double *varvalues;
875
876 dim=sim->dim;
877 nvar=sim->nvar;
878 varnames=sim->varnames;
879 varvalues=sim->varvalues;
880
881 if(!strcmp(word,"name")) { // name, got[0]
882 itct=sscanf(line2,"%s",nm);
883 CHECKS(itct==1,"error reading compartment name");
884 cmpt=compartaddcompart(sim,nm);
885 CHECKS(cmpt,"failed to add compartment");
886 CHECKS(!strnword(line2,2),"unexpected text following name"); }
887
888 else if(!strcmp(word,"surface")) { // surface
889 CHECKS(cmpt,"name has to be entered before surface");
890 CHECKS(sim->srfss,"surfaces need to be entered before compartment surfaces");
891 itct=sscanf(line2,"%s",nm);
892 CHECKS(itct==1,"error reading surface name");
893 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
894 CHECKS(s>=0,"surface name '%s' not recognized",nm);
895 er=compartaddsurf(cmpt,sim->srfss->srflist[s]);
896 CHECKS(er!=1,"out of memory adding surface to compartment");
897 CHECKS(er!=2,"cannot add surface to compartment more than once");
898 CHECKS(!strnword(line2,2),"unexpected text following surface"); }
899
900 else if(!strcmp(word,"point")) { // point
901 CHECKS(cmpt,"name has to be entered before point");
902 if(dim==1) itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[0]);
903 else if(dim==2) itct=strmathsscanf(line2,"%mlg %mlg",varnames,varvalues,nvar,&v1[0],&v1[1]);
904 else itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&v1[0],&v1[1],&v1[2]);
905 CHECKS(itct==dim,"unable to read all point values");
906 er=compartaddpoint(cmpt,dim,v1);
907 CHECKS(!er,"out of memory adding point to compartment");
908 CHECKS(!strnword(line2,dim+1),"unexpected text following point"); }
909
910 else if(!strcmp(word,"compartment")) { // compartment
911 CHECKS(cmpt,"name has to be entered before compartment");
912 itct=sscanf(line2,"%s %s",nm1,nm);
913 CHECKS(itct==2,"compartment format: symbol name");
914 sym=compartstring2cl(nm1);
915 CHECKS(sym!=CLnone,"unrecognized logic symbol");
916 cl=stringfind(sim->cmptss->cnames,sim->cmptss->ncmpt,nm);
917 CHECKS(cl>=0,"cmpartment name not recognized");
918 er=compartaddcmptl(cmpt,sim->cmptss->cmptlist[cl],sym);
919 CHECKS(er!=1,"out of memory adding compartment to compartment");
920 CHECKS(er!=2,"cannot a compartment to itself");
921 CHECKS(!strnword(line2,3),"unexpected text following compartment"); }
922
923 else { // unknown word
924 CHECKS(0,"syntax error within compartment block: statement not recognized"); }
925
926 return cmpt;
927
928 failure:
929 simParseError(sim,pfp);
930 return NULL; }
931
932
933 /* loadcompart */
loadcompart(simptr sim,ParseFilePtr * pfpptr,char * line2)934 int loadcompart(simptr sim,ParseFilePtr *pfpptr,char *line2) {
935 ParseFilePtr pfp;
936 char word[STRCHAR],errstring[STRCHAR];
937 int done,pfpcode,firstline2;
938 compartptr cmpt;
939
940 pfp=*pfpptr;
941 done=0;
942 cmpt=NULL;
943 firstline2=line2?1:0;
944
945 while(!done) {
946 if(pfp->lctr==0)
947 simLog(sim,2," Reading file: '%s'\n",pfp->fname);
948 if(firstline2) {
949 strcpy(word,"name");
950 pfpcode=1;
951 firstline2=0; }
952 else
953 pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
954 *pfpptr=pfp;
955 CHECKS(pfpcode!=3,"%s",errstring);
956
957 if(pfpcode==0); // already taken care of
958 else if(pfpcode==2) { // end reading
959 done=1; }
960 else if(!strcmp(word,"end_compartment")) { // end_compartment
961 CHECKS(!line2,"unexpected text following end_compartment");
962 return 0; }
963 else if(!line2) { // just word
964 CHECKS(0,"unknown word or missing parameter"); }
965 else {
966 cmpt=compartreadstring(sim,pfp,cmpt,word,line2);
967 CHECK(cmpt); }} // failed but error has already been reported
968
969 CHECKS(0,"end of file encountered before end_compartment statement"); // end of file
970
971 failure: // failure
972 if(ErrorType!=1) simParseError(sim,pfp);
973 *pfpptr=pfp=NULL;
974 return 1; }
975
976
977 /* compartsupdateparams */
compartsupdateparams(simptr sim)978 int compartsupdateparams(simptr sim) {
979 #ifdef OPTION_VCELL
980 return compartsupdateparams_volumeSample(sim);
981 #else
982 return compartsupdateparams_original(sim);
983 #endif
984 }
985
986
987 /* compartsupdateparams */
compartsupdateparams_original(simptr sim)988 int compartsupdateparams_original(simptr sim) {
989 boxssptr boxs;
990 boxptr bptr;
991 compartssptr cmptss;
992 compartptr cmpt;
993 int b,c,s,p,inbox,er,cl;
994 double pos[3];
995 surfaceptr srf;
996 enum CmptLogic clsym;
997
998 cmptss=sim->cmptss;
999 boxs=sim->boxs;
1000 if(!boxs || !boxs->nbox) return 2;
1001
1002 for(c=0;c<cmptss->ncmpt;c++) {
1003 cmpt=cmptss->cmptlist[c];
1004 cmpt->nbox=0;
1005
1006 for(b=0;b<boxs->nbox;b++) { // find boxes that are in the compartment
1007 bptr=boxs->blist[b];
1008 inbox=0;
1009 for(p=0;p<bptr->npanel && !inbox;p++) {
1010 srf=bptr->panel[p]->srf;
1011 for(s=0;s<cmpt->nsrf && !inbox;s++)
1012 if(cmpt->surflist[s]==srf) inbox=1; } // a compartment surface is in the box
1013 if(!inbox && cmpt->ncmptl==0) {
1014 boxrandpos(sim,pos,bptr);
1015 if(posincompart(sim,pos,cmpt,0)) inbox=2; } // compartment contains whole box
1016 if(inbox) {
1017 er=compartupdatebox(sim,cmpt,bptr,inbox==2?1:-1);
1018 if(er==-1) return 1; }}
1019
1020 for(cl=0;cl<cmpt->ncmptl;cl++) { // still finding boxes that are in compartment
1021 clsym=cmpt->clsym[cl];
1022 if(clsym==CLequal || clsym==CLor || clsym==CLxor)
1023 for(b=0;b<cmpt->cmptl[cl]->nbox;b++) {
1024 bptr=cmpt->cmptl[cl]->boxlist[b];
1025 er=compartupdatebox(sim,cmpt,bptr,-2);
1026 if(er==-1) return 1; }
1027 else if(clsym==CLequalnot || clsym==CLornot)
1028 for(b=0;b<boxs->nbox;b++) {
1029 bptr=boxs->blist[b];
1030 er=compartupdatebox(sim,cmpt,bptr,-2); }}}
1031
1032 return 0; }
1033
1034
1035 #ifdef OPTION_VCELL
1036 /* compartsupdateparams_volumeSample */
compartsupdateparams_volumeSample(simptr sim)1037 int compartsupdateparams_volumeSample(simptr sim) { //VCELL
1038 //indecies
1039 int b,c,d,i,j,k,cl;
1040 //varibles used to check possible boxes and its volFrac in a specific compartment
1041 double boxLow[3], boxHigh[3], sampleLow[3], sampleHigh[3];
1042
1043 int dim = sim->dim;
1044 double* min = sim->boxs->min; // get origin in the space
1045 double* size = sim->boxs->size; //get box x,y,z size
1046 VolumeSamples* volumeSample = sim->volumeSamplesPtr;
1047 compartssptr cmptss = sim->cmptss;
1048 boxssptr boxs = sim->boxs;
1049 if(!boxs || !boxs->nbox) return 2;
1050
1051 double sampleDz = 1;
1052 if(dim > 2){
1053 sampleDz = volumeSample->size[2]/volumeSample->num[2];
1054 }
1055 double sampleDy = 1;
1056 if(dim > 1) {
1057 sampleDy = volumeSample->size[1]/volumeSample->num[1];
1058 }
1059 double sampleDx = volumeSample->size[0]/volumeSample->num[0];
1060
1061 int inbox,p,s;
1062
1063 surfaceptr srf;
1064 double pos[DIMMAX];
1065 for(c=0;c<cmptss->ncmpt;c++) {
1066 compartptr cmpt=cmptss->cmptlist[c];
1067 cmpt->nbox=0;
1068 unsigned char cmptID = getCompartmentID(cmpt->cname, sim->volumeSamplesPtr);
1069
1070 for(b=0;b<boxs->nbox;b++) { // find boxes that are in the compartment
1071 boxptr bptr=boxs->blist[b];
1072
1073 inbox=0;
1074 for(p=0;p<bptr->npanel && !inbox;p++) {
1075 srf=bptr->panel[p]->srf;
1076 for(s=0;s<cmpt->nsrf && !inbox;s++){
1077 if(cmpt->surflist[s]==srf){
1078 inbox=1;
1079 for(i=0;i<100;i++) {
1080 boxrandpos(sim,pos,bptr);
1081 }
1082 }
1083 }
1084 }
1085 if(!inbox && cmpt->ncmptl==0) {
1086 boxrandpos(sim,pos,bptr);
1087 if(posincompart(sim,pos,cmpt,0)) inbox=2; }
1088
1089 //finding box low point and high point's indexes in volume sample values.
1090 for(d=0;d<dim;d++){
1091 boxLow[d] = min[d]+size[d]*bptr->indx[d];
1092 boxHigh[d] = min[d]+size[d]*(bptr->indx[d]+1);
1093 }
1094
1095 //initialize varibles used for each box
1096 double insideCmptVol = 0;
1097
1098 for(k = 0; k < volumeSample->num[2]; k++)
1099 {
1100 sampleLow[2] = volumeSample->origin[2]+ k*sampleDz;
1101 sampleHigh[2] = volumeSample->origin[2]+ (k+1)*sampleDz;
1102 double intersectZ = std::min<double>(boxHigh[2],sampleHigh[2]) - std::max<double>(boxLow[2], sampleLow[2]);
1103 if(intersectZ<=0)
1104 {
1105 continue;
1106 }
1107 for(j = 0; j < volumeSample->num[1]; j++)
1108 {
1109 sampleLow[1] = volumeSample->origin[1]+ j*sampleDy;
1110 sampleHigh[1] = volumeSample->origin[1]+ (j+1)*sampleDy;
1111 double intersectY = std::min<double>(boxHigh[1],sampleHigh[1]) - std::max<double>(boxLow[1], sampleLow[1]);
1112 if(intersectY<=0)
1113 {
1114 continue;
1115 }
1116 for(i = 0; i < volumeSample->num[0]; i++)
1117 {
1118 sampleLow[0] = volumeSample->origin[0]+ i*sampleDx;
1119 sampleHigh[0] = volumeSample->origin[0]+ (i+1)*sampleDx;
1120 double intersectX = std::min<double>(boxHigh[0],sampleHigh[0]) - std::max<double>(boxLow[0], sampleLow[0]);
1121 if(intersectX<=0)
1122 {
1123 continue;
1124 }
1125 //find intersect volume
1126 int sampleIdx = i + j*volumeSample->num[0] + k*volumeSample->num[0]*volumeSample->num[1];
1127 if(volumeSample->volsamples[sampleIdx] == cmptID)
1128 {
1129 insideCmptVol += intersectX*intersectY*intersectZ;
1130 }
1131 }
1132 }
1133 }
1134 double boxVolfrac = insideCmptVol/((boxs->size[0])*(boxs->size[1])*(boxs->size[2]));
1135 if(boxVolfrac <= 1e-8) boxVolfrac = 0; // volume fraction is too small, consider it as 0
1136 if((boxVolfrac + 1e-8) >= 1) boxVolfrac = 1; //if volume fraction is almost 1, consider it as 1
1137 if(boxVolfrac > 0)
1138 {
1139 int errorCode = compartupdatebox_volumeSample(sim,cmpt,bptr,boxVolfrac);
1140 if(errorCode==-1) return 1;
1141 }
1142 }
1143 enum CmptLogic clsym;
1144 for(cl=0;cl<cmpt->ncmptl;cl++) { // still finding boxes that are in compartment
1145 clsym=cmpt->clsym[cl];
1146 if(clsym==CLequal || clsym==CLor || clsym==CLxor)
1147 for(b=0;b<cmpt->cmptl[cl]->nbox;b++) {
1148 boxptr bptr=cmpt->cmptl[cl]->boxlist[b];
1149 int er=compartupdatebox_volumeSample(sim,cmpt,bptr,-2);
1150 if(er==-1) return 1; }
1151 else if(clsym==CLequalnot || CLornot)
1152 for(b=0;b<boxs->nbox;b++) {
1153 boxptr bptr=boxs->blist[b];
1154 int er=compartupdatebox_volumeSample(sim,cmpt,bptr,-2);
1155 if(er==-1) return 1; }
1156 }
1157 }
1158 //for(c=0;c<cmptss->ncmpt;c++) {
1159 // compartptr cmpt=cmptss->cmptlist[c];
1160 // using namespace std;
1161 // //write to a file the box volume for each compartment
1162 // stringstream ss;
1163 // ss << "d:\\volumeCmpt" << cmpt->cname << ".txt";
1164 // ofstream filestr;
1165 // filestr.open (ss.str());
1166 // boxptr* boxlist = cmpt->boxlist;
1167 // double* boxVol = cmpt->boxfrac;
1168 // for(int i=0; i<cmpt->nbox; i++)
1169 // {
1170 // for(b=0;b<sim->boxs->nbox;b++) {
1171 // boxptr bptr=boxlist[b];
1172 // if (boxlist[i] == sim->boxs->blist[b]) {
1173 // filestr << "box index:" << b << " " << "volFrac:" << boxVol[i] << ",nneigh=" << bptr->nneigh <<", midneigh="<<bptr->midneigh<<", nwall=" <<bptr->nwall <<",maxpanel="<<bptr->maxpanel<<",npanel="<< bptr->npanel << endl;
1174 // break;
1175 // }
1176 // }
1177 // }
1178 // filestr.close();
1179 //}
1180
1181 return 0;
1182 }
1183 #endif
1184
1185
1186 /* compartsupdatelists */
compartsupdatelists(simptr sim)1187 int compartsupdatelists(simptr sim) {
1188 (void)sim;
1189 return 0; }
1190
1191
1192 /* compartsupdate */
compartsupdate(simptr sim)1193 int compartsupdate(simptr sim) {
1194 int er;
1195 compartssptr cmptss;
1196
1197 cmptss=sim->cmptss;
1198 if(cmptss) {
1199 if(cmptss->condition<=SClists) {
1200 er=compartsupdatelists(sim);
1201 if(er) return er;
1202 compartsetcondition(cmptss,SCparams,1); }
1203 if(cmptss->condition==SCparams) {
1204 er=compartsupdateparams(sim);
1205 if(er) return er;
1206 compartsetcondition(cmptss,SCok,1); }}
1207 return 0; }
1208
1209
1210 /******************************************************************************/
1211 /************************* core simulation functions **************************/
1212 /******************************************************************************/
1213
1214
comparttranslate(simptr sim,compartptr cmpt,int code,double * translate)1215 void comparttranslate(simptr sim,compartptr cmpt,int code,double *translate) {
1216 int s,dim,ll,m,d,pt,p,lxp;
1217 surfaceptr srf;
1218 moleculeptr mptr;
1219 molssptr mols;
1220 double newpos[DIMMAX],crsspt[DIMMAX],cross,epsilon;
1221 enum PanelShape ps;
1222 panelptr pnl;
1223 enum PanelFace face1,face2;
1224 enum MolecState ms;
1225
1226 dim=sim->dim;
1227 mols=sim->mols;
1228 epsilon=sim->srfss->epsilon;
1229
1230 if(code&1) { // translate surfaces and interior-defining points
1231 for(s=0;s<cmpt->nsrf;s++) {
1232 srf=cmpt->surflist[s];
1233 surfupdateoldpos(srf,dim);
1234 surftranslatesurf(srf,dim,translate); }
1235 for(pt=0;pt<cmpt->npts;pt++)
1236 for(d=0;d<dim;d++)
1237 cmpt->points[pt][d]+=translate[d]; }
1238
1239 if(code&2) { // translate surface-bound molecules
1240 for(s=0;s<cmpt->nsrf;s++) {
1241 srf=cmpt->surflist[s];
1242 for(ll=0;ll<srf->nmollist;ll++) {
1243 for(m=0;m<srf->nmol[ll];m++) {
1244 mptr=srf->mol[ll][m];
1245 if(mptr->ident) {
1246 for(d=0;d<dim;d++)
1247 mptr->pos[d]+=translate[d];
1248 ms=mptr->mstate;
1249 if(ms==MSfront) face1=PFfront;
1250 else if(ms==MSback) face1=PFback;
1251 else face1=PFnone;
1252 fixpt2panel(mptr->pos,mptr->pnl,dim,face1,epsilon); }}}}}
1253
1254 if(code&4 || code&8) { // translate molecules
1255 for(ll=0;ll<mols->nlist;ll++)
1256 for(m=0;m<mols->nl[ll];m++) {
1257 mptr=mols->live[ll][m];
1258 if(mptr->ident && mptr->mstate==MSsoln) {
1259 if(posincompart(sim,mptr->pos,cmpt,1)) { // translate molecules within compartment
1260 if(code&4)
1261 for(d=0;d<dim;d++)
1262 mptr->pos[d]+=translate[d]; }
1263 else if(code&8) { // translate external molecules
1264 for(d=0;d<dim;d++)
1265 newpos[d]=mptr->pos[d]-translate[d]; // in frame of reference of moving surface
1266 for(s=0;s<cmpt->nsrf;s++) {
1267 srf=cmpt->surflist[s];
1268 if(srf->action[mptr->ident][MSsoln][PFfront]!=SAtrans || srf->action[mptr->ident][MSsoln][PFback]!=SAtrans) {
1269 for(ps=(enum PanelShape)0;ps<PSMAX;ps=(enum PanelShape)(ps+1))
1270 for(p=0;p<srf->npanel[ps];p++) {
1271 pnl=srf->panels[ps][p];
1272 lxp=lineXpanel(mptr->pos,newpos,pnl,dim,crsspt,&face1,&face2,&cross,NULL,NULL,1);
1273 if(lxp && srf->action[mptr->ident][MSsoln][face1]!=SAtrans) {
1274 for(d=0;d<dim;d++) {
1275 mptr->via[d]=mptr->pos[d];
1276 mptr->pos[d]+=translate[d]; }
1277 checksurfaces1mol(sim,mptr,1.0-cross);
1278 }}}}}}}}
1279 sim->mols->touch++;
1280 return; }
1281
1282
1283