axisym_linear_elasticity_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 apply surface loads to
27// the equations of axisymmetric linear elasticity
28
29#ifndef OOMPH_AXISYMMETRIC_LINEAR_ELASTICITY_TRACTION_ELEMENTS_HEADER
30#define OOMPH_AXISYMMETRIC_LINEAR_ELASTICITY_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/Qelements.h"
42
43
44namespace oomph
45{
46 //=======================================================================
47 /// Namespace containing the zero traction function for
48 /// axisymmetric linear elasticity traction elements
49 //=======================================================================
50 namespace AxisymmetricLinearElasticityTractionElementHelper
51 {
52 //=======================================================================
53 /// Default load function (zero traction)
54 //=======================================================================
55 void Zero_traction_fct(const double& time,
56 const Vector<double>& x,
57 const Vector<double>& N,
58 Vector<double>& load)
59 {
60 unsigned n_dim = load.size();
61 for (unsigned i = 0; i < n_dim; i++)
62 {
63 load[i] = 0.0;
64 }
65 }
66
67 } // namespace AxisymmetricLinearElasticityTractionElementHelper
68
69
70 //======================================================================
71 /// A class for elements that allow the imposition of an applied traction
72 /// in the equations of axisymmetric linear elasticity.
73 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
74 /// class and and thus, we can be generic enough without the need to have
75 /// a separate equations class.
76 //======================================================================
77 template<class ELEMENT>
79 : public virtual FaceGeometry<ELEMENT>,
80 public virtual FaceElement
81 {
82 protected:
83 /// Index at which the i-th displacement component is stored
85
86 /// Pointer to an imposed traction function. Arguments:
87 /// Eulerian coordinate; outer unit normal;
88 /// applied traction. (Not all of the input arguments will be
89 /// required for all specific load functions but the list should
90 /// cover all cases)
91 void (*Traction_fct_pt)(const double& time,
92 const Vector<double>& x,
93 const Vector<double>& n,
94 Vector<double>& result);
95
96
97 /// Get the traction vector: Pass number of integration point
98 /// (dummy), Eulerian coordinate and normal vector and return the load
99 /// vector (not all of the input arguments will be required for all specific
100 /// load functions but the list should cover all cases). This function is
101 /// virtual so it can be overloaded for FSI.
102 virtual void get_traction(const double& time,
103 const unsigned& intpt,
104 const Vector<double>& x,
105 const Vector<double>& n,
107 {
108 Traction_fct_pt(time, x, n, traction);
109 }
110
111
112 /// Helper function that actually calculates the residuals
113 // This small level of indirection is required to avoid calling
114 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
115 // which causes all kinds of pain if overloading later on
117 Vector<double>& residuals);
118
119
120 public:
121 /// Constructor, which takes a "bulk" element and the
122 /// value of the index and its limit
124 FiniteElement* const& element_pt, const int& face_index)
125 : FaceGeometry<ELEMENT>(), FaceElement()
126 {
127 // Attach the geometrical information to the element. N.B. This function
128 // also assigns nbulk_value from the required_nvalue of the bulk element
129 element_pt->build_face_element(face_index, this);
130
131 // Find the dimension of the problem
132 unsigned n_dim = element_pt->nodal_dimension();
133
134 // Find the index at which the displacement unknowns are stored
135 ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
136 this->U_index_axisymmetric_linear_elasticity_traction.resize(n_dim + 1);
137 for (unsigned i = 0; i < n_dim + 1; i++)
138 {
139 this->U_index_axisymmetric_linear_elasticity_traction[i] =
140 cast_element_pt->u_index_axisymmetric_linear_elasticity(i);
141 }
142
143 // Zero traction
146 }
147
148
149 /// Reference to the traction function pointer
150 void (*&traction_fct_pt())(const double& time,
151 const Vector<double>& x,
152 const Vector<double>& n,
154 {
155 return Traction_fct_pt;
156 }
157
158
159 /// Return the residuals
161 {
163 residuals);
164 }
165
166
167 /// Fill in contribution from Jacobian
169 DenseMatrix<double>& jacobian)
170 {
171 // Call the residuals
173 residuals);
174 }
175
176 /// Specify the value of nodal zeta from the face geometry
177 /// The "global" intrinsic coordinate of the element when
178 /// viewed as part of a geometric object should be given by
179 /// the FaceElement representation, by default (needed to break
180 /// indeterminacy if bulk element is SolidElement)
181 double zeta_nodal(const unsigned& n,
182 const unsigned& k,
183 const unsigned& i) const
184 {
185 return FaceElement::zeta_nodal(n, k, i);
186 }
187
188 /// Output function
189 void output(std::ostream& outfile)
190 {
191 unsigned nplot = 5;
192 output(outfile, nplot);
193 }
194
195 /// Output function
196 void output(std::ostream& outfile, const unsigned& n_plot)
197 {
198 // Number of dimensions
199 unsigned n_dim = this->nodal_dimension();
200
201 // Find out how many nodes there are
202 const unsigned n_node = nnode();
203
204 // Get continuous time from timestepper of first node
205 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
206
207 // Set up memory for the shape functions
208 Shape psi(n_node);
209
210 // Local and global coordinates
211 Vector<double> s(n_dim - 1);
213
214 // Tecplot header info
215 outfile << this->tecplot_zone_string(n_plot);
216
217 // Loop over plot points
218 unsigned num_plot_points = this->nplot_points(n_plot);
219 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
220 {
221 // Get local coordinates of plot point
222 this->get_s_plot(iplot, n_plot, s);
223
224 // Outer unit normal
225 Vector<double> unit_normal(n_dim);
226 outer_unit_normal(s, unit_normal);
227
228 // Find the shape functions
229 shape(s, psi);
230
231 // Initialise to zero
232 for (unsigned i = 0; i < n_dim; i++)
233 {
234 interpolated_x[i] = 0.0;
235 }
236
237 // Calculate stuff
238 for (unsigned l = 0; l < n_node; l++)
239 {
240 // Loop over directions
241 for (unsigned i = 0; i < n_dim; i++)
242 {
243 interpolated_x[i] += nodal_position(l, i) * psi[l];
244 }
245 }
246
247 // Get the imposed flux
249
250 // Dummy integration point
251 unsigned ipt = 0;
252
253 get_traction(time, ipt, interpolated_x, unit_normal, traction);
254
255 // Output the x,y,..
256 for (unsigned i = 0; i < n_dim; i++)
257 {
258 outfile << interpolated_x[i] << " ";
259 }
260
261 // Output the traction components
262 for (unsigned i = 0; i < n_dim + 1; i++)
263 {
264 outfile << traction[i] << " ";
265 }
266
267 // Output normal
268 for (unsigned i = 0; i < n_dim; i++)
269 {
270 outfile << unit_normal[i] << " ";
271 }
272 outfile << std::endl;
273 }
274 }
275
276 /// C_style output function
277 void output(FILE* file_pt)
278 {
279 FiniteElement::output(file_pt);
280 }
281
282 /// C-style output function
283 void output(FILE* file_pt, const unsigned& n_plot)
284 {
285 FiniteElement::output(file_pt, n_plot);
286 }
287
288
289 /// Compute traction vector at specified local coordinate
290 /// Should only be used for post-processing; ignores dependence
291 /// on integration point!
292 void traction(const double& time,
293 const Vector<double>& s,
295 };
296
297 /// ////////////////////////////////////////////////////////////////////
298 /// ////////////////////////////////////////////////////////////////////
299 /// ////////////////////////////////////////////////////////////////////
300
301 //=====================================================================
302 /// Compute traction vector at specified local coordinate
303 /// Should only be used for post-processing; ignores dependence
304 /// on integration point!
305 //=====================================================================
306 template<class ELEMENT>
308 const double& time, const Vector<double>& s, Vector<double>& traction)
309 {
310 unsigned n_dim = this->nodal_dimension();
311
312 // Position vector
313 Vector<double> x(n_dim);
314 interpolated_x(s, x);
315
316 // Outer unit normal (only in r and z direction!)
317 Vector<double> unit_normal(n_dim);
318 outer_unit_normal(s, unit_normal);
319
320 // Dummy
321 unsigned ipt = 0;
322
323 // Traction vector
324 get_traction(time, ipt, x, unit_normal, traction);
325 }
326
327
328 //=====================================================================
329 /// Return the residuals for the
330 /// AxisymmetricLinearElasticityTractionElement equations
331 //=====================================================================
332 template<class ELEMENT>
335 Vector<double>& residuals)
336 {
337 // Find out how many nodes there are
338 unsigned n_node = nnode();
339
340 // Get continuous time from timestepper of first node
341 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
342
343#ifdef PARANOID
344 // Find out how many positional dofs there are
345 unsigned n_position_type = this->nnodal_position_type();
346 if (n_position_type != 1)
347 {
348 throw OomphLibError("AxisymmetricLinearElasticity is not yet implemented "
349 "for more than one position type",
350 OOMPH_CURRENT_FUNCTION,
351 OOMPH_EXCEPTION_LOCATION);
352 }
353#endif
354
355 // Find out the dimension of the node
356 unsigned n_dim = this->nodal_dimension();
357
358 // Cache the nodal indices at which the displacement components are stored
359 unsigned u_nodal_index[n_dim + 1];
360 for (unsigned i = 0; i < n_dim + 1; i++)
361 {
362 u_nodal_index[i] =
363 this->U_index_axisymmetric_linear_elasticity_traction[i];
364 }
365
366 // Integer to hold the local equation number
367 int local_eqn = 0;
368
369 // Set up memory for the shape functions
370 // Note that in this case, the number of lagrangian coordinates is always
371 // equal to the dimension of the nodes
372 Shape psi(n_node);
373 DShape dpsids(n_node, n_dim - 1);
374
375 // Set the value of n_intpt
376 unsigned n_intpt = integral_pt()->nweight();
377
378 // Loop over the integration points
379 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
380 {
381 // Get the integral weight
382 double w = integral_pt()->weight(ipt);
383
384 // Only need to call the local derivatives
385 dshape_local_at_knot(ipt, psi, dpsids);
386
387 // Calculate the Eulerian and Lagrangian coordinates
388 Vector<double> interpolated_x(n_dim, 0.0);
389
390 // Also calculate the surface Vectors (derivatives wrt local coordinates)
391 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
392
393 // Calculate displacements and derivatives
394 for (unsigned l = 0; l < n_node; l++)
395 {
396 // Loop over directions
397 for (unsigned i = 0; i < n_dim; i++)
398 {
399 // Calculate the Eulerian coords
400 const double x_local = nodal_position(l, i);
401 interpolated_x[i] += x_local * psi(l);
402
403 // Loop over LOCAL derivative directions, to calculate the tangent(s)
404 for (unsigned j = 0; j < n_dim - 1; j++)
405 {
406 interpolated_A(j, i) += x_local * dpsids(l, j);
407 }
408 }
409 }
410
411 // Now find the local metric tensor from the tangent Vectors
412 DenseMatrix<double> A(n_dim - 1);
413 for (unsigned i = 0; i < n_dim - 1; i++)
414 {
415 for (unsigned j = 0; j < n_dim - 1; j++)
416 {
417 // Initialise surface metric tensor to zero
418 A(i, j) = 0.0;
419
420 // Take the dot product
421 for (unsigned k = 0; k < n_dim; k++)
422 {
423 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
424 }
425 }
426 }
427
428 // Get the outer unit normal
429 Vector<double> interpolated_normal(n_dim);
430 outer_unit_normal(ipt, interpolated_normal);
431
432 // Find the determinant of the metric tensor
433 double Adet = 0.0;
434 switch (n_dim)
435 {
436 case 2:
437 Adet = A(0, 0);
438 break;
439 case 3:
440 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
441 break;
442 default:
443 throw OomphLibError(
444 "Wrong dimension in AxisymmetricLinearElasticityTractionElement",
445 "AxisymmetricLinearElasticityTractionElement::fill_in_contribution_"
446 "to_residuals()",
447 OOMPH_EXCEPTION_LOCATION);
448 }
449
450 // Premultiply the weights and the square-root of the determinant of
451 // the metric tensor
452 double W = w * sqrt(Adet);
453
454 // Now calculate the load
455 Vector<double> traction(n_dim + 1);
456 get_traction(time, ipt, interpolated_x, interpolated_normal, traction);
457
458 // Loop over the test functions, nodes of the element
459 for (unsigned l = 0; l < n_node; l++)
460 {
461 // Loop over the displacement components
462 for (unsigned i = 0; i < n_dim + 1; i++)
463 {
464 // Equation number
465 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
466 /*IF it's not a boundary condition*/
467 if (local_eqn >= 0)
468 {
469 // Add the loading terms to the residuals
470 residuals[local_eqn] -=
471 traction[i] * psi(l) * interpolated_x[0] * W;
472 }
473 }
474 } // End of loop over shape functions
475 } // End of loop over integration points
476 }
477
478
479 /// ////////////////////////////////////////////////////////////////////
480 /// ////////////////////////////////////////////////////////////////////
481 /// ////////////////////////////////////////////////////////////////////
482
483
484 //======================================================================
485 /// A class for elements that allow the imposition of an applied traction
486 /// in the equations of axisymmetric linear elasticity from an adjacent
487 /// axisymmetric Navier Stokes element in a linearised FSI problem.
488 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
489 /// class and and thus, we can be generic enough without the need to have
490 /// a separate equations class.
491 //======================================================================
492 template<class ELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
494 : public virtual FaceGeometry<ELASTICITY_BULK_ELEMENT>,
495 public virtual FaceElement,
496 public virtual ElementWithExternalElement
497 {
498 protected:
499 /// Pointer to the ratio, \f$ Q \f$ , of the stress used to
500 /// non-dimensionalise the fluid stresses to the stress used to
501 /// non-dimensionalise the solid stresses.
502 double* Q_pt;
503
504 /// Static default value for the ratio of stress scales
505 /// used in the fluid and solid equations (default is 1.0)
506 static double Default_Q_Value;
507
508 /// Index at which the i-th displacement component is stored
510
511 /// Helper function that actually calculates the residuals
512 // This small level of indirection is required to avoid calling
513 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
514 // which causes all kinds of pain if overloading later on
516 Vector<double>& residuals);
517
518 public:
519 /// Constructor, which takes a "bulk" element and the
520 /// value of the index and its limit
522 FiniteElement* const& element_pt, const int& face_index)
523 : FaceGeometry<ELASTICITY_BULK_ELEMENT>(),
524 FaceElement(),
526 {
527 // Set source element storage: one interaction with an external
528 // element -- the Navier Stokes bulk element that provides the traction
529 this->set_ninteraction(1);
530
531 // Attach the geometrical information to the element. N.B. This function
532 // also assigns nbulk_value from the required_nvalue of the bulk element
533 element_pt->build_face_element(face_index, this);
534
535 // Find the dimension of the problem
536 unsigned n_dim = element_pt->nodal_dimension();
537
538 // Find the index at which the displacement unknowns are stored
539 ELASTICITY_BULK_ELEMENT* cast_element_pt =
540 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
541 this->U_index_axisym_fsi_traction.resize(n_dim + 1);
542 for (unsigned i = 0; i < n_dim + 1; i++)
543 {
544 this->U_index_axisym_fsi_traction[i] =
545 cast_element_pt->u_index_axisymmetric_linear_elasticity(i);
546 }
547 }
548
549
550 /// Return the residuals
552 {
554 }
555
556
557 /// Fill in contribution from Jacobian
559 DenseMatrix<double>& jacobian)
560 {
561 // Call the residuals
563
564 // Derivatives w.r.t. external data
566 }
567
568 /// Return the ratio of the stress scales used to non-dimensionalise
569 /// the fluid and elasticity equations. E.g.
570 /// \f$ Q = (\omega a)^2 \rho/E \f$, i.e. the
571 /// ratio between the inertial fluid stress and the solid's elastic
572 /// modulus E.
573 const double& q() const
574 {
575 return *Q_pt;
576 }
577
578 /// Return a pointer the ratio of stress scales used to
579 /// non-dimensionalise the fluid and solid equations.
580 double*& q_pt()
581 {
582 return Q_pt;
583 }
584
585
586 /// Output function
587 void output(std::ostream& outfile)
588 {
589 /// Dummy
590 unsigned nplot = 0;
591 output(outfile, nplot);
592 }
593
594 /// Output function: Plot traction etc at Gauss points
595 /// nplot is ignored.
596 void output(std::ostream& outfile, const unsigned& n_plot)
597 {
598 // Dimension
599 unsigned n_dim = this->nodal_dimension();
600
601 // Get FSI parameter
602 const double q_value = q();
603
604 // Storage for traction (includes swirl!)
605 Vector<double> traction(n_dim + 1);
606
607 outfile << "ZONE\n";
608
609 // Set the value of n_intpt
610 unsigned n_intpt = integral_pt()->nweight();
611
612 // Loop over the integration points
613 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
614 {
615 Vector<double> s_int(n_dim - 1);
616 s_int[0] = integral_pt()->knot(ipt, 0);
617
618 // Get the outer unit normal
619 Vector<double> interpolated_normal(n_dim);
620 outer_unit_normal(ipt, interpolated_normal);
621
622 // Boundary coordinate
623 Vector<double> zeta(1);
624 interpolated_zeta(s_int, zeta);
625
626 // Get bulk element for traction
627 NAVIER_STOKES_BULK_ELEMENT* ext_el_pt =
628 dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>(
629 this->external_element_pt(0, ipt));
630 Vector<double> s_ext(this->external_element_local_coord(0, ipt));
631
632 // Get traction from bulk element (on fluid scale)
633 ext_el_pt->traction(s_ext, interpolated_normal, traction);
634
635 outfile << ext_el_pt->interpolated_x(s_ext, 0) << " "
636 << ext_el_pt->interpolated_x(s_ext, 1) << " "
637 << q_value * traction[0] << " " << q_value * traction[1] << " "
638 << q_value * traction[2] << " " << interpolated_normal[0] << " "
639 << interpolated_normal[1] << " " << zeta[0] << std::endl;
640 }
641 }
642
643 /// C_style output function
644 void output(FILE* file_pt)
645 {
647 }
648
649 /// C-style output function
650 void output(FILE* file_pt, const unsigned& n_plot)
651 {
653 }
654 };
655
656 /// ////////////////////////////////////////////////////////////////////
657 /// ////////////////////////////////////////////////////////////////////
658 /// ////////////////////////////////////////////////////////////////////
659
660
661 //=================================================================
662 /// Static default value for the ratio of stress scales
663 /// used in the fluid and solid equations (default is 1.0)
664 //=================================================================
665 template<class ELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
666 double FSIAxisymmetricLinearElasticityTractionElement<
667 ELASTICITY_BULK_ELEMENT,
668 NAVIER_STOKES_BULK_ELEMENT>::Default_Q_Value = 1.0;
669
670
671 //=====================================================================
672 /// Return the residuals
673 //=====================================================================
674 template<class ELASTICITY_BULK_ELEMENT, class NAVIER_STOKES_BULK_ELEMENT>
675 void FSIAxisymmetricLinearElasticityTractionElement<
676 ELASTICITY_BULK_ELEMENT,
677 NAVIER_STOKES_BULK_ELEMENT>::
678 fill_in_contribution_to_residuals_axisym_fsi_traction(
679 Vector<double>& residuals)
680 {
681 // Find out how many nodes there are
682 unsigned n_node = nnode();
683
684#ifdef PARANOID
685 // Find out how many positional dofs there are
686 unsigned n_position_type = this->nnodal_position_type();
687 if (n_position_type != 1)
688 {
689 throw OomphLibError("LinearElasticity is not yet implemented for more "
690 "than one position type.",
691 OOMPH_CURRENT_FUNCTION,
692 OOMPH_EXCEPTION_LOCATION);
693 }
694#endif
695
696 // Find out the dimension of the node
697 unsigned n_dim = this->nodal_dimension();
698
699 // Cache the nodal indices at which the displacement components are stored
700 unsigned u_nodal_index[n_dim + 1];
701 for (unsigned i = 0; i < n_dim + 1; i++)
702 {
703 u_nodal_index[i] = this->U_index_axisym_fsi_traction[i];
704 }
705
706 // Integer to hold the local equation number
707 int local_eqn = 0;
708
709 // Set up memory for the shape functions
710 Shape psi(n_node);
711 DShape dpsids(n_node, n_dim - 1);
712
713 // Get FSI parameter
714 const double q_value = q();
715
716 // Storage for traction
717 Vector<double> traction(3);
718
719 // Set the value of n_intpt
720 unsigned n_intpt = integral_pt()->nweight();
721
722 // Loop over the integration points
723 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
724 {
725 // Get the integral weight
726 double w = integral_pt()->weight(ipt);
727
728 // Only need to call the local derivatives
729 dshape_local_at_knot(ipt, psi, dpsids);
730
731 // Calculate the coordinates
732 Vector<double> interpolated_x(n_dim, 0.0);
733
734 // Also calculate the surface tangent vectors
735 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
736
737 // Calculate displacements and derivatives
738 for (unsigned l = 0; l < n_node; l++)
739 {
740 // Loop over directions
741 for (unsigned i = 0; i < n_dim; i++)
742 {
743 // Calculate the Eulerian coords
744 const double x_local = nodal_position(l, i);
745 interpolated_x[i] += x_local * psi(l);
746
747 // Loop over LOCAL derivative directions, to calculate the tangent(s)
748 for (unsigned j = 0; j < n_dim - 1; j++)
749 {
750 interpolated_A(j, i) += x_local * dpsids(l, j);
751 }
752 }
753 }
754
755 // Now find the local metric tensor from the tangent Vectors
756 DenseMatrix<double> A(n_dim - 1);
757 for (unsigned i = 0; i < n_dim - 1; i++)
758 {
759 for (unsigned j = 0; j < n_dim - 1; j++)
760 {
761 // Initialise surface metric tensor to zero
762 A(i, j) = 0.0;
763
764 // Take the dot product
765 for (unsigned k = 0; k < n_dim; k++)
766 {
767 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
768 }
769 }
770 }
771
772 // Get the outer unit normal
773 Vector<double> interpolated_normal(n_dim);
774 outer_unit_normal(ipt, interpolated_normal);
775
776 // Find the determinant of the metric tensor
777 double Adet = 0.0;
778 switch (n_dim)
779 {
780 case 2:
781 Adet = A(0, 0);
782 break;
783 case 3:
784 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
785 break;
786 default:
787 throw OomphLibError(
788 "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
789 "TimeHarmonicLinElastLoadedByPressureElement::fill_in_contribution_"
790 "to_residuals()",
791 OOMPH_EXCEPTION_LOCATION);
792 }
793
794 // Premultiply the weights and the square-root of the determinant of
795 // the metric tensor
796 double W = w * sqrt(Adet);
797
798 // Get bulk element for traction
799 NAVIER_STOKES_BULK_ELEMENT* ext_el_pt =
800 dynamic_cast<NAVIER_STOKES_BULK_ELEMENT*>(
801 this->external_element_pt(0, ipt));
802 Vector<double> s_ext(this->external_element_local_coord(0, ipt));
803
804 // Get traction from bulk element
805 ext_el_pt->traction(s_ext, interpolated_normal, traction);
806
807 // Loop over the test functions, nodes of the element
808 for (unsigned l = 0; l < n_node; l++)
809 {
810 // Loop over the displacement components
811 for (unsigned i = 0; i < n_dim + 1; i++)
812 {
813 // Equation number
814 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
815 /*IF it's not a boundary condition*/
816 if (local_eqn >= 0)
817 {
818 // Add the loading terms (multiplied by fsi scaling factor)
819 // to the residuals
820 residuals[local_eqn] -=
821 q_value * traction[i] * psi(l) * interpolated_x[0] * W;
822 }
823 }
824 } // End of loop over shape functions
825
826 } // End of loop over integration points
827 }
828
829
830 /// ////////////////////////////////////////////////////////////////////
831 /// ////////////////////////////////////////////////////////////////////
832 /// ////////////////////////////////////////////////////////////////////
833
834
835} // namespace oomph
836
837#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 equations of axisymmetri...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
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 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...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
Vector< unsigned > U_index_axisymmetric_linear_elasticity_traction
Index at which the i-th displacement component is stored.
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)
Get the traction vector: Pass number of integration point (dummy), Eulerian coordinate and normal vec...
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
AxisymmetricLinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void fill_in_contribution_to_residuals_axisymmetric_linear_elasticity_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
FSIAxisymmetricLinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations....
Vector< unsigned > U_index_axisym_fsi_traction
Index at which the i-th displacement component is stored.
void fill_in_contribution_to_residuals_axisym_fsi_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations.
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and solid equations (default is...
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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
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 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.
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
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...