axisym_fluid_traction_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 elements that are used to integrate fluid tractions
27// This includes the guts (i.e. equations) because we want to inline them
28// for faster operation, although it slows down the compilation!
29#ifndef OOMPH_AXISYM_FLUID_TRACTION_ELEMENTS_HEADER
30#define OOMPH_AXISYM_FLUID_TRACTION_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37
38// OOMPH-LIB headers
39#include "../generic/shape.h"
40#include "../generic/elements.h"
41#include "../generic/element_with_external_element.h"
42
43
44namespace oomph
45{
46 //=======================================================================
47 /// Namespace containing the zero traction function for axisymmetric
48 /// Navier Stokes traction elements
49 //=======================================================================
50 namespace AxisymmetricNavierStokesTractionElementHelper
51 {
52 //=======================================================================
53 /// Default load function (zero traction)
54 //=======================================================================
55 extern void Zero_traction_fct(const double& time,
56 const Vector<double>& x,
57 const Vector<double>& N,
58 Vector<double>& load);
59
60 } // namespace AxisymmetricNavierStokesTractionElementHelper
61
62
63 //======================================================================
64 /// A class for elements that allow the imposition of an applied traction
65 /// in the axisym Navier Stokes eqns.
66 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
67 /// class and and thus, we can be generic enough without the need to have
68 /// a separate equations class.
69 //======================================================================
70 template<class ELEMENT>
72 : public virtual FaceGeometry<ELEMENT>,
73 public virtual FaceElement
74 {
75 protected:
76 /// Index at which the i-th velocity component is stored
78
79 /// Pointer to an imposed traction function. Arguments:
80 /// Eulerian coordinate; outer unit normal;
81 /// applied traction. (Not all of the input arguments will be
82 /// required for all specific load functions but the list should
83 /// cover all cases)
84 void (*Traction_fct_pt)(const double& time,
85 const Vector<double>& x,
86 const Vector<double>& n,
87 Vector<double>& result);
88
89
90 /// Get the traction vector: Pass number of integration point
91 /// (dummy), Eulerian coordinate and normal vector and return the load
92 /// vector (not all of the input arguments will be required for all specific
93 /// load functions but the list should cover all cases). This function is
94 /// virtual so it can be overloaded for FSI.
95 virtual void get_traction(const double& time,
96 const unsigned& intpt,
97 const Vector<double>& x,
98 const Vector<double>& n,
100 {
101 Traction_fct_pt(time, x, n, traction);
102 }
103
104
105 /// Helper function that actually calculates the residuals
106 // This small level of indirection is required to avoid calling
107 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
108 // which causes all kinds of pain if overloading later on
110 Vector<double>& residuals);
111
112
113 public:
114 /// Constructor, which takes a "bulk" element and the
115 /// value of the index and its limit
117 const int& face_index)
118 : FaceGeometry<ELEMENT>(), FaceElement()
119 {
120 // Attach the geometrical information to the element. N.B. This function
121 // also assigns nbulk_value from the required_nvalue of the bulk element
122 element_pt->build_face_element(face_index, this);
123
124 // Find the dimension of the problem
125 unsigned n_dim = element_pt->nodal_dimension();
126
127 // Find the index at which the velocity unknowns are stored
128 ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
129 this->U_index_axisymmetric_nst_traction.resize(n_dim + 1);
130 for (unsigned i = 0; i < n_dim + 1; i++)
131 {
132 this->U_index_axisymmetric_nst_traction[i] =
133 cast_element_pt->u_index_axi_nst(i);
134 }
135
136 // Zero traction
139 }
140
141
142 /// Reference to the traction function pointer
143 void (*&traction_fct_pt())(const double& time,
144 const Vector<double>& x,
145 const Vector<double>& n,
147 {
148 return Traction_fct_pt;
149 }
150
151
152 /// Return the residuals
154 {
156 }
157
158
159 /// Fill in contribution from Jacobian
161 DenseMatrix<double>& jacobian)
162 {
163 // Call the residuals
165 }
166
167 /// Specify the value of nodal zeta from the face geometry
168 /// The "global" intrinsic coordinate of the element when
169 /// viewed as part of a geometric object should be given by
170 /// the FaceElement representation, by default (needed to break
171 /// indeterminacy if bulk element is SolidElement)
172 double zeta_nodal(const unsigned& n,
173 const unsigned& k,
174 const unsigned& i) const
175 {
176 return FaceElement::zeta_nodal(n, k, i);
177 }
178
179 /// Output function
180 void output(std::ostream& outfile)
181 {
182 unsigned nplot = 5;
183 output(outfile, nplot);
184 }
185
186 /// Number of scalars/fields output by this element. Reimplements
187 /// broken virtual function in base class.
188 unsigned nscalar_paraview() const
189 {
190 // Number of dimensions
191 unsigned n_dim = this->nodal_dimension();
192
193 return 2 * (n_dim + 1);
194 }
195
196 /// Write values of the k-th scalar field at the plot points. Needs
197 /// to be implemented for each new specific element type.
198 void scalar_value_paraview(std::ofstream& file_out,
199 const unsigned& k,
200 const unsigned& nplot) const
201 {
202 // Number of dimensions
203 unsigned n_dim = this->nodal_dimension();
204
205 // Find out how many nodes there are
206 const unsigned n_node = nnode();
207
208 // Get continuous time from timestepper of first node
209 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
210
211 // Set up memory for the shape functions
212 Shape psi(n_node);
213
214 // Local and global coordinates
215 Vector<double> s(n_dim - 1);
217
218 // Loop over plot points
219 unsigned num_plot_points = this->nplot_points_paraview(nplot);
220 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
221 {
222 // Get local coordinates of plot point
223 this->get_s_plot(iplot, nplot, s);
224
225 // Outer unit normal
226 Vector<double> unit_normal(n_dim);
227 outer_unit_normal(s, unit_normal);
228
229 // Find the shape functions
230 shape(s, psi);
231
232 // Initialise to zero
233 for (unsigned i = 0; i < n_dim; i++)
234 {
235 interpolated_x[i] = 0.0;
236 }
237
238 // Calculate stuff
239 for (unsigned l = 0; l < n_node; l++)
240 {
241 // Loop over directions
242 for (unsigned i = 0; i < n_dim; i++)
243 {
244 interpolated_x[i] += this->nodal_position(l, i) * psi[l];
245 }
246 }
247
248 // Get the imposed traction
250
251 // Dummy integration point
252 unsigned ipt = 0;
253 get_traction(time, ipt, interpolated_x, unit_normal, traction);
254
255 // Traction components
256 if (k < n_dim + 1)
257 {
258 file_out << traction[k] << std::endl;
259 }
260 // Advection Diffusion
261 else if (k < 2 * n_dim + 1 && k >= n_dim + 1)
262 {
263 file_out << unit_normal[k] << std::endl;
264 }
265 // Never get here
266 else
267 {
268 std::stringstream error_stream;
269 error_stream
270 << "Axisymmetric Fluid Traction Navier-Stokes Elements only store "
271 << 2 * (n_dim + 1) << " fields " << std::endl;
272 throw OomphLibError(error_stream.str(),
273 OOMPH_CURRENT_FUNCTION,
274 OOMPH_EXCEPTION_LOCATION);
275 }
276 }
277 }
278
279 /// Name of the i-th scalar field. Default implementation
280 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
281 /// overloaded with more meaningful names in specific elements.
282 std::string scalar_name_paraview(const unsigned& i) const
283 {
284 // Number of dimensions
285 unsigned n_dim = this->nodal_dimension();
286
287 // Traction components
288 if (i < n_dim + 1)
289 {
290 return "Traction component " + StringConversion::to_string(i);
291 }
292 // Normals
293 else if (i < 2 * n_dim + 1 && i >= n_dim + 1)
294 {
295 return "Normal " + StringConversion::to_string(i % (n_dim + 1));
296 }
297 // Never get here
298 else
299 {
300 std::stringstream error_stream;
301 error_stream
302 << "Axisymmetric Fluid Traction Navier-Stokes Elements only store "
303 << 2 * (n_dim + 1) << " fields " << std::endl;
304 throw OomphLibError(
305 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
306 }
307 }
308
309 /// Output function
310 void output(std::ostream& outfile, const unsigned& n_plot)
311 {
312 // Number of dimensions
313 unsigned n_dim = this->nodal_dimension();
314
315 // Find out how many nodes there are
316 const unsigned n_node = nnode();
317
318 // Get continuous time from timestepper of first node
319 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
320
321 // Set up memory for the shape functions
322 Shape psi(n_node);
323
324 // Local and global coordinates
325 Vector<double> s(n_dim - 1);
327
328 // Tecplot header info
329 outfile << this->tecplot_zone_string(n_plot);
330
331 // Loop over plot points
332 unsigned num_plot_points = this->nplot_points(n_plot);
333 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
334 {
335 // Get local coordinates of plot point
336 this->get_s_plot(iplot, n_plot, s);
337
338 // Outer unit normal
339 Vector<double> unit_normal(n_dim);
340 outer_unit_normal(s, unit_normal);
341
342 // Find the shape functions
343 shape(s, psi);
344
345 // Initialise to zero
346 for (unsigned i = 0; i < n_dim; i++)
347 {
348 interpolated_x[i] = 0.0;
349 }
350
351 // Calculate stuff
352 for (unsigned l = 0; l < n_node; l++)
353 {
354 // Loop over directions
355 for (unsigned i = 0; i < n_dim; i++)
356 {
357 interpolated_x[i] += this->nodal_position(l, i) * psi[l];
358 }
359 }
360
361 // Get the imposed traction
363
364 // Dummy integration point
365 unsigned ipt = 0;
366 get_traction(time, ipt, interpolated_x, unit_normal, traction);
367
368 // Output the x,y,..
369 for (unsigned i = 0; i < n_dim; i++)
370 {
371 outfile << interpolated_x[i] << " ";
372 }
373
374 // Output the traction components
375 for (unsigned i = 0; i < n_dim + 1; i++)
376 {
377 outfile << traction[i] << " ";
378 }
379
380 // Output normal
381 for (unsigned i = 0; i < n_dim; i++)
382 {
383 outfile << unit_normal[i] << " ";
384 }
385 outfile << std::endl;
386 }
387 }
388
389 /// C_style output function
390 void output(FILE* file_pt)
391 {
392 FiniteElement::output(file_pt);
393 }
394
395 /// C-style output function
396 void output(FILE* file_pt, const unsigned& n_plot)
397 {
398 FiniteElement::output(file_pt, n_plot);
399 }
400
401
402 /// Compute traction vector at specified local coordinate
403 /// Should only be used for post-processing; ignores dependence
404 /// on integration point!
405 void traction(const double& time,
406 const Vector<double>& s,
408 };
409
410 /// ////////////////////////////////////////////////////////////////////
411 /// ////////////////////////////////////////////////////////////////////
412 /// ////////////////////////////////////////////////////////////////////
413
414 //=====================================================================
415 /// Compute traction vector at specified local coordinate
416 /// Should only be used for post-processing; ignores dependence
417 /// on integration point!
418 //=====================================================================
419 template<class ELEMENT>
421 const double& time, const Vector<double>& s, Vector<double>& traction)
422 {
423 unsigned n_dim = this->nodal_dimension();
424
425 // Position vector
426 Vector<double> x(n_dim);
427 interpolated_x(s, x);
428
429 // Outer unit normal (only in r and z direction!)
430 Vector<double> unit_normal(n_dim);
431 outer_unit_normal(s, unit_normal);
432
433 // Dummy
434 unsigned ipt = 0;
435
436 // Traction vector
437 get_traction(time, ipt, x, unit_normal, traction);
438 }
439
440
441 //=====================================================================
442 /// Return the residuals for the
443 /// AxisymmetricNavierStokesTractionElement equations
444 //=====================================================================
445 template<class ELEMENT>
448 Vector<double>& residuals)
449 {
450 // Find out how many nodes there are
451 unsigned n_node = nnode();
452
453 // Get continuous time from timestepper of first node
454 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
455
456#ifdef PARANOID
457 // Find out how many positional dofs there are
458 unsigned n_position_type = this->nnodal_position_type();
459 if (n_position_type != 1)
460 {
461 throw OomphLibError("AxisymmetricNavierStokes is not yet implemented for "
462 "more than one position type",
463 OOMPH_CURRENT_FUNCTION,
464 OOMPH_EXCEPTION_LOCATION);
465 }
466#endif
467
468 // Find out the dimension of the node
469 unsigned n_dim = this->nodal_dimension();
470
471 // Cache the nodal indices at which the velocity components are stored
472 unsigned u_nodal_index[n_dim + 1];
473 for (unsigned i = 0; i < n_dim + 1; i++)
474 {
475 u_nodal_index[i] = this->U_index_axisymmetric_nst_traction[i];
476 }
477
478 // Integer to hold the local equation number
479 int local_eqn = 0;
480
481 // Set up memory for the shape functions
482 // Note that in this case, the number of lagrangian coordinates is always
483 // equal to the dimension of the nodes
484 Shape psi(n_node);
485 DShape dpsids(n_node, n_dim - 1);
486
487 // Set the value of n_intpt
488 unsigned n_intpt = integral_pt()->nweight();
489
490 // Loop over the integration points
491 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
492 {
493 // Get the integral weight
494 double w = integral_pt()->weight(ipt);
495
496 // Only need to call the local derivatives
497 dshape_local_at_knot(ipt, psi, dpsids);
498
499 // Calculate the Eulerian and Lagrangian coordinates
500 Vector<double> interpolated_x(n_dim, 0.0);
501
502 // Also calculate the surface Vectors (derivatives wrt local coordinates)
503 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
504
505 // Calculate positions and derivatives
506 for (unsigned l = 0; l < n_node; l++)
507 {
508 // Loop over directions
509 for (unsigned i = 0; i < n_dim; i++)
510 {
511 // Calculate the Eulerian coords
512 const double x_local = nodal_position(l, i);
513 interpolated_x[i] += x_local * psi(l);
514
515 // Loop over LOCAL derivative directions, to calculate the tangent(s)
516 for (unsigned j = 0; j < n_dim - 1; j++)
517 {
518 interpolated_A(j, i) += x_local * dpsids(l, j);
519 }
520 }
521 }
522
523 // Now find the local metric tensor from the tangent Vectors
524 DenseMatrix<double> A(n_dim - 1);
525 for (unsigned i = 0; i < n_dim - 1; i++)
526 {
527 for (unsigned j = 0; j < n_dim - 1; j++)
528 {
529 // Initialise surface metric tensor to zero
530 A(i, j) = 0.0;
531
532 // Take the dot product
533 for (unsigned k = 0; k < n_dim; k++)
534 {
535 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
536 }
537 }
538 }
539
540 // Get the outer unit normal
541 Vector<double> interpolated_normal(n_dim);
542 outer_unit_normal(ipt, interpolated_normal);
543
544 // Find the determinant of the metric tensor
545 double Adet = 0.0;
546 switch (n_dim)
547 {
548 case 2:
549 Adet = A(0, 0);
550 break;
551 case 3:
552 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
553 break;
554 default:
555 throw OomphLibError(
556 "Wrong dimension in AxisymmetricNavierStokesTractionElement",
557 "AxisymmetricNavierStokesTractionElement::fill_in_contribution_to_"
558 "residuals()",
559 OOMPH_EXCEPTION_LOCATION);
560 }
561
562 // Premultiply the weights and the square-root of the determinant of
563 // the metric tensor
564 double W = w * sqrt(Adet);
565
566 // Now calculate the load
567 Vector<double> traction(n_dim + 1);
568 get_traction(time, ipt, interpolated_x, interpolated_normal, traction);
569
570 // Loop over the test functions, nodes of the element
571 for (unsigned l = 0; l < n_node; l++)
572 {
573 // Loop over the velocity components
574 for (unsigned i = 0; i < n_dim + 1; i++)
575 {
576 // Equation number
577 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
578 /*IF it's not a boundary condition*/
579 if (local_eqn >= 0)
580 {
581 // Add the loading terms to the residuals
582 residuals[local_eqn] -=
583 traction[i] * psi(l) * interpolated_x[0] * W;
584 }
585 }
586 } // End of loop over shape functions
587 } // End of loop over integration points
588 }
589
590
591 /// //////////////////////////////////////////////////////////////////////
592 /// //////////////////////////////////////////////////////////////////////
593 /// //////////////////////////////////////////////////////////////////////
594
595
596 //=======================================================================
597 /// Namespace containing the default Strouhal number of axisymmetric
598 /// linearised FSI.
599 //=======================================================================
600 namespace LinearisedFSIAxisymmetricNStNoSlipBCHelper
601 {
602 /// Default for fluid Strouhal number
603 extern double Default_strouhal_number;
604
605 } // namespace LinearisedFSIAxisymmetricNStNoSlipBCHelper
606
607
608 //======================================================================
609 /// A class for elements that allow the imposition of the linearised
610 /// FSI no slip condition from an adjacent linearly elastic axisymmetric
611 /// solid. The element geometry is obtained from the FaceGeometry<ELEMENT>
612 /// policy class.
613 //======================================================================
614 template<class FLUID_BULK_ELEMENT, class SOLID_BULK_ELEMENT>
616 : public virtual FaceGeometry<FLUID_BULK_ELEMENT>,
617 public virtual FaceElement,
618 public virtual ElementWithExternalElement
619 {
620 public:
621 /// Constructor, takes the pointer to the "bulk" element and the
622 /// face index identifying the face to which the element is attached.
623 /// The optional identifier can be used
624 /// to distinguish the additional nodal values created by
625 /// this element from thos created by other FaceElements.
627 FiniteElement* const& bulk_el_pt,
628 const int& face_index,
629 const unsigned& id = 0);
630
631 /// Broken copy constructor
634
635 /// Broken assignment operator
636 // Commented out broken assignment operator because this can lead to a
637 // conflict warning when used in the virtual inheritence hierarchy.
638 // Essentially the compiler doesn't realise that two separate
639 // implementations of the broken function are the same and so, quite
640 // rightly, it shouts.
641 /*void operator=(const LinearisedFSIAxisymmetricNStNoSlipBCElementElement&)
642 = delete;*/
643
644
645 /// Access function for the pointer to the fluid Strouhal number
646 /// (if not set, St defaults to 1)
647 double*& st_pt()
648 {
649 return St_pt;
650 }
651
652 /// Access function for the fluid Strouhal number
653 double st() const
654 {
655 return *St_pt;
656 }
657
658 /// Add the element's contribution to its residual vector
660 {
661 // Call the generic residuals function with flag set to 0
662 // using a dummy matrix argument
665 }
666
667
668 /// Add the element's contribution to its residual vector and its
669 /// Jacobian matrix
671 DenseMatrix<double>& jacobian)
672 {
673 // Call the generic routine with the flag set to 1
675 residuals, jacobian, 1);
676
677 // Derivatives w.r.t. external data
679 }
680
681 /// Output function
682 void output(std::ostream& outfile)
683 {
684 // Dummy
685 unsigned nplot = 0;
686 output(outfile, nplot);
687 }
688
689 /// Output function: Output at Gauss points; n_plot is ignored.
690 void output(std::ostream& outfile, const unsigned& n_plot)
691 {
692 outfile << "ZONE\n";
693
694 // Get the value of Nintpt
695 const unsigned n_intpt = integral_pt()->nweight();
696
697 // Set the Vector to hold local coordinates
698 Vector<double> s_int(Dim - 1);
699
700 // Loop over the integration points
701 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
702 {
703 // Assign values of s
704 for (unsigned i = 0; i < (Dim - 1); i++)
705 {
706 s_int[i] = integral_pt()->knot(ipt, i);
707 }
708
709 // Get boundary coordinate
710 Vector<double> zeta(1);
711 interpolated_zeta(s_int, zeta);
712
713 // Get velocity from adjacent solid
714 SOLID_BULK_ELEMENT* ext_el_pt =
715 dynamic_cast<SOLID_BULK_ELEMENT*>(external_element_pt(0, ipt));
717 Vector<double> dudt(3);
718 ext_el_pt->interpolated_du_dt_axisymmetric_linear_elasticity(s_ext,
719 dudt);
720
721 // Output
722 outfile << ext_el_pt->interpolated_x(s_ext, 0) << " "
723 << ext_el_pt->interpolated_x(s_ext, 1) << " " << dudt[0] << " "
724 << dudt[1] << " " << dudt[2] << " " << zeta[0] << std::endl;
725 }
726 }
727
728
729 /// C-style output function
730 void output(FILE* file_pt)
731 {
733 }
734
735 /// C-style output function
736 void output(FILE* file_pt, const unsigned& n_plot)
737 {
739 }
740
741
742 protected:
743 /// Function to compute the shape and test functions and to return
744 /// the Jacobian of mapping between local and global (Eulerian)
745 /// coordinates
746 inline double shape_and_test(const Vector<double>& s,
747 Shape& psi,
748 Shape& test) const
749 {
750 // Find number of nodes
751 unsigned n_node = nnode();
752
753 // Get the shape functions
754 shape(s, psi);
755
756 // Set the test functions to be the same as the shape functions
757 for (unsigned i = 0; i < n_node; i++)
758 {
759 test[i] = psi[i];
760 }
761
762 // Return the value of the jacobian
763 return J_eulerian(s);
764 }
765
766
767 /// Function to compute the shape and test functions and to return
768 /// the Jacobian of mapping between local and global (Eulerian)
769 /// coordinates
770 inline double shape_and_test_at_knot(const unsigned& ipt,
771 Shape& psi,
772 Shape& test) const
773 {
774 // Find number of nodes
775 unsigned n_node = nnode();
776
777 // Get the shape functions
778 shape_at_knot(ipt, psi);
779
780 // Set the test functions to be the same as the shape functions
781 for (unsigned i = 0; i < n_node; i++)
782 {
783 test[i] = psi[i];
784 }
785
786 // Return the value of the jacobian
787 return J_eulerian_at_knot(ipt);
788 }
789
790
791 private:
792 /// Add the element's contribution to its residual vector.
793 /// flag=1(or 0): do (or don't) compute the contribution to the
794 /// Jacobian as well.
796 Vector<double>& residuals,
797 DenseMatrix<double>& jacobian,
798 const unsigned& flag);
799
800 /// The spatial dimension of the problem
801 unsigned Dim;
802
803 /// The index at which the unknowns are stored at the nodes
805
806 /// Lagrange Id
807 unsigned Id;
808
809 /// Pointer to fluid Strouhal number
810 double* St_pt;
811 };
812
813 /// ///////////////////////////////////////////////////////////////////
814 /// ///////////////////////////////////////////////////////////////////
815 /// ///////////////////////////////////////////////////////////////////
816
817
818 //===========================================================================
819 /// Constructor, takes the pointer to the "bulk" element, and the
820 /// face index that identifies the face of the bulk element to which
821 /// this face element is to be attached.
822 /// The optional identifier can be used
823 /// to distinguish the additional nodal values created by
824 /// this element from thos created by other FaceElements.
825 //===========================================================================
826 template<class FLUID_BULK_ELEMENT, class SOLID_BULK_ELEMENT>
828 SOLID_BULK_ELEMENT>::
829 LinearisedFSIAxisymmetricNStNoSlipBCElementElement(
830 FiniteElement* const& bulk_el_pt,
831 const int& face_index,
832 const unsigned& id)
833 : FaceGeometry<FLUID_BULK_ELEMENT>(), FaceElement()
834 {
835 // Set source element storage: one interaction with an external element
836 // that provides the velocity of the adjacent linear elasticity
837 // element
838 this->set_ninteraction(1);
839
840 // Store the ID of the FaceElement -- this is used to distinguish
841 // it from any others
842 Id = id;
843
844 // Initialise pointer to fluid Strouhal number. Defaults to 1
845 St_pt =
847
848 // Let the bulk element build the FaceElement, i.e. setup the pointers
849 // to its nodes (by referring to the appropriate nodes in the bulk
850 // element), etc.
851 bulk_el_pt->build_face_element(face_index, this);
852
853 // Extract the dimension of the problem from the dimension of
854 // the first node
855 Dim = this->node_pt(0)->ndim();
856
857 // Read the index from the (cast) bulk element.
859 for (unsigned i = 0; i < 3; i++)
860 {
862 dynamic_cast<FLUID_BULK_ELEMENT*>(bulk_el_pt)->u_index_axi_nst(i);
863 }
864
865 // We need Dim+1 additional values for each FaceElement node
866 // to store the Lagrange multipliers.
867 Vector<unsigned> n_additional_values(nnode(), Dim + 1);
868
869 // Now add storage for Lagrange multipliers and set the map containing
870 // the position of the first entry of this face element's
871 // additional values.
872 add_additional_values(n_additional_values, id);
873 }
874
875
876 //===========================================================================
877 /// Helper function to compute the element's residual vector and
878 /// the Jacobian matrix.
879 //===========================================================================
880 template<class FLUID_BULK_ELEMENT, class SOLID_BULK_ELEMENT>
882 SOLID_BULK_ELEMENT>::
883 fill_in_generic_residual_contribution_fsi_no_slip_axisym(
884 Vector<double>& residuals,
885 DenseMatrix<double>& jacobian,
886 const unsigned& flag)
887 {
888 // Find out how many nodes there are
889 const unsigned n_node = nnode();
890
891 // Set up memory for the shape and test functions
892 Shape psif(n_node), testf(n_node);
893
894 // Set the value of Nintpt
895 const unsigned n_intpt = integral_pt()->nweight();
896
897 // Set the Vector to hold local coordinates
898 Vector<double> s(Dim - 1);
899
900 // Cache the Strouhal number
901 const double local_st = st();
902
903 // Integers to hold the local equation and unknown numbers
904 int local_eqn = 0;
905
906 // Loop over the integration points
907 //--------------------------------
908 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
909 {
910 // Assign values of s
911 for (unsigned i = 0; i < (Dim - 1); i++)
912 {
913 s[i] = integral_pt()->knot(ipt, i);
914 }
915
916 // Get the integral weight
917 double w = integral_pt()->weight(ipt);
918
919 // Find the shape and test functions and return the Jacobian
920 // of the mapping
921 double J = shape_and_test(s, psif, testf);
922
923 // Premultiply the weights and the Jacobian
924 double W = w * J;
925
926 // Calculate the Lagrange multiplier and the fluid veloc
927 Vector<double> lambda(Dim + 1, 0.0);
928 Vector<double> fluid_veloc(Dim + 1, 0.0);
929
930 // Loop over nodes
931 for (unsigned j = 0; j < n_node; j++)
932 {
933 Node* nod_pt = node_pt(j);
934
935 // Cast to a boundary node
936 BoundaryNodeBase* bnod_pt = dynamic_cast<BoundaryNodeBase*>(node_pt(j));
937
938 // Get the index of the first nodal value associated with
939 // this FaceElement
940 unsigned first_index =
942
943 // Assemble
944 for (unsigned i = 0; i < Dim + 1; i++)
945 {
946 lambda[i] += nod_pt->value(first_index + i) * psif(j);
947 fluid_veloc[i] +=
948 nod_pt->value(U_index_fsi_no_slip_axisym[i]) * psif(j);
949 }
950 }
951
952 // Get velocity from adjacent solid
953 SOLID_BULK_ELEMENT* ext_el_pt =
954 dynamic_cast<SOLID_BULK_ELEMENT*>(external_element_pt(0, ipt));
955 Vector<double> s_ext(external_element_local_coord(0, ipt));
956 Vector<double> dudt(3);
957 ext_el_pt->interpolated_du_dt_axisymmetric_linear_elasticity(s_ext, dudt);
958
959
960 // Now add to the appropriate equations
961
962 // Loop over the test functions
963 for (unsigned l = 0; l < n_node; l++)
964 {
965 // Loop over directions
966 for (unsigned i = 0; i < Dim + 1; i++)
967 {
968 // Add contribution to bulk Navier Stokes equations where
969 //-------------------------------------------------------
970 // the Lagrange multiplier acts as a traction
971 //-------------------------------------------
972 local_eqn = nodal_local_eqn(l, U_index_fsi_no_slip_axisym[i]);
973
974 /*IF it's not a boundary condition*/
975 if (local_eqn >= 0)
976 {
977 // Add the Lagrange multiplier "traction" to the bulk
978 residuals[local_eqn] -= lambda[i] * testf[l] * W;
979
980
981 // Jacobian entries
982 if (flag)
983 {
984 // Loop over the lagrange multiplier unknowns
985 for (unsigned l2 = 0; l2 < n_node; l2++)
986 {
987 // Cast to a boundary node
988 BoundaryNodeBase* bnod_pt =
989 dynamic_cast<BoundaryNodeBase*>(node_pt(l2));
990
991 // Local unknown
992 int local_unknown = nodal_local_eqn(
993 l2,
995 i);
996
997 // If it's not pinned
998 if (local_unknown >= 0)
999 {
1000 jacobian(local_eqn, local_unknown) -= psif[l2] * testf[l] * W;
1001 }
1002 }
1003 }
1004 }
1005
1006 // Now do the Lagrange multiplier equations
1007 //-----------------------------------------
1008 // Cast to a boundary node
1009 BoundaryNodeBase* bnod_pt =
1010 dynamic_cast<BoundaryNodeBase*>(node_pt(l));
1011
1012 // Local eqn number:
1013 int local_eqn = nodal_local_eqn(
1015
1016 // If it's not pinned
1017 if (local_eqn >= 0)
1018 {
1019 residuals[local_eqn] +=
1020 (local_st * dudt[i] - fluid_veloc[i]) * testf(l) * W;
1021
1022 // Jacobian entries
1023 if (flag)
1024 {
1025 // Loop over the velocity unknowns [derivs w.r.t. to
1026 // wall velocity taken care of by fd-ing
1027 for (unsigned l2 = 0; l2 < n_node; l2++)
1028 {
1029 int local_unknown =
1030 nodal_local_eqn(l2, U_index_fsi_no_slip_axisym[i]);
1031
1032 /*IF it's not a boundary condition*/
1033 if (local_unknown >= 0)
1034 {
1035 jacobian(local_eqn, local_unknown) -= psif[l2] * testf[l] * W;
1036 }
1037 }
1038 }
1039 }
1040 }
1041 }
1042 }
1043 }
1044
1045} // namespace oomph
1046
1047
1048#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class for elements that allow the imposition of an applied traction in the axisym Navier Stokes eqn...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &k, const unsigned &nplot) const
Write values of the k-th scalar field at the plot points. Needs to be implemented for each new specif...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class.
void traction(const double &time, const Vector< double > &s, Vector< double > &traction)
Compute traction vector at specified local coordinate Should only be used for post-processing; ignore...
Vector< unsigned > U_index_axisymmetric_nst_traction
Index at which the i-th velocity component is stored.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function. Arguments: Eulerian coordinate; outer unit normal; applied t...
virtual void get_traction(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) const
Get the traction vector: Pass number of integration point (dummy), Eulerian coordinate and normal vec...
AxisymmetricNavierStokesTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void output(FILE *file_pt)
C_style output function.
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
void fill_in_contribution_to_residuals_axisymmetric_nst_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void output(std::ostream &outfile)
Output function.
A class that contains the information required by Nodes that are located on Mesh boundaries....
Definition: nodes.h:1996
unsigned index_of_first_value_assigned_by_face_element(const unsigned &face_id=0) const
Return the index of the first value associated with the i-th face element value. If no argument is sp...
Definition: nodes.h:2061
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
void fill_in_jacobian_from_external_interaction_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from all external interaction degrees of freedom (geometr...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
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
void add_additional_values(const Vector< unsigned > &nadditional_values, const unsigned &id)
Helper function adding additional values for the unknowns associated with the FaceElement....
Definition: elements.h:4428
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
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
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5328
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2862
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4675
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
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 nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3220
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
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.
A class for elements that allow the imposition of the linearised FSI no slip condition from an adjace...
Vector< unsigned > U_index_fsi_no_slip_axisym
The index at which the unknowns are stored at the nodes.
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Output at Gauss points; n_plot is ignored.
LinearisedFSIAxisymmetricNStNoSlipBCElementElement(const LinearisedFSIAxisymmetricNStNoSlipBCElementElement &dummy)=delete
Broken copy constructor.
LinearisedFSIAxisymmetricNStNoSlipBCElementElement(FiniteElement *const &bulk_el_pt, const int &face_index, const unsigned &id=0)
Constructor, takes the pointer to the "bulk" element and the face index identifying the face to which...
void fill_in_generic_residual_contribution_fsi_no_slip_axisym(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the contri...
double st() const
Access function for the fluid Strouhal number.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
Definition: nodes.cc:2408
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
void Zero_traction_fct(const double &time, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
void output()
Doc the command line arguments.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double Default_strouhal_number
Default for fluid Strouhal number.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...