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 #include <stdio.h>
9 #include <math.h>
10 #include <string.h>
11 #include "string2.h"
12 #include "Zn.h"
13
14 #include "smoldyn.h"
15 #include "smoldynfuncs.h"
16
17 /******************************************************************************/
18 /************************************* Ports **********************************/
19 /******************************************************************************/
20
21
22 /******************************************************************************/
23 /****************************** Local declarations ****************************/
24 /******************************************************************************/
25
26 // memory management
27 portptr portalloc(void);
28 void portfree(portptr port);
29 portssptr portssalloc(portssptr portss,int maxport);
30
31 // data structure output
32
33 // structure set up
34 int portsupdateparams(simptr sim);
35 int portsupdatelists(simptr sim);
36
37 // core simulation functions
38
39
40 /******************************************************************************/
41 /******************************* memory management ****************************/
42 /******************************************************************************/
43
44 /* portalloc */
portalloc(void)45 portptr portalloc(void) {
46 portptr port;
47
48 port=(portptr)malloc(sizeof(struct portstruct));
49 CHECKMEM(port);
50 port->portname=NULL;
51 port->srf=NULL;
52 port->face=PFnone;
53 port->llport=-1;
54 return port;
55 failure:
56 simLog(NULL,10,"Unable to allocate memory in portalloc");
57 return NULL; }
58
59
60 /* portfree */
portfree(portptr port)61 void portfree(portptr port) {
62 if(!port) return;
63 free(port);
64 return; }
65
66
67 /* portssalloc */
portssalloc(portssptr portss,int maxport)68 portssptr portssalloc(portssptr portss,int maxport) {
69 int p;
70 char **newnames;
71 portptr *newportlist;
72
73 if(maxport<1) return NULL;
74
75 newnames=NULL;
76 newportlist=NULL;
77
78 if(!portss) { // new allocation
79 portss=(portssptr) malloc(sizeof(struct portsuperstruct));
80 CHECKMEM(portss);
81 portss->condition=SCinit;
82 portss->sim=NULL;
83 portss->maxport=0;
84 portss->nport=0;
85 portss->portnames=NULL;
86 portss->portlist=NULL; }
87 else { // minor check
88 if(maxport<portss->maxport) return NULL; }
89
90 if(maxport>portss->maxport) { // allocate new port names and ports
91 CHECKMEM(newnames=(char**) calloc(maxport,sizeof(char*)));
92 for(p=0;p<maxport;p++) newnames[p]=NULL;
93 for(p=0;p<portss->maxport;p++) newnames[p]=portss->portnames[p];
94 for(;p<maxport;p++)
95 CHECKMEM(newnames[p]=EmptyString());
96
97 CHECKMEM(newportlist=(portptr*) calloc(maxport,sizeof(portptr))); // port list
98 for(p=0;p<maxport;p++) newportlist[p]=NULL;
99 for(p=0;p<portss->maxport;p++) newportlist[p]=portss->portlist[p];
100 for(;p<maxport;p++) {
101 CHECKMEM(newportlist[p]=portalloc());
102 newportlist[p]->portss=portss;
103 newportlist[p]->portname=newnames[p]; }}
104
105 portss->maxport=maxport;
106 free(portss->portnames);
107 portss->portnames=newnames;
108 free(portss->portlist);
109 portss->portlist=newportlist;
110
111 return portss;
112
113 failure:
114 portssfree(portss);
115 simLog(NULL,10,"Unable to allocate memory in portssalloc");
116 return NULL; }
117
118
119 /* portssfree */
portssfree(portssptr portss)120 void portssfree(portssptr portss) {
121 int prt;
122
123 if(!portss) return;
124 if(portss->maxport && portss->portlist)
125 for(prt=0;prt<portss->maxport;prt++) portfree(portss->portlist[prt]);
126 free(portss->portlist);
127 if(portss->maxport && portss->portnames)
128 for(prt=0;prt<portss->maxport;prt++) free(portss->portnames[prt]);
129 free(portss->portnames);
130 free(portss);
131 return; }
132
133
134 /******************************************************************************/
135 /***************************** data structure output **************************/
136 /******************************************************************************/
137
138 /* portoutput */
portoutput(simptr sim)139 void portoutput(simptr sim) {
140 portssptr portss;
141 portptr port;
142 int prt;
143 char string[STRCHAR];
144
145 portss=sim->portss;
146 if(!portss) return;
147 simLog(sim,2,"PORT PARAMETERS\n");
148 simLog(sim,2," Ports allocated: %i, ports defined: %i\n",portss->maxport,portss->nport);
149 for(prt=0;prt<portss->nport;prt++) {
150 port=portss->portlist[prt];
151 simLog(sim,2," Port: %s\n",portss->portnames[prt]);
152 if(port->srf) simLog(sim,2," surface: %s, %s\n",port->srf->sname,surfface2string(port->face,string));
153 else simLog(sim,2," no surface assigned\n");
154 if(port->llport>=0) simLog(sim,2," molecule list: %s\n",sim->mols->listname[port->llport]);
155 else simLog(sim,2," no molecule list assigned"); }
156 simLog(sim,2,"\n");
157 return; }
158
159
160 /* writeports */
writeports(simptr sim,FILE * fptr)161 void writeports(simptr sim,FILE *fptr) {
162 portssptr portss;
163 portptr port;
164 int prt;
165 char string[STRCHAR];
166
167 portss=sim->portss;
168 if(!portss) return;
169 fprintf(fptr,"# Port parameters\n");
170 fprintf(fptr,"max_port %i\n",portss->maxport);
171 for(prt=0;prt<portss->nport;prt++) {
172 port=portss->portlist[prt];
173 fprintf(fptr,"start_port %s\n",port->portname);
174 if(port->srf) fprintf(fptr,"surface %s\n",port->srf->sname);
175 if(port->srf) fprintf(fptr,"face %s\n",surfface2string(port->face,string));
176 fprintf(fptr,"end_port\n\n"); }
177 return; }
178
179
180 /* checkportparams */
checkportparams(simptr sim,int * warnptr)181 int checkportparams(simptr sim,int *warnptr) {
182 int error,warn,prt,er,i;
183 portssptr portss;
184 portptr port;
185 char string[STRCHAR];
186
187 error=warn=0;
188 portss=sim->portss;
189 if(!portss) {
190 if(warnptr) *warnptr=warn;
191 return 0; }
192
193 if(portss->condition!=SCok) {
194 warn++;
195 simLog(sim,7," WARNING: port structure %s\n",simsc2string(portss->condition,string)); }
196
197 for(prt=0;prt<portss->nport;prt++) {
198 port=portss->portlist[prt]; // check for porting surface
199 if(!port->srf) {warn++;simLog(sim,5," WARNING: there is no porting surface assigned to port %s\n",port->portname);}
200 if(!(port->face==PFfront || port->face==PFback))
201 {warn++;simLog(sim,5," WARNING: no surface face has been assigned to port %s\n",port->portname);}
202
203 if(port->srf && port->srf->port[port->face]!=port)
204 {error++;simLog(sim,10," ERROR: port %s is not registered by surface %s\n",port->portname,port->srf->sname);}
205
206 if(sim->mols && port->srf && port->srf->action) {
207 er=1; // make sure surface action is set to port
208 for(i=0;i<sim->mols->nspecies && er;i++)
209 if(port->srf->action[i][MSsoln][port->face]==SAport) er=0;
210 if(er) {warn++;simLog(sim,5," WARNING: port %s is nonfunctional because no molecule actions at the surface %s are set to port\n",port->portname,port->srf->sname);}
211 if(!port->llport) {error++;simLog(sim,10," BUG: port %s has no molecule buffer\n",port->portname);} }}
212
213 if(warnptr) *warnptr=warn;
214 return error; }
215
216
217 /******************************************************************************/
218 /******************************** structure set up ****************************/
219 /******************************************************************************/
220
221
222 /* portsetcondition */
portsetcondition(portssptr portss,enum StructCond cond,int upgrade)223 void portsetcondition(portssptr portss,enum StructCond cond,int upgrade) {
224 if(!portss) return;
225 if(upgrade==0 && portss->condition>cond) portss->condition=cond;
226 else if(upgrade==1 && portss->condition<cond) portss->condition=cond;
227 else if(upgrade==2) portss->condition=cond;
228 if(portss->sim && portss->condition<portss->sim->condition) {
229 cond=portss->condition;
230 simsetcondition(portss->sim,cond==SCinit?SClists:cond,0); }
231 return; }
232
233
234 /* portenableports */
portenableports(simptr sim,int maxport)235 int portenableports(simptr sim,int maxport) {
236 portssptr portss;
237
238 if(sim->portss) // check for redundant function call
239 if(maxport==-1 || sim->portss->maxport==maxport)
240 return 0;
241 portss=portssalloc(sim->portss,maxport<0?5:maxport);
242 if(!portss) return 1;
243 sim->portss=portss;
244 portss->sim=sim;
245 portsetcondition(sim->portss,SClists,0);
246 return 0; }
247
248
249 /* portaddport */
portaddport(simptr sim,const char * portname,surfaceptr srf,enum PanelFace face)250 portptr portaddport(simptr sim,const char *portname,surfaceptr srf,enum PanelFace face) {
251 int er,p;
252 portssptr portss;
253 portptr port;
254
255 if(!sim->portss) {
256 er=portenableports(sim,-1);
257 if(er) return NULL; }
258 portss=sim->portss;
259
260 p=stringfind(portss->portnames,portss->nport,portname);
261 if(p<0) {
262 if(portss->nport==portss->maxport) {
263 er=portenableports(sim,portss->nport*2+1);
264 if(er) return NULL; }
265 p=portss->nport++;
266 strncpy(portss->portnames[p],portname,STRCHAR-1);
267 portss->portnames[p][STRCHAR-1]='\0';
268 port=portss->portlist[p]; }
269 else
270 port=portss->portlist[p];
271
272 if(srf) port->srf=srf;
273 if(face!=PFnone) port->face=face;
274 if(port->srf && port->face!=PFnone)
275 port->srf->port[port->face]=port;
276 portsetcondition(portss,SClists,0);
277 return port; }
278
279
280 /* portreadstring */
portreadstring(simptr sim,ParseFilePtr pfp,portptr port,const char * word,char * line2)281 portptr portreadstring(simptr sim,ParseFilePtr pfp,portptr port,const char *word,char *line2) {
282 char nm[STRCHAR];
283 int itct,s;
284 enum PanelFace face;
285
286
287 if(!strcmp(word,"name")) { // name
288 itct=sscanf(line2,"%s",nm);
289 CHECKS(itct==1,"error reading port name");
290 port=portaddport(sim,nm,NULL,PFnone);
291 CHECKS(port,"failed to add port");
292 CHECKS(!strnword(line2,2),"unexpected text following name"); }
293
294 else if(!strcmp(word,"surface")) { // surface
295 CHECKS(port,"port name has to be entered before surface");
296 itct=sscanf(line2,"%s",nm);
297 CHECKS(itct==1,"error reading surface name");
298 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
299 CHECKS(s>=0,"surface '%s' not recognized",nm);
300 port=portaddport(sim,port->portname,sim->srfss->srflist[s],PFnone);
301 CHECKBUG(port,"SMOLDYN BUG adding surface to port");
302 CHECKS(!strnword(line2,2),"unexpected text following surface"); }
303
304 else if(!strcmp(word,"face")) { // face
305 CHECKS(port,"port name has to be entered before face");
306 itct=sscanf(line2,"%s",nm);
307 CHECKS(itct==1,"error reading face name");
308 face=surfstring2face(nm);
309 CHECKS(face==PFfront || face==PFback,"face needs to be either front or back");
310 port=portaddport(sim,port->portname,NULL,face);
311 CHECKBUG(port,"SMOLDYN BUG adding face to port");
312 CHECKS(!strnword(line2,2),"unexpected text following face"); }
313
314 else { // unknown word
315 CHECKS(0,"syntax error within port block: statement not recognized"); }
316
317 return port;
318
319 failure:
320 simParseError(sim,pfp);
321 return NULL; }
322
323
324 /* loadport */
loadport(simptr sim,ParseFilePtr * pfpptr,char * line2)325 int loadport(simptr sim,ParseFilePtr *pfpptr,char* line2) {
326 ParseFilePtr pfp;
327 char word[STRCHAR],errstring[STRCHAR];
328 int done,pfpcode,firstline2;
329 portptr port;
330
331 pfp=*pfpptr;
332 done=0;
333 port=NULL;
334 firstline2=line2?1:0;
335
336 while(!done) {
337 if(pfp->lctr==0)
338 simLog(sim,2," Reading file: '%s'\n",pfp->fname);
339 if(firstline2) {
340 strcpy(word,"name");
341 pfpcode=1;
342 firstline2=0; }
343 else
344 pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
345 *pfpptr=pfp;
346 CHECKS(pfpcode!=3,"%s",errstring);
347
348 if(pfpcode==0); // already taken care of
349 else if(pfpcode==2) { // end reading
350 done=1; }
351 else if(!strcmp(word,"end_port")) { // end_port
352 CHECKS(!line2,"unexpected text following end_port");
353 return 0; }
354 else if(!line2) { // just word
355 CHECKS(0,"unknown word or missing parameter"); }
356 else {
357 port=portreadstring(sim,pfp,port,word,line2);
358 CHECK(port); }} // failed but error has already been reported
359
360 CHECKS(0,"end of file encountered before end_port statement"); // end of file
361
362 failure: // failure
363 if(ErrorType!=1) simParseError(sim,pfp);
364 *pfpptr=pfp=NULL;
365 return 1; }
366
367
368 /* portsupdateparams */
portsupdateparams(simptr sim)369 int portsupdateparams(simptr sim) {
370 (void)sim;
371 return 0; }
372
373
374 /* portsupdatelists */
portsupdatelists(simptr sim)375 int portsupdatelists(simptr sim) {
376 portssptr portss;
377 portptr port;
378 int prt,ll;
379
380 portss=sim->portss;
381 if(sim->mols && sim->mols->condition<SCparams) return 2;
382
383 if(sim->mols) {
384 for(prt=0;prt<portss->nport;prt++) {
385 port=portss->portlist[prt]; // create port molecule list if none
386 if(port->llport<0) {
387 ll=addmollist(sim,port->portname,MLTport);
388 if(ll<0) return 1;
389 port->llport=ll; }}}
390
391 return 0; }
392
393
394 /* portsupdate */
portsupdate(simptr sim)395 int portsupdate(simptr sim) {
396 int er;
397 portssptr portss;
398
399 portss=sim->portss;
400 if(portss) {
401 if(portss->condition<=SClists) {
402 er=portsupdatelists(sim);
403 if(er) return er;
404 portsetcondition(portss,SCparams,1); }
405 if(portss->condition==SCparams) {
406 er=portsupdateparams(sim);
407 if(er) return er;
408 portsetcondition(portss,SCok,1); }}
409 return 0; }
410
411
412 /******************************************************************************/
413 /*************************** core simulation functions ************************/
414 /******************************************************************************/
415
416
417 /* portgetmols */
portgetmols(simptr sim,portptr port,int ident,enum MolecState ms,int remove)418 int portgetmols(simptr sim,portptr port,int ident,enum MolecState ms,int remove) {
419 int ll,nmol,count,m;
420 moleculeptr *mlist;
421
422 ll=port->llport;
423 mlist=sim->mols->live[ll];
424 nmol=sim->mols->nl[ll];
425 if(ident<0 && ms==MSall && !remove) return nmol;
426 count=0;
427 for(m=0;m<nmol;m++) {
428 if((ident==-1 || mlist[m]->ident==ident) && (ms==MSall || mlist[m]->mstate==ms)) {
429 count++;
430 if(remove) molkill(sim,mlist[m],ll,m); }}
431 sim->eventcount[ETexport]+=count;
432 return count; }
433
434
435 /* portgetmols2 */
portgetmols2(simptr sim,portptr port,int ident,enum MolecState ms,int remove,double ** positions)436 int portgetmols2(simptr sim,portptr port,int ident,enum MolecState ms,int remove, double **positions) {
437 int ll,nmol,count,m;
438 moleculeptr *mlist;
439
440 ll=port->llport;
441 mlist=sim->mols->live[ll];
442 nmol=sim->mols->nl[ll];
443 if(ident<0 && ms==MSall && !remove && !positions) return nmol;
444 count=0;
445 for(m=0;m<nmol;m++) {
446 if((ident==-1 || mlist[m]->ident==ident) && (ms==MSall || mlist[m]->mstate==ms)) {
447 count++;
448 if(positions) positions[count] = mlist[m]->pos;
449 if(remove) molkill(sim,mlist[m],ll,m); }}
450 sim->eventcount[ETexport]+=count;
451 return count; }
452
453
454 /* portputmols */
portputmols(simptr sim,portptr port,int nmol,int ident,int * species,double ** positions,double ** positionsx)455 int portputmols(simptr sim,portptr port,int nmol,int ident,int *species,double **positions,double **positionsx) {
456 moleculeptr mptr;
457 int dim,m,d;
458 panelptr pnl;
459
460 if(!nmol) return 0;
461 if(!positionsx) { // no positionsx is ok, but then the port has to have a surface
462 if(!port->srf) return 2;
463 if(port->face==PFnone) return 3;
464 if(port->srf->totpanel==0) return 4; }
465 dim=sim->dim;
466
467 for(m=0;m<nmol;m++) {
468 mptr=getnextmol(sim->mols);
469 if(!mptr) return 1;
470 if(species) mptr->ident=species[m];
471 else mptr->ident=ident;
472 mptr->mstate=MSsoln;
473 mptr->list=sim->mols->listlookup[mptr->ident][MSsoln];
474 sim->mols->expand[mptr->ident]|=1;
475 if(positionsx) {
476 for(d=0;d<dim;d++) {
477 mptr->pos[d]=positions[m][d];
478 mptr->posx[d]=positionsx[m][d];}}
479 else if(positions) {
480 closestsurfacept(port->srf,sim->dim,positions[m],mptr->posx,&pnl,NULL);
481 fixpt2panel(mptr->posx,pnl,dim,port->face,sim->srfss->epsilon);
482 for(d=0;d<dim;d++) mptr->pos[d]=positions[m][d]; }
483 else {
484 pnl=surfrandpos(port->srf,mptr->posx,dim);
485 fixpt2panel(mptr->posx,pnl,dim,port->face,sim->srfss->epsilon);
486 for(d=0;d<dim;d++) mptr->pos[d]=mptr->posx[d]; }
487 mptr->box=pos2box(sim,mptr->pos); }
488 sim->eventcount[ETimport]+=nmol;
489 return 0; }
490
491
492 /* porttransport */
porttransport(simptr sim1,portptr port1,simptr sim2,portptr port2)493 int porttransport(simptr sim1,portptr port1,simptr sim2,portptr port2) {
494 int i,nmol,er;
495
496 if(!portgetmols(sim1,port1,-1,MSall,0)) return 0;
497 er=0;
498 for(i=1;i<sim1->mols->nspecies && !er;i++) {
499 nmol=portgetmols(sim1,port1,i,MSall,1);
500 er=portputmols(sim2,port2,nmol,i,NULL,NULL,NULL); }
501 return er; }
502
503