refineable_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//====================================================================
27
28
29namespace oomph
30{
31 //========================================================================
32 /// Validate against exact flux.
33 /// Flux is provided via function pointer.
34 /// Plot error at a given number of plot points.
35 //========================================================================
36 template<unsigned DIM>
38 std::ostream& outfile,
40 double& error,
41 double& norm)
42 {
43 // Initialise
44 error = 0.0;
45 norm = 0.0;
46
47 // Vector of local coordinates
48 Vector<double> s(DIM);
49
50 // Vector of global coordinates
51 Vector<double> x(DIM);
52
53 // Find out how many nodes there are in the element
54 const unsigned n_node = nnode();
55
56 // Allocate storage for shape functions
57 Shape psi(n_node);
58
59 // Set the value of n_intpt
60 const unsigned n_intpt = integral_pt()->nweight();
61
62 // Tecplot
63 outfile << "ZONE" << std::endl;
64
65 // Exact flux Vector (size is DIM)
66 Vector<double> exact_flux(DIM);
67
68 // Loop over the integration points
69 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
70 {
71 // Assign values of s
72 for (unsigned i = 0; i < DIM; i++)
73 {
74 s[i] = integral_pt()->knot(ipt, i);
75 }
76
77 // Get the integral weight
78 double w = integral_pt()->weight(ipt);
79
80 // Get jacobian of mapping
81 double J = J_eulerian(s);
82
83 // Premultiply the weights and the Jacobian
84 double W = w * J;
85
86 // Get x position as Vector
87 interpolated_x(s, x);
88
89 // Allocate storage for FE flux Vector
90 Vector<double> fe_flux(DIM);
91
92 // Get FE flux as Vector
93 get_Z2_flux(s, fe_flux);
94
95 // Get exact flux at this point
96 (*exact_flux_pt)(x, exact_flux);
97
98 // Output x,y,...
99 for (unsigned i = 0; i < DIM; i++)
100 {
101 outfile << x[i] << " ";
102 }
103
104 // Output exact flux
105 for (unsigned i = 0; i < DIM; i++)
106 {
107 outfile << exact_flux[i] << " ";
108 }
109
110 // Output error
111 for (unsigned i = 0; i < DIM; i++)
112 {
113 outfile << exact_flux[i] - fe_flux[i] << " ";
114 }
115 outfile << std::endl;
116
117 // Add to RMS error:
118 double sum = 0.0;
119 double sum2 = 0.0;
120 for (unsigned i = 0; i < DIM; i++)
121 {
122 sum += (fe_flux[i] - exact_flux[i]) * (fe_flux[i] - exact_flux[i]);
123 sum2 += exact_flux[i] * exact_flux[i];
124 }
125 error += sum * W;
126 norm += sum2 * W;
127
128 } // End of loop over the integration points
129 }
130
131
132 //========================================================================
133 /// Add element's contribution to the elemental
134 /// residual vector and/or Jacobian matrix.
135 /// flag=1: compute both
136 /// flag=0: compute only residual vector
137 //========================================================================
138 template<unsigned DIM>
141 DenseMatrix<double>& jacobian,
142 const unsigned& flag)
143 {
144 // Find out how many nodes there are in the element
145 unsigned n_node = nnode();
146
147 // Set up memory for the shape and test functions
148 Shape psi(n_node), test(n_node);
149 DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
150
151 // Set the value of n_intpt
152 unsigned n_intpt = integral_pt()->nweight();
153
154 // The local index at which the poisson variable is stored
155 unsigned u_nodal_index = this->u_index_poisson();
156
157 // Integers to store the local equation and unknown numbers
158 int local_eqn = 0, local_unknown = 0;
159
160 // Local storage for pointers to hang_info objects
161 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
162
163 // Loop over the integration points
164 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
165 {
166 // Get the integral weight
167 double w = integral_pt()->weight(ipt);
168
169 // Call the derivatives of the shape and test functions
170 double J = this->dshape_and_dtest_eulerian_at_knot_poisson(
171 ipt, psi, dpsidx, test, dtestdx);
172
173 // Premultiply the weights and the Jacobian
174 double W = w * J;
175
176 // Position and gradient
177 Vector<double> interpolated_x(DIM, 0.0);
178 Vector<double> interpolated_dudx(DIM, 0.0);
179
180 // Calculate function value and derivatives:
181 //-----------------------------------------
182
183 // Loop over nodes
184 for (unsigned l = 0; l < n_node; l++)
185 {
186 // Get the poisson value from the node
187 //(hanging-ness will be taken into account
188 double u_value = this->nodal_value(l, u_nodal_index);
189
190 // Loop over directions
191 for (unsigned j = 0; j < DIM; j++)
192 {
193 interpolated_x[j] += nodal_position(l, j) * psi(l);
194 interpolated_dudx[j] += u_value * dpsidx(l, j);
195 }
196 }
197
198 // Get body force
199 double source;
200 this->get_source_poisson(ipt, interpolated_x, source);
201
202
203 // Assemble residuals and Jacobian
204
205 // Loop over the nodes for the test functions
206 for (unsigned l = 0; l < n_node; l++)
207 {
208 // Local variables used to store the number of master nodes and the
209 // weight associated with the shape function if the node is hanging
210 unsigned n_master = 1;
211 double hang_weight = 1.0;
212
213 // Local bool (is the node hanging)
214 bool is_node_hanging = this->node_pt(l)->is_hanging();
215
216 // If the node is hanging, get the number of master nodes
217 if (is_node_hanging)
218 {
219 hang_info_pt = this->node_pt(l)->hanging_pt();
220 n_master = hang_info_pt->nmaster();
221 }
222 // Otherwise there is just one master node, the node itself
223 else
224 {
225 n_master = 1;
226 }
227
228 // Loop over the master nodes
229 for (unsigned m = 0; m < n_master; m++)
230 {
231 // Get the local equation number and hang_weight
232 // If the node is hanging
233 if (is_node_hanging)
234 {
235 // Read out the local equation number from the m-th master node
236 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
237 u_nodal_index);
238
239 // Read out the weight from the master node
240 hang_weight = hang_info_pt->master_weight(m);
241 }
242 // If the node is not hanging
243 else
244 {
245 // The local equation number comes from the node itself
246 local_eqn = this->nodal_local_eqn(l, u_nodal_index);
247
248 // The hang weight is one
249 hang_weight = 1.0;
250 }
251
252 // If the nodal equation is not a boundary condition
253 if (local_eqn >= 0)
254 {
255 // Add body force/source term here
256 residuals[local_eqn] += source * test(l) * W * hang_weight;
257
258 // The Poisson bit itself
259 for (unsigned k = 0; k < DIM; k++)
260 {
261 residuals[local_eqn] +=
262 interpolated_dudx[k] * dtestdx(l, k) * W * hang_weight;
263 }
264
265 // Calculate the Jacobian
266 if (flag)
267 {
268 // Local variables to store the number of master nodes
269 // and the weights associated with each hanging node
270 unsigned n_master2 = 1;
271 double hang_weight2 = 1.0;
272
273 // Loop over the nodes for the variables
274 for (unsigned l2 = 0; l2 < n_node; l2++)
275 {
276 // Local bool (is the node hanging)
277 bool is_node2_hanging = this->node_pt(l2)->is_hanging();
278
279 // If the node is hanging, get the number of master nodes
280 if (is_node2_hanging)
281 {
282 hang_info2_pt = this->node_pt(l2)->hanging_pt();
283 n_master2 = hang_info2_pt->nmaster();
284 }
285 // Otherwise there is one master node, the node itself
286 else
287 {
288 n_master2 = 1;
289 }
290
291 // Loop over the master nodes
292 for (unsigned m2 = 0; m2 < n_master2; m2++)
293 {
294 // Get the local unknown and weight
295 // If the node is hanging
296 if (is_node2_hanging)
297 {
298 // Read out the local unknown from the master node
299 local_unknown = this->local_hang_eqn(
300 hang_info2_pt->master_node_pt(m2), u_nodal_index);
301
302 // Read out the hanging weight from the master node
303 hang_weight2 = hang_info2_pt->master_weight(m2);
304 }
305 // If the node is not hanging
306 else
307 {
308 // The local unknown number comes from the node
309 local_unknown = this->nodal_local_eqn(l2, u_nodal_index);
310
311 // The hang weight is one
312 hang_weight2 = 1.0;
313 }
314
315 // If the unknown is not pinned
316 if (local_unknown >= 0)
317 {
318 // Add contribution to Elemental Matrix
319 for (unsigned i = 0; i < DIM; i++)
320 {
321 jacobian(local_eqn, local_unknown) +=
322 dpsidx(l2, i) * dtestdx(l, i) * W * hang_weight *
323 hang_weight2;
324 }
325 }
326 } // End of loop over master nodes
327 } // End of loop over nodes
328 } // End of Jacobian calculation
329
330 } // End of case when residual equation is not pinned
331 } // End of loop over master nodes for residual
332 } // End of loop over nodes
333
334 } // End of loop over integration points
335 }
336
337
338 //======================================================================
339 /// Compute derivatives of elemental residual vector with respect
340 /// to nodal coordinates (fully analytical).
341 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
342 /// Overloads the FD-based version in the FE base class.
343 //======================================================================
344 template<unsigned DIM>
346 RankThreeTensor<double>& dresidual_dnodal_coordinates)
347 {
348 // Determine number of nodes in element
349 const unsigned n_node = nnode();
350
351 // Set up memory for the shape and test functions
352 Shape psi(n_node), test(n_node);
353 DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
354
355 // Get number of shape controlling nodes
356 const unsigned n_shape_controlling_node = nshape_controlling_nodes();
357
358 // Deriatives of shape fct derivatives w.r.t. nodal coords
359 RankFourTensor<double> d_dpsidx_dX(
360 DIM, n_shape_controlling_node, n_node, DIM);
361 RankFourTensor<double> d_dtestdx_dX(
362 DIM, n_shape_controlling_node, n_node, DIM);
363
364 // Derivative of Jacobian of mapping w.r.t. to nodal coords
365 DenseMatrix<double> dJ_dX(DIM, n_shape_controlling_node);
366
367 // Derivatives of derivative of u w.r.t. nodal coords
368 RankThreeTensor<double> d_dudx_dX(DIM, n_shape_controlling_node, DIM);
369
370 // Source function and its gradient
371 double source;
372 Vector<double> d_source_dx(DIM);
373
374 // Index at which the poisson unknown is stored
375 const unsigned u_nodal_index = this->u_index_poisson();
376
377 // Determine the number of integration points
378 const unsigned n_intpt = integral_pt()->nweight();
379
380 // Integer to store the local equation number
381 int local_eqn = 0;
382
383 // Local storage for pointers to hang_info object
384 HangInfo* hang_info_pt = 0;
385
386 // Loop over the integration points
387 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
388 {
389 // Get the integral weight
390 double w = integral_pt()->weight(ipt);
391
392 // Call the derivatives of the shape/test functions, as well as the
393 // derivatives of these w.r.t. nodal coordinates and the derivative
394 // of the jacobian of the mapping w.r.t. nodal coordinates
395 const double J = this->dshape_and_dtest_eulerian_at_knot_poisson(
396 ipt, psi, dpsidx, d_dpsidx_dX, test, dtestdx, d_dtestdx_dX, dJ_dX);
397
398 // Calculate local values
399 // Allocate and initialise to zero
400 Vector<double> interpolated_x(DIM, 0.0);
401 Vector<double> interpolated_dudx(DIM, 0.0);
402
403 // Calculate function value and derivatives:
404 // -----------------------------------------
405 // Loop over nodes
406 for (unsigned l = 0; l < n_node; l++)
407 {
408 // Get the nodal value of the Poisson unknown
409 double u_value = nodal_value(l, u_nodal_index);
410
411 // Loop over directions
412 for (unsigned i = 0; i < DIM; i++)
413 {
414 interpolated_x[i] += nodal_position(l, i) * psi(l);
415 interpolated_dudx[i] += u_value * dpsidx(l, i);
416 }
417 }
418
419 // Calculate derivative of du/dx_i w.r.t. nodal positions X_{pq}
420
421 // Loop over shape controlling nodes
422 for (unsigned q = 0; q < n_shape_controlling_node; q++)
423 {
424 // Loop over coordinate directions
425 for (unsigned p = 0; p < DIM; p++)
426 {
427 for (unsigned i = 0; i < DIM; i++)
428 {
429 double aux = 0.0;
430 for (unsigned j = 0; j < n_node; j++)
431 {
432 aux += nodal_value(j, u_nodal_index) * d_dpsidx_dX(p, q, j, i);
433 }
434 d_dudx_dX(p, q, i) = aux;
435 }
436 }
437 }
438
439 // Get source function
440 this->get_source_poisson(ipt, interpolated_x, source);
441
442 // Get gradient of source function
443 this->get_source_gradient_poisson(ipt, interpolated_x, d_source_dx);
444
445 // std::map<Node*,unsigned>
446 // local_shape_controlling_node_lookup=shape_controlling_node_lookup();
447
448 // Assemble d res_{local_eqn} / d X_{pq}
449 // -------------------------------------
450
451 // Loop over the nodes for the test functions
452 for (unsigned l = 0; l < n_node; l++)
453 {
454 // Local variables used to store the number of master nodes and the
455 // weight associated with the shape function if the node is hanging
456 unsigned n_master = 1;
457 double hang_weight = 1.0;
458
459 // Local bool (is the node hanging)
460 bool is_node_hanging = this->node_pt(l)->is_hanging();
461
462 // If the node is hanging, get the number of master nodes
463 if (is_node_hanging)
464 {
465 hang_info_pt = this->node_pt(l)->hanging_pt();
466 n_master = hang_info_pt->nmaster();
467 }
468 // Otherwise there is just one master node, the node itself
469 else
470 {
471 n_master = 1;
472 }
473
474 // Loop over the master nodes
475 for (unsigned m = 0; m < n_master; m++)
476 {
477 // Get the local equation number and hang_weight
478 // If the node is hanging
479 if (is_node_hanging)
480 {
481 // Read out the local equation number from the m-th master node
482 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
483 u_nodal_index);
484
485 // Read out the weight from the master node
486 hang_weight = hang_info_pt->master_weight(m);
487 }
488 // If the node is not hanging
489 else
490 {
491 // The local equation number comes from the node itself
492 local_eqn = this->nodal_local_eqn(l, u_nodal_index);
493 // The hang weight is one
494 hang_weight = 1.0;
495 }
496
497 // If the nodal equation is not a boundary condition
498 if (local_eqn >= 0)
499 {
500 // Loop over coordinate directions
501 for (unsigned p = 0; p < DIM; p++)
502 {
503 // Loop over shape controlling nodes
504 for (unsigned q = 0; q < n_shape_controlling_node; q++)
505 {
506 double sum = source * test(l) * dJ_dX(p, q) +
507 d_source_dx[p] * test(l) * psi(q) * J;
508
509 for (unsigned i = 0; i < DIM; i++)
510 {
511 sum += interpolated_dudx[i] * (dtestdx(l, i) * dJ_dX(p, q) +
512 d_dtestdx_dX(p, q, l, i) * J) +
513 d_dudx_dX(p, q, i) * dtestdx(l, i) * J;
514 }
515
516 // Multiply through by integration weight
517 dresidual_dnodal_coordinates(local_eqn, p, q) +=
518 sum * w * hang_weight;
519 }
520 }
521 }
522 }
523 }
524 } // End of loop over integration points
525 }
526
527 /// Get error against and norm of exact solution
528 template<unsigned DIM>
530 std::ostream& outfile,
532 double& error,
533 double& norm)
534 {
535 // Initialise
536 error = 0.0;
537 norm = 0.0;
538
539 // Vector of local coordinates
540 Vector<double> s(DIM);
541
542 // Vector for coordintes
543 Vector<double> x(DIM);
544
545 // Set the value of n_intpt
546 unsigned n_intpt = this->integral_pt()->nweight();
547
548 // Setup output structure: Conversion is fishy but it's only output...
549 unsigned nplot;
550 if (DIM == 1)
551 {
552 nplot = n_intpt;
553 }
554 else
555 {
556 nplot = unsigned(pow(n_intpt, 1.0 / double(DIM)));
557 }
558
559 // Tecplot header info
560 outfile << this->tecplot_zone_string(nplot);
561
562 // Exact gradient Vector
563 Vector<double> exact_grad(DIM);
564
565 // Loop over the integration points
566 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
567 {
568 // Assign values of s
569 for (unsigned i = 0; i < DIM; i++)
570 {
571 s[i] = this->integral_pt()->knot(ipt, i);
572 }
573
574 // Get the integral weight
575 double w = this->integral_pt()->weight(ipt);
576
577 // Get jacobian of mapping
578 double J = this->J_eulerian(s);
579
580 // Premultiply the weights and the Jacobian
581 double W = w * J;
582
583 // Get x position as Vector
584 this->interpolated_x(s, x);
585
586 // Get FE du/dx
587 Vector<double> dudx_fe(DIM);
589
590 // Get exact gradient at this point
591 (*exact_grad_pt)(x, exact_grad);
592
593 // Output x,y,...,error
594 for (unsigned i = 0; i < DIM; i++)
595 {
596 outfile << x[i] << " ";
597 }
598 for (unsigned i = 0; i < DIM; i++)
599 {
600 outfile << exact_grad[i] << " " << exact_grad[i] - dudx_fe[i]
601 << std::endl;
602 }
603
604 // Add to error and norm
605 for (unsigned i = 0; i < DIM; i++)
606 {
607 norm += exact_grad[i] * exact_grad[i] * W;
608 error +=
609 (exact_grad[i] - dudx_fe[i]) * (exact_grad[i] - dudx_fe[i]) * W;
610 }
611 }
612 }
613
614 template<unsigned DIM>
616 {
617 if (this->tree_pt()->father_pt() != 0)
618 {
619 // Needed to set the source function pointer (if there is a father)
621 }
622 // Now do the PRefineableQElement further_build()
624 }
625
626
627 //====================================================================
628 // Force build of templates
629 //====================================================================
630 template class RefineableQPoissonElement<1, 2>;
631 template class RefineableQPoissonElement<1, 3>;
632 template class RefineableQPoissonElement<1, 4>;
633
634 template class RefineableQPoissonElement<2, 2>;
635 template class RefineableQPoissonElement<2, 3>;
636 template class RefineableQPoissonElement<2, 4>;
637
638 template class RefineableQPoissonElement<3, 2>;
639 template class RefineableQPoissonElement<3, 3>;
640 template class RefineableQPoissonElement<3, 4>;
641
642 template class PRefineableQPoissonElement<1>;
643 template class PRefineableQPoissonElement<2>;
644 template class PRefineableQPoissonElement<3>;
645
646} // 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
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
Definition: Qelements.h:2274
p-refineable version of 2D QPoissonElement elements
void further_build()
Further build: Copy source function pointer from father element.
void compute_energy_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_grad_pt, double &error, double &norm)
Get error against and norm of exact solution.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: flux[i] = du/dx_i.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1701
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
void fill_in_generic_residual_contribution_poisson(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
void compute_exact_Z2_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
Get error against and norm of exact flux.
void further_build()
Further build: Copy source function pointer from father element.
Refineable version of 2D QPoissonElement elements.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...