Flow123d
ppfcs.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  * @ingroup io
27  * @brief Cross-sections computation
28  *
29  */
30 
31 #include "transport/transport.h"
32 
33 #include "system/system.hh"
34 #include "xio.h"
35 #include "system/math_fce.h"
36 //#include "problem.h"
37 #include "mesh/mesh.h"
38 #include "ppfcs.h"
39 #include "io/read_ini.h"
40 #include "materials.hh"
41 
42 
43 static void flow_cs(struct Transport *transport);
44 static void output_AGE(struct Transport *transport,double time);
45 static double particle_test(struct Transport *transport);
46 static void clear_tbc(struct Transport *transport);
47 
48 static int create_flow_section(struct Transport *transport);
49 static void sfsec(struct Transport *transport);
50 static double node_pos(Node* node, double *eqn);
51 static void compute_cutplane(struct Transport *trannsport);
52 static double *cut_point(Node* node0, Node* node1, double eqn[4]);
53 static double *assign_node_coordinates(Node* node);
54 static struct ElementCut *new_ec(struct Transport *transport,ElementIter elm,struct ElementCut *prev_ec);
55 static void compute_flux(struct ElementCut *ec,struct FSection *fs);
56 static void output_FCS(struct Transport *transport);
57 
58 //==============================================================================
59 // Flow crosssection
60 //==============================================================================
61 void flow_cs(struct Transport *transport)
62 {
63  struct ElementCut *ec;
64 
65  if (create_flow_section(transport)) {
66  sfsec(transport);
67  if(transport->fsec->elc != NULL){
68  compute_cutplane(transport);
69  FOR_ELEMENTCUT(ec)
70  compute_flux(ec,transport->fsec);
71  output_FCS(transport);
72  }
73  }
74 }
75 //==============================================================================
76 // Create flow crosssection
77 //==============================================================================
78 int create_flow_section(struct Transport *transport)
79 {
80 
81  // This is ancient initialization from problem.c
82  // TODO: Proper implementation of cross section
83 #if 0
84  problem->ftrans_out = get_b( "Output", "Write_ftrans_out", false );
85  problem->cross_section = get_b( "Output", "Cross_section", false ); //jh
86  problem->cs_params = get_s( "Output", "Cs_params", "0 0 0 0 0 0 0" ); //jh
87 // problem->res_run = get_b( "Output", "Cs_results_run", false ); //jh
88 // problem->res_fin = get_b( "Output", "Cs_results_final", false );
89  problem->specify_elm_output = get_b( "Output", "Specify_elm_type", false ); //jh temp
90  problem->output_elm_type = get_i( "Output", "Output_elm_type", 1 ); //jh temp
91  problem->fsec_params = get_s( "Output", "FCs_params", "0 0 0 0 0" );
92 // problem->CF_params = get_s( "Output", "ConfFlow_params", "0");
93 #endif
94 
95 
96  double norm;
97  int i;
98 
99  struct FSection *fsec;//=problem->fsec;
100 
101  transport->fsec =(struct FSection*)xmalloc(sizeof(struct FSection));
102  fsec = transport->fsec;
103  fsec->n_elm = 0;
104  fsec->elc = NULL;
105 
106  fsec->fcs_params = OptGetStr( "Output", "FCs_params", "0 0 0 0 0");
107 
108  if( sscanf( transport->fsec->fcs_params, "%lf %lf %lf %lf %d",
109  &fsec->eqn[0],
110  &fsec->eqn[1],
111  &fsec->eqn[2],
112  &fsec->eqn[3],
113  &fsec->axis_output) == 5 )
114  xprintf( Msg, "O.K.\n")/*orig verb 2*/; //118
115  else
116  xprintf(UsrErr,"Invalid flow cross section params.\n");
117  norm = vector_length(transport->fsec->eqn);
118 
119  if(norm > ZERO)
120  for(i=0;i<4;i++)
121  fsec->eqn[i] /= norm;
122  else{
123  xfree(fsec);
124  return false;
125  }
126  return true;
127 }
128 //==============================================================================
129 // Select elements of flow crosssection
130 //==============================================================================
131 void sfsec(struct Transport *transport)
132 {
133  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
134 
135  ElementIter elm;
136  NodeIter nod;
137  struct FSection *fs;
138  struct ElementCut *ec;
139  int i,j;
140  int b1,b2;
141 
142  fs = transport->fsec;
143  ec = NULL;
144 
145  FOR_NODES( nod ){
146  nod->faux = node_pos(nod,fs->eqn);
147  nod->aux = SGN(nod->faux);
148  }
149 
150  FOR_ELEMENTS(elm){
151  b1 = b2 = 0;
152  elm->aux = 0;
153  FOR_ELEMENT_NODES(elm,i){ //element cut
154  if(elm->node[i]->aux == -1){
155  b1 = 1;
156  elm->aux = -1;
157  }
158  if(elm->node[i]->aux == 1)
159  b2 = 1;
160  }
161  if((b1 * b2) != 0){
162  ec = new_ec(transport,elm,ec);
163  ec->type = 2;
164  }
165  else
166  FOR_ELEMENT_SIDES(elm,i){ //side cut
167  b1 = 1;
168  FOR_SIDE_NODES(elm->side[i],j)
169  if(elm->side[i]->node[j]->aux != 0)
170  b1 = 0;
171  if((b1 == 1) && (elm->aux == -1)){
172  ec = new_ec(transport,elm,ec);
173  ec->type = 1;
174  ec->sid = i;
175  }
176  }
177  }
178 }
179 //==============================================================================
180 // Create new element cut
181 //==============================================================================
182 struct ElementCut *new_ec(struct Transport *transport,ElementIter elm,struct ElementCut *prev_ec)
183 {
184  struct ElementCut *ec;
185  ec = (struct ElementCut*)xmalloc(sizeof(struct ElementCut));
186 
187  ec->n_point = 0;
188  ec->point[0] = NULL;
189  ec->point[1] = NULL;
190  ec->point[2] = NULL;
191  ec->point[3] = NULL;
192 
193  if(transport->fsec->elc == NULL)
194  transport->fsec->elc = ec;
195  else
196  prev_ec->next = ec;
197  ec->prev = prev_ec;
198  ec->next = NULL;
199  transport->fsec->n_elm++;
200  ec->element = elm;
201  return ec;
202 }
203 //==============================================================================
204 // Node position
205 //==============================================================================
206 double node_pos(Node* node, double* eqn)
207 {
208  double val;
209 
210  val = eqn[0] * node->getX() +
211  eqn[1] * node->getY() +
212  eqn[2] * node->getZ() +
213  eqn[3];
214 
215  return val;
216 }
217 //==============================================================================
218 // Compute breakplane
219 //==============================================================================
220 void compute_cutplane(struct Transport *transport)
221 {
222  ElementIter elm;
223  struct ElementCut *ec;
224  int i,j,n;
225 
226  FOR_ELEMENTCUT(ec){
227  elm = ec->element;
228  n = 0;
229  FOR_ELEMENT_NODES(elm,i)
230  if(elm->node[i]->aux == 0)
231  ec->point[ec->n_point++] = assign_node_coordinates(elm->node[i]);
232  if(ec->type != 1)
233  FOR_ELEMENT_NODES(elm,i){
234  n++;
235  for(j=n;j<elm->n_nodes;j++)
236  if((elm->node[i]->aux != 0) && (elm->node[j]->aux != 0)){
237  ec->point[ec->n_point] =
238  cut_point(elm->node[i],elm->node[j],transport->fsec->eqn);
239  if(ec->point[ec->n_point] != NULL)
240  ec->n_point++;
241  }
242  }
243  }
244 
245  /* FOR_ELEMENTPLANE(ep){
246  printf("\n\[id:%d\tdim:%d\tn_point:%d\ttype:%d\]",ep->element->id,ep->element->dim,ep->n_point,ep->type);
247  for(j=0;j<ep->n_point;j++){
248  printf("\n");
249  for(k=0;k<3;k++)
250  printf("\t%f",ep->point[j][k]);
251  }
252  printf("\nBalance:%f\n",ep->element->balance);
253  }
254  printf("\nn_elm:%d",problem->fsec->n_elm);
255  printf("\n%s",problem->fsec_params);
256  getchar();
257 
258 */
259 }
260 //==============================================================================
261 // Compute break point
262 //==============================================================================
263 double *cut_point(Node* node0, Node* node1, double eqn[4])
264 {
265  double v[3],t,d,*point;
266 
267  point = NULL;
268 
269  v[0] = node0->getX() - node1->getX();
270  v[1] = node0->getY() - node1->getY();
271  v[2] = node0->getZ() - node1->getZ();
272  t = node1->faux;
273  d = scalar_product(eqn,v);
274  if (fabs(d) < ZERO)
275  return point;
276  else
277  t /= -d;
278 
279  if((t <= 1) && (t > ZERO)){
280  point = (double*)xmalloc(3*sizeof(double));
281  point[0] = node1->getX() + t * v[0];
282  point[1] = node1->getY() + t * v[1];
283  point[2] = node1->getZ() + t * v[2];
284  }
285  return point;
286 }
287 //==============================================================================
288 // Assign node coordinates
289 //==============================================================================
290 double *assign_node_coordinates(Node* node)
291 {
292  double *point;
293 
294  point = (double*)xmalloc(3*sizeof(double));
295  point[0] = node->getX();
296  point[1] = node->getY();
297  point[2] = node->getZ();
298 
299  return point;
300 }
301 //==============================================================================
302 // Compute flux
303 //==============================================================================
304 void compute_flux(struct ElementCut *ec,struct FSection *fs)
305 {
306  ElementIter elm;
307  double vector[3],u[3],v[3],en[3],flux,l;
308  int i,sg;
309 
310 
311  for(i=0;i<3;i++)
312  vector[i] = 0.0;
313  elm = ec->element;
314 
315  if(ec->type != 1)
316  switch(elm->dim){
317  case 1:
318  /*
319  l = line_length(elm->node[1]->getX(),elm->node[1]->getY(),elm->node[1]->getZ(),elm->node[0]->getX(),elm->node[0]->getY(),elm->node[0]->getZ());
320 
321  if(elm->side[0]->aux == 1){
322  ld = line_length(ep->point[0][0],ep->point[0][1],ep->point[0][2],elm->node[0]->getX(),elm->node[0]->getY(),elm->node[0]->getZ());
323  flux = elm->side[0]->flux + ld / l * (-elm->side[1]->flux - elm->side[0]->flux);
324  }
325  else{
326  ld = line_length(ep->point[0][0],ep->point[0][1],ep->point[0][2],elm->node[1]->getX(),elm->node[1]->getY(),elm->node[1]->getZ());
327  flux = elm->side[1]->flux + ld / l * (-elm->side[0]->flux - elm->side[1]->flux);
328  }
329  */
330  flux = SGN(scalar_product(fs->eqn,elm->vector)) * elm->side[0]->metrics * elm->v_length;
331  ec->cutflux = flux;
332  break;
333  case 2:
334  u[ 0 ] = elm->node[ 1 ]->getX() - elm->node[ 0 ]->getX();
335  u[ 1 ] = elm->node[ 1 ]->getY() - elm->node[ 0 ]->getY();
336  u[ 2 ] = elm->node[ 1 ]->getZ() - elm->node[ 0 ]->getZ();
337  v[ 0 ] = elm->node[ 2 ]->getX() - elm->node[ 0 ]->getX();
338  v[ 1 ] = elm->node[ 2 ]->getY() - elm->node[ 0 ]->getY();
339  v[ 2 ] = elm->node[ 2 ]->getZ() - elm->node[ 0 ]->getZ();
340  vector_product( u, v, en );
341  normalize_vector( en );
342  vector_difference(ec->point[0],ec->point[1],v);
343  vector_product( v, en, u );
344  normalize_vector(u);
345  l = vector_length(v);
346  sg = SGN(scalar_product(elm->vector,fs->eqn));
347  ec->cutflux = fabs(scalar_product(u,elm->vector)) * elm->material->size * l * sg;
348  //ec->cutflux *= elm->material->size * l * sg;
349  break;
350  case 3:
351  break;
352  }
353  else
354  ec->cutflux = elm->side[ec->sid]->flux;
355 }
356 
357 //==============================================================================
358 // OUTPUT FLOW CROSSSECTION
359 //==============================================================================
360 void output_FCS(struct Transport *transport)
361 {
362  FILE *out;
363  struct ElementCut *ec;
364  double flux = 0;
365 
366  char dbl_fmt[ 16 ],file[LINE_SIZE];
367 
368  sprintf( dbl_fmt, "%%.%dg", ConstantDB::getInstance()->getInt("Out_digit"));
369  sprintf(file,"%s.fcs", IONameHandler::get_instance()->get_output_file_name(OptGetFileName("Output", "Output_file", NULL)) );
370  out = xfopen(file,"wt");
371 
372  FOR_ELEMENTCUT(ec){
373  flux += ec->cutflux;
374  switch(ec->element->dim){
375  case 1:
376  // xfprintf(out,"%d\t",ec->element->id);
377  xfprintf(out,dbl_fmt,ec->point[0][transport->fsec->axis_output]);
378  xfprintf(out,"\t");
379  xfprintf(out,dbl_fmt,ec->cutflux);
380  break;
381  case 2:
382  if(ec->point[0][transport->fsec->axis_output] > ec->point[1][transport->fsec->axis_output]){
383  xfprintf(out,dbl_fmt,ec->point[0][transport->fsec->axis_output]);
384  xfprintf(out,":");
385  xfprintf(out,dbl_fmt,ec->point[1][transport->fsec->axis_output]);
386  }
387  else{
388  xfprintf(out,dbl_fmt,ec->point[1][transport->fsec->axis_output]);
389  xfprintf(out,":");
390  xfprintf(out,dbl_fmt,ec->point[0][transport->fsec->axis_output]);
391  }
392  xfprintf(out,"\t");
393  xfprintf(out,dbl_fmt,ec->cutflux);
394  break;
395  }
396  xfprintf(out,"\n");
397  }
398 
399  xfprintf(out,"Total flux of cut : ");
400  xfprintf(out,dbl_fmt,flux);
401  // printf("\nTotal flux of cut: %f\n",flux);
402  xfprintf(out,"\n");
403  xfclose( out );
404 
405 }
406 //==============================================================================
407 
408 // DECOVALEX
409 
410 //==============================================================================
411 /*
412 //==============================================================================
413 // OUTPUT AGE
414 //==============================================================================
415 void output_AGE(struct Transport *transport,double time)
416 {
417  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
418 
419  FILE *out;
420  struct ElementCut *ec;
421  double Np = 0.0;
422  double pc = 0.0;
423  double flux;
424  char dbl_fmt[ 16 ],file[LINE_SIZE];
425  int i,id;
426 
427  sprintf( dbl_fmt, "%%.%dg", ConstantDB::getInstance()->getInt("Out_digit"));
428  sprintf(file,"%s.age", OptGetFileName("Output", "Output_file", NULL) );
429  out = xfopen(file,(time==0.0) ? "wt":"at");
430 
431  // TODO: remove epos_id and find_id
432  for(i=0;i<mesh->n_elements();i++)
433  pc += transport->conc[MOBILE][0][i] * mesh->element.find_id(mesh->epos_id[i])->volume;
434 
435  FOR_ELEMENTCUT(ec){
436  id = id2pos(mesh,ELEMENT_FULL_ITER(ec->element).id(),mesh->epos_id,ELM);
437  flux = ec->cutflux;
438  Np += transport->conc[MOBILE][0][id] * flux * transport->time_step;
439  }
440 
441 
442  if((Np > 0.0) || (time == 0.0)){
443  xfprintf(out,dbl_fmt,time);
444  xfprintf(out,"\t");
445  xfprintf(out,dbl_fmt,Np);
446  xfprintf(out,"\t");
447  xfprintf(out,dbl_fmt,pc);
448  xfprintf(out,"\n");
449  }
450 
451  xfclose( out );
452  // getchar();
453 } */
454 //==============================================================================
455 // PARTICLE TEST
456 //==============================================================================
457 /*
458 double particle_test(struct Transport *transport)
459 {
460  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
461 
462  int i;
463  double pc = 0.0;
464 
465  for(i=0;i < mesh->n_elements();i++){
466  pc += transport->conc[MOBILE][0][i] * mesh->element.find_id(mesh->epos_id[i])->volume;
467  }
468  return pc;
469 } */
470 //==============================================================================
471 // CLEAR TRANSPORT BC
472 //==============================================================================
473 void clear_tbc(struct Transport *transport)
474 {
475  Mesh* mesh = (Mesh*) ConstantDB::getInstance()->getObject(MESH::MAIN_INSTANCE);
476 
477  int start,stop,i;
478  start = mesh->n_elements();
479  stop = mesh->n_elements() + mesh->n_boundaries(); //-1
480  for(i=start;i<stop;i++){
481  transport->conc[MOBILE][0][i] = 0.0;
482  transport->pconc[MOBILE][0][i] = 0.0;
483  }
484 }
485 //------------------------------------------------------------------------------
Generated on Thu May 29 2014 23:14:47 for Flow123d by   doxygen 1.8.4