my_taylor_hood_elements.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// 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 #ifndef OOMPH_MY_TAYLOR_HOOD_ELEMENTS
31 #define OOMPH_MY_TAYLOR_HOOD_ELEMENTS
32 
33 namespace oomph
34 {
35 
36 //==============================================================
37 /// Overload TaylorHood element to modify output
38 //==============================================================
40  public virtual PseudoSolidNodeUpdateElement<TTaylorHoodElement<2>,
41  TPVDElement<2,3> >
42  {
43 
44  private:
45 
46  /// Storage for elemental error estimate -- used for post-processing
47  double Error;
48 
49  public:
50 
51  /// Constructor initialise error
53  {
54  Error=0.0;
55  }
56 
57  /// Set error value for post-processing
58  void set_error(const double& error){Error=error;}
59 
60  /// Return variable identifier
61  std::string variable_identifier()
62  {
63  std::string txt="VARIABLES=";
64  txt+="\"x\",";
65  txt+="\"y\",";
66  txt+="\"u\",";
67  txt+="\"v\",";
68  txt+="\"p\",";
69  txt+="\"du/dt\",";
70  txt+="\"dv/dt\",";
71  txt+="\"u_m\",";
72  txt+="\"v_m\",";
73  txt+="\"x_h1\",";
74  txt+="\"y_h1\",";
75  txt+="\"x_h2\",";
76  txt+="\"y_h2\",";
77  txt+="\"u_h1\",";
78  txt+="\"v_h1\",";
79  txt+="\"u_h2\",";
80  txt+="\"v_h2\",";
81  txt+="\"error\",";
82  txt+="\"size\",";
83  txt+="\n";
84  return txt;
85  }
86 
87 
88  /// Overload output function
89  void output(std::ostream &outfile,
90  const unsigned &nplot)
91  {
92 
93  // Assign dimension
94  unsigned el_dim=2;
95 
96  // Vector of local coordinates
97  Vector<double> s(el_dim);
98 
99  // Acceleration
100  Vector<double> dudt(el_dim);
101 
102  // Mesh elocity
103  Vector<double> mesh_veloc(el_dim,0.0);
104 
105  // Tecplot header info
106  outfile << tecplot_zone_string(nplot);
107 
108  // Find out how many nodes there are
109  unsigned n_node = nnode();
110 
111  //Set up memory for the shape functions
112  Shape psif(n_node);
113  DShape dpsifdx(n_node,el_dim);
114 
115  // Loop over plot points
116  unsigned num_plot_points=nplot_points(nplot);
117  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
118  {
119 
120  // Get local coordinates of plot point
121  get_s_plot(iplot,nplot,s);
122 
123  //Call the derivatives of the shape and test functions
124  dshape_eulerian(s,psif,dpsifdx);
125 
126  //Allocate storage
127  Vector<double> mesh_veloc(el_dim);
128  Vector<double> dudt(el_dim);
129  Vector<double> dudt_ALE(el_dim);
130  DenseMatrix<double> interpolated_dudx(el_dim,el_dim);
131 
132  //Initialise everything to zero
133  for(unsigned i=0;i<el_dim;i++)
134  {
135  mesh_veloc[i]=0.0;
136  dudt[i]=0.0;
137  dudt_ALE[i]=0.0;
138  for(unsigned j=0;j<el_dim;j++)
139  {
140  interpolated_dudx(i,j) = 0.0;
141  }
142  }
143 
144  //Calculate velocities and derivatives
145 
146  //Loop over directions
147  for(unsigned i=0;i<el_dim;i++)
148  {
149  //Get the index at which velocity i is stored
150  unsigned u_nodal_index = u_index_nst(i);
151  // Loop over nodes
152  for(unsigned l=0;l<n_node;l++)
153  {
154  dudt[i]+=du_dt_nst(l,u_nodal_index)*psif[l];
155  mesh_veloc[i]+=dnodal_position_dt(l,i)*psif[l];
156 
157  //Loop over derivative directions for velocity gradients
158  for(unsigned j=0;j<el_dim;j++)
159  {
160  interpolated_dudx(i,j) += nodal_value(l,u_nodal_index)*
161  dpsifdx(l,j);
162  }
163  }
164  }
165 
166 
167  // Get dudt in ALE form (incl mesh veloc)
168  for(unsigned i=0;i<el_dim;i++)
169  {
170  dudt_ALE[i]=dudt[i];
171  for (unsigned k=0;k<el_dim;k++)
172  {
173  dudt_ALE[i]-=mesh_veloc[k]*interpolated_dudx(i,k);
174  }
175  }
176 
177 
178  // Coordinates
179  for(unsigned i=0;i<el_dim;i++)
180  {
181  outfile << interpolated_x(s,i) << " ";
182  }
183 
184  // Velocities
185  for(unsigned i=0;i<el_dim;i++)
186  {
187  outfile << interpolated_u_nst(s,i) << " ";
188  }
189 
190  // Pressure
191  outfile << interpolated_p_nst(s) << " ";
192 
193  // Accelerations
194  for(unsigned i=0;i<el_dim;i++)
195  {
196  outfile << dudt_ALE[i] << " ";
197  }
198 
199  // Mesh velocity
200  for(unsigned i=0;i<el_dim;i++)
201  {
202  outfile << mesh_veloc[i] << " ";
203  }
204 
205  // History values of coordinates
206  unsigned n_prev=node_pt(0)->position_time_stepper_pt()->ntstorage();
207  for (unsigned t=1;t<n_prev;t++)
208  {
209  for(unsigned i=0;i<el_dim;i++)
210  {
211  outfile << interpolated_x(t,s,i) << " ";
212  }
213  }
214 
215  // History values of velocities
216  n_prev=node_pt(0)->time_stepper_pt()->ntstorage();
217  for (unsigned t=1;t<n_prev;t++)
218  {
219  for(unsigned i=0;i<el_dim;i++)
220  {
221  outfile << interpolated_u_nst(t,s,i) << " ";
222  }
223  }
224 
225  outfile << Error << " "
226  << size() << std::endl;
227  }
228 
229  // Write tecplot footer (e.g. FE connectivity lists)
230  write_tecplot_zone_footer(outfile,nplot);
231  }
232 
233 
234  /// Get square of L2 norm of velocity
236  {
237 
238  // Assign dimension
239  unsigned el_dim=2;
240  // Initalise
241  double sum=0.0;
242 
243  //Find out how many nodes there are
244  unsigned n_node = nnode();
245 
246  //Find the indices at which the local velocities are stored
247  unsigned u_nodal_index[el_dim];
248  for(unsigned i=0;i<el_dim;i++) {u_nodal_index[i] = u_index_nst(i);}
249 
250  //Set up memory for the velocity shape fcts
251  Shape psif(n_node);
252  DShape dpsidx(n_node,el_dim);
253 
254  //Number of integration points
255  unsigned n_intpt = integral_pt()->nweight();
256 
257  //Set the Vector to hold local coordinates
258  Vector<double> s(el_dim);
259 
260  //Loop over the integration points
261  for(unsigned ipt=0;ipt<n_intpt;ipt++)
262  {
263  //Assign values of s
264  for(unsigned i=0;i<el_dim;i++) s[i] = integral_pt()->knot(ipt,i);
265 
266  //Get the integral weight
267  double w = integral_pt()->weight(ipt);
268 
269  // Call the derivatives of the veloc shape functions
270  // (Derivs not needed but they are free)
271  double J = this->dshape_eulerian_at_knot(ipt,psif,dpsidx);
272 
273  //Premultiply the weights and the Jacobian
274  double W = w*J;
275 
276  //Calculate velocities
277  Vector<double> interpolated_u(el_dim,0.0);
278 
279  // Loop over nodes
280  for(unsigned l=0;l<n_node;l++)
281  {
282  //Loop over directions
283  for(unsigned i=0;i<el_dim;i++)
284  {
285  //Get the nodal value
286  double u_value = raw_nodal_value(l,u_nodal_index[i]);
287  interpolated_u[i] += u_value*psif[l];
288  }
289  }
290 
291  //Assemble square of L2 norm
292  for(unsigned i=0;i<el_dim;i++)
293  {
294  sum+=interpolated_u[i]*interpolated_u[i]*W;
295  }
296  }
297 
298  return sum;
299 
300  }
301 
302  };
303 
304 
305 
306 
307 
308 //=======================================================================
309 /// Face geometry for element is the same as that for the underlying
310 /// wrapped element
311 //=======================================================================
312  template<>
313  class FaceGeometry<MyTaylorHoodElement>
314  : public virtual SolidTElement<1,3>
315  {
316  public:
317  FaceGeometry() : SolidTElement<1,3>() {}
318  };
319 
320 //=======================================================================
321 /// Face geometry for element is the same as that for the underlying
322 /// wrapped element
323 //=======================================================================
324  template<>
325  class FaceGeometry<FaceGeometry<MyTaylorHoodElement> >
326  : public virtual PointElement
327  {
328  public:
329  FaceGeometry() : PointElement() {}
330  };
331 
332 
333 }
334 #endif
double Error
Storage for elemental error estimate – used for post-processing.
void set_error(const double &error)
Set error value for post-processing.
MyTaylorHoodElement()
Constructor initialise error.
Overload TaylorHood element to modify output.
double square_of_l2_norm()
Get square of L2 norm of velocity.
std::string variable_identifier()
Return variable identifier.
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.