Laboratoire RDP » Biophysics Team » Vascular Networks in Growing Leaves

---------------------------------------------------------------------------

# "Cell division rule" project

(project ID : "leaf")

---------------------------------------------------------------------------

Biophysics Team - RDP ENS Lyon

Developers :

--------------

Jean-Daniel Julien

Download the repository / Rapatrier le dépot

--------------------------------------------

cd "directory that will contain the project's directory"

svn checkout http://forge.cbp.ens-lyon.fr/svn/leaf

Update (download changes from the repository)

--------------------------------------------

cd leaf

svn up

Commit (upload your changes into the repository)

-----------------------------------------------

cd leaf

svn commit

Changing the file tree

----------------------

if you want to add, delete or move individual files in the file-tree of the repository, you have to do that by svn

... cd in the right subdirectory

svn add toto

svn mv toto titi

svn rm titi

... and then make a commit (!)

Useful links

-------------

svn useful commands : http://wiki.greenstone.org/wiki/index.php/Useful_SVN_Commands

formats in a wiki : https://forge.cbp.ens-lyon.fr/redmine/help/wiki_syntax.html

#include "optimize_width.h"

#include "compute_energy_width.h"

#include <nlopt.h>

#include <stdio.h>

#include <stdlib.h>

#include <cmath>

void optimize_width(int vertex_number, double *vertices, int edge_number, int *edges_vertices,

                   double *lengths, double *rest_lengths, double *edges_width, int *periodizer, double *width, double *height){

nlopt_opt opt;

opt=nlopt_create(NLOPT_LN_BOBYQA,edge_number);

nlopt_set_stopval(opt, -HUGE_VAL);

nlopt_set_maxeval(opt, -1.);

nlopt_set_maxtime(opt, -1.);

nlopt_set_xtol_rel(opt, -1.);

nlopt_set_xtol_abs1(opt, 1.e-6);

nlopt_set_ftol_rel(opt, -1.);

nlopt_set_ftol_abs(opt, -1.);

struct parameters_list3 {int vertex_number; double *vertices; int edge_number; int *edges_vertices;

                        double *lengths; double *rest_lengths; int *periodizer; double *width; double *height;};

struct parameters_list3 parameters;

parameters.vertex_number=vertex_number;

parameters.vertices=vertices;

parameters.edge_number=edge_number;

parameters.edges_vertices=edges_vertices;

parameters.lengths=lengths;

parameters.rest_lengths=rest_lengths;

parameters.periodizer=periodizer;

parameters.width=width;

parameters.height=height;

void* pparameters=&parameters;

double width_optimize[edge_number]; double energy;

for (int i=0; i<edge_number;i++){ width_optimize[i]=edges_width[i];}

/// first resolution: find the reference width and height without the anisotropic stress

int reso=nlopt_set_min_objective(opt, compute_energy_width, pparameters);

nlopt_optimize(opt, width_optimize, &energy);

printf("%d\n",reso);

for (int i=0; i<edge_number;i++){edges_width[i]=width_optimize[i];}

nlopt_destroy(opt);

}

#include <math.h>

#include "constants.h"

#include "compute_lengths.h"

#include <stdio.h>

/* this functions computes the lengths of each wall */

void compute_lengths(double *vertices, int *edges_vertices, int edge_number, double *lengths, int *periodizer, double width, double height){

for (int i=0;i<edge_number;i++){

    lengths[i]=sqrt(pow(vertices[2*edges_vertices[2*i]]+width*(double)periodizer[4*i]-vertices[2*edges_vertices[2*i+1]]-width*(double)periodizer[4*i+2],2.)+

                    pow(vertices[2*edges_vertices[2*i]+1]+height*(double)periodizer[4*i+1]-vertices[2*edges_vertices[2*i+1]+1]-height*(double)periodizer[4*i+3],2.));

}

}

#ifndef OPTIMIZE_WIDTH_H_INCLUDED

#define OPTIMIZE_WIDTH_H_INCLUDED

void optimize_width(int vertex_number, double *vertices, int edge_number, int *edges_vertices,

                   double *lengths, double *rest_lengths, double *edges_width, int *periodizer, double *width, double *height);

