refineable_time_harmonic_linear_elasticity_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 member functions and static member data for refineable
31 //linear elasticity elements
32 
33 
35 
36 namespace oomph
37 {
38 
39 //====================================================================
40 /// Residuals for Refineable QTimeHarmonicLinearElasticityElements
41 //====================================================================
42  template<unsigned DIM>
45  Vector<double> &residuals,DenseMatrix<double> &jacobian,unsigned flag)
46  {
47 
48 #ifdef PARANOID
49 
50  //Find out how many positional dofs there are
51  unsigned n_position_type = this->nnodal_position_type();
52  if(n_position_type != 1)
53  {
54  throw OomphLibError(
55  "TimeHarmonicLinearElasticity is not yet implemented for more than one position type",
56  OOMPH_CURRENT_FUNCTION,
57  OOMPH_EXCEPTION_LOCATION);
58  }
59 
60  //Throw and error if an elasticity tensor has not been set
61  if(this->Elasticity_tensor_pt==0)
62  {
63  throw OomphLibError(
64  "No elasticity tensor set",
65  OOMPH_CURRENT_FUNCTION,
66  OOMPH_EXCEPTION_LOCATION);
67  }
68 
69 #endif
70 
71  //Find out how many nodes there are
72  unsigned n_node = this->nnode();
73 
74  //Find the indices at which the local displacements are stored
75  std::complex<unsigned> u_nodal_index[DIM];
76  for(unsigned i=0;i<DIM;i++)
77  {
78  u_nodal_index[i] = this->u_index_time_harmonic_linear_elasticity(i);
79  }
80 
81  // Square of non-dimensional frequency
82  const double omega_sq_local = this->omega_sq();
83 
84  //Set up memory for the shape functions
85  Shape psi(n_node);
86  DShape dpsidx(n_node,DIM);
87 
88  //Set the value of Nintpt -- the number of integration points
89  unsigned n_intpt = this->integral_pt()->nweight();
90 
91  //Set the vector to hold the local coordinates in the element
92  Vector<double> s(DIM);
93 
94  //Integer to store the local equation number
95  int local_eqn=0, local_unknown=0;
96 
97  // Local storage for pointers to hang_info objects
98  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
99 
100  //Loop over the integration points
101  for(unsigned ipt=0;ipt<n_intpt;ipt++)
102  {
103  //Assign the values of s
104  for(unsigned i=0;i<DIM;++i) {s[i] = this->integral_pt()->knot(ipt,i);}
105 
106  //Get the integral weight
107  double w = this->integral_pt()->weight(ipt);
108 
109  //Call the derivatives of the shape functions (and get Jacobian)
110  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
111 
112  //Storage for Eulerian coordinates (initialised to zero)
113  Vector<double> interpolated_x(DIM,0.0);
114 
115  // Displacement
117  interpolated_u(DIM,std::complex<double>(0.0,0.0));
118 
119  //Calculate interpolated values of the derivative of displacements
120  DenseMatrix<std::complex<double> > interpolated_dudx(
121  DIM,DIM,std::complex<double>(0.0,0.0));
122 
123  //Calculate displacements and derivatives
124  for(unsigned l=0;l<n_node;l++)
125  {
126  //Loop over displacement components
127  for(unsigned i=0;i<DIM;i++)
128  {
129  //Calculate the coordinates and the accelerations
130  interpolated_x[i] += this->nodal_position(l,i)*psi(l);
131 
132  //Get the nodal displacements
133  const std::complex<double> u_value =
134  std::complex<double>(nodal_value(l,u_nodal_index[i].real()),
135  nodal_value(l,u_nodal_index[i].imag()));
136 
137  interpolated_u[i]+=u_value*psi(l);
138 
139  //Loop over derivative directions
140  for(unsigned j=0;j<DIM;j++)
141  {
142  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
143  }
144  }
145  }
146 
147  //Get body force at current time
149  this->body_force(interpolated_x,b);
150 
151  //Premultiply the weights and the Jacobian
152  double W = w*J;
153 
154  //Number of master nodes and storage for the weight of the shape function
155  unsigned n_master=1; double hang_weight=1.0;
156 
157  //Loop over the test functions, nodes of the element
158  for(unsigned l=0;l<n_node;l++)
159  {
160  //Local boolean to indicate whether the node is hanging
161  bool is_node_hanging = node_pt(l)->is_hanging();
162 
163  //If the node is hanging
164  if(is_node_hanging)
165  {
166  hang_info_pt = node_pt(l)->hanging_pt();
167  //Read out number of master nodes from hanging data
168  n_master = hang_info_pt->nmaster();
169  }
170  //Otherwise the node is its own master
171  else
172  {
173  n_master = 1;
174  }
175 
176  //Loop over the master nodes
177  for(unsigned m=0;m<n_master;m++)
178  {
179  //Loop over the displacement components
180  for(unsigned a=0;a<DIM;a++)
181  {
182 
183 
184  // Real
185  //-----
186 
187  //Get the equation number
188  if(is_node_hanging)
189  {
190  //Get the equation number from the master node
191  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
192  u_nodal_index[a].real());
193  //Get the hang weight from the master node
194  hang_weight = hang_info_pt->master_weight(m);
195  }
196  //Otherwise the node is not hanging
197  else
198  {
199  local_eqn = this->nodal_local_eqn(l,u_nodal_index[a].real());
200  hang_weight = 1.0;
201  }
202 
203  /*IF it's not a boundary condition*/
204  if(local_eqn >= 0)
205  {
206  // Acceleration and body force
207  residuals[local_eqn] +=
208  (-omega_sq_local*interpolated_u[a].real()
209  -b[a].real())*psi(l)*W*hang_weight;
210 
211  // Stress term
212  for(unsigned b=0;b<DIM;b++)
213  {
214  for(unsigned c=0;c<DIM;c++)
215  {
216  for(unsigned d=0;d<DIM;d++)
217  {
218  //Add the stress terms to the residuals
219  residuals[local_eqn] +=
220  this->E(a,b,c,d)*interpolated_dudx(c,d).real()*dpsidx(l,b)*W*
221  hang_weight;
222  }
223  }
224  }
225 
226  //Jacobian entries
227  if(flag)
228  {
229  //Number of master nodes and weights
230  unsigned n_master2=1; double hang_weight2=1.0;
231  //Loop over the displacement basis functions again
232  for(unsigned l2=0;l2<n_node;l2++)
233  {
234  //Local boolean to indicate whether the node is hanging
235  bool is_node2_hanging = node_pt(l2)->is_hanging();
236 
237  //If the node is hanging
238  if(is_node2_hanging)
239  {
240  hang_info2_pt = node_pt(l2)->hanging_pt();
241  //Read out number of master nodes from hanging data
242  n_master2 = hang_info2_pt->nmaster();
243  }
244  //Otherwise the node is its own master
245  else
246  {
247  n_master2 = 1;
248  }
249 
250  //Loop over the master nodes
251  for(unsigned m2=0;m2<n_master2;m2++)
252  {
253  //Loop over the displacement components again
254  for(unsigned c=0;c<DIM;c++)
255  {
256 
257  //Get the number of the unknown
258  //If the node is hanging
259  if(is_node2_hanging)
260  {
261  //Get the equation number from the master node
262  local_unknown =
263  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
264  u_nodal_index[c].real());
265  //Get the hang weights
266  hang_weight2 = hang_info2_pt->master_weight(m2);
267  }
268  else
269  {
270  local_unknown =
271  this->nodal_local_eqn(l2,u_nodal_index[c].real());
272  hang_weight2 = 1.0;
273  }
274 
275  //If it's not pinned
276  if(local_unknown >= 0)
277  {
278  // Inertial term
279  if (a==c)
280  {
281  jacobian(local_eqn,local_unknown) -=
282  omega_sq_local*psi(l)*psi(l2)*W*
283  hang_weight*hang_weight2;
284  }
285 
286  // Stress term
287  for(unsigned b=0;b<DIM;b++)
288  {
289  for(unsigned d=0;d<DIM;d++)
290  {
291  //Add the contribution to the Jacobian matrix
292  jacobian(local_eqn,local_unknown) +=
293  this->E(a,b,c,d)*dpsidx(l2,d)*dpsidx(l,b)*W
294  *hang_weight*hang_weight2;
295  }
296  }
297  } //End of if not boundary condition
298  }
299  }
300  }
301  } //End of jacobian calculation
302 
303  } //End of if not boundary condition
304 
305 
306 
307  // Imag
308  //-----
309 
310  //Get the equation number
311  if(is_node_hanging)
312  {
313  //Get the equation number from the master node
314  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
315  u_nodal_index[a].imag());
316  //Get the hang weight from the master node
317  hang_weight = hang_info_pt->master_weight(m);
318  }
319  //Otherwise the node is not hanging
320  else
321  {
322  local_eqn = this->nodal_local_eqn(l,u_nodal_index[a].imag());
323  hang_weight = 1.0;
324  }
325 
326  /*IF it's not a boundary condition*/
327  if(local_eqn >= 0)
328  {
329  // Acceleration and body force
330  residuals[local_eqn] +=
331  (-omega_sq_local*interpolated_u[a].imag()
332  -b[a].imag())*psi(l)*W*hang_weight;
333 
334  // Stress term
335  for(unsigned b=0;b<DIM;b++)
336  {
337  for(unsigned c=0;c<DIM;c++)
338  {
339  for(unsigned d=0;d<DIM;d++)
340  {
341  //Add the stress terms to the residuals
342  residuals[local_eqn] +=
343  this->E(a,b,c,d)*interpolated_dudx(c,d).imag()*dpsidx(l,b)*W*
344  hang_weight;
345  }
346  }
347  }
348 
349  //Jacobian entries
350  if(flag)
351  {
352  //Number of master nodes and weights
353  unsigned n_master2=1; double hang_weight2=1.0;
354  //Loop over the displacement basis functions again
355  for(unsigned l2=0;l2<n_node;l2++)
356  {
357  //Local boolean to indicate whether the node is hanging
358  bool is_node2_hanging = node_pt(l2)->is_hanging();
359 
360  //If the node is hanging
361  if(is_node2_hanging)
362  {
363  hang_info2_pt = node_pt(l2)->hanging_pt();
364  //Read out number of master nodes from hanging data
365  n_master2 = hang_info2_pt->nmaster();
366  }
367  //Otherwise the node is its own master
368  else
369  {
370  n_master2 = 1;
371  }
372 
373  //Loop over the master nodes
374  for(unsigned m2=0;m2<n_master2;m2++)
375  {
376  //Loop over the displacement components again
377  for(unsigned c=0;c<DIM;c++)
378  {
379 
380  //Get the number of the unknown
381  //If the node is hanging
382  if(is_node2_hanging)
383  {
384  //Get the equation number from the master node
385  local_unknown =
386  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
387  u_nodal_index[c].imag());
388  //Get the hang weights
389  hang_weight2 = hang_info2_pt->master_weight(m2);
390  }
391  else
392  {
393  local_unknown =
394  this->nodal_local_eqn(l2,u_nodal_index[c].imag());
395  hang_weight2 = 1.0;
396  }
397 
398  //If it's not pinned
399  if(local_unknown >= 0)
400  {
401 
402  // Inertial term
403  if (a==c)
404  {
405  jacobian(local_eqn,local_unknown) -=
406  omega_sq_local*psi(l)*psi(l2)*W*
407  hang_weight*hang_weight2;
408  }
409 
410  // Stress term
411  for(unsigned b=0;b<DIM;b++)
412  {
413  for(unsigned d=0;d<DIM;d++)
414  {
415  //Add the contribution to the Jacobian matrix
416  jacobian(local_eqn,local_unknown) +=
417  this->E(a,b,c,d)*dpsidx(l2,d)*dpsidx(l,b)*W
418  *hang_weight*hang_weight2;
419  }
420  }
421  } //End of if not boundary condition
422  }
423  }
424  }
425  } //End of jacobian calculation
426 
427  } //End of if not boundary condition
428 
429 
430 
431 
432  } //End of loop over coordinate directions
433  }
434  } //End of loop over shape functions
435  } //End of loop over integration points
436  }
437 
438 
439 
440 
441 //====================================================================
442 /// Force building of required templates
443 //====================================================================
446 
447 }
cstr elem_len * i
Definition: cfortran.h:607
void fill_in_generic_contribution_to_residuals_time_harmonic_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overloaded helper function to take hanging nodes into account.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
static char t char * s
Definition: cfortran.h:572
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