root / bin / compute_energy_and_forces.cpp @ 1
Historique | Voir | Annoter | Télécharger (9,47 ko)
| 1 | 1 | akiss | #include "compute_energy_and_forces.h" |
|---|---|---|---|
| 2 | 1 | akiss | #include "constants.h" |
| 3 | 1 | akiss | #include "compute_areas.h" |
| 4 | 1 | akiss | #include "compute_cell_areas.h" |
| 5 | 1 | akiss | #include "compute_lengths.h" |
| 6 | 1 | akiss | #include <math.h> |
| 7 | 1 | akiss | #include <stdio.h> |
| 8 | 1 | akiss | #include <nlopt.h> |
| 9 | 1 | akiss | #include <omp.h> |
| 10 | 1 | akiss | |
| 11 | 1 | akiss | /* this function computes the mechanical energy and its derivatives */
|
| 12 | 1 | akiss | |
| 13 | 1 | akiss | double compute_energy_and_forces(unsigned n, const double *vertices, double *forces, void *parameters){ |
| 14 | 1 | akiss | |
| 15 | 1 | akiss | struct parameters_list {int vertex_number; int edge_number; int *edges_vertices; |
| 16 | 1 | akiss | double *lengths; double *rest_lengths; double *edges_width; int *periodizer; int anisotropy; double ref_width; double ref_height;}; |
| 17 | 1 | akiss | double energy=0; |
| 18 | 1 | akiss | |
| 19 | 1 | akiss | compute_lengths((double*)vertices, ((parameters_list*)parameters)->edges_vertices, ((parameters_list*)parameters)->edge_number,
|
| 20 | 1 | akiss | ((parameters_list*)parameters)->lengths, ((parameters_list*)parameters)->periodizer, (double)vertices[2*(((parameters_list*)parameters)->vertex_number)], (double)vertices[2*(((parameters_list*)parameters)->vertex_number)+1]); |
| 21 | 1 | akiss | |
| 22 | 1 | akiss | if (forces){ for (int i=0;i<2*((parameters_list*)parameters)->vertex_number+2;i++){forces[i]=0.;}} |
| 23 | 1 | akiss | |
| 24 | 1 | akiss | |
| 25 | 1 | akiss | if (forces){
|
| 26 | 1 | akiss | forces[2*(((parameters_list*)parameters)->vertex_number)]-=vertices[2*(((parameters_list*)parameters)->vertex_number)+1]; |
| 27 | 1 | akiss | forces[2*(((parameters_list*)parameters)->vertex_number)+1]-=vertices[2*(((parameters_list*)parameters)->vertex_number)]; |
| 28 | 1 | akiss | } |
| 29 | 1 | akiss | /* pressure energy */
|
| 30 | 1 | akiss | energy-=vertices[2*(((parameters_list*)parameters)->vertex_number)]*vertices[2*(((parameters_list*)parameters)->vertex_number)+1]; |
| 31 | 1 | akiss | |
| 32 | 1 | akiss | |
| 33 | 1 | akiss | |
| 34 | 1 | akiss | for (int i=0;i<((parameters_list*)parameters)->edge_number;i++){ |
| 35 | 1 | akiss | /* derivatives of the lengths of the walls */
|
| 36 | 1 | akiss | if (forces){
|
| 37 | 1 | akiss | forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])* |
| 38 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])* |
| 39 | 1 | akiss | (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i] |
| 40 | 1 | akiss | -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i+2])/ |
| 41 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]); |
| 42 | 1 | akiss | forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])* |
| 43 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])* |
| 44 | 1 | akiss | (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i] |
| 45 | 1 | akiss | -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i+2])/ |
| 46 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]); |
| 47 | 1 | akiss | forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])* |
| 48 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])* |
| 49 | 1 | akiss | (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+1] |
| 50 | 1 | akiss | -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+3])/ |
| 51 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]); |
| 52 | 1 | akiss | forces[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])* |
| 53 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])* |
| 54 | 1 | akiss | (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+1] |
| 55 | 1 | akiss | -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+3])/ |
| 56 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]); |
| 57 | 1 | akiss | |
| 58 | 1 | akiss | forces[2*((parameters_list*)parameters)->vertex_number]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])*
|
| 59 | 1 | akiss | ((double)(((parameters_list*)parameters)->periodizer)[4*i]-(double)(((parameters_list*)parameters)->periodizer)[4*i+2])* |
| 60 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])* |
| 61 | 1 | akiss | (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]]+vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i] |
| 62 | 1 | akiss | -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]]-vertices[2*((parameters_list*)parameters)->vertex_number]*(double)(((parameters_list*)parameters)->periodizer)[4*i+2])/ |
| 63 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]); |
| 64 | 1 | akiss | |
| 65 | 1 | akiss | forces[2*((parameters_list*)parameters)->vertex_number+1]+=youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])* |
| 66 | 1 | akiss | ((double)(((parameters_list*)parameters)->periodizer)[4*i+1]-(double)(((parameters_list*)parameters)->periodizer)[4*i+3])* |
| 67 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]-(((parameters_list*)parameters)->rest_lengths)[i])* |
| 68 | 1 | akiss | (vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i]+1]+vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+1] |
| 69 | 1 | akiss | -vertices[2*(((parameters_list*)parameters)->edges_vertices)[2*i+1]+1]-vertices[2*((parameters_list*)parameters)->vertex_number+1]*(double)(((parameters_list*)parameters)->periodizer)[4*i+3])/ |
| 70 | 1 | akiss | ((((parameters_list*)parameters)->lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]*(((parameters_list*)parameters)->rest_lengths)[i]); |
| 71 | 1 | akiss | |
| 72 | 1 | akiss | } |
| 73 | 1 | akiss | /* wall elastic energy */
|
| 74 | 1 | akiss | energy+=0.5*youngs_modulus*((((parameters_list*)parameters)->edges_width)[i])*((((parameters_list*)parameters)->rest_lengths)[i])* |
| 75 | 1 | akiss | pow((((parameters_list*)parameters)->lengths)[i]/(((parameters_list*)parameters)->rest_lengths)[i]-1.,2.); |
| 76 | 1 | akiss | } |
| 77 | 1 | akiss | |
| 78 | 1 | akiss | /* external, anisotropic stretching */
|
| 79 | 1 | akiss | if ((((parameters_list*)parameters)->anisotropy)==1){ |
| 80 | 1 | akiss | |
| 81 | 1 | akiss | if (forces){
|
| 82 | 1 | akiss | forces[2*(((parameters_list*)parameters)->vertex_number)]-=stress_x*(((parameters_list*)parameters)->ref_height);
|
| 83 | 1 | akiss | forces[2*(((parameters_list*)parameters)->vertex_number)+1]-=stress_y*(((parameters_list*)parameters)->ref_width); |
| 84 | 1 | akiss | } |
| 85 | 1 | akiss | |
| 86 | 1 | akiss | energy-=stress_x*(((parameters_list*)parameters)->ref_height)*(vertices[2*(((parameters_list*)parameters)->vertex_number)]-(((parameters_list*)parameters)->ref_width))+
|
| 87 | 1 | akiss | stress_y*(((parameters_list*)parameters)->ref_width)*(vertices[2*(((parameters_list*)parameters)->vertex_number)+1]-(((parameters_list*)parameters)->ref_height)); |
| 88 | 1 | akiss | |
| 89 | 1 | akiss | } |
| 90 | 1 | akiss | |
| 91 | 1 | akiss | return(energy);
|
| 92 | 1 | akiss | |
| 93 | 1 | akiss | } |