refineable_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 QLinearElasticityElements
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  "LinearElasticity 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 
72  //Find out how many nodes there are
73  unsigned n_node = this->nnode();
74 
75  //Find the indices at which the local displacements are stored
76  unsigned u_nodal_index[DIM];
77  for(unsigned i=0;i<DIM;i++)
78  {u_nodal_index[i] = this->u_index_linear_elasticity(i);}
79 
80  // Timescale ratio (non-dim density)
81 // hierher double Lambda_sq = this->lambda_sq();
82 
83  //Set up memory for the shape functions
84  Shape psi(n_node);
85  DShape dpsidx(n_node,DIM);
86 
87  //Set the value of Nintpt -- the number of integration points
88  unsigned n_intpt = this->integral_pt()->nweight();
89 
90  //Set the vector to hold the local coordinates in the element
91  Vector<double> s(DIM);
92 
93  //Integer to store the local equation number
94  int local_eqn=0, local_unknown=0;
95 
96  // Local storage for pointers to hang_info objects
97  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
98 
99  //Loop over the integration points
100  for(unsigned ipt=0;ipt<n_intpt;ipt++)
101  {
102  //Assign the values of s
103  for(unsigned i=0;i<DIM;++i) {s[i] = this->integral_pt()->knot(ipt,i);}
104 
105  //Get the integral weight
106  double w = this->integral_pt()->weight(ipt);
107 
108  //Call the derivatives of the shape functions (and get Jacobian)
109  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
110 
111  //Storage for Eulerian coordinates (initialised to zero)
112  Vector<double> interpolated_x(DIM,0.0);
113 
114  //Calculate interpolated values of the derivative of displacements
115  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
116 
117  // Setup memory for accelerations (initialised to zero)
118  // hierher Vector<double> accel(DIM,0.0);
119 
120 
121  //Calculate displacements and derivatives
122  for(unsigned l=0;l<n_node;l++)
123  {
124  //Loop over displacement components
125  for(unsigned i=0;i<DIM;i++)
126  {
127  //Calculate the coordinates and the accelerations
128  interpolated_x[i] += this->nodal_position(l,i)*psi(l);
129 
130  // Only compute accelerations if inertia is switched on
131  // otherwise the timestepper might not be able to
132  // work out dx_gen_dt(2,...)
133  //if (this->Unsteady)
134  // {
135  // // hierher accel[i] += this->dnodal_position_dt(2,l,i)*psi(l);
136  // }
137 
138  //Get the nodal displacements
139  const double u_value = nodal_value(l,u_nodal_index[i]);
140 
141  //Loop over derivative directions
142  for(unsigned j=0;j<DIM;j++)
143  {
144  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
145  }
146  }
147  }
148 
149  //Get body force at current time
150  Vector<double> b(DIM);
151  this->body_force(interpolated_x,b);
152 
153  //Premultiply the weights and the Jacobian
154  double W = w*J;
155 
156  //Number of master nodes and storage for the weight of the shape function
157  unsigned n_master=1; double hang_weight=1.0;
158 
159  //Loop over the test functions, nodes of the element
160  for(unsigned l=0;l<n_node;l++)
161  {
162  //Local boolean to indicate whether the node is hanging
163  bool is_node_hanging = node_pt(l)->is_hanging();
164 
165  //If the node is hanging
166  if(is_node_hanging)
167  {
168  hang_info_pt = node_pt(l)->hanging_pt();
169  //Read out number of master nodes from hanging data
170  n_master = hang_info_pt->nmaster();
171  }
172  //Otherwise the node is its own master
173  else
174  {
175  n_master = 1;
176  }
177 
178  //Loop over the master nodes
179  for(unsigned m=0;m<n_master;m++)
180  {
181  //Loop over the displacement components
182  for(unsigned a=0;a<DIM;a++)
183  {
184  //Get the equation number
185  if(is_node_hanging)
186  {
187  //Get the equation number from the master node
188  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
189  u_nodal_index[a]);
190  //Get the hang weight from the master node
191  hang_weight = hang_info_pt->master_weight(m);
192  }
193  //Otherwise the node is not hanging
194  else
195  {
196  local_eqn = this->nodal_local_eqn(l,u_nodal_index[a]);
197  hang_weight = 1.0;
198  }
199 
200  /*IF it's not a boundary condition*/
201  if(local_eqn >= 0)
202  {
203  // Acceleration and body force
204  residuals[local_eqn] +=
205  (/*hierher Lambda_sq*accel[a]*/-b[a])*psi(l)*W*hang_weight;
206 
207  // Stress term
208  for(unsigned b=0;b<DIM;b++)
209  {
210  for(unsigned c=0;c<DIM;c++)
211  {
212  for(unsigned d=0;d<DIM;d++)
213  {
214  //Add the stress terms to the residuals
215  residuals[local_eqn] +=
216  this->E(a,b,c,d)*interpolated_dudx(c,d)*dpsidx(l,b)*W*
217  hang_weight;
218  }
219  }
220  }
221 
222  //Jacobian entries
223  if(flag)
224  {
225  //Number of master nodes and weights
226  unsigned n_master2=1; double hang_weight2=1.0;
227  //Loop over the displacement basis functions again
228  for(unsigned l2=0;l2<n_node;l2++)
229  {
230  //Local boolean to indicate whether the node is hanging
231  bool is_node2_hanging = node_pt(l2)->is_hanging();
232 
233  //If the node is hanging
234  if(is_node2_hanging)
235  {
236  hang_info2_pt = node_pt(l2)->hanging_pt();
237  //Read out number of master nodes from hanging data
238  n_master2 = hang_info2_pt->nmaster();
239  }
240  //Otherwise the node is its own master
241  else
242  {
243  n_master2 = 1;
244  }
245 
246  //Loop over the master nodes
247  for(unsigned m2=0;m2<n_master2;m2++)
248  {
249  //Loop over the displacement components again
250  for(unsigned c=0;c<DIM;c++)
251  {
252 
253  //Get the number of the unknown
254  //If the node is hanging
255  if(is_node2_hanging)
256  {
257  //Get the equation number from the master node
258  local_unknown =
259  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
260  u_nodal_index[c]);
261  //Get the hang weights
262  hang_weight2 = hang_info2_pt->master_weight(m2);
263  }
264  else
265  {
266  local_unknown =
267  this->nodal_local_eqn(l2,u_nodal_index[c]);
268  hang_weight2 = 1.0;
269  }
270 
271  //If it's not pinned
272  if(local_unknown >= 0)
273  {
274  for(unsigned b=0;b<DIM;b++)
275  {
276  for(unsigned d=0;d<DIM;d++)
277  {
278  //Add the contribution to the Jacobian matrix
279  jacobian(local_eqn,local_unknown) +=
280  this->E(a,b,c,d)*dpsidx(l2,d)*dpsidx(l,b)*W
281  *hang_weight*hang_weight2;
282  }
283  }
284  } //End of if not boundary condition
285  }
286  }
287  }
288  } //End of jacobian calculation
289 
290  } //End of if not boundary condition
291  } //End of loop over coordinate directions
292  }
293  } //End of loop over shape functions
294  } //End of loop over integration points
295  }
296 
297 /// Get error against and norm of exact solution
298 template<unsigned DIM>
300  std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_grad_pt,
301  double& error, double& norm)
302  {
303  // Initialise
304  error=0.0;
305  norm=0.0;
306 
307  //Vector of local coordinates
308  Vector<double> s(DIM);
309 
310  // Vector for coordintes
311  Vector<double> x(DIM);
312 
313  //Set the value of n_intpt
314  unsigned n_intpt = this->integral_pt()->nweight();
315 
316  // Setup output structure: Conversion is fishy but it's only output...
317  unsigned nplot;
318  if (DIM==1)
319  {
320  nplot=n_intpt;
321  }
322  else
323  {
324  nplot=unsigned(pow(n_intpt,1.0/double(DIM)));
325  }
326 
327  // Tecplot header info
328  outfile << this->tecplot_zone_string(nplot);
329 
330  // Exact gradient Vector
331  Vector<double> exact_grad(DIM);
332 
333  //Loop over the integration points
334  for(unsigned ipt=0;ipt<n_intpt;ipt++)
335  {
336 
337  //Assign values of s
338  for(unsigned i=0;i<DIM;i++)
339  {
340  s[i] = this->integral_pt()->knot(ipt,i);
341  }
342 
343  //Get the integral weight
344  double w = this->integral_pt()->weight(ipt);
345 
346  // Get jacobian of mapping
347  double J=this->J_eulerian(s);
348 
349  //Premultiply the weights and the Jacobian
350  double W = w*J;
351 
352  // Get x position as Vector
353  this->interpolated_x(s,x);
354 
355  // Get FE du/dx
356  Vector<double> dudx_fe(DIM);
357 
358  // Get exact gradient at this point
359  (*exact_grad_pt)(x,exact_grad);
360 
361  //Output x,y,...,error
362  for(unsigned i=0;i<DIM;i++)
363  {
364  outfile << x[i] << " ";
365  }
366  for(unsigned i=0;i<DIM;i++)
367  {
368  outfile << exact_grad[i] << " " << exact_grad[i]-dudx_fe[i] << std::endl;
369  }
370 
371  // Add to error and norm
372  for(unsigned i=0;i<DIM;i++)
373  {
374  //LinearElasticityEquations<DIM>::get_flux(s,i,dudx_fe);
375  norm+=exact_grad[i]*exact_grad[i]*W;
376  error+=(exact_grad[i]-dudx_fe[i])*(exact_grad[i]-dudx_fe[i])*W;
377  }
378 
379  }
380  }
381 
382 template<unsigned DIM>
384  {
385  if(this->tree_pt()->father_pt()!=0)
386  {
387  // Needed to set the source function pointer (if there is a father)
389  }
390  // Now do the PRefineableQElement further_build()
391  //PRefineableQElement<DIM>::further_build();
392  }
393 
394 
395 
396 
397 //====================================================================
398 /// Force building of required templates
399 //====================================================================
402 
404 //template class PRefineableQLinearElasticityElement<3>;
405 
406 }
Class for Refineable LinearElasticity equations.
cstr elem_len * i
Definition: cfortran.h:607
p-refineable version of 2D QLinearElasticityElement elements
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
void compute_energy_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_grad_pt, double &error, double &norm)
Get error against and norm of exact solution.
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
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
Class that contains data for hanging nodes.
Definition: nodes.h:684
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
void fill_in_generic_contribution_to_residuals_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overloaded helper function to take hanging nodes into account.
void further_build()
Further build function, pass the pointers down to the sons.