| | |
| | | #include "general.h" |
| | | #include "energy.h" |
| | | #include "vertex.h" |
| | | #include "bond.h" |
| | | #include<math.h> |
| | | #include<stdio.h> |
| | | |
| | |
| | | * @returns TS_SUCCESS on successful calculation. |
| | | */ |
| | | inline ts_bool energy_vertex(ts_vertex *vtx){ |
| | | ts_uint jj; |
| | | ts_uint jjp,jjm; |
| | | ts_vertex *j,*jp, *jm; |
| | | ts_triangle *jt; |
| | | ts_double s=0.0,xh=0.0,yh=0.0,zh=0.0,txn=0.0,tyn=0.0,tzn=0.0; |
| | | ts_double x1,x2,x3,ctp,ctm,tot,xlen; |
| | | ts_double h,ht; |
| | | for(jj=1; jj<=vtx->neigh_no;jj++){ |
| | | jjp=jj+1; |
| | | if(jjp>vtx->neigh_no) jjp=1; |
| | | jjm=jj-1; |
| | | if(jjm<1) jjm=vtx->neigh_no; |
| | | j=vtx->neigh[jj-1]; |
| | | jp=vtx->neigh[jjp-1]; |
| | | jm=vtx->neigh[jjm-1]; |
| | | jt=vtx->tristar[jj-1]; |
| | | x1=vtx_distance_sq(vtx,jp); //shouldn't be zero! |
| | | x2=vtx_distance_sq(j,jp); // shouldn't be zero! |
| | | x3=(j->x-jp->x)*(vtx->x-jp->x)+ |
| | | (j->y-jp->y)*(vtx->y-jp->y)+ |
| | | (j->z-jp->z)*(vtx->z-jp->z); |
| | | |
| | | #ifdef TS_DOUBLE_DOUBLE |
| | | ctp=x3/sqrt(x1*x2-x3*x3); |
| | | #endif |
| | | #ifdef TS_DOUBLE_FLOAT |
| | | ctp=x3/sqrtf(x1*x2-x3*x3); |
| | | #endif |
| | | #ifdef TS_DOUBLE_LONGDOUBLE |
| | | ctp=x3/sqrtl(x1*x2-x3*x3); |
| | | #endif |
| | | x1=vtx_distance_sq(vtx,jm); |
| | | x2=vtx_distance_sq(j,jm); |
| | | x3=(j->x-jm->x)*(vtx->x-jm->x)+ |
| | | (j->y-jm->y)*(vtx->y-jm->y)+ |
| | | (j->z-jm->z)*(vtx->z-jm->z); |
| | | #ifdef TS_DOUBLE_DOUBLE |
| | | ctm=x3/sqrt(x1*x2-x3*x3); |
| | | #endif |
| | | #ifdef TS_DOUBLE_FLOAT |
| | | ctm=x3/sqrtf(x1*x2-x3*x3); |
| | | #endif |
| | | #ifdef TS_DOUBLE_LONGDOUBLE |
| | | ctm=x3/sqrtl(x1*x2-x3*x3); |
| | | #endif |
| | | tot=ctp+ctm; |
| | | tot=0.5*tot; |
| | | ts_uint jj, i, j, cnt=0; |
| | | ts_double edge_vector_x[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_vector_y[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_vector_z[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_normal_x[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_normal_y[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_normal_z[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_binormal_x[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_binormal_y[7]={0,0,0,0,0,0,0}; |
| | | ts_double edge_binormal_z[7]={0,0,0,0,0,0,0}; |
| | | ts_double vertex_normal_x=0.0; |
| | | ts_double vertex_normal_y=0.0; |
| | | ts_double vertex_normal_z=0.0; |
| | | ts_triangle *triedge[2]={NULL,NULL}; |
| | | |
| | | xlen=vtx_distance_sq(j,vtx); |
| | | /* |
| | | #ifdef TS_DOUBLE_DOUBLE |
| | | vtx->bond[jj-1]->bond_length=sqrt(xlen); |
| | | #endif |
| | | #ifdef TS_DOUBLE_FLOAT |
| | | vtx->bond[jj-1]->bond_length=sqrtf(xlen); |
| | | #endif |
| | | #ifdef TS_DOUBLE_LONGDOUBLE |
| | | vtx->bond[jj-1]->bond_length=sqrtl(xlen); |
| | | #endif |
| | | ts_double sumnorm; |
| | | |
| | | vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length; |
| | | */ |
| | | s+=tot*xlen; |
| | | xh+=tot*(j->x - vtx->x); |
| | | yh+=tot*(j->y - vtx->y); |
| | | zh+=tot*(j->z - vtx->z); |
| | | txn+=jt->xnorm; |
| | | tyn+=jt->ynorm; |
| | | tzn+=jt->znorm; |
| | | } |
| | | |
| | | h=xh*xh+yh*yh+zh*zh; |
| | | ht=txn*xh+tyn*yh + tzn*zh; |
| | | s=s/4.0; |
| | | #ifdef TS_DOUBLE_DOUBLE |
| | | if(ht>=0.0) { |
| | | vtx->curvature=sqrt(h); |
| | | } else { |
| | | vtx->curvature=-sqrt(h); |
| | | } |
| | | #endif |
| | | #ifdef TS_DOUBLE_FLOAT |
| | | if(ht>=0.0) { |
| | | vtx->curvature=sqrtf(h); |
| | | } else { |
| | | vtx->curvature=-sqrtf(h); |
| | | } |
| | | #endif |
| | | #ifdef TS_DOUBLE_LONGDOUBLE |
| | | if(ht>=0.0) { |
| | | vtx->curvature=sqrtl(h); |
| | | } else { |
| | | vtx->curvature=-sqrtl(h); |
| | | } |
| | | #endif |
| | | // c is forced curvature energy for each vertex. Should be set to zero for |
| | | // normal circumstances. |
| | | /* the following statement is an expression for $\frac{1}{2}\int(c_1+c_2-c_0^\prime)^2\mathrm{d}A$, where $c_0^\prime=2c_0$ (twice the spontaneous curvature) */ |
| | | vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c); |
| | | // Here edge vector is calculated |
| | | // fprintf(stderr, "Vertex has neighbours=%d\n", vtx->neigh_no); |
| | | for(jj=0;jj<vtx->neigh_no;jj++){ |
| | | edge_vector_x[jj]=vtx->neigh[jj]->x-vtx->x; |
| | | edge_vector_y[jj]=vtx->neigh[jj]->y-vtx->y; |
| | | edge_vector_z[jj]=vtx->neigh[jj]->z-vtx->z; |
| | | // We find lm and lp from k->tristar ! |
| | | cnt=0; |
| | | for(i=0;i<vtx->tristar_no;i++){ |
| | | for(j=0;j<vtx->neigh[jj]->tristar_no;j++){ |
| | | if((vtx->tristar[i] == vtx->neigh[jj]->tristar[j])){ //ce gre za skupen trikotnik |
| | | triedge[cnt]=vtx->tristar[i]; |
| | | cnt++; |
| | | } |
| | | } |
| | | } |
| | | if(cnt!=2) fatal("ts_energy_vertex: both triangles not found!", 133); |
| | | sumnorm=sqrt( pow((triedge[0]->xnorm + triedge[1]->xnorm),2) + pow((triedge[0]->ynorm + triedge[1]->ynorm), 2) + pow((triedge[0]->znorm + triedge[1]->znorm), 2)); |
| | | |
| | | return TS_SUCCESS; |
| | | edge_normal_x[jj]=(triedge[0]->xnorm+ triedge[1]->xnorm)/sumnorm; |
| | | edge_normal_y[jj]=(triedge[0]->ynorm+ triedge[1]->ynorm)/sumnorm; |
| | | edge_normal_z[jj]=(triedge[0]->znorm+ triedge[1]->znorm)/sumnorm; |
| | | |
| | | |
| | | edge_binormal_x[jj]=(edge_normal_y[jj]*edge_vector_z[jj])-(edge_normal_z[jj]*edge_vector_y[jj]); |
| | | edge_binormal_y[jj]=-(edge_normal_x[jj]*edge_vector_z[jj])+(edge_normal_z[jj]*edge_vector_x[jj]); |
| | | edge_binormal_z[jj]=(edge_normal_x[jj]*edge_vector_y[jj])-(edge_normal_y[jj]*edge_vector_x[jj]); |
| | | |
| | | printf("(%f %f %f); (%f %f %f); (%f %f %f), %d\n", edge_vector_x[jj], edge_vector_y[jj], edge_vector_z[jj], edge_normal_x[jj], edge_normal_y[jj], edge_normal_z[jj], edge_binormal_x[jj], edge_binormal_y[jj], edge_binormal_z[jj],triedge[0]->idx); |
| | | |
| | | } |
| | | for(i=0; i<vtx->tristar_no; i++){ |
| | | vertex_normal_x=vertex_normal_x + vtx->tristar[i]->xnorm*vtx->tristar[i]->area; |
| | | vertex_normal_y=vertex_normal_y + vtx->tristar[i]->ynorm*vtx->tristar[i]->area; |
| | | vertex_normal_z=vertex_normal_z + vtx->tristar[i]->znorm*vtx->tristar[i]->area; |
| | | } |
| | | printf("(%f %f %f)\n", vertex_normal_x, vertex_normal_y, vertex_normal_z); |
| | | vtx->energy=0.0; |
| | | return TS_SUCCESS; |
| | | } |
| | | |
| | | |
| | | |
| | | ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle){ |
| | | int i; |