#endif // OPTIMIZE_WIDTH_H_INCLUDED

#ifndef COMPUTE_LENGTHS_H_INCLUDED

#define COMPUTE_LENGTHS_H_INCLUDED

void compute_lengths(double *vertices, int *edges_vertices, int edge_number, double *lengths, int *periodizer, double width, double height);

#endif // COMPUTE_LENGTHS_H_INCLUDED

#include <string.h>

#include <stdlib.h>

#include <stdio.h>

#include "save_data_2.h"

#include <math.h>

#include "constants.h"

#include <gsl/gsl_matrix.h>

#include <gsl/gsl_blas.h>

#include <gsl/gsl_eigen.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_cblas.h>

#include <gsl/gsl_linalg.h>

#include "gaussian.h"

/* this function saves the tissue and its chemical and mechanical state */

void save_data_2(FILE *f, int vertex_number, double *vertices, int edge_number, int *edges_vertices,

               double *lengths, double *rest_lengths, double *edges_width, int *edges_cells, int *periodizer, double width, double height, int cell_number){

fprintf(f,"edge_number= %i;\n",edge_number);

fprintf(f,"vertex_number= %i;\n",vertex_number);

fprintf(f,"cell_number= %i;\n",cell_number);

fprintf(f,"width_after= %f;\n",width);

fprintf(f,"height_after= %f;\n",height);

fprintf(f,"stress_x= %f;\n",stress_x);

fprintf(f,"stress_y= %f;\n",stress_y);

fprintf(f,"sigma= %f;\n",sigma);

fprintf(f,"vertices=[\n");

for (int i=0;i<vertex_number;i++){

    fprintf(f,"%f %f\n",vertices[2*i],vertices[2*i+1]);}

fprintf(f,"];\n");

fprintf(f,"edges_vertices=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%i %i\n",edges_vertices[2*i],edges_vertices[2*i+1]);

}

fprintf(f,"];\n");

fprintf(f,"lengths=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%f\n",lengths[i]);

}

fprintf(f,"];\n");

fprintf(f,"rest_lengths=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%f\n",rest_lengths[i]);

}

fprintf(f,"];\n");

fprintf(f,"edges_width=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%f\n",edges_width[i]);

}

fprintf(f,"];\n");

fprintf(f,"periodizer=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%i %i %i %i\n",periodizer[4*i],periodizer[4*i+1],periodizer[4*i+2],periodizer[4*i+3]);

}

fprintf(f,"];\n");

fprintf(f,"edges_cells=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%i %i\n",edges_cells[2*i],edges_cells[2*i+1]);

}

fprintf(f,"];\n");

double barycenters[2*cell_number]; for (int i=0; i<2*cell_number;i++){ barycenters[i]=0.;}

double count[cell_number]; for (int i=0; i<cell_number;i++){ count[i]=0.;}

double texture[3*cell_number]; for (int i=0; i<3*cell_number;i++){ texture[i]=0.;}

int box[18]={-1,-1,-1,0,-1,1,0,-1,0,0,0,1,1,-1,1,0,1,1};

double distance=0.; int j_min; double distance_min;

for (int i=0; i<edge_number; i++){

    if (count[edges_cells[2*i]]==0){

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i]]+width*(double)periodizer[4*i])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)periodizer[4*i+1])*lengths[i];

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)periodizer[4*i+2])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)periodizer[4*i+3])*lengths[i];

        count[edges_cells[2*i]]+=2.*lengths[i];

    }

    else {

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i]]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i]+1]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i]]+=1.*lengths[i];

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i]]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i+1]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i]+1]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i+1]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i]]+=1.*lengths[i];

    }

    if (count[edges_cells[2*i+1]]==0){

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i]]+width*(double)periodizer[4*i])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)periodizer[4*i+1])*lengths[i];

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)periodizer[4*i+2])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)periodizer[4*i+3])*lengths[i];

        count[edges_cells[2*i+1]]+=2.*lengths[i];

    }

    else {

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i+1]]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i+1]+1]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i+1]]+=1.*lengths[i];

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i+1]]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i+1]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i+1]+1]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i+1]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i+1]]+=1.*lengths[i];

    }

