From b4e2e6dbd5049ae279f0f74838cd21b625148c72 Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@gmail.com> Date: Thu, 11 Oct 2018 10:34:25 +0000 Subject: [PATCH] Changes in plates movement --- src/energy.c | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 files changed, 56 insertions(+), 2 deletions(-) diff --git a/src/energy.c b/src/energy.c index 3173f67..4d27995 100644 --- a/src/energy.c +++ b/src/energy.c @@ -1,3 +1,4 @@ +/* vim: set ts=4 sts=4 sw=4 noet : */ #include<stdlib.h> #include "general.h" #include "energy.h" @@ -43,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}{ @@ -180,7 +184,57 @@ #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; } + + + +ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle){ + int i; + for(i=0;i<vesicle->blist->n;i++){ + attraction_bond_energy(vesicle->blist->bond[i], vesicle->tape->w); + } + return TS_SUCCESS; +} + + +inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w){ + + if(fabs(bond->vtx1->c)>1e-16 && fabs(bond->vtx2->c)>1e-16){ + bond->energy=-w; + } + else { + bond->energy=0.0; + } + 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; + +} -- Gitblit v1.9.3