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