biharmonic_problem.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: 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//====================================================================
30 #ifndef OOMPH_BIHARMONIC_PROBLEM_HEADER
31 #define OOMPH_BIHARMONIC_PROBLEM_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35  #include <oomph-lib-config.h>
36 #endif
37 
38 #ifdef OOMPH_HAS_MPI
39 //mpi headers
40 #include "mpi.h"
41 #endif
42 
43 // Generic C++ headers
44 #include <iostream>
45 #include <math.h>
46 
47 // oomph-lib headers
48 #include "../generic/problem.h"
49 #include "../generic/hijacked_elements.h"
50 #include "../meshes/hermite_element_quad_mesh.template.h"
51 #include "../meshes/hermite_element_quad_mesh.template.cc"
52 #include "biharmonic_elements.h"
54 
55 
56 namespace oomph
57 {
58 
59 
60 //=============================================================================
61 /// \short Biharmonic Plate Problem Class - for problems where the load can be
62 /// assumed to be acting normal to the surface of the plate and the deflections
63 /// are small relative to the thickness of the plate. Developed for the
64 /// topologically rectangular Hermite Element Mesh. Contains functions
65 /// allowing the following boundary conditions to be applied (on a given edge):
66 /// + clamped : u and du/dn imposed
67 /// + simply supported : u and laplacian(u) imposed
68 /// + free : laplacian(u) and dlaplacian(u)/dn imposed
69 //=============================================================================
70 template<unsigned DIM> class BiharmonicProblem : public Problem
71 {
72 
73 public:
74 
75 
76  /// \short Definition of a dirichlet boundary condition function pointer.
77  /// Takes the position along a boundary (s) in the macro element coordinate
78  /// scheme and returns the value of the boundary condition at that point (u).
79  typedef void (*DirichletBCFctPt)(const double& s, double& u);
80 
81 
82  /// \short Definition of the Source Function.
83  typedef void (*BiharmonicSourceFctPt)(const Vector<double>& x, double& f);
84 
85  /// \short Constructor
87  {
90  }
91 
92  /// Destructor. Delete the meshes
94  {
95  delete Bulk_element_mesh_pt;
96  delete Face_element_mesh_pt;
97  };
98 
99  /// actions before solve, performs self test
101  {
102 #ifdef PARANOID
103  if (0 == self_test())
104  {
105  oomph_info << "self test passed" << std::endl;
106  }
107  else
108  {
109  oomph_info << "self test failed" << std::endl;
110  }
111 #endif
112  }
113 
114  /// action after solve
116  {}
117 
118  /// \short documents the solution, and if an exact solution is provided, then
119  /// the error between the numerical and exact solution is presented
120  void doc_solution(DocInfo& doc_info,
121  FiniteElement::SteadyExactSolutionFctPt exact_soln_pt = 0);
122 
123  /// \short Access function to the bulk element mesh pt
125  {
126  return Bulk_element_mesh_pt;
127  }
128 
129 protected:
130 
131  /// \short builds the bulk mesh on a prescribed domain with a node spacing
132  /// defined by spacing fn with n_x by n_y elements
133  void build_bulk_mesh(const unsigned n_x, const unsigned n_y,
135  HermiteQuadMesh<BiharmonicElement<2> >::MeshSpacingFnPtr
136  spacing_fn = 0)
137  {
138  if (spacing_fn == 0)
139  {
140  Bulk_element_mesh_pt = new
141  HermiteQuadMesh<Hijacked<BiharmonicElement<2> > >(n_x, n_y,domain_pt);
142  }
143  else
144  {
145  Bulk_element_mesh_pt = new
146  HermiteQuadMesh<Hijacked<BiharmonicElement<2> > >(n_x, n_y,domain_pt,spacing_fn);
147  }
148 // add_sub_mesh(Bulk_element_mesh_pt);
149 // build_global_mesh();
150  }
151 
152 
153  /// \short
155  {
157  if (Face_element_mesh_pt != 0)
158  {
160  }
163  }
164 
165  /// \short Impose a load to the surface of the plate.
166  /// note : MUST be called before neumann boundary conditions are imposed,
167  /// i.e. a free edge or a simply supported edge with laplacian(u) imposed
169  {
170  //number of elements in mesh
171  unsigned n_bulk_element = Bulk_element_mesh_pt->nelement();
172 
173  //loop over bulk elements
174  for (unsigned i = 0; i < n_bulk_element; i++)
175  {
176  // upcast from generalised element to specific element
177  BiharmonicElement<2> *element_pt =
178  dynamic_cast<BiharmonicElement<2> *>
180 
181  // set the source function pointer
182  element_pt->source_fct_pt() = source_fct_pt;
183  }
184  }
185 
186  /// \short Imposes the prescribed dirichlet BCs u (u_fn) and
187  /// du/dn (dudn_fn) dirichlet BCs by 'pinning'
188  void set_dirichlet_boundary_condition(const unsigned& b,
189  DirichletBCFctPt u_fn = 0,
190  DirichletBCFctPt dudn_fn = 0);
191 
192  /// \short Imposes the prescribed Neumann BCs laplacian(u) (flux0_fct_pt) and
193  /// dlaplacian(u)/dn (flux1_fct_pt) with flux edge elements
194  void set_neumann_boundary_condition(const unsigned &b,
196  flux0_fct_pt,
198  flux1_fct_pt = 0);
199 
200  public:
201 
202  // NOTE: these two private meshes are required for the block preconditioners.
203 
204  /// \short Mesh for BiharmonicElement<DIM> only - the block preconditioner
205  /// assemble the global equation number to block number mapping from elements
206  /// in this mesh only
208 
209  /// \short mesh for face elements
211 };
212 
213 
214 
215 
216 //=============================================================================
217 /// \short Biharmonic Fluid Problem Class - describes stokes flow in 2D.
218 /// Developed for the topologically rectangular Hermite Element Mesh. Contains
219 /// functions allowing the following boundary conditions to be applied (on a
220 /// given edge):
221 /// + wall : v_n = 0 and v_t = 0 (psi must also be prescribed)
222 /// + traction free : v_t = 0
223 /// + flow : v_n and v_t are prescribed
224 /// NOTE 1 : psi is the stream function
225 /// + fluid velocity normal to boundary v_n = +/- dpsi/dt
226 /// + fluid velocity tangential to boundary v_t = -/+ dpsi/dn
227 /// NOTE 2 : when a solid wall boundary condition is applied to ensure that
228 /// v_n = 0 the the streamfunction psi must also be prescribed (and
229 /// constant)
230 //=============================================================================
231 template<unsigned DIM> class BiharmonicFluidProblem : public Problem
232 {
233 
234 public:
235 
236 
237  /// \short Definition of a dirichlet boundary condition function pointer.
238  /// Takes the position along a boundary (s) in the macro element coordinate
239  /// scheme and returns the fluid velocity normal (dpsi/dt) to the boundary
240  /// (u[0]) and the fluid velocity tangential (dpsidn) to the boundary (u[1]).
241  typedef void (*FluidBCFctPt)(const double& s, Vector<double>& u);
242 
243 
244  /// constructor
246  {
247  // initialise the number of non bulk elements
248  Npoint_element = 0;
249  }
250 
251 
252  /// actions before solve, performs self test
254  {
255 #ifdef PARANOID
256  if (0 == self_test())
257  {
258  oomph_info << "self test passed" << std::endl;
259  }
260  else
261  {
262  oomph_info << "self test failed" << std::endl;
263  }
264 #endif
265  }
266 
267 
268  /// action after solve
270  {}
271 
272 
273  /// \short documents the solution, and if an exact solution is provided, then
274  /// the error between the numerical and exact solution is presented
275  void doc_solution(DocInfo& doc_info,
276  FiniteElement::SteadyExactSolutionFctPt exact_soln_pt = 0);
277 
278 
279 protected:
280 
281 
282  /// \short Imposes a solid boundary on boundary b - no flow into boundary
283  /// or along boundary v_n = 0 and v_t = 0. User must presribe the
284  /// streamfunction psi to ensure dpsi/dt = 0 is imposed at all points on the
285  /// boundary and not just at the nodes
286  void impose_solid_boundary_on_edge(const unsigned& b,const double& psi = 0);
287 
288  /// \short Impose a traction free edge - i.e. v_t = 0 or dpsi/dn = 0. In
289  /// general dpsi/dn = 0 can only be imposed using equation elements to set
290  /// the DOFs dpsi/ds_n, however in the special case of dt/ds_n = 0, then
291  /// dpsi/ds_n = 0 and can be imposed using pinning - this is handled
292  /// automatically in this function. For a more detailed description of the
293  /// equations see the description of the class BiharmonicFluidBoundaryElement
294  void impose_traction_free_edge(const unsigned& b);
295 
296 
297  /// \short Impose a prescribed fluid flow comprising the velocity normal to
298  /// the boundary (u_imposed_fn[0]) and the velocity tangential to the
299  /// boundary (u_imposed_fn[1])
300  void impose_fluid_flow_on_edge(const unsigned& b, FluidBCFctPt u_imposed_fn);
301 
302 
303  private:
304 
305  // number of non-bulk elements - i.e. biharmonic fluid boundary elements
306  unsigned Npoint_element;
307 
308 };
309 
310 
311 
312 
313 
314 //=============================================================================
315 /// \short Point equation element used to impose the traction free edge (i.e.
316 /// du/dn = 0) on the boundary when dt/ds_n != 0. The following equation is
317 /// implemented : du/ds_n = dt/ds_n * ds_t/dt * du/dt.
318 /// The bulk biharmonic elements on the boundary must be hijackable and the
319 /// du/ds_n and d2u/ds_nds_t boundary DOFs hijacked when these elements are
320 /// applied. At any node where dt/ds_n = 0 we can impose du/ds_n = 0 and
321 /// d2u/ds_nds_t = 0 using pinning - see
322 /// BiharmonicFluidProblem::impose_traction_free_edge()
323 //=============================================================================
325 {
326 
327 
328 public :
329 
330  // constructor
331  BiharmonicFluidBoundaryElement(Node* node_pt, const unsigned s_fixed_index)
332  {
333 
334  // set the node pt
335  this->node_pt(0) = node_pt;
336 
337  // store fixed index on the boundary
338  S_fixed_index = s_fixed_index;
339 
340  }
341 
342  /// Output function -- does nothing
343  void output(std::ostream &outfile) {}
344 
345 
346  /// \short Output function -- does nothing
347  void output(std::ostream &outfile, const unsigned &n_plot) {}
348 
349 
350  /// \short Output function -- does nothing
351  void output_fluid_velocity(std::ostream &outfile, const unsigned &n_plot) {}
352 
353 
354  /// C-style output function -- does nothing
355  void output(FILE* file_pt) {}
356 
357 
358  /// \short C-style output function -- does nothing
359  void output(FILE* file_pt, const unsigned &n_plot) {}
360 
361 
362  /// compute_error -- does nothing
363  void compute_error(std::ostream &outfile,
365  exact_soln_pt, double& error, double& norm) {}
366 
367 
368  /// \short Compute the elemental residual vector - wrapper function called by
369  /// get_residuals in GeneralisedElement
371  {
372 
373  // create a dummy matrix
374  DenseDoubleMatrix dummy(1);
375 
376  // call the generic residuals functions with flag set to zero
378  dummy, 0);
379  }
380 
381 
382  /// \short Compute the elemental residual vector and jacobian matrix - wrapper
383  /// function called by get_jacobian in GeneralisedElement
385  DenseMatrix<double> &jacobian)
386  {
387 
388  // call generic routine with flag set to 1
390  jacobian,1);
391  }
392 
393 
394  /// \short Computes the elemental residual vector and the elemental jacobian
395  /// matrix if JFLAG = 0
396  /// Imposes the equations : du/ds_n = dt/ds_n * ds_t/dt * du/dt
398  Vector<double> &residuals,
399  DenseMatrix<double>& jacobian,
400  unsigned JFLAG);
401 
402 private :
403 
404  // fixed local coordinate index on boundary
405  unsigned S_fixed_index;
406 
407 };
408 
409 
410 
411 }
412 #endif
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – does nothing.
GeneralisedAxisymAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
void output_fluid_velocity(std::ostream &outfile, const unsigned &n_plot)
Output function – does nothing.
Information for documentation of results: Directory and file number to enable output in the form RESL...
void(* BiharmonicSourceFctPt)(const Vector< double > &x, double &f)
Definition of the Source Function.
cstr elem_len * i
Definition: cfortran.h:607
A two dimensional Hermite bicubic element quadrilateral mesh for a topologically rectangular domain...
The Problem class.
Definition: problem.h:152
SourceFctPt & source_fct_pt()
Access functions for the source function pointer.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the elemental residual vector and jacobian matrix - wrapper function called by get_jacobian i...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
OomphInfo oomph_info
void set_source_function(const BiharmonicSourceFctPt source_fct_pt)
Impose a load to the surface of the plate. note : MUST be called before neumann boundary conditions a...
void doc_solution(DocInfo &doc_info, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt=0)
documents the solution, and if an exact solution is provided, then the error between the numerical an...
void doc_solution(DocInfo &doc_info, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt=0)
documents the solution, and if an exact solution is provided, then the error between the numerical an...
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
Mesh * bulk_element_mesh_pt()
Access function to the bulk element mesh pt.
void(* FluidBCFctPt)(const double &s, Vector< double > &u)
Definition of a dirichlet boundary condition function pointer. Takes the position along a boundary (s...
void set_neumann_boundary_condition(const unsigned &b, BiharmonicFluxElement< 2 >::FluxFctPt flux0_fct_pt, BiharmonicFluxElement< 2 >::FluxFctPt flux1_fct_pt=0)
Imposes the prescribed Neumann BCs laplacian(u) (flux0_fct_pt) and dlaplacian(u)/dn (flux1_fct_pt) wi...
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
Mesh * Face_element_mesh_pt
mesh for face elements
void actions_after_newton_solve()
action after solve
void impose_solid_boundary_on_edge(const unsigned &b, const double &psi=0)
Imposes a solid boundary on boundary b - no flow into boundary or along boundary v_n = 0 and v_t = 0...
static char t char * s
Definition: cfortran.h:572
void actions_before_newton_solve()
actions before solve, performs self test
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
compute_error – does nothing
BiharmonicProblem()
Constructor.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void output(FILE *file_pt)
C-style output function – does nothing.
Topologically Rectangular Domain - a domain dexcribing a topologically rectangular problem - primaril...
void impose_traction_free_edge(const unsigned &b)
Impose a traction free edge - i.e. v_t = 0 or dpsi/dn = 0. In general dpsi/dn = 0 can only be imposed...
virtual void fill_in_generic_residual_contribution_biharmonic_boundary(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned JFLAG)
Computes the elemental residual vector and the elemental jacobian matrix if JFLAG = 0 Imposes the equ...
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Add a submesh to the problem and return its number, i, by which it can be accessed via mesh_pt(i)...
Definition: problem.h:1293
void actions_after_newton_solve()
action after solve
biharmonic element class
void build_bulk_mesh(const unsigned n_x, const unsigned n_y, TopologicallyRectangularDomain *domain_pt, HermiteQuadMesh< BiharmonicElement< 2 > >::MeshSpacingFnPtr spacing_fn=0)
builds the bulk mesh on a prescribed domain with a node spacing defined by spacing fn with n_x by n_y...
Point equation element used to impose the traction free edge (i.e. du/dn = 0) on the boundary when dt...
BiharmonicFluidBoundaryElement(Node *node_pt, const unsigned s_fixed_index)
void set_dirichlet_boundary_condition(const unsigned &b, DirichletBCFctPt u_fn=0, DirichletBCFctPt dudn_fn=0)
Imposes the prescribed dirichlet BCs u (u_fn) and du/dn (dudn_fn) dirichlet BCs by 'pinning'...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – does nothing.
void(* DirichletBCFctPt)(const double &s, double &u)
Definition of a dirichlet boundary condition function pointer. Takes the position along a boundary (s...
void build_global_mesh()
Build the global mesh by combining the all the submeshes. Note: The nodes boundary information refers...
Definition: problem.cc:1472
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the elemental residual vector - wrapper function called by get_residuals in GeneralisedElemen...
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1965
virtual ~BiharmonicProblem()
Destructor. Delete the meshes.
unsigned self_test()
Self-test: Check meshes and global data. Return 0 for OK.
Definition: problem.cc:13046
void actions_before_newton_solve()
actions before solve, performs self test
Biharmonic Fluid Problem Class - describes stokes flow in 2D. Developed for the topologically rectang...
Biharmonic Plate Problem Class - for problems where the load can be assumed to be acting normal to th...
void output(std::ostream &outfile)
Output function – does nothing.
A general mesh class.
Definition: mesh.h:74
Mesh * Bulk_element_mesh_pt
Mesh for BiharmonicElement<DIM> only - the block preconditioner assemble the global equation number t...
void impose_fluid_flow_on_edge(const unsigned &b, FluidBCFctPt u_imposed_fn)
Impose a prescribed fluid flow comprising the velocity normal to the boundary (u_imposed_fn[0]) and t...