//    count[edges_cells[2*i]]+=1.;

//    count[edges_cells[2*i+1]]+=1.;

}

for (int i=0; i<cell_number; i++){

    barycenters[2*i]/=count[i];

    barycenters[2*i+1]/=count[i];

}

for (int i=0; i<cell_number;i++){ count[i]=0.;}

for (int i=0; i<edge_number; i++){

    distance_min=100000.;

    for (int j=0; j<9; j++){

        distance=pow(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j]-barycenters[2*edges_cells[2*i+1]],2.)+

                 pow(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j+1]-barycenters[2*edges_cells[2*i+1]+1],2.);

        if (distance<distance_min){

            distance_min=distance;

            j_min=j;

        }

    }

    texture[3*edges_cells[2*i]]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]]);

    texture[3*edges_cells[2*i+1]]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]]);

    texture[3*edges_cells[2*i]+1]+=(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    texture[3*edges_cells[2*i+1]+1]+=(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    texture[3*edges_cells[2*i]+2]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    texture[3*edges_cells[2*i+1]+2]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    count[edges_cells[2*i]]+=1.;

    count[edges_cells[2*i+1]]+=1.;

}

for (int i=0; i<cell_number; i++){

    texture[3*i]/=count[i];

    texture[3*i+1]/=count[i];

    texture[3*i+2]/=count[i];

}

double average_texture[3*cell_number]; for (int i=0; i<3*cell_number; i++){ average_texture[i]=0.;}

for (int i=0; i<cell_number; i++){ count[i]=0.;}

for (int i=0; i<cell_number; i++){

    for (int k=i; k<cell_number; k++){

        distance_min=100000.;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*i]+width*(double)box[2*j]-barycenters[2*k],2.)+

                     pow(barycenters[2*i+1]+height*(double)box[2*j+1]-barycenters[2*k+1],2.);

            if (distance<distance_min){

                distance_min=distance;

            }

        }

        distance_min=sqrt(distance_min);

        average_texture[3*i]+=texture[3*k]*gaussian(distance_min,sigma);

        average_texture[3*i+1]+=texture[3*k+1]*gaussian(distance_min,sigma);

        average_texture[3*i+2]+=texture[3*k+2]*gaussian(distance_min,sigma);

        average_texture[3*k]+=texture[3*i]*gaussian(distance_min,sigma);

        average_texture[3*k+1]+=texture[3*i+1]*gaussian(distance_min,sigma);

        average_texture[3*k+2]+=texture[3*i+2]*gaussian(distance_min,sigma);

        count[i]+=gaussian(distance_min,sigma);

        count[k]+=gaussian(distance_min,sigma);

    }

}

for (int i=0; i<cell_number; i++){

    average_texture[3*i]/=count[i];

    average_texture[3*i+1]/=count[i];

    average_texture[3*i+2]/=count[i];

}

fprintf(f,"barycenters_after=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%f %f\n",barycenters[2*i],barycenters[2*i+1]);

}

fprintf(f,"];\n");

fprintf(f,"texture_after=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%f %f %f\n",texture[3*i],texture[3*i+1],texture[3*i+2]);

}

fprintf(f,"];\n");

fprintf(f,"average_texture_after=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%f %f %f\n",average_texture[3*i],average_texture[3*i+1],average_texture[3*i+2]);

}

fprintf(f,"];\n");

int boundary[cell_number]; for (int i=0; i<cell_number; i++){ boundary[i]=0;}

for (int i=0;i<edge_number;i++){

    if (periodizer[4*i]!=0 || periodizer[4*i+1]!=0 || periodizer[4*i+2]!=0 || periodizer[4*i+3]!=0){

        boundary[edges_cells[2*i]]=1; boundary[edges_cells[2*i+1]]=1;

    }

}

fprintf(f,"boundary=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%i\n",boundary[i]);

}

fprintf(f,"];\n");

}

#include "compute_energy_and_forces.h"

#include "constants.h"

#include "compute_areas.h"

#include "compute_cell_areas.h"

#include "compute_lengths.h"

#include <math.h>

#include <stdio.h>

#include <nlopt.h>

