Revision 3

src/compute_energy_and_forces.cpp (revision 3)
1
#include "compute_energy_and_forces.h"
2
#include "constants.h"
3
#include "compute_areas.h"
4
#include "compute_cell_areas.h"
5
#include "compute_shapes.h"
6
#include "compute_barycenters.h"
7
#include "compute_lengths.h"
8
#include <math.h>
9
#include <stdio.h>
10
#include <nlopt.h>
11
#include <omp.h>
12

  
13
/* this function computes the mechanical energy and its derivatives (basically only math equations) */
14

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

  
17
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;
18
                        int *wall_vertices; int wall_number; double *lengths; double *target_lengths; int *boundary; double *tissue_target; int anisotropy; double *relative_pressure;};
19

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

  
22
compute_lengths((double*)vertices, ((parameters_list*)parameters)->wall_vertices, ((parameters_list*)parameters)->wall_number, ((parameters_list*)parameters)->lengths);
23
compute_barycenters((double*)vertices, ((parameters_list*)parameters)->vertices_number, ((parameters_list*)parameters)->cells_vertices,
24
                     ((parameters_list*)parameters)->vertices_number_in_each_cell, ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->barycenters);
25
compute_areas((double*)vertices, ((parameters_list*)parameters)->cells_vertices, ((parameters_list*)parameters)->vertices_number_in_each_cell,
26
              ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->areas);
27
compute_cell_areas((double*)vertices, ((parameters_list*)parameters)->cells_vertices, ((parameters_list*)parameters)->vertices_number_in_each_cell,
28
              ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->areas, ((parameters_list*)parameters)->cell_areas);
29
compute_shapes((double*)vertices, ((parameters_list*)parameters)->cells_vertices, ((parameters_list*)parameters)->vertices_number_in_each_cell,
30
               ((parameters_list*)parameters)->cell_number, ((parameters_list*)parameters)->barycenters, ((parameters_list*)parameters)->shapes);
31

  
32

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

  
35
//for (int i=0;i<((parameters_list*)parameters)->vertices_number;i++){
36
//    if (vertices[2*i+1]>7.){ forces[2*i+1]+=0.01*sinh(vertices[2*i+1]-7.);}
37
//    if (vertices[2*i+1]<-7.){ forces[2*i+1]+=0.01*sinh(vertices[2*i+1]+7.);}
38
//    if (vertices[2*i+1]>7.){ energy+=0.01*(cosh(vertices[2*i+1]-7.)-1.);}
39
//    if (vertices[2*i+1]<-7.){ energy+=0.01*(cosh(vertices[2*i+1]+7.)-1.);}
40
//}
41

  
42
double ablation;
43
for (int i=0;i<((parameters_list*)parameters)->cell_number;i++){
44
    ablation=(double)(i!=ablated_cell);
45
    if (forces){
46
    for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){
47
        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;}
48
        /* derivatives of the area in the pressure energy */
49
        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]);
50
        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]]);
51
        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]);
52
        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]]);
53
    }
54
    }
55
    /* pressure energy */
56
    energy+=-pressure*(((parameters_list*)parameters)->relative_pressure)[i]*ablation*(((parameters_list*)parameters)->cell_areas)[i];
57
}
58

  
59
double boundary;
60
for (int i=0;i<((parameters_list*)parameters)->cell_number;i++){
61
//    for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){
62
//        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;}
63
//        boundary=(double)((((parameters_list*)parameters)->boundary)[j]);
64
//        if (forces){
65
//            forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]-=boundary*0.5*pressure*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]);
66
//            forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-=boundary*0.5*pressure*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]]);
67
//            forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]]-=boundary*0.5*pressure*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]);
68
//            forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]-=boundary*0.5*pressure*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]);
69
//        }
70
//
71
//        energy-=boundary*0.5*pressure*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-
72
//                                                vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]);
73
//    }
74
}
75

  
76
double trace;
77
for (int i=0;i<((parameters_list*)parameters)->cell_number;i++){
78
    trace=(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->targets)[3*i+1];
79
    N=(double)((((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1])-(double)((((parameters_list*)parameters)->vertices_number_in_each_cell)[i]);
80
    ablation=(double)(i!=ablated_cell);
81
    if (forces){
82
    for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){
83
        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;}
84
        /* derivatives of the shape matrix */
85
        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])
86
                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]-(((parameters_list*)parameters)->barycenters)[2*i])*
87
                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);
88
        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])
89
                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-(((parameters_list*)parameters)->barycenters)[2*i+1])*
90
                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);
91
        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])
92
                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-(((parameters_list*)parameters)->barycenters)[2*i+1])*
93
                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);
94
        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])
95
                                                                        *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]-(((parameters_list*)parameters)->barycenters)[2*i])*
96
                                                                        (((parameters_list*)parameters)->cell_areas)[i]/(trace*trace*N);
97
        /* derivatives of the area */
98
        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]+=ablation*shear_modulus*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1])*
99
                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+
100
                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)
101
                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);
102
        forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]+=ablation*shear_modulus*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]])*
103
                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+
104
                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)
105
                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);
106
        forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]]+=ablation*shear_modulus*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1])*
107
                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+
108
                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)
109
                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);
110
        forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]+=ablation*shear_modulus*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]])*
111
                                    (pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)+
112
                                     pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)
113
                                    +2.*pow((((parameters_list*)parameters)->shapes)[3*i+2]-(((parameters_list*)parameters)->targets)[3*i+2],2.))/(2.*trace*trace);
114
    }
115
    }
116
    /* shape energy, contribution from the shear modulus */
117
    energy+=ablation*shear_modulus*(((parameters_list*)parameters)->cell_areas)[i]*(pow((((parameters_list*)parameters)->shapes)[3*i]-(((parameters_list*)parameters)->targets)[3*i],2.)
118
                        +pow((((parameters_list*)parameters)->shapes)[3*i+1]-(((parameters_list*)parameters)->targets)[3*i+1],2.)
119
                        +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];
120
}
121

  
122
    /* contribution of the first lame parameter, this term was neglected in the final version of the paper as it has no impact */
123
for (int i=0;i<((parameters_list*)parameters)->cell_number;i++){
124
//    min=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i]; max=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];
125
//    trace=(((parameters_list*)parameters)->targets)[3*i]+(((parameters_list*)parameters)->targets)[3*i+1];
126
//    N=(double)(max)-(double)(min);
127
//    if (forces){
128
//    for (j=0;j<max-min;j++){ previous=j-1; if (j==0){ previous=max-min-1;}
129
//        /* derivatives of the shape matrix */
130
//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]]+=lame_first_periclinal*4.0*(((parameters_list*)parameters)->cell_areas)[i]*
131
//                                                    ((((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])
132
//                                                    *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]]-(((parameters_list*)parameters)->barycenters)[2*i])
133
//                                                    /(trace*trace*N);
134
//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1]+=lame_first_periclinal*4.0*(((parameters_list*)parameters)->cell_areas)[i]*
135
//                                                    ((((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])
136
//                                                    *(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1]-(((parameters_list*)parameters)->barycenters)[2*i+1])
137
//                                                    /(trace*trace*N);
138
//        /* derivatives of the area */
139
//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]]+=lame_first_periclinal*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]+1])*
140
//                                    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.)
141
//                                    /(2.*trace*trace);
142
//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1]+=lame_first_periclinal*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]])*
143
//                                    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.)
144
//                                    /(2.*trace*trace);
145
//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]]+=lame_first_periclinal*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]+1])*
146
//                                    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.)
147
//                                    /(2.*trace*trace);
148
//        forces[2*(((parameters_list*)parameters)->cells_vertices)[min+previous]+1]+=lame_first_periclinal*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[min+j]])*
149
//                                    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.)
150
//                                    /(2.*trace*trace);
151
//    }
152
//    }
153
//    /* area energy */
154
//    energy+=lame_first_periclinal*(((parameters_list*)parameters)->cell_areas)[i]*
155
//                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.)
156
//                /(trace*trace);
157
}
158

  
159
    /* energy of the anticlinal walls */
160
for (int i=0;i<((parameters_list*)parameters)->wall_number;i++){
161
//    /* derivatives of the lengths of the walls */
162
//    if (forces){
163
//    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])*
164
//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]])/
165
//                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->target_lengths)[i]);
166
//    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])*
167
//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]])/
168
//                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->target_lengths)[i]);
169
//    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])*
170
//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]+1]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]+1])/
171
//                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->target_lengths)[i]);
172
//    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])*
173
//                                                                        (vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i+1]+1]-vertices[2*(((parameters_list*)parameters)->wall_vertices)[2*i]+1])/
174
//                                                                        ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->target_lengths)[i]);
175
//    }
176
//    /* walls energy */
177
//    /* the 2. factor comes from the fact that each wall is counted twice in the sum over the cells */
178
//    energy+=lame_first_anticlinal*pow((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->target_lengths)[i],2.)/(((parameters_list*)parameters)->target_lengths)[i];
179
}
180

  
181
    /* contribution from the anisotropic stress that mimics curvature */
