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$
7 //LIC//
8 //LIC// $LastChangedDate$
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>
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
261  get_source_gradient_poisson(ipt,interpolated_x,d_source_dx);
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 /// Compute norm of the solution
464 //=======================================================================
465 template <unsigned DIM>
467 {
468 
469  // Initialise
470  norm=0.0;
471 
472  //Vector of local coordinates
473  Vector<double> s(DIM);
474 
475  // Solution
476  double u=0.0;
477 
478  //Find out how many nodes there are in the element
479  unsigned n_node = this->nnode();
480 
481  Shape psi(n_node);
482 
483  //Set the value of n_intpt
484  unsigned n_intpt = this->integral_pt()->nweight();
485 
486  //Loop over the integration points
487  for(unsigned ipt=0;ipt<n_intpt;ipt++)
488  {
489 
490  //Assign values of s
491  for(unsigned i=0;i<DIM;i++)
492  {
493  s[i] = this->integral_pt()->knot(ipt,i);
494  }
495 
496  //Get the integral weight
497  double w = this->integral_pt()->weight(ipt);
498 
499  // Get jacobian of mapping
500  double J=this->J_eulerian(s);
501 
502  //Premultiply the weights and the Jacobian
503  double W = w*J;
504 
505  // Get FE function value
506  u=this->interpolated_u_poisson(s);
507 
508  // Add to norm
509  norm+=u*u*W;
510  }
511 }
512 
513 //======================================================================
514  /// Validate against exact solution
515  ///
516  /// Solution is provided via function pointer.
517  /// Plot error at a given number of plot points.
518  ///
519 //======================================================================
520 template <unsigned DIM>
521 void PoissonEquations<DIM>::compute_error(std::ostream &outfile,
523  double& error, double& norm)
524 {
525 
526  // Initialise
527  error=0.0;
528  norm=0.0;
529 
530  //Vector of local coordinates
531  Vector<double> s(DIM);
532 
533  // Vector for coordintes
534  Vector<double> x(DIM);
535 
536  //Find out how many nodes there are in the element
537  unsigned n_node = nnode();
538 
539  Shape psi(n_node);
540 
541  //Set the value of n_intpt
542  unsigned n_intpt = integral_pt()->nweight();
543 
544  // Tecplot
545  outfile << "ZONE" << std::endl;
546 
547  // Exact solution Vector (here a scalar)
548  Vector<double> exact_soln(1);
549 
550  //Loop over the integration points
551  for(unsigned ipt=0;ipt<n_intpt;ipt++)
552  {
553 
554  //Assign values of s
555  for(unsigned i=0;i<DIM;i++)
556  {
557  s[i] = integral_pt()->knot(ipt,i);
558  }
559 
560  //Get the integral weight
561  double w = integral_pt()->weight(ipt);
562 
563  // Get jacobian of mapping
564  double J=J_eulerian(s);
565 
566  //Premultiply the weights and the Jacobian
567  double W = w*J;
568 
569  // Get x position as Vector
570  interpolated_x(s,x);
571 
572  // Get FE function value
573  double u_fe=interpolated_u_poisson(s);
574 
575  // Get exact solution at this point
576  (*exact_soln_pt)(x,exact_soln);
577 
578  //Output x,y,...,error
579  for(unsigned i=0;i<DIM;i++)
580  {
581  outfile << x[i] << " ";
582  }
583  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
584 
585  // Add to error and norm
586  norm+=exact_soln[0]*exact_soln[0]*W;
587  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
588 
589  }
590 }
591 
592 
593 
594 
595 
596 //====================================================================
597 // Force build of templates
598 //====================================================================
599 template class QPoissonElement<1,2>;
600 template class QPoissonElement<1,3>;
601 template class QPoissonElement<1,4>;
602 
603 template class QPoissonElement<2,2>;
604 template class QPoissonElement<2,3>;
605 template class QPoissonElement<2,4>;
606 
607 template class QPoissonElement<3,2>;
608 template class QPoissonElement<3,3>;
609 template class QPoissonElement<3,4>;
610 
611 }
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
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:4333
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.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1723
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 compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void compute_norm(double &norm)
Compute norm of solution: square of the L2 norm.
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.