From 579fbdf994b239ddcf68b1364fa3c84b2837d5c8 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Sat, 19 Oct 2019 16:05:49 +0000
Subject: [PATCH] Restore seems to be working, however...

---
 src/poly.c |  222 ++++++++++++++++++++++++++++---------------------------
 1 files changed, 113 insertions(+), 109 deletions(-)

diff --git a/src/poly.c b/src/poly.c
index cfaa354..cf416b5 100644
--- a/src/poly.c
+++ b/src/poly.c
@@ -7,6 +7,7 @@
 #include<math.h>
 #include"energy.h"
 #include"cell.h"
+#include"frame.h"
 ts_bool poly_assign_filament_xi(ts_vesicle *vesicle, ts_tape *tape){
 	ts_uint i;
 
@@ -47,49 +48,20 @@
 	poly->blist->bond[i]->bond_length=sqrt(vtx_distance_sq(poly->blist->bond[i]->vtx1,poly->blist->bond[i]->vtx2));
 	bond_energy(poly->blist->bond[i],poly);
 	}
-
+	vertex_list_assign_id(poly->vlist,TS_ID_FILAMENT);
 	return poly;
 }
 
 
-ts_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist, ts_vesicle *vesicle){
-	ts_poly_list *poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
-	poly_list->poly	= (ts_poly **)calloc(n_poly,sizeof(ts_poly *));
-	ts_uint i=0,j=0,k; //idx;
-	ts_uint gvtxi;
-	ts_double xnorm,ynorm,znorm,normlength;
-	ts_double dphi,dh;
 
-	// Grafting polymers:
-	if (vlist!=NULL){
-		if (n_poly > vlist->n){fatal("Number of polymers larger than numbero f vertices on a vesicle.",310);}
-	
-		while(i<n_poly){
-			gvtxi = rand() % vlist->n;
-			if (vlist->vtx[gvtxi]->grafted_poly == NULL){
-			poly_list->poly[i] = init_poly(n_mono, vlist->vtx[gvtxi]);
-			i++;
-			}
-		}
-	}
-	else
-	{
-		for(i=0;i<n_poly;i++){
-			poly_list->poly[i] = init_poly(n_mono, NULL);
-		}
-	}
-
-	poly_list->n = n_poly;
-
-	if (vlist!=NULL){
+ts_bool poly_initial_distribution(ts_poly_list *poly_list, ts_int i, ts_vesicle *vesicle){
 	/* Make straight grafted poylmers normal to membrane (polymer brush). Dist. between poly vertices put to 1*/
+		ts_double xnorm,ynorm,znorm,normlength;
 		ts_int intpoly=vesicle->tape->internal_poly;
 		ts_int cellidx;
-		ts_double posX,posY,posZ,prevPosX,prevPosY,prevPosZ, phi,theta;
+		ts_double posX,posY,posZ,prevPosX,prevPosY,prevPosZ, phi,costheta,sintheta;
 		ts_bool retval;
-		ts_int l;
-		for (i=0;i<poly_list->n;i++){
-	
+		ts_int j,k,l,m;
 			xnorm=0.0;
 			ynorm=0.0;
 			znorm=0.0;
@@ -105,104 +77,123 @@
 			xnorm=xnorm/normlength;
 			ynorm=ynorm/normlength;
 			znorm=znorm/normlength;
-
 			//prepare starting position for building the polymeres
 			prevPosX=poly_list->poly[i]->grafted_vtx->x;
 			prevPosY=poly_list->poly[i]->grafted_vtx->y;
 			prevPosZ=poly_list->poly[i]->grafted_vtx->z;
 			for (j=0;j<poly_list->poly[i]->vlist->n;j++){
+				//if(j==0){
+					posX=prevPosX+xnorm*(vesicle->clist->dmin_interspecies);
+					posY=prevPosY+ynorm*(vesicle->clist->dmin_interspecies);
+					posZ=prevPosZ+znorm*(vesicle->clist->dmin_interspecies);
+				//}else{
+				//	posX=prevPosX+xnorm;
+				//	posY=prevPosY+ynorm;
+				//	posZ=prevPosZ+znorm;
+				//}
 				//trying to go towards normal
-				posX=prevPosX+(ts_double)xnorm;
-				posY=prevPosY+(ts_double)ynorm;
-				posZ=prevPosZ+(ts_double)znorm;
 				k=0;
-				l=0;
 				while(1){
 					poly_list->poly[i]->vlist->vtx[j]->x = posX;
 					poly_list->poly[i]->vlist->vtx[j]->y = posY;
 					poly_list->poly[i]->vlist->vtx[j]->z = posZ;
 					cellidx=vertex_self_avoidance(vesicle, poly_list->poly[i]->vlist->vtx[j]);
-					//distance from neighbors
-					//retval=TS_SUCCESS;
-					/*for(k=0;k<poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
-						dist=vtx_distance_sq(poly_list->poly[i]->vlist->vtx[j],poly_list->poly[i]->vlist->vtx[j]->neigh[k]);
-						if(dist<1.0 || dist>vesicle->dmax){
-							retval=TS_FAIL;
-							printf("dist_fail! %e\n", dist);
-						}
-					}*/
-					//if(retval!=TS_FAIL){
-						//check occupation number
-						retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,poly_list->poly[i]->vlist->vtx[j]);
-					//}
+					retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,poly_list->poly[i]->vlist->vtx[j]);
 					if(retval==TS_SUCCESS){
 						retval=cell_add_vertex(vesicle->clist->cell[cellidx],poly_list->poly[i]->vlist->vtx[j]);
 						break;
 					}
 					else{
-				//		printf("%d %d Cannot put the vesicle here. Finding another position\n",i,j);
-						theta=drand48()*M_PI-M_PI/2;
-						phi=drand48()*2*M_PI;
-						posX=prevPosX+sin(phi)*cos(theta);
-						posY=prevPosY+sin(phi)*sin(theta);
-						posZ=prevPosZ+cos(theta);
-						//randomly change the normal.
-					}
-					k++;
-					if(k>1000){
-						//lets choose another grafting vertex;
-						while(1){
-							gvtxi = rand() % vesicle->vlist->n;
-							if (vesicle->vlist->vtx[gvtxi]->grafted_poly == NULL){
-								ts_fprintf(stdout,"Found new potential grafting vertex %d for poly %d\n",gvtxi,i);
-								poly_list->poly[i]->grafted_vtx->grafted_poly=NULL;
-								poly_list->poly[i]->grafted_vtx = vesicle->vlist->vtx[gvtxi];
-								vesicle->vlist->vtx[gvtxi]->grafted_poly = poly_list->poly[i];
-								l++;
-								k=0;
+					//	printf("%d %d Cannot put the vertex here. Finding another position\n",i,j);
+						//randomly change the direction.
+						m=0;
+						//we must move first vertex into the vesicle if the normal is in or out of the vesicle if the normal is out
+						do{
+							costheta=2.0*drand48()-1.0;
+							sintheta=sqrt(1-pow(costheta,2));
+							phi=drand48()*2.0*M_PI;
+							if(j==0){
+						//for special cases, when we are on the edge of bipyramid  the distance od dmin_interspecies is not enough
+								posX=prevPosX+vesicle->dmax*sintheta*cos(phi);
+								posY=prevPosY+vesicle->dmax*sintheta*sin(phi);
+								posZ=prevPosZ+vesicle->dmax*costheta;
 
-
-
-
-			xnorm=0.0;
-			ynorm=0.0;
-			znorm=0.0;
-			int o;
-			for (o=0;o<poly_list->poly[i]->grafted_vtx->tristar_no;o++){
-				xnorm-=poly_list->poly[i]->grafted_vtx->tristar[o]->xnorm;
-				ynorm-=poly_list->poly[i]->grafted_vtx->tristar[o]->ynorm;
-				znorm-=poly_list->poly[i]->grafted_vtx->tristar[o]->znorm;	
-			}
-			normlength=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm);
-			if(intpoly && i%2){
-				normlength=-normlength;
-			}
-			xnorm=xnorm/normlength;
-			ynorm=ynorm/normlength;
-			znorm=znorm/normlength;
-
-			//prepare starting position for building the polymeres
-			posX=poly_list->poly[i]->grafted_vtx->x+xnorm;
-			posY=poly_list->poly[i]->grafted_vtx->y+ynorm;
-			posZ=poly_list->poly[i]->grafted_vtx->z+znorm;
-			
-
-
-
-
-
+							} else {
+								posX=prevPosX+vesicle->clist->dmin_interspecies*sintheta*cos(phi);
+								posY=prevPosY+vesicle->clist->dmin_interspecies*sintheta*sin(phi);
+								posZ=prevPosZ+vesicle->clist->dmin_interspecies*costheta;
+							}
+							m++;
+							if(m>1000) {
+								k=9999; //break also ot of the outer loop
 								break;
 							}
 						}
-						if(l>1000)
-							fatal("Cannot make internal polymeres. No space inside the vesicle?",1001);
+						while((xnorm*(poly_list->poly[i]->grafted_vtx->x-posX)+ynorm*(poly_list->poly[i]->grafted_vtx->y-posY)+znorm*(poly_list->poly[i]->grafted_vtx->z-posZ))>0.0 && j==0);
+					}
+					k++;
+					if(k>1000){
+						//undo changes to the cell
+						for(l=0;l<j;l++){
+							cellidx=vertex_self_avoidance(vesicle, poly_list->poly[i]->vlist->vtx[l]);
+							cell_remove_vertex(vesicle->clist->cell[cellidx],poly_list->poly[i]->vlist->vtx[l]);
+						}
+						return TS_FAIL;
 					}
 				}
 				prevPosX=posX;
 				prevPosY=posY;
 				prevPosZ=posZ;
 			}
+	printf("did it\n");
+	return TS_SUCCESS;
+
+}
+
+
+ts_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist, ts_vesicle *vesicle){
+	ts_poly_list *poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
+	poly_list->poly	= (ts_poly **)calloc(n_poly,sizeof(ts_poly *));
+	ts_uint i=0,j=0; //idx;
+	ts_uint gvtxi;
+	ts_bool retval;
+	ts_double dphi,dh;
+
+	cell_occupation(vesicle); //needed for evading the membrane
+	// Grafting polymers:
+	int tries=0;
+	if (vlist!=NULL){
+		if (n_poly > vlist->n){fatal("Number of polymers larger than numbero f vertices on a vesicle.",310);}
+		while(i<n_poly){
+			gvtxi = rand() % vlist->n;
+			if (vlist->vtx[gvtxi]->grafted_poly == NULL){
+				poly_list->poly[i] = init_poly(n_mono, vlist->vtx[gvtxi]);
+				retval=poly_initial_distribution(poly_list, i, vesicle);
+				if(retval==TS_FAIL){
+					ts_fprintf(stdout,"Found new potential grafting vertex %d for poly %d\n",gvtxi,i);
+					poly_free(poly_list->poly[i]);
+					tries++;
+				}
+				else {
+					tries=0;
+					i++;
+				}
+				if(tries>5000){
+					fatal("Cannot find space for inner polymeres",1001);
+				}
+			}
 		}
+	}
+	else
+	{
+		for(i=0;i<n_poly;i++){
+			poly_list->poly[i] = init_poly(n_mono, NULL);
+		}
+	}
+
+	poly_list->n = n_poly;
+
+	if (vlist!=NULL){
 	}
 	else
 	{
@@ -244,16 +235,16 @@
 		}
 	}
 */
-
+	
 	return poly_list;
 }
 
 
 ts_bool poly_free(ts_poly *poly){
 
-	if (poly->grafted_vtx!=NULL){
-		poly->grafted_vtx->grafted_poly=NULL;
-	}
+//	if (poly->grafted_vtx!=NULL){
+//		poly->grafted_vtx->grafted_poly=NULL;
+//	}
 	vtx_list_free(poly->vlist);
 	bond_list_free(poly->blist);
 	free(poly);
@@ -263,8 +254,9 @@
 
 ts_bool poly_list_free(ts_poly_list *poly_list){
 	ts_uint i;
-
+	//fprintf(stderr,"no. of polys=%d\n", poly_list->n);
 	for(i=0;i<poly_list->n;i++){
+	//	fprintf(stderr,"%d poly address in mem=%ld\n",i+1,(long)&(poly_list->poly[i]));
 		poly_free(poly_list->poly[i]);
 	}
 	free(poly_list->poly);
@@ -311,3 +303,15 @@
 	return TS_SUCCESS;
 }
 
+
+
+ts_poly_list *init_empty_poly_list(ts_uint n_poly, ts_uint n_mono){
+	ts_poly_list *poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
+	poly_list->poly=(ts_poly **)calloc(n_poly,sizeof(ts_poly *));
+	ts_uint i;
+	for (i=0; i<n_poly;i++){
+		poly_list->poly[i]=init_poly(n_mono, NULL);
+	}
+	poly_list->n = n_poly;
+	return poly_list;
+}

--
Gitblit v1.9.3