fourier_decomposed_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 Fourier decomposed Helmholtz equations
28#ifndef OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_BC_ELEMENTS_HEADER
29#define OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_BC_ELEMENTS_HEADER
30
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37#include "math.h"
38#include <complex>
39
40// Get the Bessel functions
41#include "oomph_crbond_bessel.h"
42
44
45namespace oomph
46{
47 /// //////////////////////////////////////////////////////////////////
48 /// //////////////////////////////////////////////////////////////////
49 /// //////////////////////////////////////////////////////////////////
50
51
52 //======================================================================
53 /// A class for elements that allow the approximation of the
54 /// Sommerfeld radiation BC for Fourier decomposed Helmholtz equations.
55 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
56 /// policy class.
57 //======================================================================
58 template<class ELEMENT>
60 : public virtual FaceGeometry<ELEMENT>,
61 public virtual FaceElement
62 {
63 public:
64 /// Constructor, takes the pointer to the "bulk" element and the
65 /// index of the face to which the element is attached.
67 const int& face_index);
68
69 /// Broken empty constructor
71 {
72 throw OomphLibError("Don't call empty constructor for "
73 "FourierDecomposedHelmholtzBCElementBase",
74 OOMPH_CURRENT_FUNCTION,
75 OOMPH_EXCEPTION_LOCATION);
76 }
77
78 /// Broken copy constructor
80 const FourierDecomposedHelmholtzBCElementBase& dummy) = delete;
81
82 /// Broken assignment operator
83 // Commented out broken assignment operator because this can lead to a
84 // conflict warning when used in the virtual inheritence hierarchy.
85 // Essentially the compiler doesn't realise that two separate
86 // implementations of the broken function are the same and so, quite
87 // rightly, it shouts.
88 /*void operator=(const FourierDecomposedHelmholtzBCElementBase&) = delete;*/
89
90
91 /// Specify the value of nodal zeta from the face geometry
92 /// The "global" intrinsic coordinate of the element when
93 /// viewed as part of a geometric object should be given by
94 /// the FaceElement representation, by default (needed to break
95 /// indeterminacy if bulk element is SolidElement)
96 double zeta_nodal(const unsigned& n,
97 const unsigned& k,
98 const unsigned& i) const
99 {
100 return FaceElement::zeta_nodal(n, k, i);
101 }
102
103
104 /// Output function -- forward to broken version in FiniteElement
105 /// until somebody decides what exactly they want to plot here...
106 void output(std::ostream& outfile)
107 {
108 FiniteElement::output(outfile);
109 }
110
111 /// Output function -- forward to broken version in FiniteElement
112 /// until somebody decides what exactly they want to plot here...
113 void output(std::ostream& outfile, const unsigned& n_plot)
114 {
115 FiniteElement::output(outfile, n_plot);
116 }
117
118 /// C-style output function -- forward to broken version in FiniteElement
119 /// until somebody decides what exactly they want to plot here...
120 void output(FILE* file_pt)
121 {
122 FiniteElement::output(file_pt);
123 }
124
125 /// C-style output function -- forward to broken version in
126 /// FiniteElement until somebody decides what exactly they want to plot
127 /// here...
128 void output(FILE* file_pt, const unsigned& n_plot)
129 {
130 FiniteElement::output(file_pt, n_plot);
131 }
132
133 /// Return the index at which the real/imag unknown value
134 /// is stored.
135 virtual inline std::complex<unsigned> u_index_fourier_decomposed_helmholtz()
136 const
137 {
138 return std::complex<unsigned>(
141 }
142
143 /// Compute the element's contribution to the time-averaged
144 /// radiated power over the artificial boundary
146 {
147 // Dummy output file
148 std::ofstream outfile;
149 return global_power_contribution(outfile);
150 }
151
152 /// Compute the element's contribution to the time-averaged
153 /// radiated power over the artificial boundary. Also output the
154 /// power density as a fct of the zenith angle in the specified
155 /// output file if it's open.
156 double global_power_contribution(std::ofstream& outfile)
157 {
158 // pointer to the corresponding bulk element
159 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
160
161 // Number of nodes in bulk element
162 unsigned nnode_bulk = bulk_elem_pt->nnode();
163 const unsigned n_node_local = this->nnode();
164
165 // get the dim of the bulk and local nodes
166 const unsigned bulk_dim = bulk_elem_pt->dim();
167 const unsigned local_dim = this->node_pt(0)->ndim();
168
169 // Set up memory for the shape and test functions
170 Shape psi(n_node_local);
171
172 // Set up memory for the shape functions
173 Shape psi_bulk(nnode_bulk);
174 DShape dpsi_bulk_dx(nnode_bulk, bulk_dim);
175
176 // Set up memory for the outer unit normal
177 Vector<double> unit_normal(bulk_dim);
178
179 // Set the value of n_intpt
180 const unsigned n_intpt = integral_pt()->nweight();
181
182 // Set the Vector to hold local coordinates
183 Vector<double> s(local_dim - 1);
184 double power = 0.0;
185
186 // Output?
187 if (outfile.is_open())
188 {
189 outfile << "ZONE\n";
190 }
191
192 // Loop over the integration points
193 //--------------------------------
194 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
195 {
196 // Assign values of s
197 for (unsigned i = 0; i < (local_dim - 1); i++)
198 {
199 s[i] = integral_pt()->knot(ipt, i);
200 }
201 // get the outer_unit_ext vector
202 this->outer_unit_normal(s, unit_normal);
203
204 // Get the integral weight
205 double w = integral_pt()->weight(ipt);
206
207 // Get jacobian of mapping
208 double J = J_eulerian(s);
209
210 // Premultiply the weights and the Jacobian
211 double W = w * J;
212
213 // Get local coordinates in bulk element by copy construction
215
216 // Call the derivatives of the shape functions
217 // in the bulk -- must do this via s because this point
218 // is not an integration point the bulk element!
219 (void)bulk_elem_pt->dshape_eulerian(s_bulk, psi_bulk, dpsi_bulk_dx);
220 this->shape(s, psi);
221
222 // Derivs of Eulerian coordinates w.r.t. local coordinates
223 std::complex<double> dphi_dn(0.0, 0.0);
224 Vector<std::complex<double>> interpolated_dphidx(bulk_dim);
225 std::complex<double> interpolated_phi(0.0, 0.0);
226 Vector<double> x(bulk_dim);
227
228 // Calculate function value and derivatives:
229 //-----------------------------------------
230 // Loop over nodes
231 for (unsigned l = 0; l < nnode_bulk; l++)
232 {
233 // Get the nodal value of the helmholtz unknown
234 const std::complex<double> phi_value(
235 bulk_elem_pt->nodal_value(
236 l, bulk_elem_pt->u_index_fourier_decomposed_helmholtz().real()),
237 bulk_elem_pt->nodal_value(
238 l, bulk_elem_pt->u_index_fourier_decomposed_helmholtz().imag()));
239
240 // Loop over directions
241 for (unsigned i = 0; i < bulk_dim; i++)
242 {
243 interpolated_dphidx[i] += phi_value * dpsi_bulk_dx(l, i);
244 }
245 } // End of loop over the bulk_nodes
246
247 for (unsigned l = 0; l < n_node_local; l++)
248 {
249 // Get the nodal value of the Helmholtz unknown
250 const std::complex<double> phi_value(
253
254 interpolated_phi += phi_value * psi(l);
255 }
256
257 // define dphi_dn
258 for (unsigned i = 0; i < bulk_dim; i++)
259 {
260 dphi_dn += interpolated_dphidx[i] * unit_normal[i];
261 }
262
263 // Power density
264 double integrand = (interpolated_phi.real() * dphi_dn.imag() -
265 interpolated_phi.imag() * dphi_dn.real());
266
267 interpolated_x(s, x);
268 double theta = atan2(x[0], x[1]);
269 // Output?
270 if (outfile.is_open())
271 {
272 outfile << x[0] << " " << x[1] << " " << theta << " " << integrand
273 << "\n";
274 }
275
276 // ...add to integral
277 power += MathematicalConstants::Pi * x[0] * integrand * W;
278 }
279
280 return power;
281 }
282
283 protected:
284 /// Function to compute the test functions and to return
285 /// the Jacobian of mapping between local and global (Eulerian)
286 /// coordinates
287 inline double shape_and_test(const Vector<double>& s,
288 Shape& psi,
289 Shape& test) const
290 {
291 // Get the shape functions
292 shape(s, test);
293
294 unsigned nnod = nnode();
295 for (unsigned i = 0; i < nnod; i++)
296 {
297 psi[i] = test[i];
298 }
299
300 // Return the value of the jacobian
301 return J_eulerian(s);
302 }
303
304 /// Function to compute the shape, test functions and derivates
305 /// and to return
306 /// the Jacobian of mapping between local and global (Eulerian)
307 /// coordinates
309 Shape& psi,
310 Shape& test,
311 DShape& dpsi_ds,
312 DShape& dtest_ds) const
313 {
314 // Find number of nodes
315 unsigned n_node = nnode();
316
317 // Get the shape functions
318 dshape_local(s, psi, dpsi_ds);
319
320 // Set the test functions to be the same as the shape functions
321 for (unsigned i = 0; i < n_node; i++)
322 {
323 for (unsigned j = 0; j < 1; j++)
324 {
325 test[i] = psi[i];
326 dtest_ds(i, j) = dpsi_ds(i, j);
327 }
328 }
329 // Return the value of the jacobian
330 return J_eulerian(s);
331 }
332
333 /// The index at which the real and imag part of the unknown is
334 /// stored at the nodes
336 };
337
338
339 /// ///////////////////////////////////////////////////////////////////
340 /// ///////////////////////////////////////////////////////////////////
341 /// ///////////////////////////////////////////////////////////////////
342
343
344 /// =================================================================
345 /// Mesh for DtN boundary condition elements -- provides
346 /// functionality to apply Sommerfeld radiation condtion
347 /// via DtN BC.
348 /// =================================================================
349 template<class ELEMENT>
351 {
352 public:
353 /// Constructor: Specify radius of outer boundary and number of
354 /// terms used in the computation of the gamma integral
356 const unsigned& n_terms)
358 {
359 }
360
361 /// Compute and store the gamma integral at all integration
362 /// points of the constituent elements.
363 void setup_gamma();
364
365 /// Gamma integral evaluated at Gauss points
366 /// for specified element
368 {
369 return Gamma_at_gauss_point[el_pt];
370 }
371
372 /// Derivative of Gamma integral w.r.t global unknown, evaluated
373 /// at Gauss points for specified element
375 FiniteElement* el_pt)
376 {
377 return D_Gamma_at_gauss_point[el_pt];
378 }
379
380 /// The outer radius
381 double& outer_radius()
382 {
383 return Outer_radius;
384 }
385
386 /// Number of terms used in the computation of the
387 /// gamma integral
388 unsigned& n_terms()
389 {
390 return N_terms;
391 }
392
393 private:
394 /// Outer radius
396
397 /// Number of terms used in the Gamma computation
398 unsigned N_terms;
399
400
401 /// Container to store the gamma integral for given Gauss point
402 /// and element
403 std::map<FiniteElement*, Vector<std::complex<double>>> Gamma_at_gauss_point;
404
405
406 /// Container to store the derivate of Gamma integral w.r.t
407 /// global unknown evaluated at Gauss points for specified element
408 std::map<FiniteElement*, Vector<std::map<unsigned, std::complex<double>>>>
410 };
411
412 /// //////////////////////////////////////////////////////////////////
413 /// //////////////////////////////////////////////////////////////////
414
415 //=============================================================
416 /// FaceElement used to apply Sommerfeld radiation conditon
417 /// via Dirichlet to Neumann map.
418 //==============================================================
419 template<class ELEMENT>
422 {
423 public:
424 /// Construct element from specification of bulk element and
425 /// face index.
427 FiniteElement* const& bulk_el_pt, const int& face_index)
429 {
430 }
431
432 /// Add the element's contribution to its residual vector
434 {
435 // Call the generic residuals function with flag set to 0
436 // using a dummy matrix argument
439 }
440
441
442 /// Add the element's contribution to its residual vector and its
443 /// Jacobian matrix
445 DenseMatrix<double>& jacobian)
446 {
447 // Call the generic routine with the flag set to 1
449 residuals, jacobian, 1);
450 }
451
452 /// Compute the contribution of the element
453 /// to the Gamma integral and its derivates w.r.t
454 /// to global unknows; the function takes the wavenumber (for gamma
455 /// integral, not the one from the Fourier decomposition of the Helmholtz
456 /// equations!) and the polar angle theta as input.
458 const double& theta,
459 const unsigned& n,
460 std::complex<double>& gamma_con,
461 std::map<unsigned, std::complex<double>>& d_gamma_con);
462
463
464 /// Access function to mesh of all DtN boundary condition elements
465 /// (needed to get access to gamma values)
467 {
469 }
470
471 /// Set mesh of all DtN boundary condition elements
474 {
475 Outer_boundary_mesh_pt = mesh_pt;
476 }
477
478
479 /// Complete the setup of additional dependencies arising
480 /// through the far-away interaction with other nodes in
481 /// Outer_boundary_mesh_pt.
483 {
484 // Create a set of all nodes
485 std::set<Node*> node_set;
486 unsigned nel = Outer_boundary_mesh_pt->nelement();
487 for (unsigned e = 0; e < nel; e++)
488 {
489 FiniteElement* el_pt = Outer_boundary_mesh_pt->finite_element_pt(e);
490 unsigned nnod = el_pt->nnode();
491 for (unsigned j = 0; j < nnod; j++)
492 {
493 Node* nod_pt = el_pt->node_pt(j);
494
495 // Don't add copied nodes
496 if (!(nod_pt->is_a_copy()))
497 {
498 node_set.insert(nod_pt);
499 }
500 }
501 }
502 // Now erase the current element's own nodes
503 unsigned nnod = this->nnode();
504 for (unsigned j = 0; j < nnod; j++)
505 {
506 Node* nod_pt = this->node_pt(j);
507 node_set.erase(nod_pt);
508
509 // If the element's node is a copy then its "master" will
510 // already have been added in the set above -- remove the
511 // master to avoid double counting eqn numbers
512 if (nod_pt->is_a_copy())
513 {
514 node_set.erase(nod_pt->copied_node_pt());
515 }
516 }
517
518 // Now declare these nodes to be the element's external Data
519 for (std::set<Node*>::iterator it = node_set.begin();
520 it != node_set.end();
521 it++)
522 {
523 this->add_external_data(*it);
524 }
525 }
526
527
528 private:
529 /// Compute the element's residual vector
530 /// Jacobian matrix.
531 /// Overloaded version, using the gamma computed in the mesh
533 Vector<double>& residuals,
534 DenseMatrix<double>& jacobian,
535 const unsigned& flag)
536 {
537 // Find out how many nodes there are
538 const unsigned n_node = this->nnode();
539
540 // Set up memory for the shape and test functions
541 Shape test(n_node);
542 Shape psi(n_node);
543
544 // Set the value of Nintpt
545 const unsigned n_intpt = this->integral_pt()->nweight();
546
547 // Set the Vector to hold local coordinates
548 Vector<double> s(1);
549
550 // Integers to hold the local equation and unknown numbers
551 int local_eqn_real = 0, local_unknown_real = 0, global_unk_real = 0,
552 local_eqn_imag = 0, local_unknown_imag = 0, global_unk_imag = 0;
553 int external_global_unk_real = 0, external_unknown_real = 0,
554 external_global_unk_imag = 0, external_unknown_imag = 0;
555
556
557 // Get the gamma value for the current integration point
558 // from the mesh
560 Outer_boundary_mesh_pt->gamma_at_gauss_point(this));
561
563 Outer_boundary_mesh_pt->d_gamma_at_gauss_point(this));
564
565 // Loop over the integration points
566 //--------------------------------
567 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
568 {
569 // Assign values of s
570 for (unsigned i = 0; i < 1; i++)
571 {
572 s[i] = this->integral_pt()->knot(ipt, i);
573 }
574
575 // Get the integral weight
576 double w = this->integral_pt()->weight(ipt);
577
578 // Find the shape test functions and derivates; return the Jacobian
579 // of the mapping between local and global (Eulerian)
580 // coordinates
581 double J = this->shape_and_test(s, psi, test);
582
583 // Premultiply the weights and the Jacobian
584 double W = w * J;
585
586 // Build up radius
587 double r = 0.0;
588 for (unsigned j = 0; j < n_node; j++)
589 {
590 r += this->node_pt(j)->x(0) * psi(j);
591 }
592
593 // Now add to the appropriate equations
594 // Loop over the test functions:loop over the nodes
595 for (unsigned l = 0; l < n_node; l++)
596 {
597 local_eqn_real = this->nodal_local_eqn(
599
600 local_eqn_imag = this->nodal_local_eqn(
602
603 // IF it's not a boundary condition
604 if (local_eqn_real >= 0)
605 {
606 // Add the gamma contribution in this int_point to the res
607 residuals[local_eqn_real] += gamma[ipt].real() * test[l] * r * W;
608
609 // Calculate the jacobian
610 //-----------------------
611 if (flag)
612 {
613 // Loop over the shape functions again
614 for (unsigned l2 = 0; l2 < n_node; l2++)
615 {
616 // Add the contribution of the local data
617 local_unknown_real = this->nodal_local_eqn(
618 l2, this->U_index_fourier_decomposed_helmholtz.real());
619
620 local_unknown_imag = this->nodal_local_eqn(
621 l2, this->U_index_fourier_decomposed_helmholtz.imag());
622
623 // If at a non-zero degree of freedom add in the entry
624 if (local_unknown_real >= 0)
625 {
626 // Add the contribution
627 global_unk_real = this->eqn_number(local_unknown_real);
628 jacobian(local_eqn_real, local_unknown_real) +=
629 d_gamma[ipt][global_unk_real].real() * test[l] * r * W;
630 }
631 if (local_unknown_imag >= 0)
632 {
633 // Add the contribution
634 global_unk_imag = this->eqn_number(local_unknown_imag);
635 jacobian(local_eqn_real, local_unknown_imag) +=
636 d_gamma[ipt][global_unk_imag].real() * test[l] * r * W;
637 }
638 } // End of loop over nodes l2
639
640 // Add the contribution of the external data
641 unsigned n_ext_data = this->nexternal_data();
642 // Loop over the shape functions again
643 for (unsigned l2 = 0; l2 < n_ext_data; l2++)
644 {
645 // Add the contribution of the local data
646 external_unknown_real = this->external_local_eqn(
647 l2, this->U_index_fourier_decomposed_helmholtz.real());
648
649 external_unknown_imag = this->external_local_eqn(
650 l2, this->U_index_fourier_decomposed_helmholtz.imag());
651
652 // If at a non-zero degree of freedom add in the entry
653 if (external_unknown_real >= 0)
654 {
655 // Add the contribution
656 external_global_unk_real =
657 this->eqn_number(external_unknown_real);
658 jacobian(local_eqn_real, external_unknown_real) +=
659 d_gamma[ipt][external_global_unk_real].real() * test[l] *
660 r * W;
661 }
662 if (external_unknown_imag >= 0)
663 {
664 // Add the contribution
665 external_global_unk_imag =
666 this->eqn_number(external_unknown_imag);
667 jacobian(local_eqn_real, external_unknown_imag) +=
668 d_gamma[ipt][external_global_unk_imag].real() * test[l] *
669 r * W;
670 }
671 } // End of loop over external data
672 } // End of flag
673 } // end of local_eqn_real
674
675 if (local_eqn_imag >= 0)
676 {
677 // Add the gamma contribution in this int_point to the res
678 residuals[local_eqn_imag] += gamma[ipt].imag() * test[l] * r * W;
679
680 // Calculate the jacobian
681 //-----------------------
682 if (flag)
683 {
684 // Loop over the shape functions again
685 for (unsigned l2 = 0; l2 < n_node; l2++)
686 {
687 // Add the contribution of the local data
688 local_unknown_real = this->nodal_local_eqn(
689 l2, this->U_index_fourier_decomposed_helmholtz.real());
690
691 local_unknown_imag = this->nodal_local_eqn(
692 l2, this->U_index_fourier_decomposed_helmholtz.imag());
693
694 // If at a non-zero degree of freedom add in the entry
695 if (local_unknown_real >= 0)
696 {
697 // Add the contribution
698 global_unk_real = this->eqn_number(local_unknown_real);
699 jacobian(local_eqn_imag, local_unknown_real) +=
700 d_gamma[ipt][global_unk_real].imag() * test[l] * r * W;
701 }
702 if (local_unknown_imag >= 0)
703 {
704 // Add the contribution
705 global_unk_imag = this->eqn_number(local_unknown_imag);
706 jacobian(local_eqn_imag, local_unknown_imag) +=
707 d_gamma[ipt][global_unk_imag].imag() * test[l] * r * W;
708 }
709 } // End of loop over nodes l2
710
711 // Add the contribution of the external data
712 unsigned n_ext_data = this->nexternal_data();
713
714 // Loop over the shape functions again
715 for (unsigned l2 = 0; l2 < n_ext_data; l2++)
716 {
717 // Add the contribution of the local data
718 external_unknown_real = this->external_local_eqn(
719 l2, this->U_index_fourier_decomposed_helmholtz.real());
720
721 external_unknown_imag = this->external_local_eqn(
722 l2, this->U_index_fourier_decomposed_helmholtz.imag());
723
724 // If at a non-zero degree of freedom add in the entry
725 if (external_unknown_real >= 0)
726 {
727 // Add the contribution
728 external_global_unk_real =
729 this->eqn_number(external_unknown_real);
730 jacobian(local_eqn_imag, external_unknown_real) +=
731 d_gamma[ipt][external_global_unk_real].imag() * test[l] *
732 r * W;
733 }
734 if (external_unknown_imag >= 0)
735 {
736 // Add the contribution
737 external_global_unk_imag =
738 this->eqn_number(external_unknown_imag);
739 jacobian(local_eqn_imag, external_unknown_imag) +=
740 d_gamma[ipt][external_global_unk_imag].imag() * test[l] *
741 r * W;
742 }
743 } // End of loop over external data
744 } // End of flag
745 } // end of local_eqn_imag
746 } // end of loop over the nodes
747 } // End of loop over int_pt
748 }
749
750
751 /// Pointer to mesh of all DtN boundary condition elements
752 /// (needed to get access to gamma values)
754 };
755
756
757 /// /////////////////////////////////////////////////////////////
758 /// /////////////////////////////////////////////////////////////
759 /// /////////////////////////////////////////////////////////////
760
761
762 //===========start_compute_gamma_contribution==================
763 /// compute the contribution of the element
764 /// to the Gamma integral and its derivates w.r.t
765 /// to global unknows; the function takes wavenumber n
766 /// (from the computation of the gamma integral, not the one
767 /// from the Fourier decomposition of the Helmholtz equations!)
768 /// and polar angle theta as input.
769 //==============================================================
770 template<class ELEMENT>
773 const double& theta,
774 const unsigned& n,
775 std::complex<double>& gamma_con,
776 std::map<unsigned, std::complex<double>>& d_gamma_con)
777 {
778 // Parameters
779 int n_fourier_helmholtz =
780 dynamic_cast<ELEMENT*>(this->bulk_element_pt())->fourier_wavenumber();
781
782 // define the imaginary number
783 const std::complex<double> I(0.0, 1.0);
784
785 // Find out how many nodes there are
786 const unsigned n_node = this->nnode();
787
788 // Set up memory for the shape functions
789 Shape psi(n_node);
790 DShape dpsi(n_node, 1);
791
792 // initialise the variable
793 int local_unknown_real = 0, local_unknown_imag = 0;
794 int global_eqn_real = 0, global_eqn_imag = 0;
795
796 // Set the value of n_intpt
797 const unsigned n_intpt = this->integral_pt()->nweight();
798
799 // Set the Vector to hold local coordinates
800 Vector<double> s(1);
801
802 // Initialise
803 gamma_con = std::complex<double>(0.0, 0.0);
804 d_gamma_con.clear();
805
806 // Loop over the integration points
807 //--------------------------------
808 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
809 {
810 // Assign values of s
811 for (unsigned i = 0; i < 1; i++)
812 {
813 s[i] = this->integral_pt()->knot(ipt, i);
814 }
815
816 // Get the integral weight
817 double w = this->integral_pt()->weight(ipt);
818
819 // Get the shape functions
820 this->dshape_local(s, psi, dpsi);
821
822 // Eulerian coordinates at Gauss point
823 Vector<double> interpolated_x(2, 0.0);
824
825 // Derivs of Eulerian coordinates w.r.t. local coordinates
826 Vector<double> interpolated_dxds(2);
827 std::complex<double> interpolated_u(0.0, 0.0);
828
829 // Assemble x and its derivs
830 for (unsigned l = 0; l < n_node; l++)
831 {
832 // Loop over directions
833 for (unsigned i = 0; i < 2; i++)
834 {
835 interpolated_x[i] += this->nodal_position(l, i) * psi[l];
836 interpolated_dxds[i] += this->nodal_position(l, i) * dpsi(l, 0);
837 }
838
839 // Get the nodal value of the helmholtz unknown
840 std::complex<double> u_value(
841 this->nodal_value(l,
842 this->U_index_fourier_decomposed_helmholtz.real()),
843 this->nodal_value(l,
844 this->U_index_fourier_decomposed_helmholtz.imag()));
845
846 interpolated_u += u_value * psi(l);
847
848 } // End of loop over the nodes
849
850
851 // calculate the integral
852 //-----------------------
853 // define the variable theta
854 double phi = atan2(interpolated_x[0], interpolated_x[1]);
855
856 // define dphi_ds=-z'/r
857 double dphi_ds = -std::fabs(-interpolated_dxds[1] / interpolated_x[0]);
858
859 // define the associated legendre functions
860 double p_theta =
861 Legendre_functions_helper::plgndr2(n, n_fourier_helmholtz, cos(theta));
862
863 double p_phi =
864 Legendre_functions_helper::plgndr2(n, n_fourier_helmholtz, cos(phi));
865
866 gamma_con += interpolated_u * p_phi * p_theta * sin(phi) * w * dphi_ds;
867
868 // compute the contribution to each node to the map
869 for (unsigned l = 0; l < n_node; l++)
870 {
871 // Add the contribution of the real local data
872 local_unknown_real = this->nodal_local_eqn(
873 l, this->U_index_fourier_decomposed_helmholtz.real());
874 if (local_unknown_real >= 0)
875 {
876 global_eqn_real = this->eqn_number(local_unknown_real);
877 d_gamma_con[global_eqn_real] +=
878 psi(l) * p_phi * p_theta * sin(phi) * w * dphi_ds;
879 }
880
881 // Add the contribution of the imag local data
882 local_unknown_imag = this->nodal_local_eqn(
883 l, this->U_index_fourier_decomposed_helmholtz.imag());
884 if (local_unknown_imag >= 0)
885 {
886 global_eqn_imag = this->eqn_number(local_unknown_imag);
887 d_gamma_con[global_eqn_imag] +=
888 I * psi(l) * p_phi * p_theta * sin(phi) * w * dphi_ds;
889 }
890 } // end of loop over the node
891 } // End of loop over integration points
892 }
893
894
895 /// //////////////////////////////////////////////////////////////////
896 /// //////////////////////////////////////////////////////////////////
897 /// //////////////////////////////////////////////////////////////////
898
899
900 //===========================================================================
901 /// Namespace for checking radius of nodes on (assumed to be circular)
902 /// DtN boundary
903 //===========================================================================
904 namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
905 {
906 /// Relative tolerance to within radius of points on DtN boundary
907 /// are allowed to deviate from specified value
908 extern double Tol;
909
910 } // namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
911
912
913 //===========================================================================
914 /// Namespace for checking radius of nodes on (assumed to be circular)
915 /// DtN boundary
916 //===========================================================================
917 namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
918 {
919 /// Relative tolerance to within radius of points on DtN boundary
920 /// are allowed to deviate from specified value
921 double Tol = 1.0e-3;
922
923 } // namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
924
925 /// //////////////////////////////////////////////////////////////////
926 /// //////////////////////////////////////////////////////////////////
927 /// //////////////////////////////////////////////////////////////////
928
929
930 /// ================================================================
931 /// Compute and store the gamma integral and derivates
932 // /w.r.t global unknows at all integration points
933 /// of the mesh's constituent elements
934 //================================================================
935 template<class ELEMENT>
937 {
938#ifdef PARANOID
939 {
940 // Loop over elements e
941 unsigned nel = this->nelement();
942 for (unsigned e = 0; e < nel; e++)
943 {
944 FiniteElement* fe_pt = finite_element_pt(e);
945 unsigned nnod = fe_pt->nnode();
946 for (unsigned j = 0; j < nnod; j++)
947 {
948 Node* nod_pt = fe_pt->node_pt(j);
949
950 // Extract nodal coordinates from node:
951 Vector<double> x(2);
952 x[0] = nod_pt->x(0);
953 x[1] = nod_pt->x(1);
954
955 // Evaluate the radial distance
956 double r = sqrt(x[0] * x[0] + x[1] * x[1]);
957
958 // Check
959 if (Outer_radius == 0.0)
960 {
961 throw OomphLibError("Outer radius for DtN BC must not be zero!",
962 OOMPH_CURRENT_FUNCTION,
963 OOMPH_EXCEPTION_LOCATION);
964 }
965
966 if (std::fabs((r - this->Outer_radius) / Outer_radius) >
968 {
969 std::ostringstream error_stream;
970 error_stream
971 << "Node at " << x[0] << " " << x[1] << " has radius " << r
972 << " which does not "
973 << " agree with \nspecified outer radius " << this->Outer_radius
974 << " within relative tolerance "
976 << ".\nYou can adjust the tolerance via\n"
977 << "ToleranceForFourierDecomposedHelmholtzOuterBoundary::Tol\n"
978 << "or recompile without PARANOID.\n";
979 throw OomphLibError(error_stream.str(),
980 OOMPH_CURRENT_FUNCTION,
981 OOMPH_EXCEPTION_LOCATION);
982 }
983 }
984 }
985 }
986#endif
987
988
989 // Get parameters from first element
992 this->element_pt(0));
993 double k =
994 sqrt(dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->k_squared());
995 int n_fourier_decomposed =
996 dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->fourier_wavenumber();
997 double n_hankel_order_max = double(N_terms) + 0.5;
998 double n_hankel_order_tmp = 0.0;
999
1000 /// Imaginary unit
1001 std::complex<double> I(0.0, 1.0);
1002
1003 // Precompute factors in sum
1004 Vector<std::complex<double>> h_a(N_terms + 1), hp_a(N_terms + 1),
1005 q(N_terms + 1, std::complex<double>(0.0, 0.0));
1006 Vector<double> jv(N_terms + 1), yv(N_terms + 1), djv(N_terms + 1),
1007 dyv(N_terms + 1);
1008
1009 // Get Bessel functions
1010 CRBond_Bessel::bessjyv(n_hankel_order_max,
1011 Outer_radius * k,
1012 n_hankel_order_tmp,
1013 &jv[0],
1014 &yv[0],
1015 &djv[0],
1016 &dyv[0]);
1017
1018 // Assemble Hankel functions and their derivatives
1019 for (unsigned j = 0; j < N_terms + 1; j++)
1020 {
1021 h_a[j] = jv[j] + I * yv[j];
1022 hp_a[j] = djv[j] + I * dyv[j];
1023 }
1024
1025 // Precompute relevant terms in sum
1026 for (unsigned i = n_fourier_decomposed; i < N_terms; i++)
1027 {
1028 double n_1 =
1029 Legendre_functions_helper::factorial(i - n_fourier_decomposed);
1030 double n_2 =
1031 Legendre_functions_helper::factorial(i + n_fourier_decomposed);
1032
1033 q[i] = k * sqrt(MathematicalConstants::Pi / (2.0 * k * Outer_radius)) *
1034 (hp_a[i] - h_a[i] / (2.0 * k * Outer_radius)) *
1035 (2.0 * double(i) + 1.0) /
1036 (2.0 * sqrt(MathematicalConstants::Pi / (2.0 * k * Outer_radius)) *
1037 h_a[i]) *
1038 (n_1 / n_2);
1039 }
1040
1041 // first loop over elements e
1042 unsigned nel = this->nelement();
1043 for (unsigned e = 0; e < nel; e++)
1044 {
1045 // Get a pointer to element
1048 this->element_pt(e));
1049
1050 // Set the value of n_intpt
1051 const unsigned n_intpt = el_pt->integral_pt()->nweight();
1052
1053 // initialise gamma integral and its derivatives
1054 Vector<std::complex<double>> gamma_vector(n_intpt,
1055 std::complex<double>(0.0, 0.0));
1056 Vector<std::map<unsigned, std::complex<double>>> d_gamma_vector(n_intpt);
1057
1058 // Loop over the integration points
1059 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1060 {
1061 // Allocate and initialise coordinate
1062 Vector<double> x(el_pt->dim() + 1, 0.0);
1063
1064 // Set the Vector to hold local coordinates
1065 Vector<double> s(el_pt->dim(), 0.0);
1066 for (unsigned i = 0; i < el_pt->dim(); i++)
1067 {
1068 s[i] = el_pt->integral_pt()->knot(ipt, i);
1069 }
1070
1071 // Get the coordinates of the integration point
1072 el_pt->interpolated_x(s, x);
1073
1074 // Polar angle
1075 double theta = atan2(x[0], x[1]);
1076
1077 // Elemental contribution to gamma integral and its derivative
1078 std::complex<double> gamma_con(0.0, 0.0);
1079 std::map<unsigned, std::complex<double>> d_gamma_con;
1080
1081 // loop over the Fourier terms
1082 for (unsigned nn = n_fourier_decomposed; nn < N_terms; nn++)
1083 {
1084 // Second loop over the elements
1085 // to evaluate the complete integral
1086 for (unsigned ee = 0; ee < nel; ee++)
1087 {
1089 dynamic_cast<
1091 this->element_pt(ee));
1092
1093 // contribution of the positive term in the sum
1095 theta, nn, gamma_con, d_gamma_con);
1096
1097 unsigned n_node = eel_pt->nnode();
1098
1099 gamma_vector[ipt] += q[nn] * gamma_con;
1100 for (unsigned l = 0; l < n_node; l++)
1101 {
1102 // Add the contribution of the real local data
1103 int local_unknown_real = eel_pt->nodal_local_eqn(
1104 l, eel_pt->u_index_fourier_decomposed_helmholtz().real());
1105
1106 if (local_unknown_real >= 0)
1107 {
1108 int global_eqn_real = eel_pt->eqn_number(local_unknown_real);
1109 d_gamma_vector[ipt][global_eqn_real] +=
1110 q[nn] * d_gamma_con[global_eqn_real];
1111 }
1112
1113 // Add the contribution of the imag local data
1114 int local_unknown_imag = eel_pt->nodal_local_eqn(
1115 l, eel_pt->u_index_fourier_decomposed_helmholtz().imag());
1116
1117 if (local_unknown_imag >= 0)
1118 {
1119 int global_eqn_imag = eel_pt->eqn_number(local_unknown_imag);
1120 d_gamma_vector[ipt][global_eqn_imag] +=
1121 q[nn] * d_gamma_con[global_eqn_imag];
1122 }
1123 } // end of loop over the node
1124 } // End of second loop over the elements
1125 } // End of loop over Fourier terms
1126 } // end of loop over integration point
1127
1128 // Store it in map
1129 Gamma_at_gauss_point[el_pt] = gamma_vector;
1130 D_Gamma_at_gauss_point[el_pt] = d_gamma_vector;
1131
1132 } // end of first loop over element
1133 }
1134
1135 //===========================================================================
1136 /// Constructor, takes the pointer to the "bulk" element and the face index.
1137 //===========================================================================
1138 template<class ELEMENT>
1141 const int& face_index)
1142 : FaceGeometry<ELEMENT>(), FaceElement()
1143 {
1144 // Let the bulk element build the FaceElement, i.e. setup the pointers
1145 // to its nodes (by referring to the appropriate nodes in the bulk
1146 // element), etc.
1147 bulk_el_pt->build_face_element(face_index, this);
1148
1149 // Set up U_index_fourier_decomposedhelmholtz.
1150 U_index_fourier_decomposed_helmholtz = std::complex<unsigned>(0, 1);
1151
1152 // Cast to the appropriate FourierDecomposedHelmholtzEquation so that we can
1153 // find the index at which the variable is stored
1154 // We assume that the dimension of the full problem is the same
1155 // as the dimension of the node, if this is not the case you will have
1156 // to write custom elements, sorry
1158 dynamic_cast<FourierDecomposedHelmholtzEquations*>(bulk_el_pt);
1159 if (eqn_pt == 0)
1160 {
1161 std::string error_string =
1162 "Bulk element must inherit from FourierDecomposedHelmholtzEquations.";
1163 throw OomphLibError(
1164 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1165 }
1166 // Otherwise read out the value
1167 else
1168 {
1169 // Read the index from the (cast) bulk element
1172 }
1173 }
1174} // namespace oomph
1175
1176#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
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
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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...
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...
double global_power_contribution(std::ofstream &outfile)
Compute the element's contribution to the time-averaged radiated power over the artificial boundary....
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the test functions and to return the Jacobian of mapping between local and global...
FourierDecomposedHelmholtzBCElementBase(const FourierDecomposedHelmholtzBCElementBase &dummy)=delete
Broken copy constructor.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Broken assignment operator.
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Return the index at which the real/imag unknown value is stored.
std::complex< unsigned > U_index_fourier_decomposed_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
double global_power_contribution()
Compute the element's contribution to the time-averaged radiated power over the artificial boundary.
//////////////////////////////////////////////////////////////////
FourierDecomposedHelmholtzDtNMesh< ELEMENT > * outer_boundary_mesh_pt() const
Access function to mesh of all DtN boundary condition elements (needed to get access to gamma values)
void fill_in_generic_residual_contribution_fourier_decomposed_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...
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.
FourierDecomposedHelmholtzDtNBoundaryElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Construct element from specification of bulk element and face index.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void complete_setup_of_dependencies()
Complete the setup of additional dependencies arising through the far-away interaction with other nod...
FourierDecomposedHelmholtzDtNMesh< ELEMENT > * Outer_boundary_mesh_pt
Pointer to mesh of all DtN boundary condition elements (needed to get access to gamma values)
void set_outer_boundary_mesh_pt(FourierDecomposedHelmholtzDtNMesh< ELEMENT > *mesh_pt)
Set mesh of all DtN boundary condition elements.
void compute_gamma_contribution(const double &theta, const unsigned &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....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
std::map< FiniteElement *, Vector< std::complex< double > > > Gamma_at_gauss_point
Container to store the gamma integral for given Gauss point and element.
unsigned & n_terms()
Number of terms used in the computation of the gamma integral.
FourierDecomposedHelmholtzDtNMesh(const double &outer_radius, const unsigned &n_terms)
Constructor: Specify radius of outer boundary and number of terms used in the computation of the gamm...
Vector< std::complex< double > > & gamma_at_gauss_point(FiniteElement *el_pt)
Gamma integral evaluated at Gauss points for specified element.
void setup_gamma()
Compute and store the gamma integral at all integration points of the constituent elements.
unsigned N_terms
Number of terms used in the Gamma computation.
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 ...
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Broken assignment operator.
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
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....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double plgndr2(const unsigned &l, const unsigned &m, const double &x)
Legendre polynomials depending on two parameters.
const std::complex< double > I(0.0, 1.0)
The imaginary unit.
const double Pi
50 digits from maple
double Tol
Relative tolerance to within radius of points on DtN boundary are allowed to deviate from specified v...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...