refineable_pml_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-2024 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 #include "refineable_pml_helmholtz_elements.h"
27 
28 
29 namespace oomph
30 {
31  //======================================================================
32  /// Compute element residual Vector and/or element Jacobian matrix
33  ///
34  /// flag=1: compute both
35  /// flag=0: compute only residual Vector
36  ///
37  /// Pure version without hanging nodes
38  //======================================================================
39  template<unsigned DIM>
40  void RefineablePMLHelmholtzEquations<DIM>::
41  fill_in_generic_residual_contribution_helmholtz(
42  Vector<double>& residuals,
43  DenseMatrix<double>& jacobian,
44  const unsigned& flag)
45  {
46  // Find out how many nodes there are
47  const unsigned n_node = nnode();
48 
49  // Set up memory for the shape and test functions
50  Shape psi(n_node), test(n_node);
51  DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
52 
53  // Local storage for pointers to hang_info objects
54  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
55 
56  // Set the value of n_intpt
57  const unsigned n_intpt = integral_pt()->nweight();
58 
59  // Integers to store the local equation and unknown numbers
60  int local_eqn_real = 0, local_unknown_real = 0;
61  int local_eqn_imag = 0, local_unknown_imag = 0;
62 
63  // Loop over the integration points
64  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
65  {
66  // Get the integral weight
67  double w = integral_pt()->weight(ipt);
68 
69  // Call the derivatives of the shape and test functions
70  double J = this->dshape_and_dtest_eulerian_at_knot_helmholtz(
71  ipt, psi, dpsidx, test, dtestdx);
72 
73  // Premultiply the weights and the Jacobian
74  double W = w * J;
75 
76  // Calculate local values of unknown
77  // Allocate and initialise to zero
78  std::complex<double> interpolated_u(0.0, 0.0);
79  Vector<double> interpolated_x(DIM, 0.0);
80  Vector<std::complex<double>> interpolated_dudx(DIM);
81 
82  // Calculate function value and derivatives:
83  //-----------------------------------------
84  // Loop over nodes
85  for (unsigned l = 0; l < n_node; l++)
86  {
87  // Loop over directions
88  for (unsigned j = 0; j < DIM; j++)
89  {
90  interpolated_x[j] += nodal_position(l, j) * psi(l);
91  }
92 
93  // Get the nodal value of the helmholtz unknown
94  const std::complex<double> u_value(
95  this->nodal_value(l, this->u_index_helmholtz().real()),
96  this->nodal_value(l, this->u_index_helmholtz().imag()));
97 
98  // Add to the interpolated value
99  interpolated_u += u_value * psi(l);
100 
101  // Loop over directions
102  for (unsigned j = 0; j < DIM; j++)
103  {
104  interpolated_dudx[j] += u_value * dpsidx(l, j);
105  }
106  }
107 
108  // Get source function
109  //-------------------
110  std::complex<double> source(0.0, 0.0);
111  this->get_source_helmholtz(ipt, interpolated_x, source);
112 
113 
114  // Declare a vector of complex numbers for pml weights on the Laplace bit
115  Vector<std::complex<double>> pml_laplace_factor(DIM);
116  // Declare a complex number for pml weights on the mass matrix bit
117  std::complex<double> pml_k_squared_factor =
118  std::complex<double>(1.0, 0.0);
119 
120  // All the PML weights that participate in the assemby process
121  // are computed here. pml_laplace_factor will contain the entries
122  // for the Laplace bit, while pml_k_squared_factor contains the
123  // contributions to the Helmholtz bit. Both default to 1.0, should the PML
124  // not be enabled via enable_pml.
125  this->compute_pml_coefficients(
126  ipt, interpolated_x, pml_laplace_factor, pml_k_squared_factor);
127 
128  // Alpha adjusts the pml factors, the imaginary part produces cross terms
129  std::complex<double> alpha_pml_k_squared_factor =
130  std::complex<double>(pml_k_squared_factor.real() -
131  this->alpha() * pml_k_squared_factor.imag(),
132  this->alpha() * pml_k_squared_factor.real() +
133  pml_k_squared_factor.imag());
134 
135 
136  // std::complex<double> alpha_pml_k_squared_factor
137  // if(alpha_pt() == 0)
138  // {
139  // std::complex<double> alpha_pml_k_squared_factor =
140  // std::complex<double>(
141  // pml_k_squared_factor.real() - alpha() *
142  // pml_k_squared_factor.imag(), alpha() * pml_k_squared_factor.real() +
143  // pml_k_squared_factor.imag()
144  // );
145  // }
146  // Assemble residuals and Jacobian
147  //--------------------------------
148  // Loop over the test functions
149  for (unsigned l = 0; l < n_node; l++)
150  {
151  // Local variables used to store the number of master nodes and the
152  // weight associated with the shape function if the node is hanging
153  unsigned n_master = 1;
154  double hang_weight = 1.0;
155 
156  // Local bool (is the node hanging)
157  bool is_node_hanging = this->node_pt(l)->is_hanging();
158 
159  // If the node is hanging, get the number of master nodes
160  if (is_node_hanging)
161  {
162  hang_info_pt = this->node_pt(l)->hanging_pt();
163  n_master = hang_info_pt->nmaster();
164  }
165  // Otherwise there is just one master node, the node itself
166  else
167  {
168  n_master = 1;
169  }
170 
171  // Loop over the master nodes
172  for (unsigned m = 0; m < n_master; m++)
173  {
174  // Get the local equation number and hang_weight
175  // If the node is hanging
176  if (is_node_hanging)
177  {
178  // Read out the local equation number from the m-th master node
179  local_eqn_real =
180  this->local_hang_eqn(hang_info_pt->master_node_pt(m),
181  this->u_index_helmholtz().real());
182 
183  local_eqn_imag =
184  this->local_hang_eqn(hang_info_pt->master_node_pt(m),
185  this->u_index_helmholtz().imag());
186 
187  // Read out the weight from the master node
188  hang_weight = hang_info_pt->master_weight(m);
189  }
190  // If the node is not hanging
191  else
192  {
193  // The local equation number comes from the node itself
194  local_eqn_real =
195  this->nodal_local_eqn(l, this->u_index_helmholtz().real());
196  local_eqn_imag =
197  this->nodal_local_eqn(l, this->u_index_helmholtz().imag());
198 
199  // The hang weight is one
200  hang_weight = 1.0;
201  }
202 
203  // first, compute the real part contribution
204  //-------------------------------------------
205 
206  /*IF it's not a boundary condition*/
207  if (local_eqn_real >= 0)
208  {
209  // Add body force/source term and Helmholtz bit
210  residuals[local_eqn_real] +=
211  (source.real() - (alpha_pml_k_squared_factor.real() *
212  this->k_squared() * interpolated_u.real() -
213  alpha_pml_k_squared_factor.imag() *
214  this->k_squared() * interpolated_u.imag())) *
215  test(l) * W * hang_weight;
216 
217  // The Laplace bit
218  for (unsigned k = 0; k < DIM; k++)
219  {
220  residuals[local_eqn_real] +=
221  (pml_laplace_factor[k].real() * interpolated_dudx[k].real() -
222  pml_laplace_factor[k].imag() * interpolated_dudx[k].imag()) *
223  dtestdx(l, k) * W * hang_weight;
224  }
225 
226  // Calculate the jacobian
227  //-----------------------
228  if (flag)
229  {
230  // Local variables to store the number of master nodes
231  // and the weights associated with each hanging node
232  unsigned n_master2 = 1;
233  double hang_weight2 = 1.0;
234 
235  // Loop over the nodes for the variables
236  for (unsigned l2 = 0; l2 < n_node; l2++)
237  {
238  // Local bool (is the node hanging)
239  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
240 
241  // If the node is hanging, get the number of master nodes
242  if (is_node2_hanging)
243  {
244  hang_info2_pt = this->node_pt(l2)->hanging_pt();
245  n_master2 = hang_info2_pt->nmaster();
246  }
247  // Otherwise there is one master node, the node itself
248  else
249  {
250  n_master2 = 1;
251  }
252 
253  // Loop over the master nodes
254  for (unsigned m2 = 0; m2 < n_master2; m2++)
255  {
256  // Get the local unknown and weight
257  // If the node is hanging
258  if (is_node2_hanging)
259  {
260  // Read out the local unknown from the master node
261  local_unknown_real =
262  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
263  this->u_index_helmholtz().real());
264  local_unknown_imag =
265  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
266  this->u_index_helmholtz().imag());
267 
268  // Read out the hanging weight from the master node
269  hang_weight2 = hang_info2_pt->master_weight(m2);
270  }
271  // If the node is not hanging
272  else
273  {
274  // The local unknown number comes from the node
275  local_unknown_real = this->nodal_local_eqn(
276  l2, this->u_index_helmholtz().real());
277 
278  local_unknown_imag = this->nodal_local_eqn(
279  l2, this->u_index_helmholtz().imag());
280 
281  // The hang weight is one
282  hang_weight2 = 1.0;
283  }
284 
285 
286  // If at a non-zero degree of freedom add in the entry
287  if (local_unknown_real >= 0)
288  {
289  // Add contribution to Elemental Matrix
290  for (unsigned i = 0; i < DIM; i++)
291  {
292  jacobian(local_eqn_real, local_unknown_real) +=
293  pml_laplace_factor[i].real() * dpsidx(l2, i) *
294  dtestdx(l, i) * W * hang_weight * hang_weight2;
295  }
296  // Add the helmholtz contribution
297  jacobian(local_eqn_real, local_unknown_real) +=
298  -alpha_pml_k_squared_factor.real() * this->k_squared() *
299  psi(l2) * test(l) * W * hang_weight * hang_weight2;
300  }
301  // If at a non-zero degree of freedom add in the entry
302  if (local_unknown_imag >= 0)
303  {
304  // Add contribution to Elemental Matrix
305  for (unsigned i = 0; i < DIM; i++)
306  {
307  jacobian(local_eqn_real, local_unknown_imag) -=
308  pml_laplace_factor[i].imag() * dpsidx(l2, i) *
309  dtestdx(l, i) * W * hang_weight * hang_weight2;
310  }
311  // Add the helmholtz contribution
312  jacobian(local_eqn_real, local_unknown_imag) +=
313  alpha_pml_k_squared_factor.imag() * this->k_squared() *
314  psi(l2) * test(l) * W * hang_weight * hang_weight2;
315  }
316  }
317  }
318  }
319  }
320 
321  // Second, compute the imaginary part contribution
322  //------------------------------------------------
323 
324  /*IF it's not a boundary condition*/
325  if (local_eqn_imag >= 0)
326  {
327  // Add body force/source term and Helmholtz bit
328  residuals[local_eqn_imag] +=
329  (source.imag() - (alpha_pml_k_squared_factor.imag() *
330  this->k_squared() * interpolated_u.real() +
331  alpha_pml_k_squared_factor.real() *
332  this->k_squared() * interpolated_u.imag())) *
333  test(l) * W * hang_weight;
334 
335  // The Laplace bit
336  for (unsigned k = 0; k < DIM; k++)
337  {
338  residuals[local_eqn_imag] +=
339  (pml_laplace_factor[k].imag() * interpolated_dudx[k].real() +
340  pml_laplace_factor[k].real() * interpolated_dudx[k].imag()) *
341  dtestdx(l, k) * W * hang_weight;
342  }
343 
344  // Calculate the jacobian
345  //-----------------------
346  if (flag)
347  {
348  // Local variables to store the number of master nodes
349  // and the weights associated with each hanging node
350  unsigned n_master2 = 1;
351  double hang_weight2 = 1.0;
352 
353  // Loop over the nodes for the variables
354  for (unsigned l2 = 0; l2 < n_node; l2++)
355  {
356  // Local bool (is the node hanging)
357  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
358 
359  // If the node is hanging, get the number of master nodes
360  if (is_node2_hanging)
361  {
362  hang_info2_pt = this->node_pt(l2)->hanging_pt();
363  n_master2 = hang_info2_pt->nmaster();
364  }
365  // Otherwise there is one master node, the node itself
366  else
367  {
368  n_master2 = 1;
369  }
370 
371  // Loop over the master nodes
372  for (unsigned m2 = 0; m2 < n_master2; m2++)
373  {
374  // Get the local unknown and weight
375  // If the node is hanging
376  if (is_node2_hanging)
377  {
378  // Read out the local unknown from the master node
379  local_unknown_real =
380  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
381  this->u_index_helmholtz().real());
382  local_unknown_imag =
383  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
384  this->u_index_helmholtz().imag());
385 
386  // Read out the hanging weight from the master node
387  hang_weight2 = hang_info2_pt->master_weight(m2);
388  }
389  // If the node is not hanging
390  else
391  {
392  // The local unknown number comes from the node
393  local_unknown_real = this->nodal_local_eqn(
394  l2, this->u_index_helmholtz().real());
395 
396  local_unknown_imag = this->nodal_local_eqn(
397  l2, this->u_index_helmholtz().imag());
398 
399  // The hang weight is one
400  hang_weight2 = 1.0;
401  }
402 
403  // If at a non-zero degree of freedom add in the entry
404  if (local_unknown_imag >= 0)
405  {
406  // Add contribution to Elemental Matrix
407  for (unsigned i = 0; i < DIM; i++)
408  {
409  jacobian(local_eqn_imag, local_unknown_imag) +=
410  pml_laplace_factor[i].real() * dpsidx(l2, i) *
411  dtestdx(l, i) * W * hang_weight * hang_weight2;
412  }
413  // Add the helmholtz contribution
414  jacobian(local_eqn_imag, local_unknown_imag) +=
415  -alpha_pml_k_squared_factor.real() * this->k_squared() *
416  psi(l2) * test(l) * W * hang_weight * hang_weight2;
417  }
418  if (local_unknown_real >= 0)
419  {
420  // Add contribution to Elemental Matrix
421  for (unsigned i = 0; i < DIM; i++)
422  {
423  jacobian(local_eqn_imag, local_unknown_real) +=
424  pml_laplace_factor[i].imag() * dpsidx(l2, i) *
425  dtestdx(l, i) * W * hang_weight * hang_weight2;
426  }
427  // Add the helmholtz contribution
428  jacobian(local_eqn_imag, local_unknown_real) +=
429  -alpha_pml_k_squared_factor.imag() * this->k_squared() *
430  psi(l2) * test(l) * W * hang_weight * hang_weight2;
431  }
432  }
433  }
434  }
435  }
436  }
437  }
438 
439  } // End of loop over integration points
440  }
441 
442 
443  //====================================================================
444  // Force build of templates
445  //====================================================================
446  template class RefineableQPMLHelmholtzElement<1, 2>;
447  template class RefineableQPMLHelmholtzElement<1, 3>;
448  template class RefineableQPMLHelmholtzElement<1, 4>;
449 
450  template class RefineableQPMLHelmholtzElement<2, 2>;
451  template class RefineableQPMLHelmholtzElement<2, 3>;
452  template class RefineableQPMLHelmholtzElement<2, 4>;
453 
454  template class RefineableQPMLHelmholtzElement<3, 2>;
455  template class RefineableQPMLHelmholtzElement<3, 3>;
456  template class RefineableQPMLHelmholtzElement<3, 4>;
457 
458 } // namespace oomph
i
cstr elem_len * i
Definition: cfortran.h:603
oomph::DShape
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
oomph::HangInfo
Class that contains data for hanging nodes.
Definition: nodes.h:742
oomph::HangInfo::master_weight
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
oomph::HangInfo::master_node_pt
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::RefineablePMLHelmholtzEquations::fill_in_generic_residual_contribution_helmholtz
void fill_in_generic_residual_contribution_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
Definition: refineable_pml_helmholtz_elements.cc:41
oomph::RefineableQPMLHelmholtzElement
Refineable version of QPMLHelmholtzElement elements.
Definition: refineable_pml_helmholtz_elements.h:199
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30
refineable_pml_helmholtz_elements.h