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