fixed up kenotaxis aligner. Next going to change it to affect the whole tissue

Author: dluke

.gitignore

@@ -1,2 +1,3 @@
 build/
 *.pyc
+.*.swp

src/aligner/external/external_kenotaxis_aligner.cpp

@@ -51,10 +51,15 @@ void ExternalKenotaxisAlign::compute()
     Particle& pi = m_system->get_particle(i);
     if (pi.boundary)
     {
+      // Compute the direction perpendicular to the boundary.
+      // We cheat by using the centre of mass position to deterimine the outside and inside of the tissue.
+
+      // vector (X, Y, Z) points from centre of mass to particle i
       double force_sign = 1.0;
       double X = pi.x - xcm,  Y = pi.y - ycm,  Z = pi.z - zcm;
       double len_R = sqrt(X*X + Y*Y + Z*Z);
-      X /= len_R;  Y /= len_R;  Y /= len_R;
+      X /= len_R;  Y /= len_R;  Z /= len_R;
+      // get relative vectors of neighbouring boundary particles
       Particle& pj = m_system->get_particle(pi.boundary_neigh[0]);
       Particle& pk = m_system->get_particle(pi.boundary_neigh[1]);
       double xji = pj.x - pi.x, yji = pj.y - pi.y, zji = pj.z - pi.z;
@@ -63,10 +68,11 @@ void ExternalKenotaxisAlign::compute()
       double xki = pk.x - pi.x, yki = pk.y - pi.y, zki = pk.z - pi.z;
       double len_ki = sqrt(xki*xki + yki*yki + zki*zki);
       xki /= len_ki;  yki /= len_ki;  zki /= len_ki;
-      double x = -(xji+xki), y = -(yji+yki), z = -(zji+zki);
+      // compute vector perpendicular to r_kj
+      double x = -(yji-yki), y = (xji-xki), z = 0.0; // give up with three dimensions, z=0
       double len_r = sqrt(x*x + y*y + z*z);
       x /= len_r;  y /= len_r;  z /= len_r;
-      // compute dot product with radius the vector connecting center of mass and pi
+      // compute dot product with vector connecting center of mass and (x,y,z)
       if ((x*X + y*Y + z*Z) < 0.0)
         force_sign = -1.0;
       double J_factor = force_sign*J;

src/aligner/external/external_shape_aligner.cpp

@@ -27,6 +27,8 @@
  * \brief Declaration of ExternalShapeAlign class
  */ 
 
+// Problem with shape alignment is that orientation can flip like a nematic
+
 
 #include "external_shape_aligner.hpp"
 
@@ -61,8 +63,8 @@ void ExternalShapeAlign::compute()
         }
         A /= V.dual.size();  B /= V.dual.size();  C /= V.dual.size();
         // Maximum eigenvalue is 
-        double l1 = 0.5*(A+C + sqrt((A-C)*(A-C)+4*B*B));
-        double l2 = 0.5*(A+C - sqrt((A-C)*(A-C)+4*B*B));
+        double l1 = 0.5*(A+C + sqrt((A+C)*(A+C)+4*B*B));
+        double l2 = 0.5*(A+C - sqrt((A+C)*(A+C)+4*B*B));
         double lambda = (l1 > l2) ? l1 : l2;
         // compute corresponding eigen vector
         ax = B;

utils/pvmodel/analyse_cells.py

@@ -78,7 +78,7 @@ def process_commands():
             help='Mirrors exclude_boundary flag in samso')
 
 
-    # Want a sub parser for simple operations on mesh which don't involve stress calcualtion
+    # Want a sub parser for simple operations on mesh which don't involve stress calculation
     # Examples. Calculating the structure tensor, calculating cell neighbour topology
     # Associated with the 'Property' senario
     meshparser = subparsers.add_parser('mesh', help='Analysis on the mesh')

utils/pvmodel/cellmesh.py

@@ -620,7 +620,6 @@ def _parts(self):
         # still ignoring kinetic part for now
         stress = dict([(i, self.stress[i]) for i in self.clist])
         evalues, evectors = {}, {}
-        #print stress
         for i in self.clist:
             st = stress[i]
 
@@ -1290,8 +1289,6 @@ def makecellparts(self):
             stress += struct 
             self.cellparts[vhi] = stress
 
-            #print vhi
-            #print stress
 
 
         self.structure = structure
@@ -1322,7 +1319,6 @@ def vvext(rings, last_ring, adj):
 
                 rings.append(last_ring)
                 adj -= 1
-                #print adj
                 vvext(rings, last_ring, adj)
 
         vvext(rings, last_ring, adj)

utils/pvmodel/senario.py

@@ -407,7 +407,7 @@ def _operate(self):
         pv.calculate_forces()
 
         pv.makecellparts()
-        pv.on_centres()
+        pv.on_centres(args.adj)
 
         if self.fid == self.ref_index:
             self.ref_structure = {}
@@ -416,9 +416,10 @@ def _operate(self):
         if args.test:
             self.test_convergence(pv)
 
-        # for real
-        adjlist = [0, 4]
-        #self.calculate_U(pv, adjlist)
+        # for real 
+        print 'stress averaging adj = ', args.adj
+        adjlist = [0, args.adj]
+        self.calculate_U(pv, adjlist)
 
         # pickle useful datas
         self.pkr.update(pv, outnum)
@@ -428,7 +429,7 @@ def _operate(self):
 
     def calculate_U(self, pv, adjlist):
         self.adjavg_structure(pv, adjlist)
-        pv.structure_xx
+        #pv.structure_xx
 
 
     def test_convergence(self, pv):