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