poisson_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 Poisson elements
31 #include "poisson_elements.h"
32
33
34 namespace oomph
35 {
36
37
38 //======================================================================
39 /// Set the data for the number of Variables at each node, always one
40 /// in every case
41 //======================================================================
42  template<unsigned DIM, unsigned NNODE_1D>
43  const unsigned QPoissonElement<DIM,NNODE_1D>::Initial_Nvalue = 1;
44
45
46 //======================================================================
47 /// Compute element residual Vector and/or element Jacobian matrix
48 ///
49 /// flag=1: compute both
50 /// flag=0: compute only residual Vector
51 ///
52 /// Pure version without hanging nodes
53 //======================================================================
54 template <unsigned DIM>
57  DenseMatrix<double> &jacobian,
58  const unsigned& flag)
59 {
60  //Find out how many nodes there are
61  const unsigned n_node = nnode();
62
63  //Set up memory for the shape and test functions
64  Shape psi(n_node), test(n_node);
65  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
66
67  //Index at which the poisson unknown is stored
68  const unsigned u_nodal_index = u_index_poisson();
69
70  //Set the value of n_intpt
71  const unsigned n_intpt = integral_pt()->nweight();
72
73  //Integers to store the local equation and unknown numbers
74  int local_eqn=0, local_unknown=0;
75
76  //Loop over the integration points
77  for(unsigned ipt=0;ipt<n_intpt;ipt++)
78  {
79  //Get the integral weight
80  double w = integral_pt()->weight(ipt);
81
82  //Call the derivatives of the shape and test functions
83  double J = dshape_and_dtest_eulerian_at_knot_poisson(ipt,psi,dpsidx,
84  test,dtestdx);
85
86  //Premultiply the weights and the Jacobian
87  double W = w*J;
88
89  //Calculate local values of unknown
90  //Allocate and initialise to zero
91  double interpolated_u=0.0;
92  Vector<double> interpolated_x(DIM,0.0);
93  Vector<double> interpolated_dudx(DIM,0.0);
94
95  //Calculate function value and derivatives:
96  //-----------------------------------------
97  // Loop over nodes
98  for(unsigned l=0;l<n_node;l++)
99  {
100  //Get the nodal value of the poisson unknown
101  double u_value = raw_nodal_value(l,u_nodal_index);
102  interpolated_u += u_value*psi(l);
103  // Loop over directions
104  for(unsigned j=0;j<DIM;j++)
105  {
106  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
107  interpolated_dudx[j] += u_value*dpsidx(l,j);
108  }
109  }
110
111
112  //Get source function
113  //-------------------
114  double source;
115  get_source_poisson(ipt,interpolated_x,source);
116
117  // Assemble residuals and Jacobian
118  //--------------------------------
119
120  // Loop over the test functions
121  for(unsigned l=0;l<n_node;l++)
122  {
123  //Get the local equation
124  local_eqn = nodal_local_eqn(l,u_nodal_index);
125  // IF it's not a boundary condition
126  if(local_eqn >= 0)
127  {
128  // Add body force/source term here
129  residuals[local_eqn] += source*test(l)*W;
130
131  // The Poisson bit itself
132  for(unsigned k=0;k<DIM;k++)
133  {
134  residuals[local_eqn] += interpolated_dudx[k]*dtestdx(l,k)*W;
135  }
136
137  // Calculate the jacobian
138  //-----------------------
139  if(flag)
140  {
141  //Loop over the velocity shape functions again
142  for(unsigned l2=0;l2<n_node;l2++)
143  {
144  local_unknown = nodal_local_eqn(l2,u_nodal_index);
145  //If at a non-zero degree of freedom add in the entry
146  if(local_unknown >= 0)
147  {
148  //Add contribution to Elemental Matrix
149  for(unsigned i=0;i<DIM;i++)
150  {
151  jacobian(local_eqn,local_unknown)
152  += dpsidx(l2,i)*dtestdx(l,i)*W;
153  }
154  }
155  }
156  }
157  }
158  }
159
160  } // End of loop over integration points
161 }
162
163
164
165 //======================================================================
166 /// Compute derivatives of elemental residual vector with respect
167 /// to nodal coordinates (fully analytical).
168 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
169 /// Overloads the FD-based version in the FE base class.
170 //======================================================================
171 template <unsigned DIM>
174  dresidual_dnodal_coordinates)
175 {
176  // Determine number of nodes in element
177  const unsigned n_node = nnode();
178
179  // Set up memory for the shape and test functions
180  Shape psi(n_node), test(n_node);
181  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
182
183  // Deriatives of shape fct derivatives w.r.t. nodal coords
184  RankFourTensor<double> d_dpsidx_dX(DIM,n_node,n_node,DIM);
185  RankFourTensor<double> d_dtestdx_dX(DIM,n_node,n_node,DIM);
186
187  // Derivative of Jacobian of mapping w.r.t. to nodal coords
188  DenseMatrix<double> dJ_dX(DIM,n_node);
189
190  // Derivatives of derivative of u w.r.t. nodal coords
191  RankThreeTensor<double> d_dudx_dX(DIM,n_node,DIM);
192
193  // Source function and its gradient
194  double source;
195  Vector<double> d_source_dx(DIM);
196
197  // Index at which the poisson unknown is stored
198  const unsigned u_nodal_index = u_index_poisson();
199
200  // Determine the number of integration points
201  const unsigned n_intpt = integral_pt()->nweight();
202
203  // Integer to store the local equation number
204  int local_eqn=0;
205
206  // Loop over the integration points
207  for(unsigned ipt=0;ipt<n_intpt;ipt++)
208  {
209  // Get the integral weight
210  double w = integral_pt()->weight(ipt);
211
212  // Call the derivatives of the shape/test functions, as well as the
213  // derivatives of these w.r.t. nodal coordinates and the derivative
214  // of the jacobian of the mapping w.r.t. nodal coordinates
215  const double J = dshape_and_dtest_eulerian_at_knot_poisson(
216  ipt,psi,dpsidx,d_dpsidx_dX,test,dtestdx,d_dtestdx_dX,dJ_dX);
217
218  // Calculate local values
219  // Allocate and initialise to zero
220  Vector<double> interpolated_x(DIM,0.0);
221  Vector<double> interpolated_dudx(DIM,0.0);
222
223  // Calculate function value and derivatives:
224  // -----------------------------------------
225  // Loop over nodes
226  for(unsigned l=0;l<n_node;l++)
227  {
228  // Get the nodal value of the Poisson unknown
229  double u_value = raw_nodal_value(l,u_nodal_index);
230
231  // Loop over directions
232  for(unsigned i=0;i<DIM;i++)
233  {
234  interpolated_x[i] += raw_nodal_position(l,i)*psi(l);
235  interpolated_dudx[i] += u_value*dpsidx(l,i);
236  }
237  }
238
239  // Calculate derivative of du/dx_i w.r.t. nodal positions X_{pq}
240  for(unsigned q=0;q<n_node;q++)
241  {
242  // Loop over coordinate directions
243  for(unsigned p=0;p<DIM;p++)
244  {
245  for(unsigned i=0;i<DIM;i++)
246  {
247  double aux=0.0;
248  for(unsigned j=0;j<n_node;j++)
249  {
250  aux += raw_nodal_value(j,u_nodal_index)*d_dpsidx_dX(p,q,j,i);
251  }
252  d_dudx_dX(p,q,i) = aux;
253  }
254  }
255  }
256
257  // Get source function
258  get_source_poisson(ipt,interpolated_x,source);
259
260  // Get gradient of source function
262
263  // Assemble d res_{local_eqn} / d X_{pq}
264  // -------------------------------------
265
266  // Loop over the test functions
267  for(unsigned l=0;l<n_node;l++)
268  {
269  // Get the local equation
270  local_eqn = nodal_local_eqn(l,u_nodal_index);
271
272  // IF it's not a boundary condition
273  if(local_eqn >= 0)
274  {
275  // Loop over coordinate directions
276  for(unsigned p=0;p<DIM;p++)
277  {
278  // Loop over nodes
279  for(unsigned q=0;q<n_node;q++)
280  {
281  double sum = source*test(l)*dJ_dX(p,q)
282  + d_source_dx[p]*test(l)*psi(q)*J;
283
284  for(unsigned i=0;i<DIM;i++)
285  {
286  sum += interpolated_dudx[i]*(dtestdx(l,i)*dJ_dX(p,q) +
287  d_dtestdx_dX(p,q,l,i)*J)
288  + d_dudx_dX(p,q,i)*dtestdx(l,i)*J;
289  }
290
291  // Multiply through by integration weight
292  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*w;
293  }
294  }
295  }
296  }
297  } // End of loop over integration points
298 }
299
300
301
302 //======================================================================
303 /// Self-test: Return 0 for OK
304 //======================================================================
305 template <unsigned DIM>
307 {
308
309  bool passed=true;
310
311  // Check lower-level stuff
312  if (FiniteElement::self_test()!=0)
313  {
314  passed=false;
315  }
316
317  // Return verdict
318  if (passed)
319  {
320  return 0;
321  }
322  else
323  {
324  return 1;
325  }
326
327 }
328
329
330
331 //======================================================================
332 /// Output function:
333 ///
334 /// x,y,u or x,y,z,u
335 ///
336 /// nplot points in each coordinate direction
337 //======================================================================
338 template <unsigned DIM>
339 void PoissonEquations<DIM>::output(std::ostream &outfile,
340  const unsigned &nplot)
341 {
342
343  //Vector of local coordinates
344  Vector<double> s(DIM);
345
346  // Tecplot header info
347  outfile << tecplot_zone_string(nplot);
348
349  // Loop over plot points
350  unsigned num_plot_points=nplot_points(nplot);
351  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
352  {
353
354  // Get local coordinates of plot point
355  get_s_plot(iplot,nplot,s);
356
357  for(unsigned i=0;i<DIM;i++)
358  {
359  outfile << interpolated_x(s,i) << " ";
360  }
361  outfile << interpolated_u_poisson(s) << std::endl;
362
363  }
364
365  // Write tecplot footer (e.g. FE connectivity lists)
366  write_tecplot_zone_footer(outfile,nplot);
367
368 }
369
370
371 //======================================================================
372 /// C-style output function:
373 ///
374 /// x,y,u or x,y,z,u
375 ///
376 /// nplot points in each coordinate direction
377 //======================================================================
378 template <unsigned DIM>
379 void PoissonEquations<DIM>::output(FILE* file_pt,
380  const unsigned &nplot)
381 {
382  //Vector of local coordinates
383  Vector<double> s(DIM);
384
385  // Tecplot header info
386  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
387
388  // Loop over plot points
389  unsigned num_plot_points=nplot_points(nplot);
390  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
391  {
392  // Get local coordinates of plot point
393  get_s_plot(iplot,nplot,s);
394
395  for(unsigned i=0;i<DIM;i++)
396  {
397  fprintf(file_pt,"%g ",interpolated_x(s,i));
398  }
399  fprintf(file_pt,"%g \n",interpolated_u_poisson(s));
400  }
401
402  // Write tecplot footer (e.g. FE connectivity lists)
403  write_tecplot_zone_footer(file_pt,nplot);
404 }
405
406
407
408 //======================================================================
409  /// Output exact solution
410  ///
411  /// Solution is provided via function pointer.
412  /// Plot at a given number of plot points.
413  ///
414  /// x,y,u_exact or x,y,z,u_exact
415 //======================================================================
416 template <unsigned DIM>
417 void PoissonEquations<DIM>::output_fct(std::ostream &outfile,
418  const unsigned &nplot,
420 {
421  //Vector of local coordinates
422  Vector<double> s(DIM);
423
424  // Vector for coordintes
425  Vector<double> x(DIM);
426
427  // Tecplot header info
428  outfile << tecplot_zone_string(nplot);
429
430  // Exact solution Vector (here a scalar)
431  Vector<double> exact_soln(1);
432
433  // Loop over plot points
434  unsigned num_plot_points=nplot_points(nplot);
435  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
436  {
437
438  // Get local coordinates of plot point
439  get_s_plot(iplot,nplot,s);
440
441  // Get x position as Vector
442  interpolated_x(s,x);
443
444  // Get exact solution at this point
445  (*exact_soln_pt)(x,exact_soln);
446
447  //Output x,y,...,u_exact
448  for(unsigned i=0;i<DIM;i++)
449  {
450  outfile << x[i] << " ";
451  }
452  outfile << exact_soln[0] << std::endl;
453  }
454
455  // Write tecplot footer (e.g. FE connectivity lists)
456  write_tecplot_zone_footer(outfile,nplot);
457 }
458
459
460
461
462 //======================================================================
463  /// Validate against exact solution
464  ///
465  /// Solution is provided via function pointer.
466  /// Plot error at a given number of plot points.
467  ///
468 //======================================================================
469 template <unsigned DIM>
470 void PoissonEquations<DIM>::compute_error(std::ostream &outfile,
472  double& error, double& norm)
473 {
474
475  // Initialise
476  error=0.0;
477  norm=0.0;
478
479  //Vector of local coordinates
480  Vector<double> s(DIM);
481
482  // Vector for coordintes
483  Vector<double> x(DIM);
484
485  //Find out how many nodes there are in the element
486  unsigned n_node = nnode();
487
488  Shape psi(n_node);
489
490  //Set the value of n_intpt
491  unsigned n_intpt = integral_pt()->nweight();
492
493  // Tecplot
494  outfile << "ZONE" << std::endl;
495
496  // Exact solution Vector (here a scalar)
497  Vector<double> exact_soln(1);
498
499  //Loop over the integration points
500  for(unsigned ipt=0;ipt<n_intpt;ipt++)
501  {
502
503  //Assign values of s
504  for(unsigned i=0;i<DIM;i++)
505  {
506  s[i] = integral_pt()->knot(ipt,i);
507  }
508
509  //Get the integral weight
510  double w = integral_pt()->weight(ipt);
511
512  // Get jacobian of mapping
513  double J=J_eulerian(s);
514
515  //Premultiply the weights and the Jacobian
516  double W = w*J;
517
518  // Get x position as Vector
519  interpolated_x(s,x);
520
521  // Get FE function value
522  double u_fe=interpolated_u_poisson(s);
523
524  // Get exact solution at this point
525  (*exact_soln_pt)(x,exact_soln);
526
527  //Output x,y,...,error
528  for(unsigned i=0;i<DIM;i++)
529  {
530  outfile << x[i] << " ";
531  }
532  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
533
534  // Add to error and norm
535  norm+=exact_soln[0]*exact_soln[0]*W;
536  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
537
538  }
539 }
540
541
542
543
544
545 //====================================================================
546 // Force build of templates
547 //====================================================================
548 template class QPoissonElement<1,2>;
549 template class QPoissonElement<1,3>;
550 template class QPoissonElement<1,4>;
551
552 template class QPoissonElement<2,2>;
553 template class QPoissonElement<2,3>;
554 template class QPoissonElement<2,4>;
555
556 template class QPoissonElement<3,2>;
557 template class QPoissonElement<3,3>;
558 template class QPoissonElement<3,4>;
559
560 }
cstr elem_len * i
Definition: cfortran.h:607
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
virtual void fill_in_generic_residual_contribution_poisson(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
A Rank 4 Tensor class.
Definition: matrices.h:1625
void output(std::ostream &outfile)
Output with default number of plot points.
A Rank 3 Tensor class.
Definition: matrices.h:1337
static char t char * s
Definition: cfortran.h:572
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Overwrites default implementation in FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points.
unsigned self_test()
Self-test: Return 0 for OK.