pml_helmholtz_flux_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: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
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 prescribed flux
31 // boundary conditions to the Pml Helmholtz equations
32 #ifndef OOMPH_PML_HELMHOLTZ_FLUX_ELEMENTS_HEADER
33 #define OOMPH_PML_HELMHOLTZ_FLUX_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 // oomph-lib ncludes
42 #include "../generic/Qelements.h"
43 #include "math.h"
44 #include <complex>
45 
46 namespace oomph
47 {
48 
49 
50 
51 //======================================================================
52 /// \short A class for elements that allow the post-processing of
53 /// radiated power and flux on the boundaries of PMLHelmholtz elements.
54 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
55 /// policy class.
56 //======================================================================
57 template <class ELEMENT>
58 class PMLHelmholtzPowerElement : public virtual FaceGeometry<ELEMENT>,
59 public virtual FaceElement
60 {
61 
62 public:
63 
64 
65  /// \short Constructor, takes the pointer to the "bulk" element and the
66  /// index of the face to which the element is attached.
67  PMLHelmholtzPowerElement(FiniteElement* const &bulk_el_pt,
68  const int& face_index);
69 
70  ///\short Broken empty constructor
72  {
73  throw OomphLibError(
74  "Don't call empty constructor for PMLHelmholtzPowerElement",
75  OOMPH_CURRENT_FUNCTION,
76  OOMPH_EXCEPTION_LOCATION);
77  }
78 
79  /// Broken copy constructor
81  {
82  BrokenCopy::broken_copy("PMLHelmholtzPowerElement");
83  }
84 
85  /// Broken assignment operator
86 //Commented out broken assignment operator because this can lead to a conflict warning
87 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
88 //realise that two separate implementations of the broken function are the same and so,
89 //quite rightly, it shouts.
90  /*void operator=(const PMLHelmholtzPowerElement&)
91  {
92  BrokenCopy::broken_assign("PMLHelmholtzPowerElement");
93  }*/
94 
95 
96  /// \short Specify the value of nodal zeta from the face geometry
97  /// The "global" intrinsic coordinate of the element when
98  /// viewed as part of a geometric object should be given by
99  /// the FaceElement representation, by default (needed to break
100  /// indeterminacy if bulk element is SolidElement)
101  double zeta_nodal(const unsigned &n, const unsigned &k,
102  const unsigned &i) const
103  {return FaceElement::zeta_nodal(n,k,i);}
104 
105 
106 
107 
108 
109  /// \short Return the index at which the unknown value
110  /// is stored.
111  virtual inline std::complex<unsigned> u_index_helmholtz() const
112  {return std::complex<unsigned>(U_index_helmholtz.real(),
113  U_index_helmholtz.imag());}
114 
115 
116  /// \short Compute the element's contribution to the time-averaged
117  /// radiated power over the artificial boundary
119  {
120  // Dummy output file
121  std::ofstream outfile;
122  return global_power_contribution(outfile);
123  }
124 
125  /// \short Compute the element's contribution to the time-averaged
126  /// radiated power over the artificial boundary. Also output the
127  /// power density in the specified
128  ///output file if it's open.
129  double global_power_contribution(std::ofstream& outfile)
130  {
131  // pointer to the corresponding bulk element
132  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
133 
134  // Number of nodes in bulk element
135  unsigned nnode_bulk=bulk_elem_pt->nnode();
136  const unsigned n_node_local = nnode();
137 
138  //get the dim of the bulk and local nodes
139  const unsigned bulk_dim= bulk_elem_pt->dim();
140  const unsigned local_dim=this->dim();
141 
142  //Set up memory for the shape and test functions
143  Shape psi(n_node_local);
144 
145  //Set up memory for the shape functions
146  Shape psi_bulk(nnode_bulk);
147  DShape dpsi_bulk_dx(nnode_bulk,bulk_dim);
148 
149  //Set up memory for the outer unit normal
150  Vector< double > unit_normal(bulk_dim);
151 
152  //Set the value of n_intpt
153  const unsigned n_intpt = integral_pt()->nweight();
154 
155  //Set the Vector to hold local coordinates
156  Vector<double> s(local_dim);
157  double power=0.0;
158 
159  // Output?
160  if (outfile.is_open())
161  {
162  outfile << "ZONE\n";
163  }
164 
165  //Loop over the integration points
166  //--------------------------------
167  for(unsigned ipt=0;ipt<n_intpt;ipt++)
168  {
169  //Assign values of s
170  for(unsigned i=0;i<local_dim;i++)
171  {
172  s[i] = integral_pt()->knot(ipt,i);
173  }
174  //get the outer_unit_ext vector
175  this->outer_unit_normal(s,unit_normal);
176 
177  //Get the integral weight
178  double w = integral_pt()->weight(ipt);
179 
180  // Get jacobian of mapping
181  double J=J_eulerian(s);
182 
183  //Premultiply the weights and the Jacobian
184  double W = w*J;
185 
186  // Get local coordinates in bulk element by copy construction
188 
189  //Call the derivatives of the shape functions
190  //in the bulk -- must do this via s because this point
191  //is not an integration point the bulk element!
192  (void)bulk_elem_pt->dshape_eulerian(s_bulk,psi_bulk,dpsi_bulk_dx);
193  this->shape(s,psi);
194 
195  // Derivs of Eulerian coordinates w.r.t. local coordinates
196  std::complex<double> dphi_dn(0.0,0.0);
198  interpolated_dphidx(bulk_dim,std::complex<double>(0.0,0.0));
199  std::complex<double> interpolated_phi(0.0,0.0);
200  Vector<double> x(bulk_dim);
201 
202  //Calculate function value and derivatives:
203  //-----------------------------------------
204  // Loop over nodes
205  for(unsigned l=0;l<nnode_bulk;l++)
206  {
207  //Get the nodal value of the helmholtz unknown
208  const std::complex<double> phi_value(
209  bulk_elem_pt->nodal_value(l,bulk_elem_pt->u_index_helmholtz().real()),
210  bulk_elem_pt->nodal_value(l,bulk_elem_pt->u_index_helmholtz().imag()));
211 
212  //Loop over directions
213  for(unsigned i=0;i<bulk_dim;i++)
214  {
215  interpolated_dphidx[i] += phi_value*dpsi_bulk_dx(l,i);
216  }
217  } // End of loop over the bulk_nodes
218 
219  for(unsigned l=0;l<n_node_local;l++)
220  {
221  //Get the nodal value of the helmholtz unknown
222  const std::complex<double> phi_value(
223  nodal_value(l,u_index_helmholtz().real()),
224  nodal_value(l,u_index_helmholtz().imag()));
225 
226  interpolated_phi += phi_value*psi(l);
227  }
228 
229  //define dphi_dn
230  for(unsigned i=0;i<bulk_dim;i++)
231  {
232  dphi_dn += interpolated_dphidx[i]*unit_normal[i];
233  }
234 
235  // Power density
236  double integrand=0.5*
237  (interpolated_phi.real()*dphi_dn.imag()-
238  interpolated_phi.imag()*dphi_dn.real());
239 
240  // Output?
241  if (outfile.is_open())
242  {
243  interpolated_x(s,x);
244  double phi=atan2(x[1],x[0]);
245  outfile << x[0] << " "
246  << x[1] << " "
247  << phi << " "
248  << integrand << "\n";
249  }
250 
251  // ...add to integral
252  power+=integrand*W;
253  }
254 
255  return power;
256  }
257 
258 
259 
260 
261  /// \short Compute the element's contribution to the time-averaged
262  /// flux
263  std::complex<double> global_flux_contribution()
264  {
265  // Dummy output file
266  std::ofstream outfile;
267  return global_flux_contribution(outfile);
268  }
269 
270  /// \short Compute the element's contribution to the integral of
271  /// the flux over the artificial boundary. Also output the
272  /// flux in the specified
273  ///output file if it's open.
274  std::complex<double> global_flux_contribution(std::ofstream& outfile)
275  {
276  // pointer to the corresponding bulk element
277  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
278 
279  // Number of nodes in bulk element
280  unsigned nnode_bulk=bulk_elem_pt->nnode();
281  const unsigned n_node_local = nnode();
282 
283  //get the dim of the bulk and local nodes
284  const unsigned bulk_dim= bulk_elem_pt->dim();
285  const unsigned local_dim=this->dim();
286 
287  //Set up memory for the shape and test functions
288  Shape psi(n_node_local);
289 
290  //Set up memory for the shape functions
291  Shape psi_bulk(nnode_bulk);
292  DShape dpsi_bulk_dx(nnode_bulk,bulk_dim);
293 
294  //Set up memory for the outer unit normal
295  Vector<double> unit_normal(bulk_dim);
296 
297  //Set the value of n_intpt
298  const unsigned n_intpt = integral_pt()->nweight();
299 
300  //Set the Vector to hold local coordinates
301  Vector<double> s(local_dim);
302  std::complex<double> flux (0.0,0.0);
303 
304  // Output?
305  if (outfile.is_open())
306  {
307  outfile << "ZONE\n";
308  }
309 
310  //Loop over the integration points
311  //--------------------------------
312  for(unsigned ipt=0;ipt<n_intpt;ipt++)
313  {
314  //Assign values of s
315  for(unsigned i=0;i<local_dim;i++)
316  {
317  s[i] = integral_pt()->knot(ipt,i);
318  }
319  //get the outer_unit_ext vector
320  this->outer_unit_normal(s,unit_normal);
321 
322  //Get the integral weight
323  double w = integral_pt()->weight(ipt);
324 
325  // Get jacobian of mapping
326  double J=J_eulerian(s);
327 
328  //Premultiply the weights and the Jacobian
329  double W = w*J;
330 
331  // Get local coordinates in bulk element by copy construction
333 
334  //Call the derivatives of the shape functions
335  //in the bulk -- must do this via s because this point
336  //is not an integration point the bulk element!
337  (void)bulk_elem_pt->dshape_eulerian(s_bulk,psi_bulk,dpsi_bulk_dx);
338  this->shape(s,psi);
339 
340  // Derivs of Eulerian coordinates w.r.t. local coordinates
341  std::complex<double> dphi_dn(0.0,0.0);
343  interpolated_dphidx(bulk_dim,std::complex<double>(0.0,0.0));
344  Vector<double> x(bulk_dim);
345 
346  //Calculate function value and derivatives:
347  //-----------------------------------------
348  // Loop over nodes
349  for(unsigned l=0;l<nnode_bulk;l++)
350  {
351  //Get the nodal value of the helmholtz unknown
352  const std::complex<double> phi_value(
353  bulk_elem_pt->nodal_value(l,bulk_elem_pt->u_index_helmholtz().real()),
354  bulk_elem_pt->nodal_value(l,bulk_elem_pt->u_index_helmholtz().imag()));
355 
356  //Loop over directions
357  for(unsigned i=0;i<bulk_dim;i++)
358  {
359  interpolated_dphidx[i] += phi_value*dpsi_bulk_dx(l,i);
360  }
361  } // End of loop over the bulk_nodes
362 
363 
364  //define dphi_dn
365  for(unsigned i=0;i<bulk_dim;i++)
366  {
367  dphi_dn += interpolated_dphidx[i]*unit_normal[i];
368  }
369 
370  // Output?
371  if (outfile.is_open())
372  {
373  interpolated_x(s,x);
374  outfile << x[0] << " "
375  << x[1] << " "
376  << dphi_dn.real() << " "
377  << dphi_dn.imag() << "\n";
378  }
379 
380  // ...add to integral
381  flux+=dphi_dn*W;
382  }
383 
384  return flux;
385  }
386 
387 
388 protected:
389 
390 
391 
392 
393  /// \short The index at which the real and imag part of the unknown is stored
394  /// at the nodes
395  std::complex<unsigned> U_index_helmholtz;
396 
397  ///The spatial dimension of the problem
398  unsigned Dim;
399 
400 
401 };
402 
403 //////////////////////////////////////////////////////////////////////
404 //////////////////////////////////////////////////////////////////////
405 //////////////////////////////////////////////////////////////////////
406 
407 
408 
409 //===========================================================================
410 /// Constructor, takes the pointer to the "bulk" element, the
411 /// index of the fixed local coordinate and its value represented
412 /// by an integer (+/- 1), indicating that the face is located
413 /// at the max. or min. value of the "fixed" local coordinate
414 /// in the bulk element.
415 //===========================================================================
416 template<class ELEMENT>
419  const int &face_index) :
420  FaceGeometry<ELEMENT>(), FaceElement()
421  {
422 #ifdef PARANOID
423  {
424  //Check that the element is not a refineable 3d element
425  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
426  //If it's three-d
427  if(elem_pt->dim()==3)
428  {
429  //Is it refineable
430  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
431  if(ref_el_pt!=0)
432  {
433  if (this->has_hanging_nodes())
434  {
435  throw OomphLibError(
436  "This flux element will not work correctly if nodes are hanging\n",
437  OOMPH_CURRENT_FUNCTION,
438  OOMPH_EXCEPTION_LOCATION);
439  }
440  }
441  }
442  }
443 #endif
444 
445  // Let the bulk element build the FaceElement, i.e. setup the pointers
446  // to its nodes (by referring to the appropriate nodes in the bulk
447  // element), etc.
448  bulk_el_pt->build_face_element(face_index,this);
449 
450  // Extract the dimension of the problem from the dimension of
451  // the first node
452  Dim = this->node_pt(0)->ndim();
453 
454  //Set up U_index_helmholtz. Initialise to zero, which probably won't change
455  //in most cases, oh well, the price we pay for generality
456  U_index_helmholtz = std::complex<unsigned>(0,1);
457 
458  //Cast to the appropriate PMLHelmholtzEquation so that we can
459  //find the index at which the variable is stored
460  //We assume that the dimension of the full problem is the same
461  //as the dimension of the node, if this is not the case you will have
462  //to write custom elements, sorry
463  switch(Dim)
464  {
465  //One dimensional problem
466  case 1:
467  {
468  PMLHelmholtzEquations<1>* eqn_pt =
469  dynamic_cast<PMLHelmholtzEquations<1>*>(bulk_el_pt);
470  //If the cast has failed die
471  if(eqn_pt==0)
472  {
473  std::string error_string =
474  "Bulk element must inherit from PMLHelmholtzEquations.";
475  error_string +=
476  "Nodes are one dimensional, but cannot cast the bulk element to\n";
477  error_string += "PMLHelmholtzEquations<1>\n.";
478  error_string +=
479  "If you desire this functionality, you must implement it yourself\n";
480 
481  throw OomphLibError(error_string,
482  OOMPH_CURRENT_FUNCTION,
483  OOMPH_EXCEPTION_LOCATION);
484  }
485  //Otherwise read out the value
486  else
487  {
488  //Read the index from the (cast) bulk element
490  }
491  }
492  break;
493 
494  //Two dimensional problem
495  case 2:
496  {
497  PMLHelmholtzEquations<2>* eqn_pt =
498  dynamic_cast<PMLHelmholtzEquations<2>*>(bulk_el_pt);
499  //If the cast has failed die
500  if(eqn_pt==0)
501  {
502  std::string error_string =
503  "Bulk element must inherit from PMLHelmholtzEquations.";
504  error_string +=
505  "Nodes are two dimensional, but cannot cast the bulk element to\n";
506  error_string += "PMLHelmholtzEquations<2>\n.";
507  error_string +=
508  "If you desire this functionality, you must implement it yourself\n";
509 
510  throw OomphLibError(error_string,
511  OOMPH_CURRENT_FUNCTION,
512  OOMPH_EXCEPTION_LOCATION);
513  }
514  else
515  {
516  //Read the index from the (cast) bulk element
518  }
519  }
520 
521  break;
522 
523  //Three dimensional problem
524  case 3:
525  {
526  PMLHelmholtzEquations<3>* eqn_pt =
527  dynamic_cast<PMLHelmholtzEquations<3>*>(bulk_el_pt);
528  //If the cast has failed die
529  if(eqn_pt==0)
530  {
531  std::string error_string =
532  "Bulk element must inherit from PMLHelmholtzEquations.";
533  error_string +=
534  "Nodes are three dimensional, but cannot cast the bulk element to\n";
535  error_string += "PMLHelmholtzEquations<3>\n.";
536  error_string +=
537  "If you desire this functionality, you must implement it yourself\n";
538 
539  throw OomphLibError(error_string,
540  OOMPH_CURRENT_FUNCTION,
541  OOMPH_EXCEPTION_LOCATION);
542  }
543  else
544  {
545  //Read the index from the (cast) bulk element
547  }
548  }
549  break;
550 
551  //Any other case is an error
552  default:
553  std::ostringstream error_stream;
554  error_stream << "Dimension of node is " << Dim
555  << ". It should be 1,2, or 3!" << std::endl;
556 
557  throw OomphLibError(error_stream.str(),
558  OOMPH_CURRENT_FUNCTION,
559  OOMPH_EXCEPTION_LOCATION);
560  break;
561  }
562  }
563 
564 
565 
566 //////////////////////////////////////////////////////////////////////////
567 //////////////////////////////////////////////////////////////////////////
568 //////////////////////////////////////////////////////////////////////////
569 
570 
571 
572 
573 //======================================================================
574 /// \short A class for elements that allow the imposition of an
575 /// applied flux on the boundaries of PMLHelmholtz elements.
576 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
577 /// policy class.
578 //======================================================================
579 template <class ELEMENT>
580 class PMLHelmholtzFluxElement : public virtual FaceGeometry<ELEMENT>,
581 public virtual FaceElement
582 {
583 
584 public:
585 
586  /// \short Function pointer to the prescribed-flux function fct(x,f(x)) --
587  /// x is a Vector and the flux is a complex. NOTE THAT f(x) represents
588  /// c^2 du/dn!
589  typedef void (*PMLHelmholtzPrescribedFluxFctPt)
590  (const Vector<double>& x, std::complex<double>& flux);
591 
592  /// \short Constructor, takes the pointer to the "bulk" element and the
593  /// index of the face to which the element is attached.
594  PMLHelmholtzFluxElement(FiniteElement* const &bulk_el_pt,
595  const int& face_index);
596 
597  ///\short Broken empty constructor
599  {
600  throw OomphLibError(
601  "Don't call empty constructor for PMLHelmholtzFluxElement",
602  OOMPH_CURRENT_FUNCTION,
603  OOMPH_EXCEPTION_LOCATION);
604  }
605 
606  /// Broken copy constructor
608  {
609  BrokenCopy::broken_copy("PMLHelmholtzFluxElement");
610  }
611 
612  /// Broken assignment operator
613  /*void operator=(const PMLHelmholtzFluxElement&)
614  {
615  BrokenCopy::broken_assign("PMLHelmholtzFluxElement");
616  }*/
617 
618 
619  /// Access function for the prescribed-flux function pointer
621 
622 
623  /// Add the element's contribution to its residual vector
625  {
626  //Call the generic residuals function with flag set to 0
627  //using a dummy matrix argument
630  }
631 
632 
633  /// \short Add the element's contribution to its residual vector and its
634  /// Jacobian matrix
636  DenseMatrix<double> &jacobian)
637  {
638  //Call the generic routine with the flag set to 1
640  }
641 
642 
643  /// \short Specify the value of nodal zeta from the face geometry
644  /// The "global" intrinsic coordinate of the element when
645  /// viewed as part of a geometric object should be given by
646  /// the FaceElement representation, by default (needed to break
647  /// indeterminacy if bulk element is SolidElement)
648  double zeta_nodal(const unsigned &n, const unsigned &k,
649  const unsigned &i) const
650  {return FaceElement::zeta_nodal(n,k,i);}
651 
652 
653  /// Output function -- forward to broken version in FiniteElement
654  /// until somebody decides what exactly they want to plot here...
655  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
656 
657  /// \short Output function -- forward to broken version in FiniteElement
658  /// until somebody decides what exactly they want to plot here...
659  void output(std::ostream &outfile, const unsigned &n_plot)
660  {FiniteElement::output(outfile,n_plot);}
661 
662 
663  /// C-style output function -- forward to broken version in FiniteElement
664  /// until somebody decides what exactly they want to plot here...
665  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
666 
667  /// \short C-style output function -- forward to broken version in
668  /// FiniteElement until somebody decides what exactly they want to plot
669  /// here...
670  void output(FILE* file_pt, const unsigned &n_plot)
671  {FiniteElement::output(file_pt,n_plot);}
672 
673 
674  /// \short Return the index at which the unknown value
675  /// is stored.
676  virtual inline std::complex<unsigned> u_index_helmholtz() const
677  {return std::complex<unsigned>(U_index_helmholtz.real(),U_index_helmholtz.imag());}
678 
679 
680  /// \short The number of "DOF types" that degrees of freedom in this element
681  /// are sub-divided into: real and imaginary part
682  unsigned ndof_types() const
683  {
684  return 2;
685  }
686 
687  /// \short Create a list of pairs for all unknowns in this element,
688  /// so that the first entry in each pair contains the global equation
689  /// number of the unknown, while the second one contains the number
690  /// of the "DOF type" that this unknown is associated with.
691  /// (Function can obviously only be called if the equation numbering
692  /// scheme has been set up.) Real=0; Imag=1
694  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
695  {
696  // temporary pair (used to store dof lookup prior to being added to list)
697  std::pair<unsigned,unsigned> dof_lookup;
698 
699  // number of nodes
700  unsigned n_node = this->nnode();
701 
702  // loop over the nodes
703  for (unsigned n = 0; n < n_node; n++)
704  {
705  // determine local eqn number for real part
706  int local_eqn_number =
707  this->nodal_local_eqn(n,this->U_index_helmholtz.real());
708 
709  // ignore pinned values
710  if (local_eqn_number >= 0)
711  {
712  // store dof lookup in temporary pair: First entry in pair
713  // is global equation number; second entry is dof type
714  dof_lookup.first = this->eqn_number(local_eqn_number);
715  dof_lookup.second = 0;
716 
717  // add to list
718  dof_lookup_list.push_front(dof_lookup);
719  }
720 
721  // determine local eqn number for imag part
722  local_eqn_number =
723  this->nodal_local_eqn(n, this->U_index_helmholtz.imag());
724 
725  // ignore pinned values
726  if (local_eqn_number >= 0)
727  {
728  // store dof lookup in temporary pair: First entry in pair
729  // is global equation number; second entry is dof type
730  dof_lookup.first = this->eqn_number(local_eqn_number);
731  dof_lookup.second = 1;
732 
733  // add to list
734  dof_lookup_list.push_front(dof_lookup);
735  }
736  }
737  }
738 protected:
739 
740  /// \short Function to compute the shape and test functions and to return
741  /// the Jacobian of mapping between local and global (Eulerian)
742  /// coordinates
743  inline double shape_and_test(const Vector<double> &s, Shape &psi, Shape &test)
744  const
745  {
746  //Find number of nodes
747  unsigned n_node = nnode();
748 
749  //Get the shape functions
750  shape(s,psi);
751 
752  //Set the test functions to be the same as the shape functions
753  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
754 
755  //Return the value of the jacobian
756  return J_eulerian(s);
757  }
758 
759 
760  /// \short Function to compute the shape and test functions and to return
761  /// the Jacobian of mapping between local and global (Eulerian)
762  /// coordinates
763  inline double shape_and_test_at_knot(const unsigned &ipt,
764  Shape &psi, Shape &test)
765  const
766  {
767  //Find number of nodes
768  unsigned n_node = nnode();
769 
770  //Get the shape functions
771  shape_at_knot(ipt,psi);
772 
773  //Set the test functions to be the same as the shape functions
774  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
775 
776  //Return the value of the jacobian
777  return J_eulerian_at_knot(ipt);
778  }
779 
780 
781  /// Function to calculate the prescribed flux at a given spatial
782  /// position
783  void get_flux(const Vector<double>& x, std::complex<double>& flux)
784  {
785  //If the function pointer is zero return zero
786  if(Flux_fct_pt == 0)
787  {
788  flux = std::complex<double>(0.0,0.0);
789  }
790  //Otherwise call the function
791  else
792  {
793  (*Flux_fct_pt)(x,flux);
794  }
795  }
796 
797 
798  /// \short The index at which the real and imag part of the unknown is stored
799  /// at the nodes
800  std::complex<unsigned> U_index_helmholtz;
801 
802 
803  /// \short Add the element's contribution to its residual vector.
804  /// flag=1(or 0): do (or don't) compute the contribution to the
805  /// Jacobian as well.
807  Vector<double> &residuals, DenseMatrix<double> &jacobian,
808  const unsigned& flag);
809 
810 
811  /// Function pointer to the (global) prescribed-flux function
813 
814  ///The spatial dimension of the problem
815  unsigned Dim;
816 
817 
818 };
819 
820 //////////////////////////////////////////////////////////////////////
821 //////////////////////////////////////////////////////////////////////
822 //////////////////////////////////////////////////////////////////////
823 
824 
825 
826 //===========================================================================
827 /// Constructor, takes the pointer to the "bulk" element, the
828 /// index of the fixed local coordinate and its value represented
829 /// by an integer (+/- 1), indicating that the face is located
830 /// at the max. or min. value of the "fixed" local coordinate
831 /// in the bulk element.
832 //===========================================================================
833 template<class ELEMENT>
836  const int &face_index) :
837  FaceGeometry<ELEMENT>(), FaceElement()
838  {
839 #ifdef PARANOID
840  {
841  //Check that the element is not a refineable 3d element
842  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
843  //If it's three-d
844  if(elem_pt->dim()==3)
845  {
846  //Is it refineable
847  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
848  if(ref_el_pt!=0)
849  {
850  if (this->has_hanging_nodes())
851  {
852  throw OomphLibError(
853  "This flux element will not work correctly if nodes are hanging\n",
854  OOMPH_CURRENT_FUNCTION,
855  OOMPH_EXCEPTION_LOCATION);
856  }
857  }
858  }
859  }
860 #endif
861 
862  // Let the bulk element build the FaceElement, i.e. setup the pointers
863  // to its nodes (by referring to the appropriate nodes in the bulk
864  // element), etc.
865  bulk_el_pt->build_face_element(face_index,this);
866 
867  // Initialise the prescribed-flux function pointer to zero
868  Flux_fct_pt = 0;
869 
870  // Extract the dimension of the problem from the dimension of
871  // the first node
872  Dim = this->node_pt(0)->ndim();
873 
874  //Set up U_index_helmholtz. Initialise to zero, which probably won't change
875  //in most cases, oh well, the price we pay for generality
876  U_index_helmholtz = std::complex<unsigned>(0,1);
877 
878  //Cast to the appropriate PMLHelmholtzEquation so that we can
879  //find the index at which the variable is stored
880  //We assume that the dimension of the full problem is the same
881  //as the dimension of the node, if this is not the case you will have
882  //to write custom elements, sorry
883  switch(Dim)
884  {
885  //One dimensional problem
886  case 1:
887  {
888  PMLHelmholtzEquations<1>* eqn_pt =
889  dynamic_cast<PMLHelmholtzEquations<1>*>(bulk_el_pt);
890  //If the cast has failed die
891  if(eqn_pt==0)
892  {
893  std::string error_string =
894  "Bulk element must inherit from PMLHelmholtzEquations.";
895  error_string +=
896  "Nodes are one dimensional, but cannot cast the bulk element to\n";
897  error_string += "PMLHelmholtzEquations<1>\n.";
898  error_string +=
899  "If you desire this functionality, you must implement it yourself\n";
900 
901  throw OomphLibError(error_string,
902  OOMPH_CURRENT_FUNCTION,
903  OOMPH_EXCEPTION_LOCATION);
904  }
905  //Otherwise read out the value
906  else
907  {
908  //Read the index from the (cast) bulk element
910  }
911  }
912  break;
913 
914  //Two dimensional problem
915  case 2:
916  {
917  PMLHelmholtzEquations<2>* eqn_pt =
918  dynamic_cast<PMLHelmholtzEquations<2>*>(bulk_el_pt);
919  //If the cast has failed die
920  if(eqn_pt==0)
921  {
922  std::string error_string =
923  "Bulk element must inherit from PMLHelmholtzEquations.";
924  error_string +=
925  "Nodes are two dimensional, but cannot cast the bulk element to\n";
926  error_string += "PMLHelmholtzEquations<2>\n.";
927  error_string +=
928  "If you desire this functionality, you must implement it yourself\n";
929 
930  throw OomphLibError(error_string,
931  OOMPH_CURRENT_FUNCTION,
932  OOMPH_EXCEPTION_LOCATION);
933  }
934  else
935  {
936  //Read the index from the (cast) bulk element
938  }
939  }
940 
941  break;
942 
943  //Three dimensional problem
944  case 3:
945  {
946  PMLHelmholtzEquations<3>* eqn_pt =
947  dynamic_cast<PMLHelmholtzEquations<3>*>(bulk_el_pt);
948  //If the cast has failed die
949  if(eqn_pt==0)
950  {
951  std::string error_string =
952  "Bulk element must inherit from PMLHelmholtzEquations.";
953  error_string +=
954  "Nodes are three dimensional, but cannot cast the bulk element to\n";
955  error_string += "PMLHelmholtzEquations<3>\n.";
956  error_string +=
957  "If you desire this functionality, you must implement it yourself\n";
958 
959  throw OomphLibError(error_string,
960  OOMPH_CURRENT_FUNCTION,
961  OOMPH_EXCEPTION_LOCATION);
962  }
963  else
964  {
965  //Read the index from the (cast) bulk element
967  }
968  }
969  break;
970 
971  //Any other case is an error
972  default:
973  std::ostringstream error_stream;
974  error_stream << "Dimension of node is " << Dim
975  << ". It should be 1,2, or 3!" << std::endl;
976 
977  throw OomphLibError(error_stream.str(),
978  OOMPH_CURRENT_FUNCTION,
979  OOMPH_EXCEPTION_LOCATION);
980  break;
981  }
982  }
983 
984 
985 //===========================================================================
986 /// Compute the element's residual vector and the (zero) Jacobian matrix.
987 //===========================================================================
988 template<class ELEMENT>
991  Vector<double> &residuals, DenseMatrix<double> &jacobian,
992  const unsigned& flag)
993 {
994  //Find out how many nodes there are
995  const unsigned n_node = nnode();
996 
997  //Set up memory for the shape and test functions
998  Shape psif(n_node), testf(n_node);
999 
1000  //Set the value of Nintpt
1001  const unsigned n_intpt = integral_pt()->nweight();
1002 
1003  //Set the Vector to hold local coordinates
1004  Vector<double> s(Dim-1);
1005 
1006  //Integers to hold the local equation and unknown numbers
1007  int local_eqn_real=0 ,local_eqn_imag=0;
1008 
1009  //Loop over the integration points
1010  //--------------------------------
1011  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1012  {
1013 
1014  //Assign values of s
1015  for(unsigned i=0;i<(Dim-1);i++) {s[i] = integral_pt()->knot(ipt,i);}
1016 
1017  //Get the integral weight
1018  double w = integral_pt()->weight(ipt);
1019 
1020  //Find the shape and test functions and return the Jacobian
1021  //of the mapping
1022  double J = shape_and_test(s,psif,testf);
1023 
1024  //Premultiply the weights and the Jacobian
1025  double W = w*J;
1026 
1027  //Need to find position to feed into flux function, initialise to zero
1028  Vector<double> interpolated_x(Dim,0.0);
1029 
1030  //Calculate Eulerian position of integration point
1031  for(unsigned l=0;l<n_node;l++)
1032  {
1033  for(unsigned i=0;i<Dim;i++)
1034  {
1035  interpolated_x[i] += nodal_position(l,i)*psif[l];
1036  }
1037  }
1038 
1039  //Get the imposed flux
1040  std::complex<double> flux(0.0,0.0);
1041  get_flux(interpolated_x,flux);
1042 
1043  //Now add to the appropriate equations
1044  //Loop over the test functions
1045  for(unsigned l=0;l<n_node;l++)
1046  {
1047  local_eqn_real = nodal_local_eqn(l,U_index_helmholtz.real());
1048  /*IF it's not a boundary condition*/
1049  if(local_eqn_real >= 0)
1050  {
1051  //Add the prescribed flux terms
1052  residuals[local_eqn_real] -= flux.real()*testf[l]*W;
1053 
1054  // Imposed traction doesn't depend upon the solution,
1055  // --> the Jacobian is always zero, so no Jacobian
1056  // terms are required
1057  }
1058  local_eqn_imag = nodal_local_eqn(l,U_index_helmholtz.imag());
1059  /*IF it's not a boundary condition*/
1060  if(local_eqn_imag >= 0)
1061  {
1062  //Add the prescribed flux terms
1063  residuals[local_eqn_imag] -= flux.imag()*testf[l]*W;
1064 
1065  // Imposed traction doesn't depend upon the solution,
1066  // --> the Jacobian is always zero, so no Jacobian
1067  // terms are required
1068  }
1069  }
1070  }
1071 }
1072 
1073 
1074 }
1075 
1076 #endif
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual void fill_in_generic_residual_contribution_helmholtz_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the contri...
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
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
PMLHelmholtzFluxElement(const PMLHelmholtzFluxElement &dummy)
Broken copy constructor.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
std::complex< unsigned > U_index_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
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
PMLHelmholtzPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
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
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
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
A general Finite Element class.
Definition: elements.h:1271
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
A class for elements that allow the imposition of an applied flux on the boundaries of PMLHelmholtz e...
unsigned Dim
The spatial dimension of the problem.
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
void get_flux(const Vector< double > &x, std::complex< double > &flux)
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
PMLHelmholtzPowerElement(const PMLHelmholtzPowerElement &dummy)
Broken copy constructor.
double global_power_contribution(std::ofstream &outfile)
Compute the element's contribution to the time-averaged radiated power over the artificial boundary...
std::complex< double > global_flux_contribution()
Compute the element's contribution to the time-averaged flux.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
std::complex< double > global_flux_contribution(std::ofstream &outfile)
Compute the element's contribution to the integral of the flux over the artificial boundary...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Broken assignment operator.
virtual std::complex< unsigned > u_index_helmholtz() const
Return the index at which the unknown value is stored.
double global_power_contribution()
Compute the element's contribution to the time-averaged radiated power over the artificial boundary...
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.
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3158
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
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
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
PMLHelmholtzPowerElement()
Broken empty constructor.
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
PMLHelmholtzFluxElement()
Broken empty constructor.
std::complex< unsigned > U_index_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
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
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
PMLHelmholtzPrescribedFluxFctPt & flux_fct_pt()
Broken assignment operator.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
A class for elements that allow the post-processing of radiated power and flux on the boundaries of P...
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
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 ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: real and imag...
unsigned Dim
The spatial dimension of the problem.
virtual std::complex< unsigned > u_index_helmholtz() const
Return the index at which the unknown value is stored.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
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...
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
void(* PMLHelmholtzPrescribedFluxFctPt)(const Vector< double > &x, std::complex< double > &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector and the flux is a comple...