182
if (((parameters_list*)parameters)->anisotropy==1){
183
    double tissue_area=0.; int boundary_vertices_number=0;
184
    double tissue_shape[2]={0.,0.}; double tissue_barycenter[2]={0.,0.};
185
    for (int i=0;i<((parameters_list*)parameters)->cell_number;i++){
186
        for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){
187
            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;}
188
            boundary_vertices_number+=(((parameters_list*)parameters)->boundary)[j];
189
            boundary=(double)((((parameters_list*)parameters)->boundary)[j]);
190
            tissue_area+=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-
191
                                   vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]]*vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1])/2.;
192
            tissue_barycenter[0]+=boundary*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]];
193
            tissue_barycenter[1]+=boundary*vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1];
194
        }
195
    }
196
    tissue_barycenter[0]/=boundary_vertices_number;
197
    tissue_barycenter[1]/=boundary_vertices_number;
198
    for (int i=0;i<(((parameters_list*)parameters)->vertices_number_in_each_cell)[((parameters_list*)parameters)->cell_number];i++){
199
        boundary=(double)((((parameters_list*)parameters)->boundary)[i]);
200
        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]);
201
        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]);
202
    }
203
    tissue_shape[0]/=boundary_vertices_number;
204
    tissue_shape[1]/=boundary_vertices_number;
205
    trace=(((parameters_list*)parameters)->tissue_target)[0]+(((parameters_list*)parameters)->tissue_target)[1];
206
    if (forces){
207
        for (int i=0;i<((parameters_list*)parameters)->cell_number;i++){
208
            for (j=(((parameters_list*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list*)parameters)->vertices_number_in_each_cell)[i+1];j++){
209
                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;}
210
                boundary=(double)((((parameters_list*)parameters)->boundary)[j]);
211
                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])/
212
                                                                               (trace*boundary_vertices_number);
213
                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])/
214
                                                                                 (trace*boundary_vertices_number);
215

  
216
                forces[2*(((parameters_list*)parameters)->cells_vertices)[j]]-=boundary*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1])*
217
                                                                               (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+
218
                                                                               second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);
219
                forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]]-=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]+1])*
220
                                                                                      (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+
221
                                                                                       second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);
222
                forces[2*(((parameters_list*)parameters)->cells_vertices)[j]+1]-=boundary*(vertices[2*(((parameters_list*)parameters)->cells_vertices)[previous]])*
223
                                                                                 (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+
224
                                                                                 second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);
225
                forces[2*(((parameters_list*)parameters)->cells_vertices)[previous]+1]-=boundary*(-vertices[2*(((parameters_list*)parameters)->cells_vertices)[j]])*
226
                                                                                        (first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+
227
                                                                                        second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/(2.*trace);
228
            }
229
        }
230
    }
231
    energy-=tissue_area*(first_shape_derived_stress*(tissue_shape[0]-(((parameters_list*)parameters)->tissue_target)[0])+
232
                         second_shape_derived_stress*(tissue_shape[1]-(((parameters_list*)parameters)->tissue_target)[1]))/trace;
233
}
234

  
235
return(energy);
236

  
237
}
src/step.cpp (revision 3)
1
#include "compute_cell_areas.h"
2
#include "compute_areas.h"
3
#include "compute_lengths.h"
4
#include "compute_barycenters.h"
5
#include "step.h"
6
#include <stdio.h>
7
#include <math.h>
8
#include "minimize_BFGS.h"
9
#include "minimize_steepest_descent.h"
10
#include "growth_and_chemistry.h"
11
#include "compute_stress_periclinal.h"
12
#include "compute_stress_anticlinal.h"
13
#include "compute_cell_walls.h"
14
#include <gsl/gsl_odeiv2.h>
15
#include <gsl/gsl_errno.h>
16
#include "compute_energy_and_forces.h"
17
#include "compute_center_of_mass.h"
18
#include "save_data.h"
19
#include "compute_shapes.h"
20
#include <cstdlib>
21
#include "constants.h"
22

  
23

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

  
26
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,
27
          double *barycenters, double *shapes, double *targets, double *lengths, double *target_lengths, double *areas, double *cell_areas, double *forces, double *stress_periclinal, double *stress_anticlinal,
28
          int *boundary, double *growth_rate, double *relative_pressure){
29

  
30
/* minimization of energy: move the vertices in order to make the energy as small as possible, i.e. in order to make cells
31
   as close as possible of their targets */
32
//printf("step1\n");
33
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);
34

  
35
compute_cell_walls(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, cell_walls, wall_vertices, wall_number);
36
compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas);
37
compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas); //printf("step2\n");int bug[1]={0};
38
compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters);
39
compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes);
40
compute_lengths(vertices, wall_vertices, wall_number, lengths);
41

  
42
for (int i=0; i<cell_number;i++){
43
    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);
44
}
45
for (int i=0; i<wall_number;i++){
46
    compute_stress_anticlinal(i, vertices, vertices_number_in_each_cell, cell_number, stress_anticlinal, wall_vertices, cell_walls, lengths, target_lengths);
47
}
48

  
49
//for (int i=0; i<cell_number;i++){ stress_periclinal[3*i]=0; stress_periclinal[3*i+1]=0; stress_periclinal[3*i+2]=0;
50
//    for (int j=vertices_number_in_each_cell[i]; j<vertices_number_in_each_cell[i+1];j++){
51
//        stress_periclinal[3*i]+=stress_anticlinal[3*cell_walls[j]]; stress_periclinal[3*i+1]+=stress_anticlinal[3*cell_walls[j]+1]; stress_periclinal[3*i+2]+=stress_anticlinal[3*cell_walls[j]+2];
52
//    }
53
//}
54

  
55
        /* save the data */
56
//if ((iteration%100)==0 /*|| iteration==51*/){
57
//    FILE *f=NULL; char file_name[256];
58
//    sprintf(file_name,"results/boundary/results_%i/%i_%i_%i.txt",m,(int)threshold_division,((m/4)%4),iteration); f=fopen(file_name,"w+");
59
//    //sprintf(file_name,"results/ablation/%i_%i_%i_%i.txt",(int)threshold_division,(m/4)%300,m%4,iteration); f=fopen(file_name,"w+");
60
//    //sprintf(file_name,"results/isotropic/%i_%i_%i.txt",(int)threshold_division,m/2,iteration); f=fopen(file_name,"w+");
61
//    //sprintf(file_name,"results/results_%i/%i.txt",m,iteration); f=fopen(file_name,"w+");
62
//    save_data(f, vertices_number, vertices, cell_number, vertices_number_in_each_cell, cells_vertices, wall_number, barycenters, stress_periclinal, stress_anticlinal,  boundary, cell_areas, lengths,shapes, targets, growth_rate, cell_walls);
63
//    fclose(f);
64
//}
65

  
66
/* center the tissue */
67
//double center[2]={0,0};
68
//for (int i=0;i<cell_number;i++){center[0]+=barycenters[2*i]; center[1]+=barycenters[2*i+1];}
69
//center[0]/=cell_number; center[1]/=cell_number;
70
//for (int i=0;i<cell_number;i++){barycenters[2*i]-=center[0]; barycenters[2*i+1]-=center[1];}
71
//for (int i=0;i<vertices_number;i++){vertices[2*i]-=center[0]; vertices[2*i+1]-=center[1];}
72

  
73
    /* resolution of growth */
74
double noise[cell_number];
75
for (int i=0; i<cell_number;i++){
76
    noise[i]=1.+noise_amplitude*(2.*(double)(rand())/(double)RAND_MAX-1.);
77
}
78

  
79
struct parameters_list2 {int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *shapes; double *lengths; double *cell_areas; double *barycenters;
80
                         int *cell_walls; int wall_number; int *boundary; int *wall_vertices; double *vertices; double *growth_rate; double *noise;};
