Changes in code after some debuggin

Samo Penic

2021-04-19 efb615c065452b58799e32735647b0559698163f

Changes in code after some debuggin

 src/energy.c

@@ -144,9 +144,9 @@

    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;
        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));
@@ -161,9 +161,10 @@
    Pv31=vertex_normal_x*vertex_normal_z;
    Pv32=vertex_normal_y*vertex_normal_z;




/*    if(vtx->idx==0){
        printf("Vertex normal for vertex %d: %f, %f, %f\n",vtx->idx,vertex_normal_x, 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;
@@ -178,8 +179,11 @@

    //end normalization
//    printf("(%f %f %f)\n", vertex_normal_x, vertex_normal_y, vertex_normal_z);


/*
    if(vtx->idx==0){
        printf("Edge vector for vertex %d (vector %d): %f, %f, %f\n",vtx->idx,jj,edge_vector_x[jj], edge_vector_y[jj], edge_vector_z[jj]);
    }
*/
    it=vtx;
    k=vtx->neigh[jj];
    nei=0;
@@ -235,19 +239,30 @@
    edge_binormal_z[jj]=(edge_normal_x[jj]*edge_vector_y[jj])-(edge_normal_y[jj]*edge_vector_x[jj]);


    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;
    //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]);
    mprod[jj]=lm->xnorm*(lp->ynorm*edge_vector_z[jj]-lp->znorm*edge_vector_y[jj]) - lm->ynorm*(lp->xnorm*edge_vector_z[jj]-lp->znorm*edge_vector_z[jj])+ lm->znorm*(lp->xnorm*edge_vector_y[jj]-lp->ynorm*edge_vector_x[jj]);

    phi[jj]=copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm-1e-10),-mprod[jj])+M_PI;
/*    if(vtx->idx==0){
        printf("Angle PHI vertex %d (angle %d): %f\n",vtx->idx,jj,phi[jj]);
    }
*/
//    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);
//    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);
    he[jj]=temp_length*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];
/*
    if(vtx->idx==0){
        printf("H operator of edge vertex %d (edge %d): %f\n",vtx->idx,jj,he[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;
@@ -266,7 +281,6 @@
//    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]);
@@ -293,8 +307,10 @@
    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] );

    if(vtx->idx==0){
    printf("Eigenvalues: %f, %f, %f\n", eigenval[0], eigenval[1], eigenval[2] );
//    exit(0);
    }

    vtx->energy=(pow(eigenval[0]+eigenval[1],2))*Av;


 src/tape

@@ -25,7 +25,7 @@


#Stretching
stretchswitch=1
stretchswitch=0
xkA0=1.0

####### Polymer (brush) definitions ###########