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// Copyright (C) 2006-2025 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#ifndef OOMPH_BIHARMONIC_PROBLEM_HEADER
27#define OOMPH_BIHARMONIC_PROBLEM_HEADER
28
29// Config header
30#ifdef HAVE_CONFIG_H
31#include <oomph-lib-config.h>
32#endif
33
34#ifdef OOMPH_HAS_MPI
35// mpi headers
36#include "mpi.h"
37#endif
38
39// Generic C++ headers
40#include <iostream>
41#include <math.h>
42
43// oomph-lib headers
44#include "generic/problem.h"
46#include "../meshes/hermite_element_quad_mesh.h"
47#include "biharmonic_elements.h"
49
50
51namespace oomph
52{
53 //=============================================================================
54 /// Biharmonic Plate Problem Class - for problems where the load can be
55 /// assumed to be acting normal to the surface of the plate and the
56 /// deflections are small relative to the thickness of the plate. Developed
57 /// for the topologically rectangular Hermite Element Mesh. Contains functions
58 /// allowing the following boundary conditions to be applied (on a given
59 /// edge):
60 /// + clamped : u and du/dn imposed
61 /// + simply supported : u and laplacian(u) imposed
62 /// + free : laplacian(u) and dlaplacian(u)/dn imposed
63 //=============================================================================
64 template<unsigned DIM>
66 {
67 public:
68 /// Definition of a dirichlet boundary condition function pointer.
69 /// Takes the position along a boundary (s) in the macro element coordinate
70 /// scheme and returns the value of the boundary condition at that point
71 /// (u).
72 typedef void (*DirichletBCFctPt)(const double& s, double& u);
73
74
75 /// Definition of the Source Function.
76 typedef void (*BiharmonicSourceFctPt)(const Vector<double>& x, double& f);
77
78 /// Constructor
84
85 /// Destructor. Delete the meshes
87 {
90 };
91
92 /// actions before solve, performs self test
94 {
95#ifdef PARANOID
96 if (0 == self_test())
97 {
98 oomph_info << "self test passed" << std::endl;
99 }
100 else
101 {
102 oomph_info << "self test failed" << std::endl;
103 }
104#endif
105 }
106
107 /// action after solve
109
110 /// documents the solution, and if an exact solution is provided,
111 /// then the error between the numerical and exact solution is presented
112 void doc_solution(
113 DocInfo& doc_info,
115
116 /// Access function to the bulk element mesh pt
121
122 protected:
123 /// builds the bulk mesh on a prescribed domain with a node spacing
124 /// defined by spacing fn with n_x by n_y elements
126 const unsigned n_x,
127 const unsigned n_y,
129 HermiteQuadMesh<BiharmonicElement<2>>::MeshSpacingFnPtr spacing_fn = 0)
130 {
131 if (spacing_fn == 0)
132 {
135 n_x, n_y, domain_pt);
136 }
137 else
138 {
141 n_x, n_y, domain_pt, spacing_fn);
142 }
143 // add_sub_mesh(Bulk_element_mesh_pt);
144 // build_global_mesh();
145 }
146
147
148 /// Build global mesh and assign equation numbers
159
160 /// Impose a load to the surface of the plate.
161 /// note : MUST be called before neumann boundary conditions are imposed,
162 /// i.e. a free edge or a simply supported edge with laplacian(u) imposed
164 {
165 // number of elements in mesh
167
168 // loop over bulk elements
169 for (unsigned i = 0; i < n_bulk_element; i++)
170 {
171 // upcast from generalised element to specific element
172 BiharmonicElement<2>* element_pt = dynamic_cast<BiharmonicElement<2>*>(
174
175 // set the source function pointer
176 element_pt->source_fct_pt() = source_fct_pt;
177 }
178 }
179
180 /// Imposes the prescribed dirichlet BCs u (u_fn) and
181 /// du/dn (dudn_fn) dirichlet BCs by 'pinning'
182 void set_dirichlet_boundary_condition(const unsigned& b,
185
186 /// Imposes the prescribed Neumann BCs laplacian(u) (flux0_fct_pt)
187 /// and dlaplacian(u)/dn (flux1_fct_pt) with flux edge elements
189 const unsigned& b,
191 BiharmonicFluxElement<2>::FluxFctPt flux1_fct_pt = 0);
192
193 public:
194 // NOTE: these two private meshes are required for the block
195 // preconditioners.
196
197 /// Mesh for BiharmonicElement<DIM> only - the block preconditioner
198 /// assemble the global equation number to block number mapping from
199 /// elements in this mesh only
201
202 /// mesh for face elements
204 };
205
206
207 //=============================================================================
208 /// Biharmonic Fluid Problem Class - describes stokes flow in 2D.
209 /// Developed for the topologically rectangular Hermite Element Mesh. Contains
210 /// functions allowing the following boundary conditions to be applied (on a
211 /// given edge):
212 /// + wall : v_n = 0 and v_t = 0 (psi must also be prescribed)
213 /// + traction free : v_t = 0
214 /// + flow : v_n and v_t are prescribed
215 /// NOTE 1 : psi is the stream function
216 /// + fluid velocity normal to boundary v_n = +/- dpsi/dt
217 /// + fluid velocity tangential to boundary v_t = -/+ dpsi/dn
218 /// NOTE 2 : when a solid wall boundary condition is applied to ensure that
219 /// v_n = 0 the the streamfunction psi must also be prescribed (and
220 /// constant)
221 //=============================================================================
222 template<unsigned DIM>
224 {
225 public:
226 /// Definition of a dirichlet boundary condition function pointer.
227 /// Takes the position along a boundary (s) in the macro element coordinate
228 /// scheme and returns the fluid velocity normal (dpsi/dt) to the boundary
229 /// (u[0]) and the fluid velocity tangential (dpsidn) to the boundary
230 /// (u[1]).
231 typedef void (*FluidBCFctPt)(const double& s, Vector<double>& u);
232
233
234 /// constructor
236 {
237 // initialise the number of non bulk elements
238 Npoint_element = 0;
239 }
240
241
242 /// actions before solve, performs self test
244 {
245#ifdef PARANOID
246 if (0 == self_test())
247 {
248 oomph_info << "self test passed" << std::endl;
249 }
250 else
251 {
252 oomph_info << "self test failed" << std::endl;
253 }
254#endif
255 }
256
257
258 /// action after solve
260
261
262 /// documents the solution, and if an exact solution is provided,
263 /// then the error between the numerical and exact solution is presented
264 void doc_solution(
265 DocInfo& doc_info,
267
268
269 protected:
270 /// Imposes a solid boundary on boundary b - no flow into boundary
271 /// or along boundary v_n = 0 and v_t = 0. User must presribe the
272 /// streamfunction psi to ensure dpsi/dt = 0 is imposed at all points on the
273 /// boundary and not just at the nodes
274 void impose_solid_boundary_on_edge(const unsigned& b,
275 const double& psi = 0);
276
277 /// Impose a traction free edge - i.e. v_t = 0 or dpsi/dn = 0. In
278 /// general dpsi/dn = 0 can only be imposed using equation elements to set
279 /// the DOFs dpsi/ds_n, however in the special case of dt/ds_n = 0, then
280 /// dpsi/ds_n = 0 and can be imposed using pinning - this is handled
281 /// automatically in this function. For a more detailed description of the
282 /// equations see the description of the class
283 /// BiharmonicFluidBoundaryElement
284 void impose_traction_free_edge(const unsigned& b);
285
286
287 /// Impose a prescribed fluid flow comprising the velocity normal to
288 /// the boundary (u_imposed_fn[0]) and the velocity tangential to the
289 /// boundary (u_imposed_fn[1])
290 void impose_fluid_flow_on_edge(const unsigned& b,
292
293
294 private:
295 // number of non-bulk elements - i.e. biharmonic fluid boundary elements
297 };
298
299
300 //=============================================================================
301 /// Point equation element used to impose the traction free edge (i.e.
302 /// du/dn = 0) on the boundary when dt/ds_n != 0. The following equation is
303 /// implemented : du/ds_n = dt/ds_n * ds_t/dt * du/dt.
304 /// The bulk biharmonic elements on the boundary must be hijackable and the
305 /// du/ds_n and d2u/ds_nds_t boundary DOFs hijacked when these elements are
306 /// applied. At any node where dt/ds_n = 0 we can impose du/ds_n = 0 and
307 /// d2u/ds_nds_t = 0 using pinning - see
308 /// BiharmonicFluidProblem::impose_traction_free_edge()
309 //=============================================================================
311 {
312 public:
313 // constructor
315 {
316 // set the node pt
317 this->node_pt(0) = node_pt;
318
319 // store fixed index on the boundary
320 S_fixed_index = s_fixed_index;
321 }
322
323 /// Output function -- does nothing
324 void output(std::ostream& outfile) {}
325
326
327 /// Output function -- does nothing
328 void output(std::ostream& outfile, const unsigned& n_plot) {}
329
330
331 /// Output function -- does nothing
332 void output_fluid_velocity(std::ostream& outfile, const unsigned& n_plot) {}
333
334
335 /// C-style output function -- does nothing
337
338
339 /// C-style output function -- does nothing
340 void output(FILE* file_pt, const unsigned& n_plot) {}
341
342
343 /// compute_error -- does nothing
344 void compute_error(std::ostream& outfile,
346 double& error,
347 double& norm)
348 {
349 }
350
351
352 /// Compute the elemental residual vector - wrapper function called
353 /// by get_residuals in GeneralisedElement
355 {
356 // create a dummy matrix
358
359 // call the generic residuals functions with flag set to zero
361 residuals, dummy, 0);
362 }
363
364
365 /// Compute the elemental residual vector and jacobian matrix -
366 /// wrapper function called by get_jacobian in GeneralisedElement
368 DenseMatrix<double>& jacobian)
369 {
370 // call generic routine with flag set to 1
372 residuals, jacobian, 1);
373 }
374
375
376 /// Computes the elemental residual vector and the elemental jacobian
377 /// matrix if JFLAG = 0
378 /// Imposes the equations : du/ds_n = dt/ds_n * ds_t/dt * du/dt
380 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned JFLAG);
381
382 private:
383 // fixed local coordinate index on boundary
385 };
386
387
388} // namespace oomph
389#endif
static char t char * s
Definition cfortran.h:568
cstr elem_len * i
Definition cfortran.h:603
AdvectionDiffusionReactionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
Point equation element used to impose the traction free edge (i.e. du/dn = 0) on the boundary when dt...
void output(std::ostream &outfile)
Output function – does nothing.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – does nothing.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the elemental residual vector - wrapper function called by get_residuals in GeneralisedElemen...
void output(FILE *file_pt)
C-style output function – does nothing.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – does nothing.
BiharmonicFluidBoundaryElement(Node *node_pt, const unsigned s_fixed_index)
void output_fluid_velocity(std::ostream &outfile, const unsigned &n_plot)
Output function – does nothing.
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...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
compute_error – does nothing
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...
Biharmonic Fluid Problem Class - describes stokes flow in 2D. Developed for the topologically rectang...
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 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...
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....
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...
void actions_after_newton_solve()
action after solve
void(* FluidBCFctPt)(const double &s, Vector< double > &u)
Definition of a dirichlet boundary condition function pointer. Takes the position along a boundary (s...
void actions_before_newton_solve()
actions before solve, performs self test
Biharmonic Plate Problem Class - for problems where the load can be assumed to be acting normal to th...
void(* BiharmonicSourceFctPt)(const Vector< double > &x, double &f)
Definition of the Source Function.
void actions_before_newton_solve()
actions before solve, performs self test
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...
void build_global_mesh_and_assign_eqn_numbers()
Build global mesh and assign equation numbers.
Mesh * bulk_element_mesh_pt()
Access function to the bulk element mesh pt.
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...
Mesh * Bulk_element_mesh_pt
Mesh for BiharmonicElement<DIM> only - the block preconditioner assemble the global equation number t...
virtual ~BiharmonicProblem()
Destructor. Delete the meshes.
void(* DirichletBCFctPt)(const double &s, double &u)
Definition of a dirichlet boundary condition function pointer. Takes the position along a boundary (s...
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 actions_after_newton_solve()
action after solve
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 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...
Mesh * Face_element_mesh_pt
mesh for face elements
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition matrices.h:1271
Information for documentation of results: Directory and file number to enable output in the form RESL...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition elements.h:1763
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
A two dimensional Hermite bicubic element quadrilateral mesh for a topologically rectangular domain....
A general mesh class.
Definition mesh.h:67
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition mesh.h:452
unsigned long nelement() const
Return number of elements in the mesh.
Definition mesh.h:598
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
Point element has just a single node and a single shape function which is identically equal to one.
Definition elements.h:3443
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition problem.h:154
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:1350
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:2085
void build_global_mesh()
Build the global mesh by combining the all the submeshes. Note: The nodes boundary information refers...
Definition problem.cc:1589
unsigned self_test()
Self-test: Check meshes and global data. Return 0 for OK.
Definition problem.cc:13354
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Topologically Rectangular Domain - a domain dexcribing a topologically rectangular problem - primaril...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...