two_d_poisson_flux_bc_adapt.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 //Driver for a simple 2D Poisson problem with flux boundary conditions
31 //uses two separate meshes for bulk and surface mesh
32 
33 #ifdef OOMPH_HAS_MPI
34 #include "mpi.h"
35 #endif
36 
37 //Generic routines
38 #include "generic.h"
39 
40 // The Poisson equations
41 #include "poisson.h"
42 
43 // The mesh
44 #include "meshes/simple_rectangular_quadmesh.h"
45 
46 // We need to include the refineable_quad_mesh's
47 // templated source here to allow the build of
48 // all templates.
49 //#include "generic/refineable_quad_mesh.template.cc"
50 
51 using namespace std;
52 
53 using namespace oomph;
54 
55 
56 //==============================start_of_mesh======================
57 /// Refineable equivalent of the SimpleRectangularQuadMesh.
58 /// Refinement is performed by the QuadTree-based procedures
59 /// implemented in the RefineableQuadMesh base class.
60 //=================================================================
61 template<class ELEMENT>
63  public virtual SimpleRectangularQuadMesh<ELEMENT>,
64  public RefineableQuadMesh<ELEMENT>
65 {
66 
67 public:
68 
69  /// \short Pass number of elements in the horizontal
70  /// and vertical directions, and the corresponding dimensions.
71  /// Timestepper defaults to Static.
73  const unsigned &Ny,
74  const double &Lx, const double &Ly,
75  TimeStepper* time_stepper_pt=
76  &Mesh::Default_TimeStepper) :
77  SimpleRectangularQuadMesh<ELEMENT>(Nx,Ny,Lx,Ly,time_stepper_pt)
78  {
79  // Nodal positions etc. were created in constructor for
80  // SimpleRectangularQuadMesh<...> --> We only need to set up
81  // adaptivity information: Associate finite elements with their
82  // QuadTrees and plant them in a QuadTreeForest:
83  this->setup_quadtree_forest();
84 
85  } // end of constructor
86 
87 
88  /// Destructor: Empty
90 
91 }; // end of mesh
92 
93 
94 
95 //===== start_of_namespace=============================================
96 /// Namespace for exact solution for Poisson equation with "sharp step"
97 //=====================================================================
98 namespace TanhSolnForPoisson
99 {
100 
101  /// Parameter for steepness of "step"
102  double Alpha=1.0;
103 
104  /// Parameter for angle Phi of "step"
105  double TanPhi=0.0;
106 
107  /// Exact solution as a Vector
108  void get_exact_u(const Vector<double>& x, Vector<double>& u)
109  {
110  u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
111  }
112 
113  /// Source function required to make the solution above an exact solution
114  void source_function(const Vector<double>& x, double& source)
115  {
116  source = 2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
117  (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*
118  Alpha*Alpha*TanPhi*TanPhi+2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
119  (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*Alpha;
120  }
121 
122  /// Flux required by the exact solution on a boundary on which x is fixed
123  void prescribed_flux_on_fixed_x_boundary(const Vector<double>& x,
124  double& flux)
125  {
126  //The outer unit normal to the boundary is (1,0)
127  double N[2] = {1.0, 0.0};
128  //The flux in terms of the normal is
129  flux =
130  -(1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*TanPhi*N[0]+(
131  1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*N[1];
132  }
133 
134 } // end of namespace
135 
136 
137 
138 //========= start_of_problem_class=====================================
139 /// 2D Poisson problem on rectangular domain, discretised with
140 /// 2D QPoisson elements. Flux boundary conditions are applied
141 /// along boundary 1 (the boundary where x=L). The specific type of
142 /// element is specified via the template parameter.
143 //====================================================================
144 template<class ELEMENT>
146 {
147 
148 public:
149 
150  /// Constructor: Pass pointer to source function
151  RefineableTwoMeshFluxPoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt);
152 
153  /// Destructor (empty)
155 
156  /// Doc the solution. DocInfo object stores flags/labels for where the
157  /// output gets written to
158  void doc_solution(DocInfo& doc_info);
159 
160 
161 private:
162 
163  /// \short Update the problem specs before solve: Reset boundary conditions
164  /// to the values from the exact solution.
165  void actions_before_newton_solve();
166 
167  /// Update the problem specs after solve (empty)
169 
170  /// Actions before adapt: Wipe the mesh of prescribed flux elements
171  void actions_before_adapt();
172 
173  /// Actions after adapt: Rebuild the mesh of prescribed flux elements
174  void actions_after_adapt();
175 
176  /// \short Create Poisson flux elements on boundary b of the Mesh pointed
177  /// to by bulk_mesh_pt and add them to the Mesh object pointed to by
178  /// surface_mesh_pt
179  void create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
180  Mesh* const &surface_mesh_pt);
181 
182  /// \short Delete Poisson flux elements and wipe the surface mesh
183  void delete_flux_elements(Mesh* const &surface_mesh_pt);
184 
185  /// \short Set pointer to prescribed-flux function for all
186  /// elements in the surface mesh
187  void set_prescribed_flux_pt();
188 
189  /// Pointer to the "bulk" mesh
191 
192  /// Pointer to the "surface" mesh
194 
195  /// Pointer to source function
196  PoissonEquations<2>::PoissonSourceFctPt Source_fct_pt;
197 
198 }; // end of problem class
199 
200 
201 
202 
203 //=======start_of_constructor=============================================
204 /// Constructor for Poisson problem: Pass pointer to source function.
205 //========================================================================
206 template<class ELEMENT>
208 RefineableTwoMeshFluxPoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt)
209  : Source_fct_pt(source_fct_pt)
210 {
211 
212  // Setup "bulk" mesh
213 
214  // # of elements in x-direction
215  unsigned n_x=4;
216 
217  // # of elements in y-direction
218  unsigned n_y=4;
219 
220  // Domain length in x-direction
221  double l_x=1.0;
222 
223  // Domain length in y-direction
224  double l_y=2.0;
225 
226  // Build "bulk" mesh
227  Bulk_mesh_pt=new
229 
230  // Create/set error estimator
231  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
232 
233  // Create "surface mesh" that will contain only the prescribed-flux
234  // elements. The constructor just creates the mesh without
235  // giving it any elements, nodes, etc.
236  Surface_mesh_pt = new Mesh;
237 
238  // Create prescribed-flux elements from all elements that are
239  // adjacent to boundary 1, but add them to a separate mesh.
240  // Note that this is exactly the same function as used in the
241  // single mesh version of the problem, we merely pass different Mesh pointers.
243 
244  // Add the two sub meshes to the problem
245  add_sub_mesh(Bulk_mesh_pt);
246  add_sub_mesh(Surface_mesh_pt);
247 
248  // Rebuild the Problem's global mesh from its various sub-meshes
249  build_global_mesh();
250 
251 
252  // Set the boundary conditions for this problem: All nodes are
253  // free by default -- just pin the ones that have Dirichlet conditions
254  // here.
255  unsigned n_bound = Bulk_mesh_pt->nboundary();
256  for(unsigned b=0;b<n_bound;b++)
257  {
258  //Leave nodes on boundary 1 free
259  if (b!=1)
260  {
261  unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
262  for (unsigned n=0;n<n_node;n++)
263  {
264  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
265  }
266  }
267  }
268 
269  // Complete the build of all elements so they are fully functional
270 
271  // Loop over the Poisson bulk elements to set up element-specific
272  // things that cannot be handled by constructor: Pass pointer to
273  // source function
274  unsigned n_element = Bulk_mesh_pt->nelement();
275  for(unsigned e=0;e<n_element;e++)
276  {
277  // Upcast from GeneralisedElement to Poisson bulk element
278  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
279 
280  //Set the source function pointer
281  el_pt->source_fct_pt() = Source_fct_pt;
282  }
283 
284  // Set pointer to prescribed flux function for flux elements
286 
287  // Setup equation numbering scheme
288  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
289 
290 } // end of constructor
291 
292 
293 
294 
295 //====================start_of_actions_before_newton_solve=======================
296 /// Update the problem specs before solve: Reset boundary conditions
297 /// to the values from the exact solution.
298 //========================================================================
299 template<class ELEMENT>
301 {
302  // How many boundaries are in the bulk mesh?
303  unsigned n_bound = Bulk_mesh_pt->nboundary();
304 
305  //Loop over the boundaries in the bulk mesh
306  for(unsigned i=0;i<n_bound;i++)
307  {
308  // Only update Dirichlet nodes
309  if (i!=1)
310  {
311  // How many nodes are there on this boundary?
312  unsigned n_node = Bulk_mesh_pt->nboundary_node(i);
313 
314  // Loop over the nodes on boundary
315  for (unsigned n=0;n<n_node;n++)
316  {
317  // Get pointer to node
318  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(i,n);
319 
320  // Extract nodal coordinates from node:
321  Vector<double> x(2);
322  x[0]=nod_pt->x(0);
323  x[1]=nod_pt->x(1);
324 
325  // Compute the value of the exact solution at the nodal point
326  Vector<double> u(1);
328 
329  // Assign the value to the one (and only) nodal value at this node
330  nod_pt->set_value(0,u[0]);
331  }
332  }
333  }
334 } // end of actions before solve
335 
336 
337 
338 //=====================start_of_actions_before_adapt======================
339 /// Actions before adapt: Wipe the mesh of prescribed flux elements
340 //========================================================================
341 template<class ELEMENT>
343 {
344  // Kill the flux elements and wipe surface mesh
345  delete_flux_elements(Surface_mesh_pt);
346 
347  // Rebuild the Problem's global mesh from its various sub-meshes
348  rebuild_global_mesh();
349 
350 }// end of actions_before_adapt
351 
352 
353 
354 //=====================start_of_actions_after_adapt=======================
355 /// Actions after adapt: Rebuild the mesh of prescribed flux elements
356 //========================================================================
357 template<class ELEMENT>
359 {
360  // Create prescribed-flux elements from all elements that are
361  // adjacent to boundary 1 and add them to surfac mesh
362  create_flux_elements(1,Bulk_mesh_pt,Surface_mesh_pt);
363 
364  // Rebuild the Problem's global mesh from its various sub-meshes
365  rebuild_global_mesh();
366 
367  // Set pointer to prescribed flux function for flux elements
368  set_prescribed_flux_pt();
369 
370  // Doc refinement levels in bulk mesh
371  unsigned min_refinement_level;
372  unsigned max_refinement_level;
373  Bulk_mesh_pt->get_refinement_levels(min_refinement_level,
374  max_refinement_level);
375  cout << "Min/max. refinement levels in bulk mesh: "
376  << min_refinement_level << " "
377  << max_refinement_level << std::endl;
378 
379 }// end of actions_after_adapt
380 
381 
382 
383 //==================start_of_set_prescribed_flux_pt=======================
384 /// Set pointer to prescribed-flux function for all
385 /// elements in the surface mesh
386 //========================================================================
387 template<class ELEMENT>
389 {
390  // Loop over the flux elements to pass pointer to prescribed flux function
391  unsigned n_element=Surface_mesh_pt->nelement();
392  for(unsigned e=0;e<n_element;e++)
393  {
394  // Upcast from GeneralisedElement to Poisson flux element
395  PoissonFluxElement<ELEMENT> *el_pt =
396  dynamic_cast< PoissonFluxElement<ELEMENT>*>(
397  Surface_mesh_pt->element_pt(e));
398 
399  // Set the pointer to the prescribed flux function
400  el_pt->flux_fct_pt() =
402  }
403 
404 }// end of set prescribed flux pt
405 
406 
407 
408 
409 //=====================start_of_doc=======================================
410 /// Doc the solution: doc_info contains labels/output directory etc.
411 //========================================================================
412 template<class ELEMENT>
414 {
415 
416  // Doc refinement levels in bulk mesh
417  unsigned min_refinement_level;
418  unsigned max_refinement_level;
419  Bulk_mesh_pt->get_refinement_levels(min_refinement_level,
420  max_refinement_level);
421  cout << "Ultimate min/max. refinement levels in bulk mesh : "
422  << min_refinement_level << " "
423  << max_refinement_level << std::endl;
424 
425 
426  ofstream some_file;
427  char filename[100];
428 
429  // Number of plot points
430  unsigned npts;
431  npts=5;
432 
433  // Output solution
434  //-----------------
435  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
436  doc_info.number());
437  some_file.open(filename);
438  Bulk_mesh_pt->output(some_file,npts);
439  some_file.close();
440 
441  // Output exact solution
442  //----------------------
443  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
444  doc_info.number());
445  some_file.open(filename);
446  Bulk_mesh_pt->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
447  some_file.close();
448 
449 
450  // Doc error and return of the square of the L2 error
451  //---------------------------------------------------
452  double error,norm;
453  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
454  doc_info.number());
455  some_file.open(filename);
456  Bulk_mesh_pt->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
457  error,norm);
458  some_file.close();
459 
460  // Doc L2 error and norm of solution
461  cout << "\nNorm of error : " << sqrt(error) << std::endl;
462  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
463 
464 
465 } // end of doc
466 
467 
468 //============start_of_create_flux_elements==============================
469 /// Create Poisson Flux Elements on the b-th boundary of the Mesh object
470 /// pointed to by bulk_mesh_pt and add the elements to the Mesh object
471 /// pointeed to by surface_mesh_pt.
472 //=======================================================================
473 template<class ELEMENT>
475 create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
476  Mesh* const &surface_mesh_pt)
477 {
478  // How many bulk elements are adjacent to boundary b?
479  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
480 
481  // Loop over the bulk elements adjacent to boundary b?
482  for(unsigned e=0;e<n_element;e++)
483  {
484  // Get pointer to the bulk element that is adjacent to boundary b
485  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
486  bulk_mesh_pt->boundary_element_pt(b,e));
487 
488  //Find the index of the face of element e along boundary b
489  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
490 
491  // Build the corresponding prescribed-flux element
492  PoissonFluxElement<ELEMENT>* flux_element_pt = new
493  PoissonFluxElement<ELEMENT>(bulk_elem_pt,face_index);
494 
495  //Add the prescribed-flux element to the surface mesh
496  surface_mesh_pt->add_element_pt(flux_element_pt);
497 
498  } //end of loop over bulk elements adjacent to boundary b
499 
500 } // end of create_flux_elements
501 
502 
503 
504 //============start_of_delete_flux_elements==============================
505 /// Delete Poisson Flux Elements and wipe the surface mesh
506 //=======================================================================
507 template<class ELEMENT>
509 delete_flux_elements(Mesh* const &surface_mesh_pt)
510 {
511  // How many surface elements are in the surface mesh
512  unsigned n_element = surface_mesh_pt->nelement();
513 
514  // Loop over the surface elements
515  for(unsigned e=0;e<n_element;e++)
516  {
517  // Kill surface element
518  delete surface_mesh_pt->element_pt(e);
519  }
520 
521  // Wipe the mesh
522  surface_mesh_pt->flush_element_and_node_storage();
523 
524 } // end of delete_flux_elements
525 
526 
527 
528 //==========start_of_main=================================================
529 /// Demonstrate how to solve 2D Poisson problem with flux boundary
530 /// conditions, using two meshes.
531 //========================================================================
532 int main()
533 {
534 
535  //Set up the problem
536  //------------------
537 
538  //Set up the problem with 2D nine-node elements from the
539  //QPoissonElement family. Pass pointer to source function.
542 
543 
544  // Create label for output
545  //------------------------
546  DocInfo doc_info;
547 
548  // Set output directory
549  doc_info.set_directory("RESLT");
550 
551  // Step number
552  doc_info.number()=0;
553 
554 
555 
556  // Check if we're ready to go:
557  //----------------------------
558  cout << "\n\n\nProblem self-test ";
559  if (problem.self_test()==0)
560  {
561  cout << "passed: Problem can be solved." << std::endl;
562  }
563  else
564  {
565  throw OomphLibError("Self test failed",
566  OOMPH_CURRENT_FUNCTION,
567  OOMPH_EXCEPTION_LOCATION);
568  }
569 
570 
571  // Set the orientation of the "step" to 45 degrees
573 
574  // Initial value for the steepness of the "step"
576 
577  // Do a couple of solutions for different forcing functions
578  //---------------------------------------------------------
579  unsigned nstep=4;
580  for (unsigned istep=0;istep<nstep;istep++)
581  {
582  // Increase the steepness of the step:
584 
585  cout << "\n\nSolving for TanhSolnForPoisson::Alpha="
586  << TanhSolnForPoisson::Alpha << std::endl << std::endl;
587 
588  // Solve the problem
589  problem.newton_solve(3);
590 
591  //Output solution
592  problem.doc_solution(doc_info);
593 
594  //Increment counter for solutions
595  doc_info.number()++;
596  }
597 
598 } //end of main
599 
600 
601 
602 
603 
604 
605 
606 
607 
SimpleRefineableRectangularQuadMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
virtual ~SimpleRefineableRectangularQuadMesh()
Destructor: Empty.
SimpleRefineableRectangularQuadMesh(const unsigned &Nx, const unsigned &Ny, const double &Lx, const double &Ly, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Pass number of elements in the horizontal and vertical directions, and the corresponding dimensions...
void create_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &surface_mesh_pt)
Create Poisson flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to the...
double Alpha
Parameter for steepness of "step".
void doc_solution(DocInfo &doc_info)
Doc the solution: doc_info contains labels/output directory etc.
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the exact solutio...
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void source_function(const Vector< double > &x, double &source)
Source function required to make the solution above an exact solution.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed flux elements.
double TanPhi
Parameter for angle Phi of "step".
void set_prescribed_flux_pt()
Set pointer to prescribed-flux function for all elements in the surface mesh.
void prescribed_flux_on_fixed_x_boundary(const Vector< double > &x, double &flux)
Flux required by the exact solution on a boundary on which x is fixed.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed flux elements.
PoissonEquations< 2 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
RefineableTwoMeshFluxPoissonProblem(PoissonEquations< 2 >::PoissonSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function.
void delete_flux_elements(Mesh *const &surface_mesh_pt)
Delete Poisson flux elements and wipe the surface mesh.
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.