Laboratoire RDP » Biophysics Team » Cell Division Rule

#include "compute_areas.h"

#include "compute_cell_areas.h"

#include "compute_barycenters.h"

#include "compute_shapes.h"

#include "compute_lengths.h"

#include "initialize.h"

#include "constants.h"

#include <stdio.h>

#include "compute_cell_walls.h"

/* this function initializes the various lists used in the program */

void initialize(double *vertices, int vertices_number, int *cells_vertices, int *vertices_number_in_each_cell, int cell_number, int *wall_vertices, int wall_number, int *cell_walls,

                double *barycenters, double *shapes, double *targets, double *lengths, double *target_lengths, double *areas, double *cell_areas, double *forces, double *stress_periclinal, double *stress_anticlinal){

int i;

compute_barycenters(vertices,vertices_number,cells_vertices,vertices_number_in_each_cell,cell_number,barycenters);

compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes);

compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas);

compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas);

compute_cell_walls(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, cell_walls, wall_vertices, wall_number);

compute_lengths(vertices, wall_vertices, wall_number, lengths);

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

    targets[3*i]=initial_growth*shapes[3*i]; targets[3*i+1]=initial_growth*shapes[3*i+1]; targets[3*i+2]=initial_growth*shapes[3*i+2];

    stress_periclinal[3*i]=0.; stress_periclinal[3*i+1]=0.; stress_periclinal[3*i+2]=0.;

}

for (i=0;i<wall_number;i++){

    target_lengths[i]=lengths[i]; stress_anticlinal[3*i]=0.; stress_anticlinal[3*i+1]=0.; stress_anticlinal[3*i+2]=0.;

}

}

#include "compute_energy_and_forces.h"

#include "minimize_steepest_descent.h"

#include <nlopt.h>

#include <stdio.h>

#include <stdlib.h>

#include <cmath>

/* this function minimizes the mechanical energy. It is an alternative to the BFGS algorithm, just for debugging */

void minimize_steepest_descent(int vertices_number, int cell_number, double *vertices, int *cells_vertices, int *vertices_number_in_each_cell, double *barycenters, double *shapes, double *targets, double *areas, double *cell_areas, double *target_areas,

                   int *wall_vertices, int wall_number, double *lengths, double *target_lengths){

double step_size=1.e-2;

double precision=1.e-2;

struct parameters_list {int vertices_number; int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *barycenters; double *shapes; double *targets; double *areas; double *cell_areas; double *target_areas;

                        int *wall_vertices; int wall_number; double *lengths; double *target_lengths;};

struct parameters_list parameters;

parameters.vertices_number=vertices_number;

parameters.cell_number=cell_number;

parameters.cells_vertices=cells_vertices;

parameters.vertices_number_in_each_cell=vertices_number_in_each_cell;

parameters.barycenters=barycenters;

parameters.shapes=shapes;

parameters.targets=targets;

parameters.areas=areas;

parameters.cell_areas=cell_areas;

parameters.target_areas=target_areas;

parameters.wall_vertices=wall_vertices;

parameters.wall_number=wall_number;

parameters.lengths=lengths;

parameters.target_lengths=target_lengths;

void* pparameters=&parameters;

double forces[2*vertices_number];

double energy=compute_energy_and_forces(0, vertices, forces, pparameters);

double max_step=0.; for (int i=0;i<2*vertices_number;i++){ if (fabs(forces[i])>max_step){ max_step=fabs(forces[i]);}}

int step_number=0;

while (max_step>precision){

    step_number++;

    for (int i=0;i<2*vertices_number;i++){ vertices[i]-=step_size*forces[i];}

    energy=compute_energy_and_forces(0, vertices, forces, pparameters);

    max_step=0;

    for (int i=0;i<2*vertices_number;i++){ if (fabs(forces[i])>max_step){ max_step=fabs(forces[i]);}}

}

printf("opt %f %d\n",energy,step_number);

}

#ifndef COMPUTE_LENGTHS_H_INCLUDED

#define COMPUTE_LENGTHS_H_INCLUDED

void compute_lengths(double *vertices, int *wall_vertices, int wall_number, double *lengths);

#endif // COMPUTE_LENGTHS_H_INCLUDED

#include "compute_energy_and_forces.h"

#include "constants.h"

#include "compute_areas.h"

#include "compute_cell_areas.h"

#include "compute_shapes.h"

#include "compute_barycenters.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 (basically only math equations) */

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

struct parameters_list {int vertices_number; int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *barycenters; double *shapes; double *targets; double *areas; double *cell_areas;

                        int *wall_vertices; int wall_number; double *lengths; double *target_lengths; int *boundary; double *tissue_target; int anisotropy; double *relative_pressure;};

int j,previous; double N,energy=0;

compute_lengths((double*)vertices, ((parameters_list*)parameters)->wall_vertices, ((parameters_list*)parameters)->wall_number, ((parameters_list*)parameters)->lengths);

compute_barycenters((double*)vertices, ((parameters_list*)parameters)->vertices_number, ((parameters_list*)parameters)->cells_vertices,

                     ((parameters_list*)parameters)->vertices_number_in_each_cell, ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->barycenters);

compute_areas((double*)vertices, ((parameters_list*)parameters)->cells_vertices, ((parameters_list*)parameters)->vertices_number_in_each_cell,

              ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->areas);

compute_cell_areas((double*)vertices, ((parameters_list*)parameters)->cells_vertices, ((parameters_list*)parameters)->vertices_number_in_each_cell,

              ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->areas, ((parameters_list*)parameters)->cell_areas);

compute_shapes((double*)vertices, ((parameters_list*)parameters)->cells_vertices, ((parameters_list*)parameters)->vertices_number_in_each_cell,

               ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->barycenters, ((parameters_list*)parameters)->shapes);

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

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

//    if (vertices[2*i+1]>7.){ forces[2*i+1]+=0.01*sinh(vertices[2*i+1]-7.);}

//    if (vertices[2*i+1]<-7.){ forces[2*i+1]+=0.01*sinh(vertices[2*i+1]+7.);}

//    if (vertices[2*i+1]>7.){ energy+=0.01*(cosh(vertices[2*i+1]-7.)-1.);}

//    if (vertices[2*i+1]<-7.){ energy+=0.01*(cosh(vertices[2*i+1]+7.)-1.);}

//}

double ablation;

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

    ablation=(double)(i!=ablated_cell);

    if (forces){

    for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){

        previous=j-1; if (j==(((parameters_list*)parameters)->vertices_number_in_each_cell)[i]){ previous=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1]-1;}

        /* derivatives of the area in the pressure energy */

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]-=0.5*ablation*pressure*(((parameters_list*)parameters)->relative_pressure)[i]*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-=0.5*ablation*pressure*(((parameters_list*)parameters)->relative_pressure)[i]*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]]);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]]-=0.5*ablation*pressure*(((parameters_list*)parameters)->relative_pressure)[i]*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]-=0.5*ablation*pressure*(((parameters_list*)parameters)->relative_pressure)[i]*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]);

    }

    }

    /* pressure energy */

    energy+=-pressure*(((parameters_list*)parameters)->relative_pressure)[i]*ablation*(((parameters_list*)parameters)->cell_areas)[i];

}

double boundary;

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

//    for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){

//        previous=j-1; if (j==(((parameters_list*)parameters)->vertices_number_in_each_cell)[i]){ previous=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1]-1;}

//        boundary=(double)((((parameters_list*)parameters)->boundary)[j]);

//        if (forces){

//            forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]-=boundary*0.5*pressure*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]);

//            forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-=boundary*0.5*pressure*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]]);

//            forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]]-=boundary*0.5*pressure*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]);

//            forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]-=boundary*0.5*pressure*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]);

//        }

//

//        energy-=boundary*0.5*pressure*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-

//                                                vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]);

//    }

}