#include <omp.h>

/* this function computes the mechanical energy and its derivatives */

double compute_energy_and_forces(unsigned n, const double *vertices, double *forces, void *parameters){

struct parameters_list {int vertex_number; int edge_number; int *edges_vertices;

                        double *lengths; double *rest_lengths; double *edges_width; int *periodizer; int anisotropy; double ref_width; double ref_height;};

double energy=0;

compute_lengths((double*)vertices, ((parameters_list*)parameters)->edges_vertices, ((parameters_list*)parameters)->edge_number,

                ((parameters_list*)parameters)->lengths, ((parameters_list*)parameters)->periodizer, (double)vertices[2*(((parameters_list*)parameters)->vertex_number)], (double)vertices[2*(((parameters_list*)parameters)->vertex_number)+1]);

if (forces){ for (int i=0;i<2*((parameters_list*)parameters)->vertex_number+2;i++){forces[i]=0.;}}

if (forces){

    forces[2*(((parameters_list*)parameters)->vertex_number)]-=vertices[2*(((parameters_list*)parameters)->vertex_number)+1];

    forces[2*(((parameters_list*)parameters)->vertex_number)+1]-=vertices[2*(((parameters_list*)parameters)->vertex_number)];

}

/* pressure energy */

energy-=vertices[2*(((parameters_list*)parameters)->vertex_number)]*vertices[2*(((parameters_list*)parameters)->vertex_number)+1];

for (int i=0;i<((parameters_list*)parameters)->edge_number;i++){

    /* derivatives of the lengths of the walls */

    if (forces){

        forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i]

                                                                         -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i+2])/

                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]);

        forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i]

                                                                         -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i+2])/

                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]);

        forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+1]

                                                                         -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+3])/

                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]);

        forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+1]

                                                                         -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+3])/

                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]);

        forces[2*((parameters_list*)parameters)->vertex_number]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        ((double)(((parameters_list*)parameters)->periodizer)[4*i]-(double)(((parameters_list*)parameters)->periodizer)[4*i+2])*

                                                                        ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i]

                                                                         -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i+2])/

                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]);

        forces[2*((parameters_list*)parameters)->vertex_number+1]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        ((double)(((parameters_list*)parameters)->periodizer)[4*i+1]-(double)(((parameters_list*)parameters)->periodizer)[4*i+3])*

                                                                        ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])*

                                                                        (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+1]

                                                                         -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+3])/

                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]);

    }

    /* wall elastic energy */

    energy+=0.5*youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*

            pow((((parameters_list*)parameters)->lengths)[i]/(((parameters_list*)parameters)->rest_lengths)[i]-1.,2.);

}

    /* external, anisotropic stretching */

if ((((parameters_list*)parameters)->anisotropy)==1){

if (forces){

    forces[2*(((parameters_list*)parameters)->vertex_number)]-=stress_x*(((parameters_list*)parameters)->ref_height);

    forces[2*(((parameters_list*)parameters)->vertex_number)+1]-=stress_y*(((parameters_list*)parameters)->ref_width);

}

energy-=stress_x*(((parameters_list*)parameters)->ref_height)*(vertices[2*(((parameters_list*)parameters)->vertex_number)]-(((parameters_list*)parameters)->ref_width))+

        stress_y*(((parameters_list*)parameters)->ref_width)*(vertices[2*(((parameters_list*)parameters)->vertex_number)+1]-(((parameters_list*)parameters)->ref_height));

}

return(energy);

}

#include "gaussian.h"

#include <math.h>

double gaussian(double distance, double sigma){

    double result=distance*distance/(2.*sigma*sigma);

    result=exp(-result);

    return result;

}

#include "step.h"

#include <stdio.h>

#include <math.h>

#include "minimize_BFGS.h"

#include "growth_and_chemistry.h"

#include <gsl/gsl_odeiv2.h>

#include <gsl/gsl_errno.h>

#include "compute_energy_and_forces.h"

#include "save_data.h"

#include <cstdlib>

#include "constants.h"

#include <algorithm>

#include "compute_lengths.h"

/* this function performs a time step (growth and chemistry then mechanics) */