81

  
82
struct parameters_list2 parameters2;
83
parameters2.cell_number=cell_number;
84
parameters2.cells_vertices=cells_vertices;
85
parameters2.vertices_number_in_each_cell=vertices_number_in_each_cell;
86
parameters2.shapes=shapes;
87
parameters2.lengths=lengths;
88
parameters2.cell_areas=cell_areas;
89
parameters2.barycenters=barycenters;
90
parameters2.cell_walls=cell_walls;
91
parameters2.wall_number=wall_number;
92
parameters2.boundary=boundary;
93
parameters2.wall_vertices=wall_vertices;
94
parameters2.vertices=vertices;
95
parameters2.growth_rate=growth_rate;
96
parameters2.noise=noise;
97
void *pparameters2=&parameters2;
98
int system_size=3*cell_number+wall_number;
99

  
100
double x[system_size];
101
for (int i=0;i<3*cell_number;i++){ x[i]=targets[i];}
102
for (int i=0;i<wall_number;i++){ x[3*cell_number+i]=target_lengths[i];}
103

  
104
gsl_odeiv2_system sys={growth_and_chemistry, NULL, system_size, pparameters2};
105
gsl_odeiv2_driver *d=gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk2, 1.e-3, 1.e-3, 0.);
106
int s; double t=0.0;
107
s=gsl_odeiv2_driver_apply(d, &t, 1.e1, x);
108
//printf("résolution croissance/chimie: %d\n",s);
109
//printf("test\n");
110
for (int i=0;i<3*cell_number;i++){ targets[i]=x[i];}
111
for (int i=0;i<wall_number;i++){ target_lengths[i]=x[3*cell_number+i];}
112
gsl_odeiv2_driver_free(d);
113

  
114
}
115

  
116

  
117

  
118

  
119

  
src/constants.cpp (revision 3)
1
#include <math.h>
2

  
3
int step_number=1;//51;
4
int max_cell_number=200;//400; //1000000
5
double tau1=10.;
6

  
7
double pi=3.1415926;
8

  
9
double initial_growth=1.;
10
double threshold_area=3.;
11
double threshold_division=0.;
12

  
13
double shear_modulus=1.;
14
double lame_first_periclinal=0.;//0.1;
15
double lame_first_anticlinal=0.0;//0.02;
16
double pressure=1.e-2;//6.e-2;//12.e-2;
17

  
18
double first_shape_derived_stress=3.e-2;// 3.e-2
19
double second_shape_derived_stress=3.e-2;
20

  
21
double noise_amplitude=0.1;
22

  
23
int ablated_cell=-1;
24
int iteration=0;
25
int iteration_ablation=1000000; //50
26

  
27
double growth_differential=1.; //>=0
28

  
29
#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)
src/compute_stress_anticlinal.cpp (revision 3)
1
#include "compute_stress_anticlinal.h"
2
#include "constants.h"
3
#include <stdio.h>
4

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

  
7
void compute_stress_anticlinal(int wall_index, double *vertices, int *vertices_number_in_each_cell, int cell_number, double *stress_anticlinal,
8
                                int *wall_vertices, int *cell_walls, double *lengths, double *target_lengths){
9

  
10
double projection[3];
11
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]]);
12
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]);
13
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]]);
14
projection[0]/=lengths[wall_index]*lengths[wall_index];
15
projection[1]/=lengths[wall_index]*lengths[wall_index];
16
projection[2]/=lengths[wall_index]*lengths[wall_index];
17
stress_anticlinal[3*wall_index]=lame_first_anticlinal*(lengths[wall_index]-target_lengths[wall_index])*projection[0]/target_lengths[wall_index];
18
stress_anticlinal[3*wall_index+1]=lame_first_anticlinal*(lengths[wall_index]-target_lengths[wall_index])*projection[1]/target_lengths[wall_index];
19
stress_anticlinal[3*wall_index+2]=lame_first_anticlinal*(lengths[wall_index]-target_lengths[wall_index])*projection[2]/target_lengths[wall_index];
20

  
21

  
22
}
src/compute_shapes.cpp (revision 3)
1
#include "constants.h"
2
#include "compute_shapes.h"
3
#include <stdio.h>
4

  
5
/* this function computes the shape of each cell, as its covariance matrix */
6

  
7
void compute_shapes(double *vertices, int *cells_vertices, int *vertices_number_in_each_cell,
8
                    int cell_number, double *barycenters, double *shapes){
9

  
10
int i,j,min,max; double N;
11

  
12
for (i=0;i<cell_number;i++){ shapes[3*i]=0; shapes[3*i+1]=0; shapes[3*i+2]=0;}
13

  
14
for (i=0;i<cell_number;i++){ min=vertices_number_in_each_cell[i]; max=vertices_number_in_each_cell[i+1]; N=(double)(max-min);
15
    for (j=0;j<max-min;j++){
16
        shapes[3*i]+=(vertices[2*cells_vertices[min+j]]-barycenters[2*i])*(vertices[2*cells_vertices[min+j]]-barycenters[2*i]);
17
        shapes[3*i+1]+=(vertices[2*cells_vertices[min+j]+1]-barycenters[2*i+1])*(vertices[2*cells_vertices[min+j]+1]-barycenters[2*i+1]);
18
        shapes[3*i+2]+=(vertices[2*cells_vertices[min+j]+1]-barycenters[2*i+1])*(vertices[2*cells_vertices[min+j]]-barycenters[2*i]);}
19

  
20
    shapes[3*i]/=N; shapes[3*i+1]/=N; shapes[3*i+2]/=N;}
21
}
src/compute_barycenters.h (revision 3)
1
#ifndef COMPUTE_barycenters_H_INCLUDED
2
#define COMPUTE_barycenters_H_INCLUDED
3

  
4
void compute_barycenters(double *vertices, int vertices_number, int *cells_vertices,
5
                     int *vertices_number_in_each_cell, int cell_number, double *barycenters);
6

  
7
#endif // COMPUTE_barycenters_H_INCLUDED
src/compute_cell_walls.h (revision 3)
1
#ifndef COMPUTE_CELL_WALLS_H_INCLUDED
2
#define COMPUTE_CELL_WALLS_H_INCLUDED
3

  
4
void compute_cell_walls(double *vertices, int *cells_vertices, int *vertices_number_in_each_cell,
5
                    int cell_number, int *cell_walls, int *wall_vertices, int wall_number);
6

  
7
#endif // COMPUTE_CELL_WALLS_H_INCLUDED
src/minimize_BFGS.cpp (revision 3)
1
#include "compute_energy_and_forces.h"
2
#include "minimize_BFGS.h"
3
#include <nlopt.h>
4
#include <stdio.h>
5
#include <stdlib.h>
6
#include <cmath>
7

  
8
/* this function minimizes the mechanical energy */
9

  
10
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,
11
                   int *wall_vertices, int wall_number, double *lengths, double *target_lengths, int *boundary, double *relative_pressure){
12

  
13
nlopt_opt opt;
14
opt=nlopt_create(NLOPT_LD_LBFGS,2*vertices_number);
15

  
16
nlopt_set_stopval(opt, -HUGE_VAL);
17
nlopt_set_maxeval(opt, -1.);
18
nlopt_set_maxtime(opt, -1.);
19
nlopt_set_xtol_rel(opt, -1.);
20
nlopt_set_xtol_abs1(opt, 1.e-6);
21
nlopt_set_ftol_rel(opt, -1.);
22
nlopt_set_ftol_abs(opt, -1.);
23

  
24
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;
25
                        int *wall_vertices; int wall_number; double *lengths; double *target_lengths; int *boundary; double *tissue_target; int anisotropy; double *relative_pressure;};
26

  
27
double min_energy; double tissue_target[2]={1.,1.};
28
struct parameters_list parameters;
29
parameters.vertices_number=vertices_number;
30
parameters.cell_number=cell_number;
31
parameters.cells_vertices=cells_vertices;
32
parameters.vertices_number_in_each_cell=vertices_number_in_each_cell;
33
parameters.barycenters=barycenters;
34
parameters.shapes=shapes;
35
parameters.targets=targets;
36
parameters.areas=areas;
37
parameters.cell_areas=cell_areas;
38
parameters.wall_vertices=wall_vertices;
39
parameters.wall_number=wall_number;
40
parameters.lengths=lengths;
41
parameters.target_lengths=target_lengths;
42
parameters.boundary=boundary;
43
parameters.tissue_target=tissue_target;
44
parameters.anisotropy=0;
45
parameters.relative_pressure=relative_pressure;
46
void* pparameters=&parameters;
47

  
48
//double forces[2*vertices_number];
49
//double test=compute_energy_and_forces(0, vertices, forces, pparameters);
50
////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");
51
//double test2;
52
//for (int i=0;i<2*vertices_number;i++){vertices[i]+=1.e-8;
53
//test2=compute_energy_and_forces(0, vertices, forces, pparameters);
54
//test2-=test;
55
//test2*=1.e8;
56
//vertices[i]-=1.e-8;
57
//printf("%f %f %f\n",test2,forces[i],(test2-forces[i])/forces[i]);
58
//}printf("\n");
59

  
60
int reso=nlopt_set_min_objective(opt, compute_energy_and_forces, pparameters);
61
nlopt_optimize(opt, vertices, &min_energy);
62

  
63
tissue_target[0]=0.; tissue_target[1]=0.;
64
int boundary_vertices_number=0.; double tissue_barycenter[2]={0.,0.};
65
for (int i=0;i<vertices_number_in_each_cell[cell_number];i++){
66
    tissue_barycenter[0]+=((double)boundary[i])*vertices[2*cells_vertices[i]];
67
    tissue_barycenter[1]+=((double)boundary[i])*vertices[2*cells_vertices[i]+1];
68
    boundary_vertices_number+=boundary[i];
69
}
70
tissue_barycenter[0]/=boundary_vertices_number;
71
tissue_barycenter[1]/=boundary_vertices_number;
72

  
73
for (int i=0;i<vertices_number_in_each_cell[cell_number];i++){
74
    tissue_target[0]+=((double)boundary[i])*(vertices[2*cells_vertices[i]]-tissue_barycenter[0])*(vertices[2*cells_vertices[i]]-tissue_barycenter[0]);
75
    tissue_target[1]+=((double)boundary[i])*(vertices[2*cells_vertices[i]+1]-tissue_barycenter[1])*(vertices[2*cells_vertices[i]+1]-tissue_barycenter[1]);
76
}
77
tissue_target[0]/=boundary_vertices_number;
78
tissue_target[1]/=boundary_vertices_number;
79

  
80
parameters.tissue_target=tissue_target;
81
parameters.anisotropy=1;
82
pparameters=&parameters;
83

  
84
reso=nlopt_set_min_objective(opt, compute_energy_and_forces, pparameters);
85
nlopt_optimize(opt, vertices, &min_energy);
86
nlopt_destroy(opt);
87

  
88
}
src/cell_division.cpp (revision 3)
1
#include <math.h>
2
#include "cell_division.h"
3
#include "constants.h"
4
#include <stdio.h>
5
#include <stdlib.h>
6
#include <omp.h>
7

  
8
/* this function:   finds the division direction
9
                    computes the position of the new vertices
10
                    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) */