double trace;

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

    trace=(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->targets)[3*i+1];

    N=(double)((((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1])-(double)((((parameters_list*)parameters)->vertices_number_in_each_cell)[i]);

    ablation=(double)(i!=ablated_cell);

    if (forces){

    for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){

        previous=j-1; if (j==(((parameters_list*)parameters)->vertices_number_in_each_cell)[i]){ previous=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1]-1;}

        /* derivatives of the shape matrix */

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]+=ablation*shear_modulus*4.0*((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i])

                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]-(((parameters_list*)parameters)->barycenters)[2*i])*

                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]+=ablation*shear_modulus*4.0*((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2])

                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-(((parameters_list*)parameters)->barycenters)[2*i+1])*

                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]+=ablation*shear_modulus*4.0*((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1])

                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-(((parameters_list*)parameters)->barycenters)[2*i+1])*

                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]+=ablation*shear_modulus*4.0*((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2])

                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]-(((parameters_list*)parameters)->barycenters)[2*i])*

                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);

        /* derivatives of the area */

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]+=ablation*shear_modulus*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1])*

                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+

                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]+=ablation*shear_modulus*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]])*

                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+

                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]]+=ablation*shear_modulus*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1])*

                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+

                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);

        forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]+=ablation*shear_modulus*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]])*

                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+

                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);

    }

    }

    /* shape energy, contribution from the shear modulus */

    energy+=ablation*shear_modulus*(((parameters_list*)parameters)->cell_areas)[i]*(pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)

                        +pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

                        +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(trace*trace);//(((parameters_list*)parameters)->cell_areas)[i];

}

    /* contribution of the first lame parameter, this term was neglected in the final version of the paper as it has no impact */

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

//    min=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i]; max=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];

//    trace=(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->targets)[3*i+1];

//    N=(double)(max)-(double)(min);

//    if (forces){

//    for (j=0;j<max-min;j++){ previous=j-1; if (j==0){ previous=max-min-1;}

//        /* derivatives of the shape matrix */

//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]]+=lame_first_periclinal*4.0*(((parameters_list*)parameters)->cell_areas)[i]*

//                                                    ((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1])

//                                                    *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]]-(((parameters_list*)parameters)->barycenters)[2*i])

//                                                    /(trace*trace*N);

//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1]+=lame_first_periclinal*4.0*(((parameters_list*)parameters)->cell_areas)[i]*

//                                                    ((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1])

//                                                    *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1]-(((parameters_list*)parameters)->barycenters)[2*i+1])

//                                                    /(trace*trace*N);

//        /* derivatives of the area */

//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]]+=lame_first_periclinal*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]+1])*

//                                    pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

//                                    /(2.*trace*trace);

//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1]+=lame_first_periclinal*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]])*

//                                    pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

//                                    /(2.*trace*trace);

//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]]+=lame_first_periclinal*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1])*

//                                    pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

//                                    /(2.*trace*trace);

//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]+1]+=lame_first_periclinal*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]])*

//                                    pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

//                                    /(2.*trace*trace);

//    }

//    }

//    /* area energy */

//    energy+=lame_first_periclinal*(((parameters_list*)parameters)->cell_areas)[i]*

//                pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)

//                /(trace*trace);

}

    /* energy of the anticlinal walls */

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

//    /* derivatives of the lengths of the walls */

//    if (forces){

//    forces[2*(((parameters_list*)parameters)->wall_vertices)[2*i]]+=lame_first_anticlinal/*(((parameters_list*)parameters)->wall_out)[i]*/*2.*((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->target_lengths)[i])*

//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]])/

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

//    forces[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]]+=lame_first_anticlinal/*(((parameters_list*)parameters)->wall_out)[i]*/*2.*((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->target_lengths)[i])*

//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]])/

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

//    forces[2*(((parameters_list*)parameters)->wall_vertices)[2*i]+1]+=lame_first_anticlinal/*(((parameters_list*)parameters)->wall_out)[i]*/*2.*((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->target_lengths)[i])*

//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]+1]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]+1])/

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

//    forces[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]+1]+=lame_first_anticlinal/*(((parameters_list*)parameters)->wall_out)[i]*/*2.*((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->target_lengths)[i])*

//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]+1]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]+1])/

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

//    }

//    /* walls energy */

//    /* the 2. factor comes from the fact that each wall is counted twice in the sum over the cells */

//    energy+=lame_first_anticlinal*pow((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->target_lengths)[i],2.)/(((parameters_list*)parameters)->target_lengths)[i];

}

    /* contribution from the anisotropic stress that mimics curvature */

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

    double tissue_area=0.; int boundary_vertices_number=0;

    double tissue_shape[2]={0.,0.}; double tissue_barycenter[2]={0.,0.};

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

        for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){

            previous=j-1; if (j==(((parameters_list*)parameters)->vertices_number_in_each_cell)[i]){ previous=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1]-1;}

            boundary_vertices_number+=(((parameters_list*)parameters)->boundary)[j];

            boundary=(double)((((parameters_list*)parameters)->boundary)[j]);

            tissue_area+=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-

                                   vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1])/2.;

            tissue_barycenter[0]+=boundary*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]];

            tissue_barycenter[1]+=boundary*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1];

        }

    }

    tissue_barycenter[0]/=boundary_vertices_number;

    tissue_barycenter[1]/=boundary_vertices_number;

    for (int i=0;i<(((parameters_list*)parameters)->vertices_number_in_each_cell)[((parameters_list*)parameters)->cell_number];i++){

        boundary=(double)((((parameters_list*)parameters)->boundary)[i]);

        tissue_shape[0]+=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[i]]-tissue_barycenter[0])*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[i]]-tissue_barycenter[0]);

        tissue_shape[1]+=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[i]+1]-tissue_barycenter[1])*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[i]+1]-tissue_barycenter[1]);

    }

    tissue_shape[0]/=boundary_vertices_number;

    tissue_shape[1]/=boundary_vertices_number;

    trace=(((parameters_list*)parameters)->tissue_target)[0]+(((parameters_list*)parameters)->tissue_target)[1];

    if (forces){

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

            for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){

                previous=j-1; if (j==(((parameters_list*)parameters)->vertices_number_in_each_cell)[i]){ previous=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1]-1;}

                boundary=(double)((((parameters_list*)parameters)->boundary)[j]);

                forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]-=boundary*tissue_area*first_shape_derived_stress*2.*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]-tissue_barycenter[0])/

                                                                               (trace*boundary_vertices_number);

                forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-=boundary*tissue_area*second_shape_derived_stress*2.*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-tissue_barycenter[1])/

                                                                                 (trace*boundary_vertices_number);

                forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]-=boundary*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1])*

                                                                               (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+

                                                                               second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);

                forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]]-=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1])*

                                                                                      (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+

                                                                                       second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);

                forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]])*

                                                                                 (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+

                                                                                 second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);

                forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]-=boundary*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]])*

                                                                                        (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+

                                                                                        second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);

            }

        }

    }

    energy-=tissue_area*(first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+

                         second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/trace;

}

return(energy);

}

#ifndef CELL_DIVISION2_H_INCLUDED

#define CELL_DIVISION2_H_INCLUDED

void cell_division2(double *vertices, int vertices_number, int *cells_vertices, int *vertices_number_in_each_cell, int cell_number, int *wall_vertices, int wall_number,

          double *barycenters, double *areas, double *cell_areas, double *lengths, double *target_lengths, double *shapes, double *targets,

          int dividing_cell, int *boundary, int *border, int adjacent_cells_number, int adjacent_cell1, int adjacent_cell2,

          int *daughter_cells_vertices_number, int *cell1, int *cell2, int *wall1, int *wall2);

#endif // CELL_DIVISION2_H_INCLUDED

#include <math.h>

int step_number=1;//51;

int max_cell_number=200;//400; //1000000

double tau1=10.;

double pi=3.1415926;

double initial_growth=1.;

double threshold_area=3.;

double threshold_division=0.;

double shear_modulus=1.;

double lame_first_periclinal=0.;//0.1;

double lame_first_anticlinal=0.0;//0.02;

double pressure=1.e-2;//6.e-2;//12.e-2;

double first_shape_derived_stress=3.e-2;// 3.e-2

double second_shape_derived_stress=3.e-2;

double noise_amplitude=0.1;

int ablated_cell=-1;

int iteration=0;

int iteration_ablation=1000000; //50

double growth_differential=1.; //>=0

#pragma omp threadprivate(threshold_area,threshold_division,shear_modulus,lame_first_periclinal,lame_first_anticlinal,ablated_cell,iteration,growth_differential,first_shape_derived_stress,second_shape_derived_stress)

#include <math.h>

#include "constants.h"

#include "compute_cell_walls.h"

#include <stdio.h>

/* this function computes the list of the indices of the wall of each cell */

