navier_stokes_surface_power_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for specific surface elements
31 
32 #ifndef OOMPH_NAVIER_STOKES_SURFACE_POWER_ELEMENTS_HEADER
33 #define OOMPH_NAVIER_STOKES_SURFACE_POWER_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 
41 //OOMPH-LIB headers
42 #include "../generic/Qelements.h"
43 
44 namespace oomph
45 {
46 
47 //======================================================================
48 /// A class of elements that allow the determination of the power
49 /// input and various other fluxes over the domain boundaries.
50 /// The element operates as a FaceElement and attaches itself
51 /// to a bulk element of the type specified by the template
52 /// argument.
53 //======================================================================
54 template <class ELEMENT>
55 class NavierStokesSurfacePowerElement : public virtual FaceGeometry<ELEMENT>,
56 public virtual FaceElement
57 {
58 
59 public:
60 
61  ///Constructor, which takes a "bulk" element and the value of the index
62  ///and its limit
64  const int &face_index) :
65  FaceGeometry<ELEMENT>(), FaceElement()
66  {
67  //Attach the geometrical information to the element. N.B. This function
68  //also assigns nbulk_value from the required_nvalue of the bulk element
69  element_pt->build_face_element(face_index,this);
70 
71  //Set the dimension from the dimension of the first node
72  Dim = node_pt(0)->ndim();
73  }
74 
75 
76  /// \short The "global" intrinsic coordinate of the element when
77  /// viewed as part of a geometric object should be given by
78  /// the FaceElement representation, by default
79  /// This final over-ride is required for cases where the
80  /// FaceElement is a SolidFiniteElement because both SolidFiniteElements
81  /// and FaceElements overload zeta_nodal.
82  double zeta_nodal(const unsigned &n, const unsigned &k,
83  const unsigned &i) const
84  {return FaceElement::zeta_nodal(n,k,i);}
85 
86 
87 
88  /// \short Get drag force (traction acting on fluid)
90  {
91  std::ofstream dummy_file;
92  return drag_force(dummy_file);
93  }
94 
95 
96 
97  /// \short Get drag force (traction acting on fluid)
98  /// Doc in outfile.
99  Vector<double> drag_force(std::ofstream& outfile)
100  {
101 
102  // Spatial dimension of element
103  unsigned ndim=dim();
104 
105  // Initialise
106  Vector<double> drag(ndim+1,0.0);
107 
108  //Vector of local coordinates in face element
109  Vector<double> s(ndim);
110 
111  // Vector for global Eulerian coordinates
112  Vector<double> x(ndim+1);
113 
114  // Vector for local coordinates in bulk element
115  Vector<double> s_bulk(ndim+1);
116 
117  //Set the value of n_intpt
118  unsigned n_intpt = integral_pt()->nweight();
119 
120 
121  // Get pointer to assocated bulk element
122  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
123 
124  // Hacky: This is only appropriate for 3 point integration of
125  // 1D line elements
126  if (outfile.is_open()) outfile << "ZONE I=3" << std::endl;
127 
128  // Loop over the integration points
129  for (unsigned ipt=0;ipt<n_intpt;ipt++)
130  {
131 
132  //Assign values of s in FaceElement and local coordinates in bulk element
133  for(unsigned i=0;i<ndim;i++)
134  {
135  s[i] = integral_pt()->knot(ipt,i);
136  }
137 
138  //Get the bulk coordinates
139  this->get_local_coordinate_in_bulk(s,s_bulk);
140 
141  //Get the integral weight
142  double w = integral_pt()->weight(ipt);
143 
144  // Get jacobian of mapping
145  double J=J_eulerian(s);
146 
147  //Premultiply the weights and the Jacobian
148  double W = w*J;
149 
150  // Get x position as Vector
151  interpolated_x(s,x);
152 
153 #ifdef PARANOID
154 
155  // Get x position as Vector from bulk element
156  Vector<double> x_bulk(ndim+1);
157  bulk_el_pt->interpolated_x(s_bulk,x_bulk);
158 
159  double max_legal_error=1.0e-5;
160  double error=0.0;
161  for(unsigned i=0;i<ndim+1;i++)
162  {
163  error+=fabs(x[i]- x_bulk[i]);
164  }
165  if (error>max_legal_error)
166  {
167  std::ostringstream error_stream;
168  error_stream << "difference in Eulerian posn from bulk and face: "
169  << error << " exceeds threshold " << max_legal_error
170  << std::endl;
171  throw OomphLibError(
172  error_stream.str(),
173  OOMPH_CURRENT_FUNCTION,
174  OOMPH_EXCEPTION_LOCATION);
175  }
176 #endif
177 
178  // Outer unit normal
179  Vector<double> normal(ndim+1);
180  outer_unit_normal(s,normal);
181 
182  // Get velocity from bulk element
183  Vector<double> veloc(ndim+1);
184  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
185 
186  // Get traction from bulk element
187  Vector<double> traction(ndim+1);
188  bulk_el_pt->get_traction(s_bulk,normal,traction);
189 
190  // Integrate
191  for (unsigned i=0;i<ndim+1;i++)
192  {
193  drag[i]+=traction[i]*W;
194  }
195 
196  if (outfile.is_open())
197  {
198  //Output x,y,...,
199  for(unsigned i=0;i<ndim+1;i++)
200  {
201  outfile << x[i] << " ";
202  }
203 
204  //Output traction
205  for(unsigned i=0;i<ndim+1;i++)
206  {
207  outfile << traction[i] << " ";
208  }
209 
210 
211  //Output normal
212  for(unsigned i=0;i<ndim+1;i++)
213  {
214  outfile << normal[i] << " ";
215  }
216 
217  outfile << std::endl;
218  }
219  }
220  return drag;
221  }
222 
223  /// \short Get integral of instantaneous rate of work done by
224  /// the traction that's exerted onto the fluid.
226  {
227  std::ofstream dummy_file;
228  return get_rate_of_traction_work(dummy_file);
229  }
230 
231 
232 
233  /// \short Get integral of instantaneous rate of work done by
234  /// the traction that's exerted onto the fluid. Doc in outfile.
235  double get_rate_of_traction_work(std::ofstream& outfile)
236  {
237 
238  // Initialise
239  double rate_of_work_integral=0.0;
240 
241  // Spatial dimension of element
242  unsigned ndim=dim();
243 
244  //Vector of local coordinates in face element
245  Vector<double> s(ndim);
246 
247  // Vector for global Eulerian coordinates
248  Vector<double> x(ndim+1);
249 
250  // Vector for local coordinates in bulk element
251  Vector<double> s_bulk(ndim+1);
252 
253  //Set the value of n_intpt
254  unsigned n_intpt = integral_pt()->nweight();
255 
256 
257  // Get pointer to assocated bulk element
258  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
259 
260  // Hacky: This is only appropriate for 3x3 integration of
261  // 2D quad elements
262  if (outfile.is_open()) outfile << "ZONE I=3, J=3" << std::endl;
263 
264  //Loop over the integration points
265  for (unsigned ipt=0;ipt<n_intpt;ipt++)
266  {
267 
268  //Assign values of s in FaceElement and local coordinates in bulk element
269  for(unsigned i=0;i<ndim;i++)
270  {
271  s[i] = integral_pt()->knot(ipt,i);
272  }
273 
274  //Get the bulk coordinates
275  this->get_local_coordinate_in_bulk(s,s_bulk);
276 
277  //Get the integral weight
278  double w = integral_pt()->weight(ipt);
279 
280  // Get jacobian of mapping
281  double J=J_eulerian(s);
282 
283  //Premultiply the weights and the Jacobian
284  double W = w*J;
285 
286  // Get x position as Vector
287  interpolated_x(s,x);
288 
289 #ifdef PARANOID
290 
291  // Get x position as Vector from bulk element
292  Vector<double> x_bulk(ndim+1);
293  bulk_el_pt->interpolated_x(s_bulk,x_bulk);
294 
295  double max_legal_error=1.0e-5;
296  double error=0.0;
297  for(unsigned i=0;i<ndim+1;i++)
298  {
299  error+=fabs(x[i]- x_bulk[i]);
300  }
301  if (error>max_legal_error)
302  {
303  std::ostringstream error_stream;
304  error_stream << "difference in Eulerian posn from bulk and face: "
305  << error << " exceeds threshold " << max_legal_error
306  << std::endl;
307  throw OomphLibError(
308  error_stream.str(),
309  OOMPH_CURRENT_FUNCTION,
310  OOMPH_EXCEPTION_LOCATION);
311  }
312 #endif
313 
314  // Outer unit normal
315  Vector<double> normal(ndim+1);
316  outer_unit_normal(s,normal);
317 
318 
319  // Get velocity from bulk element
320  Vector<double> veloc(ndim+1);
321  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
322 
323  // Get traction from bulk element
324  Vector<double> traction(ndim+1);
325  bulk_el_pt->get_traction(s_bulk,normal,traction);
326 
327 
328  // Local rate of work:
329  double rate_of_work=0.0;
330  for (unsigned i=0;i<ndim+1;i++)
331  {
332  rate_of_work+=traction[i]*veloc[i];
333  }
334 
335  // Add rate of work
336  rate_of_work_integral+=rate_of_work*W;
337 
338  if (outfile.is_open())
339  {
340  //Output x,y,...,
341  for(unsigned i=0;i<ndim+1;i++)
342  {
343  outfile << x[i] << " ";
344  }
345 
346  //Output traction
347  for(unsigned i=0;i<ndim+1;i++)
348  {
349  outfile << traction[i] << " ";
350  }
351 
352  //Output veloc
353  for(unsigned i=0;i<ndim+1;i++)
354  {
355  outfile << veloc[i] << " ";
356  }
357 
358  //Output normal
359  for(unsigned i=0;i<ndim+1;i++)
360  {
361  outfile << normal[i] << " ";
362  }
363 
364  //Output local rate of work
365  for(unsigned i=0;i<ndim+1;i++)
366  {
367  outfile << rate_of_work << " ";
368  }
369 
370  outfile << std::endl;
371  }
372  }
373 
374  return rate_of_work_integral;
375 
376  }
377 
378 
379 
380  /// \short Get integral of instantaneous rate of work done by
381  /// the traction that's exerted onto the fluid, decomposed into pressure
382  /// and normal and tangential viscous components.
383  void get_rate_of_traction_work_components(double& rate_of_work_integral_p,
384  double& rate_of_work_integral_n,
385  double& rate_of_work_integral_t)
386  {
387  std::ofstream dummy_file;
389  rate_of_work_integral_p,
390  rate_of_work_integral_n,
391  rate_of_work_integral_t);
392  }
393 
394 
395 
396  /// \short Get integral of instantaneous rate of work done by
397  /// the traction that's exerted onto the fluid, decomposed into pressure
398  /// and normal and tangential viscous components. Doc in outfile.
399  void get_rate_of_traction_work_components(std::ofstream& outfile,
400  double& rate_of_work_integral_p,
401  double& rate_of_work_integral_n,
402  double& rate_of_work_integral_t)
403  {
404 
405  // Initialise
406  rate_of_work_integral_p=0;
407  rate_of_work_integral_n=0;
408  rate_of_work_integral_t=0;
409 
410  // Spatial dimension of element
411  unsigned ndim=dim();
412 
413  //Vector of local coordinates in face element
414  Vector<double> s(ndim);
415 
416  // Vector for global Eulerian coordinates
417  Vector<double> x(ndim+1);
418 
419  // Vector for local coordinates in bulk element
420  Vector<double> s_bulk(ndim+1);
421 
422  //Set the value of n_intpt
423  unsigned n_intpt = integral_pt()->nweight();
424 
425 
426  // Get pointer to assocated bulk element
427  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
428 
429  // Hacky: This is only appropriate for 3x3 integration of
430  // 2D quad elements
431  if (outfile.is_open()) outfile << "ZONE I=3, J=3" << std::endl;
432 
433  //Loop over the integration points
434  for (unsigned ipt=0;ipt<n_intpt;ipt++)
435  {
436 
437  //Assign values of s in FaceElement and local coordinates in bulk element
438  for(unsigned i=0;i<ndim;i++)
439  {
440  s[i] = integral_pt()->knot(ipt,i);
441  }
442 
443  //Get the bulk coordinates
444  this->get_local_coordinate_in_bulk(s,s_bulk);
445 
446  //Get the integral weight
447  double w = integral_pt()->weight(ipt);
448 
449  // Get jacobian of mapping
450  double J=J_eulerian(s);
451 
452  //Premultiply the weights and the Jacobian
453  double W = w*J;
454 
455  // Get x position as Vector
456  interpolated_x(s,x);
457 
458 #ifdef PARANOID
459 
460  // Get x position as Vector from bulk element
461  Vector<double> x_bulk(ndim+1);
462  bulk_el_pt->interpolated_x(s_bulk,x_bulk);
463 
464  double max_legal_error=1.0e-5;
465  double error=0.0;
466  for(unsigned i=0;i<ndim+1;i++)
467  {
468  error+=fabs(x[i]- x_bulk[i]);
469  }
470  if (error>max_legal_error)
471  {
472  std::ostringstream error_stream;
473  error_stream << "difference in Eulerian posn from bulk and face: "
474  << error << " exceeds threshold " << max_legal_error
475  << std::endl;
476  throw OomphLibError(
477  error_stream.str(),
478  OOMPH_CURRENT_FUNCTION,
479  OOMPH_EXCEPTION_LOCATION);
480  }
481 #endif
482 
483  // Outer unit normal
484  Vector<double> normal(ndim+1);
485  outer_unit_normal(s,normal);
486 
487 
488  // Get velocity from bulk element
489  Vector<double> veloc(ndim+1);
490  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
491 
492  // Get traction from bulk element
493  Vector<double> traction_p(ndim+1);
494  Vector<double> traction_n(ndim+1);
495  Vector<double> traction_t(ndim+1);
496  bulk_el_pt->get_traction(s_bulk,normal,traction_p,traction_n,traction_t);
497 
498 
499  // Local rate of work:
500  double rate_of_work_p=0.0;
501  double rate_of_work_n=0.0;
502  double rate_of_work_t=0.0;
503  for (unsigned i=0;i<ndim+1;i++)
504  {
505  rate_of_work_p+=traction_p[i]*veloc[i];
506  rate_of_work_n+=traction_n[i]*veloc[i];
507  rate_of_work_t+=traction_t[i]*veloc[i];
508  }
509 
510  // Add rate of work
511  rate_of_work_integral_p+=rate_of_work_p*W;
512  rate_of_work_integral_n+=rate_of_work_n*W;
513  rate_of_work_integral_t+=rate_of_work_t*W;
514 
515  if (outfile.is_open())
516  {
517  //Output x,y,...,
518  for(unsigned i=0;i<ndim+1;i++)
519  {
520  outfile << x[i] << " ";
521  }
522 
523  //Output traction due to pressure
524  for(unsigned i=0;i<ndim+1;i++)
525  {
526  outfile << traction_p[i] << " ";
527  }
528 
529  //Output traction due to viscous normal stress
530  for(unsigned i=0;i<ndim+1;i++)
531  {
532  outfile << traction_n[i] << " ";
533  }
534 
535  //Output traction due to viscous tangential stress
536  for(unsigned i=0;i<ndim+1;i++)
537  {
538  outfile << traction_t[i] << " ";
539  }
540 
541  //Output veloc
542  for(unsigned i=0;i<ndim+1;i++)
543  {
544  outfile << veloc[i] << " ";
545  }
546 
547  //Output normal
548  for(unsigned i=0;i<ndim+1;i++)
549  {
550  outfile << normal[i] << " ";
551  }
552 
553  //Output local rate of work due to pressure
554  for(unsigned i=0;i<ndim+1;i++)
555  {
556  outfile << rate_of_work_p << " ";
557  }
558 
559  //Output local rate of work due to viscous normal stress
560  for(unsigned i=0;i<ndim+1;i++)
561  {
562  outfile << rate_of_work_n << " ";
563  }
564 
565  //Output local rate of work due to viscous tangential stress
566  for(unsigned i=0;i<ndim+1;i++)
567  {
568  outfile << rate_of_work_t << " ";
569  }
570 
571  outfile << std::endl;
572  }
573  }
574 
575  }
576 
577 
578 
579 
580  /// \short Get integral of kinetic energy flux
582  {
583  std::ofstream dummy_file;
584  return get_kinetic_energy_flux(dummy_file);
585  }
586 
587 
588 
589  /// \short Get integral of kinetic energy flux and doc
590  double get_kinetic_energy_flux(std::ofstream& outfile)
591  {
592  // Initialise
593  double kinetic_energy_flux_integral=0.0;
594 
595  // Spatial dimension of element
596  unsigned ndim=dim();
597 
598  //Vector of local coordinates in face element
599  Vector<double> s(ndim);
600 
601  // Vector for global Eulerian coordinates
602  Vector<double> x(ndim+1);
603 
604  // Vector for local coordinates in bulk element
605  Vector<double> s_bulk(ndim+1);
606 
607  //Set the value of n_intpt
608  unsigned n_intpt = integral_pt()->nweight();
609 
610  // Get pointer to assocated bulk element
611  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
612 
613  // Hacky: This is only appropriate for 3x3 integration of
614  // 2D quad elements
615  if (outfile.is_open()) outfile << "ZONE I=3, J=3" << std::endl;
616 
617  //Loop over the integration points
618  for (unsigned ipt=0;ipt<n_intpt;ipt++)
619  {
620 
621  //Assign values of s in FaceElement and local coordinates in bulk element
622  for(unsigned i=0;i<ndim;i++)
623  {
624  s[i] = integral_pt()->knot(ipt,i);
625  }
626 
627  //Get the bulk coordinates
628  this->get_local_coordinate_in_bulk(s,s_bulk);
629 
630  //Get the integral weight
631  double w = integral_pt()->weight(ipt);
632 
633  // Get jacobian of mapping
634  double J=J_eulerian(s);
635 
636  //Premultiply the weights and the Jacobian
637  double W = w*J;
638 
639  // Get x position as Vector
640  interpolated_x(s,x);
641 
642 #ifdef PARANOID
643 
644  // Get x position as Vector from bulk element
645  Vector<double> x_bulk(ndim+1);
646  bulk_el_pt->interpolated_x(s_bulk,x_bulk);
647 
648  double max_legal_error=1.0e-5;
649  double error=0.0;
650  for(unsigned i=0;i<ndim+1;i++)
651  {
652  error+=fabs(x[i]- x_bulk[i]);
653  }
654  if (error>max_legal_error)
655  {
656  std::ostringstream error_stream;
657  error_stream << "difference in Eulerian posn from bulk and face: "
658  << error << " exceeds threshold " << max_legal_error
659  << std::endl;
660  throw OomphLibError(
661  error_stream.str(),
662  OOMPH_CURRENT_FUNCTION,
663  OOMPH_EXCEPTION_LOCATION);
664  }
665 #endif
666 
667  // Outer unit normal
668  Vector<double> normal(ndim+1);
669  outer_unit_normal(s,normal);
670 
671  // Get velocity from bulk element
672  Vector<double> veloc(ndim+1);
673  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
674 
675  double kin_energy=0.0;
676  for (unsigned i=0;i<ndim+1;i++)
677  {
678  kin_energy+=veloc[i]*veloc[i];
679  }
680  kin_energy*=0.5;
681 
682  // Kinetic energy flux
683  double kin_energy_flux=0.0;
684  for (unsigned i=0;i<ndim+1;i++)
685  {
686  kin_energy_flux+=kin_energy*normal[i]*veloc[i];
687  }
688 
689  // Add to integral
690  kinetic_energy_flux_integral+=kin_energy_flux*W;
691 
692  if (outfile.is_open())
693  {
694  //Output x,y,...,
695  for(unsigned i=0;i<ndim+1;i++)
696  {
697  outfile << x[i] << " ";
698  }
699 
700  //Output veloc
701  for(unsigned i=0;i<ndim+1;i++)
702  {
703  outfile << veloc[i] << " ";
704  }
705 
706  //Output normal
707  for(unsigned i=0;i<ndim+1;i++)
708  {
709  outfile << normal[i] << " ";
710  }
711 
712  //Output local kin energy flux
713  outfile << kin_energy_flux << " ";
714 
715  outfile << std::endl;
716  }
717 
718  }
719 
720  return kinetic_energy_flux_integral;
721 
722  }
723 
724 
725  /// \short Get integral of volume flux
727  {
728  std::ofstream dummy_file;
729  return get_volume_flux(dummy_file);
730  }
731 
732 
733  /// \short Get integral of volume flux and doc
734  double get_volume_flux(std::ofstream& outfile)
735  {
736  // Initialise
737  double volume_flux_integral=0.0;
738 
739  // Spatial dimension of element
740  unsigned ndim=dim();
741 
742  //Vector of local coordinates in face element
743  Vector<double> s(ndim);
744 
745  // Vector for global Eulerian coordinates
746  Vector<double> x(ndim+1);
747 
748  // Vector for local coordinates in bulk element
749  Vector<double> s_bulk(ndim+1);
750 
751  //Set the value of n_intpt
752  unsigned n_intpt = integral_pt()->nweight();
753 
754  // Get pointer to assocated bulk element
755  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
756 
757  // Hacky: This is only appropriate for 3x3 integration of
758  // 2D quad elements
759  if (outfile.is_open()) outfile << "ZONE I=3, J=3" << std::endl;
760 
761  //Loop over the integration points
762  for (unsigned ipt=0;ipt<n_intpt;ipt++)
763  {
764 
765  //Assign values of s in FaceElement and local coordinates in bulk element
766  for(unsigned i=0;i<ndim;i++)
767  {
768  s[i] = integral_pt()->knot(ipt,i);
769  }
770 
771 
772  //Get the bulk coordinates
773  this->get_local_coordinate_in_bulk(s,s_bulk);
774 
775  //Get the integral weight
776  double w = integral_pt()->weight(ipt);
777 
778  // Get jacobian of mapping
779  double J=J_eulerian(s);
780 
781  //Premultiply the weights and the Jacobian
782  double W = w*J;
783 
784  // Get x position as Vector
785  interpolated_x(s,x);
786 
787 #ifdef PARANOID
788 
789  // Get x position as Vector from bulk element
790  Vector<double> x_bulk(ndim+1);
791  bulk_el_pt->interpolated_x(s_bulk,x_bulk);
792 
793  double max_legal_error=1.0e-5;
794  double error=0.0;
795  for(unsigned i=0;i<ndim+1;i++)
796  {
797  error+=fabs(x[i]- x_bulk[i]);
798  }
799  if (error>max_legal_error)
800  {
801  std::ostringstream error_stream;
802  error_stream << "difference in Eulerian posn from bulk and face: "
803  << error << " exceeds threshold " << max_legal_error
804  << std::endl;
805  throw OomphLibError(
806  error_stream.str(),
807  OOMPH_CURRENT_FUNCTION,
808  OOMPH_EXCEPTION_LOCATION);
809  }
810 #endif
811 
812  // Outer unit normal
813  Vector<double> normal(ndim+1);
814  outer_unit_normal(s,normal);
815 
816  // Get velocity from bulk element
817  Vector<double> veloc(ndim+1);
818  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
819 
820  // Volume flux
821  double volume_flux=0.0;
822  for (unsigned i=0;i<ndim+1;i++)
823  {
824  volume_flux+=normal[i]*veloc[i];
825  }
826 
827  // Add to integral
828  volume_flux_integral+=volume_flux*W;
829 
830  if (outfile.is_open())
831  {
832  //Output x,y,...,
833  for(unsigned i=0;i<ndim+1;i++)
834  {
835  outfile << x[i] << " ";
836  }
837 
838  //Output veloc
839  for(unsigned i=0;i<ndim+1;i++)
840  {
841  outfile << veloc[i] << " ";
842  }
843 
844  //Output normal
845  for(unsigned i=0;i<ndim+1;i++)
846  {
847  outfile << normal[i] << " ";
848  }
849 
850  //Output local volume flux
851  outfile << volume_flux << " ";
852 
853  outfile << std::endl;
854  }
855 
856  }
857 
858  return volume_flux_integral;
859 
860  }
861 
862 
863 private:
864 
865  ///The highest dimension of the problem
866  unsigned Dim;
867 
868 
869 };
870 
871 
872 }
873 
874 #endif
unsigned ndim() const
Access function to # of Eulerian coordinates.
Definition: geom_objects.h:185
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
A general Finite Element class.
Definition: elements.h:1271
Vector< double > drag_force(std::ofstream &outfile)
Get drag force (traction acting on fluid) Doc in outfile.
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:5695
double get_volume_flux(std::ofstream &outfile)
Get integral of volume flux and doc.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
double get_kinetic_energy_flux(std::ofstream &outfile)
Get integral of kinetic energy flux and doc.
double get_rate_of_traction_work()
Get integral of instantaneous rate of work done by the traction that's exerted onto the fluid...
double get_volume_flux()
Get integral of volume flux.
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
NavierStokesSurfacePowerElement(FiniteElement *const &element_pt, const int &face_index)
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
Vector< double > drag_force()
Get drag force (traction acting on fluid)
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Calculate the vector of local coordinate in the bulk element given the local coordinates in this Face...
Definition: elements.cc:6064
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:5033
unsigned Dim
The highest dimension of the problem.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
double get_rate_of_traction_work(std::ofstream &outfile)
Get integral of instantaneous rate of work done by the traction that's exerted onto the fluid...
void get_rate_of_traction_work_components(std::ofstream &outfile, double &rate_of_work_integral_p, double &rate_of_work_integral_n, double &rate_of_work_integral_t)
Get integral of instantaneous rate of work done by the traction that's exerted onto the fluid...
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
double get_kinetic_energy_flux()
Get integral of kinetic energy flux.
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
void get_rate_of_traction_work_components(double &rate_of_work_integral_p, double &rate_of_work_integral_n, double &rate_of_work_integral_t)
Get integral of instantaneous rate of work done by the traction that's exerted onto the fluid...