two_d_poisson_flux_bc.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
28//Generic routines
29#include "generic.h"
30
31// The Poisson equations
32#include "poisson.h"
33
34// The mesh
35#include "meshes/simple_rectangular_quadmesh.h"
36
37using namespace std;
38
39using namespace oomph;
40
41//===== start_of_namespace=============================================
42/// Namespace for exact solution for Poisson equation with "sharp step"
43//=====================================================================
45{
46
47 /// Parameter for steepness of "step"
48 double Alpha=1.0;
49
50 /// Parameter for angle Phi of "step"
51 double TanPhi=0.0;
52
53 /// Exact solution as a Vector
54 void get_exact_u(const Vector<double>& x, Vector<double>& u)
55 {
56 u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
57 }
58
59 /// Source function required to make the solution above an exact solution
60 void source_function(const Vector<double>& x, double& source)
61 {
62 source = 2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
63 (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*
64 Alpha*Alpha*TanPhi*TanPhi+2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
65 (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*Alpha;
66 }
67
68 /// Flux required by the exact solution on a boundary on which x is fixed
69 void prescribed_flux_on_fixed_x_boundary(const Vector<double>& x,
70 double& flux)
71 {
72 //The outer unit normal to the boundary is (1,0)
73 double N[2] = {1.0, 0.0};
74 //The flux in terms of the normal is
75 flux =
76 -(1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*TanPhi*N[0]+(
77 1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*N[1];
78 }
79
80} // end of namespace
81
82
83
84
85
86//========= start_of_problem_class=====================================
87/// 2D Poisson problem on rectangular domain, discretised with
88/// 2D QPoisson elements. Flux boundary conditions are applied
89/// along boundary 1 (the boundary where x=L). The specific type of
90/// element is specified via the template parameter.
91//====================================================================
92template<class ELEMENT>
93class FluxPoissonProblem : public Problem
94{
95
96public:
97
98 /// Constructor: Pass pointer to source function
99 FluxPoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt);
100
101 /// Destructor (empty)
103
104 /// Doc the solution. DocInfo object stores flags/labels for where the
105 /// output gets written to
106 void doc_solution(DocInfo& doc_info);
107
108private:
109
110 /// Update the problem specs before solve: Reset boundary conditions
111 /// to the values from the exact solution.
113
114 /// Update the problem specs after solve (empty)
116
117 /// Create Poisson flux elements on the b-th boundary of the
118 /// problem's mesh
119 void create_flux_elements(const unsigned &b);
120
121 /// Number of Poisson "bulk" elements (We're attaching the flux
122 /// elements to the bulk mesh --> only the first Npoisson_elements elements
123 /// in the mesh are bulk elements!)
125
126 /// Pointer to source function
127 PoissonEquations<2>::PoissonSourceFctPt Source_fct_pt;
128
129}; // end of problem class
130
131
132
133//=======start_of_constructor=============================================
134/// Constructor for Poisson problem: Pass pointer to source function.
135//========================================================================
136template<class ELEMENT>
138 FluxPoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt)
139 : Source_fct_pt(source_fct_pt)
140{
141
142 // Setup mesh
143
144 // # of elements in x-direction
145 unsigned n_x=4;
146
147 // # of elements in y-direction
148 unsigned n_y=4;
149
150 // Domain length in x-direction
151 double l_x=1.0;
152
153 // Domain length in y-direction
154 double l_y=2.0;
155
156 // Build and assign mesh
157 Problem::mesh_pt()=new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
158
159 // Store number of Poisson bulk elements (= number of elements so far).
160 Npoisson_elements=mesh_pt()->nelement();
161
162 // Create prescribed-flux elements from all elements that are
163 // adjacent to boundary 1 and add them to the (single) global mesh
165
166 // Set the boundary conditions for this problem: All nodes are
167 // free by default -- just pin the ones that have Dirichlet conditions
168 // here.
169 unsigned n_bound = mesh_pt()->nboundary();
170 for(unsigned b=0;b<n_bound;b++)
171 {
172 //Leave nodes on boundary 1 free
173 if(b!=1)
174 {
175 unsigned n_node= mesh_pt()->nboundary_node(b);
176 for (unsigned n=0;n<n_node;n++)
177 {
178 mesh_pt()->boundary_node_pt(b,n)->pin(0);
179 }
180 }
181 }
182
183 // Loop over the Poisson bulk elements to set up element-specific
184 // things that cannot be handled by constructor: Pass pointer to source
185 // function
186 for(unsigned e=0;e<Npoisson_elements;e++)
187 {
188 // Upcast from GeneralisedElement to Poisson bulk element
189 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
190
191 //Set the source function pointer
192 el_pt->source_fct_pt() = Source_fct_pt;
193 }
194
195 // Total number of elements:
196 unsigned n_element=mesh_pt()->nelement();
197
198 // Loop over the flux elements (located at the "end" of the
199 // mesh) to pass function pointer to prescribed-flux function.
200 for (unsigned e=Npoisson_elements;e<n_element;e++)
201 {
202 // Upcast from GeneralisedElement to Poisson flux element
203 PoissonFluxElement<ELEMENT> *el_pt =
204 dynamic_cast<PoissonFluxElement<ELEMENT>*>(mesh_pt()->element_pt(e));
205
206 // Set the pointer to the prescribed flux function
207 el_pt->flux_fct_pt() =
209 }
210
211 // Setup equation numbering scheme
212 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
213
214} // end of constructor
215
216
217//============start_of_create_flux_elements==============================
218/// Create Poisson Flux Elements on the b-th boundary of the Mesh
219//=======================================================================
220template<class ELEMENT>
222{
223 // How many bulk elements are adjacent to boundary b?
224 unsigned n_element = mesh_pt()->nboundary_element(b);
225
226 // Loop over the bulk elements adjacent to boundary b?
227 for(unsigned e=0;e<n_element;e++)
228 {
229 // Get pointer to the bulk element that is adjacent to boundary b
230 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
231 mesh_pt()->boundary_element_pt(b,e));
232
233 // What is the index of the face of the bulk element at the boundary
234 int face_index = mesh_pt()->face_index_at_boundary(b,e);
235
236 // Build the corresponding prescribed-flux element
237 PoissonFluxElement<ELEMENT>* flux_element_pt = new
238 PoissonFluxElement<ELEMENT>(bulk_elem_pt,face_index);
239
240 //Add the prescribed-flux element to the mesh
241 mesh_pt()->add_element_pt(flux_element_pt);
242
243 } //end of loop over bulk elements adjacent to boundary b
244
245} // end of create_flux_elements
246
247
248
249//====================start_of_actions_before_newton_solve================
250/// Update the problem specs before solve: Reset boundary conditions
251/// to the values from the exact solution.
252//========================================================================
253template<class ELEMENT>
255{
256 // How many boundaries are there?
257 unsigned n_bound = mesh_pt()->nboundary();
258
259 //Loop over the boundaries
260 for(unsigned i=0;i<n_bound;i++)
261 {
262 // Only update Dirichlet nodes
263 if (i!=1)
264 {
265 // How many nodes are there on this boundary?
266 unsigned n_node = mesh_pt()->nboundary_node(i);
267
268 // Loop over the nodes on boundary
269 for (unsigned n=0;n<n_node;n++)
270 {
271 // Get pointer to node
272 Node* nod_pt = mesh_pt()->boundary_node_pt(i,n);
273
274 // Extract nodal coordinates from node:
275 Vector<double> x(2);
276 x[0]=nod_pt->x(0);
277 x[1]=nod_pt->x(1);
278
279 // Compute the value of the exact solution at the nodal point
280 Vector<double> u(1);
282
283 // Assign the value to the one (and only) nodal value at this node
284 nod_pt->set_value(0,u[0]);
285 }
286 }
287 }
288} // end of actions before solve
289
290
291
292//=====================start_of_doc=======================================
293/// Doc the solution: doc_info contains labels/output directory etc.
294//========================================================================
295template<class ELEMENT>
297{
298
299 ofstream some_file;
300 char filename[100];
301
302 // Number of plot points
303 unsigned npts;
304 npts=5;
305
306 // Output solution
307 //-----------------
308 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
309 doc_info.number());
310 some_file.open(filename);
311 mesh_pt()->output(some_file,npts);
312 some_file.close();
313
314 // Output exact solution
315 //----------------------
316 sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
317 doc_info.number());
318 some_file.open(filename);
319 for(unsigned e=0;e<Npoisson_elements;e++)
320 {
321 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
322 el_pt->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
323 }
324 some_file.close();
325
326 // Can't use black-box error computatation routines because
327 // the mesh contains two different types of elements.
328 // error function hasn't been implemented for the prescribed
329 // flux elements...
330
331} // end of doc
332
333
334
335//==========start_of_main=================================================
336/// Demonstrate how to solve 2D Poisson problem with flux boundary
337/// conditions
338//========================================================================
339int main()
340{
341
342 //Set up the problem
343 //------------------
344
345 //Set up the problem with 2D nine-node elements from the
346 //QPoissonElement family. Pass pointer to source function.
349
350
351 // Create label for output
352 //------------------------
353 DocInfo doc_info;
354
355 // Set output directory
356 doc_info.set_directory("RESLT");
357
358 // Step number
359 doc_info.number()=0;
360
361
362
363 // Check that we're ready to go:
364 //----------------------------
365 cout << "\n\n\nProblem self-test ";
366 if (problem.self_test()==0)
367 {
368 cout << "passed: Problem can be solved." << std::endl;
369 }
370 else
371 {
372 throw OomphLibError("Self test failed",
373 OOMPH_CURRENT_FUNCTION,
374 OOMPH_EXCEPTION_LOCATION);
375 }
376
377
378 // Set the orientation of the "step" to 45 degrees
380
381 // Initial value for the steepness of the "step"
383
384 // Do a couple of solutions for different forcing functions
385 //---------------------------------------------------------
386 unsigned nstep=4;
387 for (unsigned istep=0;istep<nstep;istep++)
388 {
389 // Increase the steepness of the step:
391
392 cout << "\n\nSolving for TanhSolnForPoisson::Alpha="
393 << TanhSolnForPoisson::Alpha << std::endl << std::endl;
394
395 // Solve the problem
396 problem.newton_solve();
397
398 //Output solution
399 problem.doc_solution(doc_info);
400
401 //Increment counter for solutions
402 doc_info.number()++;
403 }
404
405} //end of main
406
407
408
409
410
411
412
413
414
2D Poisson problem on rectangular domain, discretised with 2D QPoisson elements. Flux boundary condit...
FluxPoissonProblem(PoissonEquations< 2 >::PoissonSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function.
unsigned Npoisson_elements
Number of Poisson "bulk" elements (We're attaching the flux elements to the bulk mesh --> only the fi...
void create_flux_elements(const unsigned &b)
Create Poisson flux elements on the b-th boundary of the problem's mesh.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the exact solutio...
PoissonEquations< 2 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to.
~FluxPoissonProblem()
Destructor (empty)
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.