root / bin / main.cpp @ 1
Historique | Voir | Annoter | Télécharger (37,05 ko)
| 1 | 1 | akiss | di#include <iostream> |
|---|---|---|---|
| 2 | 1 | akiss | #include <string> |
| 3 | 1 | akiss | #include <stdlib.h> |
| 4 | 1 | akiss | #include <stdio.h> |
| 5 | 1 | akiss | #include <time.h> |
| 6 | 1 | akiss | #include "initialize.h" |
| 7 | 1 | akiss | #include "step.h" |
| 8 | 1 | akiss | #include "constants.h" |
| 9 | 1 | akiss | #include "cell_division.h" |
| 10 | 1 | akiss | #include "cell_division2.h" |
| 11 | 1 | akiss | #include <math.h> |
| 12 | 1 | akiss | #include "load_initial_state.h" |
| 13 | 1 | akiss | #include <fstream> |
| 14 | 1 | akiss | #include "compute_barycenters.h" |
| 15 | 1 | akiss | #include "compute_areas.h" |
| 16 | 1 | akiss | #include "compute_cell_areas.h" |
| 17 | 1 | akiss | #include "compute_shapes.h" |
| 18 | 1 | akiss | #include "compute_stress_periclinal.h" |
| 19 | 1 | akiss | #include "compute_stress_anticlinal.h" |
| 20 | 1 | akiss | #include "find_boundary_walls.h" |
| 21 | 1 | akiss | #include "minimize_BFGS.h" |
| 22 | 1 | akiss | #include <omp.h> |
| 23 | 1 | akiss | #include "compute_cell_areas.h" |
| 24 | 1 | akiss | #include "compute_shapes.h" |
| 25 | 1 | akiss | #include "compute_wall_vertices.h" |
| 26 | 1 | akiss | #include "compute_cell_walls.h" |
| 27 | 1 | akiss | #include <dirent.h> |
| 28 | 1 | akiss | #include "save_data.h" |
| 29 | 1 | akiss | #include <cmath> |
| 30 | 1 | akiss | #include "compute_lengths.h" |
| 31 | 1 | akiss | |
| 32 | 1 | akiss | using namespace std; |
| 33 | 1 | akiss | |
| 34 | 1 | akiss | int main(){
|
| 35 | 1 | akiss | |
| 36 | 1 | akiss | omp_set_dynamic(0);
|
| 37 | 1 | akiss | |
| 38 | 1 | akiss | srand(time(NULL));
|
| 39 | 1 | akiss | double global_barycenter[2]; |
| 40 | 1 | akiss | int init_vertices_number,init_cell_number,cumulated_vertices_number,init_wall_number;
|
| 41 | 1 | akiss | |
| 42 | 1 | akiss | // Pick the initial conditions, obtained from previous simulations
|
| 43 | 1 | akiss | //init_cell_number=5; init_vertices_number=8; cumulated_vertices_number=20;
|
| 44 | 1 | akiss | //const char* file_name="init/5_8_20.txt";
|
| 45 | 1 | akiss | init_cell_number=101; init_vertices_number=200; cumulated_vertices_number=567; |
| 46 | 1 | akiss | const char* file_name="init/101_200_567.txt"; |
| 47 | 1 | akiss | //init_cell_number=300; init_vertices_number=598; cumulated_vertices_number=1736;
|
| 48 | 1 | akiss | //const char* file_name="init/300_598_1736.txt";
|
| 49 | 1 | akiss | //init_cell_number=100; init_vertices_number=200; cumulated_vertices_number=561;
|
| 50 | 1 | akiss | //const char* file_name="init/100_200_ablation.txt";
|
| 51 | 1 | akiss | float * init_vertices; init_vertices=(float*)calloc(2*init_vertices_number,sizeof(float)); |
| 52 | 1 | akiss | int * init_vertices_number_in_each_cell; init_vertices_number_in_each_cell=(int*)calloc(init_cell_number+1,sizeof(int)); |
| 53 | 1 | akiss | int * init_cells_vertices; init_cells_vertices=(int*)calloc(cumulated_vertices_number,sizeof(int)); |
| 54 | 1 | akiss | load_initial_state(file_name, init_cell_number, init_vertices_number, init_vertices, init_vertices_number_in_each_cell, init_cells_vertices); |
| 55 | 1 | akiss | // initial tissue loaded
|
| 56 | 1 | akiss | |
| 57 | 1 | akiss | /* center the tissue */
|
| 58 | 1 | akiss | global_barycenter[1]=0; global_barycenter[2]=0; |
| 59 | 1 | akiss | //#pragma omp parallel shared(global_barycenter, init_vertices, init_vertices_number)
|
| 60 | 1 | akiss | //{
|
| 61 | 1 | akiss | //#pragma omp for
|
| 62 | 1 | akiss | for (int i=0;i<init_vertices_number;i++){global_barycenter[1]+=init_vertices[2*i]/init_vertices_number; |
| 63 | 1 | akiss | global_barycenter[2]+=init_vertices[2*i+1]/init_vertices_number;} |
| 64 | 1 | akiss | //#pragma omp for
|
| 65 | 1 | akiss | for (int i=0;i<init_vertices_number;i++){init_vertices[2*i]-=global_barycenter[1]; |
| 66 | 1 | akiss | init_vertices[2*i+1]-=global_barycenter[2];} |
| 67 | 1 | akiss | //}
|
| 68 | 1 | akiss | // tissue centered
|
| 69 | 1 | akiss | |
| 70 | 1 | akiss | /* find the walls which are on the boundary of the tissue */
|
| 71 | 1 | akiss | int* init_boundary; init_boundary=(int*)calloc(init_vertices_number_in_each_cell[init_cell_number],sizeof(int)); |
| 72 | 1 | akiss | find_boundary_walls(init_cells_vertices, init_vertices_number_in_each_cell, init_cell_number, init_boundary, &init_wall_number); |
| 73 | 1 | akiss | int * init_wall_vertices; init_wall_vertices=(int*)calloc(2*init_wall_number,sizeof(int)); |
| 74 | 1 | akiss | compute_wall_vertices( init_cells_vertices, init_vertices_number_in_each_cell, init_cell_number, init_wall_vertices); |
| 75 | 1 | akiss | // boundary identified
|
| 76 | 1 | akiss | |
| 77 | 1 | akiss | // define the parameters over witch there is a loop
|
| 78 | 1 | akiss | // threshold=0 -> division following the highest tension
|
| 79 | 1 | akiss | // threshold<0 -> division following the shortest axis of the cell
|
| 80 | 1 | akiss | |
| 81 | 1 | akiss | // initialize the files where the data will be saved. Yes that's a lot of them
|
| 82 | 1 | akiss | |
| 83 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 84 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 85 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary.txt",1.-((double)m)*0.1);
|
| 86 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 87 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 88 | 1 | akiss | // fclose(f);
|
| 89 | 1 | akiss | // }
|
| 90 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 91 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 92 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem.txt",1.-((double)m)*0.1);
|
| 93 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 94 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 95 | 1 | akiss | // fclose(f);
|
| 96 | 1 | akiss | // }
|
| 97 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 98 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 99 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary.txt",1.-((double)m)*0.1);
|
| 100 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 101 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 102 | 1 | akiss | // fclose(f);
|
| 103 | 1 | akiss | // }
|
| 104 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 105 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 106 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem.txt",1.-((double)m)*0.1);
|
| 107 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 108 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 109 | 1 | akiss | // fclose(f);
|
| 110 | 1 | akiss | // }
|
| 111 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 112 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 113 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f.txt",1.-((double)m)*0.1);
|
| 114 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 115 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 116 | 1 | akiss | // fclose(f);
|
| 117 | 1 | akiss | // }
|
| 118 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 119 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 120 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f.txt",1.-((double)m)*0.1);
|
| 121 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 122 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 123 | 1 | akiss | // fclose(f);
|
| 124 | 1 | akiss | // }
|
| 125 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 126 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 127 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary_filter.txt",1.-((double)m)*0.1);
|
| 128 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 129 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 130 | 1 | akiss | // fclose(f);
|
| 131 | 1 | akiss | // }
|
| 132 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 133 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 134 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem_filter.txt",1.-((double)m)*0.1);
|
| 135 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 136 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 137 | 1 | akiss | // fclose(f);
|
| 138 | 1 | akiss | // }
|
| 139 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 140 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 141 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary_filter.txt",1.-((double)m)*0.1);
|
| 142 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 143 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 144 | 1 | akiss | // fclose(f);
|
| 145 | 1 | akiss | // }
|
| 146 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 147 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 148 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem_filter.txt",1.-((double)m)*0.1);
|
| 149 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 150 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 151 | 1 | akiss | // fclose(f);
|
| 152 | 1 | akiss | // }
|
| 153 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 154 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 155 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f_filter.txt",1.-((double)m)*0.1);
|
| 156 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 157 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 158 | 1 | akiss | // fclose(f);
|
| 159 | 1 | akiss | // }
|
| 160 | 1 | akiss | // for (int m=7; m<=7; m++){
|
| 161 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 162 | 1 | akiss | // sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f_filter.txt",1.-((double)m)*0.1);
|
| 163 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 164 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 165 | 1 | akiss | // fclose(f);
|
| 166 | 1 | akiss | // }
|
| 167 | 1 | akiss | |
| 168 | 1 | akiss | // for (int m=5; m<=10; m++){
|
| 169 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 170 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/stress_%1.1f.txt",1.+((double)m)*0.1);
|
| 171 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 172 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 173 | 1 | akiss | // fclose(f);
|
| 174 | 1 | akiss | // }
|
| 175 | 1 | akiss | // for (int m=5; m<=10; m++){
|
| 176 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 177 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/shape_%1.1f.txt",1.+((double)m)*0.1);
|
| 178 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 179 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 180 | 1 | akiss | // fclose(f);
|
| 181 | 1 | akiss | // }
|
| 182 | 1 | akiss | // for (int m=5; m<=10; m++){
|
| 183 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 184 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/stress_%1.1f_filter.txt",1.+((double)m)*0.1);
|
| 185 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 186 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 187 | 1 | akiss | // fclose(f);
|
| 188 | 1 | akiss | // }
|
| 189 | 1 | akiss | // for (int m=5; m<=10; m++){
|
| 190 | 1 | akiss | // FILE *f=NULL; char file_name[256];
|
| 191 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/shape_%1.1f_filter.txt",1.+((double)m)*0.1);
|
| 192 | 1 | akiss | // f=fopen(file_name,"a+");
|
| 193 | 1 | akiss | // fprintf(f,"mother vertex daughter x y z\n");
|
| 194 | 1 | akiss | // fclose(f);
|
| 195 | 1 | akiss | // }
|
| 196 | 1 | akiss | |
| 197 | 1 | akiss | for (int m=0; m<=0; m++){ |
| 198 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 199 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/stress.txt");
|
| 200 | 1 | akiss | f=fopen(file_name,"a+");
|
| 201 | 1 | akiss | fprintf(f,"mother vertex daughter x y z\n");
|
| 202 | 1 | akiss | fclose(f); |
| 203 | 1 | akiss | } |
| 204 | 1 | akiss | for (int m=0; m<=0; m++){ |
| 205 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 206 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/shape.txt");
|
| 207 | 1 | akiss | f=fopen(file_name,"a+");
|
| 208 | 1 | akiss | fprintf(f,"mother vertex daughter x y z\n");
|
| 209 | 1 | akiss | fclose(f); |
| 210 | 1 | akiss | } |
| 211 | 1 | akiss | for (int m=0; m<=0; m++){ |
| 212 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 213 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/stress_filter.txt");
|
| 214 | 1 | akiss | f=fopen(file_name,"a+");
|
| 215 | 1 | akiss | fprintf(f,"mother vertex daughter x y z\n");
|
| 216 | 1 | akiss | fclose(f); |
| 217 | 1 | akiss | } |
| 218 | 1 | akiss | for (int m=0; m<=0; m++){ |
| 219 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 220 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/shape_filter.txt");
|
| 221 | 1 | akiss | f=fopen(file_name,"a+");
|
| 222 | 1 | akiss | fprintf(f,"mother vertex daughter x y z\n");
|
| 223 | 1 | akiss | fclose(f); |
| 224 | 1 | akiss | } |
| 225 | 1 | akiss | //
|
| 226 | 1 | akiss | |
| 227 | 1 | akiss | |
| 228 | 1 | akiss | // parallel loop over the parameters
|
| 229 | 1 | akiss | #pragma omp parallel for schedule(dynamic,1) collapse(3) |
| 230 | 1 | akiss | for (int meristem=0; meristem<=19; meristem++){ // several meristem, just the initial conditions are different |
| 231 | 1 | akiss | for (int anisotropy=0; anisotropy<=0; anisotropy++){ // amplitude of the stress anisotropy |
| 232 | 1 | akiss | //for (int stress_source=1; stress_source<=1; stress_source++){ // origin of the stress: 0=growth, 1=curvature
|
| 233 | 1 | akiss | for (int rule=-1; rule<=0; rule++){ |
| 234 | 1 | akiss | |
| 235 | 1 | akiss | |
| 236 | 1 | akiss | ablated_cell=-1; // -1=no ablation, any other number n = the cell number n is killed |
| 237 | 1 | akiss | |
| 238 | 1 | akiss | iteration=0;
|
| 239 | 1 | akiss | |
| 240 | 1 | akiss | threshold_division=(double)rule; // rule used for the cell division |
| 241 | 1 | akiss | // negative number = cell divided along the short axis
|
| 242 | 1 | akiss | // 0 = cell divided along the maximal stress
|
| 243 | 1 | akiss | // positive number = probabilistic rule based on the mdeviatoric stress
|
| 244 | 1 | akiss | |
| 245 | 1 | akiss | // define the stress anisotropy
|
| 246 | 1 | akiss | // if (stress_source==0){ // anisotropy comes from growth
|
| 247 | 1 | akiss | // growth_differential=0.999999-((double)anisotropy)*0.1;
|
| 248 | 1 | akiss | // first_shape_derived_stress=3.e-2;
|
| 249 | 1 | akiss | // second_shape_derived_stress=3.e-2;
|
| 250 | 1 | akiss | // }
|
| 251 | 1 | akiss | // if (stress_source==1){ // anisotropy comes from curvature
|
| 252 | 1 | akiss | // growth_differential=0.999999;
|
| 253 | 1 | akiss | // first_shape_derived_stress=3.e-2;
|
| 254 | 1 | akiss | // second_shape_derived_stress=(3.e-2)*(1.-((double)anisotropy)*0.1);
|
| 255 | 1 | akiss | // }
|
| 256 | 1 | akiss | |
| 257 | 1 | akiss | |
| 258 | 1 | akiss | /* allocation of memory */
|
| 259 | 1 | akiss | int vertices_number=init_vertices_number;
|
| 260 | 1 | akiss | int cell_number=init_cell_number;
|
| 261 | 1 | akiss | int wall_number=init_wall_number;
|
| 262 | 1 | akiss | int i;
|
| 263 | 1 | akiss | double * vertices; vertices=(double*)calloc(2*vertices_number,sizeof(double)); /* coordinates of the vertices */ |
| 264 | 1 | akiss | for (i=0;i<2*vertices_number;i++){vertices[i]=(double)init_vertices[i];} |
| 265 | 1 | akiss | double r; double theta; |
| 266 | 1 | akiss | for (i=0;i<vertices_number;i++){r=0.3*(double)(rand())/(double)RAND_MAX; theta=2.*pi*(double)(rand())/(double)RAND_MAX; |
| 267 | 1 | akiss | vertices[2*i]+=r*cos(theta); vertices[2*i+1]+=r*sin(theta); |
| 268 | 1 | akiss | } |
| 269 | 1 | akiss | 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 */ |
| 270 | 1 | akiss | for (i=0;i<cell_number+1;i++){vertices_number_in_each_cell[i]=init_vertices_number_in_each_cell[i];} |
| 271 | 1 | akiss | 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 */ |
| 272 | 1 | akiss | for (i=0;i<vertices_number_in_each_cell[cell_number];i++){cells_vertices[i]=init_cells_vertices[i];} |
| 273 | 1 | akiss | int * wall_vertices; wall_vertices=(int*)calloc(2*wall_number,sizeof(int)); /* list of the vertices of all the walls, wall-per-wall */ |
| 274 | 1 | akiss | for (i=0;i<2*wall_number;i++){wall_vertices[i]=init_wall_vertices[i];} |
| 275 | 1 | akiss | double * lengths; lengths=(double*)calloc(wall_number,sizeof(double)); /* lengths of the anticlinal walls */ |
| 276 | 1 | akiss | double * target_lengths; target_lengths=(double*)calloc(wall_number,sizeof(double)); /* target lengths of the anticlinal walls */ |
| 277 | 1 | akiss | 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 |
| 278 | 1 | akiss | double * barycenters; barycenters=(double*)calloc(cell_number*2,sizeof(double)); /* positions of centers of cells */
|
| 279 | 1 | akiss | double * shapes; shapes=(double*)calloc(cell_number*3,sizeof(double)); /* shapes of cells */ |
| 280 | 1 | akiss | double * targets; targets=(double*)calloc(cell_number*3,sizeof(double)); /* targets of cells */ |
| 281 | 1 | akiss | 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 */ |
| 282 | 1 | akiss | double * cell_areas; cell_areas=(double*)calloc(cell_number, sizeof(double)); /* areas of cells */ |
| 283 | 1 | akiss | double * forces; forces=(double*)calloc(2*vertices_number, sizeof(double)); /* forces on vertices */ |
| 284 | 1 | akiss | double * stress_periclinal; stress_periclinal=(double*)calloc(3*cell_number, sizeof(double)); /* stress on the in-the-plane walls */ |
| 285 | 1 | akiss | double * stress_anticlinal; stress_anticlinal=(double*)calloc(3*wall_number, sizeof(double)); /* stress on the out-the-plane walls */ |
| 286 | 1 | akiss | 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 */ |
| 287 | 1 | akiss | for (i=0;i<vertices_number_in_each_cell[cell_number];i++){boundary[i]=init_boundary[i];} |
| 288 | 1 | akiss | double * growth_rate; growth_rate=(double*)calloc(cell_number,sizeof(double)); /* list of the growth rate of each cell*/ |
| 289 | 1 | akiss | int * border; border=(int*)calloc(cell_number,sizeof(int)); /* indicates if a cell belongs to the boundary of the tissue */ |
| 290 | 1 | akiss | for (i=0; i<cell_number;i++){ growth_rate[i]=1.; border[i]=0; |
| 291 | 1 | akiss | for (int j=vertices_number_in_each_cell[i]; j<vertices_number_in_each_cell[i+1];j++){ |
| 292 | 1 | akiss | if (boundary[j]==1.){ growth_rate[i]=growth_differential; border[i]=1; } |
| 293 | 1 | akiss | } |
| 294 | 1 | akiss | } |
| 295 | 1 | akiss | double * relative_pressure; relative_pressure=(double*)calloc(cell_number,sizeof(double)); /* llist of the pressure in each cell */ |
| 296 | 1 | akiss | for (int i=0; i<cell_number; i++){ relative_pressure[i]=1.;} |
| 297 | 1 | akiss | |
| 298 | 1 | akiss | /* from the vertices, initialize the program by calculating everything that has to be so: lengths, areas, stresses, ... */
|
| 299 | 1 | akiss | initialize(vertices,vertices_number,cells_vertices,vertices_number_in_each_cell,cell_number, wall_vertices, wall_number, cell_walls, |
| 300 | 1 | akiss | barycenters,shapes,targets,lengths,target_lengths,areas,cell_areas,forces,stress_periclinal, stress_anticlinal); |
| 301 | 1 | akiss | |
| 302 | 1 | akiss | /* loop until the tissue reaches the limit size */
|
| 303 | 1 | akiss | int check_division; int total_cell_number=1; |
| 304 | 1 | akiss | while (cell_number<max_cell_number /*&& iteration<step_number*/){check_division=0; |
| 305 | 1 | akiss | //while (iteration<step_number || ablated_cell==-1){check_division=0;
|
| 306 | 1 | akiss | |
| 307 | 1 | akiss | compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters); |
| 308 | 1 | akiss | compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas); |
| 309 | 1 | akiss | compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas); |
| 310 | 1 | akiss | compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes); |
| 311 | 1 | akiss | /* one step of energy minimization and growth */
|
| 312 | 1 | akiss | step(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, wall_vertices, wall_number, cell_walls, |
| 313 | 1 | akiss | barycenters, shapes, targets, lengths, target_lengths, areas, cell_areas, forces, stress_periclinal, stress_anticlinal, |
| 314 | 1 | akiss | boundary, growth_rate, relative_pressure); |
| 315 | 1 | akiss | |
| 316 | 1 | akiss | /* compute the stress in the pericinal walls if it is necesssary for the divisions */
|
| 317 | 1 | akiss | if (threshold_division>=0.){ |
| 318 | 1 | akiss | for (i=0; i<cell_number;i++){ |
| 319 | 1 | akiss | 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); |
| 320 | 1 | akiss | } |
| 321 | 1 | akiss | } |
| 322 | 1 | akiss | |
| 323 | 1 | akiss | /* check if division is required, i.e. if a cell is too big */
|
| 324 | 1 | akiss | for (i=0;i<cell_number;i++){ |
| 325 | 1 | akiss | 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];//} |
| 326 | 1 | akiss | |
| 327 | 1 | akiss | |
| 328 | 1 | akiss | /* test: increase the pressure in the cell just before division, see how it modifies the divisions */
|
| 329 | 1 | akiss | relative_pressure[i]=1.+(double)anisotropy*0.1;//} |
| 330 | 1 | akiss | 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); |
| 331 | 1 | akiss | compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas); |
| 332 | 1 | akiss | compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas); //printf("step2\n");int bug[1]={0};
|
| 333 | 1 | akiss | compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters); |
| 334 | 1 | akiss | compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes); |
| 335 | 1 | akiss | compute_lengths(vertices, wall_vertices, wall_number, lengths); |
| 336 | 1 | akiss | for (int k=0; k<cell_number;k++){ |
| 337 | 1 | akiss | 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); |
| 338 | 1 | akiss | } |
| 339 | 1 | akiss | relative_pressure[i]=1.;
|
| 340 | 1 | akiss | /* end of the test, pressure back to normal */
|
| 341 | 1 | akiss | |
| 342 | 1 | akiss | /* temporal variables to identify the neighbours of the dividing cell */
|
| 343 | 1 | akiss | int adjacent_cells_number[1]; int adjacent_cell1[1]; int adjacent_cell2[1]; |
| 344 | 1 | akiss | int daughter_cells_vertices_number[2]={0,0}; int cell1[10],cell2[10],wall1[2],wall2[2]; |
| 345 | 1 | akiss | |
| 346 | 1 | akiss | /* division, with reallocation of the memory */
|
| 347 | 1 | akiss | cell_number+=1; vertices_number+=2; wall_number+=3; |
| 348 | 1 | akiss | vertices=(double*)realloc(vertices, sizeof(double)*2*vertices_number); |
| 349 | 1 | akiss | forces=(double*)realloc(forces, sizeof(double)*2*vertices_number); |
| 350 | 1 | akiss | vertices_number_in_each_cell=(int*)realloc(vertices_number_in_each_cell, sizeof(int)*(cell_number+1)); |
| 351 | 1 | akiss | barycenters=(double*)realloc(barycenters, sizeof(double)*2*cell_number); |
| 352 | 1 | akiss | shapes=(double*)realloc(shapes, sizeof(double)*3*cell_number); |
| 353 | 1 | akiss | targets=(double*)realloc(targets, sizeof(double)*3*cell_number); |
| 354 | 1 | akiss | stress_periclinal=(double*)realloc(stress_periclinal, sizeof(double)*3*cell_number); |
| 355 | 1 | akiss | stress_anticlinal=(double*)realloc(stress_anticlinal, sizeof(double)*3*wall_number); |
| 356 | 1 | akiss | cell_areas=(double*)realloc(cell_areas, sizeof(double)*cell_number); |
| 357 | 1 | akiss | growth_rate=(double*)realloc(growth_rate,sizeof(double)*cell_number); growth_rate[cell_number-1]=growth_rate[i]; |
| 358 | 1 | akiss | border=(int*)realloc(border,sizeof(int)*cell_number); border[cell_number-1]=0; |
| 359 | 1 | akiss | relative_pressure=(double*)realloc(relative_pressure,sizeof(double)*cell_number); relative_pressure[cell_number-1]=1.; |
| 360 | 1 | akiss | |
| 361 | 1 | akiss | /* find the direction of division and the position of the new vertices */
|
| 362 | 1 | akiss | cell_division(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, |
| 363 | 1 | akiss | barycenters, i, boundary, stress_periclinal, shapes, |
| 364 | 1 | akiss | adjacent_cells_number, adjacent_cell1, adjacent_cell2, daughter_cells_vertices_number, cell1, cell2, wall1, wall2); |
| 365 | 1 | akiss | |
| 366 | 1 | akiss | cells_vertices=(int*)realloc(cells_vertices, sizeof(int)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0])); |
| 367 | 1 | akiss | wall_vertices=(int*)realloc(wall_vertices, sizeof(int)*2*wall_number); |
| 368 | 1 | akiss | lengths=(double*)realloc(lengths, sizeof(double)*(wall_number)); |
| 369 | 1 | akiss | target_lengths=(double*)realloc(target_lengths, sizeof(double)*(wall_number)); |
| 370 | 1 | akiss | cell_walls=(int*)realloc(cell_walls, sizeof(int)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0])); |
| 371 | 1 | akiss | areas=(double*)realloc(areas, sizeof(double)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0])); |
| 372 | 1 | akiss | boundary=(int*)realloc(boundary, sizeof(int)*(vertices_number_in_each_cell[cell_number-1]+4+adjacent_cells_number[0])); |
| 373 | 1 | akiss | |
| 374 | 1 | akiss | /* actually performs the division based on the division identified above */
|
| 375 | 1 | akiss | cell_division2(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, wall_vertices, wall_number, |
| 376 | 1 | akiss | barycenters, areas, cell_areas, lengths, target_lengths, shapes, targets, i, boundary, border, |
| 377 | 1 | akiss | adjacent_cells_number[0], adjacent_cell1[0], adjacent_cell2[0], daughter_cells_vertices_number, cell1, cell2, wall1, wall2); |
| 378 | 1 | akiss | |
| 379 | 1 | akiss | /* update what is modified by the division */
|
| 380 | 1 | akiss | compute_cell_walls(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, cell_walls, wall_vertices, wall_number); |
| 381 | 1 | akiss | compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas); |
| 382 | 1 | akiss | compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas); |
| 383 | 1 | akiss | compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters); |
| 384 | 1 | akiss | compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes); |
| 385 | 1 | akiss | compute_lengths(vertices, wall_vertices, wall_number, lengths); |
| 386 | 1 | akiss | |
| 387 | 1 | akiss | |
| 388 | 1 | akiss | /* minimize the energy after division */
|
| 389 | 1 | akiss | 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); |
| 390 | 1 | akiss | compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas); |
| 391 | 1 | akiss | compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas); |
| 392 | 1 | akiss | compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters); |
| 393 | 1 | akiss | compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes); |
| 394 | 1 | akiss | compute_lengths(vertices, wall_vertices, wall_number, lengths); |
| 395 | 1 | akiss | |
| 396 | 1 | akiss | |
| 397 | 1 | akiss | /* save the history of the divisions, used as by Marion in her statistical analysis */
|
| 398 | 1 | akiss | #pragma omp critical
|
| 399 | 1 | akiss | {
|
| 400 | 1 | akiss | if (save_cell==1){ |
| 401 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 402 | 1 | akiss | |
| 403 | 1 | akiss | /*if (stress_source==1){
|
| 404 | 1 | akiss | if (rule==-1){
|
| 405 | 1 | akiss | sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f.txt",1.-((double)anisotropy)*0.1);
|
| 406 | 1 | akiss | }
|
| 407 | 1 | akiss | if (rule==0){
|
| 408 | 1 | akiss | sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f.txt",1.-((double)anisotropy)*0.1);
|
| 409 | 1 | akiss | }
|
| 410 | 1 | akiss | }
|
| 411 | 1 | akiss | |
| 412 | 1 | akiss | if (stress_source==0){
|
| 413 | 1 | akiss | if (growth_rate[i]==1.){
|
| 414 | 1 | akiss | if (rule==-1){
|
| 415 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);
|
| 416 | 1 | akiss | }
|
| 417 | 1 | akiss | if (rule==0){
|
| 418 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);
|
| 419 | 1 | akiss | }
|
| 420 | 1 | akiss | }
|
| 421 | 1 | akiss | if (growth_rate[i]!=1.){
|
| 422 | 1 | akiss | if (rule==-1){
|
| 423 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);
|
| 424 | 1 | akiss | }
|
| 425 | 1 | akiss | if (rule==0){
|
| 426 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);
|
| 427 | 1 | akiss | }
|
| 428 | 1 | akiss | }
|
| 429 | 1 | akiss | }*/
|
| 430 | 1 | akiss | // if (rule==-1){
|
| 431 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/shape_%1.1f.txt",1.+((double)anisotropy)*0.1);
|
| 432 | 1 | akiss | // }
|
| 433 | 1 | akiss | // if (rule==0){
|
| 434 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/stress_%1.1f.txt",1.+((double)anisotropy)*0.1);
|
| 435 | 1 | akiss | // }
|
| 436 | 1 | akiss | |
| 437 | 1 | akiss | if (rule==-1){ |
| 438 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/shape.txt");
|
| 439 | 1 | akiss | } |
| 440 | 1 | akiss | if (rule==0){ |
| 441 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/stress.txt");
|
| 442 | 1 | akiss | } |
| 443 | 1 | akiss | |
| 444 | 1 | akiss | f=fopen(file_name,"a");
|
| 445 | 1 | akiss | |
| 446 | 1 | akiss | int write=0,count=0,count2=0; |
| 447 | 1 | akiss | while (count<(daughter_cells_vertices_number[1]-1)){ |
| 448 | 1 | akiss | if (write){
|
| 449 | 1 | akiss | 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*/, |
| 450 | 1 | akiss | vertices[2*cell2[count2%daughter_cells_vertices_number[1]]], |
| 451 | 1 | akiss | vertices[2*cell2[count2%daughter_cells_vertices_number[1]]+1]); |
| 452 | 1 | akiss | count++; |
| 453 | 1 | akiss | } |
| 454 | 1 | akiss | write+=((cell2[count2%daughter_cells_vertices_number[1]])==(vertices_number-2)); |
| 455 | 1 | akiss | count2++; |
| 456 | 1 | akiss | } |
| 457 | 1 | akiss | write=0; count=0; count2=0; |
| 458 | 1 | akiss | while (count<(daughter_cells_vertices_number[0]-1)){ |
| 459 | 1 | akiss | if (write){
|
| 460 | 1 | akiss | 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*/, |
| 461 | 1 | akiss | vertices[2*cell1[count2%daughter_cells_vertices_number[0]]], |
| 462 | 1 | akiss | vertices[2*cell1[count2%daughter_cells_vertices_number[0]]+1]); |
| 463 | 1 | akiss | count++; |
| 464 | 1 | akiss | } |
| 465 | 1 | akiss | write+=((cell1[count2%daughter_cells_vertices_number[0]])==(vertices_number-1)); |
| 466 | 1 | akiss | count2++; |
| 467 | 1 | akiss | } |
| 468 | 1 | akiss | fclose(f); |
| 469 | 1 | akiss | } |
| 470 | 1 | akiss | } |
| 471 | 1 | akiss | /* division saved */
|
| 472 | 1 | akiss | |
| 473 | 1 | akiss | |
| 474 | 1 | akiss | /* same as above, but filter out the asymmetric divisions */
|
| 475 | 1 | akiss | #pragma omp critical
|
| 476 | 1 | akiss | {
|
| 477 | 1 | akiss | double area_ratio=cell_areas[i]/(cell_areas[cell_number-1]+cell_areas[i]); if (area_ratio>0.5){ area_ratio=1.-area_ratio;} |
| 478 | 1 | akiss | double aspect_ratio=(shapes[3*i]+shapes[3*i+1])/2.; |
| 479 | 1 | akiss | 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]); |
| 480 | 1 | akiss | aspect_ratio=aspect_ratio/(shapes[3*i]+shapes[3*i+1]-aspect_ratio); aspect_ratio=sqrt(aspect_ratio); |
| 481 | 1 | akiss | |
| 482 | 1 | akiss | if (save_cell==1 && area_ratio>0.4 && aspect_ratio>0.5){ |
| 483 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 484 | 1 | akiss | |
| 485 | 1 | akiss | /*if (stress_source==1){
|
| 486 | 1 | akiss | if (rule==-1){
|
| 487 | 1 | akiss | sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f_filter.txt",1.-((double)anisotropy)*0.1);
|
| 488 | 1 | akiss | }
|
| 489 | 1 | akiss | if (rule==0){
|
| 490 | 1 | akiss | sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f_filter.txt",1.-((double)anisotropy)*0.1);
|
| 491 | 1 | akiss | }
|
| 492 | 1 | akiss | }
|
| 493 | 1 | akiss | |
| 494 | 1 | akiss | if (stress_source==0){
|
| 495 | 1 | akiss | if (growth_rate[i]==1.){
|
| 496 | 1 | akiss | if (rule==-1){
|
| 497 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem_filter.txt",1.-((double)anisotropy)*0.1);
|
| 498 | 1 | akiss | }
|
| 499 | 1 | akiss | if (rule==0){
|
| 500 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem_filter.txt",1.-((double)anisotropy)*0.1);
|
| 501 | 1 | akiss | }
|
| 502 | 1 | akiss | }
|
| 503 | 1 | akiss | if (growth_rate[i]!=1.){
|
| 504 | 1 | akiss | if (rule==-1){
|
| 505 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary_filter.txt",1.-((double)anisotropy)*0.1);
|
| 506 | 1 | akiss | }
|
| 507 | 1 | akiss | if (rule==0){
|
| 508 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary_filter.txt",1.-((double)anisotropy)*0.1);
|
| 509 | 1 | akiss | }
|
| 510 | 1 | akiss | }
|
| 511 | 1 | akiss | }*/
|
| 512 | 1 | akiss | // if (rule==-1){
|
| 513 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/shape_%1.1f_filter.txt",1.+((double)anisotropy)*0.1);
|
| 514 | 1 | akiss | // }
|
| 515 | 1 | akiss | // if (rule==0){
|
| 516 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/stress_%1.1f_filter.txt",1.+((double)anisotropy)*0.1);
|
| 517 | 1 | akiss | // }
|
| 518 | 1 | akiss | |
| 519 | 1 | akiss | if (rule==-1){ |
| 520 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/shape_filter.txt");
|
| 521 | 1 | akiss | } |
| 522 | 1 | akiss | if (rule==0){ |
| 523 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/stress_filter.txt");
|
| 524 | 1 | akiss | } |
| 525 | 1 | akiss | |
| 526 | 1 | akiss | f=fopen(file_name,"a");
|
| 527 | 1 | akiss | |
| 528 | 1 | akiss | int write=0,count=0,count2=0; |
| 529 | 1 | akiss | while (count<(daughter_cells_vertices_number[1]-1)){ |
| 530 | 1 | akiss | if (write){
|
| 531 | 1 | akiss | 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*/, |
| 532 | 1 | akiss | vertices[2*cell2[count2%daughter_cells_vertices_number[1]]], |
| 533 | 1 | akiss | vertices[2*cell2[count2%daughter_cells_vertices_number[1]]+1]); |
| 534 | 1 | akiss | count++; |
| 535 | 1 | akiss | } |
| 536 | 1 | akiss | write+=((cell2[count2%daughter_cells_vertices_number[1]])==(vertices_number-2)); |
| 537 | 1 | akiss | count2++; |
| 538 | 1 | akiss | } |
| 539 | 1 | akiss | write=0; count=0; count2=0; |
| 540 | 1 | akiss | while (count<(daughter_cells_vertices_number[0]-1)){ |
| 541 | 1 | akiss | if (write){
|
| 542 | 1 | akiss | 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*/, |
| 543 | 1 | akiss | vertices[2*cell1[count2%daughter_cells_vertices_number[0]]], |
| 544 | 1 | akiss | vertices[2*cell1[count2%daughter_cells_vertices_number[0]]+1]); |
| 545 | 1 | akiss | count++; |
| 546 | 1 | akiss | } |
| 547 | 1 | akiss | write+=((cell1[count2%daughter_cells_vertices_number[0]])==(vertices_number-1)); |
| 548 | 1 | akiss | count2++; |
| 549 | 1 | akiss | } |
| 550 | 1 | akiss | fclose(f); |
| 551 | 1 | akiss | } |
| 552 | 1 | akiss | } |
| 553 | 1 | akiss | |
| 554 | 1 | akiss | |
| 555 | 1 | akiss | #pragma omp critical
|
| 556 | 1 | akiss | {
|
| 557 | 1 | akiss | if (save_cell==1){ |
| 558 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 559 | 1 | akiss | |
| 560 | 1 | akiss | /*if (stress_source==1){
|
| 561 | 1 | akiss | if (rule==-1){
|
| 562 | 1 | akiss | sprintf(file_name,"results_v7/boundary_curvature/shape_%1.1f.txt",1.-((double)anisotropy)*0.1);
|
| 563 | 1 | akiss | }
|
| 564 | 1 | akiss | if (rule==0){
|
| 565 | 1 | akiss | sprintf(file_name,"results_v7/boundary_curvature/stress_%1.1f.txt",1.-((double)anisotropy)*0.1);
|
| 566 | 1 | akiss | }
|
| 567 | 1 | akiss | }
|
| 568 | 1 | akiss | |
| 569 | 1 | akiss | if (stress_source==0){
|
| 570 | 1 | akiss | if (growth_rate[i]==1.){
|
| 571 | 1 | akiss | if (rule==-1){
|
| 572 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);
|
| 573 | 1 | akiss | }
|
| 574 | 1 | akiss | if (rule==0){
|
| 575 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_meristem.txt",1.-((double)anisotropy)*0.1);
|
| 576 | 1 | akiss | }
|
| 577 | 1 | akiss | }
|
| 578 | 1 | akiss | if (growth_rate[i]!=1.){
|
| 579 | 1 | akiss | if (rule==-1){
|
| 580 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/shape_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);
|
| 581 | 1 | akiss | }
|
| 582 | 1 | akiss | if (rule==0){
|
| 583 | 1 | akiss | sprintf(file_name,"results_v7/boundary_growth/stress_%1.1f_boundary.txt",1.-((double)anisotropy)*0.1);
|
| 584 | 1 | akiss | }
|
| 585 | 1 | akiss | }
|
| 586 | 1 | akiss | }*/
|
| 587 | 1 | akiss | // if (rule==-1){
|
| 588 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/shape_%1.1f.txt",1.+((double)anisotropy)*0.1);
|
| 589 | 1 | akiss | // }
|
| 590 | 1 | akiss | // if (rule==0){
|
| 591 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/stress_%1.1f.txt",1.+((double)anisotropy)*0.1);
|
| 592 | 1 | akiss | // }
|
| 593 | 1 | akiss | |
| 594 | 1 | akiss | if (rule==-1){ |
| 595 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/shape.txt");
|
| 596 | 1 | akiss | } |
| 597 | 1 | akiss | if (rule==0){ |
| 598 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/stress.txt");
|
| 599 | 1 | akiss | } |
| 600 | 1 | akiss | |
| 601 | 1 | akiss | f=fopen(file_name,"a");
|
| 602 | 1 | akiss | |
| 603 | 1 | akiss | int write=0,count=0,count2=0; |
| 604 | 1 | akiss | while (count<(daughter_cells_vertices_number[1]-1)){ |
| 605 | 1 | akiss | if (write){
|
| 606 | 1 | akiss | 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*/, |
| 607 | 1 | akiss | vertices[2*cell2[count2%daughter_cells_vertices_number[1]]], |
| 608 | 1 | akiss | vertices[2*cell2[count2%daughter_cells_vertices_number[1]]+1]); |
| 609 | 1 | akiss | count++; |
| 610 | 1 | akiss | } |
| 611 | 1 | akiss | write+=((cell2[count2%daughter_cells_vertices_number[1]])==(vertices_number-2)); |
| 612 | 1 | akiss | count2++; |
| 613 | 1 | akiss | } |
| 614 | 1 | akiss | write=0; count=0; count2=0; |
| 615 | 1 | akiss | while (count<(daughter_cells_vertices_number[0]-1)){ |
| 616 | 1 | akiss | if (write){
|
| 617 | 1 | akiss | 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*/, |
| 618 | 1 | akiss | vertices[2*cell1[count2%daughter_cells_vertices_number[0]]], |
| 619 | 1 | akiss | vertices[2*cell1[count2%daughter_cells_vertices_number[0]]+1]); |
| 620 | 1 | akiss | count++; |
| 621 | 1 | akiss | } |
| 622 | 1 | akiss | write+=((cell1[count2%daughter_cells_vertices_number[0]])==(vertices_number-1)); |
| 623 | 1 | akiss | count2++; |
| 624 | 1 | akiss | } |
| 625 | 1 | akiss | fclose(f); |
| 626 | 1 | akiss | } |
| 627 | 1 | akiss | } |
| 628 | 1 | akiss | |
| 629 | 1 | akiss | #pragma omp critical
|
| 630 | 1 | akiss | {
|
| 631 | 1 | akiss | double area_ratio=cell_areas[i]/(cell_areas[cell_number-1]+cell_areas[i]); if (area_ratio>0.5){ area_ratio=1.-area_ratio;} |
| 632 | 1 | akiss | double aspect_ratio=(shapes[3*i]+shapes[3*i+1])/2.; |
| 633 | 1 | akiss | 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]); |
| 634 | 1 | akiss | aspect_ratio=aspect_ratio/(shapes[3*i]+shapes[3*i+1]-aspect_ratio); aspect_ratio=sqrt(aspect_ratio); |
| 635 | 1 | akiss | |
| 636 | 1 | akiss | if (save_cell==1){ |
| 637 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 638 | 1 | akiss | |
| 639 | 1 | akiss | |
| 640 | 1 | akiss | if (rule==-1){ |
| 641 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/shape_data.txt");
|
| 642 | 1 | akiss | } |
| 643 | 1 | akiss | if (rule==0){ |
| 644 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/stress_data.txt");
|
| 645 | 1 | akiss | } |
| 646 | 1 | akiss | |
| 647 | 1 | akiss | f=fopen(file_name,"a");
|
| 648 | 1 | akiss | fprintf(f,"%d %f %f %d %f %f %f\n",1000*meristem+total_cell_number,cell_areas[i],aspect_ratio, |
| 649 | 1 | akiss | 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]); |
| 650 | 1 | akiss | fclose(f); |
| 651 | 1 | akiss | } |
| 652 | 1 | akiss | } |
| 653 | 1 | akiss | #pragma omp critical
|
| 654 | 1 | akiss | {
|
| 655 | 1 | akiss | double area_ratio=cell_areas[i]/(cell_areas[cell_number-1]+cell_areas[i]); if (area_ratio>0.5){ area_ratio=1.-area_ratio;} |
| 656 | 1 | akiss | double aspect_ratio=(shapes[3*i]+shapes[3*i+1])/2.; |
| 657 | 1 | akiss | 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]); |
| 658 | 1 | akiss | aspect_ratio=aspect_ratio/(shapes[3*i]+shapes[3*i+1]-aspect_ratio); aspect_ratio=sqrt(aspect_ratio); |
| 659 | 1 | akiss | |
| 660 | 1 | akiss | if (save_cell==1 && area_ratio>0.4 && aspect_ratio>0.5){ |
| 661 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 662 | 1 | akiss | |
| 663 | 1 | akiss | |
| 664 | 1 | akiss | if (rule==-1){ |
| 665 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/shape_filter_data.txt");
|
| 666 | 1 | akiss | } |
| 667 | 1 | akiss | if (rule==0){ |
| 668 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/stress_filter_data.txt");
|
| 669 | 1 | akiss | } |
| 670 | 1 | akiss | |
| 671 | 1 | akiss | f=fopen(file_name,"a");
|
| 672 | 1 | akiss | fprintf(f,"%d %f %f %d %f %f %f\n",1000*meristem+total_cell_number,cell_areas[i],aspect_ratio, |
| 673 | 1 | akiss | 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]); |
| 674 | 1 | akiss | fclose(f); |
| 675 | 1 | akiss | } |
| 676 | 1 | akiss | } |
| 677 | 1 | akiss | |
| 678 | 1 | akiss | |
| 679 | 1 | akiss | |
| 680 | 1 | akiss | total_cell_number+=3;
|
| 681 | 1 | akiss | |
| 682 | 1 | akiss | |
| 683 | 1 | akiss | |
| 684 | 1 | akiss | }} |
| 685 | 1 | akiss | |
| 686 | 1 | akiss | if (check_division==1){ |
| 687 | 1 | akiss | /* minimize the energy after all the divisions*/
|
| 688 | 1 | akiss | 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); |
| 689 | 1 | akiss | |
| 690 | 1 | akiss | } |
| 691 | 1 | akiss | |
| 692 | 1 | akiss | compute_cell_walls(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, cell_walls, wall_vertices, wall_number); |
| 693 | 1 | akiss | compute_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas); |
| 694 | 1 | akiss | compute_cell_areas(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, areas, cell_areas); //printf("main3\n");
|
| 695 | 1 | akiss | compute_barycenters(vertices, vertices_number, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters); |
| 696 | 1 | akiss | compute_shapes(vertices, cells_vertices, vertices_number_in_each_cell, cell_number, barycenters, shapes); |
| 697 | 1 | akiss | |
| 698 | 1 | akiss | for (int i=0; i<cell_number;i++){ |
| 699 | 1 | akiss | 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); |
| 700 | 1 | akiss | } |
| 701 | 1 | akiss | |
| 702 | 1 | akiss | for (int i=0; i<wall_number;i++){ |
| 703 | 1 | akiss | compute_stress_anticlinal(i, vertices, vertices_number_in_each_cell, cell_number, stress_anticlinal, wall_vertices, cell_walls, lengths, target_lengths); |
| 704 | 1 | akiss | } |
| 705 | 1 | akiss | |
| 706 | 1 | akiss | iteration++; |
| 707 | 1 | akiss | |
| 708 | 1 | akiss | // if (cell_number>=300 && ablated_cell==-1){ ablated_cell=m%300; iteration=0;}
|
| 709 | 1 | akiss | |
| 710 | 1 | akiss | /* save all the variables, for my analysis */
|
| 711 | 1 | akiss | if (cell_number>=200){ |
| 712 | 1 | akiss | FILE *f=NULL; char file_name[256]; |
| 713 | 1 | akiss | /*if (stress_source==0){
|
| 714 | 1 | akiss | if (rule==-1){
|
| 715 | 1 | akiss | sprintf(file_name,"results_v7/data/growth_shape_%1.1f_%i_%i.txt",1.-((double)anisotropy)*0.1,meristem,iteration); f=fopen(file_name,"w+");
|
| 716 | 1 | akiss | }
|
| 717 | 1 | akiss | if (rule==0){
|
| 718 | 1 | akiss | sprintf(file_name,"results_v7/data/growth_stress_%1.1f_%i_%i.txt",1.-((double)anisotropy)*0.1,meristem,iteration); f=fopen(file_name,"w+");
|
| 719 | 1 | akiss | }
|
| 720 | 1 | akiss | }
|
| 721 | 1 | akiss | if (stress_source==1){
|
| 722 | 1 | akiss | if (rule==-1){
|
| 723 | 1 | akiss | sprintf(file_name,"results_v7/data/curvature_shape_%1.1f_%i_%i.txt",1.-((double)anisotropy)*0.1,meristem,iteration); f=fopen(file_name,"w+");
|
| 724 | 1 | akiss | }
|
| 725 | 1 | akiss | if (rule==0){
|
| 726 | 1 | akiss | sprintf(file_name,"results_v7/data/curvature_stress_%1.1f_%i_%i.txt",1.-((double)anisotropy)*0.1,meristem,iteration); f=fopen(file_name,"w+");
|
| 727 | 1 | akiss | }
|
| 728 | 1 | akiss | }*/
|
| 729 | 1 | akiss | // if (rule==-1){
|
| 730 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/data/shape_%1.1f_%i_%i.txt",1.+((double)anisotropy)*0.1,meristem,iteration); f=fopen(file_name,"w+");
|
| 731 | 1 | akiss | // }
|
| 732 | 1 | akiss | // if (rule==0){
|
| 733 | 1 | akiss | // sprintf(file_name,"results_v7/pressure/data/stress_%1.1f_%i_%i.txt",1.+((double)anisotropy)*0.1,meristem,iteration); f=fopen(file_name,"w+");
|
| 734 | 1 | akiss | // }
|
| 735 | 1 | akiss | |
| 736 | 1 | akiss | if (rule==-1){ |
| 737 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/data/shape_%i_%i.txt",meristem,iteration); f=fopen(file_name,"w+"); |
| 738 | 1 | akiss | } |
| 739 | 1 | akiss | if (rule==0){ |
| 740 | 1 | akiss | sprintf(file_name,"results_v7/isotropic_2/data/stress_%i_%i.txt",meristem,iteration); f=fopen(file_name,"w+"); |
| 741 | 1 | akiss | } |
| 742 | 1 | akiss | |
| 743 | 1 | akiss | 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); |
| 744 | 1 | akiss | fclose(f); |
| 745 | 1 | akiss | } |
| 746 | 1 | akiss | |
| 747 | 1 | akiss | } |
| 748 | 1 | akiss | |
| 749 | 1 | akiss | free(vertices); free(cells_vertices); free(vertices_number_in_each_cell); free(barycenters); free(shapes); free(targets); free(growth_rate); free(border); free(relative_pressure); |
| 750 | 1 | akiss | free(lengths); free(areas); free(cell_areas); free(forces); free(stress_periclinal); free(stress_anticlinal); free(boundary); free(target_lengths); free(cell_walls); free(wall_vertices); |
| 751 | 1 | akiss | |
| 752 | 1 | akiss | //printf("jeu de param?tre termin? %d\n",m);
|
| 753 | 1 | akiss | |
| 754 | 1 | akiss | }}}//}
|
| 755 | 1 | akiss | |
| 756 | 1 | akiss | free(init_vertices); free(init_vertices_number_in_each_cell); free(init_cells_vertices); free(init_boundary); free(init_wall_vertices); |
| 757 | 1 | akiss | |
| 758 | 1 | akiss | printf("fin du programme\n");
|
| 759 | 1 | akiss | return 0; |
| 760 | 1 | akiss | } |