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