void compute_cell_walls(double *vertices, int *cells_vertices, int *vertices_number_in_each_cell,

                    int cell_number, int *cell_walls, int *wall_vertices, int wall_number){

int i,j,k,previous;

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

    for (j=vertices_number_in_each_cell[i];j<vertices_number_in_each_cell[i+1];j++){ previous=j-1; if (j==vertices_number_in_each_cell[i]){previous=vertices_number_in_each_cell[i+1]-1;}

        for (k=0;k<wall_number;k++){

            if ((cells_vertices[previous]==wall_vertices[2*k] && cells_vertices[j]==wall_vertices[2*k+1]) || (cells_vertices[previous]==wall_vertices[2*k+1] && cells_vertices[j]==wall_vertices[2*k])){

                cell_walls[j]=k;

            }

        }

    }

}

}

#include "compute_stress_anticlinal.h"

#include "constants.h"

#include <stdio.h>

/* this function computes the stress on a cell by interpolating the forces on its vertices */

void compute_stress_anticlinal(int wall_index, double *vertices, int *vertices_number_in_each_cell, int cell_number, double *stress_anticlinal,

                                int *wall_vertices, int *cell_walls, double *lengths, double *target_lengths){

double projection[3];

projection[0]=(vertices[2*wall_vertices[2*wall_index]]-vertices[2*wall_vertices[2*wall_index+1]])*(vertices[2*wall_vertices[2*wall_index]]-vertices[2*wall_vertices[2*wall_index+1]]);

projection[1]=(vertices[2*wall_vertices[2*wall_index]+1]-vertices[2*wall_vertices[2*wall_index+1]+1])*(vertices[2*wall_vertices[2*wall_index]+1]-vertices[2*wall_vertices[2*wall_index+1]+1]);

projection[2]=(vertices[2*wall_vertices[2*wall_index]+1]-vertices[2*wall_vertices[2*wall_index+1]+1])*(vertices[2*wall_vertices[2*wall_index]]-vertices[2*wall_vertices[2*wall_index+1]]);

projection[0]/=lengths[wall_index]*lengths[wall_index];

projection[1]/=lengths[wall_index]*lengths[wall_index];

projection[2]/=lengths[wall_index]*lengths[wall_index];

stress_anticlinal[3*wall_index]=lame_first_anticlinal*(lengths[wall_index]-target_lengths[wall_index])*projection[0]/target_lengths[wall_index];

stress_anticlinal[3*wall_index+1]=lame_first_anticlinal*(lengths[wall_index]-target_lengths[wall_index])*projection[1]/target_lengths[wall_index];

stress_anticlinal[3*wall_index+2]=lame_first_anticlinal*(lengths[wall_index]-target_lengths[wall_index])*projection[2]/target_lengths[wall_index];

}

#ifndef STEP_H_INCLUDED

#define STEP_H_INCLUDED

void step(double *vertices, int vertices_number, int *cells_vertices, int *vertices_number_in_each_cell, int cell_number, int *wall_vertices, int wall_number, int *cell_walls,

          double *barycenters, double *shapes, double *targets, double *lengths, double *target_lengths, double *areas, double *cell_areas, double *forces, double *stress_periclinal, double *stress_anticlinal,

           int *boundary, double *growth_rate, double *relative_pressure);

#endif // STEP_H_INCLUDED

#include "compute_energy_and_forces.h"

#include "minimize_BFGS.h"

#include <nlopt.h>

#include <stdio.h>

#include <stdlib.h>

#include <cmath>

/* this function minimizes the mechanical energy */

void minimize_BFGS(int vertices_number, int cell_number, double *vertices, int *cells_vertices, int *vertices_number_in_each_cell, double *barycenters, double *shapes, double *targets, double *areas, double *cell_areas,

                   int *wall_vertices, int wall_number, double *lengths, double *target_lengths, int *boundary, double *relative_pressure){

nlopt_opt opt;

opt=nlopt_create(NLOPT_LD_LBFGS,2*vertices_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_list {int vertices_number; int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *barycenters; double *shapes; double *targets; double *areas; double *cell_areas;

                        int *wall_vertices; int wall_number; double *lengths; double *target_lengths; int *boundary; double *tissue_target; int anisotropy; double *relative_pressure;};

double min_energy; double tissue_target[2]={1.,1.};

struct parameters_list parameters;

parameters.vertices_number=vertices_number;

parameters.cell_number=cell_number;

parameters.cells_vertices=cells_vertices;

parameters.vertices_number_in_each_cell=vertices_number_in_each_cell;

parameters.barycenters=barycenters;

parameters.shapes=shapes;

parameters.targets=targets;

parameters.areas=areas;

parameters.cell_areas=cell_areas;

parameters.wall_vertices=wall_vertices;

parameters.wall_number=wall_number;

parameters.lengths=lengths;

parameters.target_lengths=target_lengths;

parameters.boundary=boundary;

parameters.tissue_target=tissue_target;

parameters.anisotropy=0;

parameters.relative_pressure=relative_pressure;

void* pparameters=&parameters;

//double forces[2*vertices_number];

//double test=compute_energy_and_forces(0, vertices, forces, pparameters);

////for (int i=0; i<vertices_number;i++){ printf("%f %f %d %d\n",forces[2*i],forces[2*i+1],i,vertices_number);} printf("test\n");

//double test2;

//for (int i=0;i<2*vertices_number;i++){vertices[i]+=1.e-8;

//test2=compute_energy_and_forces(0, vertices, forces, pparameters);

//test2-=test;

//test2*=1.e8;

//vertices[i]-=1.e-8;

//printf("%f %f %f\n",test2,forces[i],(test2-forces[i])/forces[i]);

//}printf("\n");

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

nlopt_optimize(opt, vertices, &min_energy);

tissue_target[0]=0.; tissue_target[1]=0.;

int boundary_vertices_number=0.; double tissue_barycenter[2]={0.,0.};

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

    tissue_barycenter[0]+=((double)boundary[i])*vertices[2*cells_vertices[i]];

    tissue_barycenter[1]+=((double)boundary[i])*vertices[2*cells_vertices[i]+1];

    boundary_vertices_number+=boundary[i];

}

tissue_barycenter[0]/=boundary_vertices_number;

tissue_barycenter[1]/=boundary_vertices_number;

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

    tissue_target[0]+=((double)boundary[i])*(vertices[2*cells_vertices[i]]-tissue_barycenter[0])*(vertices[2*cells_vertices[i]]-tissue_barycenter[0]);

    tissue_target[1]+=((double)boundary[i])*(vertices[2*cells_vertices[i]+1]-tissue_barycenter[1])*(vertices[2*cells_vertices[i]+1]-tissue_barycenter[1]);

}

tissue_target[0]/=boundary_vertices_number;

tissue_target[1]/=boundary_vertices_number;

parameters.tissue_target=tissue_target;

parameters.anisotropy=1;

pparameters=&parameters;

reso=nlopt_set_min_objective(opt, compute_energy_and_forces, pparameters);

nlopt_optimize(opt, vertices, &min_energy);

nlopt_destroy(opt);

}

#ifndef COMPUTE_barycenters_H_INCLUDED

#define COMPUTE_barycenters_H_INCLUDED

void compute_barycenters(double *vertices, int vertices_number, int *cells_vertices,

                     int *vertices_number_in_each_cell, int cell_number, double *barycenters);

#endif // COMPUTE_barycenters_H_INCLUDED

#include <math.h>

#include "cell_division.h"

#include "constants.h"

#include <stdio.h>

#include <stdlib.h>

#include <omp.h>

/* this function:   finds the division direction

                    computes the position of the new vertices

                    saves the indices and the vertices of the two daughter cells and of the two walls which are cut by the new one (these data are usefull for cell_division2.cpp) */

