overloaded_element_body.h
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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 
27  private:
28 
30  double Error;
31 
34  unsigned Nprev_for_extrapolation_of_strain_rate;
35 
38  bool Use_fixed_point_for_strain_rate;
39 
42  bool Use_aitken_extrapolation;
43 
44  public:
45 
48  MyTaylorHoodElement()
49  {
50  Error=0.0;
51  Nprev_for_extrapolation_of_strain_rate=4;
52  Use_fixed_point_for_strain_rate=false;
53  Use_aitken_extrapolation=false;
54 
55  // Make space for fixed point iteration on strain rate
56  unsigned fp_n_val = 3;
57  Fixed_point_iteration_guess_for_strain_rate.resize(fp_n_val);
58  unsigned n_intpt = integral_pt()->nweight();
59  for (unsigned val=0;val<fp_n_val;val++)
60  {
61  Fixed_point_iteration_guess_for_strain_rate[val].resize(n_intpt);
62  for (unsigned ipt=0;ipt<n_intpt;ipt++)
63  {
64  Fixed_point_iteration_guess_for_strain_rate[val][ipt].resize(3,3,0.0);
65  }
66  }
67 
68  Aitken_index=0;
69  }
70 
71 
74  unsigned& nprev_for_extrapolation_of_strain_rate()
75  {
76  return Nprev_for_extrapolation_of_strain_rate;
77  }
78 
79 
82  void enable_fixed_point_iteration_for_strain_rate()
83  {
84  Aitken_index=0;
85  update_latest_fixed_point_iteration_guess_for_strain_rate();
86  Use_fixed_point_for_strain_rate=true;
87  }
88 
90  void disable_fixed_point_iteration_for_strain_rate()
91  {
92  Use_fixed_point_for_strain_rate=false;
93  Use_aitken_extrapolation=false;
94  }
95 
97  void enable_aitken_extrapolation()
98  {
99  Use_aitken_extrapolation=true;
100  }
101 
103  void disable_aitken_extrapolation()
104  {
105  Use_aitken_extrapolation=false;
106  }
107 
110  void latest_fixed_point_iteration_guess_for_strain_rate(
111  const unsigned& ipt, DenseMatrix<double>& strainrate) const
112  {
113  for (unsigned i=0;i<3;i++)
114  {
115  for (unsigned j=0;j<3;j++)
116  {
117  strainrate(i,j)=
118  Fixed_point_iteration_guess_for_strain_rate[Aitken_index-1][ipt](i,j);
119  }
120  }
121  }
122 
123 
126  void update_latest_fixed_point_iteration_guess_for_strain_rate()
127  {
128  Vector<double> s(2);
129  DenseMatrix<double> strain_rate(3,3,0.0);
130  unsigned n_intpt = integral_pt()->nweight();
131  for (unsigned ipt=0;ipt<n_intpt;ipt++)
132  {
133  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
134  this->strain_rate(s,strain_rate);
135  for (unsigned i=0;i<3;i++)
136  {
137  for (unsigned j=0;j<3;j++)
138  {
139  Fixed_point_iteration_guess_for_strain_rate[Aitken_index][ipt](i,j)
140  =strain_rate(i,j);
141  }
142  }
143  }
144 
145  if(Aitken_index == 2)
146  {
147  Aitken_index = 0;
148  for (unsigned ipt=0;ipt<n_intpt;ipt++)
149  {
150  for (unsigned i=0;i<3;i++)
151  {
152  for (unsigned j=0;j<3;j++)
153  {
154  double v0 =
155  Fixed_point_iteration_guess_for_strain_rate[0][ipt](i,j);
156  double v1 =
157  Fixed_point_iteration_guess_for_strain_rate[1][ipt](i,j);
158  double v2 =
159  Fixed_point_iteration_guess_for_strain_rate[2][ipt](i,j);
160 
161  double new_value=v2;
162 
163  if(Use_aitken_extrapolation)
164  {
165  double max_diff=std::max(std::fabs(v1-v0),std::fabs(v2-v1));
166 
167  if(max_diff > 1.0e-16)
168  {
169  new_value=v2-std::pow(v2-v1,2.0)/(v2-2.0*v1+v0);
170  }
171  }
172 
173  Fixed_point_iteration_guess_for_strain_rate[Aitken_index][ipt](i,j)=
174  new_value;
175 
176  }
177  }
178  }
179  }
180 
181  Aitken_index++;
182  }
183 
184 
189  void extrapolated_strain_rate(const unsigned& ipt,
190  DenseMatrix<double>& strainrate) const
191  {
192  if (Use_fixed_point_for_strain_rate)
193  {
194  latest_fixed_point_iteration_guess_for_strain_rate(ipt,strainrate);
195  }
196  else
197  {
198  Vector<double> s(2);
199  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
200  extrapolated_strain_rate(s,strainrate);
201  }
202  }
203 
204 
209  void extrapolated_strain_rate(const Vector<double>& s,
210  DenseMatrix<double>& strainrate) const
211  {
212 
213 #ifdef PARANOID
214  if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
215  {
216  std::ostringstream error_message;
217  error_message << "The strain rate has incorrect dimensions "
218  << strainrate.ncol() << " x "
219  << strainrate.nrow() << " Not 3" << std::endl;
220 
221  throw OomphLibError(error_message.str(),
222  OOMPH_CURRENT_FUNCTION,
223  OOMPH_EXCEPTION_LOCATION);
224  }
225 #endif
226 
227 
228  // Get required previous strain rates
229  Vector<DenseMatrix<double> >
230  previous_strain_rate(Nprev_for_extrapolation_of_strain_rate);
231  for (unsigned it=0;it<Nprev_for_extrapolation_of_strain_rate;it++)
232  {
233  previous_strain_rate[it].resize(3,3);
234  strain_rate(it+1,s,previous_strain_rate[it]);
235  }
236 
237  // Get timestepper from first node
238  TimeStepper* time_stepper_pt=node_pt(0)->time_stepper_pt();
239 
240 // hierher #ifdef PARANOID
241  if (time_stepper_pt->nprev_values()<
242  Nprev_for_extrapolation_of_strain_rate)
243  {
244  oomph_info << "Won't work: " << time_stepper_pt->nprev_values()
245  << " < "
246  << Nprev_for_extrapolation_of_strain_rate << std::endl;
247  abort();
248  }
249 // hierher #endif
250 
251 
252  // Which extrapolation are we doing?
253  switch(Nprev_for_extrapolation_of_strain_rate)
254  {
255 
256  // Zero-th order extrapolation; use one previous value
257  case 1:
258  {
259  strainrate=previous_strain_rate[0];
260  }
261  break;
262 
263 
264  // First order extrapolation -- two history values
265  case 2:
266  {
267  // Current and previous timesteps
268  double dt=time_stepper_pt->time_pt()->dt(0);
269  double dt_minus_1=time_stepper_pt->time_pt()->dt(1);
270 
271  // Extrapolate
272  for (unsigned i=0;i<3;i++)
273  {
274  for (unsigned j=0;j<3;j++)
275  {
276  double u_minus_one=previous_strain_rate[0](i,j);
277  double u_minus_two=previous_strain_rate[1](i,j);
278 
279  // Rate of changed based on previous two solutions
280  double slope=(u_minus_one-u_minus_two)/dt_minus_1;
281 
282  // Extrapolated value from previous computed one to current one
283  strainrate(i,j)=u_minus_one+slope*dt;
284  }
285  }
286  }
287  break;
288 
289 
290  // Four history values
291  case 4:
292  {
293  // Extrapolate
294  for (unsigned i=0;i<3;i++)
295  {
296  for (unsigned j=0;j<3;j++)
297  {
298  double u_minus_one =previous_strain_rate[0](i,j);
299  double u_minus_two =previous_strain_rate[1](i,j);
300  double u_minus_three=previous_strain_rate[2](i,j);
301  double u_minus_four =previous_strain_rate[3](i,j);
302 
303  // Current and previous timesteps
304  double dt=time_stepper_pt->time_pt()->dt(0);
305  double dt_minus_1=time_stepper_pt->time_pt()->dt(1);
306  double dt_minus_2=time_stepper_pt->time_pt()->dt(2);
307  double dt_minus_3=time_stepper_pt->time_pt()->dt(3);
308 
309 
310  double MapleGenVar1=0.0;
311  double MapleGenVar2=0.0;
312  double MapleGenVar3=0.0;
313  double MapleGenVar4=0.0;
314  double MapleGenVar5=0.0;
315  double MapleGenVar6=0.0;
316  double MapleGenVar7=0.0;
317  double t0=0.0;
318 
319  MapleGenVar2 = -1.0;
320  MapleGenVar7 = dt*dt*dt*dt_minus_1*dt_minus_1*dt_minus_2*u_minus_four-dt*
321 dt*dt*dt_minus_1*dt_minus_1*dt_minus_2*u_minus_three-dt*dt*dt*dt_minus_1*
322 dt_minus_1*dt_minus_3*u_minus_three+dt*dt*dt*dt_minus_1*dt_minus_1*dt_minus_3*
323 u_minus_two+dt*dt*dt*dt_minus_1*dt_minus_2*dt_minus_2*u_minus_four-dt*dt*dt*
324 dt_minus_1*dt_minus_2*dt_minus_2*u_minus_three-2.0*dt*dt*dt*dt_minus_1*
325 dt_minus_2*dt_minus_3*u_minus_three+2.0*dt*dt*dt*dt_minus_1*dt_minus_2*
326 dt_minus_3*u_minus_two-dt*dt*dt*dt_minus_1*dt_minus_3*dt_minus_3*u_minus_three+
327 dt*dt*dt*dt_minus_1*dt_minus_3*dt_minus_3*u_minus_two;
328  MapleGenVar6 = -dt*dt*dt*dt_minus_2*dt_minus_2*dt_minus_3*u_minus_one+dt*
329 dt*dt*dt_minus_2*dt_minus_2*dt_minus_3*u_minus_two-dt*dt*dt*dt_minus_2*
330 dt_minus_3*dt_minus_3*u_minus_one+dt*dt*dt*dt_minus_2*dt_minus_3*dt_minus_3*
331 u_minus_two+2.0*dt*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*u_minus_four
332 -2.0*dt*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*u_minus_three-2.0*dt*dt*
333 dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_3*u_minus_three+2.0*dt*dt*dt_minus_1*
334 dt_minus_1*dt_minus_1*dt_minus_3*u_minus_two+3.0*dt*dt*dt_minus_1*dt_minus_1*
335 dt_minus_2*dt_minus_2*u_minus_four-3.0*dt*dt*dt_minus_1*dt_minus_1*dt_minus_2*
336 dt_minus_2*u_minus_three+MapleGenVar7;
337  MapleGenVar7 = -6.0*dt*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*
338 u_minus_three+6.0*dt*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*u_minus_two
339 -3.0*dt*dt*dt_minus_1*dt_minus_1*dt_minus_3*dt_minus_3*u_minus_three+3.0*dt*dt*
340 dt_minus_1*dt_minus_1*dt_minus_3*dt_minus_3*u_minus_two+dt*dt*dt_minus_1*
341 dt_minus_2*dt_minus_2*dt_minus_2*u_minus_four-dt*dt*dt_minus_1*dt_minus_2*
342 dt_minus_2*dt_minus_2*u_minus_three-3.0*dt*dt*dt_minus_1*dt_minus_2*dt_minus_2*
343 dt_minus_3*u_minus_one-3.0*dt*dt*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*
344 u_minus_three+6.0*dt*dt*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*u_minus_two
345 -3.0*dt*dt*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*u_minus_one;
346  MapleGenVar5 = -3.0*dt*dt*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*
347 u_minus_three+6.0*dt*dt*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*u_minus_two
348 -dt*dt*dt_minus_1*dt_minus_3*dt_minus_3*dt_minus_3*u_minus_three+dt*dt*
349 dt_minus_1*dt_minus_3*dt_minus_3*dt_minus_3*u_minus_two-2.0*dt*dt*dt_minus_2*
350 dt_minus_2*dt_minus_2*dt_minus_3*u_minus_one+2.0*dt*dt*dt_minus_2*dt_minus_2*
351 dt_minus_2*dt_minus_3*u_minus_two-3.0*dt*dt*dt_minus_2*dt_minus_2*dt_minus_3*
352 dt_minus_3*u_minus_one+3.0*dt*dt*dt_minus_2*dt_minus_2*dt_minus_3*dt_minus_3*
353 u_minus_two-dt*dt*dt_minus_2*dt_minus_3*dt_minus_3*dt_minus_3*u_minus_one+dt*dt
354 *dt_minus_2*dt_minus_3*dt_minus_3*dt_minus_3*u_minus_two+MapleGenVar6+
355 MapleGenVar7;
356  MapleGenVar7 = dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*
357 u_minus_four-dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*
358 u_minus_three-dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_3*
359 u_minus_three+dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_3*
360 u_minus_two+2.0*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*
361 u_minus_four-2.0*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*
362 u_minus_three-4.0*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*
363 u_minus_three+4.0*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*
364 u_minus_two-2.0*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_3*dt_minus_3*
365 u_minus_three+2.0*dt*dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_3*dt_minus_3*
366 u_minus_two;
367  MapleGenVar6 = dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_2*
368 u_minus_four-dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_2*
369 u_minus_three-3.0*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*
370 u_minus_one-3.0*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*
371 u_minus_three+6.0*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*
372 u_minus_two-3.0*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*
373 u_minus_one-3.0*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*
374 u_minus_three+6.0*dt*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*
375 u_minus_two-dt*dt_minus_1*dt_minus_1*dt_minus_3*dt_minus_3*dt_minus_3*
376 u_minus_three+dt*dt_minus_1*dt_minus_1*dt_minus_3*dt_minus_3*dt_minus_3*
377 u_minus_two+MapleGenVar7;
378  MapleGenVar7 = -4.0*dt*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_2*
379 dt_minus_3*u_minus_one+4.0*dt*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_2*
380 dt_minus_3*u_minus_two-6.0*dt*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*
381 dt_minus_3*u_minus_one+6.0*dt*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*
382 dt_minus_3*u_minus_two-2.0*dt*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*
383 dt_minus_3*u_minus_one+2.0*dt*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*
384 dt_minus_3*u_minus_two-dt*dt_minus_2*dt_minus_2*dt_minus_2*dt_minus_2*
385 dt_minus_3*u_minus_one+dt*dt_minus_2*dt_minus_2*dt_minus_2*dt_minus_2*
386 dt_minus_3*u_minus_two-2.0*dt*dt_minus_2*dt_minus_2*dt_minus_2*dt_minus_3*
387 dt_minus_3*u_minus_one+2.0*dt*dt_minus_2*dt_minus_2*dt_minus_2*dt_minus_3*
388 dt_minus_3*u_minus_two-dt*dt_minus_2*dt_minus_2*dt_minus_3*dt_minus_3*
389 dt_minus_3*u_minus_one;
390  MapleGenVar4 = dt*dt_minus_2*dt_minus_2*dt_minus_3*dt_minus_3*dt_minus_3*
391 u_minus_two-dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_3*
392 u_minus_one-dt_minus_1*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*dt_minus_3*
393 u_minus_one-2.0*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_2*
394 dt_minus_3*u_minus_one-3.0*dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_2*
395 dt_minus_3*dt_minus_3*u_minus_one-dt_minus_1*dt_minus_1*dt_minus_2*dt_minus_3*
396 dt_minus_3*dt_minus_3*u_minus_one-dt_minus_1*dt_minus_2*dt_minus_2*dt_minus_2*
397 dt_minus_2*dt_minus_3*u_minus_one-2.0*dt_minus_1*dt_minus_2*dt_minus_2*
398 dt_minus_2*dt_minus_3*dt_minus_3*u_minus_one-dt_minus_1*dt_minus_2*dt_minus_2*
399 dt_minus_3*dt_minus_3*dt_minus_3*u_minus_one+MapleGenVar5+MapleGenVar6+
400 MapleGenVar7;
401  MapleGenVar5 = 1/dt_minus_1;
402  MapleGenVar3 = MapleGenVar4*MapleGenVar5;
403  MapleGenVar1 = MapleGenVar2*MapleGenVar3;
404  MapleGenVar2 = 1/dt_minus_2/dt_minus_3/(dt_minus_1*dt_minus_1*dt_minus_2+
405 dt_minus_1*dt_minus_1*dt_minus_3+2.0*dt_minus_1*dt_minus_2*dt_minus_2+3.0*
406 dt_minus_1*dt_minus_2*dt_minus_3+dt_minus_1*dt_minus_3*dt_minus_3+dt_minus_2*
407 dt_minus_2*dt_minus_2+2.0*dt_minus_2*dt_minus_2*dt_minus_3+dt_minus_2*
408 dt_minus_3*dt_minus_3);
409  t0 = MapleGenVar1*MapleGenVar2;
410 
411  // Extrapolated value from previous computed ones to current one
412  strainrate(i,j)=t0;
413 
414  }
415  }
416  }
417  break;
418 
419  default:
420  {
421  oomph_info << "Never get here\n";
422  abort();
423  }
424  break;
425 
426  }
427  }
428 
429 
430 
432  void set_error(const double& error){Error=error;}
433 
435  std::string variable_identifier()
436  {
437  std::string txt="VARIABLES=";
438  txt+="\"r\",";
439  txt+="\"z\",";
440  txt+="\"u_r\",";
441  txt+="\"u_z\",";
442  //txt+="\"u_phi\",";
443  txt+="\"p\",";
444  txt+="\"du_r/dt\",";
445  txt+="\"du_z/dt\",";
446  //txt+="\"u_rm\",";
447  //txt+="\"u_zm\",";
448  //txt+="\"r_h1\",";
449  //txt+="\"z_h1\",";
450  //txt+="\"r_h2\",";
451  //txt+="\"z_h2\",";
452  //txt+="\"u_rh1\",";
453  //txt+="\"u_zh1\",";
454  //txt+="\"u_rh2\",";
455  //txt+="\"u_zh2\",";
456  //txt+="\"strain_rate(0,0)\",";
457  //txt+="\"strain_rate(1,1)\",";
458  //txt+="\"strain_rate(2,2)\",";
459  //txt+="\"strain_rate(0,1)\","; // V19
460  //txt+="\"stress(0,0)\","; // 20
461  //txt+="\"stress(1,1)\","; // 21
462  //txt+="\"stress(2,2)\","; // 22
463  //txt+="\"stress(0,1)\","; // 23
464  txt+="\"invariant_of_strain_rate\","; // 24
465  txt+="\"invariant_of_strain_rate_extrapol\","; // 25
466  txt+="\"invariant_of_stress\","; // 26
467  txt+="\"viscosity\","; // 27
468  txt+="\"yield\","; // 28
469  txt+="\"error\","; // 29
470  txt+="\"size\","; // 30
471  //txt+="\"strain_rate_m1(0,0)\","; //31
472  //txt+="\"strain_rate_m1(1,1)\","; //32
473  //txt+="\"strain_rate_m1(2,2)\","; //33
474  //txt+="\"strain_rate_m1(0,1)\","; //34
475  //txt+="\"strain_rate_m2(0,0)\","; //v35
476  //txt+="\"strain_rate_m2(1,1)\","; //v36
477  //txt+="\"strain_rate_m2(2,2)\","; //v37
478  //txt+="\"strain_rate_m2(0,1)\","; //v38
479  //txt+="\"strain_rate_extrapol(0,0)\","; //v39
480  //txt+="\"strain_rate_extrapol(1,1)\","; //v40
481  //txt+="\"strain_rate_extrapol(2,2)\","; //v41
482  //txt+="\"strain_rate_extrapol(0,1)\","; //v42
483 
484  txt+="\n";
485  return txt;
486  }
487 
488 
490  void output(std::ostream &outfile,
491  const unsigned &nplot)
492  {
493 
494  // Assign dimension
495  unsigned el_dim=2;
496 
497  // Vector of local coordinates
498  Vector<double> s(el_dim);
499 
500  // Acceleration
501  Vector<double> dudt(el_dim);
502 
503  // Mesh velocity
504  Vector<double> mesh_veloc(el_dim,0.0);
505 
506  // Tecplot header info
507  outfile << tecplot_zone_string(nplot);
508 
509  // Find out how many nodes there are
510  unsigned n_node = nnode();
511 
512  // Get continuous time from timestepper of first node
513  //double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
514 
515  //Set up memory for the shape functions
516  Shape psif(n_node);
517  DShape dpsifdx(n_node,el_dim);
518 
519  // Loop over plot points
520  unsigned num_plot_points=nplot_points(nplot);
521  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
522  {
523 
524  // Get local coordinates of plot point
525  get_s_plot(iplot,nplot,s);
526 
527  //Call the derivatives of the shape and test functions
528  dshape_eulerian(s,psif,dpsifdx);
529 
530  //Allocate storage
531  Vector<double> mesh_veloc(el_dim,0.0);
532  Vector<double> int_x(el_dim,0.0);
533  Vector<double> dudt(el_dim,0.0);
534  Vector<double> dudt_ALE(el_dim,0.0);
535  DenseMatrix<double> interpolated_dudx(el_dim+1,el_dim,0.0);
536 
537  // //Initialise everything to zero
538  // for(unsigned i=0;i<el_dim;i++)
539  // {
540  // mesh_veloc[i]=0.0;
541  // dudt[i]=0.0;
542  // dudt_ALE[i]=0.0;
543  // for(unsigned j=0;j<el_dim;j++)
544  // {
545  // interpolated_dudx(i,j) = 0.0;
546  // }
547 
548  // int_x[i] = interpolated_x(s,i);
549  // }
550 
551  //Calculate velocities and derivatives
552 
553  //Loop over directions
554  for(unsigned i=0;i<el_dim;i++)
555  {
556  //Get the index at which velocity i is stored
557  unsigned u_nodal_index = u_index_axi_nst(i);
558  // Loop over nodes
559  for(unsigned l=0;l<n_node;l++)
560  {
561  dudt[i]+=du_dt_axi_nst(l,u_nodal_index)*psif[l];
562  mesh_veloc[i]+=dnodal_position_dt(l,i)*psif[l];
563  }
564  }
565 
566  //Loop over directions
567  for(unsigned i=0;i<el_dim+1;i++)
568  {
569  //Get the index at which velocity i is stored
570  unsigned u_nodal_index = u_index_axi_nst(i);
571  // Loop over nodes
572  for(unsigned l=0;l<n_node;l++)
573  {
574  //Loop over derivative directions for velocity gradients
575  for(unsigned j=0;j<el_dim;j++)
576  {
577  interpolated_dudx(i,j) += nodal_value(l,u_nodal_index)*
578  dpsifdx(l,j);
579  }
580  }
581  }
582 
583  // Get dudt in ALE form (incl mesh veloc)
584  for(unsigned i=0;i<el_dim;i++)
585  {
586  dudt_ALE[i]=dudt[i];
587  for (unsigned k=0;k<el_dim;k++)
588  {
589  dudt_ALE[i]-=mesh_veloc[k]*interpolated_dudx(i,k);
590  }
591  }
592 
593 
594  // Actual rate of strain
595  DenseMatrix<double> rate_of_strain(3,3,0.0);
596  this->strain_rate(s,rate_of_strain);
597 
598  // Extrapolated (or recycle actual one if there's no extrapol)
599  DenseMatrix<double> rate_of_strain_extrapol(rate_of_strain);
600  if (Use_extrapolated_strainrate_to_compute_second_invariant)
601  {
602  this->extrapolated_strain_rate(s,rate_of_strain_extrapol);
603  }
604 
605  // Get associated invariants
606  double second_invariant_strain=
607  SecondInvariantHelper::second_invariant(rate_of_strain);
608  double second_invariant_strain_extrapol=
609  SecondInvariantHelper::second_invariant(rate_of_strain_extrapol);
610 
611  // Get viscosity
612  double viscosity=this->Constitutive_eqn_pt
613  ->viscosity(second_invariant_strain_extrapol);
614 
615  // Stress
616  DenseMatrix<double> stress(3,3,0.0);
617 
618  // Calculate the stress
619  for(unsigned i=0;i<3;i++)
620  {
621  for(unsigned j=0;j<3;j++)
622  {
623  stress(i,j)=viscosity*rate_of_strain(i,j);
624  }
625  }
626 
627  double second_invariant_stress=
628  SecondInvariantHelper::second_invariant(stress);
629 
630  // Flag indicating if material has yielded (1.0) or not (0.0)
631  double yield=0.0;
632 
633  // identify whether material is yielded (1) or not (0)
634  if (sqrt(std::fabs(second_invariant_stress))
635  >= Problem_Parameter::Yield_stress)
636  {
637  yield=1.0;
638  }
639 
640  // Coordinates
641  for(unsigned i=0;i<el_dim;i++)
642  {
643  outfile << interpolated_x(s,i) << " ";
644  }
645 
646  // Velocities
647  for(unsigned i=0;i<el_dim;i++)
648  {
649  outfile << interpolated_u_axi_nst(s,i) << " ";
650  }
651 
652  // Pressure
653  outfile << interpolated_p_axi_nst(s) << " ";
654 
655  // Accelerations
656  for(unsigned i=0;i<el_dim;i++)
657  {
658  outfile << dudt_ALE[i] << " ";
659  }
660 
661  // Mesh velocity
662  //for(unsigned i=0;i<el_dim;i++)
663  // {
664  // outfile << mesh_veloc[i] << " ";
665  // }
666 
667  // History values of coordinates
668  //unsigned n_prev=node_pt(0)->position_time_stepper_pt()->ntstorage();
669  //for (unsigned t=1;t<n_prev-2;t++)
670  // {
671  // for(unsigned i=0;i<el_dim;i++)
672  // {
673  // outfile << interpolated_x(t,s,i) << " ";
674  // }
675  // }
676 
677  // History values of velocities
678  //n_prev=node_pt(0)->time_stepper_pt()->ntstorage();
679  //for (unsigned t=1;t<n_prev-2;t++)
680  // {
681  // for(unsigned i=0;i<el_dim;i++)
682  // {
683  // outfile << interpolated_u_axi_nst(t,s,i) << " ";
684  // }
685  // }
686 
687  //outfile << rate_of_strain(0,0) << " ";
688 
689  //outfile << rate_of_strain(1,1) << " ";
690 
691  //outfile << rate_of_strain(2,2) << " ";
692 
693  //outfile << rate_of_strain(0,1) << " ";
694 
695  //outfile << stress(0,0) << " ";
696 
697  //outfile << stress(1,1) << " ";
698 
699  //outfile << stress(2,2) << " ";
700 
701  //outfile << stress(0,1) << " ";
702 
703  outfile << second_invariant_strain << " ";
704  outfile << second_invariant_strain_extrapol << " ";
705  outfile << second_invariant_stress << " ";
706 
707  outfile << viscosity << " ";
708 
709  outfile << yield << " ";
710 
711  outfile << Error << " ";
712 
713  outfile << size() << " ";
714 
715  // Rate of strain at previous timestep
716  //this->strain_rate(1,s,rate_of_strain);
717  //outfile << rate_of_strain(0,0) << " ";
718  //outfile << rate_of_strain(1,1) << " ";
719  //outfile << rate_of_strain(2,2) << " ";
720  //outfile << rate_of_strain(0,1) << " ";
721 
722 
723  // Rate of strain at second previous timestep
724  //this->strain_rate(2,s,rate_of_strain);
725  //outfile << rate_of_strain(0,0) << " ";
726  //outfile << rate_of_strain(1,1) << " ";
727  //outfile << rate_of_strain(2,2) << " ";
728  //outfile << rate_of_strain(0,1) << " ";
729 
730  // Extrapolated rate of strain
731  //outfile << rate_of_strain_extrapol(0,0) << " ";
732  //outfile << rate_of_strain_extrapol(1,1) << " ";
733  //outfile << rate_of_strain_extrapol(2,2) << " ";
734  //outfile << rate_of_strain_extrapol(0,1) << " ";
735 
736  outfile << std::endl;
737 
738  }
739 
740  // Write tecplot footer (e.g. FE connectivity lists)
741  write_tecplot_zone_footer(outfile,nplot);
742  }
743 
745  void get_Z2_flux(const Vector<double>& s, Vector<double>& flux)
746  {
747  const unsigned DIM=3;
748 #ifdef PARANOID
749  unsigned num_entries=DIM+(DIM*(DIM-1))/2;
750  if (flux.size() < num_entries)
751  {
752  std::ostringstream error_message;
753  error_message << "The flux vector has the wrong number of entries, "
754  << flux.size() << ", whereas it should be at least "
755  << num_entries << std::endl;
756  throw OomphLibError(error_message.str(),
757  "RefineableNavierStokesEquations::get_Z2_flux()",
758  OOMPH_EXCEPTION_LOCATION);
759  }
760 #endif
761 
762  // Get strain rate matrix
763  DenseMatrix<double> strainrate(DIM);
764  DenseMatrix<double> stress(DIM,DIM,0.0);
765 
766  this->strain_rate(s,strainrate);
767  double second_invariant_strain=
768  SecondInvariantHelper::second_invariant(strainrate);
769 
770  // get viscosity
771  double viscosity=this->Constitutive_eqn_pt
772  ->viscosity(second_invariant_strain);
773 
774  // Calculate the stress
775  for(unsigned i=0;i<DIM;i++)
776  {
777  for(unsigned j=0;j<DIM;j++)
778  {
779  stress(i,j)=viscosity*strainrate(i,j);
780  }
781  }
782 
783  // get the second invariant of the stress tensor
784  //double second_invariant_stress=
785  // SecondInvariantHelper::second_invariant(stress);
786 
787  // Pack into flux Vector
788  unsigned icount=0;
789 
790  // identify whether material is yielded (1) or not (0)
791  double yield=0.0;
792  //if(sqrt(std::fabs(second_invariant_stress))
793  // >= Problem_Parameter::Yield_stress)
794  // obacht
795  if(std::fabs(second_invariant_strain) >
796  Problem_Parameter::Critical_strain_rate)
797  {
798  yield=1.0;
799  }
800 
801  // Add bias to create lots of refinement near yield surface
802  // by creating a jump in the z2 flux (Note: only makes sense if
803  // refinement happens at every timestep, otherwise really fine
804  // region gets immediately left behind by rapidly moving
805  // yield surface
806  bool add_yield_bias=false;
807  if (!add_yield_bias) yield=0.0;
808 
809  // Start with diagonal terms
810  for(unsigned i=0;i<DIM;i++)
811  {
812  // Force more refinement at yield surface by adding yield value
813  // hierher: This does the right thing but triangle can't handle the refinement
814  // properly. Mesh seems ok though. MH will investigate triangle mesh adaptation
815  // separately.
816  flux[icount]=strainrate(i,i)+yield;
817  icount++;
818  }
819 
820  //Off diagonals row by row
821  for(unsigned i=0;i<DIM;i++)
822  {
823  for(unsigned j=i+1;j<DIM;j++)
824  {
825  // force more refinement at yield surface by adding yield value
826  flux[icount]=strainrate(i,j)+yield;
827  icount++;
828  }
829  }
830  }
831 
835  double square_of_norm_of_strain_invariant(double& norm_squared,
836  double& extrapolated_norm_squared,
837  double& error_norm_squared)
838  {
839  // Initialise
840  norm_squared=0.0;
841  extrapolated_norm_squared=0.0;
842  error_norm_squared=0.0;
843  double area=0.0;
844 
845  //Number of integration points
846  unsigned n_intpt = integral_pt()->nweight();
847 
848  //Set the Vector to hold local coordinates
849  Vector<double> s(2);
850 
851  //Loop over the integration points
852  for(unsigned ipt=0;ipt<n_intpt;ipt++)
853  {
854  //Assign values of s
855  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
856 
857  //Get the integral weight
858  double w = integral_pt()->weight(ipt);
859 
860  // Get Jacobain of mapping between local and Eulerian coords
861  double J = this->J_eulerian(s);
862 
863  //Premultiply the weights and the Jacobian
864  double W = w*J;
865 
866  // Invariant of actual strain rate tensor
867  DenseMatrix<double> strainrate(3,3,0.0);
868  this->strain_rate(s,strainrate);
869  double invariant=SecondInvariantHelper::second_invariant(strainrate);
870 
871  // Invariant of extrapolated strain rate tensor
872  // (or recycle actual one if there's no extrapol)
873  DenseMatrix<double> extrapolated_strainrate(strainrate);
874  if (Use_extrapolated_strainrate_to_compute_second_invariant)
875  {
876  this->extrapolated_strain_rate(s,extrapolated_strainrate);
877  }
878  double extrapolated_invariant=
879  SecondInvariantHelper::second_invariant(extrapolated_strainrate);
880 
881  //Assemble norms
882  norm_squared+=
883  invariant*invariant*W;
884  extrapolated_norm_squared+=
885  extrapolated_invariant*extrapolated_invariant*W;
886  error_norm_squared+=
887  (extrapolated_invariant-invariant)*
888  (extrapolated_invariant-invariant)*W;
889  area+=W;
890  }
891 
892  return area;
893 
894  }
895 
899  double square_of_norm_of_viscosity(double& norm_squared,
900  double& extrapolated_norm_squared,
901  double& error_norm_squared)
902  {
903  // Initialise
904  norm_squared=0.0;
905  extrapolated_norm_squared=0.0;
906  error_norm_squared=0.0;
907  double area=0.0;
908 
909  //Number of integration points
910  unsigned n_intpt = integral_pt()->nweight();
911 
912  //Set the Vector to hold local coordinates
913  Vector<double> s(2);
914 
915  //Loop over the integration points
916  for(unsigned ipt=0;ipt<n_intpt;ipt++)
917  {
918  //Assign values of s
919  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
920 
921  //Get the integral weight
922  double w = integral_pt()->weight(ipt);
923 
924  // Get Jacobain of mapping between local and Eulerian coords
925  double J = this->J_eulerian(s);
926 
927  //Premultiply the weights and the Jacobian
928  double W = w*J;
929 
930  // Invariant of actual strain rate tensor
931  DenseMatrix<double> strainrate(3,3,0.0);
932  this->strain_rate(s,strainrate);
933  double invariant=SecondInvariantHelper::second_invariant(strainrate);
934 
935  double viscosity=this->Constitutive_eqn_pt->viscosity(invariant);
936 
937  // Invariant of extrapolated strain rate tensor
938  // (or recycle actual one if there's no extrapol)
939  DenseMatrix<double> extrapolated_strainrate(strainrate);
940  if (Use_extrapolated_strainrate_to_compute_second_invariant)
941  {
942  this->extrapolated_strain_rate(s,extrapolated_strainrate);
943  }
944  double extrapolated_invariant=
945  SecondInvariantHelper::second_invariant(extrapolated_strainrate);
946 
947  double extrapolated_viscosity=
948  this->Constitutive_eqn_pt->viscosity(extrapolated_invariant);
949 
950  //Assemble norms
951  norm_squared+=
952  viscosity*viscosity*W;
953  extrapolated_norm_squared+=
954  extrapolated_viscosity*extrapolated_viscosity*W;
955  error_norm_squared+=
956  (extrapolated_viscosity-viscosity)*
957  (extrapolated_viscosity-viscosity)*W;
958  area+=W;
959  }
960 
961  return area;
962 
963  }
964 
968  double square_of_norm_of_fixed_point(double& norm_squared,
969  double& latest_guess_norm_squared,
970  double& error_norm_squared)
971  {
972  // Initialise
973  norm_squared=0.0;
974  latest_guess_norm_squared=0.0;
975  error_norm_squared=0.0;
976  double area=0.0;
977 
978  //Number of integration points
979  unsigned n_intpt = integral_pt()->nweight();
980 
981  //Set the Vector to hold local coordinates
982  Vector<double> s(2);
983 
984  //Loop over the integration points
985  for(unsigned ipt=0;ipt<n_intpt;ipt++)
986  {
987  //Assign values of s
988  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
989 
990  //Get the integral weight
991  double w = integral_pt()->weight(ipt);
992 
993  // Get Jacobain of mapping between local and Eulerian coords
994  double J = this->J_eulerian(s);
995 
996  //Premultiply the weights and the Jacobian
997  double W = w*J;
998 
999  // Invariant of actual strain rate tensor
1000  DenseMatrix<double> current_strainrate(3,3,0.0);
1001  this->strain_rate(s,current_strainrate);
1002  double invariant=
1003  SecondInvariantHelper::second_invariant(current_strainrate);
1004 
1005  // Invariant of extrapolated strain rate tensor
1006  // (or recycle actual one if there's no extrapol)
1007  DenseMatrix<double> latest_guess_strainrate(3,3,0.0);
1008  latest_fixed_point_iteration_guess_for_strain_rate(ipt,
1009  latest_guess_strainrate);
1010 
1011  double latest_guess_invariant=
1012  SecondInvariantHelper::second_invariant(latest_guess_strainrate);
1013 
1014  //Assemble norms
1015  norm_squared+=
1016  invariant*invariant*W;
1017  latest_guess_norm_squared+=
1018  latest_guess_invariant*latest_guess_invariant*W;
1019  error_norm_squared+=
1020  (latest_guess_invariant-invariant)*
1021  (latest_guess_invariant-invariant)*W;
1022  area+=W;
1023  }
1024 
1025  return area;
1026 
1027  }
1028 
1029 
1031  double square_of_l2_norm()
1032  {
1033 
1034  // Assign dimension
1035  unsigned el_dim=2;
1036  // Initalise
1037  double sum=0.0;
1038 
1039  //Find out how many nodes there are
1040  unsigned n_node = nnode();
1041 
1042  //Find the indices at which the local velocities are stored
1043  unsigned u_nodal_index[el_dim+1];
1044  for(unsigned i=0;i<el_dim+1;i++) {u_nodal_index[i] = u_index_axi_nst(i);}
1045 
1046  //Set up memory for the velocity shape fcts
1047  Shape psif(n_node);
1048  DShape dpsidx(n_node,el_dim);
1049 
1050  //Number of integration points
1051  unsigned n_intpt = integral_pt()->nweight();
1052 
1053  //Set the Vector to hold local coordinates
1054  Vector<double> s(el_dim);
1055 
1056  //Loop over the integration points
1057  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1058  {
1059  //Assign values of s
1060  for(unsigned i=0;i<el_dim;i++) s[i] = integral_pt()->knot(ipt,i);
1061 
1062  //Get the integral weight
1063  double w = integral_pt()->weight(ipt);
1064 
1065  // Call the derivatives of the veloc shape functions
1066  // (Derivs not needed but they are free)
1067  double J = this->dshape_eulerian_at_knot(ipt,psif,dpsidx);
1068 
1069  //Premultiply the weights and the Jacobian
1070  double W = w*J;
1071 
1072  //Calculate velocities
1073  Vector<double> interpolated_u(el_dim+1,0.0);
1074 
1075  // Loop over nodes
1076  for(unsigned l=0;l<n_node;l++)
1077  {
1078  //Loop over directions
1079  for(unsigned i=0;i<el_dim+1;i++)
1080  {
1081  //Get the nodal value
1082  double u_value = raw_nodal_value(l,u_nodal_index[i]);
1083  interpolated_u[i] += u_value*psif[l];
1084  }
1085  }
1086 
1087  //Assemble square of L2 norm
1088  for(unsigned i=0;i<el_dim+1;i++)
1089  {
1090  sum+=interpolated_u[i]*interpolated_u[i]*W;
1091  }
1092  }
1093 
1094  return sum;
1095 
1096  }
1097 
1098 
1099  private:
1100 
1103  Vector<Vector<DenseMatrix<double> > > Fixed_point_iteration_guess_for_strain_rate;
1104 
1107  unsigned Aitken_index;
1108 
1109 // hierher moved additional member functions out I don't need them.
1110 // If required they can be re-instated (but only after merging them
1111 // properly -- there's so much code duplication!
1112 //#include "leftover_bits.h"
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
J
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
v2
Map< RowVectorXf > v2(M2.data(), M2.size())
v1
M1<< 1, 2, 3, 4, 5, 6, 7, 8, 9;Map< RowVectorXf > v1(M1.data(), M1.size())
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
max
#define max(a, b)
Definition: datatypes.h:23
s
RealScalar s
Definition: level1_cplx_impl.h:130
k
char char char int int * k
Definition: level2_impl.h:374
DIM
#define DIM
Definition: linearised_navier_stokes_elements.h:44
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
ProblemParameters::flux
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
calibrate.val
val
Definition: calibrate.py:119
calibrate.error
int error
Definition: calibrate.py:297
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
oomph::Problem_Parameter::Critical_strain_rate
double Critical_strain_rate
Critical strain rate.
Definition: axisym_vibrating_shell_non_newtonian.cc:210
oomph::Problem_Parameter::Yield_stress
double Yield_stress
have a generic yield stress, required for the output function
Definition: axisym_vibrating_shell_non_newtonian.cc:149
oomph::QuadTreeNames::W
@ W
Definition: quadtree.h:63
oomph::SecondInvariantHelper::second_invariant
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
Definition: oomph_utilities.cc:163
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
el_dim
unsigned el_dim
dimension
Definition: overloaded_cartesian_element_body.h:30
Use_aitken_extrapolation
bool Use_aitken_extrapolation
Definition: overloaded_element_body.h:42
get_Z2_flux
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery.
Definition: overloaded_element_body.h:745
output
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
Definition: overloaded_element_body.h:490
Aitken_index
unsigned Aitken_index
Definition: overloaded_element_body.h:1107
square_of_norm_of_viscosity
double square_of_norm_of_viscosity(double &norm_squared, double &extrapolated_norm_squared, double &error_norm_squared)
Definition: overloaded_element_body.h:899
square_of_norm_of_strain_invariant
double square_of_norm_of_strain_invariant(double &norm_squared, double &extrapolated_norm_squared, double &error_norm_squared)
Definition: overloaded_element_body.h:835
set_error
void set_error(const double &error)
Set error value for post-processing.
Definition: overloaded_element_body.h:432
square_of_norm_of_fixed_point
double square_of_norm_of_fixed_point(double &norm_squared, double &latest_guess_norm_squared, double &error_norm_squared)
Definition: overloaded_element_body.h:968
update_latest_fixed_point_iteration_guess_for_strain_rate
void update_latest_fixed_point_iteration_guess_for_strain_rate()
Definition: overloaded_element_body.h:126
Error
double Error
Storage for elemental error estimate – used for post-processing.
Definition: overloaded_element_body.h:30
extrapolated_strain_rate
void extrapolated_strain_rate(const unsigned &ipt, DenseMatrix< double > &strainrate) const
Definition: overloaded_element_body.h:189
square_of_l2_norm
double square_of_l2_norm()
Get square of L2 norm of velocity.
Definition: overloaded_element_body.h:1031
nprev_for_extrapolation_of_strain_rate
unsigned & nprev_for_extrapolation_of_strain_rate()
Definition: overloaded_element_body.h:74
disable_aitken_extrapolation
void disable_aitken_extrapolation()
Disable use of Aitken extrapolation.
Definition: overloaded_element_body.h:103
Fixed_point_iteration_guess_for_strain_rate
Vector< Vector< DenseMatrix< double > > > Fixed_point_iteration_guess_for_strain_rate
Definition: overloaded_element_body.h:1103
variable_identifier
std::string variable_identifier()
Return variable identifier.
Definition: overloaded_element_body.h:435
disable_fixed_point_iteration_for_strain_rate
void disable_fixed_point_iteration_for_strain_rate()
Disable use of fixed point iteration.
Definition: overloaded_element_body.h:90
Nprev_for_extrapolation_of_strain_rate
unsigned Nprev_for_extrapolation_of_strain_rate
Definition: overloaded_element_body.h:34
MyTaylorHoodElement
MyTaylorHoodElement()
Definition: overloaded_element_body.h:48
enable_fixed_point_iteration_for_strain_rate
void enable_fixed_point_iteration_for_strain_rate()
Definition: overloaded_element_body.h:82
enable_aitken_extrapolation
void enable_aitken_extrapolation()
Enable use of Aitken extrapolation.
Definition: overloaded_element_body.h:97
latest_fixed_point_iteration_guess_for_strain_rate
void latest_fixed_point_iteration_guess_for_strain_rate(const unsigned &ipt, DenseMatrix< double > &strainrate) const
Definition: overloaded_element_body.h:110
Use_fixed_point_for_strain_rate
bool Use_fixed_point_for_strain_rate
Definition: overloaded_element_body.h:38
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2