biharmonic_problem.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//====================================================================
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32  #include <oomph-lib-config.h>
33 #endif
34 
35 //oomph-lib includes
36 #include "biharmonic_problem.h"
37 
38 
39 namespace oomph
40 {
41 
42 //=============================================================================
43 /// \short Impose a Clamped Edge. Imposes the prescribed dirichlet BCs u and
44 /// du/dn dirichlet BCs by pinning
45 //=============================================================================
46 template<unsigned DIM>
48 set_dirichlet_boundary_condition(const unsigned& b, DirichletBCFctPt u_fn,
49  DirichletBCFctPt dudn_fn)
50 {
51  // number of nodes on boundary b
52  unsigned n_node = Bulk_element_mesh_pt->nboundary_node(b);
53 
54  // fixed faced index for boundary
55  int face_index = Bulk_element_mesh_pt->face_index_at_boundary(b,0);
56 
57  //Need to get the s_fixed_index
58  unsigned s_fixed_index = 0;
59  switch(face_index)
60  {
61  case -1:
62  case 1:
63  s_fixed_index = 0;
64  break;
65 
66  case -2:
67  case 2:
68  s_fixed_index = 1;
69  break;
70 
71  default:
72  throw OomphLibError("Face Index not +/-1 or +/-2: Need 2D QElements",
73  OOMPH_CURRENT_FUNCTION,
74  OOMPH_EXCEPTION_LOCATION);
75  }
76 
77  // node position along edge b in macro element boundary representation
78  // [-1,1]
79  Vector<double> s(1);
80 
81  //Set the edge sign
82  int edge_sign = 0;
83  switch(face_index)
84  {
85  case -1:
86  case 2:
87  edge_sign = -1;
88  break;
89 
90  case 1:
91  case -2:
92  edge_sign = 1;
93  break;
94  }
95 
96  // finite difference step
97  const double h = 10e-8;
98 
99  // node position along edge b in macro element boundary representation
100  // [-1,1]
101  Vector<double> m(2);
102 
103  // if u is prescribed then impose
104  if (u_fn != 0)
105  {
106 
107  // loop over nodes on boundary b
108  for (unsigned n = 0; n < n_node; n++)
109  {
110 
111  // find node position along edge [-1,1] in macro element representation
112  Bulk_element_mesh_pt->boundary_node_pt(b,n)->
113  get_coordinates_on_boundary(b,m);
114 
115  // get u at node
116  double u;
117  (*u_fn)(m[0],u);
118 
119  // finite difference is used to compute dudm_t
120  // u0 and u1 store values of u left/below and right/above node n
121  double u_L, u_R;
122 
123  // if left/lower corner node
124  if (n == 0)
125  {
126  (*u_fn)(m[0],u_L);
127  (*u_fn)(m[0]+h,u_R);
128  }
129  // if right/upper corner node
130  else if (n == n_node-1)
131  {
132  (*u_fn)(m[0]-h,u_L);
133  (*u_fn)(m[0],u_R);
134  }
135  // if other node
136  else
137  {
138  (*u_fn)(m[0]-0.5*h,u_L);
139  (*u_fn)(m[0]+0.5*h,u_R);
140  }
141 
142  // compute dudm_t
143  double dudm_t = (u_R-u_L)/h;
144 
145  //compute duds_t
146  double duds_t = m[1]*dudm_t;
147 
148  // pin and set u type dof
149  Bulk_element_mesh_pt->boundary_node_pt(b,n)->pin(0);
150  Bulk_element_mesh_pt->boundary_node_pt(b,n)->set_value(0,u);
151 
152  // pin and set duds_t type degree of freedom
153  Bulk_element_mesh_pt->boundary_node_pt(b,n)->pin(2-s_fixed_index);
154  Bulk_element_mesh_pt->boundary_node_pt(b,n)
155  ->set_value(2-s_fixed_index,duds_t);
156  }
157  }
158 
159  // if dudn is prescribed then impose
160  if (dudn_fn != 0)
161  {
162 
163  // loop over nodes on boundary b
164  for (unsigned n = 0; n < n_node; n++)
165  {
166  // vectors for dx_i/ds_n and dx_i/ds_t
167  Vector<double> dxds_n(2);
168  Vector<double> dxds_t(2);
169  dxds_n[0] =
170  Bulk_element_mesh_pt->boundary_node_pt(b,n)->x_gen(1+s_fixed_index,0);
171  dxds_n[1] =
172  Bulk_element_mesh_pt->boundary_node_pt(b,n)->x_gen(1+s_fixed_index,1);
173  dxds_t[0] =
174  Bulk_element_mesh_pt->boundary_node_pt(b,n)->x_gen(2-s_fixed_index,0);
175  dxds_t[1] =
176  Bulk_element_mesh_pt->boundary_node_pt(b,n)->x_gen(2-s_fixed_index,1);
177 
178  // vector for d2xi/ds_n ds_t
179  Vector<double> d2xds_nds_t(2);
180  d2xds_nds_t[0] = Bulk_element_mesh_pt->boundary_node_pt(b,n)->x_gen(3,0);
181  d2xds_nds_t[1] = Bulk_element_mesh_pt->boundary_node_pt(b,n)->x_gen(3,1);
182 
183  // compute dn/ds_n
184  double dnds_n = ((dxds_n[0]*dxds_t[1])-(dxds_n[1]*dxds_t[0]))
185  / (sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1])*edge_sign);
186 
187  // compute dt/ds_n
188  double dtds_n = ((dxds_n[0]*dxds_t[0])+(dxds_n[1]*dxds_t[1]))
189  /sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1]);
190 
191  // compute dt/ds_t
192  double dtds_t = sqrt(pow(dxds_t[0],2)+pow(dxds_t[1],2));
193 
194  // compute ds_n/dn
195  double ds_ndn = -1.0*(edge_sign*sqrt(pow(dxds_t[0],2)+pow(dxds_t[1],2)))/
196  (dxds_t[0]*dxds_n[1]-dxds_n[0]*dxds_t[1]);
197 
198  // compute ds_t/dt
199  double ds_tdt = 1/sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1]);
200 
201  // compute d2t/ds_nds_t
202  double d2tds_nds_t = (dxds_t[0]*d2xds_nds_t[0] + dxds_t[1]*d2xds_nds_t[1])
203  /sqrt(dxds_t[0]*dxds_t[0] + dxds_t[1]*dxds_t[1]);
204 
205  // compute d2s_t/ds_ndt
206  double d2s_tds_ndt = (dxds_t[0]*d2xds_nds_t[0] + dxds_t[1]*d2xds_nds_t[1])
207  /pow(dxds_t[0]*dxds_t[0] + dxds_t[1]*dxds_t[1],3.0/2.0);
208 
209  // get m_t and dm_t/ds_t for this node
210  Vector<double> m_N(2);
211  Bulk_element_mesh_pt->boundary_node_pt(b,n)
212  ->get_coordinates_on_boundary(b,m_N);
213 
214  // compute d2u/dm_t2 and d/dm_t(dudn) by finite difference
215  double d2udm_t2 = 0;
216  double ddm_tdudn = 0;
217  double u_2L,u_N,u_2R,dudn_L,dudn_R;
218  // evaluate u_fn and dudn_fn for current node
219  if (n == 0)
220  {
221  (*u_fn)(m_N[0],u_2L);
222  (*u_fn)(m_N[0]+h,u_N);
223  (*u_fn)(m_N[0]-2*h,u_2R);
224  (*dudn_fn)(m_N[0],dudn_L);
225  (*dudn_fn)(m_N[0]+h,dudn_R);
226  }
227  else if (n == n_node-1)
228  {
229  (*u_fn)(m_N[0]-2*h,u_2L);
230  (*u_fn)(m_N[0]-h,u_N);
231  (*u_fn)(m_N[0],u_2R);
232  (*dudn_fn)(m_N[0]-h,dudn_L);
233  (*dudn_fn)(m_N[0],dudn_R);
234  }
235  else
236  {
237  (*u_fn)(m_N[0]-h,u_2L);
238  u_N = Bulk_element_mesh_pt->boundary_node_pt(b,n)->value(0);
239  (*u_fn)(m_N[0]+h,u_2R);
240  (*dudn_fn)(m_N[0]-0.5*h,dudn_L);
241  (*dudn_fn)(m_N[0]+0.5*h,dudn_R);
242  }
243  // compute
244  d2udm_t2 = (u_2L+u_2R-2*u_N)/(h*h);
245  ddm_tdudn = (dudn_R-dudn_L)/h;
246 
247  // get dudn at the node
248  double dudn;
249  (*dudn_fn)(m_N[0],dudn);
250 
251  // compute d2u/ds_t2
252  double d2uds_t2 = m_N[1]*m_N[1]*d2udm_t2;
253 
254  // get duds_t
255  double duds_t = Bulk_element_mesh_pt->
256  boundary_node_pt(b,n)->value(2-s_fixed_index);
257 
258  // get du/dt
259  double dudt = ds_tdt*duds_t;
260 
261  // compute d2u/dndt
262  double d2udndt = ds_tdt= dtds_t*m_N[1]*ddm_tdudn;
263 
264  // compute d2udt2
265  double dtds_nd2udt2 = edge_sign*(dxds_t[0]*dxds_n[1]-
266  dxds_n[0]*dxds_t[1])*
267  (ds_tdt*(d2udndt - ds_ndn*(d2s_tds_ndt*dudt+
268  ds_tdt*d2uds_t2)));
269 
270  // compute dds_n(dudt)
271  double dds_ndudt = dtds_nd2udt2 + dnds_n*d2udndt;
272 
273  // compute du/ds_n
274  double duds_n = dnds_n*dudn + dtds_n*ds_tdt*duds_t;
275 
276  // compute d2u/ds_nds_t
277  double d2uds_nds_t = d2tds_nds_t*dudt+dtds_t*dds_ndudt;
278 
279  // pin du/ds_n dof and set value
280  Bulk_element_mesh_pt->boundary_node_pt(b,n)->pin(1+s_fixed_index);
281  Bulk_element_mesh_pt->boundary_node_pt(b,n)->
282  set_value(1+s_fixed_index,duds_n);
283 
284  // pin d2u/ds_nds_t dof and set value
285  Bulk_element_mesh_pt->boundary_node_pt(b,n)->pin(3);
286  Bulk_element_mesh_pt->boundary_node_pt(b,n)->set_value(3,d2uds_nds_t);
287  }
288  }
289 }
290 
291 
292 
293 //=============================================================================
294 /// \short Imposes a 'free' edge. Imposes the prescribed Neumann BCs
295 /// laplacian(u) and dlaplacian(u)/dn with flux edge elements
296 //=============================================================================
297 template<unsigned DIM>
301  flux0_fct_pt,
303  flux1_fct_pt)
304 {
305 
306  // if the face element mesh pt does not exist then build it
307  if (Face_element_mesh_pt == 0)
308  {
309  Face_element_mesh_pt = new Mesh();
310  }
311 
312  // How many bulk elements are adjacent to boundary b?
313  unsigned n_element = Bulk_element_mesh_pt->nboundary_element(b);
314 
315  // Loop over the bulk elements adjacent to boundary b?
316  for(unsigned e=0;e<n_element;e++)
317  {
318 
319  // Get pointer to the bulk element that is adjacent to boundary b
320  FiniteElement* bulk_elem_pt =
321  dynamic_cast<FiniteElement* >(Bulk_element_mesh_pt->boundary_element_pt(b,e));
322 
323  // What is the face index along the boundary
324  int face_index = Bulk_element_mesh_pt->face_index_at_boundary(b,e);
325 
326  // Build the corresponding prescribed-flux element
327  BiharmonicFluxElement<2> *flux_element_pt = new
329  (bulk_elem_pt, face_index, b);
330 
331 
332  // pass the flux BC pointers to the flux elements
333  flux_element_pt->flux0_fct_pt() = flux0_fct_pt;
334  if (flux1_fct_pt != 0)
335  {
336  flux_element_pt->flux1_fct_pt() = flux1_fct_pt;
337  }
338 
339  //Add the prescribed-flux element to the mesh
340  Face_element_mesh_pt->add_element_pt(flux_element_pt);
341  }
342 }
343 
344 
345 //=============================================================================
346 /// \short documents the solution, and if an exact solution is provided, then
347 /// the error between the numerical and exact solution is presented
348 //=============================================================================
349 template<unsigned DIM>
353 {
354  std::ofstream some_file;
355  std::ostringstream filename;
356 
357  // Number of plot points: npts x npts
358  unsigned npts=5;
359 
360  // Output solution
361  filename << doc_info.directory() << "/soln_" << doc_info.number() << ".dat";
362  some_file.open(filename.str().c_str());
363  Bulk_element_mesh_pt->output(some_file,npts);
364  some_file.close();
365 
366  // if the exact solution is not provided
367  if (exact_soln_pt != 0)
368  {
369 
370  // Output exact solution
371  filename.str("");
372  filename << doc_info.directory() << "/exact_soln_"
373  << doc_info.number() << ".dat";
374  some_file.open(filename.str().c_str());
375  Bulk_element_mesh_pt->output_fct(some_file,npts, exact_soln_pt);
376  some_file.close();
377 
378  // Doc error and return of the square of the L2 error
379  double error,norm;
380  filename.str("");
381  filename << doc_info.directory() << "/error_"
382  << doc_info.number() << ".dat";
383  some_file.open(filename.str().c_str());
384  Bulk_element_mesh_pt->compute_error(some_file, exact_soln_pt, error, norm);
385  some_file.close();
386 
387  // Doc L2 error and norm of solution
388  oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;
389  oomph_info << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
390  }
391 }
392 
393 
394 
395 
396 //=============================================================================
397 /// \short Computes the elemental residual vector and the elemental jacobian
398 /// matrix if JFLAG = 0
399 /// Imposes the equations : du/ds_n = dt/ds_n * ds_t/dt * du/dt
400 //=============================================================================
403  Vector<double> &residual,
404  DenseMatrix<double>& jacobian,
405  unsigned JFLAG)
406 {
407 
408  // dof # corresponding to d/ds_n
409  unsigned k_normal = 1+S_fixed_index;
410 
411  // dof # corresponding to d/ds_t
412  unsigned k_tangential = 2-S_fixed_index;
413 
414  // vectors for dx_i/ds_n and dx_i/ds_t
415  Vector<double> dxds_n(2);
416  Vector<double> dxds_t(2);
417  dxds_n[0] = this->node_pt(0)->x_gen(1+S_fixed_index,0);
418  dxds_n[1] = this->node_pt(0)->x_gen(1+S_fixed_index,1);
419  dxds_t[0] = this->node_pt(0)->x_gen(2-S_fixed_index,0);
420  dxds_t[1] = this->node_pt(0)->x_gen(2-S_fixed_index,1);
421 
422  // vector for d2xi/ds_n ds_t
423  Vector<double> d2xds_nds_t(2);
424  d2xds_nds_t[0] = this->node_pt(0)->x_gen(3,0);
425  d2xds_nds_t[1] = this->node_pt(0)->x_gen(3,1);
426 
427  // compute dt/ds_n
428  double dtds_n = ((dxds_n[0]*dxds_t[0])+(dxds_n[1]*dxds_t[1]))
429  /sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1]);
430 
431  // compute ds_t/dt
432  double ds_tdt = 1/sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1]);
433 
434  // determine the local equation number
435  int local_eqn_number = this->nodal_local_eqn(0,k_normal);
436 
437  // if local equation number equal to -1 then its a boundary node(pinned)
438  if (local_eqn_number >= 0)
439  {
440 
441  // additions to residual for duds_n
442  residual[local_eqn_number] += (this->node_pt(0)->value(k_normal) -
443  dtds_n*ds_tdt *
444  this->node_pt(0)->value(k_tangential));
445 
446  // if contributions to jacobian are required
447  if (JFLAG == 1)
448  {
449 
450  // additions to jacobian for du/ds_n
451  int local_dof_number = this->nodal_local_eqn(0,k_normal);
452  if (local_dof_number >= 0)
453  {
454  jacobian(local_eqn_number,local_dof_number) += 1.0;
455  }
456  local_dof_number = this->nodal_local_eqn(0,k_tangential);
457  if (local_dof_number >= 0)
458  {
459  jacobian(local_eqn_number,local_dof_number) -= dtds_n*ds_tdt;
460  }
461  }
462  }
463 }
464 
465 
466 
467 
468 
469 
470 //=============================================================================
471 /// \short Imposes a solid boundary - no flow into boundary or along boundary
472 /// v_n = 0 and v_t = 0. User must presribe the streamfunction psi to ensure
473 /// dpsi/dt = 0 is imposed at all points on the boundary and not just at the
474 /// nodes
475 //=============================================================================
476 template<unsigned DIM>
478 impose_solid_boundary_on_edge(const unsigned& b, const double& psi)
479 {
480  // number of nodes on boundary b
481  unsigned n_node = mesh_pt()->nboundary_node(b);
482 
483  // loop over nodes on boundary b
484  for (unsigned n = 0; n < n_node; n++)
485  {
486  // loop over DOFs to be pinned du/ds_n, du/ds_t and d2u/ds_nds_t
487  for ( unsigned k = 1; k < 4; k++)
488  {
489  // pin and zero DOF k
490  mesh_pt()->boundary_node_pt(b,n)->pin(k);
491  mesh_pt()->boundary_node_pt(b,n)->set_value(k,0.0);
492  }
493  mesh_pt()->boundary_node_pt(b,n)->pin(0);
494  mesh_pt()->boundary_node_pt(b,n)->set_value(0,psi);
495  }
496 }
497 
498 
499 
500 //=============================================================================
501 /// \short Impose a traction free edge - i.e. v_t = 0 or dpsi/dn = 0. In
502 /// general dpsi/dn = 0 can only be imposed using equation elements to set
503 /// the DOFs dpsi/ds_n, however in the special case of dt/ds_n = 0, then
504 /// dpsi/ds_n = 0 and can be imposed using pinning - this is handled
505 /// automatically in this function. For a more detailed description of the
506 /// equations see the description of the class : BiharmonicFluidBoundaryElement
507 //=============================================================================
508 template<unsigned DIM>
510 impose_traction_free_edge(const unsigned& b)
511 {
512  // fixed faced index for boundary
513  int face_index = mesh_pt()->face_index_at_boundary(b,0);
514 
515  //Need to get the s_fixed_index
516  unsigned s_fixed_index = 0;
517  switch(face_index)
518  {
519  case -1:
520  case 1:
521  s_fixed_index = 0;
522  break;
523 
524  case -2:
525  case 2:
526  s_fixed_index = 1;
527  break;
528 
529  default:
530  throw OomphLibError("Face Index not +/-1 or +/-2: Need 2D QElements",
531  OOMPH_CURRENT_FUNCTION,
532  OOMPH_EXCEPTION_LOCATION);
533  }
534 
535  // create a point to a hijacked biharmonic element
536  Hijacked<BiharmonicElement<2> >* hijacked_element_pt;
537 
538  // vectors for dx_i/ds_n and dx_i/ds_t
539  Vector<double> dxds_n(2);
540  Vector<double> dxds_t(2);
541 
542  // number of nodes along edge
543  unsigned n_node = mesh_pt()->nboundary_node(b);
544 
545  // loop over boudnary nodes
546  for (unsigned n = 0; n < n_node; n++)
547  {
548 
549  // get dx_i/ds_t and dx_i/ds_n at node n
550  dxds_n[0] = mesh_pt()->boundary_node_pt(b,n)->x_gen(1+s_fixed_index,0);
551  dxds_n[1] = mesh_pt()->boundary_node_pt(b,n)->x_gen(1+s_fixed_index,1);
552  dxds_t[0] = mesh_pt()->boundary_node_pt(b,n)->x_gen(2-s_fixed_index,0);
553  dxds_t[1] = mesh_pt()->boundary_node_pt(b,n)->x_gen(2-s_fixed_index,1);
554 
555  // compute dt/ds_n
556  double dtds_n = ((dxds_n[0]*dxds_t[0])+(dxds_n[1]*dxds_t[1]))
557  /sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1]);
558 
559  // if dt/ds_n = 0 we can impose the traction free edge at this node by
560  // pinning dpsi/ds_n = 0, otherwise the equation elements
561  // BiharmonicFluidBoundaryElement are used
562  if (dtds_n == 0.0)
563  {
564  // pin dpsi/dn
565  mesh_pt()->boundary_node_pt(b,n)->pin(1+s_fixed_index);
566  mesh_pt()->boundary_node_pt(b,n)->set_value(1+s_fixed_index,0.0);
567  }
568 
569  // otherwise impose equation elements
570  else
571  {
572  // hijack DOFs in element either side of node
573  // boundary 0
574  if (b == 0)
575  {
576  // hijack DOFs in left element
577  if (n > 0)
578  {
579  hijacked_element_pt =
580  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
581  (mesh_pt()->boundary_element_pt(b,n-1));
582  delete hijacked_element_pt->hijack_nodal_value(1,1+s_fixed_index);
583  }
584  // hijack DOFs in right element
585  if (n < (n_node-1))
586  {
587  hijacked_element_pt =
588  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
589  (mesh_pt()->boundary_element_pt(b,n));
590  delete hijacked_element_pt->hijack_nodal_value(0,1+s_fixed_index);
591  }
592  }
593  // boundary 1
594  else if (b == 1)
595  {
596  // hijack DOFs in left element
597  if (n > 0)
598  {
599  hijacked_element_pt =
600  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
601  (mesh_pt()->boundary_element_pt(b,n-1));
602  delete hijacked_element_pt->hijack_nodal_value(3,1+s_fixed_index);
603  }
604  // hijack DOFs in right element
605  if (n < n_node-1)
606  {
607  hijacked_element_pt =
608  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
609  (mesh_pt()->boundary_element_pt(b,n));
610  delete hijacked_element_pt->hijack_nodal_value(1,1+s_fixed_index);
611  }
612  }
613 
614  // boundary 2
615  else if (b == 2)
616  {
617  // hijack DOFs in left element
618  if (n > 0)
619  {
620  hijacked_element_pt =
621  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
622  (mesh_pt()->boundary_element_pt(b,n-1));
623  delete hijacked_element_pt->hijack_nodal_value(3,1+s_fixed_index);
624  }
625  if (n < n_node-1)
626  {
627  // hijack DOFs in right element
628  hijacked_element_pt =
629  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
630  (mesh_pt()->boundary_element_pt(b,n));
631  delete hijacked_element_pt->hijack_nodal_value(2,1+s_fixed_index);
632  }
633  }
634  // boundary 3
635  else if (b == 3)
636  {
637  // hijack DOFs in left element
638  if (n > 0)
639  {
640  hijacked_element_pt =
641  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
642  (mesh_pt()->boundary_element_pt(b,n-1));
643  delete hijacked_element_pt->hijack_nodal_value(2,1+s_fixed_index);
644  }
645  if (n < n_node-1)
646  {
647  // hijack DOFs in right element
648  hijacked_element_pt =
649  dynamic_cast<Hijacked<BiharmonicElement<2> > *>
650  (mesh_pt()->boundary_element_pt(b,n));
651  delete hijacked_element_pt->hijack_nodal_value(0,1+s_fixed_index);
652  }
653  }
654 
655  // create the boundary point element
656  BiharmonicFluidBoundaryElement* boundary_point_element_pt =
657  new BiharmonicFluidBoundaryElement(mesh_pt()->boundary_node_pt(b,n),
658  s_fixed_index);
659 
660  // add element to mesh
661  mesh_pt()->add_element_pt(boundary_point_element_pt);
662 
663  // increment number of non bulk elements
664  Npoint_element++;
665  }
666  }
667 }
668 
669 
670 //=============================================================================
671 /// \short Impose a prescribed fluid flow comprising the velocity normal to
672 /// the boundary (u_imposed_fn[0]) and the velocity tangential to the
673 /// boundary (u_imposed_fn[1])
674 //=============================================================================
675 template<unsigned DIM>
677 impose_fluid_flow_on_edge(const unsigned& b,FluidBCFctPt u_imposed_fn)
678 {
679 
680  // number of nodes on boundary b
681  unsigned n_node = mesh_pt()->nboundary_node(b);
682 
683  // fixed faced index for boundary
684  int face_index = mesh_pt()->face_index_at_boundary(b,0);
685 
686  //Need to get the s_fixed_index
687  unsigned s_fixed_index = 0;
688  //Also set the edge sign
689  int edge_sign= 0;
690 
691  switch(face_index)
692  {
693  case -1:
694  s_fixed_index = 0;
695  edge_sign = -1;
696  break;
697 
698  case 1:
699  s_fixed_index = 0;
700  edge_sign = 1;
701  break;
702 
703  case -2:
704  s_fixed_index = 1;
705  edge_sign = 1;
706  break;
707 
708  case 2:
709  s_fixed_index = 1;
710  edge_sign = -1;
711  break;
712 
713  default:
714  throw OomphLibError("Face Index not +/-1 or +/-2: Need 2D QElements",
715  OOMPH_CURRENT_FUNCTION,
716  OOMPH_EXCEPTION_LOCATION);
717  }
718 
719 
720  // node position along edge b in macro element boundary representation
721  // [-1,1]
722  Vector<double> s(1);
723 
724  // finite difference step
725  const double h = 10e-8;
726 
727  // loop over nodes on boundary b
728  for (unsigned n = 0; n < n_node; n++)
729  {
730 
731  // get m_t and dm_t/ds_t for this node
732  Vector<double> m_N(2);
733  mesh_pt()->boundary_node_pt(b,n)->get_coordinates_on_boundary(b,m_N);
734 
735  // vectors for dx_i/ds_n and dx_i/ds_t
736  Vector<double> dxds_n(2);
737  Vector<double> dxds_t(2);
738  dxds_n[0] = mesh_pt()->boundary_node_pt(b,n)->x_gen(1+s_fixed_index,0);
739  dxds_n[1] = mesh_pt()->boundary_node_pt(b,n)->x_gen(1+s_fixed_index,1);
740  dxds_t[0] = mesh_pt()->boundary_node_pt(b,n)->x_gen(2-s_fixed_index,0);
741  dxds_t[1] = mesh_pt()->boundary_node_pt(b,n)->x_gen(2-s_fixed_index,1);
742 
743  // vector for d2xi/ds_n ds_t
744  Vector<double> d2xds_nds_t(2);
745  d2xds_nds_t[0] = mesh_pt()->boundary_node_pt(b,n)->x_gen(3,0);
746  d2xds_nds_t[1] = mesh_pt()->boundary_node_pt(b,n)->x_gen(3,1);
747 
748  // compute dn/ds_n
749  double dnds_n = ((dxds_n[0]*dxds_t[1])-(dxds_n[1]*dxds_t[0]))
750  / (sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1])*edge_sign);
751 
752  // compute dt/ds_n
753  double dtds_n = ((dxds_n[0]*dxds_t[0])+(dxds_n[1]*dxds_t[1]))
754  /sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1]);
755 
756  // compute dt/ds_t
757  double dtds_t = sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1]);
758 
759  // compute d2t/ds_nds_t
760  double d2tds_nds_t = (dxds_t[0]*d2xds_nds_t[0] + dxds_t[1]*d2xds_nds_t[1])
761  /sqrt(dxds_t[0]*dxds_t[0] + dxds_t[1]*dxds_t[1]);
762 
763  // get imposed velocities
764  Vector<double> u(2);
765  (*u_imposed_fn)(m_N[0],u);
766  u[0]*=edge_sign;
767  u[1]*=-edge_sign;
768 
769  // compute d/dm_t(dpsidn) by finite difference
770  double ddm_tdudn = 0;
771  double ddm_tdudt = 0;
772  Vector<double> u_L(2);
773  Vector<double> u_R(2);
774  // evaluate dudn_fn about current node
775  if (n == 0)
776  {
777  (*u_imposed_fn)(m_N[0],u_L);
778  (*u_imposed_fn)(m_N[0]+h,u_R);
779  }
780  else if (n == n_node-1)
781  {
782  (*u_imposed_fn)(m_N[0]-h,u_L);
783  (*u_imposed_fn)(m_N[0],u_R);
784  }
785  else
786  {
787  (*u_imposed_fn)(m_N[0]-0.5*h,u_L);
788  (*u_imposed_fn)(m_N[0]+0.5*h,u_R);
789  }
790  // compute
791  ddm_tdudn = (u_R[1]-u_L[1])/h;
792  ddm_tdudt = (u_R[0]-u_L[0])/h;
793 
794  // compute du/ds_t
795  double duds_t = dtds_t*u[0];
796 
797  // compute du/ds_n
798  double duds_n = dnds_n*u[1] + dtds_n*u[0];
799 
800  // compute d2u/ds_n ds_t
801  double d2uds_nds_t = dnds_n*m_N[1]*ddm_tdudn + d2tds_nds_t*u[0]
802  + dtds_n*m_N[1]*ddm_tdudt;
803 
804  // pin du/ds_n dof and set value
805  mesh_pt()->boundary_node_pt(b,n)->pin(1+s_fixed_index);
806  mesh_pt()->boundary_node_pt(b,n)->set_value(1+s_fixed_index,duds_n);
807 
808  // pin du/ds_t dof and set value
809  mesh_pt()->boundary_node_pt(b,n)->pin(2-s_fixed_index);
810  mesh_pt()->boundary_node_pt(b,n)->set_value(2-s_fixed_index,duds_t);
811 
812  // pin du/ds_t dof and set value
813  mesh_pt()->boundary_node_pt(b,n)->pin(3);
814  mesh_pt()->boundary_node_pt(b,n)->set_value(3,d2uds_nds_t);
815  }
816 }
817 
818 
819 //=============================================================================
820 /// \short documents the solution, and if an exact solution is provided, then
821 /// the error between the numerical and exact solution is presented
822 //=============================================================================
823 template<unsigned DIM>
827 {
828  // create an output stream
829  std::ofstream some_file;
830  std::ostringstream filename;
831 
832  // Number of plot points: npts x npts
833  unsigned npts=5;
834 
835  // Output solution
836  filename << doc_info.directory() << "/soln_"
837  << doc_info.label() << ".dat";
838  some_file.open(filename.str().c_str());
839  mesh_pt()->output(some_file,npts);
840  some_file.close();
841 
842  // Output fluid velocity solution
843  filename.str("");
844  filename << doc_info.directory() << "/soln_velocity_"
845  << doc_info.label() << ".dat";
846  some_file.open(filename.str().c_str());
847  unsigned n_element = mesh_pt()->nelement();
848  for (unsigned i = 0; i < n_element-Npoint_element; i++)
849  {
850  BiharmonicElement<2>* biharmonic_element_pt =
851  dynamic_cast<BiharmonicElement<2>* >(mesh_pt()->element_pt(i));
852  biharmonic_element_pt->output_fluid_velocity(some_file,npts);
853  }
854  some_file.close();
855 
856  // if the exact solution is not provided
857  if (exact_soln_pt != 0)
858  {
859 
860  // Output exact solution
861  filename.str("");
862  filename << doc_info.directory() << "/exact_soln_"
863  << doc_info.label() << ".dat";
864  some_file.open(filename.str().c_str());
865  mesh_pt()->output_fct(some_file,npts, exact_soln_pt);
866  some_file.close();
867 
868  // Doc error and return of the square of the L2 error
869  double error,norm;
870  filename.str("");
871  filename << doc_info.directory() << "/error_"
872  << doc_info.label() << ".dat";
873  some_file.open(filename.str().c_str());
874  mesh_pt()->compute_error(some_file, exact_soln_pt, error, norm);
875  some_file.close();
876 
877  // Doc L2 error and norm of solution
878  oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;
879  oomph_info << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
880  }
881 }
882 
883 // ensure build
884 template class BiharmonicFluidProblem<2>;
885 template class BiharmonicProblem<2>;
886 
887 }
888 
889 
890 
891 
892 
893 
894 
895 
Data * hijack_nodal_value(const unsigned &n, const unsigned &i)
Hijack the i-th value stored at node n. Return a custom-made (copied) data object that contains only ...
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
Information for documentation of results: Directory and file number to enable output in the form RESL...
cstr elem_len * i
Definition: cfortran.h:607
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
std::string & label()
String used (e.g.) for labeling output files.
A general Finite Element class.
Definition: elements.h:1271
FluxFctPt & flux0_fct_pt()
Access function for the flux0 function pointer.
void output_fluid_velocity(std::ostream &outfile, const unsigned &nplot)
output fluid velocity field
OomphInfo oomph_info
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...
e
Definition: cfortran.h:575
unsigned & number()
Number used (e.g.) for labeling output files.
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 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 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(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
biharmonic element class
std::string directory() const
Output directory.
Point equation element used to impose the traction free edge (i.e. du/dn = 0) on the boundary when dt...
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'...
Hijacked elements are elements in which one or more Data values that affect the element's residuals...
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...
A general mesh class.
Definition: mesh.h:74
FluxFctPt & flux1_fct_pt()
Access function for the flux1 function pointer.
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...