pml_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: 1301 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-10-17 08:14:01 +0100 (Tue, 17 Oct 2017) $
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 frequency
40  template <unsigned DIM>
41  double PMLTimeHarmonicLinearElasticityEquationsBase<DIM>::
42  Default_omega_sq_value=1.0;
43 
44 
45 //////////////////////////////////////////////////////////////////
46 //////////////////////////////////////////////////////////////////
47 //////////////////////////////////////////////////////////////////
48 
49 //======================================================================
50 /// Compute the strain tensor at local coordinate s
51 //======================================================================
52  template<unsigned DIM>
54  const Vector<double> &s,
55  DenseMatrix<std::complex<double> >& strain) const
56  {
57 #ifdef PARANOID
58  if ((strain.ncol()!=DIM)||(strain.nrow()!=DIM))
59  {
60  std::ostringstream error_message;
61  error_message << "Strain matrix is " << strain.ncol() << " x "
62  << strain.nrow() << ", but dimension of the equations is "
63  << DIM << std::endl;
64  throw OomphLibError(error_message.str(),
65  OOMPH_CURRENT_FUNCTION,
66  OOMPH_EXCEPTION_LOCATION);
67  }
68 
69  //Find out how many position types there are
70  unsigned n_position_type = this->nnodal_position_type();
71 
72  if(n_position_type != 1)
73  {
74  std::ostringstream error_message;
75  error_message << "PMLTimeHarmonicLinearElasticity is not yet "
76  << "implemented for more than one position type"
77  << std::endl;
78  throw OomphLibError(error_message.str(),
79  OOMPH_CURRENT_FUNCTION,
80  OOMPH_EXCEPTION_LOCATION);
81  }
82 #endif
83 
84 
85  //Find out how many nodes there are in the element
86  unsigned n_node = nnode();
87 
88  //Find the indices at which the local velocities are stored
89  std::complex<unsigned> u_nodal_index[DIM];
90  for(unsigned i=0;i<DIM;i++)
91  {
92  u_nodal_index[i] = u_index_time_harmonic_linear_elasticity(i);
93  }
94 
95  //Set up memory for the shape and derivative functions
96  Shape psi(n_node);
97  DShape dpsidx(n_node,DIM);
98 
99  //Call the derivatives of the shape functions
100  (void) dshape_eulerian(s,psi,dpsidx);
101 
102  //Calculate interpolated values of the derivative of global position
104  interpolated_dudx(DIM,DIM,std::complex<double>(0.0,0.0));
105 
106  //Loop over nodes
107  for(unsigned l=0;l<n_node;l++)
108  {
109  //Loop over velocity components
110  for(unsigned i=0;i<DIM;i++)
111  {
112  //Get the nodal value
113  const std::complex<double> u_value=
114  std::complex<double>(this->nodal_value(l,u_nodal_index[i].real()),
115  this->nodal_value(l,u_nodal_index[i].imag()));
116 
117  //Loop over derivative directions
118  for(unsigned j=0;j<DIM;j++)
119  {
120  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
121  }
122  }
123  }
124 
125  ///Now fill in the entries of the strain tensor
126  for(unsigned i=0;i<DIM;i++)
127  {
128  //Do upper half of matrix
129  //Note that j must be signed here for the comparison test to work
130  //Also i must be cast to an int
131  for(int j=(DIM-1);j>=static_cast<int>(i);j--)
132  {
133  //Off diagonal terms
134  if(static_cast<int>(i)!=j)
135  {
136  strain(i,j) =
137  0.5*(interpolated_dudx(i,j) + interpolated_dudx(j,i));
138  }
139  //Diagonal terms will including growth factor when it comes back in
140  else
141  {
142  strain(i,i) = interpolated_dudx(i,i);
143  }
144  }
145  //Matrix is symmetric so just copy lower half
146  for(int j=(i-1);j>=0;j--)
147  {
148  strain(i,j) = strain(j,i);
149  }
150  }
151  }
152 
153 
154 
155 
156 
157 //======================================================================
158 /// Compute the Cauchy stress tensor at local coordinate s for
159 /// displacement formulation.
160 //======================================================================
161 template<unsigned DIM>
164  DenseMatrix<std::complex<double> >&stress) const
165 {
166 #ifdef PARANOID
167  if ((stress.ncol()!=DIM)||(stress.nrow()!=DIM))
168  {
169  std::ostringstream error_message;
170  error_message << "Stress matrix is " << stress.ncol() << " x "
171  << stress.nrow() << ", but dimension of the equations is "
172  << DIM << std::endl;
173  throw OomphLibError(error_message.str(),
174  OOMPH_CURRENT_FUNCTION,
175  OOMPH_EXCEPTION_LOCATION);
176  }
177 #endif
178 
179  // Get strain
180  DenseMatrix<std::complex<double> > strain(DIM,DIM);
181  this->get_strain(s,strain);
182 
183  // Now fill in the entries of the stress tensor without exploiting
184  // symmetry -- sorry too lazy. This fct is only used for
185  // postprocessing anyway...
186  for(unsigned i=0;i<DIM;i++)
187  {
188  for(unsigned j=0;j<DIM;j++)
189  {
190  stress(i,j) = 0.0;
191  for (unsigned k=0;k<DIM;k++)
192  {
193  for (unsigned l=0;l<DIM;l++)
194  {
195  stress(i,j)+=this->E(i,j,k,l)*strain(k,l);
196  }
197  }
198  }
199  }
200 }
201 
202 
203 //=======================================================================
204 /// Compute the residuals for the linear elasticity equations in
205 /// cartesian coordinates. Flag indicates if we want Jacobian too.
206 //=======================================================================
207 template <unsigned DIM>
210  Vector<double> &residuals, DenseMatrix<double> &jacobian,unsigned flag)
211 {
212  //Find out how many nodes there are
213  unsigned n_node = this->nnode();
214 
215 #ifdef PARANOID
216  //Find out how many positional dofs there are
217  unsigned n_position_type = this->nnodal_position_type();
218 
219  if(n_position_type != 1)
220  {
221  std::ostringstream error_message;
222  error_message << "PMLTimeHarmonicLinearElasticity is not yet "
223  << "implemented for more than one position type"
224  << std::endl;
225  throw OomphLibError(error_message.str(),
226  OOMPH_CURRENT_FUNCTION,
227  OOMPH_EXCEPTION_LOCATION);
228  }
229 
230  //Throw and error if an elasticity tensor has not been set
231  if(this->Elasticity_tensor_pt==0)
232  {
233  throw OomphLibError(
234  "No elasticity tensor set.",
235  OOMPH_CURRENT_FUNCTION,
236  OOMPH_EXCEPTION_LOCATION);
237  }
238 #endif
239 
240  //Find the indices at which the local velocities are stored
241  std::complex<unsigned> u_nodal_index[DIM];
242  for(unsigned i=0;i<DIM;i++)
243  {
244  u_nodal_index[i] = this->u_index_time_harmonic_linear_elasticity(i);
245  }
246 
247 
248  // Square of non-dimensional frequency
249  const double omega_sq_local = this->omega_sq();
250 
251  //Set up memory for the shape functions
252  Shape psi(n_node);
253  DShape dpsidx(n_node,DIM);
254 
255  //Set the value of Nintpt -- the number of integration points
256  unsigned n_intpt = this->integral_pt()->nweight();
257 
258  //Set the vector to hold the local coordinates in the element
259  Vector<double> s(DIM);
260 
261  //Integers to store the local equation numbers
262  int local_eqn_real=0, local_unknown_real=0;
263  int local_eqn_imag=0, local_unknown_imag=0;
264 
265  //Loop over the integration points
266  for(unsigned ipt=0;ipt<n_intpt;ipt++)
267  {
268  //Assign the values of s
269  for(unsigned i=0;i<DIM;++i)
270  {
271  s[i] = this->integral_pt()->knot(ipt,i);
272  }
273 
274  //Get the integral weight
275  double w = this->integral_pt()->weight(ipt);
276 
277  //Call the derivatives of the shape functions (and get Jacobian)
278  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
279 
280  //Storage for Eulerian coordinates (initialised to zero)
281  Vector<double> interpolated_x(DIM,0.0);
282 
283  // Displacement
285  interpolated_u(DIM,std::complex<double>(0.0,0.0));
286 
287  //Calculate interpolated values of the derivative of global position
288  //wrt lagrangian coordinates
290  interpolated_dudx(DIM,DIM,std::complex<double>(0.0,0.0));
291 
292  //Calculate displacements and derivatives
293  for(unsigned l=0;l<n_node;l++)
294  {
295  //Loop over displacement components (deformed position)
296  for(unsigned i=0;i<DIM;i++)
297  {
298  //Calculate the Lagrangian coordinates and the accelerations
299  interpolated_x[i] += this->raw_nodal_position(l,i)*psi(l);
300 
301 
302  //Get the nodal displacements
303  const std::complex<double> u_value
304  = std::complex<double>(
305  this->raw_nodal_value(l,u_nodal_index[i].real()),
306  this->raw_nodal_value(l,u_nodal_index[i].imag()));
307 
308  interpolated_u[i]+=u_value*psi(l);
309 
310  //Loop over derivative directions
311  for(unsigned j=0;j<DIM;j++)
312  {
313  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
314  }
315  }
316  }
317 
318  //Get body force
319  Vector<std::complex<double> > body_force_vec(DIM);
320  this->body_force(interpolated_x,body_force_vec);
321 
322  //Premultiply the weights and the Jacobian
323  double W = w*J;
324 
325 
326  /// \short All the PML weights that participate in the assemby process
327  /// are computed here. pml_inverse_jacobian_diagonals are are used to
328  /// transform derivatives in real x to derivatives in transformed space
329  /// \f$\tilde x \f$.
330  /// pml_jacobian_det allows us to transform volume integrals in
331  /// transformed space to real space.
332  /// If the PML is not enabled via enable_pml, both default to 1.0.
333  Vector< std::complex<double> > pml_inverse_jacobian_diagonal(DIM);
334  std::complex<double> pml_jacobian_det;
335  this->compute_pml_coefficients(ipt, interpolated_x,
336  pml_inverse_jacobian_diagonal,
337  pml_jacobian_det);
338 
339  //Loop over the test functions, nodes of the element
340  for(unsigned l=0;l<n_node;l++)
341  {
342  //Loop over the displacement components
343  for(unsigned a=0;a<DIM;a++)
344  {
345  //Get real and imaginary equation numbers
346  local_eqn_real = this->nodal_local_eqn(l,u_nodal_index[a].real());
347  local_eqn_imag = this->nodal_local_eqn(l,u_nodal_index[a].imag());
348 
349  //===== EQUATIONS OF PML TIME HARMONIC LINEAR ELASTICITY ========
350  // (This is where the maths happens)
351 
352  // Calculate jacobian and residual contributions from acceleration and
353  // body force
354  std::complex<double> mass_jacobian_contribution =
355  -omega_sq_local*pml_jacobian_det;
356  std::complex<double> mass_residual_contribution =
357  (-omega_sq_local*interpolated_u[a]-body_force_vec[a])*pml_jacobian_det;
358 
359  // Calculate jacobian and residual contributions from stress term
360  std::complex<double> stress_jacobian_contributions[DIM][DIM][DIM];
361  std::complex<double> stress_residual_contributions[DIM];
362  for(unsigned b=0;b<DIM;b++)
363  {
364  stress_residual_contributions[b] = 0.0;
365  for(unsigned c=0;c<DIM;c++)
366  {
367  for(unsigned d=0;d<DIM;d++)
368  {
369  stress_jacobian_contributions[b][c][d] =
370  this->E(a,b,c,d)*pml_jacobian_det
371  *pml_inverse_jacobian_diagonal[b]
372  *pml_inverse_jacobian_diagonal[d];
373 
374  stress_residual_contributions[b] +=
375  stress_jacobian_contributions[b][c][d]*interpolated_dudx(c,d);
376  }
377  }
378  }
379 
380  //===== ADD CONTRIBUTIONS TO GLOBAL RESIDUALS AND JACOBIAN ========
381 
382  /*IF it's not a boundary condition*/
383  if(local_eqn_real >= 0)
384  {
385  // Acceleration and body force
386  residuals[local_eqn_real] += mass_residual_contribution.real()
387  *psi(l)*W;
388 
389  // Stress term
390  for(unsigned b=0;b<DIM;b++)
391  {
392  //Add the stress terms to the residuals
393  residuals[local_eqn_real] += stress_residual_contributions[b].real()
394  *dpsidx(l,b)*W;
395  }
396 
397  //Jacobian entries
398  if(flag)
399  {
400  //Loop over the displacement basis functions again
401  for(unsigned l2=0;l2<n_node;l2++)
402  {
403  //Loop over the displacement components again
404  for(unsigned c=0;c<DIM;c++)
405  {
406  local_unknown_real=
407  this->nodal_local_eqn(l2,u_nodal_index[c].real());
408  local_unknown_imag=
409  this->nodal_local_eqn(l2,u_nodal_index[c].imag());
410  //If it's not pinned
411  if(local_unknown_real >= 0)
412  {
413  // Inertial term
414  if (a==c)
415  {
416  jacobian(local_eqn_real,local_unknown_real) +=
417  mass_jacobian_contribution.real()*psi(l)*psi(l2)*W;
418  }
419 
420  // Stress term
421  for(unsigned b=0;b<DIM;b++)
422  {
423  for(unsigned d=0;d<DIM;d++)
424  {
425  //Add the contribution to the Jacobian matrix
426  jacobian(local_eqn_real,local_unknown_real) +=
427  stress_jacobian_contributions[b][c][d].real()
428  *dpsidx(l2,d)*dpsidx(l,b)*W;
429  }
430  }
431  } //End of if not boundary condition
432 
433  if(local_unknown_imag >= 0)
434  {
435  // Inertial term
436  if (a==c)
437  {
438  jacobian(local_eqn_real,local_unknown_imag) -=
439  mass_jacobian_contribution.imag()*psi(l)*psi(l2)*W;
440  }
441 
442  // Stress term
443  for(unsigned b=0;b<DIM;b++)
444  {
445  for(unsigned d=0;d<DIM;d++)
446  {
447  //Add the contribution to the Jacobian matrix
448  jacobian(local_eqn_real,local_unknown_imag) -=
449  stress_jacobian_contributions[b][c][d].imag()
450  *dpsidx(l2,d)*dpsidx(l,b)*W;
451  }
452  }
453  } //End of if not boundary condition
454 
455  }
456  }
457  } //End of jacobian calculation
458 
459  } //End of if not boundary condition for real eqn
460 
461 
462  /*IF it's not a boundary condition*/
463  if(local_eqn_imag >= 0)
464  {
465  // Acceleration and body force
466  residuals[local_eqn_imag] += mass_residual_contribution.imag()
467  *psi(l)*W;
468 
469  // Stress term
470  for(unsigned b=0;b<DIM;b++)
471  {
472  //Add the stress terms to the residuals
473  residuals[local_eqn_imag] += stress_residual_contributions[b].imag()
474  *dpsidx(l,b)*W;
475  }
476 
477  //Jacobian entries
478  if(flag)
479  {
480  //Loop over the displacement basis functions again
481  for(unsigned l2=0;l2<n_node;l2++)
482  {
483  //Loop over the displacement components again
484  for(unsigned c=0;c<DIM;c++)
485  {
486  local_unknown_imag=
487  this->nodal_local_eqn(l2,u_nodal_index[c].imag());
488  local_unknown_real=
489  this->nodal_local_eqn(l2,u_nodal_index[c].real());
490  //If it's not pinned
491  if(local_unknown_imag >= 0)
492  {
493  // Inertial term
494  if (a==c)
495  {
496  jacobian(local_eqn_imag,local_unknown_imag) +=
497  mass_jacobian_contribution.real()*psi(l)*psi(l2)*W;
498  }
499 
500  // Stress term
501  for(unsigned b=0;b<DIM;b++)
502  {
503  for(unsigned d=0;d<DIM;d++)
504  {
505  //Add the contribution to the Jacobian matrix
506  jacobian(local_eqn_imag,local_unknown_imag) +=
507  stress_jacobian_contributions[b][c][d].real()
508  *dpsidx(l2,d)*dpsidx(l,b)*W;
509  }
510  }
511  } //End of if not boundary condition
512 
513  if(local_unknown_real >= 0)
514  {
515  // Inertial term
516  if (a==c)
517  {
518  jacobian(local_eqn_imag,local_unknown_real) +=
519  mass_jacobian_contribution.imag()*psi(l)*psi(l2)*W;
520  }
521 
522  // Stress term
523  for(unsigned b=0;b<DIM;b++)
524  {
525  for(unsigned d=0;d<DIM;d++)
526  {
527  //Add the contribution to the Jacobian matrix
528  jacobian(local_eqn_imag,local_unknown_real) +=
529  stress_jacobian_contributions[b][c][d].imag()
530  *dpsidx(l2,d)*dpsidx(l,b)*W;
531  }
532  }
533  } //End of if not boundary condition
534 
535  }
536  }
537  } //End of jacobian calculation
538 
539  } //End of if not boundary condition for imag eqn
540 
541  } //End of loop over coordinate directions
542  } //End of loop over shape functions
543  } //End of loop over integration points
544 
545  }
546 
547 //=======================================================================
548 /// Output exact solution x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
549 //=======================================================================
550  template <unsigned DIM>
552  std::ostream &outfile,
553  const unsigned &nplot,
555  {
556  //Vector of local coordinates
557  Vector<double> s(DIM);
558 
559  // Vector for coordintes
560  Vector<double> x(DIM);
561 
562  // Tecplot header info
563  outfile << this->tecplot_zone_string(nplot);
564 
565  // Exact solution Vector
566  Vector<double> exact_soln(2*DIM);
567 
568  // Loop over plot points
569  unsigned num_plot_points=this->nplot_points(nplot);
570  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
571  {
572 
573  // Get local coordinates of plot point
574  this->get_s_plot(iplot,nplot,s);
575 
576  // Get x position as Vector
577  this->interpolated_x(s,x);
578 
579  // Get exact solution at this point
580  (*exact_soln_pt)(x,exact_soln);
581 
582  //Output x,y,...,u_exact,...
583  for(unsigned i=0;i<DIM;i++)
584  {
585  outfile << x[i] << " ";
586  }
587  for(unsigned i=0;i<2*DIM;i++)
588  {
589  outfile << exact_soln[i] << " ";
590  }
591  outfile << std::endl;
592  }
593 
594  // Write tecplot footer (e.g. FE connectivity lists)
595  this->write_tecplot_zone_footer(outfile,nplot);
596 
597  }
598 
599 
600 
601 //=======================================================================
602 /// Output: x,y,[z],u,v,[w]
603 //=======================================================================
604  template <unsigned DIM>
606  std::ostream &outfile, const unsigned &nplot)
607  {
608  // Initialise local coord, global coord and solution vectors
609  Vector<double> s(DIM);
610  Vector<double> x(DIM);
612 
613  // Tecplot header info
614  outfile << this->tecplot_zone_string(nplot);
615 
616  // Loop over plot points
617  unsigned num_plot_points=this->nplot_points(nplot);
618  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
619  {
620 
621  // Get local coordinates of plot point
622  this->get_s_plot(iplot,nplot,s);
623 
624  // Get Eulerian coordinates and displacements
625  this->interpolated_x(s,x);
626  this->interpolated_u_time_harmonic_linear_elasticity(s,u);
627 
628  //Output the x,y,..
629  for(unsigned i=0;i<DIM;i++)
630  {outfile << x[i] << " ";}
631 
632  // Output u,v,..
633  for(unsigned i=0;i<DIM;i++)
634  {outfile << u[i].real() << " ";}
635 
636  // Output u,v,..
637  for(unsigned i=0;i<DIM;i++)
638  {outfile << u[i].imag() << " ";}
639 
640  outfile << std::endl;
641  }
642 
643  // Write tecplot footer (e.g. FE connectivity lists)
644  this->write_tecplot_zone_footer(outfile,nplot);
645  }
646 
647 
648 //======================================================================
649 /// Output function for real part of full time-dependent solution
650 /// constructed by adding the scattered field
651 ///
652 /// u = Re( (u_r +i u_i) exp(-i omega t)
653 ///
654 /// at phase angle omega t = phi computed here, to the corresponding
655 /// incoming wave specified via the function pointer.
656 ///
657 /// x,y,u or x,y,z,u
658 ///
659 /// Output at nplot points in each coordinate direction
660 //======================================================================
661 template <unsigned DIM>
663  std::ostream &outfile,
664  FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt,
665  const double& phi,
666  const unsigned &nplot)
667 {
668 
669  // Initialise local coord, global coord and solution vectors
670  Vector<double> s(DIM);
671  Vector<double> x(DIM);
673 
674  // Real and imag part of incoming wave
675  Vector<double> incoming_soln(2*DIM);
676 
677  // Tecplot header info
678  outfile << this->tecplot_zone_string(nplot);
679 
680  // Loop over plot points
681  unsigned num_plot_points=this->nplot_points(nplot);
682  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
683  {
684  // Get local coordinates of plot point
685  this->get_s_plot(iplot,nplot,s);
686 
687  // Get Eulerian coordinates and displacements
688  this->interpolated_x(s,x);
689  this->interpolated_u_time_harmonic_linear_elasticity(s,u);
690 
691  // Get exact solution at this point
692  (*incoming_wave_fct_pt)(x,incoming_soln);
693 
694  //Output the x,y,..
695  for(unsigned i=0;i<DIM;i++)
696  {outfile << x[i] << " ";}
697 
698  // Output u,v,..
699  for(unsigned i=0;i<DIM;i++)
700  {
701  outfile << (u[i].real()+incoming_soln[2*i])*cos(phi)+
702  (u[i].imag()+incoming_soln[2*i+1])*sin(phi) << " ";
703  }
704 
705  outfile << std::endl;
706 
707  }
708 
709  // Write tecplot footer (e.g. FE connectivity lists)
710  this->write_tecplot_zone_footer(outfile,nplot);
711 
712 }
713 
714 //======================================================================
715 /// Output function for imaginary part of full time-dependent solution
716 ///
717 /// u = Im( (u_r +i u_i) exp(-i omega t))
718 ///
719 /// at phase angle omega t = phi.
720 ///
721 /// x,y,u or x,y,z,u
722 ///
723 /// Output at nplot points in each coordinate direction
724 //======================================================================
725 template <unsigned DIM>
727  ::output_real(std::ostream &outfile,
728  const double& phi,
729  const unsigned &nplot)
730 {
731 
732  // Initialise local coord, global coord and solution vectors
733  Vector<double> s(DIM);
734  Vector<double> x(DIM);
736 
737  // Tecplot header info
738  outfile << this->tecplot_zone_string(nplot);
739 
740  // Loop over plot points
741  unsigned num_plot_points=this->nplot_points(nplot);
742  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
743  {
744 
745  // Get local coordinates of plot point
746  this->get_s_plot(iplot,nplot,s);
747 
748  // Get Eulerian coordinates and displacements
749  this->interpolated_x(s,x);
750  this->interpolated_u_time_harmonic_linear_elasticity(s,u);
751 
752  //Output the x,y,..
753  for(unsigned i=0;i<DIM;i++)
754  {outfile << x[i] << " ";}
755 
756  // Output u,v,..
757  for(unsigned i=0;i<DIM;i++)
758  {
759  outfile << u[i].real()*cos(phi)+ u[i].imag()*sin(phi) << " ";
760  }
761 
762  outfile << std::endl;
763 
764  }
765 
766  // Write tecplot footer (e.g. FE connectivity lists)
767  this->write_tecplot_zone_footer(outfile,nplot);
768 
769 }
770 
771 //======================================================================
772 /// Output function for imaginary part of full time-dependent solution
773 ///
774 /// u = Im( (u_r +i u_i) exp(-i omega t))
775 ///
776 /// at phase angle omega t = phi.
777 ///
778 /// x,y,u or x,y,z,u
779 ///
780 /// Output at nplot points in each coordinate direction
781 //======================================================================
782 template <unsigned DIM>
784  ::output_imag(std::ostream &outfile,
785  const double& phi,
786  const unsigned &nplot)
787 {
788 
789  // Initialise local coord, global coord and solution vectors
790  Vector<double> s(DIM);
791  Vector<double> x(DIM);
793 
794  // Tecplot header info
795  outfile << this->tecplot_zone_string(nplot);
796 
797  // Loop over plot points
798  unsigned num_plot_points=this->nplot_points(nplot);
799  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
800  {
801 
802  // Get local coordinates of plot point
803  this->get_s_plot(iplot,nplot,s);
804 
805  // Get Eulerian coordinates and displacements
806  this->interpolated_x(s,x);
807  this->interpolated_u_time_harmonic_linear_elasticity(s,u);
808 
809  //Output the x,y,..
810  for(unsigned i=0;i<DIM;i++)
811  {outfile << x[i] << " ";}
812 
813  // Output u,v,..
814  for(unsigned i=0;i<DIM;i++)
815  {
816  outfile << u[i].imag()*cos(phi)- u[i].real()*sin(phi) << " ";
817  }
818 
819  outfile << std::endl;
820 
821  }
822 
823  // Write tecplot footer (e.g. FE connectivity lists)
824  this->write_tecplot_zone_footer(outfile,nplot);
825 
826 }
827 
828 
829 //=======================================================================
830 /// C-style output: x,y,[z],u,v,[w]
831 //=======================================================================
832 template <unsigned DIM>
834  FILE* file_pt, const unsigned &nplot)
835 {
836  //Vector of local coordinates
837  Vector<double> s(DIM);
838 
839  // Tecplot header info
840  fprintf(file_pt,"%s",this->tecplot_zone_string(nplot).c_str());
841 
842  // Loop over plot points
843  unsigned num_plot_points=this->nplot_points(nplot);
844  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
845  {
846  // Get local coordinates of plot point
847  this->get_s_plot(iplot,nplot,s);
848 
849  // Coordinates
850  for(unsigned i=0;i<DIM;i++)
851  {
852  fprintf(file_pt,"%g ",this->interpolated_x(s,i));
853  }
854 
855  // Displacement
856  for(unsigned i=0;i<DIM;i++)
857  {
858  fprintf(file_pt,"%g ",
859  this->interpolated_u_time_harmonic_linear_elasticity(s,i).real());
860  }
861  for(unsigned i=0;i<DIM;i++)
862  {
863  fprintf(file_pt,"%g ",
864  this->interpolated_u_time_harmonic_linear_elasticity(s,i).imag());
865  }
866  }
867  fprintf(file_pt,"\n");
868 
869  // Write tecplot footer (e.g. FE connectivity lists)
870  this->write_tecplot_zone_footer(file_pt,nplot);
871 
872 }
873 
874 
875 //=======================================================================
876 /// Compute norm of the solution
877 //=======================================================================
878 template <unsigned DIM>
880  double& norm)
881 {
882 
883  // Initialise
884  norm=0.0;
885 
886  //Vector of local coordinates
887  Vector<double> s(DIM);
888 
889  // Vector for coordintes
890  Vector<double> x(DIM);
891 
892  // Displacement vector
893  Vector<std::complex<double> > disp(DIM);
894 
895  //Find out how many nodes there are in the element
896  unsigned n_node = this->nnode();
897 
898  Shape psi(n_node);
899 
900  //Set the value of n_intpt
901  unsigned n_intpt = this->integral_pt()->nweight();
902 
903  //Loop over the integration points
904  for(unsigned ipt=0;ipt<n_intpt;ipt++)
905  {
906 
907  //Assign values of s
908  for(unsigned i=0;i<DIM;i++)
909  {
910  s[i] = this->integral_pt()->knot(ipt,i);
911  }
912 
913  //Get the integral weight
914  double w = this->integral_pt()->weight(ipt);
915 
916  // Get jacobian of mapping
917  double J = this->J_eulerian(s);
918 
919  //Premultiply the weights and the Jacobian
920  double W = w*J;
921 
922  // Get FE function value
923  this->interpolated_u_time_harmonic_linear_elasticity(s,disp);
924 
925  // Add to norm
926  for (unsigned ii=0;ii<DIM;ii++)
927  {
928  norm+=(disp[ii].real()*disp[ii].real()+disp[ii].imag()*disp[ii].imag())*W;
929  }
930  }
931 }
932 
933 //======================================================================
934  /// Validate against exact solution
935  ///
936  /// Solution is provided via function pointer.
937  ///
938 //======================================================================
939 template <unsigned DIM>
941  std::ostream &outfile,
943  double& error, double& norm)
944 {
945 
946  // Initialise
947  error=0.0;
948  norm=0.0;
949 
950  //Vector of local coordinates
951  Vector<double> s(DIM);
952 
953  // Vector for coordintes
954  Vector<double> x(DIM);
955 
956  // Displacement vector
957  Vector<std::complex<double> > disp(DIM);
958 
959  //Find out how many nodes there are in the element
960  unsigned n_node = this->nnode();
961 
962  Shape psi(n_node);
963 
964  //Set the value of n_intpt
965  unsigned n_intpt = this->integral_pt()->nweight();
966 
967  // Exact solution Vector
968  Vector<double> exact_soln(2*DIM);
969 
970  //Loop over the integration points
971  for(unsigned ipt=0;ipt<n_intpt;ipt++)
972  {
973 
974  //Assign values of s
975  for(unsigned i=0;i<DIM;i++)
976  {
977  s[i] = this->integral_pt()->knot(ipt,i);
978  }
979 
980  //Get the integral weight
981  double w = this->integral_pt()->weight(ipt);
982 
983  // Get jacobian of mapping
984  double J = this->J_eulerian(s);
985 
986  //Premultiply the weights and the Jacobian
987  double W = w*J;
988 
989  // Get x position as Vector
990  this->interpolated_x(s,x);
991 
992  // Get exact solution at this point
993  (*exact_soln_pt)(x,exact_soln);
994 
995  // Get FE function value
996  this->interpolated_u_time_harmonic_linear_elasticity(s,disp);
997 
998  // Add to norm
999  for (unsigned ii=0;ii<DIM;ii++)
1000  {
1001  // Add to error and norm
1002  error+=((exact_soln[ii]-disp[ii].real())
1003  *(exact_soln[ii]-disp[ii].real())+
1004  (exact_soln[ii+DIM]-disp[ii].imag())
1005  *(exact_soln[ii+DIM]-disp[ii].imag()))*W;
1006  norm+=(disp[ii].real()*disp[ii].real()+disp[ii].imag()*disp[ii].imag())*W;
1007  }
1008  }
1009 }
1010 
1011 
1012 //Instantiate the required elements
1015 
1019 
1020 template<unsigned DIM> ContinuousBermudezPMLMapping
1022 
1023 
1024 }
virtual void fill_in_generic_contribution_to_residuals_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...
void compute_norm(double &norm)
Compute norm of solution: square of the L2 norm.
cstr elem_len * i
Definition: cfortran.h:607
void output_real(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for real part of full time-dependent solution u = Re( (u_r +i u_i) exp(-i omega t)) a...
void output(std::ostream &outfile)
Output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
void output_total_real(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt, const double &phi, const unsigned &nplot)
static char t char * s
Definition: cfortran.h:572
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Class for dense matrices, storing all the values of the matrix as a pointer to a pointer with assorte...
Definition: communicator.h:50
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double > > &strain) const
Return the strain tensor.
void output_imag(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for imaginary part of full time-dependent solution u = Im( (u_r +i u_i) exp(-i omega ...
void get_stress(const Vector< double > &s, DenseMatrix< std::complex< double > > &sigma) const
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...