11

  
12
void cell_division(double *vertices, int vertices_number, int *cells_vertices, int *vertices_number_in_each_cell, int cell_number,
13
          double *barycenters, int dividing_cell, int *boundary, double *stress, double *shapes,
14
          int *adjacent_cells_number, int *adjacent_cell1, int *adjacent_cell2, int *daughter_cells_vertices_number, int *cell1, int *cell2, int *wall1, int *wall2){
15

  
16
///bias division
17
//double width=3; double dilation=1.;
18

  
19

  
20
int j;
21

  
22
/* vector of the division line */
23
double division_direction[2];
24
/* positions of new vertices */
25
double new_vertex1[2]={0}; double new_vertex2[2]={0};
26

  
27
double eigenvalue=0; double norm;
28
if (threshold_division>0.0){ /* probabilistic division: the probability to divide a cell following the vector n is proportionnal to
29
                        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 */
30
    double anisotropic_stress[2];
31
    anisotropic_stress[0]=(0.5*stress[3*dividing_cell]-0.5*stress[3*dividing_cell+1])/threshold_division;
32
    anisotropic_stress[1]=0.5*stress[3*dividing_cell+2]/threshold_division;
33
    double probability[20];
34
    for (int i=0;i<20;i++){
35
        probability[i]=exp(cos(i*pi/19.0)*anisotropic_stress[0]+sin(i*pi/19.0)*anisotropic_stress[1]);
36
        if (probability[i]==HUGE_VAL){ printf("division_h\n"); goto division_highest_tension;}
37
    }
38
    probability[0]/=2.; probability[19]/=2.;
39
    for (int i=1;i<20;i++){probability[i]+=probability[i-1];}
40
    for (int i=0;i<20;i++){probability[i]/=probability[19];}
41
    double random=(double)rand()/(double)RAND_MAX; //printf("random %f\n",random);
42
    int random2=rand()%2;
43
    if (random<=probability[0]){ //printf("%i\n",0);
44
        division_direction[0]=1;
45
        division_direction[1]=0;}
46
    for (int i=1;i<20;i++){//printf("%f %f\n",probability[i-1],probability[i]);
47
        if (random<=probability[i] && random>probability[i-1]){
48
            if(random2==0){//printf("%i %i\n",i,0);
49
                division_direction[0]=cos(i*0.5*pi/19.0);
50
                division_direction[1]=sin(i*0.5*pi/19.0);}
51
            if(random2==1){//printf("%i %i\n",i,1);
52
                division_direction[0]=cos(i*0.5*pi/19.0);
53
                division_direction[1]=-sin(i*0.5*pi/19.0);}
54
        }
55
    }
56
}
57

  
58
if (threshold_division==0.0){ /* division following the maximal stress = particular case of the previous rule */
59
    division_highest_tension:; //printf("%f %d %d\n",stress[3*dividing_cell+2],dividing_cell,cell_number);
60
    if (stress[3*dividing_cell+2]!=0){
61
    eigenvalue=(stress[3*dividing_cell]+stress[3*dividing_cell+1])/2.;
62
    eigenvalue+=0.5*sqrt((stress[3*dividing_cell]-stress[3*dividing_cell+1])*(stress[3*dividing_cell]-stress[3*dividing_cell+1])+
63
                    4.*stress[3*dividing_cell+2]*stress[3*dividing_cell+2]);
64
    division_direction[0]=stress[3*dividing_cell+2];
65
    division_direction[1]=-stress[3*dividing_cell]+eigenvalue;
66
    norm=sqrt(division_direction[0]*division_direction[0]+division_direction[1]*division_direction[1]);
67
    division_direction[0]/=norm;
68
    division_direction[1]/=norm;}
69
    else if (stress[3*dividing_cell]>stress[3*dividing_cell+1]){
70
    division_direction[0]=1;
71
    division_direction[1]=0;}
72
    else if (stress[3*dividing_cell]<stress[3*dividing_cell+1]){
73
    division_direction[0]=0;
74
    division_direction[1]=1;}
75
    else if (stress[3*dividing_cell]==stress[3*dividing_cell+1]){
76
    division_direction[0]=rand()/(double)RAND_MAX;
77
    division_direction[1]=sqrt(1.-division_direction[0]*division_direction[0]);}
78
}
79

  
80
double temp_shape[3];
81
temp_shape[0]=shapes[3*dividing_cell];
82
temp_shape[1]=shapes[3*dividing_cell+1];
83
temp_shape[2]=shapes[3*dividing_cell+2];
84
//if ((barycenters[2*dividing_cell]*barycenters[2*dividing_cell])<(width*width)){ temp_shape[0]*=dilation*dilation; temp_shape[2]*=dilation;}
85

  
86
if (threshold_division<0.0){ /* geometrical rule: division following the shortest direction */
87
    double shape_direction[2];
88
    eigenvalue=(temp_shape[0]+temp_shape[1])/2.;
89
    eigenvalue+=0.5*sqrt((temp_shape[0]-temp_shape[1])*(temp_shape[0]-temp_shape[1])+
90
                    4.*temp_shape[2]*temp_shape[2]);
91
    shape_direction[0]=temp_shape[2];
92
    shape_direction[1]=eigenvalue-temp_shape[0];
93
    norm=sqrt(shape_direction[0]*shape_direction[0]+shape_direction[1]*shape_direction[1]);
94
    shape_direction[0]/=norm;
95
    shape_direction[1]/=norm;
96
    division_direction[0]=shape_direction[1];
97
    division_direction[1]=-shape_direction[0];
98
}
99

  
100
/* computation of the new vertices as the intersection between the walls of the cell and the division direction */
101
int l=0,m=0,max=vertices_number_in_each_cell[dividing_cell+1],min=vertices_number_in_each_cell[dividing_cell],previous;
102
double test;
103

  
104

  
105
for (int i=min;i<max;i++){ /* loop over the vertices of the dividing cell */
106
    previous=i-1; if (i==min){previous=max-1;}
107
    test=-division_direction[1]*(vertices[2*cells_vertices[previous]]-barycenters[2*dividing_cell]);
108
    test+=division_direction[0]*(vertices[2*cells_vertices[previous]+1]-barycenters[2*dividing_cell+1]);
109
    test/=-division_direction[1]*(vertices[2*cells_vertices[previous]]-vertices[2*cells_vertices[i]])+
110
            division_direction[0]*(vertices[2*cells_vertices[previous]+1]-vertices[2*cells_vertices[i]+1]);
111

  
112
        /* test if the considered point of the wall is on the division direction */
113
    if ((test<1. && test>=0.) && (m==0 || m==1) && l<2){
114
        //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
115
        l++;
116
        if (l==1){
117
            new_vertex1[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];
118
            new_vertex1[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];
119
            cell1[daughter_cells_vertices_number[0]]=vertices_number-2; daughter_cells_vertices_number[0]+=1;
120
            cell2[daughter_cells_vertices_number[1]]=vertices_number-2; daughter_cells_vertices_number[1]+=1;
121
            wall1[0]=cells_vertices[i];
122
            wall1[1]=cells_vertices[previous];}
123
        if (l==2){
124
            new_vertex2[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];
125
            new_vertex2[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];
126
            cell1[daughter_cells_vertices_number[0]]=vertices_number-1; daughter_cells_vertices_number[0]+=1;
127
            cell2[daughter_cells_vertices_number[1]]=vertices_number-1; daughter_cells_vertices_number[1]+=1;
128
            wall2[0]=cells_vertices[i];
129
            wall2[1]=cells_vertices[previous];}
130
            m=2;
131
        }
132

  
133
    else if ((test<1. && test>=0.) && (m==0 || m==2) && l<2){
134
        //if (test>0.9){test=0.9;} if (test<0.1){test=0.1;}
135
        l++;
136
        if (l==1){
137
            new_vertex1[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];
138
            new_vertex1[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];
139
            cell1[daughter_cells_vertices_number[0]]=vertices_number-2; daughter_cells_vertices_number[0]+=1;
140
            cell2[daughter_cells_vertices_number[1]]=vertices_number-2; daughter_cells_vertices_number[1]+=1;
141
            wall1[0]=cells_vertices[i];
142
            wall1[1]=cells_vertices[previous];}
143
        if (l==2){
144
            new_vertex2[0]=test*(vertices[2*cells_vertices[i]]-vertices[2*cells_vertices[previous]])+vertices[2*cells_vertices[previous]];
145
            new_vertex2[1]=test*(vertices[2*cells_vertices[i]+1]-vertices[2*cells_vertices[previous]+1])+vertices[2*cells_vertices[previous]+1];
146
            cell1[daughter_cells_vertices_number[0]]=vertices_number-1; daughter_cells_vertices_number[0]+=1;
147
            cell2[daughter_cells_vertices_number[1]]=vertices_number-1; daughter_cells_vertices_number[1]+=1;
148
            wall2[0]=cells_vertices[i];
149
            wall2[1]=cells_vertices[previous];}
150
            m=1;
151
        }
152

  
153
    /* write the old vertices in the two daughter cells */
154
if (l==0 || l==2){cell1[daughter_cells_vertices_number[0]]=cells_vertices[i]; daughter_cells_vertices_number[0]+=1;}
155
if (l==1){cell2[daughter_cells_vertices_number[1]]=cells_vertices[i]; daughter_cells_vertices_number[1]+=1;}
156
}
157

  
158
/* if l<2, there has been a problem in the computation of the new vertices */
159
if (l!=2){printf("\n"); printf("division de %d ok,recherche des voisines\n",dividing_cell); exit(EXIT_FAILURE);}
160

  
161
/* the new vertices are saved in the list "vertices" */
162
vertices[2*(vertices_number-2)]=new_vertex1[0]; vertices[2*(vertices_number-2)+1]=new_vertex1[1];
163
vertices[2*(vertices_number-1)]=new_vertex2[0]; vertices[2*(vertices_number-1)+1]=new_vertex2[1];
164

  
165
/* find the cells which are adjacent to the dividing one and have to be modified to include the new vertices */
166
adjacent_cell1[0]=-1; adjacent_cell2[0]=-1; adjacent_cells_number[0]=0;
167
for (int i=0;i<cell_number-1;i++){ if (i==dividing_cell){goto nextcell;}
168
    min=vertices_number_in_each_cell[i]; max=vertices_number_in_each_cell[i+1]; previous=max-min-1;
169
    for (j=0;j<max-min;j++){
170
        if ((cells_vertices[min+j]==wall1[1] && cells_vertices[min+previous]==wall1[0])
171
            ||(cells_vertices[min+j]==wall2[1] && cells_vertices[min+previous]==wall2[0])
172
            ){
173

  
174
            if (adjacent_cells_number[0]==0){adjacent_cell1[0]=i;}
175
            if (adjacent_cells_number[0]==1){adjacent_cell2[0]=i;}
176
            adjacent_cells_number[0]+=1; goto nextcell;
177
            }
178
        previous=j;}
179
    nextcell:;
180
}
181
//printf("%d voisines trouvee(s): %d %d\n",adjacent_cells_number[0],adjacent_cell1[0],adjacent_cell2[0]);
182
}
src/compute_center_of_mass.cpp (revision 3)
1
#include "compute_center_of_mass.h"
2

  
3
/* knowing the vertices in each cell, computation of the center of mass of each cell */
4

  
5
void compute_center_of_mass(double *vertices, int vertices_number, int *cells_vertices,
6
                     int *vertices_number_in_each_cell, int cell_number, double *center_of_mass, double *cell_areas){
7

  
8
int i,j,previous,min;
9

  
10
for (i=0;i<cell_number;i++){center_of_mass[2*i]=0; center_of_mass[2*i+1]=0;}
11

  
12
for (i=0;i<cell_number;i++){ previous=vertices_number_in_each_cell[i+1]-vertices_number_in_each_cell[i]-1; min=vertices_number_in_each_cell[i];
13
    for (j=0;j<vertices_number_in_each_cell[i+1]-vertices_number_in_each_cell[i];j++){
14
        center_of_mass[2*i]+=(vertices[2*cells_vertices[min+j]]+vertices[2*cells_vertices[min+previous]])*
15
        (vertices[2*cells_vertices[min+previous]]*vertices[2*cells_vertices[min+j]+1]-vertices[2*cells_vertices[min+j]]*vertices[2*cells_vertices[min+previous]+1]);
16
        center_of_mass[2*i+1]+=(vertices[2*cells_vertices[min+j]+1]+vertices[2*cells_vertices[min+previous]+1])*
17
        (vertices[2*cells_vertices[min+previous]]*vertices[2*cells_vertices[min+j]+1]-vertices[2*cells_vertices[min+j]]*vertices[2*cells_vertices[min+previous]+1]);
18
        previous=j;}
19
    center_of_mass[2*i]/=6.0*cell_areas[i];
20
    center_of_mass[2*i+1]/=6.0*cell_areas[i];
21
}
22
}
src/compute_stress_periclinal.h (revision 3)
1
#ifndef compute_stress_periclinal_H_INCLUDED
2
#define compute_stress_periclinal_H_INCLUDED
3

  
4
void compute_stress_periclinal(int cell_index, double *vertices, int *vertices_number_in_each_cell, int cell_number, double *stress,
5
                      double *shapes, double *targets, double *cell_areas, int *wall_vertices, int *cell_walls, double *lengths, double *target_lengths);