void cell_division(double *vertices, int vertices_number, int *cells_vertices, int *vertices_number_in_each_cell, int cell_number,

          double *barycenters, int dividing_cell, int *boundary, double *stress, double *shapes,

          int *adjacent_cells_number, int *adjacent_cell1, int *adjacent_cell2, int *daughter_cells_vertices_number, int *cell1, int *cell2, int *wall1, int *wall2){

///bias division

//double width=3; double dilation=1.;

int j;

/* vector of the division line */

double division_direction[2];

/* positions of new vertices */

double new_vertex1[2]={0}; double new_vertex2[2]={0};

double eigenvalue=0; double norm;

if (threshold_division>0.0){ /* probabilistic division: the probability to divide a cell following the vector n is proportionnal to

                        exp(n^t D n / D0) where D is the deviatoric stress (S-tr(S)/2) and D0 the threshold over which we start felling the stress */

    double anisotropic_stress[2];

    anisotropic_stress[0]=(0.5*stress[3*dividing_cell]-0.5*stress[3*dividing_cell+1])/threshold_division;

    anisotropic_stress[1]=0.5*stress[3*dividing_cell+2]/threshold_division;

    double probability[20];

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

        probability[i]=exp(cos(i*pi/19.0)*anisotropic_stress[0]+sin(i*pi/19.0)*anisotropic_stress[1]);

        if (probability[i]==HUGE_VAL){ printf("division_h\n"); goto division_highest_tension;}

    }

    probability[0]/=2.; probability[19]/=2.;

    for (int i=1;i<20;i++){probability[i]+=probability[i-1];}

    for (int i=0;i<20;i++){probability[i]/=probability[19];}

    double random=(double)rand()/(double)RAND_MAX; //printf("random %f\n",random);

    int random2=rand()%2;

    if (random<=probability[0]){ //printf("%i\n",0);

        division_direction[0]=1;

        division_direction[1]=0;}

    for (int i=1;i<20;i++){//printf("%f %f\n",probability[i-1],probability[i]);

        if (random<=probability[i] && random>probability[i-1]){

            if(random2==0){//printf("%i %i\n",i,0);

                division_direction[0]=cos(i*0.5*pi/19.0);

                division_direction[1]=sin(i*0.5*pi/19.0);}

            if(random2==1){//printf("%i %i\n",i,1);

                division_direction[0]=cos(i*0.5*pi/19.0);

                division_direction[1]=-sin(i*0.5*pi/19.0);}

        }

    }

}

if (threshold_division==0.0){ /* division following the maximal stress = particular case of the previous rule */

    division_highest_tension:; //printf("%f %d %d\n",stress[3*dividing_cell+2],dividing_cell,cell_number);

    if (stress[3*dividing_cell+2]!=0){

    eigenvalue=(stress[3*dividing_cell]+stress[3*dividing_cell+1])/2.;

    eigenvalue+=0.5*sqrt((stress[3*dividing_cell]-stress[3*dividing_cell+1])*(stress[3*dividing_cell]-stress[3*dividing_cell+1])+

                    4.*stress[3*dividing_cell+2]*stress[3*dividing_cell+2]);

    division_direction[0]=stress[3*dividing_cell+2];

    division_direction[1]=-stress[3*dividing_cell]+eigenvalue;

    norm=sqrt(division_direction[0]*division_direction[0]+division_direction[1]*division_direction[1]);

    division_direction[0]/=norm;

    division_direction[1]/=norm;}

    else if (stress[3*dividing_cell]>stress[3*dividing_cell+1]){

    division_direction[0]=1;

    division_direction[1]=0;}

    else if (stress[3*dividing_cell]<stress[3*dividing_cell+1]){

    division_direction[0]=0;

    division_direction[1]=1;}

    else if (stress[3*dividing_cell]==stress[3*dividing_cell+1]){

    division_direction[0]=rand()/(double)RAND_MAX;

    division_direction[1]=sqrt(1.-division_direction[0]*division_direction[0]);}

}

double temp_shape[3];

temp_shape[0]=shapes[3*dividing_cell];

temp_shape[1]=shapes[3*dividing_cell+1];

temp_shape[2]=shapes[3*dividing_cell+2];

//if ((barycenters[2*dividing_cell]*barycenters[2*dividing_cell])<(width*width)){ temp_shape[0]*=dilation*dilation; temp_shape[2]*=dilation;}

if (threshold_division<0.0){ /* geometrical rule: division following the shortest direction */

    double shape_direction[2];

    eigenvalue=(temp_shape[0]+temp_shape[1])/2.;

    eigenvalue+=0.5*sqrt((temp_shape[0]-temp_shape[1])*(temp_shape[0]-temp_shape[1])+

                    4.*temp_shape[2]*temp_shape[2]);

    shape_direction[0]=temp_shape[2];

    shape_direction[1]=eigenvalue-temp_shape[0];

    norm=sqrt(shape_direction[0]*shape_direction[0]+shape_direction[1]*shape_direction[1]);

    shape_direction[0]/=norm;

    shape_direction[1]/=norm;

    division_direction[0]=shape_direction[1];

    division_direction[1]=-shape_direction[0];

}

/* computation of the new vertices as the intersection between the walls of the cell and the division direction */

int l=0,m=0,max=vertices_number_in_each_cell[dividing_cell+1],min=vertices_number_in_each_cell[dividing_cell],previous;

double test;

