refineable_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//====================================================================
31 
32 
33 namespace oomph
34 {
35 
36 
37 
38 //========================================================================
39 /// Add element's contribution to the elemental
40 /// residual vector and/or Jacobian matrix.
41 /// flag=1: compute both
42 /// flag=0: compute only residual vector
43 //========================================================================
44 template<unsigned DIM>
47  DenseMatrix<double> &jacobian,
48  const unsigned& flag)
49 {
50 
51 //Find out how many nodes there are in the element
52 unsigned n_node = nnode();
53 
54 //Set up memory for the shape and test functions
55 Shape psi(n_node), test(n_node);
56 DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
57 
58 //Set the value of n_intpt
59 unsigned n_intpt = integral_pt()->nweight();
60 
61 //Integers to store the local equation and unknown numbers
62 int local_eqn_real=0, local_unknown_real=0;
63 int local_eqn_imag=0, local_unknown_imag=0;
64 
65 // Local storage for pointers to hang_info objects
66  HangInfo *hang_info_pt=0, *hang_info2_pt=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 =
76  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  // Position and gradient
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 
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  //Get the nodal value of the helmholtz unknown
99  const std::complex<double> u_value(
100  nodal_value(l,this->u_index_helmholtz().real()),
101  nodal_value(l,this->u_index_helmholtz().imag()));
102  //Add to the interpolated value
103  interpolated_u += u_value*psi(l);
104 
105  // Loop over directions
106  for(unsigned j=0;j<DIM;j++)
107  {
108  interpolated_dudx[j] += u_value*dpsidx(l,j);
109  }
110  }
111 
112  //Get body force
113  std::complex<double> source(0.0,0.0);
114  this->get_source_helmholtz(ipt,interpolated_x,source);
115 
116 
117  // Assemble residuals and Jacobian
118 
119  // Loop over the nodes for the test functions
120  for(unsigned l=0;l<n_node;l++)
121  {
122  //Local variables used to store the number of master nodes and the
123  //weight associated with the shape function if the node is hanging
124  unsigned n_master=1; double hang_weight=1.0;
125 
126  //Local bool (is the node hanging)
127  bool is_node_hanging = this->node_pt(l)->is_hanging();
128 
129  //If the node is hanging, get the number of master nodes
130  if(is_node_hanging)
131  {
132  hang_info_pt = this->node_pt(l)->hanging_pt();
133  n_master = hang_info_pt->nmaster();
134  }
135  //Otherwise there is just one master node, the node itself
136  else
137  {
138  n_master = 1;
139  }
140 
141  //Loop over the master nodes
142  for(unsigned m=0;m<n_master;m++)
143  {
144  //Get the local equation number and hang_weight
145  //If the node is hanging
146  if(is_node_hanging)
147  {
148  //Read out the local equation number from the m-th master node
149  local_eqn_real=this->local_hang_eqn(hang_info_pt->master_node_pt(m),
150  this->u_index_helmholtz().real());
151  local_eqn_imag=this->local_hang_eqn(hang_info_pt->master_node_pt(m),
152  this->u_index_helmholtz().imag());
153 
154  //Read out the weight from the master node
155  hang_weight = hang_info_pt->master_weight(m);
156  }
157  //If the node is not hanging
158  else
159  {
160  //The local equation number comes from the node itself
161  local_eqn_real=this->nodal_local_eqn(l,this->u_index_helmholtz().real());
162  local_eqn_imag=this->nodal_local_eqn(l,this->u_index_helmholtz().imag());
163 
164  //The hang weight is one
165  hang_weight = 1.0;
166  }
167 
168  //If the nodal equation is not a boundary condition
169  if(local_eqn_real >= 0)
170  {
171  // Add body force/source term here and Helmholtz bit
172  residuals[local_eqn_real] +=
173  (source.real()-this->k_squared()*interpolated_u.real())*
174  test(l)*W*hang_weight;
175 
176  // The Helmholtz bit itself
177  for(unsigned k=0;k<DIM;k++)
178  {
179  residuals[local_eqn_real] +=
180  interpolated_dudx[k].real()*dtestdx(l,k)*W*hang_weight;
181  }
182 
183  // Calculate the Jacobian
184  if(flag)
185  {
186  //Local variables to store the number of master nodes
187  //and the weights associated with each hanging node
188  unsigned n_master2=1; double hang_weight2=1.0;
189 
190  //Loop over the nodes for the variables
191  for(unsigned l2=0;l2<n_node;l2++)
192  {
193  //Local bool (is the node hanging)
194  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
195 
196  //If the node is hanging, get the number of master nodes
197  if(is_node2_hanging)
198  {
199  hang_info2_pt = this->node_pt(l2)->hanging_pt();
200  n_master2 = hang_info2_pt->nmaster();
201  }
202  //Otherwise there is one master node, the node itself
203  else
204  {
205  n_master2 = 1;
206  }
207 
208  //Loop over the master nodes
209  for(unsigned m2=0;m2<n_master2;m2++)
210  {
211  //Get the local unknown and weight
212  //If the node is hanging
213  if(is_node2_hanging)
214  {
215  //Read out the local unknown from the master node
216  local_unknown_real =
217  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
218  this->u_index_helmholtz().real());
219 
220  //Read out the hanging weight from the master node
221  hang_weight2 = hang_info2_pt->master_weight(m2);
222  }
223  //If the node is not hanging
224  else
225  {
226  //The local unknown number comes from the node
227  local_unknown_real =
228  this->nodal_local_eqn(l2,this->u_index_helmholtz().real());
229 
230  //The hang weight is one
231  hang_weight2 = 1.0;
232  }
233 
234  //If the unknown is not pinned
235  if(local_unknown_real >= 0)
236  {
237  //Add contribution to Elemental Matrix
238  for(unsigned i=0;i<DIM;i++)
239  {
240  jacobian(local_eqn_real,local_unknown_real) +=
241  dpsidx(l2,i)*dtestdx(l,i)*W*hang_weight*hang_weight2;
242  }
243 
244  // Add the helmholtz contribution
245  jacobian(local_eqn_real,local_unknown_real) -=
246  this->k_squared()*psi(l2)*test(l)*W*
247  hang_weight*hang_weight2;
248  }
249  } //End of loop over master nodes
250  } //End of loop over nodes
251  } //End of Jacobian calculation
252  } //End of case when residual equation is not pinned
253 
254 
255  if(local_eqn_imag >= 0)
256  {
257  // Add body force/source term here and Helmholtz bit
258  residuals[local_eqn_imag] +=
259  (source.imag()-this->k_squared()*interpolated_u.imag())*
260  test(l)*W*hang_weight;
261 
262  // The Helmholtz bit itself
263  for(unsigned k=0;k<DIM;k++)
264  {
265  residuals[local_eqn_imag] +=
266  interpolated_dudx[k].imag()*dtestdx(l,k)*W*hang_weight;
267  }
268 
269  // Calculate the Jacobian
270  if(flag)
271  {
272  //Local variables to store the number of master nodes
273  //and the weights associated with each hanging node
274  unsigned n_master2=1; double hang_weight2=1.0;
275 
276  //Loop over the nodes for the variables
277  for(unsigned l2=0;l2<n_node;l2++)
278  {
279  //Local bool (is the node hanging)
280  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
281 
282  //If the node is hanging, get the number of master nodes
283  if(is_node2_hanging)
284  {
285  hang_info2_pt = this->node_pt(l2)->hanging_pt();
286  n_master2 = hang_info2_pt->nmaster();
287  }
288  //Otherwise there is one master node, the node itself
289  else
290  {
291  n_master2 = 1;
292  }
293 
294  //Loop over the master nodes
295  for(unsigned m2=0;m2<n_master2;m2++)
296  {
297  //Get the local unknown and weight
298  //If the node is hanging
299  if(is_node2_hanging)
300  {
301  //Read out the local unknown from the master node
302  local_unknown_imag =
303  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
304  this->u_index_helmholtz().imag());
305 
306  //Read out the hanging weight from the master node
307  hang_weight2 = hang_info2_pt->master_weight(m2);
308  }
309  //If the node is not hanging
310  else
311  {
312  //The local unknown number comes from the node
313 
314  local_unknown_imag =
315  this->nodal_local_eqn(l2,this->u_index_helmholtz().imag());
316 
317  //The hang weight is one
318  hang_weight2 = 1.0;
319  }
320 
321  if(local_unknown_imag >= 0)
322  {
323  //Add contribution to Elemental Matrix
324  for(unsigned i=0;i<DIM;i++)
325  {
326  jacobian(local_eqn_imag,local_unknown_imag) +=
327  dpsidx(l2,i)*dtestdx(l,i)*W*hang_weight*hang_weight2;
328  }
329 
330  // Add the helmholtz contribution
331  jacobian(local_eqn_imag,local_unknown_imag) -=
332  this->k_squared()*psi(l2)*test(l)*W*
333  hang_weight*hang_weight2;
334  }
335  } //End of loop over master nodes
336  } //End of loop over nodes
337  } //End of Jacobian calculation
338  }
339  } //End of loop over master nodes for residual
340  } //End of loop over nodes
341 
342 } // End of loop over integration points
343 }
344 
345 
346 
347 
348 
349 //====================================================================
350 // Force build of templates
351 //====================================================================
352 template class RefineableQHelmholtzElement<1,2>;
353 template class RefineableQHelmholtzElement<1,3>;
354 template class RefineableQHelmholtzElement<1,4>;
355 
356 template class RefineableQHelmholtzElement<2,2>;
357 template class RefineableQHelmholtzElement<2,3>;
358 template class RefineableQHelmholtzElement<2,4>;
359 
360 template class RefineableQHelmholtzElement<3,2>;
361 template class RefineableQHelmholtzElement<3,3>;
362 template class RefineableQHelmholtzElement<3,4>;
363 
364 }
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...