unsteady_heat_flux_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Header file for elements that are used to apply prescribed flux
31 // boundary conditions to the UnsteadyHeat equations
32 
33 
34 #ifndef OOMPH_UNSTEADY_HEAT_FLUX_ELEMENTS_HEADER
35 #define OOMPH_UNSTEADY_HEAT_FLUX_ELEMENTS_HEADER
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39  #include <oomph-lib-config.h>
40 #endif
41 
42 //Standard libray headers
43 #include <cmath>
44 
45 // oomph-lib ncludes
46 #include "../generic/Qelements.h"
47 
48 
49 namespace oomph
50 {
51 
52 //======================================================================
53 /// \short A class for elements that allow the imposition of an
54 /// applied flux on the boundaries of UnsteadyHeat elements.
55 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
56 /// policy class.
57 //======================================================================
58 template <class ELEMENT>
59 class UnsteadyHeatFluxElement : public virtual FaceGeometry<ELEMENT>,
60 public virtual FaceElement
61 {
62 
63 public:
64 
65 
66  /// \short Function pointer to the prescribed-flux function fct(x,f(x)) --
67  /// x is a Vector!
68  typedef void (*UnsteadyHeatPrescribedFluxFctPt)
69  (const double& time, const Vector<double>& x, double& flux);
70 
71 
72  /// \short Constructor, takes the pointer to the "bulk" element and the
73  /// index of the face to be created
74  UnsteadyHeatFluxElement(FiniteElement* const &bulk_el_pt,
75  const int &face_index);
76 
77  /// Broken copy constructor
79  {
80  BrokenCopy::broken_copy("UnsteadyHeatFluxElement");
81  }
82 
83  /// Broken assignment operator
84 //Commented out broken assignment operator because this can lead to a conflict warning
85 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
86 //realise that two separate implementations of the broken function are the same and so,
87 //quite rightly, it shouts.
88  /*void operator=(const UnsteadyHeatFluxElement&)
89  {
90  BrokenCopy::broken_assign("UnsteadyHeatFluxElement");
91  }*/
92 
93  /// Access function for the prescribed-flux function pointer
95 
96  /// Compute the element residual vector
98  {
99  //Call the generic residuals function with flag set to 0
100  //using a dummy matrix argument
103  }
104 
105 
106  /// Compute the element's residual vector and its Jacobian matrix
108  DenseMatrix<double> &jacobian)
109  {
110  //Call the generic routine with the flag set to 1
112  }
113 
114 
115  /// \short Specify the value of nodal zeta from the face geometry:
116  /// The "global" intrinsic coordinate of the element when
117  /// viewed as part of a geometric object should be given by
118  /// the FaceElement representation, by default (needed to break
119  /// indeterminacy if bulk element is SolidElement)
120  double zeta_nodal(const unsigned &n, const unsigned &k,
121  const unsigned &i) const
122  {return FaceElement::zeta_nodal(n,k,i);}
123 
124  /// Output function -- forward to broken version in FiniteElement
125  /// until somebody decides what exactly they want to plot here...
126  void output(std::ostream &outfile) {FaceGeometry<ELEMENT>::output(outfile);}
127 
128  /// \short Output function -- forward to broken version in FiniteElement
129  /// until somebody decides what exactly they want to plot here...
130  void output(std::ostream &outfile, const unsigned &n_plot)
131  {FaceGeometry<ELEMENT>::output(outfile,n_plot);}
132 
133  /// C-style output function -- forward to broken version in FiniteElement
134  /// until somebody decides what exactly they want to plot here...
135  void output(FILE* file_pt) {FaceGeometry<ELEMENT>::output(file_pt);}
136 
137  /// \short C-style output function -- forward to broken version in
138  /// FiniteElement until somebody decides what exactly they want to plot
139  /// here...
140  void output(FILE* file_pt, const unsigned &n_plot)
141  {FaceGeometry<ELEMENT>::output(file_pt,n_plot);}
142 
143 protected:
144 
145  /// \short Function to compute the shape and test functions and to return
146  /// the Jacobian of mapping between local and global (Eulerian)
147  /// coordinates
148  inline double shape_and_test(const Vector<double> &s, Shape &psi, Shape &test)
149  const
150  {
151  //Find number of nodes
152  unsigned n_node = nnode();
153 
154  //Get the shape functions
155  shape(s,psi);
156 
157  //Set the test functions to be the same as the shape functions
158  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
159 
160  //Return the value of the jacobian
161  return J_eulerian(s);
162  }
163 
164  /// Function to calculate the prescribed flux at a given spatial
165  /// position and at a gien time
166  void get_flux(const double& time, const Vector<double>& x, double& flux)
167  {
168  //If the function pointer is zero return zero
169  if(Flux_fct_pt == 0)
170  {
171  flux=0.0;
172  }
173  //Otherwise call the function
174  else
175  {
176  (*Flux_fct_pt)(time,x,flux);
177  }
178  }
179 
180 
181 private:
182 
183 
184  /// \short Compute the element residual vector.
185  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
187  Vector<double> &residuals, DenseMatrix<double> &jacobian,
188  unsigned flag);
189 
190 
191  /// Function pointer to the (global) prescribed-flux function
193 
194  ///The spatial dimension of the problem
195  unsigned Dim;
196 
197  ///The index at which the unknown is stored at the nodes
199 
200 
201 };
202 
203 ////////////////////////////////////////////////////////////////////////
204 ////////////////////////////////////////////////////////////////////////
205 ////////////////////////////////////////////////////////////////////////
206 
207 //===========================================================================
208 /// Constructor, takes the pointer to the "bulk" element and the
209 /// index of the face to be created.
210 //===========================================================================
211 template<class ELEMENT>
214  const int &face_index) :
215  FaceGeometry<ELEMENT>(), FaceElement()
216 {
217 
218  // Let the bulk element build the FaceElement, i.e. setup the pointers
219  // to its nodes (by referring to the appropriate nodes in the bulk
220  // element), etc.
221  bulk_el_pt->build_face_element(face_index,this);
222 
223 #ifdef PARANOID
224  {
225  //Check that the element is not a refineable 3d element
226  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
227 
228  //If it's three-d
229  if(elem_pt->dim()==3)
230  {
231  //Is it refineable
232  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
233  if(ref_el_pt!=0)
234  {
235  if (this->has_hanging_nodes())
236  {
237  throw OomphLibError(
238  "This flux element will not work correctly if nodes are hanging\n",
239  OOMPH_CURRENT_FUNCTION,
240  OOMPH_EXCEPTION_LOCATION);
241  }
242  }
243  }
244  }
245 #endif
246 
247  // Initialise the prescribed-flux function pointer to zero
248  Flux_fct_pt = 0;
249 
250  // Extract the dimension of the problem from the dimension of
251  // the first node
252  Dim = this->node_pt(0)->ndim();
253 
254  //Set up U_index_ust_heat. Initialise to zero, which probably won't change
255  //in most cases, oh well, the price we pay for generality
256  U_index_ust_heat = 0;
257 
258  //Cast to the appropriate UnsteadyHeatEquation so that we can
259  //find the index at which the variable is stored
260  //We assume that the dimension of the full problem is the same
261  //as the dimension of the node, if this is not the case you will have
262  //to write custom elements, sorry
263  switch(Dim)
264  {
265  //One dimensional problem
266  case 1:
267  {
268  UnsteadyHeatEquations<1>* eqn_pt =
269  dynamic_cast<UnsteadyHeatEquations<1>*>(bulk_el_pt);
270  //If the cast has failed die
271  if(eqn_pt==0)
272  {
273  std::string error_string =
274  "Bulk element must inherit from UnsteadyHeatEquations.";
275  error_string +=
276  "Nodes are one dimensional, but cannot cast the bulk element to\n";
277  error_string += "UnsteadyHeatEquations<1>\n.";
278  error_string +=
279  "If you desire this functionality, you must implement it yourself\n";
280 
281  throw OomphLibError(error_string,
282  OOMPH_CURRENT_FUNCTION,
283  OOMPH_EXCEPTION_LOCATION);
284  }
285  //Otherwise read out the value
286  else
287  {
288  //Read the index from the (cast) bulk element
290  }
291  }
292  break;
293 
294  //Two dimensional problem
295  case 2:
296  {
297  UnsteadyHeatEquations<2>* eqn_pt =
298  dynamic_cast<UnsteadyHeatEquations<2>*>(bulk_el_pt);
299  //If the cast has failed die
300  if(eqn_pt==0)
301  {
302  std::string error_string =
303  "Bulk element must inherit from UnsteadyHeatEquations.";
304  error_string +=
305  "Nodes are two dimensional, but cannot cast the bulk element to\n";
306  error_string += "UnsteadyHeatEquations<2>\n.";
307  error_string +=
308  "If you desire this functionality, you must implement it yourself\n";
309 
310  throw OomphLibError(error_string,
311  OOMPH_CURRENT_FUNCTION,
312  OOMPH_EXCEPTION_LOCATION);
313  }
314  else
315  {
316  //Read the index from the (cast) bulk element.
318  }
319  }
320  break;
321 
322  //Three dimensional problem
323  case 3:
324  {
325  UnsteadyHeatEquations<3>* eqn_pt =
326  dynamic_cast<UnsteadyHeatEquations<3>*>(bulk_el_pt);
327  //If the cast has failed die
328  if(eqn_pt==0)
329  {
330  std::string error_string =
331  "Bulk element must inherit from UnsteadyHeatEquations.";
332  error_string +=
333  "Nodes are three dimensional, but cannot cast the bulk element to\n";
334  error_string += "UnsteadyHeatEquations<3>\n.";
335  error_string +=
336  "If you desire this functionality, you must implement it yourself\n";
337 
338  throw OomphLibError(error_string,
339  OOMPH_CURRENT_FUNCTION,
340  OOMPH_EXCEPTION_LOCATION);
341 
342  }
343  else
344  {
345  //Read the index from the (cast) bulk element.
347  }
348  }
349  break;
350 
351  //Any other case is an error
352  default:
353  std::ostringstream error_stream;
354  error_stream << "Dimension of node is " << Dim
355  << ". It should be 1,2, or 3!" << std::endl;
356 
357  throw OomphLibError(error_stream.str(),
358  OOMPH_CURRENT_FUNCTION,
359  OOMPH_EXCEPTION_LOCATION);
360  break;
361  }
362 
363 }
364 
365 
366 
367 //===========================================================================
368 /// Compute the element's residual vector and the (zero) Jacobian matrix.
369 //===========================================================================
370 template<class ELEMENT>
373  Vector<double> &residuals, DenseMatrix<double> &jacobian,
374  unsigned flag)
375 {
376  //Find out how many nodes there are
377  const unsigned n_node = nnode();
378 
379  // Get continuous time from timestepper of first node
380  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
381 
382  //Set up memory for the shape and test functions
383  Shape psif(n_node), testf(n_node);
384 
385  //Set the value of n_intpt
386  const unsigned n_intpt = integral_pt()->nweight();
387 
388  //Set the Vector to hold local coordinates
389  Vector<double> s(Dim-1);
390 
391  //Integer to store the local equation and unknowns
392  int local_eqn=0;
393 
394  // Locally cache the index at which the variable is stored
395  const unsigned u_index_ust_heat = U_index_ust_heat;
396 
397  //Loop over the integration points
398  //--------------------------------
399  for(unsigned ipt=0;ipt<n_intpt;ipt++)
400  {
401 
402  //Assign values of s
403  for(unsigned i=0;i<(Dim-1);i++)
404  {
405  s[i] = integral_pt()->knot(ipt,i);
406  }
407 
408  //Get the integral weight
409  double w = integral_pt()->weight(ipt);
410 
411  //Find the shape and test functions and return the Jacobian
412  //of the mapping
413  double J = shape_and_test(s,psif,testf);
414 
415  //Premultiply the weights and the Jacobian
416  double W = w*J;
417 
418  //Need to find position to feed into flux function
419  Vector<double> interpolated_x(Dim);
420 
421  //Initialise to zero
422  for(unsigned i=0;i<Dim;i++)
423  {
424  interpolated_x[i] = 0.0;
425  }
426 
427  //Calculate velocities and derivatives
428  for(unsigned l=0;l<n_node;l++)
429  {
430  //Loop over velocity components
431  for(unsigned i=0;i<Dim;i++)
432  {
433  interpolated_x[i] += nodal_position(l,i)*psif[l];
434  }
435  }
436 
437  //Get the imposed flux
438  double flux;
439  get_flux(time,interpolated_x,flux);
440 
441  //Now add to the appropriate equations
442 
443  //Loop over the test functions
444  for(unsigned l=0;l<n_node;l++)
445  {
446  local_eqn = nodal_local_eqn(l,u_index_ust_heat);
447  /*IF it's not a boundary condition*/
448  if(local_eqn >= 0)
449  {
450  //Add the prescribed flux terms
451  residuals[local_eqn] -= flux*testf[l]*W;
452 
453  // Imposed traction doesn't depend upon the solution,
454  // --> the Jacobian is always zero, so no Jacobian
455  // terms are required
456  }
457  }
458  }
459 }
460 
461 
462 
463 
464 
465 }
466 
467 #endif
virtual unsigned u_index_ust_heat() const
Broken assignment operator.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element residual vector.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void fill_in_generic_residual_contribution_ust_heat_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute the element residual vector. flag=1(or 0): do (or don't) compute the Jacobian as well...
void get_flux(const double &time, const Vector< double > &x, double &flux)
unsigned Dim
The spatial dimension of the problem.
unsigned U_index_ust_heat
The index at which the unknown is stored at the nodes.
cstr elem_len * i
Definition: cfortran.h:607
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
UnsteadyHeatPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
A general Finite Element class.
Definition: elements.h:1271
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
UnsteadyHeatFluxElement(const UnsteadyHeatFluxElement &dummy)
Broken copy constructor.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
UnsteadyHeatFluxElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Constructor, takes the pointer to the "bulk" element and the index of the face to be created...
static char t char * s
Definition: cfortran.h:572
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:5033
void output(std::ostream &outfile)
Output with default number of plot points.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
A class for elements that allow the imposition of an applied flux on the boundaries of UnsteadyHeat e...
void(* UnsteadyHeatPrescribedFluxFctPt)(const double &time, const Vector< double > &x, double &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector!
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry: The "global" intrinsic coordinate of the elem...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
UnsteadyHeatPrescribedFluxFctPt & flux_fct_pt()
Broken assignment operator.
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual vector and its Jacobian matrix.
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...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228