TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
b_utilities.c
Go to the documentation of this file.
1 
25 #include "turtle.h"
26 #include "t_alloc.h"
27 #include "rw_maps.h"
28 #include "geometry.h"
29 #include "g_raster2plvector.h"
30 #include "bigcells2.h"
31 #include "geometry2.h"
32 #include "b_utilities.h"
33 
34 #define D_INIT_VALUE -9999.0
35 //#include "rw_maps.h"
36 //#include "gridded.element.input.geotop.h"
37 
38 //#include "g_raster2plvector.h"
39 DOUBLEMATRIX *doublemap_from_longmap(LONGMATRIX *lmap, long novalue, double dnovalue){
52  long r,c;
53  DOUBLEMATRIX *M;
54  M=new_doublematrix(lmap->nrh,lmap->nch);
55 
56  for (r=lmap->nrl;r<=lmap->nrh;r++){
57  for (c=lmap->ncl;c<=lmap->nch;c++)
58  if (lmap->element[r][c]==novalue) {
59  M->element[r][c]=dnovalue;
60  } else {
61  M->element[r][c]=(double)lmap->element[r][c];
62  }
63  }
64 
65  return M;
66 
67 }
68 
69 DOUBLEVECTOR *interpolete_volume2lines(DOUBLEVECTOR *v,GRID *grid) {
79  DOUBLEVECTOR *ve;
80  long j,i,l,kl,kp,kbond;
81  kbond=grid->boundary_indicator;
82 
83  if (grid->polygons->nh!=v->nh) printf ("Error in interpolete_volume2lines v (quantity vector) size (%ld) is not equal to the numbers of polygons (%ld) !! \n",v->nh,grid->polygons->nh);
84 
85  ve=new_doublevector(grid->lines->nh);
86 
87  for (j=ve->nl;j<=ve->nh;j++) {
88  ve->element[j]=D_INIT_VALUE;
89  }
90 
91  for (i=v->nl;i<=v->nh;i++) {
92  for(l=grid->polygons->element[i]->edge_indices->nl;l<=grid->polygons->element[i]->edge_indices->nh;l++){
93  kl=grid->polygons->element[i]->edge_indices->element[l];
94  kp=grid->links->element[i]->connections->element[l];
95  if (kp==kbond) {
96  if (ve->element[kl]==D_INIT_VALUE) ve->element[kl]=v->element[i];
97  } else {
98  if (ve->element[kl]==D_INIT_VALUE) ve->element[kl]=(v->element[i]+v->element[kp])/2.0;
99  }
100  }
101  }
102  for (j=ve->nl;j<=ve->nh;j++) {
103  if (ve->element[j]==D_INIT_VALUE) printf ("Error in interpolete_volume2lines v (quantity vector) (%ld) at position %ld was not well initialized!! \n",ve->nh,j);
104  }
105 
106  return ve;
107 }
108 
109 
110 
111 DOUBLEMATRIX *get_doublemap_spav(LONGMATRIX *lmap, long novalue, double dnovalue,long dr, long dc, DOUBLEMATRIX *Msource){
126  long r,c;
127  DOUBLEMATRIX *M;
128  M=new_doublematrix(lmap->nrh,lmap->nch);
129 
130  for (r=lmap->nrl;r<=lmap->nrh;r++){
131  for (c=lmap->ncl;c<=lmap->nch;c++) {
132  if ((lmap->element[r][c]==novalue) || (r==lmap->nrl) || (r==lmap->nrh) || (c==lmap->ncl) || (c==lmap->nrh)) {
133  M->element[r][c]=dnovalue;
134  } else if ((Msource->element[r][c]!=dnovalue) && (Msource->element[r+dr][c+dc])!=dnovalue){
135  M->element[r][c]=(Msource->element[r][c]+Msource->element[r+dr][c+dc])/2.0;
136 
137  } else if ((Msource->element[r][c]==dnovalue) && (Msource->element[r+dr][c+dc]!=dnovalue)){
138  M->element[r][c]=Msource->element[r+dr][c+dc];
139 
140  } else if ((Msource->element[r][c]!=dnovalue) && (Msource->element[r+dr][c+dc]==dnovalue)) {
141  M->element[r][c]=Msource->element[r][c];
142 
143  } else {
144  printf("Error 1 in get_doublemap_spav no value : address[%ld][%ld]=%ld and ref_map[%ld][%ld]=%lf ref_map[%ld][%ld]=%lf \n",r,c,lmap->element[r][c],r,c,Msource->element[r][c],r+dr,c+dc,Msource->element[r+dr][c+dc]);
145  return NULL;
146  }
147  }
148  }
149 
150  for (r=M->nrl;r<=M->nrh;r++){
151  for (c=M->ncl;c<=M->nch;c++){
152  if ((lmap->element[r][c]==novalue) && (M->element[r][c]==dnovalue)) {
153  } else if ((lmap->element[r][c]==novalue) || (M->element[r][c]==dnovalue)) {
154  printf("Error 2 in get_doublemap_spav no value : (row=%ld col=%ld) address[%ld][%ld]=%ld and ref_map[%ld][%ld]=%lf ref_map[%ld][%ld]=%lf map[%ld][%ld]=%lf\n",M->nrh,M->nch,r,c,lmap->element[r][c],r,c,Msource->element[r][c],r+dr,c+dc,Msource->element[r+dr][c+dc],r,c,M->element[r][c]);
155  return NULL;
156  }
157  }
158  }
159 
160  return M;
161 
162 }
163 
164 
165 DOUBLEVECTOR *get_doublevector_for_lines(LONGMATRIX *h_addresses,LONGMATRIX *v_addresses,DOUBLEMATRIX *mh,DOUBLEMATRIX *mv,double novalue) {
179  //DOUBLEVECTOR *get_doublevector_from_doublematrix(LONGMATRIX *indices,DOUBLEMATRIX *M, double novalue){
180  long r,c,j;
181  long n_horizontal;
182 
183  DOUBLEVECTOR *v_horizontal,*v_vertical, *v;
184  printf ("0) function get_doublevector_for_lines (novalue=%lf)\n",novalue);
185  v_horizontal=get_doublevector_from_doublematrix(h_addresses,mh,novalue);
186  n_horizontal=v_horizontal->nh;
187  printf ("1) function get_doublevector_for_lines n_horizontal=%ld \n",n_horizontal);
188  for (r=v_addresses->nrl;r<=v_addresses->nrh;r++) {
189  for(c=v_addresses->ncl;c<=v_addresses->nch;c++) {
190  if (mh->element[r][c]!=novalue) v_addresses->element[r][c]=v_addresses->element[r][c]-n_horizontal;
191  }
192  }
193 
194 
195  v_vertical=get_doublevector_from_doublematrix(v_addresses,mv,novalue);
196  v=new_doublevector(v_horizontal->nh+v_vertical->nh);
197  for (j=v->nl;j<=n_horizontal;j++) {
198  v->element[j]=v_horizontal->element[j];
199  }
200  for (j=n_horizontal+1;j<=v->nh;j++) {
201  v->element[j]=v_vertical->element[j-n_horizontal];
202  }
203 
204  for (r=v_addresses->nrl;r<=v_addresses->nrh;r++) {
205  for(c=v_addresses->ncl;r<=v_addresses->nch;c++) {
206  if (mh->element[r][c]!=novalue) v_addresses->element[r][c]=v_addresses->element[r][c]+n_horizontal;
207  }
208 
209  }
210 
211  free_doublevector(v_horizontal);
212  free_doublevector(v_vertical);
213  printf ("function get_doublevector_for_lines n_lines=%ld n_horizontal=%ld \n",v->nh,n_horizontal);
214  return v;
215 }
216 
217 
218 
219 
220 
221 
222 
223 
224 int write_raster_from_doublevector_v2(char *filename, DOUBLEVECTOR *v, long n0, T_INIT *UVref, LONGMATRIX *indices, DOUBLEMATRIX *Mref, short float_type, short map_format) {
245  DOUBLEMATRIX *M;
246  long r,c;
247 
248  for(r=indices->nrl;r<=indices->nrh;r++){
249  for(c=indices->ncl;c<=indices->nch;c++){
250  indices->element[r][c]=indices->element[r][c]-n0;
251  }
252  }
253  M=get_doublematrix_from_doublevector(v,indices,Mref,UVref->V->element[2]);
254  write_map(filename,float_type,map_format,M,UVref);
255 
256  for(r=indices->nrl;r<=indices->nrh;r++){
257  for(c=indices->ncl;c<=indices->nch;c++){
258  indices->element[r][c]=indices->element[r][c]+n0;
259  }
260  }
261  free_doublematrix(M);
262 
263  return 0;
264 
265 }
266 
267 int write_line_map(char *filename, DOUBLEVECTOR *v, long n_horizontal, T_INIT *UVref, LONGMATRIX *indices_h, LONGMATRIX *indices_v,DOUBLEMATRIX *Mref_h,DOUBLEMATRIX *Mref_v, short float_type, short map_format) {
289  char *filename_h=join_strings(filename,"_horizontal_lines");
290  char *filename_v=join_strings(filename,"_vertical_lines");
291  DOUBLEVECTOR *v_vertical, *v_horizontal;
292 
293  long i;
294  v_horizontal=new_doublevector(n_horizontal);
295  v_vertical=new_doublevector(v->nh-n_horizontal);
296  for (i=v_horizontal->nl;i<=v_horizontal->nh;i++) {
297  v_horizontal->element[i]=v->element[i];
298  }
299 
300  for (i=v_vertical->nl;i<=v_vertical->nh;i++) {
301  v_vertical->element[i]=v->element[i+v_horizontal->nh]; /* modified by Emanuele Cordano replace v_vertical->nh with v_horizontal->nh with 9 September 2009 */
302  }
303 
304 
305  write_raster_from_doublevector_v2(filename_h,v_horizontal,0,UVref,indices_h,Mref_h,float_type,map_format);
306  write_raster_from_doublevector_v2(filename_v,v_vertical,n_horizontal,UVref,indices_v,Mref_v,float_type,map_format);
307 
308  free_doublevector(v_horizontal);
309  free_doublevector(v_vertical);
310 
311  //int write_line_map(char *filename, DOUBLEVECTOR *v, long n_horizontal, T_INIT *UVref, LONGMATRIX *indices_h, LONGMATRIX *indices_v,DOUBLEMATRIX *Mref_h,DOUBLEMATRIX *Mref_v,float_type,map_format);
312  free(filename_h);
313  free(filename_v);
314  return 0;
315 
316 
317 }
318 
319 int write_lines(char *filename,long n_lines, long n_horizontal, T_INIT *UVref, LONGMATRIX *indices_h, LONGMATRIX *indices_v,DOUBLEMATRIX *Mref_h,DOUBLEMATRIX *Mref_v, short float_type, short map_format) {
341  long i;
342  DOUBLEVECTOR *v;
343  v=new_doublevector(n_lines);
344  for (i=v->nl;i<=v->nh;i++) {
345  v->element[i]=(double)i;
346  }
347  write_line_map(filename,v,n_horizontal,UVref,indices_h,indices_v,Mref_h,Mref_v,float_type,map_format);
348  free_doublevector(v);
349  return 0;
350 }
351 
352 int write_cell(char *filename, long n_cells, T_INIT *UVref, LONGMATRIX *indices, DOUBLEMATRIX *Mref, short float_type, short map_format) {
367  long i;
368  DOUBLEVECTOR *v;
369  v=new_doublevector(n_cells);
370  for (i=v->nl;i<=v->nh;i++){
371  v->element[i]=(double)i;
372  }
373  write_raster_from_doublevector_v2(filename,v,0,UVref,indices,Mref,float_type,map_format);
374  free_doublevector(v);
375  return 0;
376 }
377 
378 int zeros(DOUBLEVECTOR *v) {
385 long i;
386  for(i=v->nl;i<=v->nh;i++){
387  v->element[i]=0.0;
388  }
389 return 0;
390 }
391 
392 int check_matrices(DOUBLEMATRIX *M1,LONGMATRIX *L1, double dnull, long lnull,short print) {
398  long r,c;
399  int s;
400  if ((M1->nrh!=L1->nrh) || (M1->nch!=L1->nch)) printf("Error:: in check_matrices [%ld,%ld] and M [%ld,%ld] has different sizes! \n",L1->nrh,L1->nch,M1->nrh,M1->nch);
401  s=0;
402  for (r=M1->nrl;r<=M1->nrh;r++) {
403  for (c=M1->ncl;c<=M1->nch;c++) {
404  if ((M1->element[r][c]==dnull) && (L1->element[r][c]==lnull)) {
405  } else if ((M1->element[r][c]==dnull) || (L1->element[r][c]==lnull)) {
406  printf ("Error Function check_matrices: dmatrix and lomatrix [%ld,%ld] have different values %lf and %ld!! \n",r,c,M1->element[r][c],L1->element[r][c]);
407  printf("Neighnourg ponts: dmatrix[%ld][%ld]=%lf and dmatrix[%ld][%ld]=%lf\n",r-1,c,M1->element[r-1][c],r,c-1,M1->element[r][c-1]);
408  return -1;
409  }
410  }
411  }
412  if (print==1) printf("Function check matrix Exit %d\n",s);
413 
414  return s;
415 }
416 
417 double get_value_from_doublevector(long k,DOUBLEVECTOR *v) {
429  double error_v=-9999.0;
430  if ((k<v->nl) || (k>v->nh)) {
431  printf("Warning in function in get_value_from_doublevector: value %ld execeeds the size of the vector %ld %ld, function returns %lf ! \n",k,v->nl,v->nh,error_v);
432  return error_v;
433  } else {
434  return v->element[k];
435  }
436  printf("Error in function in get_value_from_doublevector: value %ld (between %ld %ld) is not read correctly, function returns %lf ! \n",k,v->nl,v->nh,error_v);
437  return error_v;
438 }
439 
440 int write_doublevector_in_a_ascii_ft_format_single_file(char *filename, DOUBLEVECTOR *v) {
441  /*
442  * \author Emanuele Cordano
443  * \date 18 March 2010
444  *
445  *
446  */
447  char *function_name="write_doublevector_in_ft_format";
448  FILE *fd;
449  long l;
450 // if (!v->name) v->name="missing_name";
451  fd=t_fopen(filename,"w");
452  fprintf(fd,"index{1} \n");
453  fprintf(fd,"\n");
454  fprintf(fd,"\n");
455  fprintf(fd,"1: double array STANDARD_NAME {\n");
456  for (l=v->nl;l<v->nh;l++) {
457  fprintf(fd,"%le, \n",v->element[l]);
458  }
459  fprintf(fd,"%le} \n",v->element[v->nh]);
460 
461  t_fclose(fd);
462 
463  return 0;
464 
465 }
466 
467 DOUBLEVECTOR *read_doublevector_from_a_ft_format_single_file (char *filename,short print) {
475  DOUBLEVECTOR *vect=NULL;
476  FILE *fd;
477  char *function_name="read_doublevector_from_a_ft_format_single_file";
478  fd=fopen(filename,"r");
479 
480  if (fd==NULL) {
481  printf("Waring in %s: velocity file is missing",function_name);
482 
483  return NULL;
484  }
485  return NULL;
486  short index=read_index(fd,print);
487  vect=read_doublearray(fd,print);
488  fclose(fd);
489  return vect;
490 }