void step(double *vertices, int vertex_number, int edge_number, int *edges_vertices, double *lengths, double *rest_lengths, double *edges_width, int *periodizer, double *width, double *height){

/* minimization of energy: move the vertices in order to make the energy as small as possible, i.e. in order to make cells

   as close as possible of their targets */

minimize_BFGS(vertex_number, vertices, edge_number, edges_vertices, lengths, rest_lengths, edges_width, periodizer, width, height);

compute_lengths(vertices, edges_vertices, edge_number, lengths, periodizer, width[0], height[0]);

struct parameters_list2 {int edge_number; double *lengths; double *edges_width;};

struct parameters_list2 parameters2;

parameters2.edge_number=edge_number;

parameters2.lengths=lengths;

parameters2.edges_width=edges_width;

void *pparameters2=&parameters2;

int system_size=edge_number;

double x[system_size];

for (int i=0;i<edge_number;i++){ x[i]=rest_lengths[i];}

gsl_odeiv2_system sys={growth_and_chemistry, NULL, system_size, pparameters2};

gsl_odeiv2_driver *d=gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk2, 1.e-8, 1.e-8, 0.);

int s; double t=0.0;

s=gsl_odeiv2_driver_apply(d, &t, 1.e-7, x);

for (int i=0;i<edge_number;i++){ rest_lengths[i]=x[i];}

gsl_odeiv2_driver_free(d);

}

#include <math.h>

double pi=3.1415926;

double youngs_modulus=300.;

double stress_x=0.;//2.

double stress_y=0.;

double growth_threshold=10.;

double sigma=40;

int growth=0;

#pragma omp threadprivate(stress_x, stress_y, sigma, growth, growth_threshold)

#ifndef SAVE_DATA_2_H_INCLUDED

#define SAVE_DATA_2_H_INCLUDED

void save_data_2(FILE *f, int vertex_number, double *vertices, int edge_number, int *edges_vertices,

               double *lengths, double *rest_lengths, double *edges_width, int *edges_cells, int *periodizer, double width, double height, int cell_number);

#endif // SAVE_DATA_2_H_INCLUDED

#ifndef COMPUTE_ENERGY_AND_FORCES_H_INCLUDED

#define COMPUTE_ENERGY_AND_FORCES_H_INCLUDED

double compute_energy_and_forces(unsigned n, const double *vertices, double *forces, void *parameters);

#endif // COMPUTE_ENERGY_AND_FORCES_H_INCLUDED

#ifndef GAUSSIAN_H_INCLUDED

#define GAUSSIAN_H_INCLUDED

double gaussian(double distance, double sigma);

#endif // GAUSSIAN_H_INCLUDED

#ifndef STEP_H_INCLUDED

#define STEP_H_INCLUDED

void step(double *vertices, int vertex_number, int edge_number, int *edges_vertices, double *lengths, double *rest_lengths, double *edges_width, int *periodizer, double *width, double *height);

#endif // STEP_H_INCLUDED

#ifndef CONSTANTS_H_INCLUDED

#define CONSTANTS_H_INCLUDED

extern double pi;

extern double youngs_modulus;

extern double stress_x;

extern double stress_y;

extern double growth_threshold;

extern double sigma;

extern int growth;

#pragma omp threadprivate(stress_x, stress_y, sigma, growth, growth_threshold)

#endif // CONSTANTS_H_INCLUDED

53.73237637879983

59.62851357493969

86

129

43

41.00605427519755

43.84479514583994

11.11288381815081

6.143416160775232

24.1403013424371

5.227544003018701

30.1660648218456

40.64222965901772

2.426251365215387

3.36138988298003

44.18169811048831

57.860182871652974

44.26863245351591

57.25143624724735

30.05779385705528

28.24536383636932

27.21680454415969

4.637052593595216

33.93035403432599

38.52227721012749

24.3734025425864

22.05447126597037

20.6540676178174

20.38305183805124

39.23489968845102

0.5454952385536798

53.130412457309454

36.98236472617719

17.48416479480939

46.22713497941957

9.590597622999248

29.85658282898794

30.02835746873472

26.54115576464108

14.67166730897258

5.397306638938377

46.95224069568225

9.063338711712447

