time_harmonic_fourier_decomposed_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 functions for elements that solve the equations of linear
31 // elasticity in cartesian coordinates
32 
34 
35 
36 namespace oomph
37 {
38 
39  /// Static default value for square of angular frequency: Zero
40  std::complex<double>
43 
44  /// Static default value for Young's modulus (1.0 -- for natural
45  /// scaling, i.e. all stresses have been non-dimensionalised by
46  /// the same reference Young's modulus. Setting the "non-dimensional"
47  /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
48  /// to a number larger than one means that the material is stiffer
49  /// than assumed in the non-dimensionalisation.
50  std::complex<double>
53 
54 
55 //////////////////////////////////////////////////////////////////
56 //////////////////////////////////////////////////////////////////
57 //////////////////////////////////////////////////////////////////
58 
59 
60 //=======================================================================
61 /// Compute norm of the solution
62 //=======================================================================
64 compute_norm(double& norm)
65  {
66 
67  // Initialise
68  norm=0.0;
69 
70  //Vector of local coordinates
71  Vector<double> s(2);
72 
73  // Vector for coordintes
74  Vector<double> x(2);
75 
76  // Displacement vector
78 
79  //Find out how many nodes there are in the element
80  unsigned n_node = this->nnode();
81 
82  Shape psi(n_node);
83 
84  //Set the value of n_intpt
85  unsigned n_intpt = this->integral_pt()->nweight();
86 
87  //Loop over the integration points
88  for(unsigned ipt=0;ipt<n_intpt;ipt++)
89  {
90 
91  //Assign values of s
92  for(unsigned i=0;i<2;i++)
93  {
94  s[i] = this->integral_pt()->knot(ipt,i);
95  }
96 
97  //Get the integral weight
98  double w = this->integral_pt()->weight(ipt);
99 
100  // Get jacobian of mapping
101  double J=this->J_eulerian(s);
102 
103  //Premultiply the weights and the Jacobian
104  double W = w*J;
105 
106  // Get FE function value
108  (s,disp);
109 
110  // Add to norm
111  for (unsigned ii=0;ii<3;ii++)
112  {
113  norm+=(disp[ii].real()*disp[ii].real()+disp[ii].imag()*disp[ii].imag())*W;
114  }
115  }
116 }
117 
118 //=======================================================================
119 /// Get strain tensor
120 //=======================================================================
122  get_strain(const Vector<double>& s, DenseMatrix<std::complex<double> >& strain)
123  {
124  //Find out how many nodes there are
125  unsigned n_node = this->nnode();
126 
127 #ifdef PARANOID
128  //Find out how many positional dofs there are
129  unsigned n_position_type = this->nnodal_position_type();
130 
131  if(n_position_type != 1)
132  {
133  throw OomphLibError(
134  "TimeHarmonicFourierDecomposedLinearElasticity is not yet implemented for more than one position type",
135  OOMPH_CURRENT_FUNCTION,
136  OOMPH_EXCEPTION_LOCATION);
137  }
138 #endif
139 
140  //Find the indices at which the local displacements are stored
141  std::complex<unsigned> u_nodal_index[3];
142  for(unsigned i=0;i<3;i++)
143  {
144  u_nodal_index[i] = this->
146  }
147 
148  // Fourier wavenumber as double
149  double n = double(this->fourier_wavenumber());
150 
151  //Set up memory for the shape functions
152  Shape psi(n_node);
153 
154  // Derivs only w.r.t. r [0] and z [1]
155  DShape dpsidx(n_node,2);
156 
157  //Call the derivatives of the shape functions (ignore Jacobian)
158  this->dshape_eulerian(s,psi,dpsidx);
159 
160  //Storage for Eulerian coordinates (r,z; initialised to zero)
162 
163  // Displacements u_r,u_z,u_theta
165  interpolated_u(3,std::complex<double>(0.0,0.0));
166 
167  //Calculate interpolated values of the derivatives w.r.t.
168  //Eulerian coordinates
170  interpolated_dudx(3,2,std::complex<double>(0.0,0.0));
171 
172  //Calculate displacements and derivatives
173  for(unsigned l=0;l<n_node;l++)
174  {
175  //Calculate the coordinates
176  for(unsigned i=0;i<2;i++)
177  {
178  interpolated_x[i] += this->nodal_position(l,i)*psi(l);
179  }
180  //Get the nodal displacements
181  for(unsigned i=0;i<3;i++)
182  {
183  const std::complex<double> u_value
184  = std::complex<double>(
185  this->nodal_value(l,u_nodal_index[i].real()),
186  this->nodal_value(l,u_nodal_index[i].imag()));
187 
188  interpolated_u[i]+=u_value*psi(l);
189 
190  //Loop over derivative directions
191  for(unsigned j=0;j<2;j++)
192  {
193  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
194  }
195  }
196  }
197 
198  // define shorthand notation for regularly-occurring terms
199  double r = interpolated_x[0];
200  const std::complex<double> I(0.0,1.0);
201 
202  // r component of displacement
203  std::complex<double> ur = interpolated_u[0];
204 
205  // z component of displacement
206  std::complex<double> uz = interpolated_u[1];
207 
208  // theta component of displacement
209  std::complex<double> uth = interpolated_u[2];
210 
211  // derivatives w.r.t. r and z:
212  std::complex<double> durdr = interpolated_dudx(0,0);
213  std::complex<double> durdz = interpolated_dudx(0,1);
214  std::complex<double> duzdr = interpolated_dudx(1,0);
215  std::complex<double> duzdz = interpolated_dudx(1,1);
216  std::complex<double> duthdr = interpolated_dudx(2,0);
217  std::complex<double> duthdz = interpolated_dudx(2,1);
218 
219  strain(0,0)=durdr;
220  strain(2,2)=I*double(n)*uth/r+ur/r;
221  strain(1,1)=duzdz;
222  strain(0,2)=0.5*(I*double(n)*ur/r-uth/r+duthdr);
223  strain(2,0)=0.5*(I*double(n)*ur/r-uth/r+duthdr);
224  strain(0,1)=0.5*(durdz+duzdr);
225  strain(1,0)=0.5*(durdz+duzdr);
226  strain(2,1)=0.5*(duthdz+I*double(n)*uz/r);
227  strain(1,2)=0.5*(duthdz+I*double(n)*uz/r);
228  }
229 
230 ///////////////////////////////////////////////////////////////
231 ///////////////////////////////////////////////////////////////
232 ///////////////////////////////////////////////////////////////
233 
234 //=======================================================================
235 /// Compute the residuals for the Fourier decomposed (in cyl. polars)
236 /// time harmonic linear elasticity equations in. Flag indicates if
237 /// we want Jacobian too.
238 //=======================================================================
241  Vector<double> &residuals, DenseMatrix<double> &jacobian, unsigned flag)
242  {
243  //Find out how many nodes there are
244  unsigned n_node = this->nnode();
245 
246 #ifdef PARANOID
247  //Find out how many positional dofs there are
248  unsigned n_position_type = this->nnodal_position_type();
249 
250  if(n_position_type != 1)
251  {
252  throw OomphLibError(
253  "TimeHarmonicFourierDecomposedLinearElasticity is not yet implemented for more than one position type",
254  OOMPH_CURRENT_FUNCTION,
255  OOMPH_EXCEPTION_LOCATION);
256  }
257 #endif
258 
259  //Find the indices at which the local displacements are stored
260  std::complex<unsigned> u_nodal_index[3];
261  for(unsigned i=0;i<3;i++)
262  {
263  u_nodal_index[i] = this->
265  }
266 
267  // Get (complex) elastic parameters
268  std::complex<double> nu_local = this->nu();
269  std::complex<double> youngs_modulus_local = this->youngs_modulus();
270 
271  // Obtain Lame parameters from Young's modulus and Poisson's ratio
272  std::complex<double> lambda =
273  youngs_modulus_local*nu_local/(1.0+nu_local)/(1.0-2.0*nu_local);
274  std::complex<double> mu = youngs_modulus_local/2.0/(1.0+nu_local);
275 
276 
277  // Fourier wavenumber as double
278  double n = double(this->fourier_wavenumber());
279 
280  // Square of non-dimensional frequency
281  const std::complex<double> omega_sq = this->omega_sq();
282 
283  //Set up memory for the shape functions
284  Shape psi(n_node);
285 
286  // Derivs only w.r.t. r [0] and z [1]
287  DShape dpsidx(n_node,2);
288 
289  //Set the value of Nintpt -- the number of integration points
290  unsigned n_intpt = this->integral_pt()->nweight();
291 
292  //Set the vector to hold the local coordinates in the element
293  Vector<double> s(2);
294 
295  //Integers to store the local equation numbers
296  int local_eqn_real=0, local_eqn_imag=0,
297  local_unknown_real=0, local_unknown_imag=0;
298 
299  //Loop over the integration points
300  for(unsigned ipt=0;ipt<n_intpt;ipt++)
301  {
302  //Assign the values of s
303  for(unsigned i=0;i<2;++i)
304  {
305  s[i] = this->integral_pt()->knot(ipt,i);
306  }
307 
308  //Get the integral weight
309  double w = this->integral_pt()->weight(ipt);
310 
311  //Call the derivatives of the shape functions (and get Jacobian)
312  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
313 
314  //Storage for Eulerian coordinates (r,z; initialised to zero)
316 
317  // Displacements u_r,u_z,u_theta
319  interpolated_u(3,std::complex<double>(0.0,0.0));
320 
321  //Calculate interpolated values of the derivatives w.r.t.
322  //Eulerian coordinates
324  interpolated_dudx(3,2,std::complex<double>(0.0,0.0));
325 
326  //Calculate displacements and derivatives
327  for(unsigned l=0;l<n_node;l++)
328  {
329  //Calculate the coordinates
330  for(unsigned i=0;i<2;i++)
331  {
332  interpolated_x[i] += this->raw_nodal_position(l,i)*psi(l);
333  }
334  //Get the nodal displacements
335  for(unsigned i=0;i<3;i++)
336  {
337  const std::complex<double> u_value
338  = std::complex<double>(
339  this->raw_nodal_value(l,u_nodal_index[i].real()),
340  this->raw_nodal_value(l,u_nodal_index[i].imag()));
341 
342  interpolated_u[i]+=u_value*psi(l);
343 
344  //Loop over derivative directions
345  for(unsigned j=0;j<2;j++)
346  {
347  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
348  }
349  }
350  }
351 
352  //Get body force
354  this->body_force(interpolated_x,b);
355 
356  //Premultiply the weights and the Jacobian
357  double W = w*J;
358 
359 //=====EQUATIONS OF FOURIER DECOMPOSED TIME HARMONIC LINEAR ELASTICITY ========
360 
361  // define shorthand notation for regularly-occurring terms
362  double r = interpolated_x[0];
363  double rsq = pow(r,2);
364  double nsq = pow(n,2);
365  const std::complex<double> I(0.0,1.0);
366 
367  // r component of displacement
368  std::complex<double> ur = interpolated_u[0];
369 
370  // z component of displacement
371  std::complex<double> uz = interpolated_u[1];
372 
373  // theta component of displacement
374  std::complex<double> uth = interpolated_u[2];
375 
376  // derivatives w.r.t. r and z:
377  std::complex<double> durdr = interpolated_dudx(0,0);
378  std::complex<double> durdz = interpolated_dudx(0,1);
379  std::complex<double> duzdr = interpolated_dudx(1,0);
380  std::complex<double> duzdz = interpolated_dudx(1,1);
381  std::complex<double> duthdr = interpolated_dudx(2,0);
382  std::complex<double> duthdz = interpolated_dudx(2,1);
383 
384  // storage for (complex) terms required for analytic Jacobian
385  std::complex<double> G_r, G_z, G_theta;
386 
387  //Loop over the test functions, nodes of the element
388  for(unsigned l=0;l<n_node;l++)
389  {
390  //Loop over the displacement components
391  for(unsigned a=0;a<3;a++)
392  {
393  //Get the REAL equation number
394  local_eqn_real = this->nodal_local_eqn(l,u_nodal_index[a].real());
395 
396  /*IF it's not a boundary condition*/
397  if(local_eqn_real >= 0)
398  {
399  // Acceleration and body force
400  residuals[local_eqn_real] +=
401  (-omega_sq.real()*interpolated_u[a].real()
402  +omega_sq.imag()*interpolated_u[a].imag()
403  -b[a].real())*psi(l)*r*W;
404 
405  // Three components of the stress divergence term:
406  // a=0: r; a=1: z; a=2: theta
407 
408  // Real part of r-equation
409  if(a==0)
410  {
411  residuals[local_eqn_real] +=
412  (
413  mu.real()*(
414  2.0*durdr.real()*dpsidx(l,0)+
415  n*psi(l)/r*(n*ur.real()/r+
416  duthdr.imag()-uth.imag()/r)
417  +dpsidx(l,1)*(durdz.real()+duzdr.real())+
418  2.0*psi(l)/pow(r,2)*(ur.real()-n*uth.imag())
419  )
420  +mu.imag()*(
421  -2.0*durdr.imag()*dpsidx(l,0)+
422  n*psi(l)/r*(-n*ur.imag()/r+duthdr.real()-uth.real()/r)
423  -dpsidx(l,1)*(durdz.imag()+duzdr.imag())-2.0*psi(l)/pow(r,2)*
424  (ur.imag()+n*uth.real())
425  )
426  +lambda.real()*(
427  durdr.real()+ur.real()/r-n/r*uth.imag()+duzdz.real()
428  )*(dpsidx(l,0)+psi(l)/r)
429  -lambda.imag()*(
430  durdr.imag()+ur.imag()/r+n/r*uth.real()+duzdz.imag()
431  )*(dpsidx(l,0)+psi(l)/r)
432  )*r*W;
433  }
434  // Real part of z-equation
435  else if(a==1)
436  {
437  residuals[local_eqn_real] +=
438  (
439  mu.real()*(
440  dpsidx(l,0)*(durdz.real()+duzdr.real())
441  +n*psi(l)/r*(n*uz.real()/r+duthdz.imag())+
442  2.0*duzdz.real()*dpsidx(l,1)
443  )
444  +mu.imag()*(
445  -dpsidx(l,0)*(durdz.imag()+duzdr.imag())
446  +n*psi(l)/r*(-n*uz.imag()/r+duthdz.real())-
447  2.0*duzdz.imag()*dpsidx(l,1)
448  )
449  +lambda.real()*(
450  durdr.real()+ur.real()/r-n/r*uth.imag()+duzdz.real()
451  )*dpsidx(l,1)
452  -lambda.imag()*(
453  durdr.imag()+ur.imag()/r+n/r*uth.real()+duzdz.imag()
454  )*dpsidx(l,1)
455  )*r*W;
456  }
457  // Real part of theta-equation
458  else if(a==2)
459  {
460  residuals[local_eqn_real] +=
461  (
462  mu.real()*(
463  (duthdr.real()-uth.real()/r-n*ur.imag()/r)*
464  (dpsidx(l,0)-psi(l)/r)
465  +2.0*n*psi(l)/pow(r,2)*(n*uth.real()+ur.imag())+
466  dpsidx(l,1)*(duthdz.real()-n/r*uz.imag())
467  )
468  +mu.imag()*(
469  (-duthdr.imag()+uth.imag()/r-n*ur.real()/r)*
470  (dpsidx(l,0)-psi(l)/r)
471  +2.0*n*psi(l)/pow(r,2)*(-n*uth.imag()+ur.real())+
472  dpsidx(l,1)*(-duthdz.imag()-n/r*uz.real())
473  )
474  +lambda.real()*(
475  durdr.imag()+ur.imag()/r+n/r*uth.real()+duzdz.imag()
476  )*n*psi(l)/r
477  +lambda.imag()*(
478  durdr.real()+ur.real()/r-n/r*uth.imag()+duzdz.real()
479  )*n*psi(l)/r
480  )*r*W;
481  }
482  // error: a should be 0, 1 or 2
483  else
484  {
485  throw OomphLibError(
486  "a should equal 0, 1 or 2",
487  OOMPH_CURRENT_FUNCTION,
488  OOMPH_EXCEPTION_LOCATION);
489  }
490 
491  //Jacobian entries
492  if(flag)
493  {
494  //Loop over the displacement basis functions again
495  for(unsigned l2=0;l2<n_node;l2++)
496  {
497  // define complex terms used to obtain entries of current row in the Jacobian:
498 
499  // terms for rows of Jacobian matrix corresponding to r-equation
500  if(a==0)
501  {
502  G_r = (mu*(2.0*dpsidx(l2,0)*dpsidx(l,0)+
503  (nsq+2.0)/rsq*psi(l2)*psi(l)+
504  dpsidx(l2,1)*dpsidx(l,1))+
505  lambda*(dpsidx(l2,0)+psi(l2)/r)*
506  (dpsidx(l,0)+psi(l)/r)-omega_sq*psi(l2)*psi(l))*r*W;
507 
508  G_z=(mu*dpsidx(l2,0)*dpsidx(l,1)+
509  lambda*dpsidx(l2,1)*(dpsidx(l,0)+psi(l)/r))*r*W;
510 
511  G_theta=(-I*(mu*(n/r*dpsidx(l2,0)*psi(l)-3.0*n/rsq*psi(l2)*psi(l))-
512  lambda*n/r*psi(l2)*(dpsidx(l,0)+psi(l)/r)))*r*W;
513  }
514  // terms for rows of Jacobian matrix corresponding to z-equation
515  else if(a==1)
516  {
517  G_r=(mu*dpsidx(l2,1)*dpsidx(l,0)+
518  lambda*(dpsidx(l2,0)+psi(l2)/r)*dpsidx(l,1))*r*W;
519 
520  G_z=(mu*(dpsidx(l2,0)*dpsidx(l,0)+nsq/rsq*psi(l2)*psi(l)
521  +2.0*dpsidx(l2,1)*dpsidx(l,1))+
522  lambda*dpsidx(l2,1)*dpsidx(l,1)-omega_sq*psi(l2)*psi(l))*r*W;
523 
524  G_theta=(-I*(mu*n/r*dpsidx(l2,1)*psi(l)-
525  lambda*n/r*dpsidx(l,1)*psi(l2)))*r*W;
526  }
527  // terms for rows of Jacobian matrix corresponding to theta-equation
528  else if(a==2)
529  {
530  G_r=(-I*(mu*(-n/r*psi(l2)*dpsidx(l,0)+3.0*n/rsq*psi(l2)*psi(l))+
531  lambda*n/r*(dpsidx(l2,0)+psi(l2)/r)*psi(l)))*r*W;
532 
533  G_z=(-I*(-mu*n/r*psi(l2)*dpsidx(l,1)+
534  lambda*n/r*dpsidx(l2,1)*psi(l)))*r*W;
535 
536  G_theta=(mu*((dpsidx(l2,0)-psi(l2)/r)*(dpsidx(l,0)-psi(l)/r)+
537  2.0*nsq/rsq*psi(l2)*psi(l)+dpsidx(l2,1)*dpsidx(l,1))+
538  lambda*nsq/rsq*psi(l2)*psi(l)-omega_sq*psi(l2)*psi(l))*r*W;
539  }
540 
541  //Loop over the displacement components
542  for(unsigned c=0;c<3;c++)
543  {
544  // Get real and imaginary parts of local unknown
545  local_unknown_real=this->nodal_local_eqn(l2,u_nodal_index[c].real());
546  local_unknown_imag=this->nodal_local_eqn(l2,u_nodal_index[c].imag());
547 
548  //If the real part of the local unknown is not pinned
549  if(local_unknown_real >= 0)
550  {
551  if(c==0)
552  {
553  jacobian(local_eqn_real,local_unknown_real) += G_r.real();
554  }
555  else if(c==1)
556  {
557  jacobian(local_eqn_real,local_unknown_real) += G_z.real();
558  }
559  else if(c==2)
560  {
561  jacobian(local_eqn_real,local_unknown_real) += G_theta.real();
562  }
563  }
564 
565  //If the imaginary part of the local unknown is not pinned
566  if(local_unknown_imag >= 0)
567  {
568  if(c==0)
569  {
570  jacobian(local_eqn_real,local_unknown_imag) += -G_r.imag();
571  }
572  else if(c==1)
573  {
574  jacobian(local_eqn_real,local_unknown_imag) += -G_z.imag();
575  }
576  else if(c==2)
577  {
578  jacobian(local_eqn_real,local_unknown_imag) += -G_theta.imag();
579  }
580  }
581  }
582  }
583  } //End of jacobian calculation
584 
585  } //End of if not boundary condition for real eqn
586 
587 
588 
589  //Get the IMAG equation number
590  local_eqn_imag = this->nodal_local_eqn(l,u_nodal_index[a].imag());
591 
592  /*IF it's not a boundary condition*/
593  if(local_eqn_imag >= 0)
594  {
595  // Acceleration and body force
596  residuals[local_eqn_imag] +=
597  (-omega_sq.real()*interpolated_u[a].imag()
598  -omega_sq.imag()*interpolated_u[a].real()
599  -b[a].imag())*psi(l)*r*W;
600 
601  // Three components of the stress divergence term:
602  // a=0: r; a=1: z; a=2: theta
603 
604  // Imag part of r-equation
605  if(a==0)
606  {
607  residuals[local_eqn_imag] +=
608  (
609  mu.real()*(
610  2*durdr.imag()*dpsidx(l,0)+
611  n*psi(l)/r*(n*ur.imag()/r-duthdr.real()+uth.real()/r)
612  +dpsidx(l,1)*(durdz.imag()+duzdr.imag())+
613  2*psi(l)/pow(r,2)*(ur.imag()+n*uth.real())
614  )
615  +mu.imag()*(
616  2.0*durdr.real()*dpsidx(l,0)+
617  n*psi(l)/r*(n*ur.real()/r+duthdr.imag()-uth.imag()/r)
618  +dpsidx(l,1)*(durdz.real()+duzdr.real())+
619  2.0*psi(l)/pow(r,2)*(ur.real()-n*uth.imag())
620  )
621  +lambda.real()*(
622  durdr.imag()+ur.imag()/r+n/r*uth.real()+duzdz.imag()
623  )*(dpsidx(l,0)+psi(l)/r)
624  +lambda.imag()*(
625  durdr.real()+ur.real()/r-n/r*uth.imag()+duzdz.real()
626  )*(dpsidx(l,0)+psi(l)/r)
627  )*r*W;
628  }
629  // Imag part of z-equation
630  else if(a==1)
631  {
632  residuals[local_eqn_imag] +=
633  (
634  mu.real()*(
635  dpsidx(l,0)*(durdz.imag()+duzdr.imag())+
636  n*psi(l)/r*(n*uz.imag()/r-duthdz.real())+
637  2*duzdz.imag()*dpsidx(l,1)
638  )
639  +mu.imag()*(
640  dpsidx(l,0)*(durdz.real()+duzdr.real())+
641  n*psi(l)/r*(n*uz.real()/r+duthdz.imag())+
642  2.0*duzdz.real()*dpsidx(l,1)
643  )
644  +lambda.real()*(
645  durdr.imag()+ur.imag()/r+n/r*uth.real()+duzdz.imag()
646  )*dpsidx(l,1)
647  +lambda.imag()*(
648  durdr.real()+ur.real()/r-n/r*uth.imag()+duzdz.real()
649  )*dpsidx(l,1)
650  )*r*W;
651  }
652  // Imag part of theta-equation
653  else if(a==2)
654  {
655  residuals[local_eqn_imag] +=
656  (
657  mu.real()*(
658  (duthdr.imag()-uth.imag()/r+n*ur.real()/r)*
659  (dpsidx(l,0)-psi(l)/r)
660  +2.0*n*psi(l)/pow(r,2.)*(n*uth.imag()-ur.real())+
661  dpsidx(l,1)*(duthdz.imag()+n/r*uz.real())
662  )
663  +mu.imag()*(
664  (duthdr.real()-uth.real()/r-n*ur.imag()/r)*
665  (dpsidx(l,0)-psi(l)/r)
666  +2.0*n*psi(l)/pow(r,2)*(n*uth.real()+ur.imag())+
667  dpsidx(l,1)*(duthdz.real()-n/r*uz.imag())
668  )
669  +lambda.real()*(
670  -durdr.real()-ur.real()/r+n/r*uth.imag()-duzdz.real()
671  )*n*psi(l)/r
672  +lambda.imag()*(
673  durdr.imag()+ur.imag()/r+n/r*uth.real()+duzdz.imag()
674  )*n*psi(l)/r
675  )*r*W;
676  }
677  // error: a should be 0, 1 or 2
678  else
679  {
680  throw OomphLibError(
681  "a should equal 0, 1 or 2",
682  OOMPH_CURRENT_FUNCTION,
683  OOMPH_EXCEPTION_LOCATION);
684  }
685 
686  //Jacobian entries
687  if(flag)
688  {
689  //Loop over the displacement basis functions again
690  for(unsigned l2=0;l2<n_node;l2++)
691  {
692  // define complex terms used to obtain entries of current row in the Jacobian:
693 
694  // terms for rows of Jacobian matrix corresponding to r-equation
695  if(a==0)
696  {
697  G_r=(mu*(2.0*dpsidx(l2,0)*dpsidx(l,0)+
698  (nsq+2.0)/rsq*psi(l2)*psi(l)+
699  dpsidx(l2,1)*dpsidx(l,1))+
700  lambda*(dpsidx(l2,0)+psi(l2)/r)*
701  (dpsidx(l,0)+psi(l)/r)-omega_sq*psi(l2)*psi(l))*r*W;
702 
703  G_z=(mu*dpsidx(l2,0)*dpsidx(l,1)+
704  lambda*dpsidx(l2,1)*(dpsidx(l,0)+psi(l)/r))*r*W;
705 
706  G_theta=(-I*(mu*(n/r*dpsidx(l2,0)*psi(l)-3.0*n/rsq*psi(l2)*psi(l))-
707  lambda*n/r*psi(l2)*(dpsidx(l,0)+psi(l)/r)))*r*W;
708  }
709  // terms for rows of Jacobian matrix corresponding to z-equation
710  else if(a==1)
711  {
712  G_r=(mu*dpsidx(l2,1)*dpsidx(l,0)+
713  lambda*(dpsidx(l2,0)+psi(l2)/r)*dpsidx(l,1))*r*W;
714 
715  G_z=(mu*(dpsidx(l2,0)*dpsidx(l,0)+nsq/rsq*psi(l2)*psi(l)
716  +2.0*dpsidx(l2,1)*dpsidx(l,1))+
717  lambda*dpsidx(l2,1)*dpsidx(l,1)-omega_sq*psi(l2)*psi(l))*r*W;
718 
719  G_theta=(-I*(mu*n/r*dpsidx(l2,1)*psi(l)-
720  lambda*n/r*dpsidx(l,1)*psi(l2)))*r*W;
721  }
722  // terms for rows of Jacobian matrix corresponding to theta-equation
723  else if(a==2)
724  {
725  G_r=(-I*(mu*(-n/r*psi(l2)*dpsidx(l,0)+3.0*n/rsq*psi(l2)*psi(l))+
726  lambda*n/r*(dpsidx(l2,0)+psi(l2)/r)*psi(l)))*r*W;
727 
728  G_z=(-I*(-mu*n/r*psi(l2)*dpsidx(l,1)+
729  lambda*n/r*dpsidx(l2,1)*psi(l)))*r*W;
730 
731  G_theta=(mu*((dpsidx(l2,0)-psi(l2)/r)*(dpsidx(l,0)-psi(l)/r)+
732  2.0*nsq/rsq*psi(l2)*psi(l)+dpsidx(l2,1)*dpsidx(l,1))+
733  lambda*nsq/rsq*psi(l2)*psi(l)-omega_sq*psi(l2)*psi(l))*r*W;
734  }
735 
736  //Loop over the displacement components
737  for(unsigned c=0;c<3;c++)
738  {
739  // Get real and imaginary parts of local unknown
740  local_unknown_real=this->nodal_local_eqn(l2,u_nodal_index[c].real());
741  local_unknown_imag=this->nodal_local_eqn(l2,u_nodal_index[c].imag());
742 
743  //If the real part of the local unknown is not pinned
744  if(local_unknown_real >= 0)
745  {
746  if(c==0)
747  {
748  jacobian(local_eqn_imag,local_unknown_real) += G_r.imag();
749  }
750  else if(c==1)
751  {
752  jacobian(local_eqn_imag,local_unknown_real) += G_z.imag();
753  }
754  else if(c==2)
755  {
756  jacobian(local_eqn_imag,local_unknown_real) += G_theta.imag();
757  }
758  }
759 
760  //If the imaginary part of the local unknown is not pinned
761  if(local_unknown_imag >= 0)
762  {
763  if(c==0)
764  {
765  jacobian(local_eqn_imag,local_unknown_imag) += G_r.real();
766  }
767  else if(c==1)
768  {
769  jacobian(local_eqn_imag,local_unknown_imag) += G_z.real();
770  }
771  else if(c==2)
772  {
773  jacobian(local_eqn_imag,local_unknown_imag) +=
774  G_theta.real();
775  }
776  }
777  }
778  }
779  } //End of jacobian calculation
780 
781  } //End of if not boundary condition for imag eqn
782 
783  } //End of loop over coordinate directions
784  } //End of loop over shape functions
785  } //End of loop over integration points
786 
787  }
788 
789 //=======================================================================
790 /// Output exact solution r,z, u_r_real, u_z_real, ..., u_theta_imag
791 //=======================================================================
793  std::ostream &outfile,
794  const unsigned &nplot,
796  {
797  //Vector of local coordinates
798  Vector<double> s(2);
799 
800  // Vector for coordintes
801  Vector<double> x(2);
802 
803  // Tecplot header info
804  outfile << this->tecplot_zone_string(nplot);
805 
806  // Exact solution Vector
807  Vector<double> exact_soln(6);
808 
809  // Loop over plot points
810  unsigned num_plot_points=this->nplot_points(nplot);
811  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
812  {
813 
814  // Get local coordinates of plot point
815  this->get_s_plot(iplot,nplot,s);
816 
817  // Get x position as Vector
818  this->interpolated_x(s,x);
819 
820  // Get exact solution at this point
821  (*exact_soln_pt)(x,exact_soln);
822 
823  //Output r,z,...,u_exact,...
824  for(unsigned i=0;i<2;i++)
825  {
826  outfile << x[i] << " ";
827  }
828  for(unsigned i=0;i<6;i++)
829  {
830  outfile << exact_soln[i] << " ";
831  }
832  outfile << std::endl;
833  }
834 
835  // Write tecplot footer (e.g. FE connectivity lists)
836  this->write_tecplot_zone_footer(outfile,nplot);
837 
838  }
839 
840 //=======================================================================
841 /// Output: r,z, u_r_real, u_z_real, ..., u_theta_imag
842 //=======================================================================
844  output(std::ostream &outfile,
845  const unsigned &nplot)
846  {
847  //Set output Vector
848  Vector<double> s(2);
849  Vector<double> x(2);
851 
852  // Tecplot header info
853  outfile << this->tecplot_zone_string(nplot);
854 
855  // Loop over plot points
856  unsigned num_plot_points=this->nplot_points(nplot);
857  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
858  {
859 
860  // Get local coordinates of plot point
861  this->get_s_plot(iplot,nplot,s);
862 
863  // Get Eulerian coordinates and displacements
864  this->interpolated_x(s,x);
865  this->
867 
868  //Output the r,z,..
869  for(unsigned i=0;i<2;i++)
870  {outfile << x[i] << " ";}
871 
872  // Output real part of displacement
873  for(unsigned i=0;i<3;i++)
874  {outfile << u[i].real() << " ";}
875 
876  // Output imag part of displacement
877  for(unsigned i=0;i<3;i++)
878  {outfile << u[i].imag() << " ";}
879 
880  outfile << std::endl;
881  }
882 
883  // Write tecplot footer (e.g. FE connectivity lists)
884  this->write_tecplot_zone_footer(outfile,nplot);
885  }
886 
887 
888 //=======================================================================
889 /// C-style output:r,z, u_r_real, u_z_real, ..., u_theta_imag
890 //=======================================================================
892  FILE* file_pt,
893  const unsigned &nplot)
894 {
895  //Vector of local coordinates
896  Vector<double> s(2);
897 
898  // Tecplot header info
899  fprintf(file_pt,"%s",this->tecplot_zone_string(nplot).c_str());
900 
901  // Loop over plot points
902  unsigned num_plot_points=this->nplot_points(nplot);
903  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
904  {
905  // Get local coordinates of plot point
906  this->get_s_plot(iplot,nplot,s);
907 
908  // Coordinates
909  for(unsigned i=0;i<2;i++)
910  {
911  fprintf(file_pt,"%g ",this->interpolated_x(s,i));
912  }
913 
914  // Displacement
915  for(unsigned i=0;i<3;i++)
916  {
917  fprintf(file_pt,"%g ",
918  this->
920  }
921  for(unsigned i=0;i<3;i++)
922  {
923  fprintf(file_pt,"%g ",
924  this->
926  }
927  }
928  fprintf(file_pt,"\n");
929 
930  // Write tecplot footer (e.g. FE connectivity lists)
931  this->write_tecplot_zone_footer(file_pt,nplot);
932 
933 }
934 
935 //======================================================================
936 /// Validate against exact solution
937 /// Solution is provided via function pointer.
938 /// Plot at a given number of plot points and compute L2 error
939 /// and L2 norm of velocity solution over element.
940 //=======================================================================
942  std::ostream &outfile,
944  double& error, double& norm)
945 {
946 
947  error=0.0;
948  norm=0.0;
949 
950  //Vector of local coordinates
951  Vector<double> s(2);
952 
953  // Vector for coordinates
954  Vector<double> x(2);
955 
956  //Set the value of n_intpt
957  unsigned n_intpt = this->integral_pt()->nweight();
958 
959  outfile << "ZONE" << std::endl;
960 
961  // Exact solution Vector (u_r_real, u_z_real, ..., u_theta_imag)
962  Vector<double> exact_soln(6);
963 
964  //Loop over the integration points
965  for(unsigned ipt=0;ipt<n_intpt;ipt++)
966  {
967 
968  //Assign values of s
969  for(unsigned i=0;i<2;i++)
970  {
971  s[i] = this->integral_pt()->knot(ipt,i);
972  }
973 
974  //Get the integral weight
975  double w = this->integral_pt()->weight(ipt);
976 
977  // Get jacobian of mapping
978  double J= this->J_eulerian(s);
979 
980  //Premultiply the weights and the Jacobian
981  double W = w*J;
982 
983  // Get x position as Vector
984  this->interpolated_x(s,x);
985 
986  // Get exact solution at this point
987  (*exact_soln_pt)(x,exact_soln);
988 
989  // Displacement error
990  for(unsigned i=0;i<3;i++)
991  {
992  norm+=(exact_soln[i]*exact_soln[i]+exact_soln[i+3]*exact_soln[i+3])*W;
993  error+=(
994  (exact_soln[i]-
995  this->
997  (exact_soln[i]-
998  this->
1000  (exact_soln[i+3]-
1001  this->
1003  (exact_soln[i+3]-
1004  this->
1006  )*W;
1007  }
1008 
1009 
1010  //Output r,z coordinates
1011  for(unsigned i=0;i<2;i++)
1012  {
1013  outfile << x[i] << " ";
1014  }
1015 
1016  //Output ur_error, uz_error, uth_error
1017  for(unsigned i=0;i<3;i++)
1018  {
1019  outfile << exact_soln[i]-
1020  this->
1022  << " ";
1023  }
1024  for(unsigned i=0;i<3;i++)
1025  {
1026  outfile << exact_soln[i+3]-
1027  this->
1029  << " ";
1030  }
1031  outfile << std::endl;
1032  }
1033 }
1034 
1035 // Instantiate required elements
1039 
1040 
1041 }
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution: r,z, u_r_real, u_z_real, ..., u_theta_imag.
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
const std::complex< double > & omega_sq() const
Access function for square of non-dim frequency.
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
virtual std::complex< unsigned > u_index_time_harmonic_fourier_decomposed_linear_elasticity(const unsigned i) const
Return the index at which the i-th (i=0: r, i=1: z; i=2: theta) real or imag unknown displacement com...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
virtual void fill_in_generic_contribution_to_residuals_fourier_decomp_time_harmonic_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Private helper function to compute residuals and (if requested via flag) also the Jacobian matrix...
const std::complex< double > I(0.0, 1.0)
The imaginary unit.
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void output(std::ostream &outfile)
Output: r,z, u_r_real, u_z_real, ..., u_theta_imag.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
static std::complex< double > Default_omega_sq_value
Static default value for squared frequency.
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
void interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(const Vector< double > &s, Vector< std::complex< double > > &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
Class for dense matrices, storing all the values of the matrix as a pointer to a pointer with assorte...
Definition: communicator.h:50
void body_force(const Vector< double > &x, Vector< std::complex< double > > &b) const
Evaluate body force at Eulerian coordinate x at present time (returns zero vector if no body force fu...
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double > > &strain)
Get strain tensor.
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252