for (int i=min;i<max;i++){ /* loop over the vertices of the dividing cell */

    previous=i-1; if (i==min){previous=max-1;}

    test=-division_direction[1]*(vertices[2*cells_vertices[previous]]-barycenters[2*dividing_cell]);

    test+=division_direction[0]*(vertices[2*cells_vertices[previous]+1]-barycenters[2*dividing_cell+1]);

    test/=-division_direction[1]*(vertices[2*cells_vertices[previous]]-vertices[2*cells_vertices[i]])+

            division_direction[0]*(vertices[2*cells_vertices[previous]+1]-vertices[2*cells_vertices[i]+1]);

        /* test if the considered point of the wall is on the division direction */

    if ((test<1. && test>=0.) && (m==0 || m==1) && l<2){

        //if (test>0.9){test=0.9;} if (test<0.1){test=0.1;} // if new vertices should not be too closed of the older ones for stability reasons

        l++;

        if (l==1){

            new_vertex1[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex1[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-2; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-2; daughter_cells_vertices_number[1]+=1;

            wall1[0]=cells_vertices[i];

            wall1[1]=cells_vertices[previous];}

        if (l==2){

            new_vertex2[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex2[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-1; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-1; daughter_cells_vertices_number[1]+=1;

            wall2[0]=cells_vertices[i];

            wall2[1]=cells_vertices[previous];}

            m=2;

        }

    else if ((test<1. && test>=0.) && (m==0 || m==2) && l<2){

        //if (test>0.9){test=0.9;} if (test<0.1){test=0.1;}

        l++;

        if (l==1){

            new_vertex1[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex1[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-2; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-2; daughter_cells_vertices_number[1]+=1;

            wall1[0]=cells_vertices[i];

            wall1[1]=cells_vertices[previous];}

        if (l==2){

            new_vertex2[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];

            new_vertex2[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];

            cell1[daughter_cells_vertices_number[0]]=vertices_number-1; daughter_cells_vertices_number[0]+=1;

            cell2[daughter_cells_vertices_number[1]]=vertices_number-1; daughter_cells_vertices_number[1]+=1;

            wall2[0]=cells_vertices[i];

            wall2[1]=cells_vertices[previous];}

            m=1;

        }

    /* write the old vertices in the two daughter cells */

if (l==0 || l==2){cell1[daughter_cells_vertices_number[0]]=cells_vertices[i]; daughter_cells_vertices_number[0]+=1;}

if (l==1){cell2[daughter_cells_vertices_number[1]]=cells_vertices[i]; daughter_cells_vertices_number[1]+=1;}

}

/* if l<2, there has been a problem in the computation of the new vertices */

if (l!=2){printf("\n"); printf("division de %d ok,recherche des voisines\n",dividing_cell); exit(EXIT_FAILURE);}

/* the new vertices are saved in the list "vertices" */

vertices[2*(vertices_number-2)]=new_vertex1[0]; vertices[2*(vertices_number-2)+1]=new_vertex1[1];

vertices[2*(vertices_number-1)]=new_vertex2[0]; vertices[2*(vertices_number-1)+1]=new_vertex2[1];

/* find the cells which are adjacent to the dividing one and have to be modified to include the new vertices */

adjacent_cell1[0]=-1; adjacent_cell2[0]=-1; adjacent_cells_number[0]=0;

for (int i=0;i<cell_number-1;i++){ if (i==dividing_cell){goto nextcell;}

    min=vertices_number_in_each_cell[i]; max=vertices_number_in_each_cell[i+1]; previous=max-min-1;

    for (j=0;j<max-min;j++){

        if ((cells_vertices[min+j]==wall1[1] && cells_vertices[min+previous]==wall1[0])

            ||(cells_vertices[min+j]==wall2[1] && cells_vertices[min+previous]==wall2[0])

            ){

            if (adjacent_cells_number[0]==0){adjacent_cell1[0]=i;}

            if (adjacent_cells_number[0]==1){adjacent_cell2[0]=i;}

            adjacent_cells_number[0]+=1; goto nextcell;

            }

        previous=j;}

    nextcell:;

}

//printf("%d voisines trouvee(s): %d %d\n",adjacent_cells_number[0],adjacent_cell1[0],adjacent_cell2[0]);

}

#include "constants.h"

#include "growth_and_chemistry.h"

#include <gsl/gsl_errno.h>

#include <math.h>

#include <algorithm>

#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>

/* this function computes the equations of chemistry and growth */

int growth_and_chemistry(double t, const double x[], double dx[], void *parameters){

struct parameters_list2 {int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *shapes; double *lengths; double *cell_areas; double *barycenters;

                         int *cell_walls; int wall_number; int *boundary; int *wall_vertices; double *vertices; double *growth_rate; double *noise;};

/* diagonalize the deformation tensor in order to compute growth */

double delta_shapes[3*(((parameters_list2*)parameters)->cell_number)];

for (int i=0; i<(((parameters_list2*)parameters)->cell_number);i++){

    gsl_matrix *gsl_delta_shapes=gsl_matrix_calloc(2,2);

    gsl_matrix_set(gsl_delta_shapes,0,0,(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i]); gsl_matrix_set(gsl_delta_shapes,1,1,(((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1]);

    gsl_matrix_set(gsl_delta_shapes,1,0,(((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2]); gsl_matrix_set(gsl_delta_shapes,0,1,(((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2]);

    gsl_eigen_symmv_workspace *workspace=gsl_eigen_symmv_alloc(2); gsl_vector *eigenvalues=gsl_vector_alloc(2); gsl_matrix *eigenvectors=gsl_matrix_alloc(2,2); gsl_eigen_symmv(gsl_delta_shapes,eigenvalues,eigenvectors,workspace);

    if (gsl_vector_get(eigenvalues,0)<0.){ gsl_vector_set(eigenvalues,0,0.);}

    if (gsl_vector_get(eigenvalues,1)<0.){ gsl_vector_set(eigenvalues,1,0.);}

    delta_shapes[3*i]=gsl_vector_get(eigenvalues,0)*gsl_matrix_get(eigenvectors,0,0)*gsl_matrix_get(eigenvectors,0,0)+gsl_vector_get(eigenvalues,1)*gsl_matrix_get(eigenvectors,0,1)*gsl_matrix_get(eigenvectors,0,1);

    delta_shapes[3*i+1]=gsl_vector_get(eigenvalues,0)*gsl_matrix_get(eigenvectors,1,0)*gsl_matrix_get(eigenvectors,1,0)+gsl_vector_get(eigenvalues,1)*gsl_matrix_get(eigenvectors,1,1)*gsl_matrix_get(eigenvectors,1,1);

    delta_shapes[3*i+2]=gsl_vector_get(eigenvalues,0)*gsl_matrix_get(eigenvectors,0,0)*gsl_matrix_get(eigenvectors,1,0)+gsl_vector_get(eigenvalues,1)*gsl_matrix_get(eigenvectors,0,1)*gsl_matrix_get(eigenvectors,1,1);

    gsl_vector_free(eigenvalues); gsl_matrix_free(eigenvectors); gsl_matrix_free(gsl_delta_shapes); gsl_eigen_symmv_free(workspace);

}

    /* computation of the growth rate:

     the growth rate of periclinal walls is directly given by the list "growth_rate", and here we add the noise

     the growth of anticlinal walls is the average of the growth of the two neighbouring periclinal walls */

double cell_growth[(((parameters_list2*)parameters)->cell_number)];

double wall_growth[(((parameters_list2*)parameters)->wall_number)];

for (int i=0;i<(((parameters_list2*)parameters)->wall_number);i++){

    wall_growth[i]=0.;

}

for (int i=0;i<(((parameters_list2*)parameters)->cell_number);i++){

    //cell_growth[i]=1.;

    /* with noise */

    cell_growth[i]=(((parameters_list2*)parameters)->growth_rate)[i]*(((parameters_list2*)parameters)->noise)[i];

    for (int j=(((parameters_list2*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list2*)parameters)->vertices_number_in_each_cell)[i+1];j++){

        wall_growth[(((parameters_list2*)parameters)->cell_walls)[j]]+=0.5*cell_growth[i]*(((parameters_list2*)parameters)->noise)[i];

        if ((((parameters_list2*)parameters)->boundary)[j]==1){wall_growth[(((parameters_list2*)parameters)->cell_walls)[j]]*=2.;}

    }

}

    /* growth of the periclinal walls, proportionnal to the deformation, with possiblity to set negative growth, corresponding to compression, to zero

     this was not used in the end since it made no difference */

for (int i=0;i<(((parameters_list2*)parameters)->cell_number);i++){

    dx[3*i]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

               //std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i]);

               //(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i];

               delta_shapes[3*i];

    dx[3*i+1]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

               // std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1]);

               // (((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1];

              delta_shapes[3*i+1];

    dx[3*i+2]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

              //  std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2]);

             //   (((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2];

               delta_shapes[3*i+2];

}

    /* growth of the anticlinal walls */

for (int i=0;i<(((parameters_list2*)parameters)->wall_number);i++){

        dx[3*(((parameters_list2*)parameters)->cell_number)+i]=wall_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*

                std::max((((parameters_list2*)parameters)->lengths)[i]-x[3*(((parameters_list2*)parameters)->cell_number)+i],0.);

}

return GSL_SUCCESS;

}

di#include <iostream>

#include <string>

#include <stdlib.h>

#include <stdio.h>

#include <time.h>

#include "initialize.h"

#include "step.h"

#include "constants.h"

#include "cell_division.h"

#include "cell_division2.h"

#include <math.h>

#include "load_initial_state.h"

#include <fstream>

#include "compute_barycenters.h"

#include "compute_areas.h"

#include "compute_cell_areas.h"

#include "compute_shapes.h"

#include "compute_stress_periclinal.h"

#include "compute_stress_anticlinal.h"

#include "find_boundary_walls.h"

#include "minimize_BFGS.h"

#include <omp.h>

#include "compute_cell_areas.h"

#include "compute_shapes.h"

#include "compute_wall_vertices.h"

#include "compute_cell_walls.h"

#include <dirent.h>

#include "save_data.h"

#include <cmath>

#include "compute_lengths.h"

using namespace std;

int main(){

omp_set_dynamic(0);

srand(time(NULL));

double global_barycenter[2];

int init_vertices_number,init_cell_number,cumulated_vertices_number,init_wall_number;

// Pick the initial conditions, obtained from previous simulations

//init_cell_number=5; init_vertices_number=8; cumulated_vertices_number=20;

//const char* file_name="init/5_8_20.txt";

init_cell_number=101; init_vertices_number=200; cumulated_vertices_number=567;

const char* file_name="init/101_200_567.txt";

//init_cell_number=300; init_vertices_number=598; cumulated_vertices_number=1736;

//const char* file_name="init/300_598_1736.txt";

//init_cell_number=100; init_vertices_number=200; cumulated_vertices_number=561;

//const char* file_name="init/100_200_ablation.txt";

float * init_vertices; init_vertices=(float*)calloc(2*init_vertices_number,sizeof(float));

int * init_vertices_number_in_each_cell; init_vertices_number_in_each_cell=(int*)calloc(init_cell_number+1,sizeof(int));

int * init_cells_vertices; init_cells_vertices=(int*)calloc(cumulated_vertices_number,sizeof(int));

load_initial_state(file_name, init_cell_number, init_vertices_number, init_vertices, init_vertices_number_in_each_cell, init_cells_vertices);

// initial tissue loaded

/* center the tissue */

global_barycenter[1]=0; global_barycenter[2]=0;

//#pragma omp parallel shared(global_barycenter, init_vertices, init_vertices_number)

//{

//#pragma omp for

for (int i=0;i<init_vertices_number;i++){global_barycenter[1]+=init_vertices[2*i]/init_vertices_number;

                                global_barycenter[2]+=init_vertices[2*i+1]/init_vertices_number;}

//#pragma omp for

for (int i=0;i<init_vertices_number;i++){init_vertices[2*i]-=global_barycenter[1];

                                     init_vertices[2*i+1]-=global_barycenter[2];}

//}

// tissue centered

/* find the walls which are on the boundary of the tissue */

int* init_boundary; init_boundary=(int*)calloc(init_vertices_number_in_each_cell[init_cell_number],sizeof(int));

find_boundary_walls(init_cells_vertices, init_vertices_number_in_each_cell, init_cell_number, init_boundary, &init_wall_number);

int * init_wall_vertices; init_wall_vertices=(int*)calloc(2*init_wall_number,sizeof(int));

compute_wall_vertices( init_cells_vertices, init_vertices_number_in_each_cell, init_cell_number, init_wall_vertices);

// boundary identified

// define the parameters over witch there is a loop

// threshold=0 -> division following the highest tension

// threshold<0 -> division following the shortest axis of the cell

// initialize the files where the data will be saved. Yes that's a lot of them

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary_filter.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem_filter.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary_filter.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem_filter.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f_filter.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=7; m<=7; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f_filter.txt",1.-((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=5; m<=10; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/pressure/stress_%1.1f.txt",1.+((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=5; m<=10; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/pressure/shape_%1.1f.txt",1.+((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=5; m<=10; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/pressure/stress_%1.1f_filter.txt",1.+((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

//    for (int m=5; m<=10; m++){

//    FILE *f=NULL; char file_name[256];

//    sprintf(file_name,"results_v7/pressure/shape_%1.1f_filter.txt",1.+((double)m)*0.1);

//    f=fopen(file_name,"a+");

//    fprintf(f,"mother vertex daughter x y z\n");

//    fclose(f);

//    }

    for (int m=0; m<=0; m++){

    FILE *f=NULL; char file_name[256];

    sprintf(file_name,"results_v7/isotropic_2/stress.txt");

    f=fopen(file_name,"a+");

    fprintf(f,"mother vertex daughter x y z\n");

    fclose(f);

    }

    for (int m=0; m<=0; m++){

    FILE *f=NULL; char file_name[256];

    sprintf(file_name,"results_v7/isotropic_2/shape.txt");

    f=fopen(file_name,"a+");

    fprintf(f,"mother vertex daughter x y z\n");

    fclose(f);

    }

    for (int m=0; m<=0; m++){

    FILE *f=NULL; char file_name[256];

    sprintf(file_name,"results_v7/isotropic_2/stress_filter.txt");

    f=fopen(file_name,"a+");

    fprintf(f,"mother vertex daughter x y z\n");

    fclose(f);

    }

    for (int m=0; m<=0; m++){

    FILE *f=NULL; char file_name[256];

    sprintf(file_name,"results_v7/isotropic_2/shape_filter.txt");

    f=fopen(file_name,"a+");

    fprintf(f,"mother vertex daughter x y z\n");

    fclose(f);

    }

//

// parallel loop over the parameters

#pragma omp parallel for schedule(dynamic,1) collapse(3)

for (int meristem=0; meristem<=19; meristem++){ // several meristem, just the initial conditions are different

for (int anisotropy=0; anisotropy<=0; anisotropy++){ // amplitude of the stress anisotropy

    //for (int stress_source=1; stress_source<=1; stress_source++){ // origin of the stress: 0=growth, 1=curvature

        for (int rule=-1; rule<=0; rule++){

    ablated_cell=-1; // -1=no ablation, any other number n = the cell number n is killed

    iteration=0;

    threshold_division=(double)rule; // rule used for the cell division

            // negative number = cell divided along the short axis

            // 0 = cell divided along the maximal stress

            // positive number = probabilistic rule based on the mdeviatoric stress

            // define the stress anisotropy

//    if (stress_source==0){ // anisotropy comes from growth

//        growth_differential=0.999999-((double)anisotropy)*0.1;

//        first_shape_derived_stress=3.e-2;

//        second_shape_derived_stress=3.e-2;

//    }

//    if (stress_source==1){ // anisotropy comes from curvature

//        growth_differential=0.999999;

//        first_shape_derived_stress=3.e-2;

//        second_shape_derived_stress=(3.e-2)*(1.-((double)anisotropy)*0.1);

//    }

/* allocation of memory */

int vertices_number=init_vertices_number;

int cell_number=init_cell_number;

int wall_number=init_wall_number;

int i;

double * vertices; vertices=(double*)calloc(2*vertices_number,sizeof(double)); /* coordinates of the vertices */

for (i=0;i<2*vertices_number;i++){vertices[i]=(double)init_vertices[i];}

double r; double theta;

for (i=0;i<vertices_number;i++){r=0.3*(double)(rand())/(double)RAND_MAX; theta=2.*pi*(double)(rand())/(double)RAND_MAX;

    vertices[2*i]+=r*cos(theta); vertices[2*i+1]+=r*sin(theta);

}

int * vertices_number_in_each_cell; vertices_number_in_each_cell=(int*)calloc(cell_number+1,sizeof(int)); /* list of the number of vertices in each cell */

for (i=0;i<cell_number+1;i++){vertices_number_in_each_cell[i]=init_vertices_number_in_each_cell[i];}

int * cells_vertices; cells_vertices=(int*)calloc(vertices_number_in_each_cell[cell_number],sizeof(int)); /* list of the vertices of each cell, cell-per-cell */

for (i=0;i<vertices_number_in_each_cell[cell_number];i++){cells_vertices[i]=init_cells_vertices[i];}

int * wall_vertices; wall_vertices=(int*)calloc(2*wall_number,sizeof(int)); /* list of the vertices of all the walls, wall-per-wall */

for (i=0;i<2*wall_number;i++){wall_vertices[i]=init_wall_vertices[i];}

double * lengths; lengths=(double*)calloc(wall_number,sizeof(double)); /* lengths of the anticlinal walls */

double * target_lengths; target_lengths=(double*)calloc(wall_number,sizeof(double)); /* target lengths of the anticlinal walls */

int * cell_walls; cell_walls=(int*)calloc(vertices_number_in_each_cell[cell_number],sizeof(int)); /* list of the wall of all the cells, cell-per-cell

double * barycenters; barycenters=(double*)calloc(cell_number*2,sizeof(double)); /* positions of centers of cells */

double * shapes; shapes=(double*)calloc(cell_number*3,sizeof(double)); /* shapes of cells */

double * targets; targets=(double*)calloc(cell_number*3,sizeof(double)); /* targets of cells */

double * areas; areas=(double*)calloc(vertices_number_in_each_cell[cell_number],sizeof(double)); /* areas of all the triangles between two vertices and the barycenter */

double * cell_areas; cell_areas=(double*)calloc(cell_number, sizeof(double)); /* areas of cells */

double * forces; forces=(double*)calloc(2*vertices_number, sizeof(double)); /* forces on vertices */

double * stress_periclinal; stress_periclinal=(double*)calloc(3*cell_number, sizeof(double)); /* stress on the in-the-plane walls */

double * stress_anticlinal; stress_anticlinal=(double*)calloc(3*wall_number, sizeof(double)); /* stress on the out-the-plane walls */

int * boundary; boundary=(int*)calloc(vertices_number_in_each_cell[cell_number],sizeof(int)); /* indicates if a wall belongs to the boundary of the tissue */

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

double * growth_rate; growth_rate=(double*)calloc(cell_number,sizeof(double)); /* list of the growth rate of each cell*/

int * border; border=(int*)calloc(cell_number,sizeof(int)); /* indicates if a cell belongs to the boundary of the tissue  */

for (i=0; i<cell_number;i++){ growth_rate[i]=1.; border[i]=0;

    for (int j=vertices_number_in_each_cell[i]; j<vertices_number_in_each_cell[i+1];j++){

        if (boundary[j]==1.){ growth_rate[i]=growth_differential; border[i]=1; }

    }

}

double * relative_pressure; relative_pressure=(double*)calloc(cell_number,sizeof(double)); /* llist of the pressure in each cell */

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

/* from the vertices, initialize the program by calculating everything that has to be so: lengths, areas, stresses, ... */

initialize(vertices,vertices_number,cells_vertices,vertices_number_in_each_cell,cell_number, wall_vertices, wall_number, cell_walls,

           barycenters,shapes,targets,lengths,target_lengths,areas,cell_areas,forces,stress_periclinal, stress_anticlinal);

/* loop until the tissue reaches the limit size */

int check_division; int total_cell_number=1;

while (cell_number<max_cell_number /*&& iteration<step_number*/){check_division=0;

//while (iteration<step_number || ablated_cell==-1){check_division=0;

    compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters);

    compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas);

    compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas);

    compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes);

    /* one step of energy minimization and growth */

    step(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, wall_vertices, wall_number, cell_walls,

         barycenters, shapes, targets, lengths, target_lengths, areas, cell_areas, forces, stress_periclinal, stress_anticlinal,

         boundary, growth_rate, relative_pressure);

    /* compute the stress in the pericinal walls if it is necesssary for the divisions */

    if (threshold_division>=0.){

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

            compute_stress_periclinal(i, vertices, vertices_number_in_each_cell, cell_number, stress_periclinal, shapes, targets, cell_areas, wall_vertices, cell_walls, lengths, target_lengths);

        }

    }

    /* check if division is required, i.e. if a cell is too big */

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

        if (cell_areas[i]>threshold_area && i!=ablated_cell){check_division=1; int save_cell=1; /*if (stress_source==1){ */save_cell=1-border[i];//}

            /* test: increase the pressure in the cell just before division, see how it modifies the divisions */

                relative_pressure[i]=1.+(double)anisotropy*0.1;//}

            minimize_BFGS(vertices_number, cell_number, vertices, cells_vertices, vertices_number_in_each_cell, barycenters, shapes, targets, areas, cell_areas, wall_vertices, wall_number, lengths, target_lengths, boundary,relative_pressure);

            compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas);

            compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas); //printf("step2\n");int bug[1]={0};

            compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters);

            compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes);

            compute_lengths(vertices, wall_vertices, wall_number, lengths);

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

                compute_stress_periclinal(k, vertices, vertices_number_in_each_cell, cell_number, stress_periclinal, shapes, targets, cell_areas, wall_vertices, cell_walls, lengths, target_lengths);

            }

            relative_pressure[i]=1.;

            /* end of the test, pressure back to normal */

    /* temporal variables to identify the neighbours of the dividing cell  */

    int adjacent_cells_number[1]; int adjacent_cell1[1]; int adjacent_cell2[1];

    int daughter_cells_vertices_number[2]={0,0}; int cell1[10],cell2[10],wall1[2],wall2[2];

    /* division, with reallocation of the memory */

    cell_number+=1; vertices_number+=2; wall_number+=3;

    vertices=(double*)realloc(vertices, sizeof(double)*2*vertices_number);

    forces=(double*)realloc(forces, sizeof(double)*2*vertices_number);

    vertices_number_in_each_cell=(int*)realloc(vertices_number_in_each_cell, sizeof(int)*(cell_number+1));

    barycenters=(double*)realloc(barycenters, sizeof(double)*2*cell_number);

    shapes=(double*)realloc(shapes, sizeof(double)*3*cell_number);

    targets=(double*)realloc(targets, sizeof(double)*3*cell_number);

    stress_periclinal=(double*)realloc(stress_periclinal, sizeof(double)*3*cell_number);

    stress_anticlinal=(double*)realloc(stress_anticlinal, sizeof(double)*3*wall_number);

    cell_areas=(double*)realloc(cell_areas, sizeof(double)*cell_number);

    growth_rate=(double*)realloc(growth_rate,sizeof(double)*cell_number); growth_rate[cell_number-1]=growth_rate[i];

    border=(int*)realloc(border,sizeof(int)*cell_number); border[cell_number-1]=0;

    relative_pressure=(double*)realloc(relative_pressure,sizeof(double)*cell_number); relative_pressure[cell_number-1]=1.;

            /* find the direction of division and the position of the new vertices */

    cell_division(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number,

          barycenters, i, boundary, stress_periclinal, shapes,

          adjacent_cells_number, adjacent_cell1, adjacent_cell2, daughter_cells_vertices_number, cell1, cell2, wall1, wall2);

    cells_vertices=(int*)realloc(cells_vertices, sizeof(int)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0]));

    wall_vertices=(int*)realloc(wall_vertices, sizeof(int)*2*wall_number);

    lengths=(double*)realloc(lengths, sizeof(double)*(wall_number));

    target_lengths=(double*)realloc(target_lengths, sizeof(double)*(wall_number));

    cell_walls=(int*)realloc(cell_walls, sizeof(int)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0]));

    areas=(double*)realloc(areas, sizeof(double)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0]));

    boundary=(int*)realloc(boundary, sizeof(int)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0]));

            /* actually performs the division based on the division identified above */

    cell_division2(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, wall_vertices, wall_number,

        barycenters, areas, cell_areas, lengths, target_lengths, shapes, targets, i, boundary, border,

        adjacent_cells_number[0], adjacent_cell1[0], adjacent_cell2[0], daughter_cells_vertices_number, cell1, cell2, wall1, wall2);

            /* update what is modified by the division */

    compute_cell_walls(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, cell_walls, wall_vertices, wall_number);

    compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas);

    compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas);

    compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters);

    compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes);

    compute_lengths(vertices, wall_vertices, wall_number, lengths);

            /* minimize the energy after division */

    minimize_BFGS(vertices_number, cell_number, vertices, cells_vertices, vertices_number_in_each_cell, barycenters, shapes, targets, areas, cell_areas, wall_vertices, wall_number, lengths, target_lengths, boundary,relative_pressure);

    compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas);

    compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas);

    compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters);

    compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes);

    compute_lengths(vertices, wall_vertices, wall_number, lengths);

            /* save the history of the divisions, used as by Marion in her statistical analysis */

    #pragma omp critical

    {

    if (save_cell==1){

        FILE *f=NULL; char file_name[256];

        /*if (stress_source==1){

            if (rule==-1){

                sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f.txt",1.-((double)anisotropy)*0.1);

            }

            if (rule==0){

                sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f.txt",1.-((double)anisotropy)*0.1);

            }

        }

        if (stress_source==0){

           if (growth_rate[i]==1.){

                if (rule==-1){

                    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);

                }

                if (rule==0){

                    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);

                }

            }

            if (growth_rate[i]!=1.){

                if (rule==-1){

                    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);

                }

                if (rule==0){

                    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);

                }

            }

        }*/

//            if (rule==-1){

//                sprintf(file_name,"results_v7/pressure/shape_%1.1f.txt",1.+((double)anisotropy)*0.1);

//            }

//            if (rule==0){

//                sprintf(file_name,"results_v7/pressure/stress_%1.1f.txt",1.+((double)anisotropy)*0.1);

//            }

            if (rule==-1){

                sprintf(file_name,"results_v7/isotropic_2/shape.txt");

            }

            if (rule==0){

                sprintf(file_name,"results_v7/isotropic_2/stress.txt");

            }

        f=fopen(file_name,"a");

        int write=0,count=0,count2=0;

        while (count<(daughter_cells_vertices_number[1]-1)){

            if (write){

                fprintf(f,"%d %d %d %f %f 0.0\n",1000*meristem+total_cell_number/*i+1*/,cell2[count2%daughter_cells_vertices_number[1]]+1,/*total_cell_number+1*/0/*cell_number-1*/,

                        vertices[2*cell2[count2%daughter_cells_vertices_number[1]]],

                        vertices[2*cell2[count2%daughter_cells_vertices_number[1]]+1]);

                count++;

            }

            write+=((cell2[count2%daughter_cells_vertices_number[1]])==(vertices_number-2));

            count2++;

        }

        write=0; count=0; count2=0;

        while (count<(daughter_cells_vertices_number[0]-1)){

            if (write){

                fprintf(f,"%d %d %d %f %f 0.0\n",1000*meristem+total_cell_number/*i+1*/,cell1[count2%daughter_cells_vertices_number[0]]+1,/*total_cell_number+2*/1/*i*/,

                        vertices[2*cell1[count2%daughter_cells_vertices_number[0]]],

                        vertices[2*cell1[count2%daughter_cells_vertices_number[0]]+1]);

                count++;

            }

            write+=((cell1[count2%daughter_cells_vertices_number[0]])==(vertices_number-1));

            count2++;

        }

        fclose(f);

        }

    }

            /* division saved */

            /* same as above, but filter out the asymmetric divisions */

    #pragma omp critical

    {

    double area_ratio=cell_areas[i]/(cell_areas[cell_number-1]+cell_areas[i]); if (area_ratio>0.5){ area_ratio=1.-area_ratio;}

    double aspect_ratio=(shapes[3*i]+shapes[3*i+1])/2.;

    aspect_ratio-=0.5*sqrt((shapes[3*i]-shapes[3*i+1])*(shapes[3*i]-shapes[3*i+1])+4.*shapes[3*i+2]*shapes[3*i+2]);

    aspect_ratio=aspect_ratio/(shapes[3*i]+shapes[3*i+1]-aspect_ratio); aspect_ratio=sqrt(aspect_ratio);

    if (save_cell==1 && area_ratio>0.4 && aspect_ratio>0.5){

        FILE *f=NULL; char file_name[256];

        /*if (stress_source==1){

            if (rule==-1){

                sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f_filter.txt",1.-((double)anisotropy)*0.1);

            }

            if (rule==0){

                sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f_filter.txt",1.-((double)anisotropy)*0.1);

            }

        }

        if (stress_source==0){

           if (growth_rate[i]==1.){

                if (rule==-1){

                    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem_filter.txt",1.-((double)anisotropy)*0.1);

                }

                if (rule==0){

                    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem_filter.txt",1.-((double)anisotropy)*0.1);

                }

            }

            if (growth_rate[i]!=1.){

                if (rule==-1){

                    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary_filter.txt",1.-((double)anisotropy)*0.1);

                }

                if (rule==0){

                    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary_filter.txt",1.-((double)anisotropy)*0.1);

                }

            }

        }*/

//            if (rule==-1){

//                sprintf(file_name,"results_v7/pressure/shape_%1.1f_filter.txt",1.+((double)anisotropy)*0.1);

//            }

//            if (rule==0){

//                sprintf(file_name,"results_v7/pressure/stress_%1.1f_filter.txt",1.+((double)anisotropy)*0.1);

//            }

            if (rule==-1){

                sprintf(file_name,"results_v7/isotropic_2/shape_filter.txt");

            }

            if (rule==0){

                sprintf(file_name,"results_v7/isotropic_2/stress_filter.txt");

            }

        f=fopen(file_name,"a");

        int write=0,count=0,count2=0;

        while (count<(daughter_cells_vertices_number[1]-1)){

            if (write){

                fprintf(f,"%d %d %d %f %f 0.0\n",1000*meristem+total_cell_number/*i+1*/,cell2[count2%daughter_cells_vertices_number[1]]+1,/*total_cell_number+1*/0/*cell_number-1*/,

                        vertices[2*cell2[count2%daughter_cells_vertices_number[1]]],

                        vertices[2*cell2[count2%daughter_cells_vertices_number[1]]+1]);

                count++;

            }

            write+=((cell2[count2%daughter_cells_vertices_number[1]])==(vertices_number-2));

            count2++;

        }

        write=0; count=0; count2=0;

        while (count<(daughter_cells_vertices_number[0]-1)){

            if (write){

                fprintf(f,"%d %d %d %f %f 0.0\n",1000*meristem+total_cell_number/*i+1*/,cell1[count2%daughter_cells_vertices_number[0]]+1,/*total_cell_number+2*/1/*i*/,

                        vertices[2*cell1[count2%daughter_cells_vertices_number[0]]],

                        vertices[2*cell1[count2%daughter_cells_vertices_number[0]]+1]);

                count++;

            }

            write+=((cell1[count2%daughter_cells_vertices_number[0]])==(vertices_number-1));

            count2++;

        }

        fclose(f);

        }

    }

    #pragma omp critical

    {

    if (save_cell==1){

        FILE *f=NULL; char file_name[256];

        /*if (stress_source==1){

            if (rule==-1){

                sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f.txt",1.-((double)anisotropy)*0.1);

            }

            if (rule==0){

                sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f.txt",1.-((double)anisotropy)*0.1);

            }

        }

        if (stress_source==0){

           if (growth_rate[i]==1.){

                if (rule==-1){

                    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);

                }

                if (rule==0){

                    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);

                }

            }

            if (growth_rate[i]!=1.){

                if (rule==-1){

                    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);

                }

                if (rule==0){

                    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);

                }

            }

        }*/

//            if (rule==-1){

//                sprintf(file_name,"results_v7/pressure/shape_%1.1f.txt",1.+((double)anisotropy)*0.1);

//            }

//            if (rule==0){

//                sprintf(file_name,"results_v7/pressure/stress_%1.1f.txt",1.+((double)anisotropy)*0.1);

//            }

            if (rule==-1){

                sprintf(file_name,"results_v7/isotropic_2/shape.txt");

            }

            if (rule==0){

                sprintf(file_name,"results_v7/isotropic_2/stress.txt");

            }

        f=fopen(file_name,"a");

        int write=0,count=0,count2=0;

        while (count<(daughter_cells_vertices_number[1]-1)){

            if (write){

                fprintf(f,"%d %d %d %f %f 0.0\n",1000*meristem+total_cell_number/*i+1*/,cell2[count2%daughter_cells_vertices_number[1]]+1,/*total_cell_number+1*/0/*cell_number-1*/,

                        vertices[2*cell2[count2%daughter_cells_vertices_number[1]]],

                        vertices[2*cell2[count2%daughter_cells_vertices_number[1]]+1]);

                count++;

            }

            write+=((cell2[count2%daughter_cells_vertices_number[1]])==(vertices_number-2));

            count2++;

        }

        write=0; count=0; count2=0;

        while (count<(daughter_cells_vertices_number[0]-1)){

            if (write){

                fprintf(f,"%d %d %d %f %f 0.0\n",1000*meristem+total_cell_number/*i+1*/,cell1[count2%daughter_cells_vertices_number[0]]+1,/*total_cell_number+2*/1/*i*/,

                        vertices[2*cell1[count2%daughter_cells_vertices_number[0]]],

                        vertices[2*cell1[count2%daughter_cells_vertices_number[0]]+1]);

                count++;

            }

            write+=((cell1[count2%daughter_cells_vertices_number[0]])==(vertices_number-1));

            count2++;

        }

        fclose(f);

        }

    }

    #pragma omp critical

    {

    double area_ratio=cell_areas[i]/(cell_areas[cell_number-1]+cell_areas[i]); if (area_ratio>0.5){ area_ratio=1.-area_ratio;}

    double aspect_ratio=(shapes[3*i]+shapes[3*i+1])/2.;

    aspect_ratio-=0.5*sqrt((shapes[3*i]-shapes[3*i+1])*(shapes[3*i]-shapes[3*i+1])+4.*shapes[3*i+2]*shapes[3*i+2]);

    aspect_ratio=aspect_ratio/(shapes[3*i]+shapes[3*i+1]-aspect_ratio); aspect_ratio=sqrt(aspect_ratio);

    if (save_cell==1){

        FILE *f=NULL; char file_name[256];

            if (rule==-1){

                sprintf(file_name,"results_v7/isotropic_2/shape_data.txt");

            }

            if (rule==0){

                sprintf(file_name,"results_v7/isotropic_2/stress_data.txt");

            }

        f=fopen(file_name,"a");

        fprintf(f,"%d %f %f %d %f %f %f\n",1000*meristem+total_cell_number,cell_areas[i],aspect_ratio,

                vertices_number_in_each_cell[i+1]-vertices_number_in_each_cell[i],stress_periclinal[3*i],stress_periclinal[3*i+1],stress_periclinal[3*i+2]);

        fclose(f);

        }

    }

    #pragma omp critical

    {

    double area_ratio=cell_areas[i]/(cell_areas[cell_number-1]+cell_areas[i]); if (area_ratio>0.5){ area_ratio=1.-area_ratio;}

    double aspect_ratio=(shapes[3*i]+shapes[3*i+1])/2.;

    aspect_ratio-=0.5*sqrt((shapes[3*i]-shapes[3*i+1])*(shapes[3*i]-shapes[3*i+1])+4.*shapes[3*i+2]*shapes[3*i+2]);

    aspect_ratio=aspect_ratio/(shapes[3*i]+shapes[3*i+1]-aspect_ratio); aspect_ratio=sqrt(aspect_ratio);