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 functions for elements that solve the equations of linear
31 // elasticity in cartesian coordinates
32 
34 
35 
36 namespace oomph
37 {
38 
39 
40 /// Static default value for square of frequency
41  template <unsigned DIM>
42  double TimeHarmonicLinearElasticityEquationsBase<DIM>::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  throw OomphLibError(
75  "TimeHarmonicLinearElasticity is not yet implemented for more than one position type",
76  OOMPH_CURRENT_FUNCTION,
77  OOMPH_EXCEPTION_LOCATION);
78  }
79 #endif
80 
81 
82  //Find out how many nodes there are in the element
83  unsigned n_node = nnode();
84 
85  //Find the indices at which the local velocities are stored
86  std::complex<unsigned> u_nodal_index[DIM];
87  for(unsigned i=0;i<DIM;i++)
88  {
89  u_nodal_index[i] = u_index_time_harmonic_linear_elasticity(i);
90  }
91 
92  //Set up memory for the shape and derivative functions
93  Shape psi(n_node);
94  DShape dpsidx(n_node,DIM);
95 
96  //Call the derivatives of the shape functions
97  (void) dshape_eulerian(s,psi,dpsidx);
98 
99  //Calculate interpolated values of the derivative of global position
101  interpolated_dudx(DIM,DIM,std::complex<double>(0.0,0.0));
102 
103  //Loop over nodes
104  for(unsigned l=0;l<n_node;l++)
105  {
106  //Loop over velocity components
107  for(unsigned i=0;i<DIM;i++)
108  {
109  //Get the nodal value
110  const std::complex<double> u_value=
111  std::complex<double>(this->nodal_value(l,u_nodal_index[i].real()),
112  this->nodal_value(l,u_nodal_index[i].imag()));
113 
114  //Loop over derivative directions
115  for(unsigned j=0;j<DIM;j++)
116  {
117  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
118  }
119  }
120  }
121 
122  ///Now fill in the entries of the strain tensor
123  for(unsigned i=0;i<DIM;i++)
124  {
125  //Do upper half of matrix
126  //Note that j must be signed here for the comparison test to work
127  //Also i must be cast to an int
128  for(int j=(DIM-1);j>=static_cast<int>(i);j--)
129  {
130  //Off diagonal terms
131  if(static_cast<int>(i)!=j)
132  {
133  strain(i,j) =
134  0.5*(interpolated_dudx(i,j) + interpolated_dudx(j,i));
135  }
136  //Diagonal terms will including growth factor when it comes back in
137  else
138  {
139  strain(i,i) = interpolated_dudx(i,i);
140  }
141  }
142  //Matrix is symmetric so just copy lower half
143  for(int j=(i-1);j>=0;j--)
144  {
145  strain(i,j) = strain(j,i);
146  }
147  }
148  }
149 
150 
151 
152 
153 
154 //======================================================================
155 /// Compute the Cauchy stress tensor at local coordinate s for
156 /// displacement formulation.
157 //======================================================================
158 template<unsigned DIM>
161  DenseMatrix<std::complex<double> >&stress) const
162 {
163 #ifdef PARANOID
164  if ((stress.ncol()!=DIM)||(stress.nrow()!=DIM))
165  {
166  std::ostringstream error_message;
167  error_message << "Stress matrix is " << stress.ncol() << " x "
168  << stress.nrow() << ", but dimension of the equations is "
169  << DIM << std::endl;
170  throw OomphLibError(error_message.str(),
171  OOMPH_CURRENT_FUNCTION,
172  OOMPH_EXCEPTION_LOCATION);
173  }
174 #endif
175 
176  // Get strain
177  DenseMatrix<std::complex<double> > strain(DIM,DIM);
178  this->get_strain(s,strain);
179 
180  // Now fill in the entries of the stress tensor without exploiting
181  // symmetry -- sorry too lazy. This fct is only used for
182  // postprocessing anyway...
183  for(unsigned i=0;i<DIM;i++)
184  {
185  for(unsigned j=0;j<DIM;j++)
186  {
187  stress(i,j) = 0.0;
188  for (unsigned k=0;k<DIM;k++)
189  {
190  for (unsigned l=0;l<DIM;l++)
191  {
192  stress(i,j)+=this->E(i,j,k,l)*strain(k,l);
193  }
194  }
195  }
196  }
197 }
198 
199 
200 //=======================================================================
201 /// Compute the residuals for the linear elasticity equations in
202 /// cartesian coordinates. Flag indicates if we want Jacobian too.
203 //=======================================================================
204  template <unsigned DIM>
207  Vector<double> &residuals, DenseMatrix<double> &jacobian,unsigned flag)
208  {
209  //Find out how many nodes there are
210  unsigned n_node = this->nnode();
211 
212 #ifdef PARANOID
213  //Find out how many positional dofs there are
214  unsigned n_position_type = this->nnodal_position_type();
215 
216  if(n_position_type != 1)
217  {
218  throw OomphLibError(
219  "TimeHarmonicLinearElasticity is not yet implemented for more than one position type",
220  OOMPH_CURRENT_FUNCTION,
221  OOMPH_EXCEPTION_LOCATION);
222  }
223 
224  //Throw and error if an elasticity tensor has not been set
225  if(this->Elasticity_tensor_pt==0)
226  {
227  throw OomphLibError(
228  "No elasticity tensor set.",
229  OOMPH_CURRENT_FUNCTION,
230  OOMPH_EXCEPTION_LOCATION);
231  }
232 #endif
233 
234  //Find the indices at which the local velocities are stored
235  std::complex<unsigned> u_nodal_index[DIM];
236  for(unsigned i=0;i<DIM;i++)
237  {
238  u_nodal_index[i] = this->u_index_time_harmonic_linear_elasticity(i);
239  }
240 
241 
242  // Square of non-dimensional frequency
243  const double omega_sq_local = this->omega_sq();
244 
245  //Set up memory for the shape functions
246  Shape psi(n_node);
247  DShape dpsidx(n_node,DIM);
248 
249  //Set the value of Nintpt -- the number of integration points
250  unsigned n_intpt = this->integral_pt()->nweight();
251 
252  //Set the vector to hold the local coordinates in the element
253  Vector<double> s(DIM);
254 
255  //Integer to store the local equation number
256  int local_eqn=0, local_unknown=0;
257 
258  //Loop over the integration points
259  for(unsigned ipt=0;ipt<n_intpt;ipt++)
260  {
261  //Assign the values of s
262  for(unsigned i=0;i<DIM;++i)
263  {
264  s[i] = this->integral_pt()->knot(ipt,i);
265  }
266 
267  //Get the integral weight
268  double w = this->integral_pt()->weight(ipt);
269 
270  //Call the derivatives of the shape functions (and get Jacobian)
271  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
272 
273  //Storage for Eulerian coordinates (initialised to zero)
274  Vector<double> interpolated_x(DIM,0.0);
275 
276  // Displacement
278  interpolated_u(DIM,std::complex<double>(0.0,0.0));
279 
280  //Calculate interpolated values of the derivative of global position
281  //wrt lagrangian coordinates
283  interpolated_dudx(DIM,DIM,std::complex<double>(0.0,0.0));
284 
285  //Calculate displacements and derivatives
286  for(unsigned l=0;l<n_node;l++)
287  {
288  //Loop over displacement components (deformed position)
289  for(unsigned i=0;i<DIM;i++)
290  {
291  //Calculate the Lagrangian coordinates and the accelerations
292  interpolated_x[i] += this->raw_nodal_position(l,i)*psi(l);
293 
294 
295  //Get the nodal displacements
296  const std::complex<double> u_value
297  = std::complex<double>(
298  this->raw_nodal_value(l,u_nodal_index[i].real()),
299  this->raw_nodal_value(l,u_nodal_index[i].imag()));
300 
301  interpolated_u[i]+=u_value*psi(l);
302 
303  //Loop over derivative directions
304  for(unsigned j=0;j<DIM;j++)
305  {
306  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
307  }
308  }
309  }
310 
311  //Get body force at current time
313  this->body_force(interpolated_x,b);
314 
315  //Premultiply the weights and the Jacobian
316  double W = w*J;
317 
318 //=====EQUATIONS OF LINEAR ELASTICITY ========
319 
320  //Loop over the test functions, nodes of the element
321  for(unsigned l=0;l<n_node;l++)
322  {
323  //Loop over the displacement components
324  for(unsigned a=0;a<DIM;a++)
325  {
326  //Get the REAL equation number
327  local_eqn = this->nodal_local_eqn(l,u_nodal_index[a].real());
328 
329  /*IF it's not a boundary condition*/
330  if(local_eqn >= 0)
331  {
332  // Acceleration and body force
333  residuals[local_eqn] +=
334  (-omega_sq_local*interpolated_u[a].real()-b[a].real())*psi(l)*W;
335 
336  // Stress term
337  for(unsigned b=0;b<DIM;b++)
338  {
339  for(unsigned c=0;c<DIM;c++)
340  {
341  for(unsigned d=0;d<DIM;d++)
342  {
343  //Add the stress terms to the residuals
344  residuals[local_eqn] +=
345  this->E(a,b,c,d)*interpolated_dudx(c,d).real()*dpsidx(l,b)*W;
346  }
347  }
348  }
349 
350  //Jacobian entries
351  if(flag)
352  {
353  //Loop over the displacement basis functions again
354  for(unsigned l2=0;l2<n_node;l2++)
355  {
356  //Loop over the displacement components again
357  for(unsigned c=0;c<DIM;c++)
358  {
359  local_unknown=this->nodal_local_eqn(l2,u_nodal_index[c].real());
360  //If it's not pinned
361  if(local_unknown >= 0)
362  {
363  // Inertial term
364  if (a==c)
365  {
366  jacobian(local_eqn,local_unknown) -=
367  omega_sq_local*psi(l)*psi(l2)*W;
368  }
369 
370  // Stress term
371  for(unsigned b=0;b<DIM;b++)
372  {
373  for(unsigned d=0;d<DIM;d++)
374  {
375  //Add the contribution to the Jacobian matrix
376  jacobian(local_eqn,local_unknown) +=
377  this->E(a,b,c,d)*dpsidx(l2,d)*dpsidx(l,b)*W;
378  }
379  }
380  } //End of if not boundary condition
381  }
382  }
383  } //End of jacobian calculation
384 
385  } //End of if not boundary condition for real eqn
386 
387 
388 
389  //Get the IMAG equation number
390  local_eqn = this->nodal_local_eqn(l,u_nodal_index[a].imag());
391 
392  /*IF it's not a boundary condition*/
393  if(local_eqn >= 0)
394  {
395  // Acceleration and body force
396  residuals[local_eqn] +=
397  (-omega_sq_local*interpolated_u[a].imag()-b[a].imag())*psi(l)*W;
398 
399  // Stress term
400  for(unsigned b=0;b<DIM;b++)
401  {
402  for(unsigned c=0;c<DIM;c++)
403  {
404  for(unsigned d=0;d<DIM;d++)
405  {
406  //Add the stress terms to the residuals
407  residuals[local_eqn] +=
408  this->E(a,b,c,d)*interpolated_dudx(c,d).imag()*dpsidx(l,b)*W;
409  }
410  }
411  }
412 
413  //Jacobian entries
414  if(flag)
415  {
416  //Loop over the displacement basis functions again
417  for(unsigned l2=0;l2<n_node;l2++)
418  {
419  //Loop over the displacement components again
420  for(unsigned c=0;c<DIM;c++)
421  {
422  local_unknown=this->nodal_local_eqn(l2,u_nodal_index[c].imag());
423  //If it's not pinned
424  if(local_unknown >= 0)
425  {
426  // Inertial term
427  if (a==c)
428  {
429  jacobian(local_eqn,local_unknown) -=
430  omega_sq_local*psi(l)*psi(l2)*W;
431  }
432 
433  // Stress term
434  for(unsigned b=0;b<DIM;b++)
435  {
436  for(unsigned d=0;d<DIM;d++)
437  {
438  //Add the contribution to the Jacobian matrix
439  jacobian(local_eqn,local_unknown) +=
440  this->E(a,b,c,d)*dpsidx(l2,d)*dpsidx(l,b)*W;
441  }
442  }
443  } //End of if not boundary condition
444  }
445  }
446  } //End of jacobian calculation
447 
448  } //End of if not boundary condition for imag eqn
449 
450  } //End of loop over coordinate directions
451  } //End of loop over shape functions
452  } //End of loop over integration points
453 
454  }
455 
456 //=======================================================================
457 /// Output exact solution x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
458 //=======================================================================
459  template <unsigned DIM>
461  std::ostream &outfile,
462  const unsigned &nplot,
464  {
465  //Vector of local coordinates
466  Vector<double> s(DIM);
467 
468  // Vector for coordintes
469  Vector<double> x(DIM);
470 
471  // Tecplot header info
472  outfile << this->tecplot_zone_string(nplot);
473 
474  // Exact solution Vector
475  Vector<double> exact_soln(2*DIM);
476 
477  // Loop over plot points
478  unsigned num_plot_points=this->nplot_points(nplot);
479  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
480  {
481 
482  // Get local coordinates of plot point
483  this->get_s_plot(iplot,nplot,s);
484 
485  // Get x position as Vector
486  this->interpolated_x(s,x);
487 
488  // Get exact solution at this point
489  (*exact_soln_pt)(x,exact_soln);
490 
491  //Output x,y,...,u_exact,...
492  for(unsigned i=0;i<DIM;i++)
493  {
494  outfile << x[i] << " ";
495  }
496  for(unsigned i=0;i<2*DIM;i++)
497  {
498  outfile << exact_soln[i] << " ";
499  }
500  outfile << std::endl;
501  }
502 
503  // Write tecplot footer (e.g. FE connectivity lists)
504  this->write_tecplot_zone_footer(outfile,nplot);
505 
506  }
507 
508 
509 
510 //=======================================================================
511 /// Output: x,y,[z],u,v,[w]
512 //=======================================================================
513  template <unsigned DIM>
515  const unsigned &nplot)
516  {
517  //Set output Vector
518  Vector<double> s(DIM);
519  Vector<double> x(DIM);
521 
522  // Tecplot header info
523  outfile << this->tecplot_zone_string(nplot);
524 
525  // Loop over plot points
526  unsigned num_plot_points=this->nplot_points(nplot);
527  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
528  {
529 
530  // Get local coordinates of plot point
531  this->get_s_plot(iplot,nplot,s);
532 
533  // Get Eulerian coordinates and displacements
534  this->interpolated_x(s,x);
535  this->interpolated_u_time_harmonic_linear_elasticity(s,u);
536 
537  //Output the x,y,..
538  for(unsigned i=0;i<DIM;i++)
539  {outfile << x[i] << " ";}
540 
541  // Output u,v,..
542  for(unsigned i=0;i<DIM;i++)
543  {outfile << u[i].real() << " ";}
544 
545  // Output u,v,..
546  for(unsigned i=0;i<DIM;i++)
547  {outfile << u[i].imag() << " ";}
548 
549  outfile << std::endl;
550  }
551 
552  // Write tecplot footer (e.g. FE connectivity lists)
553  this->write_tecplot_zone_footer(outfile,nplot);
554  }
555 
556 
557 
558 
559 //=======================================================================
560 /// C-style output: x,y,[z],u,v,[w]
561 //=======================================================================
562 template <unsigned DIM>
564  const unsigned &nplot)
565 {
566  //Vector of local coordinates
567  Vector<double> s(DIM);
568 
569  // Tecplot header info
570  fprintf(file_pt,"%s",this->tecplot_zone_string(nplot).c_str());
571 
572  // Loop over plot points
573  unsigned num_plot_points=this->nplot_points(nplot);
574  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
575  {
576  // Get local coordinates of plot point
577  this->get_s_plot(iplot,nplot,s);
578 
579  // Coordinates
580  for(unsigned i=0;i<DIM;i++)
581  {
582  fprintf(file_pt,"%g ",this->interpolated_x(s,i));
583  }
584 
585  // Displacement
586  for(unsigned i=0;i<DIM;i++)
587  {
588  fprintf(file_pt,"%g ",
589  this->interpolated_u_time_harmonic_linear_elasticity(s,i).real());
590  }
591  for(unsigned i=0;i<DIM;i++)
592  {
593  fprintf(file_pt,"%g ",
594  this->interpolated_u_time_harmonic_linear_elasticity(s,i).imag());
595  }
596  }
597  fprintf(file_pt,"\n");
598 
599  // Write tecplot footer (e.g. FE connectivity lists)
600  this->write_tecplot_zone_footer(file_pt,nplot);
601 
602 }
603 
604 
605 //=======================================================================
606 /// Compute norm of the solution
607 //=======================================================================
608 template <unsigned DIM>
610 {
611 
612  // Initialise
613  norm=0.0;
614 
615  //Vector of local coordinates
616  Vector<double> s(DIM);
617 
618  // Vector for coordintes
619  Vector<double> x(DIM);
620 
621  // Displacement vector
622  Vector<std::complex<double> > disp(DIM);
623 
624  //Find out how many nodes there are in the element
625  unsigned n_node = this->nnode();
626 
627  Shape psi(n_node);
628 
629  //Set the value of n_intpt
630  unsigned n_intpt = this->integral_pt()->nweight();
631 
632  //Loop over the integration points
633  for(unsigned ipt=0;ipt<n_intpt;ipt++)
634  {
635 
636  //Assign values of s
637  for(unsigned i=0;i<DIM;i++)
638  {
639  s[i] = this->integral_pt()->knot(ipt,i);
640  }
641 
642  //Get the integral weight
643  double w = this->integral_pt()->weight(ipt);
644 
645  // Get jacobian of mapping
646  double J=this->J_eulerian(s);
647 
648  //Premultiply the weights and the Jacobian
649  double W = w*J;
650 
651  // Get FE function value
652  this->interpolated_u_time_harmonic_linear_elasticity(s,disp);
653 
654  // Add to norm
655  for (unsigned ii=0;ii<DIM;ii++)
656  {
657  norm+=(disp[ii].real()*disp[ii].real()+disp[ii].imag()*disp[ii].imag())*W;
658  }
659  }
660 }
661 
662 
663 //Instantiate the required elements
666 
670 
671 
672 
673 }
cstr elem_len * i
Definition: cfortran.h:607
void output(std::ostream &outfile)
Output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
void compute_norm(double &norm)
Compute norm of solution: square of the L2 norm.
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...
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double > > &strain) const
Return the strain tensor.
static char t char * s
Definition: cfortran.h:572
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 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].
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...