48.98500227570604

2.475549594269329

45.96502204715647

39.39080378295612

0.2150610878445689

22.44754257525715

14.96796336777424

25.14631127368362

2.80872768095769

45.2358817124404

31.30690565517686

1.951626818952569

15.3424208474777

25.11667364137173

36.37739183218368

21.04400817587831

12.78716971477792

38.28566419249559

51.9718625131994

25.33595327044449

14.5153537428476

53.10917419454713

39.16460333520742

8.706121659048014

14.60912375427341

57.35499733779572

46.22372020161806

16.55464653202851

26.73641157474853

41.66851870804089

2.630386591135188

32.27528433460367

49.07363903777353

24.23969816922651

15.18793030297121

54.67809133530172

33.28650805737531

34.01725007314936

18.3069344296281

44.65626766477443

15.02572660705558

39.39806863038923

8.345942562755766

11.26386898540523

38.11511809476249

20.92351825895694

6.728035099376779

17.5125405894656

47.37846757123232

11.38910923371441

41.2571587330546

47.15701865113333

7.697565141268082

32.19655381556069

23.90021366288214

43.83491267583069

23.66447223319609

5.206701529933286

18.89873770700942

12.50564202067935

10.60583197761373

20.37447099960424

45.34099792801138

19.06395770864364

3.123270759886557

49.02228412571285

17.94994079426626

6.880624613151146

53.17355719326821

50.50223555452406

20.90321203982454

31.84187922596023

20.82343695085001

32.01750906454653

41.92271232146813

31.90706613261577

40.93760583941135

9.217592730034955

28.51015173245191

52.59199247042525

6.959803903539296

12.83996662023328

34.21782816393571

48.94842521543874

43.64335304978598

41.8750472998147

37.73169089323487

0.5851993210129022

53.60511943947328

12.88238734026699

40.16439606485149

13.87062975861086

41.98101543916131

19.46638197842913

28.12384965673833

53.34720490746081

10.73033924657607

20.86039484815298

5.855374868565534

50.26579203945766

4.191681636006709

56.63569714231279

8.919383852200541

50.89323479063303

7.147289239787876

25.37801141592199

3.235864763657833

20.84161600534828

28.60297542757392

14.66922137378002

32.94979934206445

11.70538983818954

32.99391417681286

12.4953449139785

42.69042772349333

30.49154142882092

24.9397573300001

12.40655161134378

34.43963964778495

14.74944849184517

8.55174025523531

42.41374836133627

18.90502784026141

13.3550755311257

32.44039342914015

51.0978327495856

22.17080019983781

55.64088700878654

13.81073090555738

15.58328572034042

52.26595526795432

16.86642863832872

11.11078199584252

46.14142974397822

0

1

0

0

0

9

0

0

8.854079911241516

2.7756533968765713

1

2

0

0

0

44

0

0

3.3217281643406205

3.430072695581895

2

0

0

0

0

61

0

0

3.2916943058805885

2.6668961593506504

3

4

1

0

0

4

0

0

9.121252851595218

3.3340381054359294

5

3

1

0

0

17

0

0

3.6361544732204667

5.094816436869106

4

5

1

0

0

40

0

0

5.820223452760936

1.6117704333223108

6

7

2

0

0

77

0

0

7.223384256457613

0.8902525792086653

7

8

2

0

0

8

0

0

3.1326589432640586

3.2995942249953183

8

6

2

0

0

47

0

0

0.47628536602063226

3.8693366871458186

9

1

3

0

0

60

0

0

9.241735270056887

0.8328752199743118

1

10

3

0

0

9

0

0

4.320193474921223

2.084233441930639

10

9

3

0

0

33

0

0

3.5799148882994896

1.1546833331219073

11

4

4

0

0

19

-1

0

7.228112850300355

2.5231407538546837

3

11

4

0

0

69

0

-1

6.5948981752810445

1.124841730886247

12

13

5

0

0

12

0

1

5.461190814960067

3.6614125389900516

11

12

5

0

0

6

0

0

0.6149227859844556

5.482835720588489

13

11

5

0

0

19

0

1

6.409543734430895

3.1976084732888697

2

12

