poisson_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 //Header file for Poisson elements
31 #ifndef OOMPH_POISSON_ELEMENTS_HEADER
32 #define OOMPH_POISSON_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 #include<sstream>
41 
42 //OOMPH-LIB headers
43 #include "../generic/projection.h"
44 #include "../generic/nodes.h"
45 #include "../generic/Qelements.h"
46 #include "../generic/oomph_utilities.h"
47 
48 
49 
50 namespace oomph
51 {
52 
53 //=============================================================
54 /// A class for all isoparametric elements that solve the
55 /// Poisson equations.
56 /// \f[
57 /// \frac{\partial^2 u}{\partial x_i^2} = f(x_j)
58 /// \f]
59 /// This contains the generic maths. Shape functions, geometric
60 /// mapping etc. must get implemented in derived class.
61 //=============================================================
62 template <unsigned DIM>
63 class PoissonEquations : public virtual FiniteElement
64 {
65 
66 public:
67 
68  /// \short Function pointer to source function fct(x,f(x)) --
69  /// x is a Vector!
70  typedef void (*PoissonSourceFctPt)(const Vector<double>& x, double& f);
71 
72 
73  /// \short Function pointer to gradient of source function fct(x,g(x)) --
74  /// x is a Vector!
75  typedef void (*PoissonSourceFctGradientPt)(const Vector<double>& x,
76  Vector<double>& gradient);
77 
78 
79  /// Constructor (must initialise the Source_fct_pt to null)
81  {}
82 
83  /// Broken copy constructor
85  {
86  BrokenCopy::broken_copy("PoissonEquations");
87  }
88 
89  /// Broken assignment operator
91  {
92  BrokenCopy::broken_assign("PoissonEquations");
93  }
94 
95  /// \short Return the index at which the unknown value
96  /// is stored. The default value, 0, is appropriate for single-physics
97  /// problems, when there is only one variable, the value that satisfies
98  /// the poisson equation.
99  /// In derived multi-physics elements, this function should be overloaded
100  /// to reflect the chosen storage scheme. Note that these equations require
101  /// that the unknown is always stored at the same index at each node.
102  virtual inline unsigned u_index_poisson() const {return 0;}
103 
104  /// \short Number of scalars/fields output by this element. Reimplements
105  /// broken virtual function in base class.
106  unsigned nscalar_paraview() const
107  {
108  return 1;
109  }
110 
111  /// \short Write values of the i-th scalar field at the plot points. Needs
112  /// to be implemented for each new specific element type.
113  void scalar_value_paraview(std::ofstream& file_out,
114  const unsigned& i,
115  const unsigned& nplot) const
116  {
117 #ifdef PARANOID
118  if (i!=0)
119  {
120  std::stringstream error_stream;
121  error_stream
122  << "Poisson elements only store a single field so i must be 0 rather"
123  << " than " << i << std::endl;
124  throw OomphLibError(
125  error_stream.str(),
126  OOMPH_CURRENT_FUNCTION,
127  OOMPH_EXCEPTION_LOCATION);
128  }
129 #endif
130 
131  unsigned local_loop=this->nplot_points_paraview(nplot);
132  for(unsigned j=0;j<local_loop;j++)
133  {
134  // Get the local coordinate of the required plot point
135  Vector<double> s(DIM);
136  this->get_s_plot(j,nplot,s);
137 
138  file_out << this->interpolated_u_poisson(s)
139  << std::endl;
140  }
141  }
142 
143  /// \short Name of the i-th scalar field. Default implementation
144  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
145  /// overloaded with more meaningful names in specific elements.
146  std::string scalar_name_paraview(const unsigned& i) const
147  {
148 
149 #ifdef PARANOID
150  if (i!=0)
151  {
152  std::stringstream error_stream;
153  error_stream
154  << "Poisson elements only store a single field so i must be 0 rather"
155  << " than " << i << std::endl;
156  throw OomphLibError(
157  error_stream.str(),
158  OOMPH_CURRENT_FUNCTION,
159  OOMPH_EXCEPTION_LOCATION);
160  }
161 #endif
162 
163  return "Poisson solution";
164  }
165 
166  /// Output with default number of plot points
167  void output(std::ostream &outfile)
168  {
169  const unsigned n_plot=5;
170  output(outfile,n_plot);
171  }
172 
173  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
174  /// n_plot^DIM plot points
175  void output(std::ostream &outfile, const unsigned &n_plot);
176 
177  /// C_style output with default number of plot points
178  void output(FILE* file_pt)
179  {
180  const unsigned n_plot=5;
181  output(file_pt,n_plot);
182  }
183 
184  /// \short C-style output FE representation of soln: x,y,u or x,y,z,u at
185  /// n_plot^DIM plot points
186  void output(FILE* file_pt, const unsigned &n_plot);
187 
188  /// Output exact soln: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
189  void output_fct(std::ostream &outfile, const unsigned &n_plot,
191 
192  /// \short Output exact soln: x,y,u_exact or x,y,z,u_exact at
193  /// n_plot^DIM plot points (dummy time-dependent version to
194  /// keep intel compiler happy)
195  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
196  const double& time,
198  exact_soln_pt)
199  {
200  throw OomphLibError(
201  "There is no time-dependent output_fct() for Poisson elements ",
202  OOMPH_CURRENT_FUNCTION,
203  OOMPH_EXCEPTION_LOCATION);
204  }
205 
206 
207  /// Get error against and norm of exact solution
208  void compute_error(std::ostream &outfile,
210  double& error, double& norm);
211 
212 
213  /// Dummy, time dependent error checker
214  void compute_error(std::ostream &outfile,
216  const double& time, double& error, double& norm)
217  {
218  throw OomphLibError(
219  "There is no time-dependent compute_error() for Poisson elements",
220  OOMPH_CURRENT_FUNCTION,
221  OOMPH_EXCEPTION_LOCATION);
222  }
223 
224  /// Access function: Pointer to source function
226 
227  /// Access function: Pointer to source function. Const version
229 
230  /// Access function: Pointer to gradient of source function
232  {return Source_fct_gradient_pt;}
233 
234  /// Access function: Pointer to gradient source function. Const version
236  {return Source_fct_gradient_pt;}
237 
238 
239  /// Get source term at (Eulerian) position x. This function is
240  /// virtual to allow overloading in multi-physics problems where
241  /// the strength of the source function might be determined by
242  /// another system of equations.
243  inline virtual void get_source_poisson(const unsigned& ipt,
244  const Vector<double>& x,
245  double& source) const
246  {
247  //If no source function has been set, return zero
248  if(Source_fct_pt==0) {source = 0.0;}
249  else
250  {
251  // Get source strength
252  (*Source_fct_pt)(x,source);
253  }
254  }
255 
256 
257  /// Get gradient of source term at (Eulerian) position x. This function is
258  /// virtual to allow overloading in multi-physics problems where
259  /// the strength of the source function might be determined by
260  /// another system of equations. Computed via function pointer
261  /// (if set) or by finite differencing (default)
262  inline virtual void get_source_gradient_poisson(
263  const unsigned& ipt,
264  const Vector<double>& x,
265  Vector<double>& gradient) const
266  {
267  //If no gradient function has been set, FD it
269  {
270  // Reference value
271  double source=0.0;
272  get_source_poisson(ipt,x,source);
273 
274  // FD it
276  double source_pls=0.0;
277  Vector<double> x_pls(x);
278  for (unsigned i=0;i<DIM;i++)
279  {
280  x_pls[i]+=eps_fd;
281  get_source_poisson(ipt,x_pls,source_pls);
282  gradient[i]=(source_pls-source)/eps_fd;
283  x_pls[i]=x[i];
284  }
285  }
286  else
287  {
288  // Get gradient
289  (*Source_fct_gradient_pt)(x,gradient);
290  }
291  }
292 
293 
294 
295 
296  /// Get flux: flux[i] = du/dx_i
297  void get_flux(const Vector<double>& s, Vector<double>& flux) const
298  {
299  //Find out how many nodes there are in the element
300  const unsigned n_node = nnode();
301 
302  //Get the index at which the unknown is stored
303  const unsigned u_nodal_index = u_index_poisson();
304 
305  //Set up memory for the shape and test functions
306  Shape psi(n_node);
307  DShape dpsidx(n_node,DIM);
308 
309  //Call the derivatives of the shape and test functions
310  dshape_eulerian(s,psi,dpsidx);
311 
312  //Initialise to zero
313  for(unsigned j=0;j<DIM;j++)
314  {
315  flux[j] = 0.0;
316  }
317 
318  // Loop over nodes
319  for(unsigned l=0;l<n_node;l++)
320  {
321  //Loop over derivative directions
322  for(unsigned j=0;j<DIM;j++)
323  {
324  flux[j] += this->nodal_value(l,u_nodal_index)*dpsidx(l,j);
325  }
326  }
327  }
328 
329 
330  /// Add the element's contribution to its residual vector (wrapper)
332  {
333  //Call the generic residuals function with flag set to 0
334  //using a dummy matrix argument
337  }
338 
339 
340  /// Add the element's contribution to its residual vector and
341  /// element Jacobian matrix (wrapper)
343  DenseMatrix<double> &jacobian)
344  {
345  //Call the generic routine with the flag set to 1
346  fill_in_generic_residual_contribution_poisson(residuals,jacobian,1);
347  }
348 
349 
350 
351  /// \short Return FE representation of function value u_poisson(s)
352  /// at local coordinate s
353  inline double interpolated_u_poisson(const Vector<double> &s) const
354  {
355  //Find number of nodes
356  const unsigned n_node = nnode();
357 
358  //Get the index at which the poisson unknown is stored
359  const unsigned u_nodal_index = u_index_poisson();
360 
361  //Local shape function
362  Shape psi(n_node);
363 
364  //Find values of shape function
365  shape(s,psi);
366 
367  //Initialise value of u
368  double interpolated_u = 0.0;
369 
370  //Loop over the local nodes and sum
371  for(unsigned l=0;l<n_node;l++)
372  {
373  interpolated_u += this->nodal_value(l,u_nodal_index)*psi[l];
374  }
375 
376  return(interpolated_u);
377  }
378 
379 
380  /// \short Compute derivatives of elemental residual vector with respect
381  /// to nodal coordinates. Overwrites default implementation in
382  /// FiniteElement base class.
383  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
385  dresidual_dnodal_coordinates);
386 
387  /// \short Self-test: Return 0 for OK
388  unsigned self_test();
389 
390 
391 protected:
392 
393  /// \short Shape/test functions and derivs w.r.t. to global coords at
394  /// local coord. s; return Jacobian of mapping
395  virtual double dshape_and_dtest_eulerian_poisson(const Vector<double> &s,
396  Shape &psi,
397  DShape &dpsidx, Shape &test,
398  DShape &dtestdx) const=0;
399 
400 
401  /// \short Shape/test functions and derivs w.r.t. to global coords at
402  /// integration point ipt; return Jacobian of mapping
403  virtual double dshape_and_dtest_eulerian_at_knot_poisson(const unsigned &ipt,
404  Shape &psi,
405  DShape &dpsidx,
406  Shape &test,
407  DShape &dtestdx)
408  const=0;
409 
410  /// \short Shape/test functions and derivs w.r.t. to global coords at
411  /// integration point ipt; return Jacobian of mapping (J). Also compute
412  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
414  const unsigned &ipt,
415  Shape &psi,
416  DShape &dpsidx,
417  RankFourTensor<double> &d_dpsidx_dX,
418  Shape &test,
419  DShape &dtestdx,
420  RankFourTensor<double> &d_dtestdx_dX,
421  DenseMatrix<double> &djacobian_dX) const=0;
422 
423  /// \short Compute element residual Vector only (if flag=and/or element
424  /// Jacobian matrix
426  Vector<double> &residuals, DenseMatrix<double> &jacobian,
427  const unsigned& flag);
428 
429  /// Pointer to source function:
431 
432  /// Pointer to gradient of source function
434 
435 };
436 
437 
438 
439 
440 
441 
442 ///////////////////////////////////////////////////////////////////////////
443 ///////////////////////////////////////////////////////////////////////////
444 ///////////////////////////////////////////////////////////////////////////
445 
446 
447 
448 //======================================================================
449 /// QPoissonElement elements are linear/quadrilateral/brick-shaped
450 /// Poisson elements with isoparametric interpolation for the function.
451 //======================================================================
452 template <unsigned DIM, unsigned NNODE_1D>
453  class QPoissonElement : public virtual QElement<DIM,NNODE_1D>,
454  public virtual PoissonEquations<DIM>
455 {
456 
457 private:
458 
459  /// \short Static int that holds the number of variables at
460  /// nodes: always the same
461  static const unsigned Initial_Nvalue;
462 
463  public:
464 
465 
466  ///\short Constructor: Call constructors for QElement and
467  /// Poisson equations
468  QPoissonElement() : QElement<DIM,NNODE_1D>(), PoissonEquations<DIM>()
469  {}
470 
471  /// Broken copy constructor
473  {
474  BrokenCopy::broken_copy("QPoissonElement");
475  }
476 
477  /// Broken assignment operator
479  {
480  BrokenCopy::broken_assign("QPoissonElement");
481  }
482 
483  /// \short Required # of `values' (pinned or dofs)
484  /// at node n
485  inline unsigned required_nvalue(const unsigned &n) const
486  {return Initial_Nvalue;}
487 
488  /// \short Output function:
489  /// x,y,u or x,y,z,u
490  void output(std::ostream &outfile)
492 
493 
494  /// \short Output function:
495  /// x,y,u or x,y,z,u at n_plot^DIM plot points
496  void output(std::ostream &outfile, const unsigned &n_plot)
497  {PoissonEquations<DIM>::output(outfile,n_plot);}
498 
499 
500  /// \short C-style output function:
501  /// x,y,u or x,y,z,u
502  void output(FILE* file_pt)
504 
505 
506  /// \short C-style output function:
507  /// x,y,u or x,y,z,u at n_plot^DIM plot points
508  void output(FILE* file_pt, const unsigned &n_plot)
509  {PoissonEquations<DIM>::output(file_pt,n_plot);}
510 
511 
512  /// \short Output function for an exact solution:
513  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
514  void output_fct(std::ostream &outfile, const unsigned &n_plot,
516  {PoissonEquations<DIM>::output_fct(outfile,n_plot,exact_soln_pt);}
517 
518 
519 
520  /// \short Output function for a time-dependent exact solution.
521  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
522  /// (Calls the steady version)
523  void output_fct(std::ostream &outfile, const unsigned &n_plot,
524  const double& time,
526  {PoissonEquations<DIM>::output_fct(outfile,n_plot,time,exact_soln_pt);}
527 
528 
529 protected:
530 
531 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
532  inline double dshape_and_dtest_eulerian_poisson(
533  const Vector<double> &s, Shape &psi, DShape &dpsidx,
534  Shape &test, DShape &dtestdx) const;
535 
536 
537  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
538  /// integration point ipt. Return Jacobian.
539  inline double dshape_and_dtest_eulerian_at_knot_poisson(const unsigned& ipt,
540  Shape &psi,
541  DShape &dpsidx,
542  Shape &test,
543  DShape &dtestdx)
544  const;
545 
546  /// \short Shape/test functions and derivs w.r.t. to global coords at
547  /// integration point ipt; return Jacobian of mapping (J). Also compute
548  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
550  const unsigned &ipt,
551  Shape &psi,
552  DShape &dpsidx,
553  RankFourTensor<double> &d_dpsidx_dX,
554  Shape &test,
555  DShape &dtestdx,
556  RankFourTensor<double> &d_dtestdx_dX,
557  DenseMatrix<double> &djacobian_dX) const;
558 
559 };
560 
561 
562 
563 
564 //Inline functions:
565 
566 
567 //======================================================================
568 /// Define the shape functions and test functions and derivatives
569 /// w.r.t. global coordinates and return Jacobian of mapping.
570 ///
571 /// Galerkin: Test functions = shape functions
572 //======================================================================
573 template<unsigned DIM, unsigned NNODE_1D>
575  const Vector<double> &s,
576  Shape &psi,
577  DShape &dpsidx,
578  Shape &test,
579  DShape &dtestdx) const
580 {
581  //Call the geometrical shape functions and derivatives
582  const double J = this->dshape_eulerian(s,psi,dpsidx);
583 
584  //Set the test functions equal to the shape functions
585  test = psi;
586  dtestdx= dpsidx;
587 
588  //Return the jacobian
589  return J;
590 }
591 
592 
593 
594 
595 //======================================================================
596 /// Define the shape functions and test functions and derivatives
597 /// w.r.t. global coordinates and return Jacobian of mapping.
598 ///
599 /// Galerkin: Test functions = shape functions
600 //======================================================================
601 template<unsigned DIM, unsigned NNODE_1D>
604  const unsigned &ipt,
605  Shape &psi,
606  DShape &dpsidx,
607  Shape &test,
608  DShape &dtestdx) const
609 {
610  //Call the geometrical shape functions and derivatives
611  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
612 
613  //Set the pointers of the test functions
614  test = psi;
615  dtestdx = dpsidx;
616 
617  //Return the jacobian
618  return J;
619 }
620 
621 
622 
623 //======================================================================
624 /// Define the shape functions (psi) and test functions (test) and
625 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
626 /// and return Jacobian of mapping (J). Additionally compute the
627 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
628 ///
629 /// Galerkin: Test functions = shape functions
630 //======================================================================
631 template<unsigned DIM, unsigned NNODE_1D>
634  const unsigned &ipt,
635  Shape &psi,
636  DShape &dpsidx,
637  RankFourTensor<double> &d_dpsidx_dX,
638  Shape &test,
639  DShape &dtestdx,
640  RankFourTensor<double> &d_dtestdx_dX,
641  DenseMatrix<double> &djacobian_dX) const
642  {
643  // Call the geometrical shape functions and derivatives
644  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
645  djacobian_dX,d_dpsidx_dX);
646 
647  // Set the pointers of the test functions
648  test = psi;
649  dtestdx = dpsidx;
650  d_dtestdx_dX = d_dpsidx_dX;
651 
652  //Return the jacobian
653  return J;
654 }
655 
656 
657 
658 ////////////////////////////////////////////////////////////////////////
659 ////////////////////////////////////////////////////////////////////////
660 ////////////////////////////////////////////////////////////////////////
661 
662 
663 
664 //=======================================================================
665 /// Face geometry for the QPoissonElement elements: The spatial
666 /// dimension of the face elements is one lower than that of the
667 /// bulk element but they have the same number of points
668 /// along their 1D edges.
669 //=======================================================================
670 template<unsigned DIM, unsigned NNODE_1D>
671 class FaceGeometry<QPoissonElement<DIM,NNODE_1D> >:
672  public virtual QElement<DIM-1,NNODE_1D>
673 {
674 
675  public:
676 
677  /// \short Constructor: Call the constructor for the
678  /// appropriate lower-dimensional QElement
679  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
680 
681 };
682 
683 ////////////////////////////////////////////////////////////////////////
684 ////////////////////////////////////////////////////////////////////////
685 ////////////////////////////////////////////////////////////////////////
686 
687 
688 //=======================================================================
689 /// Face geometry for the 1D QPoissonElement elements: Point elements
690 //=======================================================================
691 template<unsigned NNODE_1D>
692 class FaceGeometry<QPoissonElement<1,NNODE_1D> >:
693  public virtual PointElement
694 {
695 
696  public:
697 
698  /// \short Constructor: Call the constructor for the
699  /// appropriate lower-dimensional QElement
701 
702 };
703 
704 
705 
706 ////////////////////////////////////////////////////////////////////////
707 ////////////////////////////////////////////////////////////////////////
708 ////////////////////////////////////////////////////////////////////////
709 
710 
711 //==========================================================
712 /// Poisson upgraded to become projectable
713 //==========================================================
714  template<class POISSON_ELEMENT>
716  public virtual ProjectableElement<POISSON_ELEMENT>
717  {
718 
719  public:
720 
721  /// \short Specify the values associated with field fld.
722  /// The information is returned in a vector of pairs which comprise
723  /// the Data object and the value within it, that correspond to field fld.
725  {
726 
727 #ifdef PARANOID
728  if (fld!=0)
729  {
730  std::stringstream error_stream;
731  error_stream
732  << "Poisson elements only store a single field so fld must be 0 rather"
733  << " than " << fld << std::endl;
734  throw OomphLibError(
735  error_stream.str(),
736  OOMPH_CURRENT_FUNCTION,
737  OOMPH_EXCEPTION_LOCATION);
738  }
739 #endif
740 
741  // Create the vector
742  unsigned nnod=this->nnode();
743  Vector<std::pair<Data*,unsigned> > data_values(nnod);
744 
745  // Loop over all nodes
746  for (unsigned j=0;j<nnod;j++)
747  {
748  // Add the data value associated field: The node itself
749  data_values[j]=std::make_pair(this->node_pt(j),fld);
750  }
751 
752  // Return the vector
753  return data_values;
754  }
755 
756  /// \short Number of fields to be projected: Just one
758  {
759  return 1;
760  }
761 
762  /// \short Number of history values to be stored for fld-th field
763  /// (includes current value!)
764  unsigned nhistory_values_for_projection(const unsigned &fld)
765  {
766 #ifdef PARANOID
767  if (fld!=0)
768  {
769  std::stringstream error_stream;
770  error_stream
771  << "Poisson elements only store a single field so fld must be 0 rather"
772  << " than " << fld << std::endl;
773  throw OomphLibError(
774  error_stream.str(),
775  OOMPH_CURRENT_FUNCTION,
776  OOMPH_EXCEPTION_LOCATION);
777  }
778 #endif
779  return this->node_pt(0)->ntstorage();
780  }
781 
782  ///\short Number of positional history values
783  /// (Note: count includes current value!)
785  {
786  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
787  }
788 
789  /// \short Return Jacobian of mapping and shape functions of field fld
790  /// at local coordinate s
791  double jacobian_and_shape_of_field(const unsigned &fld,
792  const Vector<double> &s,
793  Shape &psi)
794  {
795 #ifdef PARANOID
796  if (fld!=0)
797  {
798  std::stringstream error_stream;
799  error_stream
800  << "Poisson elements only store a single field so fld must be 0 rather"
801  << " than " << fld << std::endl;
802  throw OomphLibError(
803  error_stream.str(),
804  OOMPH_CURRENT_FUNCTION,
805  OOMPH_EXCEPTION_LOCATION);
806  }
807 #endif
808  unsigned n_dim=this->dim();
809  unsigned n_node=this->nnode();
810  Shape test(n_node);
811  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
812  double J=this->dshape_and_dtest_eulerian_poisson(s,psi,dpsidx,
813  test,dtestdx);
814  return J;
815  }
816 
817 
818 
819  /// \short Return interpolated field fld at local coordinate s, at time level
820  /// t (t=0: present; t>0: history values)
821  double get_field(const unsigned &t,
822  const unsigned &fld,
823  const Vector<double>& s)
824  {
825 #ifdef PARANOID
826  if (fld!=0)
827  {
828  std::stringstream error_stream;
829  error_stream
830  << "Poisson elements only store a single field so fld must be 0 rather"
831  << " than " << fld << std::endl;
832  throw OomphLibError(
833  error_stream.str(),
834  OOMPH_CURRENT_FUNCTION,
835  OOMPH_EXCEPTION_LOCATION);
836  }
837 #endif
838  //Find the index at which the variable is stored
839  unsigned u_nodal_index = this->u_index_poisson();
840 
841  //Local shape function
842  unsigned n_node=this->nnode();
843  Shape psi(n_node);
844 
845  //Find values of shape function
846  this->shape(s,psi);
847 
848  //Initialise value of u
849  double interpolated_u = 0.0;
850 
851  //Sum over the local nodes
852  for(unsigned l=0;l<n_node;l++)
853  {
854  interpolated_u += this->nodal_value(l,u_nodal_index)*psi[l];
855  }
856  return interpolated_u;
857  }
858 
859 
860 
861 
862  ///Return number of values in field fld: One per node
863  unsigned nvalue_of_field(const unsigned &fld)
864  {
865 #ifdef PARANOID
866  if (fld!=0)
867  {
868  std::stringstream error_stream;
869  error_stream
870  << "Poisson elements only store a single field so fld must be 0 rather"
871  << " than " << fld << std::endl;
872  throw OomphLibError(
873  error_stream.str(),
874  OOMPH_CURRENT_FUNCTION,
875  OOMPH_EXCEPTION_LOCATION);
876  }
877 #endif
878  return this->nnode();
879  }
880 
881 
882  ///Return local equation number of value j in field fld.
883  int local_equation(const unsigned &fld,
884  const unsigned &j)
885  {
886 #ifdef PARANOID
887  if (fld!=0)
888  {
889  std::stringstream error_stream;
890  error_stream
891  << "Poisson elements only store a single field so fld must be 0 rather"
892  << " than " << fld << std::endl;
893  throw OomphLibError(
894  error_stream.str(),
895  OOMPH_CURRENT_FUNCTION,
896  OOMPH_EXCEPTION_LOCATION);
897  }
898 #endif
899  const unsigned u_nodal_index = this->u_index_poisson();
900  return this->nodal_local_eqn(j,u_nodal_index);
901  }
902 
903  };
904 
905 
906 //=======================================================================
907 /// Face geometry for element is the same as that for the underlying
908 /// wrapped element
909 //=======================================================================
910  template<class ELEMENT>
912  : public virtual FaceGeometry<ELEMENT>
913  {
914  public:
915  FaceGeometry() : FaceGeometry<ELEMENT>() {}
916  };
917 
918 
919 //=======================================================================
920 /// Face geometry of the Face Geometry for element is the same as
921 /// that for the underlying wrapped element
922 //=======================================================================
923  template<class ELEMENT>
925  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
926  {
927  public:
929  };
930 
931 
932 
933 
934 }
935 
936 #endif
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
PoissonEquations()
Constructor (must initialise the Source_fct_pt to null)
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual void get_source_poisson(const unsigned &ipt, const Vector< double > &x, double &source) const
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
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot ...
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
PoissonSourceFctGradientPt & source_fct_gradient_pt()
Access function: Pointer to gradient of source function.
void output(FILE *file_pt)
C_style output with default number of plot points.
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
A general Finite Element class.
Definition: elements.h:1271
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
void(* PoissonSourceFctGradientPt)(const Vector< double > &x, Vector< double > &gradient)
Function pointer to gradient of source function fct(x,g(x)) – x is a Vector!
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.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
QPoissonElement()
Constructor: Call constructors for QElement and Poisson equations.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: flux[i] = du/dx_i.
A Rank 4 Tensor class.
Definition: matrices.h:1625
double dshape_and_dtest_eulerian_poisson(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void operator=(const PoissonEquations &)
Broken assignment operator.
QPoissonElement(const QPoissonElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
void output(std::ostream &outfile)
Output with default number of plot points.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
A Rank 3 Tensor class.
Definition: matrices.h:1337
void(* PoissonSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field (includes current value!)
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
PoissonEquations(const PoissonEquations &dummy)
Broken copy constructor.
static char t char * s
Definition: cfortran.h:572
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points...
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
PoissonSourceFctGradientPt source_fct_gradient_pt() const
Access function: Pointer to gradient source function. Const version.
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 operator=(const QPoissonElement< DIM, NNODE_1D > &)
Broken assignment operator.
PoissonSourceFctGradientPt Source_fct_gradient_pt
Pointer to gradient of source function.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double interpolated_u_poisson(const Vector< double > &s) const
Return FE representation of function value u_poisson(s) at local coordinate s.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points (dummy time-dependent versi...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
PoissonSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
unsigned nfields_for_projection()
Number of fields to be projected: Just one.
virtual unsigned u_index_poisson() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
PoissonSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
virtual void get_source_gradient_poisson(const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient) const
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
Poisson upgraded to become projectable.
virtual double dshape_and_dtest_eulerian_poisson(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
virtual double dshape_and_dtest_eulerian_at_knot_poisson(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
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.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
PoissonSourceFctPt Source_fct_pt
Pointer to source function:
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
double dshape_and_dtest_eulerian_at_knot_poisson(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
unsigned self_test()
Self-test: Return 0 for OK.