6

  
7
#endif // compute_stress_periclinal_H_INCLUDED
src/compute_cell_areas.cpp (revision 3)
1
#include <math.h>
2
#include "constants.h"
3
#include "compute_cell_areas.h"
4
#include <stdio.h>
5
#include <stdlib.h>
6

  
7
/* this function computes the areas of each cell, knowing the areas of each triangles between two vertices and the barycenter */
8

  
9
void compute_cell_areas(double *vertices, int *cells_vertices, int *vertices_number_in_each_cell,
10
                    int cell_number, double *areas, double *cell_areas){
11

  
12
int i,j;
13

  
14
/* usual loop over the cells, then over the vertices of each cell */
15
for (i=0;i<cell_number;i++){cell_areas[i]=0;
16
    for (j=vertices_number_in_each_cell[i];j<vertices_number_in_each_cell[i+1];j++){
17
        cell_areas[i]+=areas[j];
18
        }
19
}
20

  
21
}
src/compute_areas.cpp (revision 3)
1
#include <math.h>
2
#include "constants.h"
3
#include "compute_areas.h"
4
#include <stdio.h>
5
#include <stdlib.h>
6

  
7
/* this function computes the areas of the triangles made with two vertices of the same cell and the barycenter of this cell
8
   the area is given by half the cross product of the centered coordinates of the vertices */
9

  
10
void compute_areas(double *vertices, int *cells_vertices, int *vertices_number_in_each_cell,
11
                    int cell_number, double *areas){
12

  
13
int i,j,previous;
14

  
15
/* usual loop over the cells, then over the vertices of each cell */
16
for (i=0;i<cell_number;i++){
17
    for (j=vertices_number_in_each_cell[i];j<vertices_number_in_each_cell[i+1];j++){
18
        previous=j-1; if (j==vertices_number_in_each_cell[i]){ previous=vertices_number_in_each_cell[i+1]-1;}
19
        areas[j]=0.5*(vertices[2*cells_vertices[previous]]*vertices[2*cells_vertices[j]+1]-
20
                    vertices[2*cells_vertices[previous]+1]*vertices[2*cells_vertices[j]]);
21
    }
22
}
23

  
24
}
src/growth_and_chemistry.h (revision 3)
1
#ifndef GROWTH_AND_CHEMISTRY_H_INCLUDED
2
#define GROWTH_AND_CHEMISTRY_H_INCLUDED
3

  
4
int growth_and_chemistry(double t, const double x[], double dx[], void *parameters);
5

  
6
#endif // GROWTH_AND_CHEMISTRY_H_INCLUDED
src/find_boundary_walls.h (revision 3)
1
#ifndef FIND_BOUNDARY_WALLS_H_INCLUDED
2
#define FIND_BOUNDARY_WALLS_H_INCLUDED
3

  
4
void find_boundary_walls(int *cells_vertices, int *vertices_number_in_each_cell, int cell_number,
5
                         int *boundary, int *wall_number);
6

  
7
#endif // FIND_BOUNDARY_WALLS_H_INCLUDED
src/save_data.h (revision 3)
1
#ifndef SAVE_DATA_H_INCLUDED
2
#define SAVE_DATA_H_INCLUDED
3

  
4
void save_data(FILE *f, int vertices_number, double *vertices, int cell_number, int *vertices_number_in_each_cell, int *cells_vertices, int wall_number,
5
               double *barycenters, double *stress_periclinal, double *stress_anticlinal, int *boundary, double *cell_areas, double *lengths, double *shapes,
6
               double *targets, double *growth_rate, int *cell_walls);//, double *growth);
