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