fourier_decomposed_helmholtz_time_harmonic_linear_elasticity_interaction.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#ifndef OOMPH_FD_HH_TIME_HARMONIC_LIN_ELAST_HEADER
27#define OOMPH_FD_HH_TIME_HARMONIC_LIN_ELAST_HEADER
28
29
30// Oomph-lib includes
31#include "generic.h"
32#include "fourier_decomposed_helmholtz.h"
33#include "time_harmonic_fourier_decomposed_linear_elasticity.h"
34
35namespace oomph
36{
37 //======================================================================
38 /// A class for elements that allow the imposition of an applied traction
39 /// in the equations of time-harmonic linear elasticity from a Helmholtz
40 /// potential (interpreted as a displacement potential for the fluid in a
41 /// quasi-steady, linearised FSI problem.)
42 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
43 /// class and and thus, we can be generic enough without the need to have
44 /// a separate equations class.
45 //======================================================================
46 template<class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
48 : public virtual FaceGeometry<ELASTICITY_BULK_ELEMENT>,
49 public virtual FaceElement,
50 public virtual ElementWithExternalElement
51 {
52 protected:
53 /// Pointer to the ratio, \f$ Q \f$ , of the stress used to
54 /// non-dimensionalise the fluid stresses to the stress used to
55 /// non-dimensionalise the solid stresses.
56 double* Q_pt;
57
58 /// Static default value for the ratio of stress scales
59 /// used in the fluid and solid equations (default is 1.0)
60 static double Default_Q_Value;
61
62 /// Index at which the i-th displacement component is stored
65
66 /// Helper function that actually calculates the residuals
67 // This small level of indirection is required to avoid calling
68 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
69 // which causes all kinds of pain if overloading later on
71 Vector<double>& residuals);
72
73 public:
74 /// Constructor, which takes a "bulk" element and the
75 /// value of the index and its limit
77 FiniteElement* const& element_pt, const int& face_index)
78 : FaceGeometry<ELASTICITY_BULK_ELEMENT>(),
81 {
82 // Attach the geometrical information to the element. N.B. This function
83 // also assigns nbulk_value from the required_nvalue of the bulk element
84 element_pt->build_face_element(face_index, this);
85
86#ifdef PARANOID
87 {
88 // Check that the element is not a refineable 3d element
89 ELASTICITY_BULK_ELEMENT* elem_pt =
90 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
91 // If it's three-d
92 if (elem_pt->dim() == 3)
93 {
94 // Is it refineable
95 RefineableElement* ref_el_pt =
96 dynamic_cast<RefineableElement*>(elem_pt);
97 if (ref_el_pt != 0)
98 {
99 if (this->has_hanging_nodes())
100 {
101 throw OomphLibError("This flux element will not work correctly "
102 "if nodes are hanging\n",
103 OOMPH_CURRENT_FUNCTION,
104 OOMPH_EXCEPTION_LOCATION);
105 }
106 }
107 }
108 }
109#endif
110
111 // Set source element storage: one interaction with an external
112 // element -- the Helmholtz bulk element that provides the pressure
113 this->set_ninteraction(1);
114
115 // Find the dimension of the problem
116 unsigned n_dim = element_pt->nodal_dimension();
117
118 // Find the index at which the displacement unknowns are stored
119 ELASTICITY_BULK_ELEMENT* cast_element_pt =
120 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
121 this
122 ->U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
123 .resize(n_dim + 1);
124 for (unsigned i = 0; i < n_dim + 1; i++)
125 {
126 this
127 ->U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
128 [i] =
129 cast_element_pt
130 ->u_index_time_harmonic_fourier_decomposed_linear_elasticity(i);
131 }
132 }
133
134
135 /// Return the residuals
137 {
139 }
140
141
142 /// Fill in contribution from Jacobian
144 DenseMatrix<double>& jacobian)
145 {
146 // Call the residuals
148
149 // Derivatives w.r.t. external data
151 }
152
153 /// Return the ratio of the stress scales used to non-dimensionalise
154 /// the fluid and elasticity equations. E.g.
155 /// \f$ Q = (\omega a)^2 \rho/E \f$, i.e. the
156 /// ratio between the inertial fluid stress and the solid's elastic
157 /// modulus E.
158 const double& q() const
159 {
160 return *Q_pt;
161 }
162
163 /// Return a pointer the ratio of stress scales used to
164 /// non-dimensionalise the fluid and solid equations.
165 double*& q_pt()
166 {
167 return Q_pt;
168 }
169
170
171 /// Output function
172 void output(std::ostream& outfile)
173 {
174 /// Dummy
175 unsigned nplot = 0;
176 output(outfile, nplot);
177 }
178
179 /// Output function: Plot traction etc at Gauss points
180 /// nplot is ignored.
181 void output(std::ostream& outfile, const unsigned& n_plot)
182 {
183 // Dimension
184 unsigned n_dim = this->nodal_dimension();
185
186 // Storage for traction
187 Vector<std::complex<double>> traction(n_dim + 1);
188
189 // Get FSI parameter
190 const double q_value = q();
191
192 outfile << "ZONE\n";
193
194 // Set the value of n_intpt
195 unsigned n_intpt = integral_pt()->nweight();
196
197 // Loop over the integration points
198 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
199 {
200 Vector<double> s_int(n_dim - 1);
201 s_int[0] = integral_pt()->knot(ipt, 0);
202
203 // Get the outer unit normal
204 Vector<double> interpolated_normal(n_dim);
205 outer_unit_normal(ipt, interpolated_normal);
206
207 // Boundary coordinate
208 Vector<double> zeta(1);
209 interpolated_zeta(s_int, zeta);
210
211
212 // Get bulk element for potential
213 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
214 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0, ipt));
216 std::complex<double> u_helmholtz =
217 ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
218
219 // Traction: Pressure is proportional to POSITIVE potential
220 ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
221 for (unsigned i = 0; i < n_dim; i++)
222 {
223 traction[i] = -q_value * interpolated_normal[i] * u_helmholtz;
224 }
225
226 outfile << ext_el_pt->interpolated_x(s_ext, 0) << " "
227 << ext_el_pt->interpolated_x(s_ext, 1) << " " << traction[0]
228 << " " << traction[1] << " " << interpolated_normal[0] << " "
229 << interpolated_normal[1] << " " << u_helmholtz << " "
230 << zeta[0] << std::endl;
231 }
232 }
233
234 /// C_style output function
235 void output(FILE* file_pt)
236 {
238 }
239
240 /// C-style output function
241 void output(FILE* file_pt, const unsigned& n_plot)
242 {
244 }
245 };
246
247 /// ////////////////////////////////////////////////////////////////////
248 /// ////////////////////////////////////////////////////////////////////
249 /// ////////////////////////////////////////////////////////////////////
250
251
252 //=================================================================
253 /// Static default value for the ratio of stress scales
254 /// used in the fluid and solid equations (default is 1.0)
255 //=================================================================
256 template<class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
257 double
258 FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
259 ELASTICITY_BULK_ELEMENT,
260 HELMHOLTZ_BULK_ELEMENT>::Default_Q_Value = 1.0;
261
262
263 //=====================================================================
264 /// Return the residuals
265 //=====================================================================
266 template<class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
267 void FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
268 ELASTICITY_BULK_ELEMENT,
269 HELMHOLTZ_BULK_ELEMENT>::
270 fill_in_contribution_to_residuals_helmholtz_traction(
271 Vector<double>& residuals)
272 {
273 // Find out how many nodes there are
274 unsigned n_node = nnode();
275
276#ifdef PARANOID
277 // Find out how many positional dofs there are
278 unsigned n_position_type = this->nnodal_position_type();
279 if (n_position_type != 1)
280 {
281 throw OomphLibError("LinearElasticity is not yet implemented for more "
282 "than one position type.",
283 OOMPH_CURRENT_FUNCTION,
284 OOMPH_EXCEPTION_LOCATION);
285 }
286#endif
287
288 // Find out the dimension of the node
289 const unsigned n_dim = this->nodal_dimension();
290
291 // Cache the nodal indices at which the displacement components are stored
292 std::vector<std::complex<unsigned>> u_nodal_index(n_dim + 1);
293 for (unsigned i = 0; i < n_dim + 1; i++)
294 {
295 u_nodal_index[i] =
296 this
297 ->U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
298 [i];
299 }
300
301 // Integer to hold the local equation number
302 int local_eqn = 0;
303
304 // Set up memory for the shape functions
305 Shape psi(n_node);
306 DShape dpsids(n_node, n_dim - 1);
307
308 // Get FSI parameter
309 const double q_value = q();
310
311 // Storage for traction
312 Vector<std::complex<double>> traction(n_dim + 1);
313
314 // Set the value of n_intpt
315 unsigned n_intpt = integral_pt()->nweight();
316
317 // Loop over the integration points
318 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
319 {
320 // Get the integral weight
321 double w = integral_pt()->weight(ipt);
322
323 // Only need to call the local derivatives
324 dshape_local_at_knot(ipt, psi, dpsids);
325
326 // Calculate the coordinates
327 Vector<double> interpolated_x(n_dim, 0.0);
328
329 // Also calculate the surface tangent vectors
330 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
331
332 // Calculate displacements and derivatives
333 for (unsigned l = 0; l < n_node; l++)
334 {
335 // Loop over directions
336 for (unsigned i = 0; i < n_dim; i++)
337 {
338 // Calculate the Eulerian coords
339 const double x_local = nodal_position(l, i);
340 interpolated_x[i] += x_local * psi(l);
341
342 // Loop over LOCAL derivative directions, to calculate the tangent(s)
343 for (unsigned j = 0; j < n_dim - 1; j++)
344 {
345 interpolated_A(j, i) += x_local * dpsids(l, j);
346 }
347 }
348 }
349
350 // First component for the residuals
351 double r = interpolated_x[0];
352
353 // Now find the local metric tensor from the tangent Vectors
354 DenseMatrix<double> A(n_dim - 1);
355 for (unsigned i = 0; i < n_dim - 1; i++)
356 {
357 for (unsigned j = 0; j < n_dim - 1; j++)
358 {
359 // Initialise surface metric tensor to zero
360 A(i, j) = 0.0;
361
362 // Take the dot product
363 for (unsigned k = 0; k < n_dim; k++)
364 {
365 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
366 }
367 }
368 }
369
370 // Get the outer unit normal
371 Vector<double> interpolated_normal(n_dim);
372 outer_unit_normal(ipt, interpolated_normal);
373
374 // Find the determinant of the metric tensor
375 double Adet = 0.0;
376 switch (n_dim)
377 {
378 case 2:
379 Adet = A(0, 0);
380 break;
381 case 3:
382 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
383 break;
384 default:
385 throw OomphLibError(
386 "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
387 "TimeHarmonicLinElastLoadedByPressureElement::fill_in_contribution_"
388 "to_residuals()",
389 OOMPH_EXCEPTION_LOCATION);
390 }
391
392 // Premultiply the weights and the square-root of the determinant of
393 // the metric tensor
394 double W = w * sqrt(Adet);
395
396 // Get bulk element for potential
397 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
398 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0, ipt));
399 Vector<double> s_ext(external_element_local_coord(0, ipt));
400
401 // Traction: Pressure is proportional to POSITIVE potential
402 std::complex<double> u_helmholtz =
403 ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
404 for (unsigned i = 0; i < n_dim; i++)
405 {
406 traction[i] = -q_value * interpolated_normal[i] * u_helmholtz;
407 }
408 // Loop over the test functions, nodes of the element
409 for (unsigned l = 0; l < n_node; l++)
410 {
411 // Loop over the displacement components
412 for (unsigned i = 0; i < n_dim + 1; i++)
413 {
414 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].real());
415 /*IF it's not a boundary condition*/
416 if (local_eqn >= 0)
417 {
418 // Add the loading terms to the residuals
419 residuals[local_eqn] -= traction[i].real() * psi(l) * r * W;
420 }
421
422 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].imag());
423 /*IF it's not a boundary condition*/
424 if (local_eqn >= 0)
425 {
426 // Add the loading terms to the residuals
427 residuals[local_eqn] -= traction[i].imag() * psi(l) * r * W;
428 }
429
430 } // End of if not boundary condition
431 } // End of loop over shape functions
432 } // End of loop over integration points
433 }
434
435
436 /// //////////////////////////////////////////////////////////////////////
437 /// //////////////////////////////////////////////////////////////////////
438 /// //////////////////////////////////////////////////////////////////////
439
440
441 //======================================================================
442 /// A class for elements that allow the imposition of an
443 /// prescribed flux (determined from the normal displacements of an
444 /// adjacent linearly elastic solid. Normal derivative for displacement
445 /// potential is given by normal displacement of adjacent solid multiplies
446 /// by FSI parameter (q = k^2 B/E).
447 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
448 /// policy class.
449 //======================================================================
450 template<class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
452 : public virtual FaceGeometry<HELMHOLTZ_BULK_ELEMENT>,
453 public virtual FaceElement,
454 public virtual ElementWithExternalElement
455 {
456 public:
457 /// Constructor, takes the pointer to the "bulk" element and the
458 /// face index identifying the face to which the element is attached.
460 FiniteElement* const& bulk_el_pt, const int& face_index);
461
462 /// Broken copy constructor
465 dummy) = delete;
466
467 /// Broken assignment operator
468 // Commented out broken assignment operator because this can lead to a
469 // conflict warning when used in the virtual inheritence hierarchy.
470 // Essentially the compiler doesn't realise that two separate
471 // implementations of the broken function are the same and so, quite
472 // rightly, it shouts.
473 /*void operator=(
474 const FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement&) =
475 delete;*/
476
477
478 /// Add the element's contribution to its residual vector
480 {
481 // Call the generic residuals function with flag set to 0
482 // using a dummy matrix argument
485 }
486
487
488 /// Add the element's contribution to its residual vector and its
489 /// Jacobian matrix
491 DenseMatrix<double>& jacobian)
492 {
493 // Call the generic routine with the flag set to 1
495 residuals, jacobian, 1);
496
497 // Derivatives w.r.t. external data
499 }
500
501 /// Output function
502 void output(std::ostream& outfile)
503 {
504 // Dummy
505 unsigned nplot = 0;
506 output(outfile, nplot);
507 }
508
509 /// Output function: flux etc at Gauss points; n_plot is ignored.
510 void output(std::ostream& outfile, const unsigned& n_plot)
511 {
512 outfile << "ZONE\n";
513
514 // Get the value of Nintpt
515 const unsigned n_intpt = integral_pt()->nweight();
516
517 // Set the Vector to hold local coordinates
518 Vector<double> s_int(Dim - 1);
519
520 // Loop over the integration points
521 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
522 {
523 // Assign values of s
524 for (unsigned i = 0; i < (Dim - 1); i++)
525 {
526 s_int[i] = integral_pt()->knot(ipt, i);
527 }
528
529 // Get the unit normal
530 Vector<double> interpolated_normal(Dim);
531 outer_unit_normal(ipt, interpolated_normal);
532
533 Vector<double> zeta(1);
534 interpolated_zeta(s_int, zeta);
535 // Get displacements
536 ELASTICITY_BULK_ELEMENT* ext_el_pt =
537 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(external_element_pt(0, ipt));
540 ext_el_pt
541 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
542 s_ext, displ);
543
544 // Convert into flux BC: This takes the dot product of the
545 // actual displacement with the flux element's own outer
546 // unit normal so the plus sign is OK.
547 std::complex<double> flux = (displ[0] * interpolated_normal[0] +
548 displ[1] * interpolated_normal[1]);
549
550 // Output
551 outfile << ext_el_pt->interpolated_x(s_ext, 0) << " "
552 << ext_el_pt->interpolated_x(s_ext, 1) << " "
553 << flux * interpolated_normal[0] << " "
554 << flux * interpolated_normal[1] << " "
555 << interpolated_normal[0] << " " << interpolated_normal[1]
556 << " " << flux << " " << zeta[0] << std::endl;
557 }
558 }
559
560
561 /// C-style output function
562 void output(FILE* file_pt)
563 {
565 }
566
567 /// C-style output function
568 void output(FILE* file_pt, const unsigned& n_plot)
569 {
571 }
572
573
574 protected:
575 /// Function to compute the shape and test functions and to return
576 /// the Jacobian of mapping between local and global (Eulerian)
577 /// coordinates
578 inline double shape_and_test(const Vector<double>& s,
579 Shape& psi,
580 Shape& test) const
581 {
582 // Find number of nodes
583 unsigned n_node = nnode();
584
585 // Get the shape functions
586 shape(s, psi);
587
588 // Set the test functions to be the same as the shape functions
589 for (unsigned i = 0; i < n_node; i++)
590 {
591 test[i] = psi[i];
592 }
593
594 // Return the value of the jacobian
595 return J_eulerian(s);
596 }
597
598
599 /// Function to compute the shape and test functions and to return
600 /// the Jacobian of mapping between local and global (Eulerian)
601 /// coordinates
602 inline double shape_and_test_at_knot(const unsigned& ipt,
603 Shape& psi,
604 Shape& test) const
605 {
606 // Find number of nodes
607 unsigned n_node = nnode();
608
609 // Get the shape functions
610 shape_at_knot(ipt, psi);
611
612 // Set the test functions to be the same as the shape functions
613 for (unsigned i = 0; i < n_node; i++)
614 {
615 test[i] = psi[i];
616 }
617
618 // Return the value of the jacobian
619 return J_eulerian_at_knot(ipt);
620 }
621
622
623 private:
624 /// Add the element's contribution to its residual vector.
625 /// flag=1(or 0): do (or don't) compute the contribution to the
626 /// Jacobian as well.
628 Vector<double>& residuals,
629 DenseMatrix<double>& jacobian,
630 const unsigned& flag);
631
632 /// The spatial dimension of the problem
633 unsigned Dim;
634
635 /// The index at which the unknown is stored at the nodes
637 };
638
639 /// ///////////////////////////////////////////////////////////////////
640 /// ///////////////////////////////////////////////////////////////////
641 /// ///////////////////////////////////////////////////////////////////
642
643
644 //===========================================================================
645 /// Constructor, takes the pointer to the "bulk" element, and the
646 /// face index that identifies the face of the bulk element to which
647 /// this face element is to be attached.
648 //===========================================================================
649 template<class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
651 HELMHOLTZ_BULK_ELEMENT,
652 ELASTICITY_BULK_ELEMENT>::
653 FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(
654 FiniteElement* const& bulk_el_pt, const int& face_index)
655 : FaceGeometry<HELMHOLTZ_BULK_ELEMENT>(), FaceElement()
656 {
657 // Let the bulk element build the FaceElement, i.e. setup the pointers
658 // to its nodes (by referring to the appropriate nodes in the bulk
659 // element), etc.
660 bulk_el_pt->build_face_element(face_index, this);
661
662#ifdef PARANOID
663 {
664 // Check that the element is not a refineable 3d element
665 HELMHOLTZ_BULK_ELEMENT* elem_pt =
666 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(bulk_el_pt);
667 // If it's three-d
668 if (elem_pt->dim() == 3)
669 {
670 // Is it refineable
671 RefineableElement* ref_el_pt =
672 dynamic_cast<RefineableElement*>(elem_pt);
673 if (ref_el_pt != 0)
674 {
675 if (this->has_hanging_nodes())
676 {
677 throw OomphLibError("This flux element will not work correctly if "
678 "nodes are hanging\n",
679 OOMPH_CURRENT_FUNCTION,
680 OOMPH_EXCEPTION_LOCATION);
681 }
682 }
683 }
684 }
685#endif
686
687 // Set source element storage: one interaction with an external element
688 // that provides the displacement of the adjacent linear elasticity
689 // element
690 this->set_ninteraction(1);
691
692
693 // Extract the dimension of the problem from the dimension of
694 // the first node
695 Dim = bulk_el_pt->nodal_dimension(); // this->node_pt(0)->ndim();// icidav
696
697 // Set up U_index_helmholtz_displacement. Initialise to zero, which
698 // probably won't change in most cases, oh well, the price we
699 // pay for generality
700 U_index_helmholtz_from_displacement = std::complex<unsigned>(0, 0);
701
702 // Cast to the appropriate HelmholtzEquation so that we can
703 // find the index at which the variable is stored
704 // We assume that the dimension of the full problem is the same
705 // as the dimension of the node, if this is not the case you will have
706 // to write custom elements, sorry
707
709 dynamic_cast<FourierDecomposedHelmholtzEquations*>(bulk_el_pt);
710 // If the cast has failed die
711 if (eqn_pt == 0)
712 {
713 std::string error_string =
714 "Bulk element must inherit from HelmholtzEquations.";
715 error_string +=
716 "Nodes are three dimensional, but cannot cast the bulk element to\n";
717 error_string += "HelmholtzEquations<3>\n.";
718 error_string +=
719 "If you desire this functionality, you must implement it yourself\n";
720
721 throw OomphLibError(
722 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
723 }
724 else
725 {
726 // Read the index from the (cast) bulk element.
729 }
730 }
731
732
733 //===========================================================================
734 /// Helper function to compute the element's residual vector and
735 /// the Jacobian matrix.
736 //===========================================================================
737 template<class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
739 HELMHOLTZ_BULK_ELEMENT,
740 ELASTICITY_BULK_ELEMENT>::
741 fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(
742 Vector<double>& residuals,
743 DenseMatrix<double>& jacobian,
744 const unsigned& flag)
745 {
746 // Find out how many nodes there are
747 const unsigned n_node = nnode();
748
749 // Set up memory for the shape and test functions
750 Shape psif(n_node), testf(n_node);
751
752 // Set the value of Nintpt
753 const unsigned n_intpt = integral_pt()->nweight();
754
755 // Set the Vector to hold local coordinates
756 Vector<double> s(Dim - 1);
757
758 // Integers to hold the local equation and unknown numbers
759 int local_eqn = 0;
760
761 // Locally cache the index at which the variable is stored
762 const std::complex<unsigned> u_index_helmholtz =
763 U_index_helmholtz_from_displacement;
764
765 // Loop over the integration points
766 //--------------------------------
767 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
768 {
769 // Assign values of s
770 for (unsigned i = 0; i < (Dim - 1); i++)
771 {
772 s[i] = integral_pt()->knot(ipt, i);
773 }
774
775 // Get the integral weight
776 double w = integral_pt()->weight(ipt);
777
778 // Find the shape and test functions and return the Jacobian
779 // of the mapping
780 double J = shape_and_test(s, psif, testf);
781
782 // Premultiply the weights and the Jacobian
783 double W = w * J;
784
785 // Need to find position to feed into flux function, initialise to zero
786 Vector<double> interpolated_x(Dim, 0.0);
787
788 // Calculate velocities and derivatives
789 for (unsigned l = 0; l < n_node; l++)
790 {
791 // Loop over velocity components
792 for (unsigned i = 0; i < Dim; i++)
793 {
794 interpolated_x[i] += nodal_position(l, i) * psif[l];
795 }
796 }
797
798 // first component
799 double r = interpolated_x[0];
800 // Get the outer unit normal
801 Vector<double> interpolated_normal(Dim);
802 outer_unit_normal(ipt, interpolated_normal);
803
804 // Get displacements
805 ELASTICITY_BULK_ELEMENT* ext_el_pt =
806 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(external_element_pt(0, ipt));
807 Vector<double> s_ext(external_element_local_coord(0, ipt));
809 ext_el_pt
810 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
811 s_ext, displ);
812
813 // Convert into flux BC: This takes the dot product of the
814 // actual displacement with the flux element's own outer
815 // unit normal so the plus sign is OK.
816 std::complex<double> flux =
817 (displ[0] * interpolated_normal[0] + displ[1] * interpolated_normal[1]);
818
819 // Now add to the appropriate equations
820
821 // Loop over the test functions
822 for (unsigned l = 0; l < n_node; l++)
823 {
824 local_eqn = nodal_local_eqn(l, u_index_helmholtz.real());
825 /*IF it's not a boundary condition*/
826 if (local_eqn >= 0)
827 {
828 // Add the prescribed flux terms
829 residuals[local_eqn] -= flux.real() * testf[l] * r * W;
830
831 // Imposed traction doesn't depend upon the solution,
832 // --> the Jacobian is always zero, so no Jacobian
833 // terms are required
834 }
835
836 local_eqn = nodal_local_eqn(l, u_index_helmholtz.imag());
837 /*IF it's not a boundary condition*/
838 if (local_eqn >= 0)
839 {
840 // Add the prescribed flux terms
841 residuals[local_eqn] -= flux.imag() * testf[l] * r * W;
842
843 // Imposed traction doesn't depend upon the solution,
844 // --> the Jacobian is always zero, so no Jacobian
845 // terms are required
846 }
847 }
848 }
849 }
850} // namespace oomph
851
852
853#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
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 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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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 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
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Broken assignment operator.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Constructor, takes the pointer to the "bulk" element and the face index identifying the face to which...
void fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(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...
FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(const FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement &dummy)=delete
Broken copy constructor.
std::complex< unsigned > U_index_helmholtz_from_displacement
The index at which the unknown is stored at the nodes.
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.
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(std::ostream &outfile, const unsigned &n_plot)
Output function: flux etc at Gauss points; n_plot is ignored.
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 ...
A class for elements that allow the imposition of an applied traction in the equations of time-harmon...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations.
Vector< std::complex< unsigned > > U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
Index at which the i-th displacement component is stored.
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and solid equations (default is...
void fill_in_contribution_to_residuals_helmholtz_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.
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations....
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
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.
An OomphLibError object which should be thrown when an run-time error is encountered....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
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.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
//////////////////////////////////////////////////////////////////// ////////////////////////////////...