foeppl_von_karman_displacement_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 FoepplvonKarman displacement elements
31 
33 
34 #include <iostream>
35 
36 namespace oomph
37 {
38 
39  // Foeppl von Karman displacement equations static data
40 
41  /// Default value for Poisson's ratio
43 
44 //======================================================================
45 /// Self-test: Return 0 for OK
46 //======================================================================
48 {
49 
50  bool passed=true;
51 
52  // Check lower-level stuff
53  if (FiniteElement::self_test()!=0)
54  {
55  passed=false;
56  }
57 
58  // Return verdict
59  if (passed)
60  {
61  return 0;
62  }
63  else
64  {
65  return 1;
66  }
67 
68 }
69 
70 //======================================================================
71 /// Output function:
72 ///
73 /// x,y,w
74 ///
75 /// nplot points in each coordinate direction
76 //======================================================================
78  const unsigned &nplot)
79 {
80 
81  // Vector of local coordinates
82  Vector<double> s(2);
83  // Vector of Eulerian coordinates
84  Vector<double> x(2);
85 
86  // Sigma_{\alpha \beta}
87  DenseMatrix<double> sigma(2,2);
88 
89  // E_{\alpha \beta}
90  DenseMatrix<double> strain(2,2);
91 
92  // Tecplot header info
93  outfile << tecplot_zone_string(nplot);
94 
95  // Loop over plot points
96  unsigned num_plot_points=nplot_points(nplot);
97  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
98  {
99 
100  // Get local coordinates of plot point
101  get_s_plot(iplot,nplot,s);
102 
103  // Get the Eulerian coordinates
104  for(unsigned i=0;i<2;i++)
105  {
106  x[i] = interpolated_x(s,i);
107  outfile << x[i] << " ";
108  }
109 
110  outfile << interpolated_w_fvk(s,0) << " "; // w
111  outfile << interpolated_w_fvk(s,1) << " "; // Laplacian W
112  outfile << interpolated_w_fvk(s,2) << " "; // Ux
113  outfile << interpolated_w_fvk(s,3) << " "; // Uy
114  // outfile << std::endl; // New line
115 
116  // Get the stresses at local coordinate s
117  this->get_stress_and_strain_for_output(s,sigma,strain);
118 
119  // Output stress
120  outfile << sigma(0,0) << " "; // sigma_xx
121  outfile << sigma(0,1) << " "; // sigma_xy
122  outfile << sigma(1,1) << " "; // sigma_yy
123 
124  // Output strain
125  outfile << strain(0,0) << " "; // E_xx
126  outfile << strain(0,1) << " "; // E_xy
127  outfile << strain(1,1) << " "; // E_yy
128  outfile << std::endl;
129 
130  }
131 
132  // Write tecplot footer (e.g. FE connectivity lists)
133  write_tecplot_zone_footer(outfile,nplot);
134 
135 }
136 
137 
138 //======================================================================
139 /// C-style output function:
140 ///
141 /// x,y,w
142 ///
143 /// nplot points in each coordinate direction
144 //======================================================================
146  const unsigned &nplot)
147 {
148  //Vector of local coordinates
149  Vector<double> s(2);
150 
151  // Tecplot header info
152  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
153 
154  // Loop over plot points
155  unsigned num_plot_points=nplot_points(nplot);
156  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
157  {
158  // Get local coordinates of plot point
159  get_s_plot(iplot,nplot,s);
160 
161  for(unsigned i=0;i<2;i++)
162  {
163  fprintf(file_pt,"%g ",interpolated_x(s,i));
164  }
165  fprintf(file_pt,"%g \n",interpolated_w_fvk(s));
166  }
167 
168  // Write tecplot footer (e.g. FE connectivity lists)
169  write_tecplot_zone_footer(file_pt,nplot);
170 }
171 
172 
173 
174 //======================================================================
175  /// Output exact solution
176  ///
177  /// Solution is provided via function pointer.
178  /// Plot at a given number of plot points.
179  ///
180  /// x,y,w_exact
181 //======================================================================
183  const unsigned &nplot,
185  {
186  //Vector of local coordinates
187  Vector<double> s(2);
188 
189  // Vector for coordintes
190  Vector<double> x(2);
191 
192  // Tecplot header info
193  outfile << tecplot_zone_string(nplot);
194 
195  // Exact solution Vector (here a scalar)
196  //Vector<double> exact_soln(1);
197  Vector<double> exact_soln(3);
198 
199  // Loop over plot points
200  unsigned num_plot_points=nplot_points(nplot);
201  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
202  {
203 
204  // Get local coordinates of plot point
205  get_s_plot(iplot,nplot,s);
206 
207  // Get x position as Vector
208  interpolated_x(s,x);
209 
210  // Get exact solution at this point
211  (*exact_soln_pt)(x,exact_soln);
212 
213  //Output x,y,w_exact
214  for(unsigned i=0;i<2;i++)
215  {
216  outfile << x[i] << " ";
217  }
218  outfile << exact_soln[0] << " ";
219  outfile << exact_soln[1] << " ";
220  outfile << exact_soln[2] << std::endl;
221  }
222 
223  // Write tecplot footer (e.g. FE connectivity lists)
224  write_tecplot_zone_footer(outfile,nplot);
225  }
226 
227 
228 
229 
230 //======================================================================
231 /// Validate against exact solution
232 ///
233 /// Solution is provided via function pointer.
234 /// Plot error at a given number of plot points.
235 ///
236 //======================================================================
239  double& error, double& norm)
240  {
241  // Initialise
242  error=0.0;
243  norm=0.0;
244 
245  //Vector of local coordinates
246  Vector<double> s(2);
247 
248  // Vector for coordintes
249  Vector<double> x(2);
250 
251  //Find out how many nodes there are in the element
252  unsigned n_node = nnode();
253 
254  Shape psi(n_node);
255 
256  //Set the value of n_intpt
257  unsigned n_intpt = integral_pt()->nweight();
258 
259  // Tecplot
260  outfile << "ZONE" << std::endl;
261 
262  // Exact solution Vector (here a scalar)
263  //Vector<double> exact_soln(1);
264  Vector<double> exact_soln(3);
265 
266  //Loop over the integration points
267  for(unsigned ipt=0;ipt<n_intpt;ipt++)
268  {
269 
270  //Assign values of s
271  for(unsigned i=0;i<2;i++)
272  {
273  s[i] = integral_pt()->knot(ipt,i);
274  }
275 
276  //Get the integral weight
277  double w = integral_pt()->weight(ipt);
278 
279  // Get jacobian of mapping
280  double J=J_eulerian(s);
281 
282  //Premultiply the weights and the Jacobian
283  double W = w*J;
284 
285  // Get x position as Vector
286  interpolated_x(s,x);
287 
288  // Get FE function value
289  double w_fe=interpolated_w_fvk(s);
290 
291  // Get exact solution at this point
292  (*exact_soln_pt)(x,exact_soln);
293 
294  //Output x,y,error
295  for(unsigned i=0;i<2;i++)
296  {
297  outfile << x[i] << " ";
298  }
299  outfile << exact_soln[0] << " " << exact_soln[0]-w_fe << std::endl;
300 
301  // Add to error and norm
302  norm+=exact_soln[0]*exact_soln[0]*W;
303  error+=(exact_soln[0]-w_fe)*(exact_soln[0]-w_fe)*W;
304 
305  }
306 
307  }
308 
309 } // namespace oomph
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
cstr elem_len * i
Definition: cfortran.h:607
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
double interpolated_w_fvk(const Vector< double > &s, unsigned index=0) const
Return FE representation of function value w_fvk(s) at local coordinate s (by default - if index > 0...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
void get_stress_and_strain_for_output(const Vector< double > &s, DenseMatrix< double > &sigma, DenseMatrix< double > &strain)
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,w_exact at n_plot^DIM plot points.
void output(std::ostream &outfile)
Output with default number of plot points.
static double Default_Nu_Value
Default value for Poisson's ratio.