From efb615c065452b58799e32735647b0559698163f Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Mon, 19 Apr 2021 13:42:52 +0000
Subject: [PATCH] Changes in code after some debuggin
---
src/tape | 2 +-
src/energy.c | 58 +++++++++++++++++++++++++++++++++++++---------------------
2 files changed, 38 insertions(+), 22 deletions(-)
diff --git a/src/energy.c b/src/energy.c
index 5913863..efe5911 100644
--- a/src/energy.c
+++ b/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;
diff --git a/src/tape b/src/tape
index 4e78aff..41bd4da 100644
--- a/src/tape
+++ b/src/tape
@@ -25,7 +25,7 @@
#Stretching
-stretchswitch=1
+stretchswitch=0
xkA0=1.0
####### Polymer (brush) definitions ###########
--
Gitblit v1.9.3