helmholtz_bc_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for elements that are used to apply Sommerfeld
27// boundary conditions to the Helmholtz equations
28#ifndef OOMPH_HELMHOLTZ_BC_ELEMENTS_HEADER
29#define OOMPH_HELMHOLTZ_BC_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36#include "math.h"
37#include <complex>
38
39// Get the Bessel functions
40#include "oomph_crbond_bessel.h"
41
42
43namespace oomph
44{
45 /// ////////////////////////////////////////////////////////////////////
46 /// ////////////////////////////////////////////////////////////////////
47 // Collection of the Bessel functions used in the Helmholtz problem
48 /// ////////////////////////////////////////////////////////////////////
49 /// ////////////////////////////////////////////////////////////////////
50
51 //====================================================================
52 /// Namespace to provide Hankel function of the first kind
53 /// and various orders -- needed for Helmholtz computations.
54 //====================================================================
55 namespace Hankel_functions_for_helmholtz_problem
56 {
57 //====================================================================
58 /// Compute Hankel function of the first kind of orders 0...n and
59 /// its derivates at coordinate x. The function returns the vector
60 /// then its derivative.
61 //====================================================================
62 void Hankel_first(const unsigned& n,
63 const double& x,
64 Vector<std::complex<double>>& h,
65 Vector<std::complex<double>>& hp)
66 {
67 int n_actual = 0;
68 Vector<double> jn(n + 1), yn(n + 1), jnp(n + 1), ynp(n + 1);
69 CRBond_Bessel::bessjyna(
70 int(n), x, n_actual, &jn[0], &yn[0], &jnp[0], &ynp[0]);
71#ifdef PARANOID
72 if (n_actual != int(n))
73 {
74 std::ostringstream error_stream;
75 error_stream << "CRBond_Bessel::bessjyna() only computed " << n_actual
76 << " rather than " << n << " Bessel functions.\n";
77 throw OomphLibError(
78 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
79 }
80#endif
81 for (unsigned i = 0; i < n; i++)
82 {
83 h[i] = std::complex<double>(jn[i], yn[i]);
84 hp[i] = std::complex<double>(jnp[i], ynp[i]);
85 }
86 } // End of Hankel_first
87
88 //====================================================================
89 /// Compute Hankel function of the first kind of orders 0...n and
90 /// its derivates at coordinate x. The function returns the vector
91 /// then its derivative (complex version). This functionality is only
92 /// required in the computation of the solution for the complex-
93 /// shifted Laplacian preconditioner.
94 //====================================================================
95 void CHankel_first(const unsigned& n,
96 const std::complex<double>& x,
97 Vector<std::complex<double>>& h,
98 Vector<std::complex<double>>& hp)
99 {
100 // Set the highest order actually calculated
101 int n_actual = 0;
102
103 // Create a vector for the Bessel function of the 1st kind
105
106 // Create a vector for the Bessel function of the 2nd kind
108
109 // Create a vector for the Bessel function (1st kind) derivative
110 Vector<std::complex<double>> jnp(n + 1);
111
112 // Create a vector for the Bessel function (2nd kind) derivative
113 Vector<std::complex<double>> ynp(n + 1);
114
115 // Call the (complex) Bessel function to calculate the solution
116 CRBond_Bessel::cbessjyna(
117 int(n), x, n_actual, &jn[0], &yn[0], &jnp[0], &ynp[0]);
118
119#ifdef PARANOID
120 // Tell the user if they tried to calculate higher order terms
121 if (n_actual != int(n))
122 {
123 // Create an output stream
124 std::ostringstream error_message_stream;
125
126 // Create the error message
127 error_message_stream << "CRBond_Bessel::cbessjyna() only computed "
128 << n_actual << " rather than " << n
129 << " Bessel functions.\n";
130
131 // Throw the error message
132 throw OomphLibError(error_message_stream.str(),
133 OOMPH_CURRENT_FUNCTION,
134 OOMPH_EXCEPTION_LOCATION);
135 }
136#endif
137
138 // Loop over the number of terms requested (only the first entry
139 // has *actually* been calculated)
140 for (unsigned i = 0; i < n; i++)
141 {
142 // Set the entries in the Hankel function vector
143 h[i] = std::complex<double>(jn[i].real() - yn[i].imag(),
144 jn[i].imag() + yn[i].real());
145
146 // Set the entries in the Hankel function derivative vector
147 hp[i] = std::complex<double>(jnp[i].real() - ynp[i].imag(),
148 jnp[i].imag() + ynp[i].real());
149 }
150 } // End of Hankel_first
151 } // namespace Hankel_functions_for_helmholtz_problem
152
153
154 /// //////////////////////////////////////////////////////////////////
155 /// //////////////////////////////////////////////////////////////////
156 /// //////////////////////////////////////////////////////////////////
157
158
159 //======================================================================
160 /// A class for elements that allow the approximation of the
161 /// Sommerfeld radiation BC.
162 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
163 /// policy class.
164 //======================================================================
165 template<class ELEMENT>
166 class HelmholtzBCElementBase : public virtual FaceGeometry<ELEMENT>,
167 public virtual FaceElement
168 {
169 public:
170 /// Constructor, takes the pointer to the "bulk" element and the
171 /// index of the face to which the element is attached.
172 HelmholtzBCElementBase(FiniteElement* const& bulk_el_pt,
173 const int& face_index);
174
175 /// Broken empty constructor
177 {
178 throw OomphLibError(
179 "Don't call empty constructor for HelmholtzBCElementBase",
180 OOMPH_CURRENT_FUNCTION,
181 OOMPH_EXCEPTION_LOCATION);
182 }
183
184 /// Broken copy constructor
186
187 /// Broken assignment operator
188 // Commented out broken assignment operator because this can lead to a
189 // conflict warning when used in the virtual inheritence hierarchy.
190 // Essentially the compiler doesn't realise that two separate
191 // implementations of the broken function are the same and so, quite
192 // rightly, it shouts.
193 /*void operator=(const HelmholtzBCElementBase&) = delete;*/
194
195
196 /// Specify the value of nodal zeta from the face geometry
197 /// The "global" intrinsic coordinate of the element when
198 /// viewed as part of a geometric object should be given by
199 /// the FaceElement representation, by default (needed to break
200 /// indeterminacy if bulk element is SolidElement)
201 double zeta_nodal(const unsigned& n,
202 const unsigned& k,
203 const unsigned& i) const
204 {
205 return FaceElement::zeta_nodal(n, k, i);
206 }
207
208
209 /// Output function -- forward to broken version in FiniteElement
210 /// until somebody decides what exactly they want to plot here...
211 void output(std::ostream& outfile)
212 {
213 FiniteElement::output(outfile);
214 }
215
216 /// Output function -- forward to broken version in FiniteElement
217 /// until somebody decides what exactly they want to plot here...
218 void output(std::ostream& outfile, const unsigned& n_plot)
219 {
220 FiniteElement::output(outfile, n_plot);
221 }
222
223 /// C-style output function -- forward to broken version in FiniteElement
224 /// until somebody decides what exactly they want to plot here...
225 void output(FILE* file_pt)
226 {
227 FiniteElement::output(file_pt);
228 }
229
230 /// C-style output function -- forward to broken version in
231 /// FiniteElement until somebody decides what exactly they want to plot
232 /// here...
233 void output(FILE* file_pt, const unsigned& n_plot)
234 {
235 FiniteElement::output(file_pt, n_plot);
236 }
237
238 /// Return the index at which the real/imag unknown value
239 /// is stored.
240 virtual inline std::complex<unsigned> u_index_helmholtz() const
241 {
242 return std::complex<unsigned>(U_index_helmholtz.real(),
243 U_index_helmholtz.imag());
244 }
245
246 /// Compute the element's contribution to the time-averaged
247 /// radiated power over the artificial boundary
249 {
250 // Dummy output file
251 std::ofstream outfile;
252 return global_power_contribution(outfile);
253 }
254
255 /// Compute the element's contribution to the time-averaged
256 /// radiated power over the artificial boundary. Also output the
257 /// power density as a fct of the polar angle in the specified
258 /// output file if it's open.
259 double global_power_contribution(std::ofstream& outfile)
260 {
261 // pointer to the corresponding bulk element
262 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
263
264 // Number of nodes in bulk element
265 unsigned nnode_bulk = bulk_elem_pt->nnode();
266 const unsigned n_node_local = nnode();
267
268 // get the dim of the bulk and local nodes
269 const unsigned bulk_dim = bulk_elem_pt->dim();
270 const unsigned local_dim = this->dim();
271
272 // Set up memory for the shape and test functions
273 Shape psi(n_node_local);
274
275 // Set up memory for the shape functions
276 Shape psi_bulk(nnode_bulk);
277 DShape dpsi_bulk_dx(nnode_bulk, bulk_dim);
278
279 // Set up memory for the outer unit normal
280 Vector<double> unit_normal(bulk_dim);
281
282 // Set the value of n_intpt
283 const unsigned n_intpt = integral_pt()->nweight();
284
285 // Set the Vector to hold local coordinates
286 Vector<double> s(local_dim);
287 double power = 0.0;
288
289 // Output?
290 if (outfile.is_open())
291 {
292 outfile << "ZONE\n";
293 }
294
295 // Loop over the integration points
296 //--------------------------------
297 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
298 {
299 // Assign values of s
300 for (unsigned i = 0; i < local_dim; i++)
301 {
302 s[i] = integral_pt()->knot(ipt, i);
303 }
304 // get the outer_unit_ext vector
305 this->outer_unit_normal(s, unit_normal);
306
307 // Get the integral weight
308 double w = integral_pt()->weight(ipt);
309
310 // Get jacobian of mapping
311 double J = J_eulerian(s);
312
313 // Premultiply the weights and the Jacobian
314 double W = w * J;
315
316 // Get local coordinates in bulk element by copy construction
318
319 // Call the derivatives of the shape functions
320 // in the bulk -- must do this via s because this point
321 // is not an integration point the bulk element!
322 (void)bulk_elem_pt->dshape_eulerian(s_bulk, psi_bulk, dpsi_bulk_dx);
323 this->shape(s, psi);
324
325 // Derivs of Eulerian coordinates w.r.t. local coordinates
326 std::complex<double> dphi_dn(0.0, 0.0);
327 Vector<std::complex<double>> interpolated_dphidx(bulk_dim);
328 std::complex<double> interpolated_phi(0.0, 0.0);
329 Vector<double> x(bulk_dim);
330
331 // Calculate function value and derivatives:
332 //-----------------------------------------
333 // Loop over nodes
334 for (unsigned l = 0; l < nnode_bulk; l++)
335 {
336 // Get the nodal value of the helmholtz unknown
337 const std::complex<double> phi_value(
338 bulk_elem_pt->nodal_value(l,
339 bulk_elem_pt->u_index_helmholtz().real()),
340 bulk_elem_pt->nodal_value(
341 l, bulk_elem_pt->u_index_helmholtz().imag()));
342
343 // Loop over directions
344 for (unsigned i = 0; i < bulk_dim; i++)
345 {
346 interpolated_dphidx[i] += phi_value * dpsi_bulk_dx(l, i);
347 }
348 } // End of loop over the bulk_nodes
349
350 for (unsigned l = 0; l < n_node_local; l++)
351 {
352 // Get the nodal value of the helmholtz unknown
353 const std::complex<double> phi_value(
354 nodal_value(l, u_index_helmholtz().real()),
355 nodal_value(l, u_index_helmholtz().imag()));
356
357 interpolated_phi += phi_value * psi(l);
358 }
359
360 // define dphi_dn
361 for (unsigned i = 0; i < bulk_dim; i++)
362 {
363 dphi_dn += interpolated_dphidx[i] * unit_normal[i];
364 }
365
366 // Power density
367 double integrand = 0.5 * (interpolated_phi.real() * dphi_dn.imag() -
368 interpolated_phi.imag() * dphi_dn.real());
369
370 // Output?
371 if (outfile.is_open())
372 {
373 interpolated_x(s, x);
374 double phi = atan2(x[1], x[0]);
375 outfile << x[0] << " " << x[1] << " " << phi << " " << integrand
376 << "\n";
377 }
378
379 // ...add to integral
380 power += integrand * W;
381 }
382
383 return power;
384 }
385
386
387 /// Compute element's contribution to Fourier components of the
388 /// solution -- length of vector indicates number of terms to be computed.
390 Vector<std::complex<double>>& a_coeff_pos,
391 Vector<std::complex<double>>& a_coeff_neg)
392 {
393#ifdef PARANOID
394 if (a_coeff_pos.size() != a_coeff_neg.size())
395 {
396 std::ostringstream error_stream;
397 error_stream << "a_coeff_pos and a_coeff_neg must have "
398 << "the same size. \n";
399 throw OomphLibError(
400 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
401 }
402#endif
403
404 // define the imaginary number
405 const std::complex<double> I(0.0, 1.0);
406
407 // Find out how many nodes there are
408 const unsigned n_node = this->nnode();
409
410 // Set up memory for the shape functions
411 Shape psi(n_node);
412 DShape dpsi(n_node, 1);
413
414 // Set the value of n_intpt
415 const unsigned n_intpt = this->integral_pt()->nweight();
416
417 // Set the Vector to hold local coordinates
418 Vector<double> s(this->Dim - 1);
419
420 // Initialise
421 unsigned n = a_coeff_pos.size();
422 for (unsigned i = 0; i < n; i++)
423 {
424 a_coeff_pos[i] = std::complex<double>(0.0, 0.0);
425 a_coeff_neg[i] = std::complex<double>(0.0, 0.0);
426 }
427
428 // Loop over the integration points
429 //--------------------------------
430 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
431 {
432 // Assign values of s
433 for (unsigned i = 0; i < (this->Dim - 1); i++)
434 {
435 s[i] = this->integral_pt()->knot(ipt, i);
436 }
437
438 // Get the integral weight
439 double w = this->integral_pt()->weight(ipt);
440
441 // Get the shape functions
442 this->dshape_local(s, psi, dpsi);
443
444 // Eulerian coordinates at Gauss point
446
447 // Derivs of Eulerian coordinates w.r.t. local coordinates
448 Vector<double> interpolated_dxds(this->Dim);
449 std::complex<double> interpolated_u(0.0, 0.0);
450
451 // Assemble x and its derivs
452 for (unsigned l = 0; l < n_node; l++)
453 {
454 // Loop over directions
455 for (unsigned i = 0; i < this->Dim; i++)
456 {
457 interpolated_x[i] += this->nodal_position(l, i) * psi[l];
458 interpolated_dxds[i] += this->nodal_position(l, i) * dpsi(l, 0);
459 }
460
461 // Get the nodal value of the helmholtz unknown
462 const std::complex<double> u_value(
463 this->nodal_value(l, this->U_index_helmholtz.real()),
464 this->nodal_value(l, this->U_index_helmholtz.imag()));
465
466 // Add to the interpolated value
467 interpolated_u += u_value * psi(l);
468 } // End of loop over the nodes
469
470 // calculate the integral
471 //-----------------------
472
473 // Get polar angle
474 double phi = atan2(interpolated_x[1], interpolated_x[0]);
475
476 // define dphi_ds=(-yx'+y'x)/(x^2+y^2)
477 double denom = (interpolated_x[0] * interpolated_x[0]) +
479 double nom = -interpolated_dxds[1] * interpolated_x[0] +
480 interpolated_dxds[0] * interpolated_x[1];
481 double dphi_ds = std::fabs(nom / denom);
482
483 // Positive coefficients
484 for (unsigned i = 0; i < n; i++)
485 {
486 a_coeff_pos[i] +=
487 interpolated_u * exp(-I * phi * double(i)) * dphi_ds * w;
488 }
489 // Negative coefficients
490 for (unsigned i = 1; i < n; i++)
491 {
492 a_coeff_neg[i] +=
493 interpolated_u * exp(I * phi * double(i)) * dphi_ds * w;
494 }
495
496 } // End of loop over integration points
497 }
498
499
500 protected:
501 /// Function to compute the shape and test functions and to return
502 /// the Jacobian of mapping between local and global (Eulerian)
503 /// coordinates
504 inline double test_only(const Vector<double>& s, Shape& test) const
505 {
506 // Get the shape functions
507 shape(s, test);
508
509 // Return the value of the jacobian
510 return J_eulerian(s);
511 }
512
513 /// Function to compute the shape, test functions and derivates
514 /// and to return
515 /// the Jacobian of mapping between local and global (Eulerian)
516 /// coordinates
518 Shape& psi,
519 Shape& test,
520 DShape& dpsi_ds,
521 DShape& dtest_ds) const
522 {
523 // Find number of nodes
524 unsigned n_node = nnode();
525
526 // Get the shape functions
527 dshape_local(s, psi, dpsi_ds);
528
529 // Set the test functions to be the same as the shape functions
530 for (unsigned i = 0; i < n_node; i++)
531 {
532 for (unsigned j = 0; j < (Dim - 1); j++)
533 {
534 test[i] = psi[i];
535 dtest_ds(i, j) = dpsi_ds(i, j);
536 }
537 }
538 // Return the value of the jacobian
539 return J_eulerian(s);
540 }
541
542 /// The index at which the real and imag part of the unknown is
543 /// stored at the nodes
544 std::complex<unsigned> U_index_helmholtz;
545
546 /// The spatial dimension of the problem
547 unsigned Dim;
548 };
549
550 /// ///////////////////////////////////////////////////////////////////
551 /// ///////////////////////////////////////////////////////////////////
552 /// ///////////////////////////////////////////////////////////////////
553
554
555 /// =================================================================
556 /// Mesh for DtN boundary condition elements -- provides
557 /// functionality to apply Sommerfeld radiation condtion
558 /// via DtN BC
559 /// =================================================================
560 template<class ELEMENT>
561 class HelmholtzDtNMesh : public virtual Mesh
562 {
563 public:
564 /// Constructor: Specify radius of outer boundary and number of
565 /// Fourier terms used in the computation of the gamma integral
566 HelmholtzDtNMesh(const double& outer_radius, const unsigned& nfourier_terms)
568 {
569 }
570
571 /// Compute and store the gamma integral at all integration
572 /// points of the constituent elements.
573 void setup_gamma();
574
575 /// Compute Fourier components of the solution -- length of
576 /// vector indicates number of terms to be computed.
577 void compute_fourier_components(Vector<std::complex<double>>& a_coeff_pos,
578 Vector<std::complex<double>>& a_coeff_neg);
579
580 /// Gamma integral evaluated at Gauss points
581 /// for specified element
583 {
584 return Gamma_at_gauss_point[el_pt];
585 }
586
587 /// Derivative of Gamma integral w.r.t global unknown, evaluated
588 /// at Gauss points for specified element
590 FiniteElement* el_pt)
591 {
592 return D_Gamma_at_gauss_point[el_pt];
593 }
594
595 /// The outer radius
596 double& outer_radius()
597 {
598 return Outer_radius;
599 }
600
601 /// Number of Fourier terms used in the computation of the
602 /// gamma integral
603 unsigned& nfourier_terms()
604 {
605 return Nfourier_terms;
606 }
607
608 private:
609 /// Outer radius
611
612 /// Nbr of Fourier terms used in the Gamma computation
614
615
616 /// Container to store the gamma integral for given Gauss point
617 /// and element
618 std::map<FiniteElement*, Vector<std::complex<double>>> Gamma_at_gauss_point;
619
620
621 /// Container to store the derivate of Gamma integral w.r.t
622 /// global unknown evaluated at Gauss points for specified element
623 std::map<FiniteElement*, Vector<std::map<unsigned, std::complex<double>>>>
625 };
626
627
628 /// //////////////////////////////////////////////////////////////////
629 /// //////////////////////////////////////////////////////////////////
630 /// //////////////////////////////////////////////////////////////////
631
632
633 //=============================================================
634 /// Absorbing BC element for approximation imposition of
635 /// Sommerfeld radiation condition
636 //==============================================================
637 template<class ELEMENT>
639 {
640 public:
641 /// Construct element from specification of bulk element and
642 /// face index.
644 const int& face_index)
645 : HelmholtzBCElementBase<ELEMENT>(bulk_el_pt, face_index)
646 {
647 // Initialise pointers
648 Outer_radius_pt = 0;
649
650 // Initialise order of absorbing boundary condition
651 ABC_order_pt = 0;
652 }
653
654 /// Pointer to order of absorbing boundary condition
655 unsigned*& abc_order_pt()
656 {
657 return ABC_order_pt;
658 }
659
660
661 /// Add the element's contribution to its residual vector
663 {
664 // Call the generic residuals function with flag set to 0
665 // using a dummy matrix argument
668 }
669
670 /// Add the element's contribution to its residual vector and its
671 /// Jacobian matrix
673 DenseMatrix<double>& jacobian)
674 {
675 // Call the generic routine with the flag set to 1
677 residuals, jacobian, 1);
678 }
679
680
681 /// Get pointer to radius of outer boundary (must be a cirle)
683 {
684 return Outer_radius_pt;
685 }
686
687 private:
688 /// Compute the element's residual vector and the (
689 /// Jacobian matrix.
690 /// Overloaded version, using the abc approximation
692 Vector<double>& residuals,
693 DenseMatrix<double>& jacobian,
694 const unsigned& flag)
695 {
696#ifdef PARANOID
697
698 if (ABC_order_pt == 0)
699 {
700 throw OomphLibError("Order of ABC hasn't been set!",
701 OOMPH_CURRENT_FUNCTION,
702 OOMPH_EXCEPTION_LOCATION);
703 }
704
705 if (this->Outer_radius_pt == 0)
706 {
707 throw OomphLibError("Pointer to outer radius hasn't been set!",
708 OOMPH_CURRENT_FUNCTION,
709 OOMPH_EXCEPTION_LOCATION);
710 }
711
712#endif
713
714 // Find out how many nodes there are
715 const unsigned n_node = this->nnode();
716
717 // Set up memory for the shape and test functions
718 Shape psi(n_node), test(n_node);
719 DShape dpsi_ds(n_node, this->Dim - 1), dtest_ds(n_node, this->Dim - 1),
720 dtest_dS(n_node, this->Dim - 1), dpsi_dS(n_node, this->Dim - 1);
721
722 // Set the value of Nintpt
723 const unsigned n_intpt = this->integral_pt()->nweight();
724
725 // Set the Vector to hold local coordinates
726 Vector<double> s(this->Dim - 1);
727
728 // Integers to hold the local equation and unknown numbers
729 int local_eqn_real = 0, local_unknown_real = 0;
730 int local_eqn_imag = 0, local_unknown_imag = 0;
731
732 // Define the problem parameters
733 double R = *Outer_radius_pt;
734 double k =
735 sqrt(dynamic_cast<ELEMENT*>(this->bulk_element_pt())->k_squared());
736
737 // Loop over the integration points
738 //--------------------------------
739 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
740 {
741 // Assign values of s
742 for (unsigned i = 0; i < (this->Dim - 1); i++)
743 {
744 s[i] = this->integral_pt()->knot(ipt, i);
745 }
746
747 // Get the integral weight
748 double w = this->integral_pt()->weight(ipt);
749
750 // Find the shape test functions and derivates; return the Jacobian
751 // of the mapping between local and global (Eulerian)
752 // coordinates
753 double J =
754 this->d_shape_and_test_local(s, psi, test, dpsi_ds, dtest_ds);
755
756 // Premultiply the weights and the Jacobian
757 double W = w * J;
758 // get the inverse of jacibian
759 double inv_J = 1 / J;
760
761 // Need to find position to feed into flux function,
762 // initialise to zero
763 std::complex<double> interpolated_u(0.0, 0.0);
764 std::complex<double> du_dS(0.0, 0.0);
765
766 // Calculate velocities and derivatives:loop over the nodes
767 for (unsigned l = 0; l < n_node; l++)
768 {
769 // Loop over real and imag part
770 // Get the nodal value of the helmholtz unknown
771 const std::complex<double> u_value(
772 this->nodal_value(l, this->U_index_helmholtz.real()),
773 this->nodal_value(l, this->U_index_helmholtz.imag()));
774
775 interpolated_u += u_value * psi[l];
776
777 du_dS += u_value * dpsi_ds(l, 0) * inv_J;
778
779 // Get the value of dtest_dS
780 dtest_dS(l, 0) = dtest_ds(l, 0) * inv_J;
781 // Get the value of dpsif_dS
782 dpsi_dS(l, 0) = dpsi_ds(l, 0) * inv_J;
783 }
784
785 // use ABC first order approximation
786 if (*ABC_order_pt == 1)
787 {
788 // Now add to the appropriate equations:use second order approximation
789 // Loop over the test functions
790 for (unsigned l = 0; l < n_node; l++)
791 {
792 local_eqn_real =
793 this->nodal_local_eqn(l, this->U_index_helmholtz.real());
794 local_eqn_imag =
795 this->nodal_local_eqn(l, this->U_index_helmholtz.imag());
796
797 // first, calculate the real part contrubution
798 //-----------------------
799 // IF it's not a boundary condition
800 if (local_eqn_real >= 0)
801 {
802 // Add the second order terms:compute manually real and imag part
803
804 residuals[local_eqn_real] += (k * interpolated_u.imag() +
805 (0.5 / R) * interpolated_u.real()) *
806 test[l] * W;
807
808 // Calculate the jacobian
809 //-----------------------
810 if (flag)
811 {
812 // Loop over the shape functions again
813 for (unsigned l2 = 0; l2 < n_node; l2++)
814 {
815 local_unknown_real =
816 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
817 local_unknown_imag =
818 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
819 // If at a non-zero degree of freedom add in the entry
820 if (local_unknown_real >= 0)
821 {
822 jacobian(local_eqn_real, local_unknown_real) +=
823 (0.5 / R) * psi[l2] * test[l] * W;
824 }
825
826 // If at a non-zero degree of freedom add in the entry
827 if (local_unknown_imag >= 0)
828 {
829 jacobian(local_eqn_real, local_unknown_imag) +=
830 k * psi[l2] * test[l] * W;
831 }
832 }
833 }
834 } // end of local_eqn_real
835
836 // second, calculate the imag part contrubution
837 //-----------------------
838 // IF it's not a boundary condition
839 if (local_eqn_imag >= 0)
840 {
841 // Add the second order terms contibution to the residual
842 residuals[local_eqn_imag] += (-k * interpolated_u.real() +
843 (0.5 / R) * interpolated_u.imag()) *
844 test[l] * W;
845
846
847 // Calculate the jacobian
848 //-----------------------
849 if (flag)
850 {
851 // Loop over the shape functions again
852 for (unsigned l2 = 0; l2 < n_node; l2++)
853 {
854 local_unknown_real =
855 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
856 local_unknown_imag =
857 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
858 // If at a non-zero degree of freedom add in the entry
859 if (local_unknown_real >= 0)
860 {
861 jacobian(local_eqn_imag, local_unknown_real) +=
862 (-k) * psi[l2] * test[l] * W;
863 }
864
865 // If at a non-zero degree of freedom add in the entry
866 if (local_unknown_imag >= 0)
867 {
868 jacobian(local_eqn_imag, local_unknown_imag) +=
869 (0.5 / R) * psi[l2] * test[l] * W;
870 }
871 }
872 }
873 }
874 } // End of loop over the nodes
875 }
876
877 //:use second order approximation
878 if (*ABC_order_pt == 2)
879 {
880 // Now add to the appropriate equations:use second order approximation
881 // Loop over the test functions
882 for (unsigned l = 0; l < n_node; l++)
883 {
884 local_eqn_real =
885 this->nodal_local_eqn(l, this->U_index_helmholtz.real());
886 local_eqn_imag =
887 this->nodal_local_eqn(l, this->U_index_helmholtz.imag());
888
889 // first, calculate the real part contrubution
890 //-----------------------
891 // IF it's not a boundary condition
892 if (local_eqn_real >= 0)
893 {
894 // Add the second order terms:compute manually real and imag part
895
896 residuals[local_eqn_real] +=
897 (k * interpolated_u.imag() +
898 (0.5 / R) * interpolated_u.real()) *
899 test[l] * W +
900 ((0.125 / (k * R * R)) * interpolated_u.imag()) * test[l] * W;
901
902 residuals[local_eqn_real] +=
903 (-0.5 / k) * du_dS.imag() * dtest_dS(l, 0) * W;
904
905 // Calculate the jacobian
906 //-----------------------
907 if (flag)
908 {
909 // Loop over the shape functions again
910 for (unsigned l2 = 0; l2 < n_node; l2++)
911 {
912 local_unknown_real =
913 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
914 local_unknown_imag =
915 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
916 // If at a non-zero degree of freedom add in the entry
917 if (local_unknown_real >= 0)
918 {
919 jacobian(local_eqn_real, local_unknown_real) +=
920 (0.5 / R) * psi[l2] * test[l] * W;
921 }
922
923 // If at a non-zero degree of freedom add in the entry
924 if (local_unknown_imag >= 0)
925 {
926 jacobian(local_eqn_real, local_unknown_imag) +=
927 k * psi[l2] * test[l] * W +
928 (0.125 / (k * R * R)) * psi[l2] * test[l] * W;
929
930 jacobian(local_eqn_real, local_unknown_imag) +=
931 (-0.5 / k) * dpsi_dS(l2, 0) * dtest_dS(l, 0) * W;
932 }
933 }
934 }
935 } // end of local_eqn_real
936
937 // second, calculate the imag part contrubution
938 //-----------------------
939 // IF it's not a boundary condition
940 if (local_eqn_imag >= 0)
941 {
942 // Add the second order terms contibution to the residual
943 residuals[local_eqn_imag] +=
944 (-k * interpolated_u.real() +
945 (0.5 / R) * interpolated_u.imag()) *
946 test[l] * W +
947 ((-0.125 / (k * R * R)) * interpolated_u.real()) * test[l] * W;
948
949 residuals[local_eqn_imag] +=
950 (0.5 / k) * du_dS.real() * dtest_dS(l, 0) * W;
951
952 // Calculate the jacobian
953 //-----------------------
954 if (flag)
955 {
956 // Loop over the shape functions again
957 for (unsigned l2 = 0; l2 < n_node; l2++)
958 {
959 local_unknown_real =
960 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
961 local_unknown_imag =
962 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
963 // If at a non-zero degree of freedom add in the entry
964 if (local_unknown_real >= 0)
965 {
966 jacobian(local_eqn_imag, local_unknown_real) +=
967 (-k) * psi[l2] * test[l] * W -
968 (0.125 / (k * R * R)) * psi[l2] * test[l] * W;
969
970 jacobian(local_eqn_imag, local_unknown_real) +=
971 (0.5 / k) * dpsi_dS(l2, 0) * dtest_dS(l, 0) * W;
972 }
973
974 // If at a non-zero degree of freedom add in the entry
975 if (local_unknown_imag >= 0)
976 {
977 jacobian(local_eqn_imag, local_unknown_imag) +=
978 (0.5 / R) * psi[l2] * test[l] * W;
979 }
980 }
981 }
982 }
983 } // End of loop over the nodes
984 }
985
986 if (*ABC_order_pt == 3)
987 {
988 // Now add to the appropriate equations:use second order approximation
989 // Loop over the test functions
990 for (unsigned l = 0; l < n_node; l++)
991 {
992 local_eqn_real =
993 this->nodal_local_eqn(l, this->U_index_helmholtz.real());
994 local_eqn_imag =
995 this->nodal_local_eqn(l, this->U_index_helmholtz.imag());
996
997 // first, calculate the real part contrubution
998 //-----------------------
999 // IF it's not a boundary condition
1000 if (local_eqn_real >= 0)
1001 {
1002 // Add the second order terms:compute manually real and imag part
1003 residuals[local_eqn_real] +=
1004 ((k * (1 + 0.125 / (k * k * R * R))) * interpolated_u.imag() +
1005 (0.5 / R - 0.125 / (k * k * R * R * R)) *
1006 interpolated_u.real()) *
1007 test[l] * W;
1008
1009 residuals[local_eqn_real] +=
1010 ((-0.5 / k) * du_dS.imag() +
1011 (0.5 / (k * k * R)) * du_dS.real()) *
1012 dtest_dS(l, 0) * W;
1013
1014 // Calculate the jacobian
1015 //-----------------------
1016 if (flag)
1017 {
1018 // Loop over the shape functions again
1019 for (unsigned l2 = 0; l2 < n_node; l2++)
1020 {
1021 local_unknown_real =
1022 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
1023 local_unknown_imag =
1024 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
1025 // If at a non-zero degree of freedom add in the entry
1026 if (local_unknown_real >= 0)
1027 {
1028 jacobian(local_eqn_real, local_unknown_real) +=
1029 (0.5 / R - 0.125 / (k * k * R * R * R)) * psi[l2] *
1030 test[l] * W;
1031
1032 jacobian(local_eqn_real, local_unknown_real) +=
1033 (0.5 / (k * k * R)) * dpsi_dS(l2, 0) * dtest_dS(l, 0) * W;
1034 }
1035
1036 // If at a non-zero degree of freedom add in the entry
1037 if (local_unknown_imag >= 0)
1038 {
1039 jacobian(local_eqn_real, local_unknown_imag) +=
1040 (k * (1 + 0.125 / (k * k * R * R))) * psi[l2] * test[l] *
1041 W;
1042
1043 jacobian(local_eqn_real, local_unknown_imag) +=
1044 (-0.5 / k) * dpsi_dS(l2, 0) * dtest_dS(l, 0) * W;
1045 }
1046 }
1047 }
1048 } // end of local_eqn_real
1049
1050 // second, calculate the imag part contrubution
1051 //-----------------------
1052 // IF it's not a boundary condition
1053 if (local_eqn_imag >= 0)
1054 {
1055 // Add the second order terms contibution to the residual
1056 residuals[local_eqn_imag] +=
1057 ((-k * (1 + 0.125 / (k * k * R * R))) * interpolated_u.real() +
1058 (0.5 / R - 0.125 / (k * k * R * R * R)) *
1059 interpolated_u.imag()) *
1060 test[l] * W;
1061
1062 residuals[local_eqn_imag] +=
1063 ((0.5 / k) * du_dS.real() +
1064 (0.5 / (k * k * R)) * du_dS.imag()) *
1065 dtest_dS(l, 0) * W;
1066
1067 // Calculate the jacobian
1068 //-----------------------
1069 if (flag)
1070 {
1071 // Loop over the shape functions again
1072 for (unsigned l2 = 0; l2 < n_node; l2++)
1073 {
1074 local_unknown_real =
1075 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
1076 local_unknown_imag =
1077 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
1078 // If at a non-zero degree of freedom add in the entry
1079 if (local_unknown_real >= 0)
1080 {
1081 jacobian(local_eqn_imag, local_unknown_real) +=
1082 (-k * (1 + 0.125 / (k * k * R * R))) * psi[l2] * test[l] *
1083 W;
1084
1085 jacobian(local_eqn_imag, local_unknown_real) +=
1086 (0.5 / k) * dpsi_dS(l2, 0) * dtest_dS(l, 0) * W;
1087 }
1088 // If at a non-zero degree of freedom add in the entry
1089 if (local_unknown_imag >= 0)
1090 {
1091 jacobian(local_eqn_imag, local_unknown_imag) +=
1092 (0.5 / R - 0.125 / (k * k * R * R * R)) * psi[l2] *
1093 test[l] * W;
1094
1095 jacobian(local_eqn_imag, local_unknown_imag) +=
1096 (0.5 / (k * k * R)) * dpsi_dS(l2, 0) * dtest_dS(l, 0) * W;
1097 }
1098 }
1099 }
1100 }
1101 } // End of loop over the nodes
1102 }
1103 } // End of loop over int_pt
1104
1105 } // End of fill_in_generic_residual_contribution_helmholtz_flux
1106
1107 private:
1108 /// Pointer to radius of outer boundary (must be a circle!)
1110
1111 /// Pointer to order of absorbing boundary condition
1112 unsigned* ABC_order_pt;
1113 };
1114
1115 //=============================================================
1116 /// FaceElement used to apply Sommerfeld radiation conditon
1117 /// via Dirichlet to Neumann map.
1118 //==============================================================
1119 template<class ELEMENT>
1121 {
1122 public:
1123 /// Construct element from specification of bulk element and
1124 /// face index.
1126 const int& face_index)
1127 : HelmholtzBCElementBase<ELEMENT>(bulk_el_pt, face_index)
1128 {
1129 }
1130
1131 /// Add the element's contribution to its residual vector
1133 {
1134 // Call the generic residuals function with flag set to 0
1135 // using a dummy matrix argument
1137 residuals, GeneralisedElement::Dummy_matrix, 0);
1138 }
1139
1140 /// Add the element's contribution to its residual vector and its
1141 /// Jacobian matrix
1143 DenseMatrix<double>& jacobian)
1144 {
1145 // Call the generic routine with the flag set to 1
1147 residuals, jacobian, 1);
1148 }
1149
1150 /// Compute the contribution of the element
1151 /// to the Gamma integral and its derivates w.r.t
1152 /// to global unknows; the function takes the wavenumber and the polar
1153 /// angle phi as input
1155 const double& phi,
1156 const int& n,
1157 std::complex<double>& gamma_con,
1158 std::map<unsigned, std::complex<double>>& d_gamma_con);
1159
1160
1161 /// Access function to mesh of all DtN boundary condition elements
1162 /// (needed to get access to gamma values)
1164 {
1166 }
1167
1168 /// Set mesh of all DtN boundary condition elements
1170 {
1171 Outer_boundary_mesh_pt = mesh_pt;
1172 }
1173
1174
1175 /// Complete the setup of additional dependencies arising
1176 /// through the far-away interaction with other nodes in
1177 /// Outer_boundary_mesh_pt.
1179 {
1180 // Create a set of all nodes
1181 std::set<Node*> node_set;
1182 unsigned nel = Outer_boundary_mesh_pt->nelement();
1183 for (unsigned e = 0; e < nel; e++)
1184 {
1185 FiniteElement* el_pt = Outer_boundary_mesh_pt->finite_element_pt(e);
1186 unsigned nnod = el_pt->nnode();
1187 for (unsigned j = 0; j < nnod; j++)
1188 {
1189 Node* nod_pt = el_pt->node_pt(j);
1190
1191 // Don't add copied nodes
1192 if (!(nod_pt->is_a_copy()))
1193 {
1194 node_set.insert(nod_pt);
1195 }
1196 }
1197 }
1198 // Now erase the current element's own nodes
1199 unsigned nnod = this->nnode();
1200 for (unsigned j = 0; j < nnod; j++)
1201 {
1202 Node* nod_pt = this->node_pt(j);
1203 node_set.erase(nod_pt);
1204
1205 // If the element's node is a copy then its "master" will
1206 // already have been added in the set above -- remove the
1207 // master to avoid double counting eqn numbers
1208 if (nod_pt->is_a_copy())
1209 {
1210 node_set.erase(nod_pt->copied_node_pt());
1211 }
1212 }
1213
1214 // Now declare these nodes to be the element's external Data
1215 for (std::set<Node*>::iterator it = node_set.begin();
1216 it != node_set.end();
1217 it++)
1218 {
1219 this->add_external_data(*it);
1220 }
1221 }
1222
1223
1224 private:
1225 /// Compute the element's residual vector
1226 /// Jacobian matrix.
1227 /// Overloaded version, using the gamma computed in the mesh
1229 Vector<double>& residuals,
1230 DenseMatrix<double>& jacobian,
1231 const unsigned& flag)
1232 {
1233 // Find out how many nodes there are
1234 const unsigned n_node = this->nnode();
1235
1236 // Set up memory for the shape and test functions
1237 Shape test(n_node);
1238
1239 // Set the value of Nintpt
1240 const unsigned n_intpt = this->integral_pt()->nweight();
1241
1242 // Set the Vector to hold local coordinates
1243 Vector<double> s(this->Dim - 1);
1244
1245 // Integers to hold the local equation and unknown numbers
1246 int local_eqn_real = 0, local_unknown_real = 0, global_eqn_real = 0,
1247 local_eqn_imag = 0, local_unknown_imag = 0, global_eqn_imag = 0;
1248 int external_global_eqn_real = 0, external_unknown_real = 0,
1249 external_global_eqn_imag = 0, external_unknown_imag = 0;
1250
1251
1252 // Get the gamma value for the current integration point
1253 // from the mesh
1255 Outer_boundary_mesh_pt->gamma_at_gauss_point(this));
1256
1258 Outer_boundary_mesh_pt->d_gamma_at_gauss_point(this));
1259
1260 // Loop over the integration points
1261 //--------------------------------
1262 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1263 {
1264 // Assign values of s
1265 for (unsigned i = 0; i < (this->Dim - 1); i++)
1266 {
1267 s[i] = this->integral_pt()->knot(ipt, i);
1268 }
1269
1270 // Get the integral weight
1271 double w = this->integral_pt()->weight(ipt);
1272
1273 // Find the shape test functions and derivates; return the Jacobian
1274 // of the mapping between local and global (Eulerian)
1275 // coordinates
1276 double J = this->test_only(s, test);
1277
1278 // Premultiply the weights and the Jacobian
1279 double W = w * J;
1280
1281 // Now add to the appropriate equations
1282 // Loop over the test functions:loop over the nodes
1283 for (unsigned l = 0; l < n_node; l++)
1284 {
1285 local_eqn_real =
1286 this->nodal_local_eqn(l, this->U_index_helmholtz.real());
1287 local_eqn_imag =
1288 this->nodal_local_eqn(l, this->U_index_helmholtz.imag());
1289
1290 // IF it's not a boundary condition
1291 if (local_eqn_real >= 0)
1292 {
1293 // Add the gamma contribution in this int_point to the res
1294 residuals[local_eqn_real] -= gamma[ipt].real() * test[l] * W;
1295
1296 // Calculate the jacobian
1297 //-----------------------
1298 if (flag)
1299 {
1300 // Loop over the shape functions again
1301 for (unsigned l2 = 0; l2 < n_node; l2++)
1302 {
1303 // Add the contribution of the local data
1304 local_unknown_real =
1305 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
1306
1307 local_unknown_imag =
1308 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
1309
1310 // If at a non-zero degree of freedom add in the entry
1311 if (local_unknown_real >= 0)
1312 {
1313 global_eqn_real = this->eqn_number(local_unknown_real);
1314
1315 // Add the first order terms contribution
1316 jacobian(local_eqn_real, local_unknown_real) -=
1317 d_gamma[ipt][global_eqn_real].real() * test[l] * W;
1318 }
1319 if (local_unknown_imag >= 0)
1320 {
1321 global_eqn_imag = this->eqn_number(local_unknown_imag);
1322
1323 // Add the first order terms contribution
1324 jacobian(local_eqn_real, local_unknown_imag) -=
1325 d_gamma[ipt][global_eqn_imag].real() * test[l] * W;
1326 }
1327 } // End of loop over nodes l2
1328
1329 // Add the contribution of the external data
1330 unsigned n_ext_data = this->nexternal_data();
1331 // Loop over the shape functions again
1332 for (unsigned l2 = 0; l2 < n_ext_data; l2++)
1333 {
1334 // Add the contribution of the local data
1335 external_unknown_real =
1336 this->external_local_eqn(l2, this->U_index_helmholtz.real());
1337
1338 external_unknown_imag =
1339 this->external_local_eqn(l2, this->U_index_helmholtz.imag());
1340
1341 // If at a non-zero degree of freedom add in the entry
1342 if (external_unknown_real >= 0)
1343 {
1344 external_global_eqn_real =
1345 this->eqn_number(external_unknown_real);
1346
1347 // Add the first order terms contribution
1348 jacobian(local_eqn_real, external_unknown_real) -=
1349 d_gamma[ipt][external_global_eqn_real].real() * test[l] * W;
1350 }
1351 if (external_unknown_imag >= 0)
1352 {
1353 external_global_eqn_imag =
1354 this->eqn_number(external_unknown_imag);
1355
1356 // Add the first order terms contribution
1357 jacobian(local_eqn_real, external_unknown_imag) -=
1358 d_gamma[ipt][external_global_eqn_imag].real() * test[l] * W;
1359 }
1360 } // End of loop over external data
1361 } // End of flag
1362 } // end of local_eqn_real
1363
1364 if (local_eqn_imag >= 0)
1365 {
1366 // Add the gamma contribution in this int_point to the res
1367 residuals[local_eqn_imag] -= gamma[ipt].imag() * test[l] * W;
1368
1369 // Calculate the jacobian
1370 //-----------------------
1371 if (flag)
1372 {
1373 // Loop over the shape functions again
1374 for (unsigned l2 = 0; l2 < n_node; l2++)
1375 {
1376 // Add the contribution of the local data
1377 local_unknown_real =
1378 this->nodal_local_eqn(l2, this->U_index_helmholtz.real());
1379
1380 local_unknown_imag =
1381 this->nodal_local_eqn(l2, this->U_index_helmholtz.imag());
1382
1383 // If at a non-zero degree of freedom add in the entry
1384 if (local_unknown_real >= 0)
1385 {
1386 global_eqn_real = this->eqn_number(local_unknown_real);
1387
1388 // Add the first order terms contribution
1389 jacobian(local_eqn_imag, local_unknown_real) -=
1390 d_gamma[ipt][global_eqn_real].imag() * test[l] * W;
1391 }
1392 if (local_unknown_imag >= 0)
1393 {
1394 global_eqn_imag = this->eqn_number(local_unknown_imag);
1395
1396 // Add the first order terms contribution
1397 jacobian(local_eqn_imag, local_unknown_imag) -=
1398 d_gamma[ipt][global_eqn_imag].imag() * test[l] * W;
1399 }
1400 } // End of loop over nodes l2
1401
1402 // Add the contribution of the external data
1403 unsigned n_ext_data = this->nexternal_data();
1404 // Loop over the shape functions again
1405 for (unsigned l2 = 0; l2 < n_ext_data; l2++)
1406 {
1407 // Add the contribution of the local data
1408 external_unknown_real =
1409 this->external_local_eqn(l2, this->U_index_helmholtz.real());
1410
1411 external_unknown_imag =
1412 this->external_local_eqn(l2, this->U_index_helmholtz.imag());
1413
1414 // If at a non-zero degree of freedom add in the entry
1415 if (external_unknown_real >= 0)
1416 {
1417 external_global_eqn_real =
1418 this->eqn_number(external_unknown_real);
1419
1420 // Add the first order terms contribution
1421 jacobian(local_eqn_imag, external_unknown_real) -=
1422 d_gamma[ipt][external_global_eqn_real].imag() * test[l] * W;
1423 }
1424 if (external_unknown_imag >= 0)
1425 {
1426 external_global_eqn_imag =
1427 this->eqn_number(external_unknown_imag);
1428
1429 // Add the first order terms contribution
1430 jacobian(local_eqn_imag, external_unknown_imag) -=
1431 d_gamma[ipt][external_global_eqn_imag].imag() * test[l] * W;
1432 }
1433 } // End of loop over external data
1434 } // End of flag
1435 } // end of local_eqn_imag
1436 } // end of llop over yhe node
1437 } // End of loop over int_pt
1438 } // End of fill_in_generic_residual_contribution_helmholtz_flux
1439
1440
1441 /// Pointer to mesh of all DtN boundary condition elements
1442 /// (needed to get access to gamma values)
1444 };
1445
1446
1447 /// /////////////////////////////////////////////////////////////
1448 /// /////////////////////////////////////////////////////////////
1449 /// /////////////////////////////////////////////////////////////
1450
1451
1452 //===========start_compute_gamma_contribution==================
1453 /// compute the contribution of the element
1454 /// to the Gamma integral and its derivates w.r.t
1455 /// to global unknows; the function takes wavenumber n
1456 /// and polar angle phi as input.
1457 //==============================================================
1458 template<class ELEMENT>
1460 const double& phi,
1461 const int& n,
1462 std::complex<double>& gamma_con,
1463 std::map<unsigned, std::complex<double>>& d_gamma_con)
1464 {
1465 // define the imaginary number
1466 const std::complex<double> I(0.0, 1.0);
1467
1468 // Find out how many nodes there are
1469 const unsigned n_node = this->nnode();
1470
1471 // Set up memory for the shape functions
1472 Shape psi(n_node);
1473 DShape dpsi(n_node, 1);
1474
1475 // initialise the variable
1476 int local_unknown_real = 0, local_unknown_imag = 0;
1477 int global_unknown_real = 0, global_unknown_imag = 0;
1478
1479 // Set the value of n_intpt
1480 const unsigned n_intpt = this->integral_pt()->nweight();
1481
1482 // Set the Vector to hold local coordinates
1483 Vector<double> s(this->Dim - 1);
1484
1485 // Initialise
1486 gamma_con = std::complex<double>(0.0, 0.0);
1487 d_gamma_con.clear();
1488
1489 // Loop over the integration points
1490 //--------------------------------
1491 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1492 {
1493 // Assign values of s
1494 for (unsigned i = 0; i < (this->Dim - 1); i++)
1495 {
1496 s[i] = this->integral_pt()->knot(ipt, i);
1497 }
1498
1499 // Get the integral weight
1500 double w = this->integral_pt()->weight(ipt);
1501
1502 // Get the shape functions
1503 this->dshape_local(s, psi, dpsi);
1504
1505 // Eulerian coordinates at Gauss point
1506 Vector<double> interpolated_x(this->Dim, 0.0);
1507
1508 // Derivs of Eulerian coordinates w.r.t. local coordinates
1509 Vector<double> interpolated_dxds(this->Dim);
1510 std::complex<double> interpolated_u(0.0, 0.0);
1511
1512 // Assemble x and its derivs
1513 for (unsigned l = 0; l < n_node; l++)
1514 {
1515 // Loop over directions
1516 for (unsigned i = 0; i < this->Dim; i++)
1517 {
1518 interpolated_x[i] += this->nodal_position(l, i) * psi[l];
1519 interpolated_dxds[i] += this->nodal_position(l, i) * dpsi(l, 0);
1520 }
1521
1522 // Get the nodal value of the helmholtz unknown
1523 std::complex<double> u_value(
1524 this->nodal_value(l, this->U_index_helmholtz.real()),
1525 this->nodal_value(l, this->U_index_helmholtz.imag()));
1526
1527 interpolated_u += u_value * psi(l);
1528 } // End of loop over the nodes
1529
1530 // calculate the integral
1531 //-----------------------
1532 // define the variable phi_p
1533 double phi_p = atan2(interpolated_x[1], interpolated_x[0]);
1534
1535 // define dphi_ds=(-yx'+y'x)/(x^2+y^2)
1536 double denom = (interpolated_x[0] * interpolated_x[0]) +
1537 (interpolated_x[1] * interpolated_x[1]);
1538 double nom = -interpolated_dxds[1] * interpolated_x[0] +
1539 interpolated_dxds[0] * interpolated_x[1];
1540 double dphi_ds = std::fabs(nom / denom);
1541
1542 // compute the element contribution to gamma
1543 // ALH: The awkward construction with pow and the static_cast is to
1544 // avoid a floating point error on my machine when running unoptimised
1545 // (no idea why!)
1546 gamma_con += (dphi_ds)*w *
1547 pow(exp(I * (phi - phi_p)), static_cast<double>(n)) *
1548 interpolated_u;
1549
1550 // compute the contribution to each node to the map
1551 for (unsigned l = 0; l < n_node; l++)
1552 {
1553 // Add the contribution of the real local data
1554 local_unknown_real =
1555 this->nodal_local_eqn(l, this->U_index_helmholtz.real());
1556 if (local_unknown_real >= 0)
1557 {
1558 global_unknown_real = this->eqn_number(local_unknown_real);
1559 d_gamma_con[global_unknown_real] +=
1560 (dphi_ds)*w * exp(I * (phi - phi_p) * double(n)) * psi(l);
1561 }
1562
1563 // Add the contribution of the imag local data
1564 local_unknown_imag =
1565 this->nodal_local_eqn(l, this->U_index_helmholtz.imag());
1566 if (local_unknown_imag >= 0)
1567 {
1568 global_unknown_imag = this->eqn_number(local_unknown_imag);
1569 // ALH: The awkward construction with pow and the static_cast is to
1570 // avoid a floating point error on my machine when running unoptimised
1571 // (no idea why!)
1572 d_gamma_con[global_unknown_imag] +=
1573 I * (dphi_ds)*w *
1574 pow(exp(I * (phi - phi_p)), static_cast<double>(n)) * psi(l);
1575 }
1576 } // end of loop over the node
1577 } // End of loop over integration points
1578 }
1579
1580
1581 /// //////////////////////////////////////////////////////////////////
1582 /// //////////////////////////////////////////////////////////////////
1583 /// //////////////////////////////////////////////////////////////////
1584
1585
1586 //===========================================================================
1587 /// Namespace for checking radius of nodes on (assumed to be circular)
1588 /// DtN boundary
1589 //===========================================================================
1590 namespace ToleranceForHelmholtzOuterBoundary
1591 {
1592 /// Relative tolerance to within radius of points on DtN boundary
1593 /// are allowed to deviate from specified value
1594 extern double Tol;
1595
1596 } // namespace ToleranceForHelmholtzOuterBoundary
1597
1598
1599 //===========================================================================
1600 /// Namespace for checking radius of nodes on (assumed to be circular)
1601 /// DtN boundary
1602 //===========================================================================
1603 namespace ToleranceForHelmholtzOuterBoundary
1604 {
1605 /// Relative tolerance to within radius of points on DtN boundary
1606 /// are allowed to deviate from specified value
1607 double Tol = 1.0e-3;
1608
1609 } // namespace ToleranceForHelmholtzOuterBoundary
1610
1611 /// //////////////////////////////////////////////////////////////////
1612 /// //////////////////////////////////////////////////////////////////
1613 /// //////////////////////////////////////////////////////////////////
1614
1615
1616 /// ================================================================
1617 /// Compute Fourier components of the solution -- length of
1618 /// vector indicates number of terms to be computed.
1619 /// ================================================================
1620 template<class ELEMENT>
1622 Vector<std::complex<double>>& a_coeff_pos,
1623 Vector<std::complex<double>>& a_coeff_neg)
1624 {
1625#ifdef PARANOID
1626 if (a_coeff_pos.size() != a_coeff_neg.size())
1627 {
1628 std::ostringstream error_stream;
1629 error_stream << "a_coeff_pos and a_coeff_neg must have "
1630 << "the same size. \n";
1631 throw OomphLibError(
1632 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1633 }
1634#endif
1635
1636 // Initialise
1637 unsigned n = a_coeff_pos.size();
1638 Vector<std::complex<double>> el_a_coeff_pos(n);
1639 Vector<std::complex<double>> el_a_coeff_neg(n);
1640 for (unsigned i = 0; i < n; i++)
1641 {
1642 a_coeff_pos[i] = std::complex<double>(0.0, 0.0);
1643 a_coeff_neg[i] = std::complex<double>(0.0, 0.0);
1644 }
1645
1646 // Loop over elements e
1647 unsigned nel = this->nelement();
1648 for (unsigned e = 0; e < nel; e++)
1649 {
1650 // Get a pointer to element
1652 dynamic_cast<HelmholtzBCElementBase<ELEMENT>*>(this->element_pt(e));
1653
1654 // Compute contribution
1655 el_pt->compute_contribution_to_fourier_components(el_a_coeff_pos,
1656 el_a_coeff_neg);
1657
1658 // Add to coefficients
1659 for (unsigned i = 0; i < n; i++)
1660 {
1661 a_coeff_pos[i] += el_a_coeff_pos[i];
1662 a_coeff_neg[i] += el_a_coeff_neg[i];
1663 }
1664 }
1665 }
1666
1667
1668 /// ================================================================
1669 /// Compute and store the gamma integral and derivates
1670 // /w.r.t global unknows at all integration points
1671 /// of the mesh's constituent elements
1672 //================================================================
1673 template<class ELEMENT>
1675 {
1676#ifdef PARANOID
1677 {
1678 // Loop over elements e
1679 unsigned nel = this->nelement();
1680 for (unsigned e = 0; e < nel; e++)
1681 {
1682 FiniteElement* fe_pt = finite_element_pt(e);
1683 unsigned nnod = fe_pt->nnode();
1684 for (unsigned j = 0; j < nnod; j++)
1685 {
1686 Node* nod_pt = fe_pt->node_pt(j);
1687
1688 // Extract nodal coordinates from node:
1689 Vector<double> x(2);
1690 x[0] = nod_pt->x(0);
1691 x[1] = nod_pt->x(1);
1692
1693 // Evaluate the radial distance
1694 double r = sqrt(x[0] * x[0] + x[1] * x[1]);
1695
1696 // Check
1697 if (Outer_radius == 0.0)
1698 {
1699 throw OomphLibError("Outer radius for DtN BC must not be zero!",
1700 OOMPH_CURRENT_FUNCTION,
1701 OOMPH_EXCEPTION_LOCATION);
1702 }
1703
1704 if (std::fabs((r - this->Outer_radius) / Outer_radius) >
1706 {
1707 std::ostringstream error_stream;
1708 error_stream << "Node at " << x[0] << " " << x[1] << " has radius "
1709 << r << " which does not "
1710 << " agree with \nspecified outer radius "
1711 << this->Outer_radius << " within relative tolerance "
1713 << ".\nYou can adjust the tolerance via\n"
1714 << "ToleranceForHelmholtzOuterBoundary::Tol\n"
1715 << "or recompile without PARANOID.\n";
1716 throw OomphLibError(error_stream.str(),
1717 OOMPH_CURRENT_FUNCTION,
1718 OOMPH_EXCEPTION_LOCATION);
1719 }
1720 }
1721 }
1722 }
1723#endif
1724
1725
1726 // Get Helmholtz parameter from bulk element
1728 dynamic_cast<HelmholtzDtNBoundaryElement<ELEMENT>*>(this->element_pt(0));
1729 double k =
1730 sqrt(dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->k_squared());
1731
1732 // Precompute factors in sum
1733 Vector<std::complex<double>> h_a(Nfourier_terms), hp_a(Nfourier_terms),
1734 q(Nfourier_terms);
1736 Nfourier_terms, Outer_radius * k, h_a, hp_a);
1737 for (unsigned i = 0; i < Nfourier_terms; i++)
1738 {
1739 q[i] = (k / (2.0 * MathematicalConstants::Pi)) * (hp_a[i] / h_a[i]);
1740 }
1741
1742 // first loop over elements e
1743 unsigned nel = this->nelement();
1744 for (unsigned e = 0; e < nel; e++)
1745 {
1746 // Get a pointer to element
1749 this->element_pt(e));
1750
1751 // Set the value of n_intpt
1752 const unsigned n_intpt = el_pt->integral_pt()->nweight();
1753
1754 // initialise gamma integral and its derivatives
1755 Vector<std::complex<double>> gamma_vector(n_intpt,
1756 std::complex<double>(0.0, 0.0));
1757 Vector<std::map<unsigned, std::complex<double>>> d_gamma_vector(n_intpt);
1758
1759 // Loop over the integration points
1760 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1761 {
1762 // Allocate and initialise coordinate
1763 unsigned ndim_local = el_pt->dim();
1764 Vector<double> x(ndim_local + 1, 0.0);
1765
1766 // Set the Vector to hold local coordinates
1767 Vector<double> s(ndim_local, 0.0);
1768 for (unsigned i = 0; i < ndim_local; i++)
1769 {
1770 s[i] = el_pt->integral_pt()->knot(ipt, i);
1771 }
1772
1773 // Get the coordinates of the integration point
1774 el_pt->interpolated_x(s, x);
1775
1776 // Polar angle
1777 double phi = atan2(x[1], x[0]);
1778
1779 // Elemental contribution to gamma integral and its derivative
1780 std::complex<double> gamma_con_p(0.0, 0.0), gamma_con_n(0.0, 0.0);
1781 std::map<unsigned, std::complex<double>> d_gamma_con_p, d_gamma_con_n;
1782
1783 // loop over the Fourier terms
1784 for (unsigned nn = 0; nn < Nfourier_terms; nn++)
1785 {
1786 // Second loop over the element
1787 // to evaluate the complete integral
1788 for (unsigned ee = 0; ee < nel; ee++)
1789 {
1792 this->element_pt(ee));
1793
1794 // contribution of the positive term in the sum
1796 phi, int(nn), gamma_con_p, d_gamma_con_p);
1797
1798 // contribution of the negative term in the sum
1800 phi, -int(nn), gamma_con_n, d_gamma_con_n);
1801
1802 unsigned n_node = eel_pt->nnode();
1803 if (nn == 0)
1804 {
1805 gamma_vector[ipt] += q[nn] * gamma_con_p;
1806 for (unsigned l = 0; l < n_node; l++)
1807 {
1808 // Add the contribution of the real local data
1809 int local_unknown_p_real = eel_pt->nodal_local_eqn(
1810 l, eel_pt->u_index_helmholtz().real());
1811 if (local_unknown_p_real >= 0)
1812 {
1813 int global_unknown_p_real =
1814 eel_pt->eqn_number(local_unknown_p_real);
1815 d_gamma_vector[ipt][global_unknown_p_real] +=
1816 q[nn] * d_gamma_con_p[global_unknown_p_real];
1817 }
1818
1819 // Add the contribution of the imag local data
1820 int local_unknown_p_imag = eel_pt->nodal_local_eqn(
1821 l, eel_pt->u_index_helmholtz().imag());
1822
1823 if (local_unknown_p_imag >= 0)
1824 {
1825 int global_unknown_p_imag =
1826 eel_pt->eqn_number(local_unknown_p_imag);
1827
1828 d_gamma_vector[ipt][global_unknown_p_imag] +=
1829 q[nn] * d_gamma_con_p[global_unknown_p_imag];
1830 }
1831 } // end of loop over the node
1832 } // End of if
1833 else
1834 {
1835 gamma_vector[ipt] += q[nn] * (gamma_con_p + gamma_con_n);
1836 for (unsigned l = 0; l < n_node; l++)
1837 {
1838 // Add the contribution of the real local data
1839 int local_unknown_real = eel_pt->nodal_local_eqn(
1840 l, eel_pt->u_index_helmholtz().real());
1841 if (local_unknown_real >= 0)
1842 {
1843 int global_unknown_real =
1844 eel_pt->eqn_number(local_unknown_real);
1845 d_gamma_vector[ipt][global_unknown_real] +=
1846 q[nn] * (d_gamma_con_p[global_unknown_real] +
1847 d_gamma_con_n[global_unknown_real]);
1848 }
1849 // Add the contribution of the imag local data
1850 int local_unknown_imag = eel_pt->nodal_local_eqn(
1851 l, eel_pt->u_index_helmholtz().imag());
1852 if (local_unknown_imag >= 0)
1853 {
1854 int global_unknown_imag =
1855 eel_pt->eqn_number(local_unknown_imag);
1856 d_gamma_vector[ipt][global_unknown_imag] +=
1857 q[nn] * (d_gamma_con_p[global_unknown_imag] +
1858 d_gamma_con_n[global_unknown_imag]);
1859 }
1860 } // end of loop over the node
1861 } // End of else
1862 } // End of second loop over the elements
1863 } // End of loop over Fourier terms
1864 } // end of loop over integration point
1865
1866 // Store it in map
1867 Gamma_at_gauss_point[el_pt] = gamma_vector;
1868 D_Gamma_at_gauss_point[el_pt] = d_gamma_vector;
1869
1870 } // end of first loop over element
1871 }
1872
1873 //===========================================================================
1874 /// Constructor, takes the pointer to the "bulk" element and the face index
1875 //===========================================================================
1876 template<class ELEMENT>
1878 FiniteElement* const& bulk_el_pt, const int& face_index)
1879 : FaceGeometry<ELEMENT>(), FaceElement()
1880 {
1881#ifdef PARANOID
1882 {
1883 // Check that the element is not a refineable 3d element
1884 ELEMENT* elem_pt = new ELEMENT;
1885 // If it's three-d
1886 if (elem_pt->dim() == 3)
1887 {
1888 // Is it refineable
1889 if (dynamic_cast<RefineableElement*>(elem_pt))
1890 {
1891 // Issue a warning
1892 OomphLibWarning("This flux element will not"
1893 "work correctly if nodes are hanging\n",
1894 "HelmholtzBCElementBase::Constructor",
1895 OOMPH_EXCEPTION_LOCATION);
1896 }
1897 }
1898 }
1899#endif
1900
1901 // Let the bulk element build the FaceElement, i.e. setup the pointers
1902 // to its nodes (by referring to the appropriate nodes in the bulk
1903 // element), etc.
1904 bulk_el_pt->build_face_element(face_index, this);
1905
1906 // Extract the dimension of the problem from the dimension of
1907 // the first node
1908 Dim = this->node_pt(0)->ndim();
1909
1910 // Set up U_index_helmholtz. Initialise to zero, which probably won't change
1911 // in most cases, oh well, the price we pay for generality
1912 U_index_helmholtz = std::complex<unsigned>(0, 1);
1913
1914 // Cast to the appropriate HelmholtzEquation so that we can
1915 // find the index at which the variable is stored
1916 // We assume that the dimension of the full problem is the same
1917 // as the dimension of the node, if this is not the case you will have
1918 // to write custom elements, sorry
1919 switch (Dim)
1920 {
1921 // One dimensional problem
1922 case 1:
1923 {
1924 HelmholtzEquations<1>* eqn_pt =
1925 dynamic_cast<HelmholtzEquations<1>*>(bulk_el_pt);
1926 // If the cast has failed die
1927 if (eqn_pt == 0)
1928 {
1929 std::string error_string =
1930 "Bulk element must inherit from HelmholtzEquations.";
1931 error_string +=
1932 "Nodes are one dimensional, but cannot cast the bulk element to\n";
1933 error_string += "HelmholtzEquations<1>\n.";
1934 error_string += "If you desire this functionality, you must "
1935 "implement it yourself\n";
1936
1937 throw OomphLibError(
1938 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1939 }
1940 // Otherwise read out the value
1941 else
1942 {
1943 // Read the index from the (cast) bulk element
1945 }
1946 }
1947 break;
1948
1949 // Two dimensional problem
1950 case 2:
1951 {
1952 HelmholtzEquations<2>* eqn_pt =
1953 dynamic_cast<HelmholtzEquations<2>*>(bulk_el_pt);
1954 // If the cast has failed die
1955 if (eqn_pt == 0)
1956 {
1957 std::string error_string =
1958 "Bulk element must inherit from HelmholtzEquations.";
1959 error_string +=
1960 "Nodes are two dimensional, but cannot cast the bulk element to\n";
1961 error_string += "HelmholtzEquations<2>\n.";
1962 error_string += "If you desire this functionality, you must "
1963 "implement it yourself\n";
1964
1965 throw OomphLibError(
1966 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1967 }
1968 else
1969 {
1970 // Read the index from the (cast) bulk element
1972 }
1973 }
1974
1975 break;
1976
1977 // Three dimensional problem
1978 case 3:
1979 {
1980 HelmholtzEquations<3>* eqn_pt =
1981 dynamic_cast<HelmholtzEquations<3>*>(bulk_el_pt);
1982 // If the cast has failed die
1983 if (eqn_pt == 0)
1984 {
1985 std::string error_string =
1986 "Bulk element must inherit from HelmholtzEquations.";
1987 error_string += "Nodes are three dimensional, but cannot cast the "
1988 "bulk element to\n";
1989 error_string += "HelmholtzEquations<3>\n.";
1990 error_string += "If you desire this functionality, you must "
1991 "implement it yourself\n";
1992
1993 throw OomphLibError(
1994 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1995 }
1996 else
1997 {
1998 // Read the index from the (cast) bulk element
2000 }
2001 }
2002 break;
2003
2004 // Any other case is an error
2005 default:
2006 std::ostringstream error_stream;
2007 error_stream << "Dimension of node is " << Dim
2008 << ". It should be 1,2, or 3!" << std::endl;
2009
2010 throw OomphLibError(
2011 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2012 break;
2013 }
2014 }
2015
2016
2017} // namespace oomph
2018
2019#endif
e
Definition: cfortran.h:571
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
virtual bool is_a_copy() const
Return a boolean to indicate whether the Data objact contains any copied values. A base Data object c...
Definition: nodes.h:253
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary,...
Definition: elements.h:4497
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
virtual 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...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
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
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1981
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 nexternal_data() const
Return the number of external data objects.
Definition: elements.h:829
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
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
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data.
Definition: elements.h:311
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
Definition: elements.cc:307
////////////////////////////////////////////////////////////////// //////////////////////////////////...
void fill_in_generic_residual_contribution_helmholtz_abc(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute the element's residual vector and the ( Jacobian matrix. Overloaded version,...
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.
unsigned * ABC_order_pt
Pointer to order of absorbing boundary condition.
double * Outer_radius_pt
Pointer to radius of outer boundary (must be a circle!)
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
HelmholtzAbsorbingBCElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Construct element from specification of bulk element and face index.
double *& outer_radius_pt()
Get pointer to radius of outer boundary (must be a cirle)
unsigned *& abc_order_pt()
Pointer to order of absorbing boundary condition.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
double test_only(const Vector< double > &s, 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 – forward to broken version in FiniteElement until somebody decides what exactly they...
double global_power_contribution()
Compute the element's contribution to the time-averaged radiated power over the artificial boundary.
double global_power_contribution(std::ofstream &outfile)
Compute the element's contribution to the time-averaged radiated power over the artificial boundary....
virtual std::complex< unsigned > u_index_helmholtz() const
Return the index at which the real/imag unknown value is stored.
HelmholtzBCElementBase()
Broken empty constructor.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Broken assignment operator.
void output(FILE *file_pt)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
HelmholtzBCElementBase(const HelmholtzBCElementBase &dummy)=delete
Broken copy constructor.
void compute_contribution_to_fourier_components(Vector< std::complex< double > > &a_coeff_pos, Vector< std::complex< double > > &a_coeff_neg)
Compute element's contribution to Fourier components of the solution – length of vector indicates num...
double d_shape_and_test_local(const Vector< double > &s, Shape &psi, Shape &test, DShape &dpsi_ds, DShape &dtest_ds) const
Function to compute the shape, test functions and derivates and to return the Jacobian of mapping bet...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
unsigned Dim
The spatial dimension of the problem.
std::complex< unsigned > U_index_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
FaceElement used to apply Sommerfeld radiation conditon via Dirichlet to Neumann map.
void fill_in_generic_residual_contribution_helmholtz_DtN_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute the element's residual vector Jacobian matrix. Overloaded version, using the gamma computed i...
HelmholtzDtNBoundaryElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Construct element from specification of bulk element and face index.
void compute_gamma_contribution(const double &phi, const int &n, std::complex< double > &gamma_con, std::map< unsigned, std::complex< double > > &d_gamma_con)
Compute the contribution of the element to the Gamma integral and its derivates w....
void set_outer_boundary_mesh_pt(HelmholtzDtNMesh< ELEMENT > *mesh_pt)
Set mesh of all DtN boundary condition elements.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
HelmholtzDtNMesh< ELEMENT > * outer_boundary_mesh_pt() const
Access function to mesh of all DtN boundary condition elements (needed to get access to gamma values)
HelmholtzDtNMesh< ELEMENT > * Outer_boundary_mesh_pt
Pointer to mesh of all DtN boundary condition elements (needed to get access to gamma values)
void complete_setup_of_dependencies()
Complete the setup of additional dependencies arising through the far-away interaction with other nod...
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.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Vector< std::map< unsigned, std::complex< double > > > & d_gamma_at_gauss_point(FiniteElement *el_pt)
Derivative of Gamma integral w.r.t global unknown, evaluated at Gauss points for specified element.
std::map< FiniteElement *, Vector< std::complex< double > > > Gamma_at_gauss_point
Container to store the gamma integral for given Gauss point and element.
std::map< FiniteElement *, Vector< std::map< unsigned, std::complex< double > > > > D_Gamma_at_gauss_point
Container to store the derivate of Gamma integral w.r.t global unknown evaluated at Gauss points for ...
unsigned Nfourier_terms
Nbr of Fourier terms used in the Gamma computation.
void setup_gamma()
Compute and store the gamma integral at all integration points of the constituent elements.
double & outer_radius()
The outer radius.
void compute_fourier_components(Vector< std::complex< double > > &a_coeff_pos, Vector< std::complex< double > > &a_coeff_neg)
Compute Fourier components of the solution – length of vector indicates number of terms to be compute...
double Outer_radius
Outer radius.
Vector< std::complex< double > > & gamma_at_gauss_point(FiniteElement *el_pt)
Gamma integral evaluated at Gauss points for specified element.
HelmholtzDtNMesh(const double &outer_radius, const unsigned &nfourier_terms)
Constructor: Specify radius of outer boundary and number of Fourier terms used in the computation of ...
unsigned & nfourier_terms()
Number of Fourier terms used in the computation of the gamma integral.
A class for all isoparametric elements that solve the Helmholtz equations.
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
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.
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
virtual Node * copied_node_pt() const
Return pointer to copied node (null if the current node is not a copy – always the case here; it's ov...
Definition: nodes.h:1107
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
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
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void CHankel_first(const unsigned &n, const std::complex< double > &x, Vector< std::complex< double > > &h, Vector< std::complex< double > > &hp)
Compute Hankel function of the first kind of orders 0...n and its derivates at coordinate x....
void Hankel_first(const unsigned &n, const double &x, Vector< std::complex< double > > &h, Vector< std::complex< double > > &hp)
Compute Hankel function of the first kind of orders 0...n and its derivates at coordinate x....
const std::complex< double > I(0.0, 1.0)
The imaginary unit.
const double Pi
50 digits from maple
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
double Tol
Relative tolerance to within radius of points on DtN boundary are allowed to deviate from specified v...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...