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