7

  
8
#endif // SAVE_DATA_H_INCLUDED
src/load_initial_state.h (revision 3)
1
#ifndef LOAD_INITIAL_STATE_H_INCLUDED
2
#define LOAD_INITIAL_STATE_H_INCLUDED
3

  
4
void load_initial_state(const char * file_name, int cell_number, int vertices_number, float *vertices, int *vertices_number_in_each_cell, int *cells_vertices);
5

  
6
#endif // LOAD_INITIAL_STATE_H_INCLUDED
src/compute_lengths.cpp (revision 3)
1
#include <math.h>
2
#include "constants.h"
3
#include "compute_lengths.h"
4
#include <stdio.h>
5

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

  
8
void compute_lengths(double *vertices, int *wall_vertices, int wall_number, double *lengths){
9

  
10
for (int i=0;i<wall_number;i++){
11
    lengths[i]=sqrt((vertices[2*wall_vertices[2*i]]-vertices[2*wall_vertices[2*i+1]])*(vertices[2*wall_vertices[2*i]]-vertices[2*wall_vertices[2*i+1]])+
12
                    (vertices[2*wall_vertices[2*i]+1]-vertices[2*wall_vertices[2*i+1]+1])*(vertices[2*wall_vertices[2*i]+1]-vertices[2*wall_vertices[2*i+1]+1]));
13
}
14
}
src/initialize.cpp (revision 3)
1
#include "compute_areas.h"
2
#include "compute_cell_areas.h"
3
#include "compute_barycenters.h"
4
#include "compute_shapes.h"
5
#include "compute_lengths.h"
6
#include "initialize.h"
7
#include "constants.h"
8
#include <stdio.h>
9
#include "compute_cell_walls.h"
10

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

  
13
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,
14
                double *barycenters, double *shapes, double *targets, double *lengths, double *target_lengths, double *areas, double *cell_areas, double *forces, double *stress_periclinal, double *stress_anticlinal){
15

  
16
int i;
17
compute_barycenters(vertices,vertices_number,cells_vertices,vertices_number_in_each_cell,cell_number,barycenters);
18
compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes);
19
compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas);
20
compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas);
21
compute_cell_walls(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, cell_walls, wall_vertices, wall_number);
22
compute_lengths(vertices, wall_vertices, wall_number, lengths);
23
for (i=0;i<cell_number;i++){
24
    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];
25
    stress_periclinal[3*i]=0.; stress_periclinal[3*i+1]=0.; stress_periclinal[3*i+2]=0.;
26
}
27
for (i=0;i<wall_number;i++){
28
    target_lengths[i]=lengths[i]; stress_anticlinal[3*i]=0.; stress_anticlinal[3*i+1]=0.; stress_anticlinal[3*i+2]=0.;
29
}
30
}
src/compute_wall_vertices.h (revision 3)
1
#ifndef COMPUTE_WALL_VERTICES_H_INCLUDED
2
#define COMPUTE_WALL_VERTICES_H_INCLUDED
3

  
4
void compute_wall_vertices(int *cells_vertices, int *vertices_number_in_each_cell, int cell_number, int *wall_vertices);
5

  
6
#endif // COMPUTE_WALL_VERTICES_H_INCLUDED
src/minimize_steepest_descent.h (revision 3)
1
#ifndef MINIMIZE_STEEPEST_DESCENT_H_INCLUDED
2
#define MINIMIZE_STEEPEST_DESCENT_H_INCLUDED
3

  
4
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,
5
                   int *wall_vertices, int wall_number, double *lengths, double *target_lengths);
6

  
7
#endif // MINIMIZE_STEEPEST_DESCENT_H_INCLUDED
src/cell_division2.h (revision 3)
1
#ifndef CELL_DIVISION2_H_INCLUDED
2
#define CELL_DIVISION2_H_INCLUDED
3

  
4
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,
5
          double *barycenters, double *areas, double *cell_areas, double *lengths, double *target_lengths, double *shapes, double *targets,
6
          int dividing_cell, int *boundary, int *border, int adjacent_cells_number, int adjacent_cell1, int adjacent_cell2,
7
          int *daughter_cells_vertices_number, int *cell1, int *cell2, int *wall1, int *wall2);
8

  
9
#endif // CELL_DIVISION2_H_INCLUDED
src/compute_barycenters.cpp (revision 3)
1
#include "compute_barycenters.h"
2
#include <stdio.h>
3
#include <stdlib.h>
4
#include <math.h>
5

  
6
/* knowing the vertices in each cell, computation of the center of each cell */
7

  
8
void compute_barycenters(double *vertices, int vertices_number, int *cells_vertices,
9
                     int *vertices_number_in_each_cell, int cell_number, double *barycenters){
10

  
11
int i,j,min,max; double N;
12

  
13
for (i=0;i<cell_number;i++){barycenters[2*i]=0; barycenters[2*i+1]=0;}
14

  
15
for (i=0;i<cell_number;i++){ min=vertices_number_in_each_cell[i]; max=vertices_number_in_each_cell[i+1]; N=(double)(max-min);
16
    for (j=0;j<max-min;j++){
17
        barycenters[2*i]+=vertices[2*cells_vertices[min+j]];
18
        barycenters[2*i+1]+=vertices[2*cells_vertices[min+j]+1];
19
        }
20
    barycenters[2*i]/=N;
21
    barycenters[2*i+1]/=N;
22
}
23

  
24
}
src/compute_cell_walls.cpp (revision 3)
1
#include <math.h>
2
#include "constants.h"
3
#include "compute_cell_walls.h"
4
#include <stdio.h>
5

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

  
8
void compute_cell_walls(double *vertices, int *cells_vertices, int *vertices_number_in_each_cell,
9
                    int cell_number, int *cell_walls, int *wall_vertices, int wall_number){
10

  
11
int i,j,k,previous;
12

  
13
for (i=0;i<cell_number;i++){
14
    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;}
15
        for (k=0;k<wall_number;k++){
16
            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])){
17
                cell_walls[j]=k;
18
            }
19
        }
20
    }
21
}
22

  
23
}
src/compute_energy_and_forces.h (revision 3)
1
#ifndef COMPUTE_ENERGY_AND_FORCES_H_INCLUDED
2
#define COMPUTE_ENERGY_AND_FORCES_H_INCLUDED
3

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

  
6
#endif // COMPUTE_ENERGY_AND_FORCES_H_INCLUDED
src/step.h (revision 3)
1
#ifndef STEP_H_INCLUDED
2
#define STEP_H_INCLUDED
3

  
4
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,
5
          double *barycenters, double *shapes, double *targets, double *lengths, double *target_lengths, double *areas, double *cell_areas, double *forces, double *stress_periclinal, double *stress_anticlinal,
6
           int *boundary, double *growth_rate, double *relative_pressure);
7

  
8
#endif // STEP_H_INCLUDED
src/constants.h (revision 3)
1
#ifndef CONSTANTS_H_INCLUDED
2
#define CONSTANTS_H_INCLUDED
3

  
4
extern int step_number;
5
extern int max_cell_number;
6
extern double tau1;
7

  
8
extern double pi;
9

  
10
extern double initial_growth;
11
extern double threshold_area;
12
extern double threshold_division;
13

  
14
extern double shear_modulus;
15
extern double lame_first_periclinal;
16
extern double lame_first_anticlinal;
17
extern double pressure;
18

  
19
extern double first_shape_derived_stress;
20
extern double second_shape_derived_stress;
21

  
22
extern double noise_amplitude;
23

  
24
extern int ablated_cell;
25
extern int iteration;
26
extern int iteration_ablation;
27

  
28
extern double growth_differential;
29

  
30
#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)
31

  
32
#endif // CONSTANTS_H_INCLUDED
src/compute_stress_periclinal.cpp (revision 3)
1
#include <math.h>
2
#include "compute_stress_periclinal.h"
3
#include "constants.h"
4
#include <stdio.h>
5

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

  
8
void compute_stress_periclinal(int cell_index, double *vertices, int *vertices_number_in_each_cell, int cell_number, double *stress,
9
                      double *shapes, double *targets, double *cell_areas, int *wall_vertices, int *cell_walls, double *lengths, double *target_lengths){
10

  
11
double trace=targets[3*cell_index]+targets[3*cell_index+1];
12
stress[3*cell_index]=(shear_modulus+lame_first_periclinal)*(shapes[3*cell_index]-targets[3*cell_index])+lame_first_periclinal*(shapes[3*cell_index+1]-targets[3*cell_index+1]);
13
stress[3*cell_index+1]=(shear_modulus+lame_first_periclinal)*(shapes[3*cell_index+1]-targets[3*cell_index+1])+lame_first_periclinal*(shapes[3*cell_index]-targets[3*cell_index]);
14
stress[3*cell_index+2]=shear_modulus*(shapes[3*cell_index+2]-targets[3*cell_index+2]);
15

  
16
stress[3*cell_index]*=2./trace;
17
stress[3*cell_index+1]*=2./trace;
18
stress[3*cell_index+2]*=2./trace;
19

  
20
}
src/compute_stress_anticlinal.h (revision 3)
1
#ifndef COMPUTE_STRESS_ANTICLINAL_H_INCLUDED
2
#define COMPUTE_STRESS_ANTICLINAL_H_INCLUDED
3

  
4
void compute_stress_anticlinal(int wall_index, double *vertices, int *vertices_number_in_each_cell, int cell_number, double *stress_anticlinal,
5
                                int *wall_vertices, int *cell_walls, double *lengths, double *target_lengths);