6

0

0

44

0

0

10.534051479547907

3.209462328992895

11

2

6

0

0

53

0

0

11.17360257979981

1.2558426903041493

14

10

7

0

0

54

0

0

9.83571508310077

1.024934055006501

10

15

7

0

0

37

0

0

6.6135668226516815

2.1610100622891437

15

14

7

0

0

16

0

0

1.7044622767022604

1.7269393196112666

16

8

8

0

0

24

0

0

4.892896758512175

2.5412315161415338

7

16

8

0

0

74

0

0

9.101023060959813

1.151120586641937

0

10

9

0

0

37

0

0

4.550802912338099

2.494666346896691

17

14

10

0

0

11

0

0

4.077633527750127

2.106756008073527

18

17

10

0

0

73

0

0

8.510652311046387

1.1250764035300234

14

18

10

0

0

16

0

0

7.218646203228216

2.075771900641477

19

14

11

0

0

25

0

0

7.114827249465112

1.9059591832016949

17

19

11

0

0

80

0

0

7.24234707219834

2.0723141474224427

16

13

12

0

0

30

0

0

8.16092918436066

1.0146127300272096

12

16

12

0

0

62

0

0

1.5037330534969058

2.729533053292074

20

21

13

0

0

23

1

0

8.930473663224427

1.2416386387610803

22

20

13

0

0

20

0

0

7.559325262226159

1.6451443301594082

21

22

13

0

0

34

1

0

5.710052158186178

2.078204568551323

23

24

14

0

0

85

0

0

6.373959027915761

0.7799771166031012

25

23

14

0

0

38

0

0

1.7732946715648046

2.173580731963505

24

25

14

0

0

29

0

0

7.495085244398764

1.3178970263775502

26

27

15

0

0

22

0

0

7.148616682812731

1.3102700672912762

27

28

15

0

0

71

0

0

5.101701457102572

0.9212679769043179

28

26

15

0

0

45

0

0

3.0098232830004346

1.806485244432378

15

18

16

0

0

26

0

0

8.398146757523122

1.7393652853252513

29

3

17

0

0

31

0

-1

7.671077844244745

2.2370951987414274

5

29

17

0

0

52

0

0

3.598236965098637

5.697151431678375

4

13

18

0

0

19

0

0

6.894279157179725

1.816085499909942

30

4

18

0

0

43

0

0

2.3645037260394544

1.1252827615110104

13

30

18

0

0

57

0

0

6.016612564927007

0.9212192682809177

0

20

20

0

0

61

0

0

3.400237165905518

2.110071038197222

22

0

20

0

0

56

0

0

8.505680286580033

1.193213937768918

31

28

21

0

0

28

-1

0

3.499401768981641

1.7004373258240852

32

31

21

0

0

84

-1

0

5.828912004953271

0.7994025711042237

28

32

21

0

0

72

0

0

3.421148198448322

0.9910089568830227

26

19

22

0

0

25

0

0

0.3756285310714842

2.3354412127296196

19

27

22

0

0

67

0

0

6.027150056137137

1.3479771010724726

33

21

23

0

0

79

0

0

6.398955389727695

0.8905786536815244

20

33

23

0

0

51

0

0

3.799444773086914

1.3115327544561985

34

8

24

0

0

66

0

-1

8.826838177566064

1.4769627121541855

16

34

24

0

0

62

0

0

6.56848456360322

2.5574272552811044

26

14

25

0

0

54

0

0

8.726441934709847

1.075660138674797

35

18

26

0

0

78

0

0

6.586073575659283

0.9977374871095045

15

35

26

0

0

41

0

0

1.7418984998374087

1.5291397101193123

21

23

27

0

0

79

0

0

5.914384296931883

0.8786479172264603

23

26

27

0

0

39

0

0

2.499716102562306

1.523049654651023

26

21

27

0

0

45

0

0

7.936078370151628

1.381559507270994

22

28

28

0

0

34

1

0

8.211840380753

1.6314292128691357

31

22

28

0

0

35

0

0

3.0986246239994975

1.4402180870960941

36

25

29

0

0

36

0

0

1.7070032875719003

1.7119227017909877

24

36