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