6

  
7
#endif // COMPUTE_STRESS_ANTICLINAL_H_INCLUDED
src/growth_and_chemistry.cpp (revision 3)
1
#include "constants.h"
2
#include "growth_and_chemistry.h"
3
#include <gsl/gsl_errno.h>
4
#include <math.h>
5
#include <algorithm>
6
#include <gsl/gsl_matrix.h>
7
#include <gsl/gsl_blas.h>
8
#include <gsl/gsl_eigen.h>
9
#include <gsl/gsl_math.h>
10
#include <gsl/gsl_cblas.h>
11

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

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

  
16
struct parameters_list2 {int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *shapes; double *lengths; double *cell_areas; double *barycenters;
17
                         int *cell_walls; int wall_number; int *boundary; int *wall_vertices; double *vertices; double *growth_rate; double *noise;};
18
    
19
/* diagonalize the deformation tensor in order to compute growth */
20
double delta_shapes[3*(((parameters_list2*)parameters)->cell_number)];
21
for (int i=0; i<(((parameters_list2*)parameters)->cell_number);i++){
22
    gsl_matrix *gsl_delta_shapes=gsl_matrix_calloc(2,2);
23

  
24
    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]);
25
    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]);
26

  
27
    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);
28
    if (gsl_vector_get(eigenvalues,0)<0.){ gsl_vector_set(eigenvalues,0,0.);}
29
    if (gsl_vector_get(eigenvalues,1)<0.){ gsl_vector_set(eigenvalues,1,0.);}
30

  
31
    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);
32
    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);
33
    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);
34

  
35
    gsl_vector_free(eigenvalues); gsl_matrix_free(eigenvectors); gsl_matrix_free(gsl_delta_shapes); gsl_eigen_symmv_free(workspace);
36
}
37

  
38
    /* computation of the growth rate:
39
     the growth rate of periclinal walls is directly given by the list "growth_rate", and here we add the noise
40
     the growth of anticlinal walls is the average of the growth of the two neighbouring periclinal walls */
41
double cell_growth[(((parameters_list2*)parameters)->cell_number)];
42
double wall_growth[(((parameters_list2*)parameters)->wall_number)];
43
for (int i=0;i<(((parameters_list2*)parameters)->wall_number);i++){
44
    wall_growth[i]=0.;
45
}
46
for (int i=0;i<(((parameters_list2*)parameters)->cell_number);i++){
47
    //cell_growth[i]=1.;
48
    /* with noise */
49
    cell_growth[i]=(((parameters_list2*)parameters)->growth_rate)[i]*(((parameters_list2*)parameters)->noise)[i];
50
    for (int j=(((parameters_list2*)parameters)->vertices_number_in_each_cell)[i];j<(((parameters_list2*)parameters)->vertices_number_in_each_cell)[i+1];j++){
51
        wall_growth[(((parameters_list2*)parameters)->cell_walls)[j]]+=0.5*cell_growth[i]*(((parameters_list2*)parameters)->noise)[i];
52
        if ((((parameters_list2*)parameters)->boundary)[j]==1){wall_growth[(((parameters_list2*)parameters)->cell_walls)[j]]*=2.;}
53
    }
54
}
55

  
56
    /* growth of the periclinal walls, proportionnal to the deformation, with possiblity to set negative growth, corresponding to compression, to zero
57
     this was not used in the end since it made no difference */
58
for (int i=0;i<(((parameters_list2*)parameters)->cell_number);i++){
59
    dx[3*i]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*
60
               //std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i]);
61
               //(((parameters_list2*)parameters)->shapes)[3*i]-x[3*i];
62
               delta_shapes[3*i];
63
    dx[3*i+1]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*
64
               // std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1]);
65
               // (((parameters_list2*)parameters)->shapes)[3*i+1]-x[3*i+1];
66
              delta_shapes[3*i+1];
67
    dx[3*i+2]=cell_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*
68
              //  std::max(0.,(((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2]);
69
             //   (((parameters_list2*)parameters)->shapes)[3*i+2]-x[3*i+2];
70
               delta_shapes[3*i+2];
71
}
72

  
73
    /* growth of the anticlinal walls */
