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