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-2025 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
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",
76 }
77
78 /// Broken copy constructor
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 {
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 {
116 }
117
118 /// C-style output function -- forward to broken version in FiniteElement
119 /// until somebody decides what exactly they want to plot here...
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 {
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;
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
171
172 // Set up memory for the shape functions
175
176 // Set up memory for the outer unit normal
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
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
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!
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);
225 std::complex<double> interpolated_phi(0.0, 0.0);
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(
236 l, bulk_elem_pt->u_index_fourier_decomposed_helmholtz().real()),
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 {
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
255 }
256
257 // define dphi_dn
258 for (unsigned i = 0; i < bulk_dim; i++)
259 {
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
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,
312 DShape& dtest_ds) const
313 {
314 // Find number of nodes
315 unsigned n_node = nnode();
316
317 // Get the shape functions
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
371
372 /// Derivative of Gamma integral w.r.t global unknown, evaluated
373 /// at Gauss points for specified element
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.
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)
470
471 /// Set mesh of all DtN boundary condition elements
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 {
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
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
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
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 {
599
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
619
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
629 d_gamma[ipt][global_unk_real].real() * test[l] * r * W;
630 }
631 if (local_unknown_imag >= 0)
632 {
633 // Add the contribution
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
648
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
660 r * W;
661 }
662 if (external_unknown_imag >= 0)
663 {
664 // Add the contribution
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
690
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
700 d_gamma[ipt][global_unk_real].imag() * test[l] * r * W;
701 }
702 if (local_unknown_imag >= 0)
703 {
704 // Add the contribution
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
720
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
732 r * W;
733 }
734 if (external_unknown_imag >= 0)
735 {
736 // Add the contribution
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
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
790 DShape dpsi(n_node, 1);
791
792 // initialise the variable
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
824
825 // Derivs of Eulerian coordinates w.r.t. local coordinates
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()),
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 =
862
863 double p_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
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);
878 psi(l) * p_phi * p_theta * sin(phi) * w * dphi_ds;
879 }
880
881 // Add the contribution of the imag local data
883 l, this->U_index_fourier_decomposed_helmholtz.imag());
884 if (local_unknown_imag >= 0)
885 {
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 {
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!",
964 }
965
966 if (std::fabs((r - this->Outer_radius) / Outer_radius) >
968 {
969 std::ostringstream 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(),
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());
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,
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 =
1030 double n_2 =
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
1055 std::complex<double>(0.0, 0.0));
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
1094 eel_pt->compute_gamma_contribution(
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
1104 l, eel_pt->u_index_fourier_decomposed_helmholtz().real());
1105
1106 if (local_unknown_real >= 0)
1107 {
1111 }
1112
1113 // Add the contribution of the imag local data
1115 l, eel_pt->u_index_fourier_decomposed_helmholtz().imag());
1116
1117 if (local_unknown_imag >= 0)
1118 {
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.
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
1159 if (eqn_pt == 0)
1160 {
1161 std::string error_string =
1162 "Bulk element must inherit from FourierDecomposedHelmholtzEquations.";
1163 throw OomphLibError(
1165 }
1166 // Otherwise read out the value
1167 else
1168 {
1169 // Read the index from the (cast) bulk element
1171 eqn_pt->u_index_fourier_decomposed_helmholtz();
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
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition elements.h:4342
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition elements.h:4630
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:6037
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:4501
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:6384
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:4532
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition elements.h:4739
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:5273
FaceGeometry class definition: This policy class is used to allow construction of face elements that ...
Definition elements.h:5002
A general Finite Element class.
Definition elements.h:1317
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition elements.h:1967
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:2597
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:3054
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...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition elements.cc:3992
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:1436
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition elements.h:2615
unsigned nnode() const
Return the number of nodes.
Definition elements.h:2214
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:2321
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition elements.cc:3328
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
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:1985
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:5163
A class for elements that allow the approximation of the Sommerfeld radiation BC for Fourier decompos...
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.
FaceElement used to apply Sommerfeld radiation conditon via Dirichlet to Neumann map.
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....
================================================================= Mesh for DtN boundary condition ele...
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.
A class for all isoparametric elements that solve the Helmholtz equations.
unsigned nexternal_data() const
Return the number of external data objects.
Definition elements.h:816
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:691
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:312
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
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
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
double plgndr2(const unsigned &l, const unsigned &m, const double &x)
Legendre polynomials depending on two parameters.
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...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).