two_d_poisson.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
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
42//===== start_of_namespace=============================================
43/// Namespace for exact solution for Poisson equation with "sharp step"
44//=====================================================================
46{
47
48 /// Parameter for steepness of "step"
49 double Alpha=1.0;
50
51 /// Parameter for angle Phi of "step"
52 double TanPhi=0.0;
53
54 /// Exact solution as a Vector
55 void get_exact_u(const Vector<double>& x, Vector<double>& u)
56 {
57 u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
58 }
59
60 /// Source function required to make the solution above an exact solution
61 void source_function(const Vector<double>& x, double& source)
62 {
63 source = 2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
64 (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*
65 Alpha*Alpha*TanPhi*TanPhi+2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
66 (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*Alpha;
67 }
68
69} // end of namespace
70
71
72
73
74
75
76
77
78//====== start_of_problem_class=======================================
79/// 2D Poisson problem on rectangular domain, discretised with
80/// 2D QPoisson elements. The specific type of element is
81/// specified via the template parameter.
82//====================================================================
83template<class ELEMENT>
84class PoissonProblem : public Problem
85{
86
87public:
88
89 /// Constructor: Pass pointer to source function
90 PoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt);
91
92 /// Destructor (empty)
94
95 /// Update the problem specs before solve: Reset boundary conditions
96 /// to the values from the exact solution.
98
99 /// Update the problem after solve (empty)
101
102 /// Doc the solution. DocInfo object stores flags/labels for where the
103 /// output gets written to
104 void doc_solution(DocInfo& doc_info);
105
106private:
107
108 /// Pointer to source function
109 PoissonEquations<2>::PoissonSourceFctPt Source_fct_pt;
110
111}; // end of problem class
112
113
114
115
116//=====start_of_constructor===============================================
117/// Constructor for Poisson problem: Pass pointer to source function.
118//========================================================================
119template<class ELEMENT>
121 PoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt)
122 : Source_fct_pt(source_fct_pt)
123{
124 // Setup mesh
125
126 // # of elements in x-direction
127 unsigned n_x=4;
128
129 // # of elements in y-direction
130 unsigned n_y=4;
131
132 // Domain length in x-direction
133 double l_x=1.0;
134
135 // Domain length in y-direction
136 double l_y=2.0;
137
138 // Build and assign mesh
139 Problem::mesh_pt() = new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
140
141 // Set the boundary conditions for this problem: All nodes are
142 // free by default -- only need to pin the ones that have Dirichlet conditions
143 // here.
144 unsigned n_bound = mesh_pt()->nboundary();
145 for(unsigned i=0;i<n_bound;i++)
146 {
147 unsigned n_node = mesh_pt()->nboundary_node(i);
148 for (unsigned n=0;n<n_node;n++)
149 {
150 mesh_pt()->boundary_node_pt(i,n)->pin(0);
151 }
152 }
153
154 // Complete the build of all elements so they are fully functional
155
156 // Loop over the elements to set up element-specific
157 // things that cannot be handled by the (argument-free!) ELEMENT
158 // constructor: Pass pointer to source function
159 unsigned n_element = mesh_pt()->nelement();
160 for(unsigned i=0;i<n_element;i++)
161 {
162 // Upcast from GeneralsedElement to the present element
163 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
164
165 //Set the source function pointer
166 el_pt->source_fct_pt() = Source_fct_pt;
167 }
168
169
170 // Setup equation numbering scheme
171 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
172
173} // end of constructor
174
175
176
177
178//========================================start_of_actions_before_newton_solve===
179/// Update the problem specs before solve: (Re-)set boundary conditions
180/// to the values from the exact solution.
181//========================================================================
182template<class ELEMENT>
184{
185 // How many boundaries are there?
186 unsigned n_bound = mesh_pt()->nboundary();
187
188 //Loop over the boundaries
189 for(unsigned i=0;i<n_bound;i++)
190 {
191 // How many nodes are there on this boundary?
192 unsigned n_node = mesh_pt()->nboundary_node(i);
193
194 // Loop over the nodes on boundary
195 for (unsigned n=0;n<n_node;n++)
196 {
197 // Get pointer to node
198 Node* nod_pt=mesh_pt()->boundary_node_pt(i,n);
199
200 // Extract nodal coordinates from node:
201 Vector<double> x(2);
202 x[0]=nod_pt->x(0);
203 x[1]=nod_pt->x(1);
204
205 // Compute the value of the exact solution at the nodal point
206 Vector<double> u(1);
208
209 // Assign the value to the one (and only) nodal value at this node
210 nod_pt->set_value(0,u[0]);
211 }
212 }
213} // end of actions before solve
214
215
216
217//===============start_of_doc=============================================
218/// Doc the solution: doc_info contains labels/output directory etc.
219//========================================================================
220template<class ELEMENT>
222{
223
224 ofstream some_file;
225 char filename[100];
226
227 // Number of plot points: npts x npts
228 unsigned npts=5;
229
230 // Output solution
231 //-----------------
232 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
233 doc_info.number());
234 some_file.open(filename);
235 mesh_pt()->output(some_file,npts);
236 some_file.close();
237
238
239 // Output exact solution
240 //----------------------
241 sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
242 doc_info.number());
243 some_file.open(filename);
244 mesh_pt()->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
245 some_file.close();
246
247 // Doc error and return of the square of the L2 error
248 //---------------------------------------------------
249 double error,norm;
250 sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
251 doc_info.number());
252 some_file.open(filename);
253 mesh_pt()->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
254 error,norm);
255 some_file.close();
256
257 // Doc L2 error and norm of solution
258 cout << "\nNorm of error : " << sqrt(error) << std::endl;
259 cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
260
261} // end of doc
262
263
264
265
266
267
268//===== start_of_main=====================================================
269/// Driver code for 2D Poisson problem
270//========================================================================
271int main()
272{
273
274 //Set up the problem
275 //------------------
276
277 // Create the problem with 2D nine-node elements from the
278 // QPoissonElement family. Pass pointer to source function.
281
282 // Create label for output
283 //------------------------
284 DocInfo doc_info;
285
286 // Set output directory
287 doc_info.set_directory("RESLT");
288
289 // Step number
290 doc_info.number()=0;
291
292 // Check that we're ready to go:
293 //----------------------------
294 cout << "\n\n\nProblem self-test ";
295 if (problem.self_test()==0)
296 {
297 cout << "passed: Problem can be solved." << std::endl;
298 }
299 else
300 {
301 throw OomphLibError("Self test failed",
302 OOMPH_CURRENT_FUNCTION,
303 OOMPH_EXCEPTION_LOCATION);
304 }
305
306
307 // Set the orientation of the "step" to 45 degrees
309
310 // Initial value for the steepness of the "step"
312
313
314 // Do a couple of solutions for different forcing functions
315 //---------------------------------------------------------
316 unsigned nstep=4;
317 for (unsigned istep=0;istep<nstep;istep++)
318 {
319
320 // Increase the steepness of the step:
322
323 cout << "\n\nSolving for TanhSolnForPoisson::Alpha="
324 << TanhSolnForPoisson::Alpha << std::endl << std::endl;
325
326 // Solve the problem
327 problem.newton_solve();
328
329 //Output the solution
330 problem.doc_solution(doc_info);
331
332 //Increment counter for solutions
333 doc_info.number()++;
334
335 }
336
337
338} //end of main
339
340
341
342
343
344
345
346
347
2D Poisson problem on rectangular domain, discretised with 2D QPoisson elements. The specific type of...
PoissonEquations< 2 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
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 after solve (empty)
PoissonProblem(PoissonEquations< 2 >::PoissonSourceFctPt source_fct_pt)
Constructor: Pass 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.
~PoissonProblem()
Destructor (empty)
Namespace for exact solution for Poisson equation with "sharp step".
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()
Driver code for 2D Poisson problem.