74
for (int i=0;i<(((parameters_list2*)parameters)->wall_number);i++){
75
        dx[3*(((parameters_list2*)parameters)->cell_number)+i]=wall_growth[i]*(1./(tau1*(((parameters_list2*)parameters)->cell_number)))*
76
                std::max((((parameters_list2*)parameters)->lengths)[i]-x[3*(((parameters_list2*)parameters)->cell_number)+i],0.);
77
}
78

  
79
return GSL_SUCCESS;
80
}
src/main.cpp (revision 3)
1
di#include <iostream>
2
#include <string>
3
#include <stdlib.h>
4
#include <stdio.h>
5
#include <time.h>
6
#include "initialize.h"
7
#include "step.h"
8
#include "constants.h"
9
#include "cell_division.h"
10
#include "cell_division2.h"
11
#include <math.h>
12
#include "load_initial_state.h"
13
#include <fstream>
14
#include "compute_barycenters.h"
15
#include "compute_areas.h"
16
#include "compute_cell_areas.h"
17
#include "compute_shapes.h"
18
#include "compute_stress_periclinal.h"
19
#include "compute_stress_anticlinal.h"
20
#include "find_boundary_walls.h"
21
#include "minimize_BFGS.h"
22
#include <omp.h>
23
#include "compute_cell_areas.h"
24
#include "compute_shapes.h"
25
#include "compute_wall_vertices.h"
26
#include "compute_cell_walls.h"
27
#include <dirent.h>
28
#include "save_data.h"
29
#include <cmath>
30
#include "compute_lengths.h"
31

  
32
using namespace std;
33

  
34
int main(){
35

  
36
omp_set_dynamic(0);
37

  
38
srand(time(NULL));
39
double global_barycenter[2];
40
int init_vertices_number,init_cell_number,cumulated_vertices_number,init_wall_number;
41

  
42
// Pick the initial conditions, obtained from previous simulations
43
//init_cell_number=5; init_vertices_number=8; cumulated_vertices_number=20;
44
//const char* file_name="init/5_8_20.txt";
45
init_cell_number=101; init_vertices_number=200; cumulated_vertices_number=567;
46
const char* file_name="init/101_200_567.txt";
47
//init_cell_number=300; init_vertices_number=598; cumulated_vertices_number=1736;
48
//const char* file_name="init/300_598_1736.txt";
49
//init_cell_number=100; init_vertices_number=200; cumulated_vertices_number=561;
50
//const char* file_name="init/100_200_ablation.txt";
51
float * init_vertices; init_vertices=(float*)calloc(2*init_vertices_number,sizeof(float));
52
int * init_vertices_number_in_each_cell; init_vertices_number_in_each_cell=(int*)calloc(init_cell_number+1,sizeof(int));
53
int * init_cells_vertices; init_cells_vertices=(int*)calloc(cumulated_vertices_number,sizeof(int));
54
load_initial_state(file_name, init_cell_number, init_vertices_number, init_vertices, init_vertices_number_in_each_cell, init_cells_vertices);
55
// initial tissue loaded
56

  
57
/* center the tissue */
58
global_barycenter[1]=0; global_barycenter[2]=0;
59
//#pragma omp parallel shared(global_barycenter, init_vertices, init_vertices_number)
60
//{
61
//#pragma omp for
62
for (int i=0;i<init_vertices_number;i++){global_barycenter[1]+=init_vertices[2*i]/init_vertices_number;
63
                                global_barycenter[2]+=init_vertices[2*i+1]/init_vertices_number;}
64
//#pragma omp for
65
for (int i=0;i<init_vertices_number;i++){init_vertices[2*i]-=global_barycenter[1];
66
                                     init_vertices[2*i+1]-=global_barycenter[2];}
67
//}
68
// tissue centered
69
    
70
/* find the walls which are on the boundary of the tissue */
71
int* init_boundary; init_boundary=(int*)calloc(init_vertices_number_in_each_cell[init_cell_number],sizeof(int));
72
find_boundary_walls(init_cells_vertices, init_vertices_number_in_each_cell, init_cell_number, init_boundary, &init_wall_number);
73
int * init_wall_vertices; init_wall_vertices=(int*)calloc(2*init_wall_number,sizeof(int));
74
compute_wall_vertices( init_cells_vertices, init_vertices_number_in_each_cell, init_cell_number, init_wall_vertices);
75
// boundary identified
76

  
77
// define the parameters over witch there is a loop
78
// threshold=0 -> division following the highest tension
79
// threshold<0 -> division following the shortest axis of the cell
80

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

  
83
//    for (int m=7; m<=7; m++){
84
//    FILE *f=NULL; char file_name[256];
85
//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary.txt",1.-((double)m)*0.1);
86
//    f=fopen(file_name,"a+");
87
//    fprintf(f,"mother vertex daughter x y z\n");
88
//    fclose(f);
89
//    }
90
//    for (int m=7; m<=7; m++){
91
//    FILE *f=NULL; char file_name[256];
92
//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem.txt",1.-((double)m)*0.1);
93
//    f=fopen(file_name,"a+");
94
//    fprintf(f,"mother vertex daughter x y z\n");
95
//    fclose(f);
96
//    }
97
//    for (int m=7; m<=7; m++){
98
//    FILE *f=NULL; char file_name[256];
99
//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary.txt",1.-((double)m)*0.1);
100
//    f=fopen(file_name,"a+");
101
//    fprintf(f,"mother vertex daughter x y z\n");
102
//    fclose(f);
103
//    }
104
//    for (int m=7; m<=7; m++){
105
//    FILE *f=NULL; char file_name[256];
106
//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem.txt",1.-((double)m)*0.1);
107
//    f=fopen(file_name,"a+");
108
//    fprintf(f,"mother vertex daughter x y z\n");
109
//    fclose(f);
110
//    }
111
//    for (int m=7; m<=7; m++){
112
//    FILE *f=NULL; char file_name[256];
113
//    sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f.txt",1.-((double)m)*0.1);
114
//    f=fopen(file_name,"a+");
115
//    fprintf(f,"mother vertex daughter x y z\n");
116
//    fclose(f);
117
//    }
118
//    for (int m=7; m<=7; m++){
119
//    FILE *f=NULL; char file_name[256];
120
//    sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f.txt",1.-((double)m)*0.1);
121
//    f=fopen(file_name,"a+");
122
//    fprintf(f,"mother vertex daughter x y z\n");
123
//    fclose(f);
124
//    }
125
//    for (int m=7; m<=7; m++){
126
//    FILE *f=NULL; char file_name[256];
127
//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary_filter.txt",1.-((double)m)*0.1);
128
//    f=fopen(file_name,"a+");
129
//    fprintf(f,"mother vertex daughter x y z\n");
130
//    fclose(f);
131
//    }
132
//    for (int m=7; m<=7; m++){
133
//    FILE *f=NULL; char file_name[256];
134
//    sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem_filter.txt",1.-((double)m)*0.1);
135
//    f=fopen(file_name,"a+");
136
//    fprintf(f,"mother vertex daughter x y z\n");
137
//    fclose(f);
138
//    }
139
//    for (int m=7; m<=7; m++){
140
//    FILE *f=NULL; char file_name[256];
141
//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary_filter.txt",1.-((double)m)*0.1);
142
//    f=fopen(file_name,"a+");
143
//    fprintf(f,"mother vertex daughter x y z\n");
144
//    fclose(f);
145
//    }
146
//    for (int m=7; m<=7; m++){
147
//    FILE *f=NULL; char file_name[256];
148
//    sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem_filter.txt",1.-((double)m)*0.1);
149
//    f=fopen(file_name,"a+");
150
//    fprintf(f,"mother vertex daughter x y z\n");
151
//    fclose(f);
152
//    }
153
//    for (int m=7; m<=7; m++){
154
//    FILE *f=NULL; char file_name[256];
155
//    sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f_filter.txt",1.-((double)m)*0.1);
156
//    f=fopen(file_name,"a+");
157
//    fprintf(f,"mother vertex daughter x y z\n");
158
//    fclose(f);
159
//    }
160
//    for (int m=7; m<=7; m++){
161
//    FILE *f=NULL; char file_name[256];
162
//    sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f_filter.txt",1.-((double)m)*0.1);
163
//    f=fopen(file_name,"a+");
164
//    fprintf(f,"mother vertex daughter x y z\n");
165
//    fclose(f);
166
//    }
167

  
168
//    for (int m=5; m<=10; m++){
169
//    FILE *f=NULL; char file_name[256];
170
//    sprintf(file_name,"results_v7/pressure/stress_%1.1f.txt",1.+((double)m)*0.1);
171
//    f=fopen(file_name,"a+");
172
//    fprintf(f,"mother vertex daughter x y z\n");
173
//    fclose(f);
174
//    }
175
//    for (int m=5; m<=10; m++){
176
//    FILE *f=NULL; char file_name[256];
177
//    sprintf(file_name,"results_v7/pressure/shape_%1.1f.txt",1.+((double)m)*0.1);
178
//    f=fopen(file_name,"a+");
179
//    fprintf(f,"mother vertex daughter x y z\n");
180
//    fclose(f);
181
//    }
182
//    for (int m=5; m<=10; m++){
183
//    FILE *f=NULL; char file_name[256];
184
//    sprintf(file_name,"results_v7/pressure/stress_%1.1f_filter.txt",1.+((double)m)*0.1);
185
//    f=fopen(file_name,"a+");
186
//    fprintf(f,"mother vertex daughter x y z\n");
187
//    fclose(f);
188
//    }
189
//    for (int m=5; m<=10; m++){
190
//    FILE *f=NULL; char file_name[256];
191
//    sprintf(file_name,"results_v7/pressure/shape_%1.1f_filter.txt",1.+((double)m)*0.1);
192
//    f=fopen(file_name,"a+");
193
//    fprintf(f,"mother vertex daughter x y z\n");
194
//    fclose(f);
195
//    }
196

  
197
    for (int m=0; m<=0; m++){
198
    FILE *f=NULL; char file_name[256];
199
    sprintf(file_name,"results_v7/isotropic_2/stress.txt");
200
    f=fopen(file_name,"a+");
201
    fprintf(f,"mother vertex daughter x y z\n");
202
    fclose(f);
203
    }
204
    for (int m=0; m<=0; m++){
205
    FILE *f=NULL; char file_name[256];
206
    sprintf(file_name,"results_v7/isotropic_2/shape.txt");
207
    f=fopen(file_name,"a+");
208
    fprintf(f,"mother vertex daughter x y z\n");
209
    fclose(f);
210
    }
211
    for (int m=0; m<=0; m++){
212
    FILE *f=NULL; char file_name[256];
213
    sprintf(file_name,"results_v7/isotropic_2/stress_filter.txt");
214
    f=fopen(file_name,"a+");
215
    fprintf(f,"mother vertex daughter x y z\n");
216
    fclose(f);
217
    }
218
    for (int m=0; m<=0; m++){
219
    FILE *f=NULL; char file_name[256];
220
    sprintf(file_name,"results_v7/isotropic_2/shape_filter.txt");
221
    f=fopen(file_name,"a+");
222
    fprintf(f,"mother vertex daughter x y z\n");
223
    fclose(f);
224
    }
225
//
226

  
227

  
228
// parallel loop over the parameters
229
#pragma omp parallel for schedule(dynamic,1) collapse(3)
230
for (int meristem=0; meristem<=19; meristem++){ // several meristem, just the initial conditions are different
231
for (int anisotropy=0; anisotropy<=0; anisotropy++){ // amplitude of the stress anisotropy
232
    //for (int stress_source=1; stress_source<=1; stress_source++){ // origin of the stress: 0=growth, 1=curvature
233
        for (int rule=-1; rule<=0; rule++){
234

  
235
            
236
    ablated_cell=-1; // -1=no ablation, any other number n = the cell number n is killed
237
            
238
    iteration=0;
239
            
240
    threshold_division=(double)rule; // rule used for the cell division
241
            // negative number = cell divided along the short axis
242
            // 0 = cell divided along the maximal stress
243
            // positive number = probabilistic rule based on the mdeviatoric stress
244
            
245
            // define the stress anisotropy
246
//    if (stress_source==0){ // anisotropy comes from growth
247
//        growth_differential=0.999999-((double)anisotropy)*0.1;
248
//        first_shape_derived_stress=3.e-2;
249
//        second_shape_derived_stress=3.e-2;
250
//    }
251
//    if (stress_source==1){ // anisotropy comes from curvature
252
//        growth_differential=0.999999;
253
//        first_shape_derived_stress=3.e-2;
254
//        second_shape_derived_stress=(3.e-2)*(1.-((double)anisotropy)*0.1);
255
//    }
256

  
257

  
258
/* allocation of memory */
259
int vertices_number=init_vertices_number;
260
int cell_number=init_cell_number;
261
int wall_number=init_wall_number;
262
int i;
263
double * vertices; vertices=(double*)calloc(2*vertices_number,sizeof(double)); /* coordinates of the vertices */
264
for (i=0;i<2*vertices_number;i++){vertices[i]=(double)init_vertices[i];}
265
double r; double theta;
266
for (i=0;i<vertices_number;i++){r=0.3*(double)(rand())/(double)RAND_MAX; theta=2.*pi*(double)(rand())/(double)RAND_MAX;
267
    vertices[2*i]+=r*cos(theta); vertices[2*i+1]+=r*sin(theta);
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff