| | |
| | | #include "bond.h" |
| | | #include<math.h> |
| | | #include<stdio.h> |
| | | #include <gsl/gsl_vector_complex.h> |
| | | #include <gsl/gsl_matrix.h> |
| | | #include <gsl/gsl_eigen.h> |
| | | |
| | | |
| | | |
| | | int cmpfunc(const void *x, const void *y) |
| | | { |
| | | double diff= fabs(*(double*)x) - fabs(*(double*)y); |
| | | if(diff<0) return 1; |
| | | else return -1; |
| | | } |
| | | |
| | | |
| | | |
| | | /** @brief Wrapper that calculates energy of every vertex in vesicle |
| | |
| | | ts_vertex *it, *k, *kp,*km; |
| | | ts_triangle *lm=NULL, *lp=NULL; |
| | | ts_double sumnorm; |
| | | ts_double temp_length; |
| | | |
| | | |
| | | ts_double Se11, Se21, Se22, Se31, Se32, Se33; |
| | | ts_double Pv11, Pv21, Pv22, Pv31, Pv32, Pv33; |
| | | ts_double We; |
| | | ts_double Av, We_Av; |
| | | |
| | | ts_double eigenval[3]; |
| | | |
| | | gsl_matrix *gsl_Sv=gsl_matrix_alloc(3,3); |
| | | gsl_vector *Sv_eigen=gsl_vector_alloc(3); |
| | | gsl_eigen_symm_workspace *workspace=gsl_eigen_symm_alloc(3); |
| | | |
| | | ts_double mprod[7], phi[7], he[7]; |
| | | ts_double Sv[3][3]={{0,0,0},{0,0,0},{0,0,0}}; |
| | | // Here edge vector is calculated |
| | | // fprintf(stderr, "Vertex has neighbours=%d\n", vtx->neigh_no); |
| | | |
| | | |
| | | |
| | | |
| | | Av=0; |
| | | 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; |
| | | Av+=vtx->tristar[i]->area/3; |
| | | } |
| | | temp_length=sqrt(pow(vertex_normal_x,2)+pow(vertex_normal_y,2)+pow(vertex_normal_z,2)); |
| | | vertex_normal_x=vertex_normal_x/temp_length; |
| | | vertex_normal_y=vertex_normal_y/temp_length; |
| | | vertex_normal_z=vertex_normal_z/temp_length; |
| | | |
| | | Pv11=1-vertex_normal_x*vertex_normal_x; |
| | | Pv22=1-vertex_normal_y*vertex_normal_y; |
| | | Pv33=1-vertex_normal_z*vertex_normal_z; |
| | | Pv21=vertex_normal_x*vertex_normal_y; |
| | | Pv31=vertex_normal_x*vertex_normal_z; |
| | | Pv32=vertex_normal_y*vertex_normal_z; |
| | | |
| | | |
| | | |
| | | |
| | | 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; |
| | | |
| | | //Here we calculate normalized edge vector |
| | | |
| | | temp_length=sqrt(edge_vector_x[jj]*edge_vector_x[jj]+edge_vector_y[jj]*edge_vector_y[jj]+edge_vector_z[jj]*edge_vector_z[jj]); |
| | | edge_vector_x[jj]=edge_vector_x[jj]/temp_length; |
| | | edge_vector_y[jj]=edge_vector_y[jj]/temp_length; |
| | | edge_vector_z[jj]=edge_vector_z[jj]/temp_length; |
| | | |
| | | //end normalization |
| | | // printf("(%f %f %f)\n", vertex_normal_x, vertex_normal_y, vertex_normal_z); |
| | | |
| | | |
| | | it=vtx; |
| | |
| | | kp=it->neigh[neip]; |
| | | |
| | | if(km==NULL || kp==NULL){ |
| | | fatal("In bondflip, cannot determine km and kp!",999); |
| | | fatal("energy_vertex: cannot determine km and kp!",233); |
| | | } |
| | | |
| | | for(i=0;i<it->tristar_no;i++){ |
| | |
| | | } |
| | | } |
| | | } |
| | | if(lm==NULL || lp==NULL) fatal("ts_flip_bond: Cannot find triangles lm and lp!",999); |
| | | if(lm==NULL || lp==NULL) fatal("energy_vertex: Cannot find triangles lm and lp!",233); |
| | | |
| | | |
| | | /* |
| | | // 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); |
| | | */ |
| | | //Triangle normals are NORMALIZED! |
| | | |
| | | sumnorm=sqrt( pow((lm->xnorm + lp->xnorm),2) + pow((lm->ynorm + lp->ynorm), 2) + pow((lm->znorm + lp->znorm), 2)); |
| | | |
| | |
| | | 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)\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]); |
| | | |
| | | } |
| | | 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; |
| | | mprod[jj]=it->x*(k->y*edge_vector_z[jj]-edge_vector_y[jj]*k->z)-it->y*(k->x*edge_vector_z[jj]-k->z*edge_vector_x[jj])+it->z*(k->x*edge_vector_y[jj]-k->y*edge_vector_x[jj]); |
| | | phi[jj]=copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm-1e-15),mprod[jj])+M_PI; |
| | | // printf("ACOS arg=%e\n", lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm); |
| | | //he was multiplied with 2 before... |
| | | he[jj]=sqrt( pow((edge_vector_x[jj]),2) + pow((edge_vector_y[jj]), 2) + pow((edge_vector_z[jj]), 2))*cos(phi[jj]/2.0); |
| | | // printf("phi[%d]=%f\n", jj,phi[jj]); |
| | | |
| | | Se11=edge_binormal_x[jj]*edge_binormal_x[jj]*he[jj]; |
| | | Se21=edge_binormal_x[jj]*edge_binormal_y[jj]*he[jj]; |
| | | Se22=edge_binormal_y[jj]*edge_binormal_y[jj]*he[jj]; |
| | | Se31=edge_binormal_x[jj]*edge_binormal_z[jj]*he[jj]; |
| | | Se32=edge_binormal_y[jj]*edge_binormal_z[jj]*he[jj]; |
| | | Se33=edge_binormal_z[jj]*edge_binormal_z[jj]*he[jj]; |
| | | |
| | | We=vertex_normal_x*edge_normal_x[jj]+vertex_normal_y*edge_normal_y[jj]+vertex_normal_z*edge_normal_z[jj]; |
| | | We_Av=We/Av; |
| | | |
| | | Sv[0][0]+=We_Av* ( Pv11*(Pv11*Se11+Pv21*Se21+Pv31*Se31)+Pv21*(Pv11*Se21+Pv21*Se22+Pv31*Se32)+Pv31*(Pv11*Se31+Pv21*Se32+Pv31*Se33) ); |
| | | Sv[0][1]+=We_Av* (Pv21*(Pv11*Se11+Pv21*Se21+Pv31*Se31)+Pv22*(Pv11*Se21+Pv21*Se22+Pv31*Se32)+Pv32*(Pv11*Se31+Pv21*Se32+Pv31*Se33)); |
| | | Sv[0][2]+=We_Av* (Pv31*(Pv11*Se11+Pv21*Se21+Pv31*Se31)+Pv32*(Pv11*Se21+Pv21*Se22+Pv31*Se32)+Pv33*(Pv11*Se31+Pv21*Se32+Pv31*Se33)); |
| | | |
| | | Sv[1][0]+=We_Av* (Pv11*(Pv21*Se11+Pv22*Se21+Pv32*Se31)+Pv21*(Pv21*Se21+Pv22*Se22+Pv32*Se32)+Pv31*(Pv21*Se31+Pv22*Se32+Pv32*Se33)); |
| | | Sv[1][1]+=We_Av* (Pv21*(Pv21*Se11+Pv22*Se21+Pv32*Se31)+Pv22*(Pv21*Se21+Pv22*Se22+Pv32*Se32)+Pv32*(Pv21*Se31+Pv22*Se32+Pv32*Se33)); |
| | | Sv[1][2]+=We_Av* (Pv31*(Pv21*Se11+Pv22*Se21+Pv32*Se31)+Pv32*(Pv21*Se21+Pv22*Se22+Pv32*Se32)+Pv33*(Pv21*Se31+Pv22*Se32+Pv32*Se33)); |
| | | |
| | | Sv[2][0]+=We_Av* (Pv11*(Pv31*Se11+Pv32*Se21+Pv33*Se31)+Pv21*(Pv31*Se21+Pv32*Se22+Pv33*Se32)+Pv31*(Pv31*Se31+Pv32*Se32+Pv33*Se33)); |
| | | Sv[2][1]+=We_Av* (Pv21*(Pv31*Se11+Pv32*Se21+Pv33*Se31)+Pv22*(Pv31*Se21+Pv32*Se22+Pv33*Se32)+Pv32*(Pv31*Se31+Pv32*Se32+Pv33*Se33)); |
| | | Sv[2][2]+=We_Av* (Pv31*(Pv31*Se11+Pv32*Se21+Pv33*Se31)+Pv32*(Pv31*Se21+Pv32*Se22+Pv33*Se32)+Pv33*(Pv31*Se31+Pv32*Se32+Pv33*Se33)); |
| | | // printf("(%f %f %f); (%f %f %f); (%f %f %f)\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]); |
| | | |
| | | } // END FOR JJ |
| | | |
| | | gsl_matrix_set(gsl_Sv, 0,0, Sv[0][0]); |
| | | gsl_matrix_set(gsl_Sv, 0,1, Sv[0][1]); |
| | | gsl_matrix_set(gsl_Sv, 0,2, Sv[0][2]); |
| | | gsl_matrix_set(gsl_Sv, 1,0, Sv[1][0]); |
| | | gsl_matrix_set(gsl_Sv, 1,1, Sv[1][1]); |
| | | gsl_matrix_set(gsl_Sv, 1,2, Sv[1][2]); |
| | | gsl_matrix_set(gsl_Sv, 2,0, Sv[2][0]); |
| | | gsl_matrix_set(gsl_Sv, 2,1, Sv[2][1]); |
| | | gsl_matrix_set(gsl_Sv, 2,2, Sv[2][2]); |
| | | |
| | | // printf("Se= %f, %f, %f\n %f, %f, %f\n %f, %f, %f\n", Se11, Se21, Se31, Se21, Se22, Se32, Se31, Se32, Se33); |
| | | // printf("Pv= %f, %f, %f\n %f, %f, %f\n %f, %f, %f\n", Pv11, Pv21, Pv31, Pv21, Pv22, Pv32, Pv31, Pv32, Pv33); |
| | | // printf("Sv= %f, %f, %f\n %f, %f, %f\n %f, %f, %f\n", Sv[0][0], Sv[0][1], Sv[0][2], Sv[1][0], Sv[1][1], Sv[1][2], Sv[2][0], Sv[2][1], Sv[2][2]); |
| | | |
| | | |
| | | gsl_eigen_symm(gsl_Sv, Sv_eigen, workspace); |
| | | |
| | | // printf("Eigenvalues: %f, %f, %f\n", gsl_vector_get(Sv_eigen, 0),gsl_vector_get(Sv_eigen, 1), gsl_vector_get(Sv_eigen, 2) ); |
| | | // printf("Eigenvalues: %f, %f, %f\n", gsl_matrix_get(evec, 0,0),gsl_matrix_get(evec, 0,1), gsl_matrix_get(evec, 0,2) ); |
| | | |
| | | |
| | | eigenval[0]= gsl_vector_get(Sv_eigen, 0); |
| | | eigenval[1]= gsl_vector_get(Sv_eigen, 1); |
| | | eigenval[2]= gsl_vector_get(Sv_eigen, 2); |
| | | |
| | | qsort(eigenval, 3, sizeof(ts_double), cmpfunc); |
| | | // printf("Eigenvalues: %f, %f, %f\n", eigenval[0], eigenval[1], eigenval[2] ); |
| | | |
| | | |
| | | vtx->energy=(pow(eigenval[0]+eigenval[1],2))*Av; |
| | | |
| | | gsl_matrix_free(gsl_Sv); |
| | | gsl_vector_free(Sv_eigen); |
| | | // gsl_matrix_free(evec); |
| | | gsl_eigen_symm_free(workspace); |
| | | return TS_SUCCESS; |
| | | } |
| | | |