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