helmholtz_elements.cc
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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 //Non-inline functions for Helmholtz elements
31 #include "helmholtz_elements.h"
32 
33 
34 namespace oomph
35 {
36 
37 
38 //======================================================================
39 /// Set the data for the number of Variables at each node, always two
40 /// (real and imag part) in every case
41 //======================================================================
42  template<unsigned DIM, unsigned NNODE_1D>
43  const unsigned QHelmholtzElement<DIM,NNODE_1D>::Initial_Nvalue = 2;
44 
45 
46 //======================================================================
47 /// Compute element residual Vector and/or element Jacobian matrix
48 ///
49 /// flag=1: compute both
50 /// flag=0: compute only residual Vector
51 ///
52 /// Pure version without hanging nodes
53 //======================================================================
54 template <unsigned DIM>
57  DenseMatrix<double> &jacobian,
58  const unsigned& flag)
59 {
60  //Find out how many nodes there are
61  const unsigned n_node = nnode();
62 
63  //Set up memory for the shape and test functions
64  Shape psi(n_node), test(n_node);
65  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
66 
67  //Set the value of n_intpt
68  const unsigned n_intpt = integral_pt()->nweight();
69 
70  //Integers to store the local equation and unknown numbers
71  int local_eqn_real=0, local_unknown_real=0;
72  int local_eqn_imag=0, local_unknown_imag=0;
73 
74  //Loop over the integration points
75  for(unsigned ipt=0;ipt<n_intpt;ipt++)
76  {
77  //Get the integral weight
78  double w = integral_pt()->weight(ipt);
79 
80  //Call the derivatives of the shape and test functions
81  double J = dshape_and_dtest_eulerian_at_knot_helmholtz(ipt,psi,dpsidx,
82  test,dtestdx);
83 
84  //Premultiply the weights and the Jacobian
85  double W = w*J;
86 
87  //Calculate local values of unknown
88  //Allocate and initialise to zero
89  std::complex<double> interpolated_u(0.0,0.0);
90  Vector<double> interpolated_x(DIM,0.0);
91  Vector< std::complex<double> > interpolated_dudx(DIM);
92 
93  //Calculate function value and derivatives:
94  //-----------------------------------------
95  // Loop over nodes
96  for(unsigned l=0;l<n_node;l++)
97  {
98  // Loop over directions
99  for(unsigned j=0;j<DIM;j++)
100  {
101  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
102  }
103 
104  //Get the nodal value of the helmholtz unknown
105  const std::complex<double>
106  u_value(raw_nodal_value(l,u_index_helmholtz().real()),
107  raw_nodal_value(l,u_index_helmholtz().imag()));
108 
109  //Add to the interpolated value
110  interpolated_u += u_value*psi(l);
111 
112  // Loop over directions
113  for(unsigned j=0;j<DIM;j++)
114  {
115  interpolated_dudx[j] += u_value*dpsidx(l,j);
116  }
117  }
118 
119  //Get source function
120  //-------------------
121  std::complex<double> source(0.0,0.0);
122  get_source_helmholtz(ipt,interpolated_x,source);
123 
124  // Assemble residuals and Jacobian
125  //--------------------------------
126 
127  // Loop over the test functions
128  for(unsigned l=0;l<n_node;l++)
129  {
130 
131  // first, compute the real part contribution
132  //-------------------------------------------
133 
134  //Get the local equation
135  local_eqn_real = nodal_local_eqn(l,u_index_helmholtz().real());
136 
137  /*IF it's not a boundary condition*/
138  if(local_eqn_real >= 0)
139  {
140  // Add body force/source term and Helmholtz bit
141  residuals[local_eqn_real] +=
142  (source.real()-k_squared()*interpolated_u.real())*test(l)*W;
143 
144  // The Helmholtz bit itself
145  for(unsigned k=0;k<DIM;k++)
146  {
147  residuals[local_eqn_real] += interpolated_dudx[k].real()*dtestdx(l,k)*W;
148  }
149 
150  // Calculate the jacobian
151  //-----------------------
152  if(flag)
153  {
154  //Loop over the velocity shape functions again
155  for(unsigned l2=0;l2<n_node;l2++)
156  {
157  local_unknown_real = nodal_local_eqn(l2,u_index_helmholtz().real());
158  //If at a non-zero degree of freedom add in the entry
159  if(local_unknown_real >= 0)
160  {
161  //Add contribution to Elemental Matrix
162  for(unsigned i=0;i<DIM;i++)
163  {
164  jacobian(local_eqn_real,local_unknown_real)
165  += dpsidx(l2,i)*dtestdx(l,i)*W;
166  }
167  // Add the helmholtz contribution
168  jacobian(local_eqn_real,local_unknown_real)
169  -= k_squared()*psi(l2)*test(l)*W;
170  }//end of local_unknown
171  }
172  }
173  }
174 
175  // Second, compute the imaginary part contribution
176  //------------------------------------------------
177 
178  //Get the local equation
179  local_eqn_imag = nodal_local_eqn(l,u_index_helmholtz().imag());
180 
181  /*IF it's not a boundary condition*/
182  if(local_eqn_imag >= 0)
183  {
184  // Add body force/source term and Helmholtz bit
185  residuals[local_eqn_imag] +=
186  (source.imag()-k_squared()*interpolated_u.imag())*test(l)*W;
187 
188  // The Helmholtz bit itself
189  for(unsigned k=0;k<DIM;k++)
190  {
191  residuals[local_eqn_imag] += interpolated_dudx[k].imag()*dtestdx(l,k)*W;
192  }
193 
194  // Calculate the jacobian
195  //-----------------------
196  if(flag)
197  {
198  //Loop over the velocity shape functions again
199  for(unsigned l2=0;l2<n_node;l2++)
200  {
201  local_unknown_imag = nodal_local_eqn(l2,u_index_helmholtz().imag());
202  //If at a non-zero degree of freedom add in the entry
203  if(local_unknown_imag >= 0)
204  {
205  //Add contribution to Elemental Matrix
206  for(unsigned i=0;i<DIM;i++)
207  {
208  jacobian(local_eqn_imag,local_unknown_imag)
209  += dpsidx(l2,i)*dtestdx(l,i)*W;
210  }
211  // Add the helmholtz contribution
212  jacobian(local_eqn_imag,local_unknown_imag)
213  -= k_squared()*psi(l2)*test(l)*W;
214  }
215  }
216  }
217  }
218 
219  }
220  } // End of loop over integration points
221 }
222 
223 
224 //======================================================================
225 /// Self-test: Return 0 for OK
226 //======================================================================
227 template <unsigned DIM>
229 {
230 
231  bool passed=true;
232 
233  // Check lower-level stuff
234  if (FiniteElement::self_test()!=0)
235  {
236  passed=false;
237  }
238 
239  // Return verdict
240  if (passed)
241  {
242  return 0;
243  }
244  else
245  {
246  return 1;
247  }
248 
249 }
250 
251 
252 
253 //======================================================================
254 /// Output function:
255 ///
256 /// x,y,u_re,u_imag or x,y,z,u_re,u_imag
257 ///
258 /// nplot points in each coordinate direction
259 //======================================================================
260 template <unsigned DIM>
261 void HelmholtzEquations<DIM>::output(std::ostream &outfile,
262  const unsigned &nplot)
263 {
264 
265  //Vector of local coordinates
266  Vector<double> s(DIM);
267 
268  // Tecplot header info
269  outfile << tecplot_zone_string(nplot);
270 
271  // Loop over plot points
272  unsigned num_plot_points=nplot_points(nplot);
273  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
274  {
275 
276  // Get local coordinates of plot point
277  get_s_plot(iplot,nplot,s);
278  std::complex<double> u(interpolated_u_helmholtz(s));
279  for(unsigned i=0;i<DIM;i++)
280  {
281  outfile << interpolated_x(s,i) << " ";
282  }
283  outfile << u.real() << " " << u.imag() << std::endl;
284 
285  }
286 
287  // Write tecplot footer (e.g. FE connectivity lists)
288  write_tecplot_zone_footer(outfile,nplot);
289 
290 }
291 
292 
293 
294 
295 //======================================================================
296 /// Output function for real part of full time-dependent solution
297 ///
298 /// u = Re( (u_r +i u_i) exp(-i omega t)
299 ///
300 /// at phase angle omega t = phi.
301 ///
302 /// x,y,u or x,y,z,u
303 ///
304 /// Output at nplot points in each coordinate direction
305 //======================================================================
306 template <unsigned DIM>
307 void HelmholtzEquations<DIM>::output_real(std::ostream &outfile,
308  const double& phi,
309  const unsigned &nplot)
310 {
311 
312  //Vector of local coordinates
313  Vector<double> s(DIM);
314 
315  // Tecplot header info
316  outfile << tecplot_zone_string(nplot);
317 
318  // Loop over plot points
319  unsigned num_plot_points=nplot_points(nplot);
320  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
321  {
322 
323  // Get local coordinates of plot point
324  get_s_plot(iplot,nplot,s);
325  std::complex<double> u(interpolated_u_helmholtz(s));
326  for(unsigned i=0;i<DIM;i++)
327  {
328  outfile << interpolated_x(s,i) << " ";
329  }
330  outfile << u.real()*cos(phi)+u.imag()*sin(phi) << std::endl;
331 
332  }
333 
334  // Write tecplot footer (e.g. FE connectivity lists)
335  write_tecplot_zone_footer(outfile,nplot);
336 
337 }
338 
339 
340 //======================================================================
341 /// C-style output function:
342 ///
343 /// x,y,u or x,y,z,u
344 ///
345 /// nplot points in each coordinate direction
346 //======================================================================
347 template <unsigned DIM>
349  const unsigned &nplot)
350 {
351  //Vector of local coordinates
352  Vector<double> s(DIM);
353 
354  // Tecplot header info
355  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
356 
357  // Loop over plot points
358  unsigned num_plot_points=nplot_points(nplot);
359  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
360  {
361  // Get local coordinates of plot point
362  get_s_plot(iplot,nplot,s);
363  std::complex<double> u(interpolated_u_helmholtz(s));
364 
365  for(unsigned i=0;i<DIM;i++)
366  {
367  fprintf(file_pt,"%g ",interpolated_x(s,i));
368  }
369 
370  for(unsigned i=0;i<DIM;i++)
371  {
372  fprintf(file_pt,"%g ",interpolated_x(s,i));
373  }
374  fprintf(file_pt,"%g ",u.real());
375  fprintf(file_pt,"%g \n",u.imag());
376 
377  }
378 
379  // Write tecplot footer (e.g. FE connectivity lists)
380  write_tecplot_zone_footer(file_pt,nplot);
381 }
382 
383 
384 
385 //======================================================================
386  /// Output exact solution
387  ///
388  /// Solution is provided via function pointer.
389  /// Plot at a given number of plot points.
390  ///
391  /// x,y,u_exact or x,y,z,u_exact
392 //======================================================================
393 template <unsigned DIM>
394 void HelmholtzEquations<DIM>::output_fct(std::ostream &outfile,
395  const unsigned &nplot,
397 {
398  //Vector of local coordinates
399  Vector<double> s(DIM);
400 
401  // Vector for coordintes
402  Vector<double> x(DIM);
403 
404  // Tecplot header info
405  outfile << tecplot_zone_string(nplot);
406 
407  // Exact solution Vector
408  Vector<double> exact_soln(2);
409 
410  // Loop over plot points
411  unsigned num_plot_points=nplot_points(nplot);
412  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
413  {
414 
415  // Get local coordinates of plot point
416  get_s_plot(iplot,nplot,s);
417 
418  // Get x position as Vector
419  interpolated_x(s,x);
420 
421  // Get exact solution at this point
422  (*exact_soln_pt)(x,exact_soln);
423 
424  //Output x,y,...,u_exact
425  for(unsigned i=0;i<DIM;i++)
426  {
427  outfile << x[i] << " ";
428  }
429  outfile << exact_soln[0] << " " << exact_soln[1] << std::endl;
430  }
431 
432  // Write tecplot footer (e.g. FE connectivity lists)
433  write_tecplot_zone_footer(outfile,nplot);
434 }
435 
436 
437 
438 //======================================================================
439 /// Output function for real part of full time-dependent fct
440 ///
441 /// u = Re( (u_r +i u_i) exp(-i omega t)
442 ///
443 /// at phase angle omega t = phi.
444 ///
445 /// x,y,u or x,y,z,u
446 ///
447 /// Output at nplot points in each coordinate direction
448 //======================================================================
449 template <unsigned DIM>
451  std::ostream &outfile,
452  const double& phi,
453  const unsigned &nplot,
455 {
456  //Vector of local coordinates
457  Vector<double> s(DIM);
458 
459  // Vector for coordintes
460  Vector<double> x(DIM);
461 
462  // Tecplot header info
463  outfile << tecplot_zone_string(nplot);
464 
465  // Exact solution Vector
466  Vector<double> exact_soln(2);
467 
468  // Loop over plot points
469  unsigned num_plot_points=nplot_points(nplot);
470  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
471  {
472 
473  // Get local coordinates of plot point
474  get_s_plot(iplot,nplot,s);
475 
476  // Get x position as Vector
477  interpolated_x(s,x);
478 
479  // Get exact solution at this point
480  (*exact_soln_pt)(x,exact_soln);
481 
482  //Output x,y,...,u_exact
483  for(unsigned i=0;i<DIM;i++)
484  {
485  outfile << x[i] << " ";
486  }
487  outfile << exact_soln[0]*cos(phi)+ exact_soln[1]*sin(phi) << std::endl;
488  }
489 
490  // Write tecplot footer (e.g. FE connectivity lists)
491  write_tecplot_zone_footer(outfile,nplot);
492 }
493 
494 
495 
496 
497 //======================================================================
498  /// Validate against exact solution
499  ///
500  /// Solution is provided via function pointer.
501  /// Plot error at a given number of plot points.
502  ///
503 //======================================================================
504 template <unsigned DIM>
505 void HelmholtzEquations<DIM>::compute_error(std::ostream &outfile,
507  double& error, double& norm)
508 {
509 
510  // Initialise
511  error=0.0;
512  norm=0.0;
513 
514  //Vector of local coordinates
515  Vector<double> s(DIM);
516 
517  // Vector for coordintes
518  Vector<double> x(DIM);
519 
520  //Find out how many nodes there are in the element
521  unsigned n_node = nnode();
522 
523  Shape psi(n_node);
524 
525  //Set the value of n_intpt
526  unsigned n_intpt = integral_pt()->nweight();
527 
528  // Tecplot
529  outfile << "ZONE" << std::endl;
530 
531  // Exact solution Vector
532  Vector<double> exact_soln(2);
533 
534  //Loop over the integration points
535  for(unsigned ipt=0;ipt<n_intpt;ipt++)
536  {
537 
538  //Assign values of s
539  for(unsigned i=0;i<DIM;i++)
540  {
541  s[i] = integral_pt()->knot(ipt,i);
542  }
543 
544  //Get the integral weight
545  double w = integral_pt()->weight(ipt);
546 
547  // Get jacobian of mapping
548  double J=J_eulerian(s);
549 
550  //Premultiply the weights and the Jacobian
551  double W = w*J;
552 
553  // Get x position as Vector
554  interpolated_x(s,x);
555 
556  // Get FE function value
557  std::complex<double> u_fe=interpolated_u_helmholtz(s);
558 
559  // Get exact solution at this point
560  (*exact_soln_pt)(x,exact_soln);
561 
562  //Output x,y,...,error
563  for(unsigned i=0;i<DIM;i++)
564  {
565  outfile << x[i] << " ";
566  }
567  outfile << exact_soln[0] << " " << exact_soln[1] << " "
568  << exact_soln[0]-u_fe.real() << " " << exact_soln[1]-u_fe.imag()
569  << std::endl;
570 
571  // Add to error and norm
572  norm+=(exact_soln[0]*exact_soln[0]+exact_soln[1]*exact_soln[1])*W;
573  error+=( (exact_soln[0]-u_fe.real())*(exact_soln[0]-u_fe.real())+
574  (exact_soln[1]-u_fe.imag())*(exact_soln[1]-u_fe.imag()) )*W;
575 
576  }
577 }
578 
579 
580 
581 
582 
583 //====================================================================
584 // Force build of templates
585 //====================================================================
586 template class HelmholtzEquations<1>;
587 template class HelmholtzEquations<2>;
588 template class HelmholtzEquations<3>;
589 
590 template class QHelmholtzElement<1,2>;
591 template class QHelmholtzElement<1,3>;
592 template class QHelmholtzElement<1,4>;
593 
594 template class QHelmholtzElement<2,2>;
595 template class QHelmholtzElement<2,3>;
596 template class QHelmholtzElement<2,4>;
597 
598 template class QHelmholtzElement<3,2>;
599 template class QHelmholtzElement<3,3>;
600 template class QHelmholtzElement<3,4>;
601 
602 }
cstr elem_len * i
Definition: cfortran.h:607
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
void output_real_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for real part of full time-dependent fct u = Re( (u_r +i u_i) exp(-i omega t) at phas...
virtual void fill_in_generic_residual_contribution_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
unsigned self_test()
Self-test: Return 0 for OK.
static char t char * s
Definition: cfortran.h:572
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
void output_real(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for real part of full time-dependent solution u = Re( (u_r +i u_i) exp(-i omega t) at...
void output(std::ostream &outfile)
Output with default number of plot points.