Flow123d
convert.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2007 Technical University of Liberec. All rights reserved.
4  *
5  * Please make a following refer to Flow123d on your project site if you use the program for any purpose,
6  * especially for academic research:
7  * Flow123d, Research Centre: Advanced Remedial Technologies, Technical University of Liberec, Czech Republic
8  *
9  * This program is free software; you can redistribute it and/or modify it under the terms
10  * of the GNU General Public License version 3 as published by the Free Software Foundation.
11  *
12  * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
13  * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14  * See the GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License along with this program; if not,
17  * write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 021110-1307, USA.
18  *
19  *
20  * $Id$
21  * $Revision$
22  * $LastChangedBy$
23  * $LastChangedDate$
24  *
25  * @file
26  * @brief Functions for exporting some input files into output formats.
27  * @ingroup io
28  *
29  * Should be replaced by unified output classes
30  *
31  */
32 
33 #include "transport/transport.h"
34 
35 #include "system/system.hh"
36 #include "xio.h"
37 #include "io/output.h"
38 #include "system/math_fce.h"
39 #include "mesh/mesh.h"
40 #include "mesh/boundaries.h"
41 #include "convert.h"
42 
43 #include "mesh/neighbours.h"
44 
45 static void output_convert_to_pos_source(struct Problem *problem);
46 static void output_convert_to_pos_bcd(struct Problem *problem);
47 static void output_convert_to_pos_material(struct Problem *problem);
48 static void output_convert_to_pos_concentration(struct Problem *problem);
49 static void output_convert_to_pos_transport_bcd(struct Problem *problem);
50 /*
51 static void output_transport_convert(struct Problem *problem);
52 static void write_transport_ascii_data(FILE *out, struct Problem *problem, struct TNode **nodes, struct TElement **elements, int time_steps, int ph);
53 static void write_transport_binary_data(FILE *out, struct Problem *problem, struct TNode **nodes, struct TElement **elements, int time_steps, int ph);
54 */
55 
56 //=============================================================================
57 // OUTPUT ROUTINE FOR CONVERTING TO POS
58 //=============================================================================
59 void output_convert_to_pos(struct Problem *problem)
60 {
61  F_ENTRY;
62 
63  ASSERT(!( problem == NULL ),"NULL as argument of function output_convert_to_pos()\n");
64  if( OptGetBool("Output", "Write_output_file", "no") == false )
65  return;
66  if (ConstantDB::getInstance()->getChar("Sources_fname") != NULL)
67  output_convert_to_pos_source(problem);
68  output_convert_to_pos_bcd(problem);
69  output_convert_to_pos_material(problem);
70  /*
71  if (OptGetBool("Transport", "Transport_on", "no") == true)
72  {
73  if (ConstantDB::getInstance()->getChar("Concentration_fname") != NULL)
74  output_convert_to_pos_concentration(problem);
75  if (ConstantDB::getInstance()->getChar("Transport_bcd_fname") != NULL)
76  output_convert_to_pos_transport_bcd(problem);
77  }*/
78 }
79 //=============================================================================
80 // OUTPUT ROUTINE FOR CONVERTING SOURCES TO POS
81 //=============================================================================
82 void output_convert_to_pos_source(struct Problem *problem)
83 {
84  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
85 
86  FILE *out;
87  int li;
88  char* filename = (char*)ConstantDB::getInstance()->getChar("Sources_fname");
89  char dbl_fmt[ 16 ];
90  ElementIter elm;
91  Node *nod;
92 
93  F_ENTRY;
94 
95  ASSERT(!( problem == NULL ),"NULL as argument of function output_convert_to_pos_source()\n");
96  sprintf( dbl_fmt, "%%.%dg ", ConstantDB::getInstance()->getInt("Out_digit"));
97  strcat(filename,".pos");
98  out = xfopen( filename, "wt" );
99  xfprintf( out, "View \"%s - sources\" {\n", OptGetStr("Global", "Description", "No description.") );
100  FOR_ELEMENTS(elm)
101  {
102  switch( elm->type ) {
103  case LINE:
104  xfprintf( out, "SL (" );
105  break;
106  case TRIANGLE:
107  xfprintf( out, "ST (" );
108  break;
109  case TETRAHEDRON:
110  xfprintf( out, "SS (" );
111  break;
112  }
113  for( li = 0; li < elm->n_nodes; li++ ) {
114  nod = elm->node[ li ];
115  xfprintf( out, dbl_fmt, nod->getX() );
116  xfprintf( out, ", " );
117  xfprintf( out, dbl_fmt, nod->getY() );
118  xfprintf( out, ", " );
119  xfprintf( out, dbl_fmt, nod->getZ() );
120  if( li == elm->n_nodes - 1 )
121  xfprintf( out, ") {" );
122  else
123  xfprintf( out, ", " );
124  }
125  for( li = 0; li < elm->n_nodes; li++ ) {
126  xfprintf( out, dbl_fmt, 0.0);
127  if( li == elm->n_nodes - 1 )
128  xfprintf( out, "};\n" );
129  else
130  xfprintf( out, ", " );
131  }
132  }
133  xfprintf( out, "};\n" );
134  xfclose(out);
135 }
136 //=============================================================================
137 // OUTPUT ROUTINE FOR CONVERTING BCDs TO POS
138 //=============================================================================
139 void output_convert_to_pos_bcd(struct Problem *problem)
140 {
141  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
142 
143  const int bcd_types[] = {DIRICHLET,NEUMANN,NEWTON,NEWTON};
144  const char *bcd_names[] = {"Dirichlet","Neumann","Newton P","Newton Sigma"};
145  FILE *out;
146  int li,j,test;
147  // Opravdu je to output - meni se na .pos
148  std::string filename = IONameHandler::get_instance()->get_output_file_name(OptGetStr( "Input", "Boundary", "\\" ));
149  char dbl_fmt[ 16 ];
150  ElementIter elm;
151  Node* nod;
152  struct Boundary *bcd;
153 
154  F_ENTRY;
155 
156  ASSERT(!( problem == NULL ),"NULL as argument of function output_convert_to_pos_bcd()\n");
157  sprintf( dbl_fmt, "%%.%dg ", ConstantDB::getInstance()->getInt("Out_digit"));
158  const std::string& file = filename + ".pos";
159  out = xfopen( file, "wt" );
160  xfprintf( out, "View \"%s - mesh\" {\n", OptGetStr("Global", "Description", "No description.") );
161  FOR_ELEMENTS(elm)
162  {
163  switch( elm->type ) {
164  case LINE:
165  xfprintf( out, "SL (" );
166  break;
167  case TRIANGLE:
168  xfprintf( out, "ST (" );
169  break;
170  case TETRAHEDRON:
171  xfprintf( out, "SS (" );
172  break;
173  }
174  for( li = 0; li < elm->n_nodes; li++ ) {
175  nod = elm->node[ li ];
176  xfprintf( out, dbl_fmt, nod->getX() );
177  xfprintf( out, ", " );
178  xfprintf( out, dbl_fmt, nod->getY() );
179  xfprintf( out, ", " );
180  xfprintf( out, dbl_fmt, nod->getZ() );
181  if( li == elm->n_nodes - 1 )
182  xfprintf( out, ") {" );
183  else
184  xfprintf( out, ", " );
185  }
186  for( li = 0; li < elm->n_nodes; li++ ) {
187  xfprintf( out, dbl_fmt, 0.0);
188  if( li == elm->n_nodes - 1 )
189  xfprintf( out, "};\n" );
190  else
191  xfprintf( out, ", " );
192  }
193  }
194  xfprintf( out, "};\n" );
195  for (j = 0;j < 4; j++)
196  {
197  test = false;
198  FOR_BOUNDARIES(bcd)
199  if (bcd->type == bcd_types[ j ])
200  {
201  test = true;
202  break;
203  }
204  if (test == false)
205  continue;
206  xfprintf( out, "View \"%s - bcd %s\" {\n", OptGetStr("Global", "Description", "No description."), bcd_names[j] );
207  FOR_BOUNDARIES(bcd)
208  {
209  if (bcd->type != bcd_types[ j ]){
210  continue;
211  }
212  switch( bcd->side->shape ) {
213  case xPOINT:
214  xfprintf( out, "SP (" );
215  break;
216  case LINE:
217  xfprintf( out, "SL (" );
218  break;
219  case TRIANGLE:
220  xfprintf( out, "ST (" );
221  break;
222  }
223  for( li = 0; li < bcd->side->n_nodes; li++ ) {
224  nod = bcd->side->node[ li ];
225  xfprintf( out, dbl_fmt, nod->getX() );
226  xfprintf( out, ", " );
227  xfprintf( out, dbl_fmt, nod->getY() );
228  xfprintf( out, ", " );
229  xfprintf( out, dbl_fmt, nod->getZ() );
230  if( li == bcd->side->n_nodes - 1 )
231  xfprintf( out, ") {" );
232  else
233  xfprintf( out, ", " );
234  }
235  for( li = 0; li < bcd->side->n_nodes; li++ ) {
236  switch (bcd->type)
237  {
238  case DIRICHLET: xfprintf( out, dbl_fmt, bcd->scalar);break;
239  case NEUMANN: xfprintf( out, dbl_fmt, bcd->flux);break;
240  case NEWTON:
241  switch (j)
242  {
243  case 2: xfprintf( out, dbl_fmt, bcd->scalar);break;
244  case 3: xfprintf( out, dbl_fmt, bcd->sigma);break;
245  }
246  break;
247  }
248  if( li == bcd->side->n_nodes - 1 )
249  xfprintf( out, "};\n" );
250  else
251  xfprintf( out, ", " );
252  }
253  }
254  xfprintf( out, "};\n" );
255  }
256  xfclose(out);
257 }
258 //=============================================================================
259 // OUTPUT ROUTINE FOR CONVERTING MATERIALS TO POS
260 //=============================================================================
261 void output_convert_to_pos_material(struct Problem *problem)
262 {
263  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
264 
265  FILE* out;
266  int li;
267  char* filename = (char*)ConstantDB::getInstance()->getChar("Material_fname");
268  char dbl_fmt[ 16 ];
269  char int_fmt[ 16 ];
270  ElementIter elm;
271  Node* nod;
272 
273  F_ENTRY;
274 
275  sprintf( dbl_fmt, "%%.%dg ", ConstantDB::getInstance()->getInt("Out_digit"));
276  // sprintf( int_fmt, "%%d ", problem->out_digit );
277  sprintf( int_fmt, "%d ", ConstantDB::getInstance()->getInt("Out_digit"));
278  strcat( filename, ".pos" );
279  out = xfopen( filename, "wt" );
280  xfprintf( out, "View \"%s - materials\" {\n", OptGetStr("Global", "Description", "No description.") );
281  FOR_ELEMENTS(elm)
282  {
283  switch( elm->type ) {
284  case LINE:
285  xfprintf( out, "SL (" );
286  break;
287  case TRIANGLE:
288  xfprintf( out, "ST (" );
289  break;
290  case TETRAHEDRON:
291  xfprintf( out, "SS (" );
292  break;
293  }
294  for( li = 0; li < elm->n_nodes; li++ ) {
295  nod = elm->node[ li ];
296  xfprintf( out, dbl_fmt, nod->getX() );
297  xfprintf( out, ", " );
298  xfprintf( out, dbl_fmt, nod->getY() );
299  xfprintf( out, ", " );
300  xfprintf( out, dbl_fmt, nod->getZ() );
301  if( li == elm->n_nodes - 1 )
302  xfprintf( out, ") {" );
303  else
304  xfprintf( out, ", " );
305  }
306  for( li = 0; li < elm->n_nodes; li++ ) {
307  xfprintf( out, int_fmt, elm->mid );
308  if( li == elm->n_nodes - 1 )
309  xfprintf( out, "};\n" );
310  else
311  xfprintf( out, ", " );
312  }
313  }
314  xfprintf( out, "};\n" );
315  xfclose(out);
316 }
317 //=============================================================================
318 // OUTPUT ROUTINE FOR CONVERTING CONCENTRATIONS TO POS
319 //=============================================================================
320 /*
321 void output_convert_to_pos_concentration(struct Problem *problem)
322 {
323  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
324 
325  FILE *out;
326  int li;
327  char* filename = (char*)ConstantDB::getInstance()->getChar("Concentration_fname");
328  char dbl_fmt[ 16 ];
329  ElementIter elm;
330  Node* nod;
331 
332  F_ENTRY;
333 
334  sprintf( dbl_fmt, "%%.%dg ", ConstantDB::getInstance()->getInt("Out_digit"));
335  strcat(filename,".pos");
336  out = xfopen( filename, "wt" );
337  xfprintf( out, "View \"%s - Concentrations\" {\n", OptGetStr("Global", "Description", "No description.") );
338  FOR_ELEMENTS(elm)
339  {
340  switch( elm->type ) {
341  case LINE:
342  xfprintf( out, "SL (" );
343  break;
344  case TRIANGLE:
345  xfprintf( out, "ST (" );
346  break;
347  case TETRAHEDRON:
348  xfprintf( out, "SS (" );
349  break;
350  }
351  for( li = 0; li < elm->n_nodes; li++ ) {
352  nod = elm->node[ li ];
353  xfprintf( out, dbl_fmt, nod->getX() );
354  xfprintf( out, ", " );
355  xfprintf( out, dbl_fmt, nod->getY() );
356  xfprintf( out, ", " );
357  xfprintf( out, dbl_fmt, nod->getZ() );
358  if( li == elm->n_nodes - 1 )
359  xfprintf( out, ") {" );
360  else
361  xfprintf( out, ", " );
362  }
363  for( li = 0; li < elm->n_nodes; li++ ) {
364  if (elm->start_conc != NULL)
365  xfprintf( out, dbl_fmt, elm->start_conc->conc);
366  else
367  xfprintf( out, dbl_fmt, 0.0);
368  if( li == elm->n_nodes - 1 )
369  xfprintf( out, "};\n" );
370  else
371  xfprintf( out, ", " );
372  }
373  }
374  xfprintf( out, "};\n" );
375  xfclose(out);
376 }
377 //=============================================================================
378 // OUTPUT ROUTINE FOR CONVERTING TRANSPORT_BCDs TO POS
379 //=============================================================================
380 void output_convert_to_pos_transport_bcd(struct Problem *problem)
381 {
382  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
383 
384  FILE *out;
385  int li;
386  char *filename = (char*)ConstantDB::getInstance()->getChar("Transport_bcd_fname");
387  char dbl_fmt[ 16 ];
388  ElementIter elm;
389  Node* nod;
390  struct Boundary *bcd;
391 
392  F_ENTRY;
393 
394  sprintf( dbl_fmt, "%%.%dg ", ConstantDB::getInstance()->getInt("Out_digit"));
395  strcat(filename,".pos");
396  out = xfopen( filename, "wt" );
397  xfprintf( out, "View \"%s - Mesh\" {\n", OptGetStr("Global", "Description", "No description.") );
398  FOR_ELEMENTS(elm)
399  {
400  switch( elm->type ) {
401  case LINE:
402  xfprintf( out, "SL (" );
403  break;
404  case TRIANGLE:
405  xfprintf( out, "ST (" );
406  break;
407  case TETRAHEDRON:
408  xfprintf( out, "SS (" );
409  break;
410  }
411  for( li = 0; li < elm->n_nodes; li++ ) {
412  nod = elm->node[ li ];
413  xfprintf( out, dbl_fmt, nod->getX() );
414  xfprintf( out, ", " );
415  xfprintf( out, dbl_fmt, nod->getY() );
416  xfprintf( out, ", " );
417  xfprintf( out, dbl_fmt, nod->getZ() );
418  if( li == elm->n_nodes - 1 )
419  xfprintf( out, ") {" );
420  else
421  xfprintf( out, ", " );
422  }
423  for( li = 0; li < elm->n_nodes; li++ ) {
424  xfprintf( out, dbl_fmt, 0.0);
425  if( li == elm->n_nodes - 1 )
426  xfprintf( out, "};\n" );
427  else
428  xfprintf( out, ", " );
429  }
430  }
431  xfprintf( out, "};\n" );
432  xfprintf( out, "View \"%s - Transport Boundary Conditions\" {\n", OptGetStr("Global", "Description", "No description.") );
433  FOR_BOUNDARIES(bcd)
434  {
435  switch( bcd->side->shape ) {
436  case xPOINT:
437  xfprintf( out, "SP (" );
438  break;
439  case LINE:
440  xfprintf( out, "SL (" );
441  break;
442  case TRIANGLE:
443  xfprintf( out, "ST (" );
444  break;
445  }
446  for( li = 0; li < bcd->side->n_nodes; li++ ) {
447  nod = bcd->side->node[ li ];
448  xfprintf( out, dbl_fmt, nod->getX() );
449  xfprintf( out, ", " );
450  xfprintf( out, dbl_fmt, nod->getY() );
451  xfprintf( out, ", " );
452  xfprintf( out, dbl_fmt, nod->getZ() );
453  if( li == bcd->side->n_nodes - 1 )
454  xfprintf( out, ") {" );
455  else
456  xfprintf( out, ", " );
457  }
458  for( li = 0; li < bcd->side->n_nodes; li++ ) {
459  xfprintf( out, dbl_fmt, bcd->transport_bcd->conc);
460  if( li == bcd->side->n_nodes - 1 )
461  xfprintf( out, "};\n" );
462  else
463  xfprintf( out, ", " );
464  }
465  }
466  xfprintf( out, "};\n" );
467  xfclose(out);
468 }*/
469 
470 // folowing function seems to be completly WRONG by desing and implementation
471 // - no need to recreate whole Element and Node structures HERE
472 // - wrong allocation
473 //
474 #if 0
475 
476 //==============================================================================
477 // OUTPUT FLOW CONVERT
478 //==============================================================================
479 void output_convert(struct Problem *problem){
480  /* struct TNode{
481  double *scalar;
482  };
483  struct TElement{
484  double *scalar;
485  double **vector;
486  }; */
487  Mesh* mesh;
488  ElementIter ele;
489  Node* nod;
490  struct Node *nodes;
491  struct TElement *elements;
492  FILE *out,*in;
493  char dbl_fmt[ 16 ];
494  char filenamei[255];
495  char filenameo[255];
496  char line[LINE_SIZE];
497  int li,time_steps,j,id,k,ci;
498  double value;
499  bool found;
500 
501  F_ENTRY;
502 
503  mesh = problem->mesh;
504  // nodes = (struct TNode**)xmalloc(sizeof(struct TNode*));
505  // elements = (struct TElement**)xmalloc(sizeof(struct TElement*));
506  sprintf( dbl_fmt, "%%.%dg ", problem->out_digit );
507  sprintf( filenameo,"%s",problem->out_fname);
508  sprintf( filenamei,"%s.tmp",problem->out_fname);
509  out = xfopen( filenameo, "wt" );
510  in = xfopen( filenamei, "rt" );
511 
512  found=skip_to(in, "VALUES");
513  ASSERT(found, "Can not find section: VALUES\n");
514 
515  // find number of time steps
516  time_steps = 0;
517  while(skip_to(in, "Time"))
518  time_steps++;
519 
520  nodes = (struct Node *)xmalloc((mesh->max_nod_id + 1)* sizeof(struct Node));
521  elements = (struct TElement *)xmalloc((mesh->max_elm_id + 1)* sizeof(struct TElement));
522 
523  // allocate and initialize tmp nodes scalar
524  for (li = 0; li <= mesh->max_nod_id; li++){
525  nodes[li].scalar = (double *)xmalloc(time_steps* sizeof(double));
526  for (j = 0; j < time_steps; j++)
527  nodes[li].scalar[j] = 0;
528  }
529 
530 
531  // allocate and initialize tmp element scalar & vector
532  for (li = 0; li <= mesh->max_elm_id; li++){
533  elements[li].scalar = (double *)xmalloc( time_steps * sizeof(double));
534  elements[li].vector = (double **)xmalloc(3 * sizeof(double*));
535  for( ci = 0; ci < 3; ci++ )
536  elements[li].vector[ci] = (double *)xmalloc(time_steps * sizeof(double));
537  for (j = 0; j < time_steps; j++)
538  {
539  elements[li].scalar[j] = 0;
540  for( ci = 0; ci < 3; ci++ )
541  elements[li].vector[ci][j] = 0;
542  }
543  }
544 
545  // now read values for nodes
546  xfclose(in);
547  in = xfopen( filenamei, "rt" );
548  found=skip_to(in, "VALUES");
549  ASSERT(found, "Can not find section: VALUES\n");
550 
551  j = 0;
552  while(skip_to(in, "Nodes")){
553  for (li = 0; li < mesh->n_nodes; li++){
554  xfgets( line, LINE_SIZE - 2, in );
555  id = atoi( xstrtok( line) );
556  value = atof ( xstrtok( NULL) );
557  nodes[id].scalar[j] = value;
558  }
559  j++;
560  }
561 
562  // end of TMP file reading
563  // read values for elements
564  xfclose(in);
565  in = xfopen( filenamei, "rt" );
566  found=skip_to(in, "VALUES");
567  ASSERT(found, "Can not find section: VALUES\n");
568 
569  j = 0;
570  while(skip_to(in, "Elements")){
571  for (li = 0; li < mesh->n_elements(); li++){
572  xfgets( line, LINE_SIZE - 2, in );
573  id = atoi( xstrtok( line) );
574  value = atof ( xstrtok( NULL) );
575  elements[id].scalar[j] = value;
576  for( ci = 0; ci < 3; ci++ ) {
577  value = atof ( xstrtok( NULL) );
578  elements[id].vector[ci][j] = value;
579  }
580  }
581  j++;
582  }
583 
584 
585 
586 
587  write_ascii_header(problem,out);
588 
589  for( k = 0; k < 4; k++){ // p node, p element, pz element, u element
590  switch(k){
591  case 0:
592  xfprintf( out, "View \"%s - p\" {\n", problem->description );
593  break;
594  case 1:
595  xfprintf( out, "View \"%s - pc\" {\n", problem->description );
596  break;
597  case 2:
598  xfprintf( out, "View \"%s - pz\" {\n", problem->description );
599  break;
600  case 3:
601  xfprintf( out, "View \"%s - u\" {\n", problem->description );
602  break; }
603 
604 
605  if(k != 3)
606  {
607  FOR_ELEMENTS( ele ){
608  switch( ele->type ) {
609  case LINE:
610  xfprintf( out, "SL (" );
611  break;
612  case TRIANGLE:
613  xfprintf( out, "ST (" );
614  break;
615  case TETRAHEDRON:
616  xfprintf( out, "SS (" );
617  break;
618  }
619  for( li = 0; li < ele->n_nodes; li++ ) {
620  nod = ele->node[ li ];
621  xfprintf( out, dbl_fmt, nod->x );
622  xfprintf( out, ", " );
623  xfprintf( out, dbl_fmt, nod->y );
624  xfprintf( out, ", " );
625  xfprintf( out, dbl_fmt, nod->z );
626  if( li == ele->n_nodes - 1 )
627  xfprintf( out, ") {" );
628  else
629  xfprintf( out, ", " );
630  }
631 
632  for( li = 0; li < time_steps; li++ ) {
633  for(j = 0; j < ele->n_nodes; j++){
634  switch(k){
635  case 0:
636  nod = ele->node[ j ];
637  xfprintf( out, dbl_fmt,nodes[nod->id].scalar[li]);
638  break;
639  case 1:
640  xfprintf( out, dbl_fmt,elements[ele->id].scalar[li]);
641  break;
642  case 2:
643  xfprintf( out, dbl_fmt,elements[ele->id].scalar[li] + ele->centre[2]);
644  break;
645  }
646  if( li == time_steps - 1 && j == ele->n_nodes - 1)
647  xfprintf( out, "};\n" );
648  else
649  xfprintf( out, ", " );
650  }
651  }
652 
653 
654  } // for elements
655  } // end if k != 3
656  else
657 
658  FOR_ELEMENTS( ele ){
659  xfprintf( out, "VP (" );
660  xfprintf( out, dbl_fmt, ele->centre[ 0 ] );
661  xfprintf( out, ", " );
662  xfprintf( out, dbl_fmt, ele->centre[ 1 ] );
663  xfprintf( out, ", " );
664  xfprintf( out, dbl_fmt, ele->centre[ 2 ] );
665  xfprintf( out, ") {" );
666  for( li = 0; li < time_steps; li++ ){
667  xfprintf( out, dbl_fmt, elements[ele->id].vector[ 0 ][li] );
668  xfprintf( out, ", " );
669  xfprintf( out, dbl_fmt, elements[ele->id].vector[ 1 ][li] );
670  xfprintf( out, ", " );
671  xfprintf( out, dbl_fmt, elements[ele->id].vector[ 2 ][li] );
672  if( li == time_steps - 1)
673  xfprintf( out, "};\n" );
674  else
675  xfprintf( out, ", " );
676  }
677  }
678  xfprintf( out, "};\n" );
679  } //p node, p element, pz element, u element
680 
681  xfclose(in);
682  xfclose(out);
683  xfree(nodes);
684  xfree(elements);
685 
686 }
687 #endif
688 
689 #if 0
690 //==============================================================================
691 // WRITE TRANSPORT ASCII DATA TO THE POS FILE
692 //==============================================================================
693 static void write_transport_ascii_data(FILE *out, struct Problem *problem, struct Node **nodes, struct TElement **elements, int time_steps, int ph)
694 {
695  int k,sbi,li,j,n_subst_;
696  Mesh* mesh;
697  ElementIter ele;
698  Node* nod;
699  char dbl_fmt[ 16 ];
700  double norm,vconc[3];
701 
702  n_subst_ = problem->transport->n_subst_;
703  mesh = problem->mesh;
704  sprintf( dbl_fmt, "%%.%dg ", problem->out_digit );
705 
706  write_ascii_header(problem, out); // header in POS for user view
707 
708  for( k = 0; k < 2; k++){ // nodes and elements
709  for( sbi = 0; sbi < n_subst_; sbi++ ) {
710 
711  xfprintf( out, "View \"Concentration in %s of %s\" {\n",
712  (k == 0) ? "node" : "element",problem->transport->subst_names_[ sbi ] );
713 
714 
715  FOR_ELEMENTS( ele ){
716  switch( ele->type ) {
717  case LINE:
718  xfprintf( out, "SL (" );
719  break;
720  case TRIANGLE:
721  xfprintf( out, "ST (" );
722  break;
723  case TETRAHEDRON:
724  xfprintf( out, "SS (" );
725  break;
726  }
727  FOR_ELEMENT_NODES(ele,li){
728  nod = ele->node[ li ];
729  xfprintf( out, dbl_fmt, nod->x );
730  xfprintf( out, ", " );
731  xfprintf( out, dbl_fmt, nod->y );
732  xfprintf( out, ", " );
733  xfprintf( out, dbl_fmt, nod->z );
734  if( li == ele->n_nodes - 1 )
735  xfprintf( out, ") {" );
736  else
737  xfprintf( out, ", " );
738  }
739 
740  for( li = 0; li < time_steps; li++ ) {
741  FOR_ELEMENT_NODES(ele,j){
742  switch(k){
743  case 0:
744  nod = ele->node[ j ];
745  xfprintf( out, dbl_fmt,nodes[ph][nod->id].conc[sbi][li]);
746  break;
747  case 1:
748  xfprintf( out, dbl_fmt,elements[ph][ele->id].conc[sbi][li]);
749  break;
750  }
751  if( li == time_steps - 1 && j == ele->n_nodes - 1)
752  xfprintf( out, "};\n" );
753  else
754  xfprintf( out, ", " );
755  }
756  }
757 
758  } // END FOR_ELEMENTS
759  xfprintf( out, "};\n" );
760 
761  // CONC VECTOR IN MOBILE ZONE FOR STEADY SATURATED PROBEM
762  if((ph == 0) & (k == 1) & (problem->type == STEADY_SATURATED)){
763  xfprintf( out, "View \"Concentration vector of %s\" {\n",
764  problem->transport->subst_names_[ sbi ] );
765  FOR_ELEMENTS( ele ){
766  xfprintf( out, "VP (" );
767  xfprintf( out, dbl_fmt, ele->centre[ 0 ] );
768  xfprintf( out, ", " );
769  xfprintf( out, dbl_fmt, ele->centre[ 1 ] );
770  xfprintf( out, ", " );
771  xfprintf( out, dbl_fmt, ele->centre[ 2 ] );
772  xfprintf( out, ") {" );
773  norm = vector_length(ele->vector);
774  if(norm > 0){
775  vconc[0] = ele->vector[0]/norm; vconc[1] = ele->vector[1]/norm; vconc[2] = ele->vector[2]/norm;
776  }
777  else
778  scale_vector(vconc,0);
779  for( li = 0; li < time_steps; li++ ){
780  xfprintf( out, dbl_fmt, elements[ph][ele->id].conc[sbi][li] * vconc[0]);
781  xfprintf( out, ", " );
782  xfprintf( out, dbl_fmt, elements[ph][ele->id].conc[sbi][li] * vconc[1]);
783  xfprintf( out, ", " );
784  xfprintf( out, dbl_fmt, elements[ph][ele->id].conc[sbi][li] * vconc[2]);
785  if( li == time_steps - 1)
786  xfprintf( out, "};\n" );
787  else
788  xfprintf( out, ", " );
789  }
790  }
791  xfprintf( out, "};\n" );
792  }
793  // END CONC VECTOR
794  } // end for substances
795  } // end nodes and elements
796 
797 
798 }
799 //==============================================================================
800 // WRITE TRANSPORT BINARY DATA TO THE POS FILE
801 //==============================================================================
802 static void write_transport_binary_data(FILE *out, struct Problem *problem, struct Node **nodes, struct TElement **elements, int time_steps, int ph)
803 {
804  int k,sbi,i,li,j,n_subst_;
805  double ts;
806  Mesh* mesh;
807  ElementIter ele;
808  int one = 1;
809  double norm,vconc[3],vconct[3];
810 
811  n_subst_ = problem->transport->n_subst_;
812  mesh = problem->mesh;
813 
814  xfprintf(out, "$PostFormat\n");
815  xfprintf(out, "%g %d %d\n", 1.4, 1, sizeof(double));
816  xfprintf(out, "$EndPostFormat\n");
817 
818  for( k = 0; k < 2; k++){ // nodes and elements
819  for( sbi = 0; sbi < n_subst_; sbi++ ) {
820  xfprintf( out, "$View\nConcentration^in^%s^of^%s ",
821  (k == 0) ? "node" : "element",problem->transport->subst_names_[ sbi ] );
822  //uprava
823  // mesh->n_tetrahedras = 0;
824  xfprintf(out, "%d 0 0 0 %d 0 0 %d 0 0 0 0 0 %d 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 \n",time_steps,mesh->n_lines,mesh->n_triangles,mesh->n_tetrahedras);
825  xfwrite(&one, sizeof(int), 1, out);
826  ts = 0;
827  for(li=0;li<time_steps;li++){
828  xfwrite(&ts, sizeof(double), 1, out);
829  ts += 1;//problem->time_step;
830  }
831  for(i = 1;i < 4;i++)
832  FOR_ELEMENTS( ele )
833  if(ele->dim == i){ // && ele->dim != 3){
834  write_elm_position_to_binary_output(ele,out);
835  for( j = 0; j < time_steps; j++ )
836  FOR_ELEMENT_NODES(ele,li)
837  xfwrite( (k==0) ? &nodes[ph][ele->node[li]->id].conc[sbi][j] :
838  &elements[ph][ele->id].conc[sbi][j],sizeof(double),1,out);
839  }
840  xfprintf(out, "\n$EndView\n");
841  }
842  }
843  // Vector conc
844  for( sbi = 0; sbi < n_subst_; sbi++ ) {
845  xfprintf( out, "$View\nConcentration^vector^of^%s ",problem->transport->subst_names_[ sbi ] );
846  xfprintf(out, "%d 0 %d 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",time_steps,mesh->n_elements());
847  xfwrite(&one, sizeof(int), 1, out);
848  ts = 0;
849  for(li=0;li<time_steps;li++){
850  xfwrite(&ts, sizeof(double), 1, out);
851  ts += 1;//problem->time_step;
852  }
853  FOR_ELEMENTS( ele ){
854  xfwrite(&ele->centre[0],sizeof(double),3,out);
855  norm = vector_length(ele->vector);
856  if(norm > 0){
857  vconc[0] = ele->vector[0]/norm; vconc[1] = ele->vector[1]/norm; vconc[2] = ele->vector[2]/norm;
858  }
859  else
860  scale_vector(vconc,0);
861  for( j = 0; j < time_steps; j++ ){
862  vconct[0] = vconc[0] * elements[ph][ele->id].conc[sbi][j];
863  vconct[1] = vconc[1] * elements[ph][ele->id].conc[sbi][j];
864  vconct[2] = vconc[2] * elements[ph][ele->id].conc[sbi][j];
865  xfwrite(&vconct[0],sizeof(double),3,out);
866  }
867  }
868  xfprintf(out, "\n$EndView\n");
869  }
870 }
871 
872 //==============================================================================
873 // OUTPUT TRANSPORT CONVERT
874 //==============================================================================
875 void output_transport_convert(struct Problem *problem){
876 
877  Mesh* mesh = problem->mesh;
878  struct Node** nodes;
879  struct TElement** elements;
880  FILE **out,**in;
881  char line[LINE_SIZE];
882  int li,time_steps,j,id,ph;
883  double value;
884  int sbi;
885  int n_subst_ = problem->transport->n_subst_;
886  bool found;
887 
888 
889  nodes = (struct Node**)xmalloc((4)* sizeof(struct Node*));
890  elements = (struct TElement**)xmalloc((4)* sizeof(struct TElement*));
891 
892 
893  in = open_temp_files(problem->transport, "%s.tmp", "rt" );
894  out = open_temp_files(problem->transport, "%s", (problem->pos_format_id == GMSH_MSH_ASCII) ? "wt" : "wb" );
895 
896 
897  // find number of time steps
898  found=skip_to(in[0], "VALUES");
899  ASSERT(found , "Can not find section: VALUES\n");
900  time_steps = 0;
901  while(skip_to(in[0], "Time"))
902  time_steps++;
903 
904  for(ph = 0; ph < 4; ph++){ // phase for cycle
905 
906  if (in[ph] != NULL){
907  nodes[ph] = (struct Node *)xmalloc((mesh->max_nod_id + 1)* sizeof(struct Node));
908  elements[ph] = (struct TElement *)xmalloc((mesh->max_elm_id + 1)* sizeof(struct TElement));
909  }
910  else
911  continue;
912 
913  // allocate and initialize tmp nodes conc
914  for (li = 0; li <= mesh->max_nod_id; li++){
915  nodes[ph][li].conc = (double **)xmalloc(n_subst_ * sizeof(double*));
916  for( sbi = 0; sbi < n_subst_; sbi++ )
917  nodes[ph][li].conc[sbi] = (double *)xmalloc(time_steps* sizeof(double));
918  for( sbi = 0; sbi < n_subst_; sbi++ )
919  for (j = 0; j < time_steps; j++)
920  nodes[ph][li].conc[sbi][j] = 0;
921  }
922  // allocate and initialize tmp element conc
923  for (li = 0; li <= mesh->max_elm_id; li++){
924  elements[ph][li].conc = (double **)xmalloc(n_subst_ * sizeof(double*));
925  for( sbi = 0; sbi < n_subst_; sbi++ )
926  elements[ph][li].conc[sbi] = (double *)xmalloc(time_steps* sizeof(double));
927  for( sbi = 0; sbi < n_subst_; sbi++ )
928  for (j = 0; j < time_steps; j++)
929  elements[ph][li].conc[sbi][j] = 0;
930  }
931 
932  // now read values for nodes
933  xfclose(in[ph]);
934  in = open_temp_files(problem->transport, "%s.tmp", "rt" );
935  found=skip_to(in[ph], "VALUES");
936  ASSERT(found, "Can not find section: VALUES\n");
937 
938 
939  j = 0;
940  while(skip_to(in[ph], "Nodes")){
941  for (li = 0; li < mesh->n_nodes; li++){
942  xfgets( line, LINE_SIZE - 2, in[ph] );
943  id = atoi( xstrtok( line) );
944  for( sbi = 0; sbi < n_subst_; sbi++ ) {
945  value = atof ( xstrtok( NULL) );
946  nodes[ph][id].conc[sbi][j] = value;
947  }
948  }
949  j++;
950  }
951 
952 
953  // end of TMP file reading
954  // read values for elements
955  xfclose(in[ph]);
956  in = open_temp_files(problem->transport, "%s.tmp", "rt" );
957  found=skip_to(in[ph], "VALUES");
958  ASSERT(found, "Can not find section: VALUES\n");
959 
960  j = 0;
961  while(skip_to(in[ph], "Elements")){
962  for (li = 0; li < mesh->n_elements(); li++){
963  xfgets( line, LINE_SIZE - 2, in[ph] );
964  id = atoi( xstrtok( line) );
965  for( sbi = 0; sbi < n_subst_; sbi++ ) {
966  value = atof ( xstrtok( NULL) );
967  elements[ph][id].conc[sbi][j] = value;
968  }
969  }
970  j++;
971  }
972  xfclose(in[ph]);
973 
974  xprintf( Msg, "Writing transport output files... ")/*orig verb 2*/;
975  switch(problem->pos_format_id){
976  case GMSH_MSH_ASCII:
977  write_transport_ascii_data(out[ph],problem,nodes,elements,time_steps,ph);
978  break;
979  case GMSH_MSH_BIN:
980  write_transport_binary_data(out[ph],problem,nodes,elements,time_steps,ph);
981  break;
982  }
983  xfclose(out[ph]);
984  }// end ph
985  xprintf( Msg, "O.K.\n")/*orig verb 2*/;
986 
987  xfree(nodes);
988  xfree(elements);
989  xfree(in);
990  xfree(out);
991 
992 }
993 #endif
994 
Generated on Thu May 29 2014 23:14:47 for Flow123d by   doxygen 1.8.4