From dba66b4f8a25cb7eb18de5f7ffca97e179d7a5d9 Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@gmail.com> Date: Sat, 19 Oct 2019 11:45:15 +0000 Subject: [PATCH] Moved compression/decompression routines to separate file --- src/energy.c | 39 +++++++++++++++++++++++++++++++++++++-- 1 files changed, 37 insertions(+), 2 deletions(-) diff --git a/src/energy.c b/src/energy.c index 0e57e71..1dae415 100644 --- a/src/energy.c +++ b/src/energy.c @@ -44,9 +44,12 @@ return TS_SUCCESS; }; -/** @brief Calculation of energy of the vertex +/** @brief Calculation of the bending energy of the vertex. * - * Main function that calculates energy of the vertex \f$i\f$. Nearest neighbors (NN) must be ordered in counterclockwise direction for this function to work. + * Main function that calculates energy of the vertex \f$i\f$. Function returns \f$\frac{1}{2}(c_1+c_2-c)^2 s\f$, where \f$(c_1+c_2)/2\f$ is mean curvature, + * \f$c/2\f$ is spontaneous curvature and \f$s\f$ is area per vertex \f$i\f$. + * + * Nearest neighbors (NN) must be ordered in counterclockwise direction for this function to work. * Firstly NNs that form two neighboring triangles are found (\f$j_m\f$, \f$j_p\f$ and common \f$j\f$). Later, the scalar product of vectors \f$x_1=(\mathbf{i}-\mathbf{j_p})\cdot (\mathbf{i}-\mathbf{j_p})(\mathbf{i}-\mathbf{j_p})\f$, \f$x_2=(\mathbf{j}-\mathbf{j_p})\cdot (\mathbf{j}-\mathbf{j_p})\f$ and \f$x_3=(\mathbf{j}-\mathbf{j_p})\cdot (\mathbf{i}-\mathbf{j_p})\f$ are calculated. From these three vectors the \f$c_{tp}=\frac{1}{\tan(\varphi_p)}\f$ is calculated, where \f$\varphi_p\f$ is the inner angle at vertex \f$j_p\f$. The procedure is repeated for \f$j_m\f$ instead of \f$j_p\f$ resulting in \f$c_{tn}\f$. * \begin{tikzpicture}{ @@ -181,6 +184,7 @@ #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); return TS_SUCCESS; @@ -207,3 +211,34 @@ } return TS_SUCCESS; } + +ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old){ + if(fabs(vtx->c)<1e-15) return 0.0; +// printf("was here"); + if(fabs(vesicle->tape->F)<1e-15) return 0.0; + + ts_double norml,ddp=0.0; + ts_uint i; + ts_double xnorm=0.0,ynorm=0.0,znorm=0.0; + /*find normal of the vertex as sum of all the normals of the triangles surrounding it. */ + for(i=0;i<vtx->tristar_no;i++){ + xnorm+=vtx->tristar[i]->xnorm; + ynorm+=vtx->tristar[i]->ynorm; + znorm+=vtx->tristar[i]->znorm; + } + /*normalize*/ + norml=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm); + xnorm/=norml; + ynorm/=norml; + znorm/=norml; + /*calculate ddp, perpendicular displacement*/ + ddp=xnorm*(vtx->x-vtx_old->x)+ynorm*(vtx->y-vtx_old->y)+znorm*(vtx->z-vtx_old->z); + /*calculate dE*/ +// printf("ddp=%e",ddp); + return vesicle->tape->F*ddp; + +} + +void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle){ + triangle->energy=vesicle->tape->xkA0/2.0*pow((triangle->area/vesicle->tlist->a0-1.0),2); +} -- Gitblit v1.9.3