pml_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_GEN_HH_TIME_HARMONIC_LIN_ELAST_HEADER
27#define OOMPH_GEN_HH_TIME_HARMONIC_LIN_ELAST_HEADER
28
29
30// Oomph-lib includes
31#include "generic.h"
32#include "pml_helmholtz.h"
33#include "time_harmonic_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
40 /// PMLHelmholtz potential (interpreted as a displacement potential
41 /// for the fluid in a quasi-steady, linearised FSI problem.)
42 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
43 /// class 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 PMLHelmholtz bulk element that provides
113 // the pressure
114 this->set_ninteraction(1);
115
116 // Find the dimension of the problem
117 unsigned n_dim = element_pt->nodal_dimension();
118
119 // Find the index at which the displacement unknowns are stored
120 ELASTICITY_BULK_ELEMENT* cast_element_pt =
121 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(element_pt);
123 n_dim);
124 for (unsigned i = 0; i < n_dim; i++)
125 {
127 cast_element_pt->u_index_time_harmonic_linear_elasticity(i);
128 }
129 }
130
131
132 /// Return the residuals
134 {
136 }
137
138
139 /// Fill in contribution from Jacobian
141 DenseMatrix<double>& jacobian)
142 {
143 // Call the residuals
145
146 // Derivatives w.r.t. external data
148 }
149
150 /// Return the ratio of the stress scales used to non-dimensionalise
151 /// the fluid and elasticity equations. E.g.
152 /// \f$ Q = (\omega a)^2 \rho/E \f$, i.e. the
153 /// ratio between the inertial fluid stress and the solid's elastic
154 /// modulus E.
155 const double& q() const
156 {
157 return *Q_pt;
158 }
159
160 /// Return a pointer the ratio of stress scales used to
161 /// non-dimensionalise the fluid and solid equations.
162 double*& q_pt()
163 {
164 return Q_pt;
165 }
166
167
168 /// Output function
169 void output(std::ostream& outfile)
170 {
171 /// Dummy
172 unsigned nplot = 0;
173 output(outfile, nplot);
174 }
175
176 /// Output function: Plot traction etc at Gauss points
177 /// nplot is ignored.
178 void output(std::ostream& outfile, const unsigned& n_plot)
179 {
180 // Dimension
181 unsigned n_dim = this->nodal_dimension();
182
183 // Storage for traction
184 Vector<std::complex<double>> traction(n_dim);
185
186 // Get FSI parameter
187 const double q_value = q();
188
189 outfile << "ZONE\n";
190
191 // Set the value of n_intpt
192 unsigned n_intpt = integral_pt()->nweight();
193
194 // Loop over the integration points
195 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
196 {
197 Vector<double> s_int(n_dim - 1);
198 for (unsigned i = 0; i < n_dim - 1; i++)
199 {
200 s_int[i] = integral_pt()->knot(ipt, i);
201 }
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 // Get external element for potential
212 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
213 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0, ipt));
215 std::complex<double> u_helmholtz =
216 ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
217
218 // Traction: Pressure is proportional to POSITIVE potential
219 ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
220 traction[0] = -q_value * interpolated_normal[0] * u_helmholtz;
221 traction[1] = -q_value * interpolated_normal[1] * u_helmholtz;
222
223 outfile << ext_el_pt->interpolated_x(s_ext, 0) << " "
224 << ext_el_pt->interpolated_x(s_ext, 1) << " "
225 << traction[0].real() << " " << traction[1].real() << " "
226 << traction[0].imag() << " " << traction[1].imag() << " "
227 << interpolated_normal[0] << " " << interpolated_normal[1]
228 << " " << u_helmholtz.real() << " " << u_helmholtz.imag() << " "
229 << interpolated_x(s_int, 0) << " " << interpolated_x(s_int, 1)
230 << " "
231 << sqrt(pow(ext_el_pt->interpolated_x(s_ext, 0) -
232 interpolated_x(s_int, 0),
233 2) +
234 pow(ext_el_pt->interpolated_x(s_ext, 1) -
235 interpolated_x(s_int, 1),
236 2))
237 << " " << zeta[0] << std::endl;
238 }
239 }
240
241
242 /// C_style output function
243 void output(FILE* file_pt)
244 {
246 }
247
248 /// C-style output function
249 void output(FILE* file_pt, const unsigned& n_plot)
250 {
252 }
253
254 /// Compute the global_flux_contribution through the template
255 /// elasticity elements : we compute u.n
257 {
258 // Dummy output file
259 std::ofstream outfile;
261 }
262
263 /// Compute the element's contribution to the integral of
264 /// the flux over the artificial boundary. Also output the
265 /// flux in the specified output file if it's open.
267 std::ofstream& outfile)
268 {
269 // pointer to the corresponding elasticity bulk element
270 ELASTICITY_BULK_ELEMENT* bulk_elem_pt =
271 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(this->bulk_element_pt());
272
273 // get the dim of the bulk and local nodes
274 const unsigned bulk_dim = bulk_elem_pt->dim();
275 const unsigned local_dim = this->dim();
276
277 // Set up memory for the outer unit normal
278 Vector<double> unit_normal(bulk_dim);
279
280 // Set the value of n_intpt
281 const unsigned n_intpt = integral_pt()->nweight();
282
283 // Set the Vector to hold local coordinates
284 Vector<double> s(local_dim);
285 std::complex<double> flux(0.0, 0.0);
286 Vector<std::complex<double>> u(bulk_dim);
287
288 // Output?
289 if (outfile.is_open())
290 {
291 outfile << "ZONE\n";
292 }
293
294 // Loop over the integration points
295 //--------------------------------
296 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
297 {
298 // Assign values of s
299 for (unsigned i = 0; i < local_dim; i++)
300 {
301 s[i] = integral_pt()->knot(ipt, i);
302 }
303 // get the outer_unit_ext vector
304 this->outer_unit_normal(s, unit_normal);
305
306 // Get the integral weight
307 double w = integral_pt()->weight(ipt);
308
309 // Get jacobian of mapping
310 double J = J_eulerian(s);
311
312 // Premultiply the weights and the Jacobian
313 double W = w * J;
314
315 // Get local coordinates in bulk element by copy construction
317
318 // Get the displacement from the bulk element
319 bulk_elem_pt->interpolated_u_time_harmonic_linear_elasticity(s_bulk, u);
320
321 // Compute normal displacement u.n
322 std::complex<double> u_dot_n(0.0, 0.0);
323 for (unsigned i = 0; i < bulk_dim; i++)
324 {
325 u_dot_n += u[i] * unit_normal[i];
326 }
327
328 // Output?
329 if (outfile.is_open())
330 {
331 Vector<double> x(bulk_dim);
332 interpolated_x(s, x);
333 outfile << x[0] << " " << x[1] << " " << u_dot_n.real() << " "
334 << u_dot_n.imag() << "\n";
335 }
336
337 // ...add to integral
338 flux += u_dot_n * W;
339 }
340
341 return flux;
342 }
343
344
345 /// Compute the global_flux_contribution through the helmholtz
346 /// elements : we compute dphidn.n
348 {
349 // Dummy output file
350 std::ofstream outfile;
352 }
353
354 /// Compute the element's contribution to the integral of
355 /// the flux over the artificial boundary. Also output the
356 /// flux in the specified output file if it's open.
358 std::ofstream& outfile)
359 {
360 // Get the dim of this element
361 const unsigned local_dim = this->dim();
362
363 // Set the value of n_intpt
364 const unsigned n_intpt = integral_pt()->nweight();
365
366 // Set the Vector to hold local coordinates
367 Vector<double> s(local_dim);
368 std::complex<double> flux(0.0, 0.0);
369
370 // Output?
371 if (outfile.is_open())
372 {
373 outfile << "ZONE\n";
374 }
375
376 // Loop over the integration points
377 //--------------------------------
378 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
379 {
380 // Assign values of s
381 for (unsigned i = 0; i < local_dim; i++)
382 {
383 s[i] = integral_pt()->knot(ipt, i);
384 }
385
386 // Get local coordinates in bulk element by copy construction
388
389 // Get external element for potential
390 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
391 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0, ipt));
392
393 // number of nodes in ext element
394 unsigned nnode_ext = ext_el_pt->nnode();
395
396 // get the dim of the external nodes
397 const unsigned ext_dim = ext_el_pt->dim();
398
399 // Set up memory for the external shape function
400 Shape psi_ext(nnode_ext);
401 DShape dpsi_ext_dx(nnode_ext, ext_dim);
402
403 // Call the derivatives of the shape functions
404 // in the external element -- must do this via s because this point
405 // is not an integration point the bulk element!
406 (void)ext_el_pt->dshape_eulerian(s_bulk, psi_ext, dpsi_ext_dx);
407
408 // Set up memory for the outer unit normal
409 Vector<double> unit_normal(ext_dim);
410
411 // get the outer_unit_ext vector
412 this->outer_unit_normal(s, unit_normal);
413
414 // Get the integral weight
415 double w = integral_pt()->weight(ipt);
416
417 // Get jacobian of mapping
418 double J = J_eulerian(s);
419
420 // Premultiply the weights and the Jacobian
421 double W = w * J;
422
423 // Global gradient
424 Vector<std::complex<double>> interpolated_dphidx(
425 ext_dim, std::complex<double>(0.0, 0.0));
426 for (unsigned l = 0; l < nnode_ext; l++)
427 {
428 // Get the nodal value of the helmholtz unknown
429 const std::complex<double> phi_value(
430 ext_el_pt->nodal_value(l, ext_el_pt->u_index_helmholtz().real()),
431 ext_el_pt->nodal_value(l, ext_el_pt->u_index_helmholtz().imag()));
432
433 // Loop over directions
434 for (unsigned i = 0; i < ext_dim; i++)
435 {
436 interpolated_dphidx[i] += phi_value * dpsi_ext_dx(l, i);
437 }
438 }
439
440 // define dphi_dn
441 std::complex<double> dphi_dn(0.0, 0.0);
442 for (unsigned i = 0; i < ext_dim; i++)
443 {
444 dphi_dn += interpolated_dphidx[i] * unit_normal[i];
445 }
446
447
448#ifdef PARANOID
449 double max_permitted_discrepancy = 1.0e-10;
450 double diff = sqrt(
451 pow(ext_el_pt->interpolated_x(s_bulk, 0) - interpolated_x(s, 0), 2) +
452 pow(ext_el_pt->interpolated_x(s_bulk, 1) - interpolated_x(s, 1), 2));
453 if (diff > max_permitted_discrepancy)
454 {
455 std::ostringstream error_stream;
456 error_stream
457 << "Position computed by external el and face element differ by "
458 << diff << " which is more than the acceptable tolerance "
459 << max_permitted_discrepancy << std::endl;
460 throw OomphLibError(error_stream.str(),
461 OOMPH_CURRENT_FUNCTION,
462 OOMPH_EXCEPTION_LOCATION);
463 }
464#endif
465
466 // Output?
467 if (outfile.is_open())
468 {
469 Vector<double> x(ext_dim);
470 interpolated_x(s, x);
471 outfile << x[0] << " " << x[1] << " " << dphi_dn.real() << " "
472 << dphi_dn.imag() << " "
473 << "\n";
474 }
475
476 // ...add to integral
477 flux += dphi_dn * W;
478 }
479
480 return flux;
481 }
482 };
483
484
485 /// ////////////////////////////////////////////////////////////////////
486 /// ////////////////////////////////////////////////////////////////////
487 /// ////////////////////////////////////////////////////////////////////
488
489
490 //=================================================================
491 /// Static default value for the ragoogletio of stress scales
492 /// used in the fluid and solid equations (default is 1.0)
493 //=================================================================
494 template<class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
495 double TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement<
496 ELASTICITY_BULK_ELEMENT,
497 HELMHOLTZ_BULK_ELEMENT>::Default_Q_Value = 1.0;
498
499
500 //=====================================================================
501 /// Return the residuals
502 //=====================================================================
503 template<class ELASTICITY_BULK_ELEMENT, class HELMHOLTZ_BULK_ELEMENT>
504 void TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement<
505 ELASTICITY_BULK_ELEMENT,
506 HELMHOLTZ_BULK_ELEMENT>::
507 fill_in_contribution_to_residuals_helmholtz_traction(
508 Vector<double>& residuals)
509 {
510 // Find out how many nodes there are
511 unsigned n_node = nnode();
512
513#ifdef PARANOID
514 // Find out how many positional dofs there are
515 unsigned n_position_type = this->nnodal_position_type();
516 if (n_position_type != 1)
517 {
518 throw OomphLibError("LinearElasticity is not yet implemented for more "
519 "than one position type.",
520 OOMPH_CURRENT_FUNCTION,
521 OOMPH_EXCEPTION_LOCATION);
522 }
523#endif
524
525 // Find out the dimension of the node
526 const unsigned n_dim = this->nodal_dimension();
527
528 // Cache the nodal indices at which the displacement components are stored
529 std::complex<unsigned> u_nodal_index[n_dim];
530 for (unsigned i = 0; i < n_dim; i++)
531 {
532 u_nodal_index[i] =
533 this->U_index_time_harmonic_linear_elasticity_helmholtz_traction[i];
534 }
535
536 // Integer to hold the local equation number
537 int local_eqn = 0;
538
539 // Set up memory for the shape functions
540 Shape psi(n_node);
541 DShape dpsids(n_node, n_dim - 1);
542
543 // Get FSI parameter
544 const double q_value = q();
545
546 // Storage for traction
547 Vector<std::complex<double>> traction(2);
548
549 // Set the value of n_intpt
550 unsigned n_intpt = integral_pt()->nweight();
551
552 // Loop over the integration points
553 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
554 {
555 // Get the integral weight
556 double w = integral_pt()->weight(ipt);
557
558 // Only need to call the local derivatives
559 dshape_local_at_knot(ipt, psi, dpsids);
560
561 // Calculate the coordinates
562 Vector<double> interpolated_x(n_dim, 0.0);
563
564 // Also calculate the surface tangent vectors
565 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
566
567 // Calculate displacements and derivatives
568 for (unsigned l = 0; l < n_node; l++)
569 {
570 // Loop over directions
571 for (unsigned i = 0; i < n_dim; i++)
572 {
573 // Calculate the Eulerian coords
574 const double x_local = nodal_position(l, i);
575 interpolated_x[i] += x_local * psi(l);
576
577 // Loop over LOCAL derivative directions, to calculate the tangent(s)
578 for (unsigned j = 0; j < n_dim - 1; j++)
579 {
580 interpolated_A(j, i) += x_local * dpsids(l, j);
581 }
582 }
583 }
584
585 // Now find the local metric tensor from the tangent Vectors
586 DenseMatrix<double> A(n_dim - 1);
587 for (unsigned i = 0; i < n_dim - 1; i++)
588 {
589 for (unsigned j = 0; j < n_dim - 1; j++)
590 {
591 // Initialise surface metric tensor to zero
592 A(i, j) = 0.0;
593
594 // Take the dot product
595 for (unsigned k = 0; k < n_dim; k++)
596 {
597 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
598 }
599 }
600 }
601
602 // Get the outer unit normal
603 Vector<double> interpolated_normal(n_dim);
604 outer_unit_normal(ipt, interpolated_normal);
605
606 // Find the determinant of the metric tensor
607 double Adet = 0.0;
608 switch (n_dim)
609 {
610 case 2:
611 Adet = A(0, 0);
612 break;
613 case 3:
614 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
615 break;
616 default:
617 throw OomphLibError(
618 "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
619 OOMPH_CURRENT_FUNCTION,
620 OOMPH_EXCEPTION_LOCATION);
621 }
622
623 // Premultiply the weights and the square-root of the determinant of
624 // the metric tensor
625 double W = w * sqrt(Adet);
626
627 // Get bulk element for potential
628 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
629 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(external_element_pt(0, ipt));
630 Vector<double> s_ext(external_element_local_coord(0, ipt));
631
632 // Traction: Pressure is proportional to POSITIVE potential
633 std::complex<double> u_helmholtz =
634 ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
635 traction[0] = -q_value * interpolated_normal[0] * u_helmholtz;
636 traction[1] = -q_value * interpolated_normal[1] * u_helmholtz;
637
638 // Loop over the test functions, nodes of the element
639 for (unsigned l = 0; l < n_node; l++)
640 {
641 // Loop over the displacement components
642 for (unsigned i = 0; i < n_dim; i++)
643 {
644 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].real());
645 /*IF it's not a boundary condition*/
646 if (local_eqn >= 0)
647 {
648 // Add the loading terms to the residuals
649 residuals[local_eqn] -= traction[i].real() * psi(l) * W;
650 }
651
652 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].imag());
653 /*IF it's not a boundary condition*/
654 if (local_eqn >= 0)
655 {
656 // Add the loading terms to the residuals
657 residuals[local_eqn] -= traction[i].imag() * psi(l) * W;
658 }
659
660 } // End of if not boundary condition
661 } // End of loop over shape functions
662 } // End of loop over integration points
663 }
664
665
666 /// ////////////////////////////////////////////////////////////////////
667 /// ////////////////////////////////////////////////////////////////////
668 /// ////////////////////////////////////////////////////////////////////
669
670
671 //======================================================================
672 /// A class for elements that allow the imposition of an
673 /// prescribed flux (determined from the normal displacements of an
674 /// adjacent linearly elastic solid. Normal derivative for displacement
675 /// potential is given by normal displacement of adjacent solid multiplies
676 /// by FSI parameter (q = k^2 B/E).
677 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
678 /// policy class.
679 //======================================================================
680 template<class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
682 : public virtual FaceGeometry<HELMHOLTZ_BULK_ELEMENT>,
683 public virtual FaceElement,
684 public virtual ElementWithExternalElement
685 {
686 public:
687 /// Constructor, takes the pointer to the "bulk" element and the
688 /// face index identifying the face to which the element is attached.
690 FiniteElement* const& bulk_el_pt, const int& face_index);
691
692 /// Broken copy constructor
695
696 /// Broken assignment operator
697 // Commented out broken assignment operator because this can lead to a
698 // conflict warning when used in the virtual inheritence hierarchy.
699 // Essentially the compiler doesn't realise that two separate
700 // implementations of the broken function are the same and so, quite
701 // rightly, it shouts.
702 /*void operator=(const
703 PMLHelmholtzFluxFromNormalDisplacementBCElement&) = delete;*/
704
705
706 /// Add the element's contribution to its residual vector
708 {
709 // Call the generic residuals function with flag set to 0
710 // using a dummy matrix argument
713 }
714
715
716 /// Add the element's contribution to its residual vector and its
717 /// Jacobian matrix
719 DenseMatrix<double>& jacobian)
720 {
721 // Call the generic routine with the flag set to 1
723 residuals, jacobian, 1);
724
725 // Derivatives w.r.t. external data
727 }
728
729 /// Output function
730 void output(std::ostream& outfile)
731 {
732 // Dummy
733 unsigned nplot = 0;
734 output(outfile, nplot);
735 }
736
737 /// Output function: flux etc at Gauss points; n_plot is ignored.
738 void output(std::ostream& outfile, const unsigned& n_plot)
739 {
740 outfile << "ZONE\n";
741
742 // Get the value of Nintpt
743 const unsigned n_intpt = integral_pt()->nweight();
744
745 // Set the Vector to hold local coordinates
746 Vector<double> s_int(Dim - 1);
747
748 // Loop over the integration points
749 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
750 {
751 // Assign values of s
752 for (unsigned i = 0; i < (Dim - 1); i++)
753 {
754 s_int[i] = integral_pt()->knot(ipt, i);
755 }
756
757 // Get the unit normal
758 Vector<double> interpolated_normal(Dim);
759 outer_unit_normal(ipt, interpolated_normal);
760
761 Vector<double> zeta(1);
762 interpolated_zeta(s_int, zeta);
763
764 // Get displacements
765 ELASTICITY_BULK_ELEMENT* ext_el_pt =
766 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(external_element_pt(0, ipt));
769 ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext, displ);
770
771 // Convert into flux BC: This takes the dot product of the
772 // actual displacement with the flux element's own outer
773 // unit normal so the plus sign is OK.
774 std::complex<double> flux = (displ[0] * interpolated_normal[0] +
775 displ[1] * interpolated_normal[1]);
776
777 // Output
778 outfile << ext_el_pt->interpolated_x(s_ext, 0) << " "
779 << ext_el_pt->interpolated_x(s_ext, 1) << " "
780 << flux.real() * interpolated_normal[0] << " "
781 << flux.real() * interpolated_normal[1] << " "
782 << flux.imag() * interpolated_normal[0] << " "
783 << flux.imag() * interpolated_normal[1] << " "
784 << interpolated_normal[0] << " " << interpolated_normal[1]
785 << " " << flux.real() << " " << flux.imag() << " "
786 << interpolated_x(s_int, 0) << " " << interpolated_x(s_int, 1)
787 << " "
788 << sqrt(pow(ext_el_pt->interpolated_x(s_ext, 0) -
789 interpolated_x(s_int, 0),
790 2) +
791 pow(ext_el_pt->interpolated_x(s_ext, 1) -
792 interpolated_x(s_int, 1),
793 2))
794 << " " << zeta[0] << std::endl;
795 }
796 }
797
798
799 /// C-style output function
800 void output(FILE* file_pt)
801 {
803 }
804
805 /// C-style output function
806 void output(FILE* file_pt, const unsigned& n_plot)
807 {
809 }
810
811
812 protected:
813 /// Function to compute the shape and test functions and to return
814 /// the Jacobian of mapping between local and global (Eulerian)
815 /// coordinates
816 inline double shape_and_test(const Vector<double>& s,
817 Shape& psi,
818 Shape& test) const
819 {
820 // Find number of nodes
821 unsigned n_node = nnode();
822
823 // Get the shape functions
824 shape(s, psi);
825
826 // Set the test functions to be the same as the shape functions
827 for (unsigned i = 0; i < n_node; i++)
828 {
829 test[i] = psi[i];
830 }
831
832 // Return the value of the jacobian
833 return J_eulerian(s);
834 }
835
836
837 /// Function to compute the shape and test functions and to return
838 /// the Jacobian of mapping between local and global (Eulerian)
839 /// coordinates
840 inline double shape_and_test_at_knot(const unsigned& ipt,
841 Shape& psi,
842 Shape& test) const
843 {
844 // Find number of nodes
845 unsigned n_node = nnode();
846
847 // Get the shape functions
848 shape_at_knot(ipt, psi);
849
850 // Set the test functions to be the same as the shape functions
851 for (unsigned i = 0; i < n_node; i++)
852 {
853 test[i] = psi[i];
854 }
855
856 // Return the value of the jacobian
857 return J_eulerian_at_knot(ipt);
858 }
859
860
861 private:
862 /// Add the element's contribution to its residual vector.
863 /// flag=1(or 0): do (or don't) compute the contribution to the
864 /// Jacobian as well.
866 Vector<double>& residuals,
867 DenseMatrix<double>& jacobian,
868 const unsigned& flag);
869
870 /// The spatial dimension of the problem
871 unsigned Dim;
872
873 /// The index at which the unknown is stored at the nodes
875 };
876
877 /// ///////////////////////////////////////////////////////////////////
878 /// ///////////////////////////////////////////////////////////////////
879 /// ///////////////////////////////////////////////////////////////////
880
881
882 //===========================================================================
883 /// Constructor, takes the pointer to the "bulk" element, and the
884 /// face index that identifies the face of the bulk element to which
885 /// this face element is to be attached.
886 //===========================================================================
887 template<class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
889 ELASTICITY_BULK_ELEMENT>::
890 PMLHelmholtzFluxFromNormalDisplacementBCElement(
891 FiniteElement* const& bulk_el_pt, const int& face_index)
892 : FaceGeometry<HELMHOLTZ_BULK_ELEMENT>(), FaceElement()
893 {
894 // Let the bulk element build the FaceElement, i.e. setup the pointers
895 // to its nodes (by referring to the appropriate nodes in the bulk
896 // element), etc.
897 bulk_el_pt->build_face_element(face_index, this);
898
899#ifdef PARANOID
900 {
901 // Check that the element is not a refineable 3d element
902 HELMHOLTZ_BULK_ELEMENT* elem_pt =
903 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*>(bulk_el_pt);
904 // If it's three-d
905 if (elem_pt->dim() == 3)
906 {
907 // Is it refineable
908 RefineableElement* ref_el_pt =
909 dynamic_cast<RefineableElement*>(elem_pt);
910 if (ref_el_pt != 0)
911 {
912 if (this->has_hanging_nodes())
913 {
914 throw OomphLibError("This flux element will not work correctly if "
915 "nodes are hanging\n",
916 OOMPH_CURRENT_FUNCTION,
917 OOMPH_EXCEPTION_LOCATION);
918 }
919 }
920 }
921 }
922#endif
923
924 // Set source element storage: one interaction with an external element
925 // that provides the displacement of the adjacent linear elasticity
926 // element
927 this->set_ninteraction(1);
928
929 // Extract the dimension of the problem from the dimension of
930 // the first node
931 Dim = this->node_pt(0)->ndim();
932
933 // Set up U_index_helmholtz_displacement. Initialise to zero, which
934 // probably won't change in most cases, oh well, the price we
935 // pay for generality
936 U_index_helmholtz_from_displacement = std::complex<unsigned>(0, 0);
937
938 // Cast to the appropriate PMLHelmholtzEquation so that we can
939 // find the index at which the variable is stored
940 // We assume that the dimension of the full problem is the same
941 // as the dimension of the node, if this is not the case you will have
942 // to write custom elements, sorry
943 switch (Dim)
944 {
945 // One dimensional problem
946 case 1:
947 {
949 dynamic_cast<PMLHelmholtzEquations<1>*>(bulk_el_pt);
950 // If the cast has failed die
951 if (eqn_pt == 0)
952 {
953 std::string error_string =
954 "Bulk element must inherit from PMLHelmholtzEquations.";
955 error_string +=
956 "Nodes are one dimensional, but cannot cast the bulk element to\n";
957 error_string += "PMLHelmholtzEquations<1>\n.";
958 error_string += "If you desire this functionality, you must "
959 "implement it yourself\n";
960
961 throw OomphLibError(
962 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
963 }
964 // Otherwise read out the value
965 else
966 {
967 // Read the index from the (cast) bulk element
969 }
970 }
971 break;
972
973 // Two dimensional problem
974 case 2:
975 {
977 dynamic_cast<PMLHelmholtzEquations<2>*>(bulk_el_pt);
978 // If the cast has failed die
979 if (eqn_pt == 0)
980 {
981 std::string error_string =
982 "Bulk element must inherit from PMLHelmholtzEquations.";
983 error_string +=
984 "Nodes are two dimensional, but cannot cast the bulk element to\n";
985 error_string += "PMLHelmholtzEquations<2>\n.";
986 error_string += "If you desire this functionality, you must "
987 "implement it yourself\n";
988
989 throw OomphLibError(
990 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
991 }
992 else
993 {
994 // Read the index from the (cast) bulk element.
996 }
997 }
998 break;
999
1000 // Three dimensional problem
1001 case 3:
1002 {
1003 PMLHelmholtzEquations<3>* eqn_pt =
1004 dynamic_cast<PMLHelmholtzEquations<3>*>(bulk_el_pt);
1005 // If the cast has failed die
1006 if (eqn_pt == 0)
1007 {
1008 std::string error_string =
1009 "Bulk element must inherit from PMLHelmholtzEquations.";
1010 error_string += "Nodes are three dimensional, but cannot cast the "
1011 "bulk element to\n";
1012 error_string += "PMLHelmholtzEquations<3>\n.";
1013 error_string += "If you desire this functionality, you must "
1014 "implement it yourself\n";
1015
1016 throw OomphLibError(
1017 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1018 }
1019 else
1020 {
1021 // Read the index from the (cast) bulk element.
1023 }
1024 }
1025 break;
1026
1027 // Any other case is an error
1028 default:
1029 std::ostringstream error_stream;
1030 error_stream << "Dimension of node is " << Dim
1031 << ". It should be 1,2, or 3!" << std::endl;
1032
1033 throw OomphLibError(
1034 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1035 break;
1036 }
1037 }
1038
1039
1040 //===========================================================================
1041 /// Helper function to compute the element's residual vector and
1042 /// the Jacobian matrix.
1043 //===========================================================================
1044 template<class HELMHOLTZ_BULK_ELEMENT, class ELASTICITY_BULK_ELEMENT>
1046 HELMHOLTZ_BULK_ELEMENT,
1047 ELASTICITY_BULK_ELEMENT>::
1048 fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(
1049 Vector<double>& residuals,
1050 DenseMatrix<double>& jacobian,
1051 const unsigned& flag)
1052 {
1053 // Find out how many nodes there are
1054 const unsigned n_node = nnode();
1055
1056 // Set up memory for the shape and test functions
1057 Shape psif(n_node), testf(n_node);
1058
1059 // Set the value of Nintpt
1060 const unsigned n_intpt = integral_pt()->nweight();
1061
1062 // Set the Vector to hold local coordinates
1063 Vector<double> s(Dim - 1);
1064
1065 // Integers to hold the local equation and unknown numbers
1066 int local_eqn = 0;
1067
1068 // Locally cache the index at which the variable is stored
1069 const std::complex<unsigned> u_index_helmholtz =
1070 U_index_helmholtz_from_displacement;
1071
1072 // Loop over the integration points
1073 //--------------------------------
1074 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1075 {
1076 // Assign values of s
1077 for (unsigned i = 0; i < (Dim - 1); i++)
1078 {
1079 s[i] = integral_pt()->knot(ipt, i);
1080 }
1081
1082 // Get the integral weight
1083 double w = integral_pt()->weight(ipt);
1084
1085 // Find the shape and test functions and return the Jacobian
1086 // of the mapping
1087 double J = shape_and_test(s, psif, testf);
1088
1089 // Premultiply the weights and the Jacobian
1090 double W = w * J;
1091
1092 // Need to find position to feed into flux function, initialise to zero
1093 Vector<double> interpolated_x(Dim, 0.0);
1094
1095 // Calculate velocities and derivatives
1096 for (unsigned l = 0; l < n_node; l++)
1097 {
1098 // Loop over velocity components
1099 for (unsigned i = 0; i < Dim; i++)
1100 {
1101 interpolated_x[i] += nodal_position(l, i) * psif[l];
1102 }
1103 }
1104
1105 // Get the outer unit normal
1106 Vector<double> interpolated_normal(2);
1107 outer_unit_normal(ipt, interpolated_normal);
1108
1109 // Get displacements
1110 ELASTICITY_BULK_ELEMENT* ext_el_pt =
1111 dynamic_cast<ELASTICITY_BULK_ELEMENT*>(external_element_pt(0, ipt));
1112 Vector<double> s_ext(external_element_local_coord(0, ipt));
1114 ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext, displ);
1115
1116 // Convert into flux BC: This takes the dot product of the
1117 // actual displacement with the flux element's own outer
1118 // unit normal so the plus sign is OK.
1119 std::complex<double> flux =
1120 (displ[0] * interpolated_normal[0] + displ[1] * interpolated_normal[1]);
1121
1122 // Now add to the appropriate equations
1123
1124 // Loop over the test functions
1125 for (unsigned l = 0; l < n_node; l++)
1126 {
1127 local_eqn = nodal_local_eqn(l, u_index_helmholtz.real());
1128 /*IF it's not a boundary condition*/
1129 if (local_eqn >= 0)
1130 {
1131 // Add the prescribed flux terms
1132 residuals[local_eqn] -= flux.real() * testf[l] * W;
1133
1134 // Imposed traction doesn't depend upon the solution,
1135 // --> the Jacobian is always zero, so no Jacobian
1136 // terms are required
1137 }
1138
1139 local_eqn = nodal_local_eqn(l, u_index_helmholtz.imag());
1140 /*IF it's not a boundary condition*/
1141 if (local_eqn >= 0)
1142 {
1143 // Add the prescribed flux terms
1144 residuals[local_eqn] -= flux.imag() * testf[l] * W;
1145
1146 // Imposed traction doesn't depend upon the solution,
1147 // --> the Jacobian is always zero, so no Jacobian
1148 // terms are required
1149 }
1150 }
1151 }
1152 }
1153} // namespace oomph
1154
1155#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
Vector< double > local_coordinate_in_bulk(const Vector< double > &s) const
Return vector of local coordinates in bulk element, given the local coordinates in this FaceElement.
Definition: elements.cc:6353
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s....
Definition: elements.cc:5242
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 dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
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
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
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....
A class for all isoparametric elements that solve the Helmholtz equations with pml capabilities....
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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...
PMLHelmholtzFluxFromNormalDisplacementBCElement(const PMLHelmholtzFluxFromNormalDisplacementBCElement &dummy)=delete
Broken copy constructor.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
std::complex< unsigned > U_index_helmholtz_from_displacement
The index at which the unknown is stored at the nodes.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Broken assignment operator.
PMLHelmholtzFluxFromNormalDisplacementBCElement(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...
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.
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...
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...
std::complex< double > global_flux_contribution_from_solid(std::ofstream &outfile)
Compute the element's contribution to the integral of the flux over the artificial boundary....
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
std::complex< double > global_flux_contribution_from_helmholtz()
Compute the global_flux_contribution through the helmholtz elements : we compute dphidn....
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....
TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
std::complex< double > global_flux_contribution_from_helmholtz(std::ofstream &outfile)
Compute the element's contribution to the integral of the flux over the artificial boundary....
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
std::complex< double > global_flux_contribution_from_solid()
Compute the global_flux_contribution through the template elasticity elements : we compute u....
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...