two_d_poisson_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 adaptive mesh refinement
31 
32 //Generic routines
33 #include "generic.h"
34 
35 // The Poisson equations
36 #include "poisson.h"
37 
38 // The mesh
39 #include "meshes/simple_rectangular_quadmesh.h"
40 
41 // We need to include the refineable_quad_mesh's
42 // templated source here to allow the build of
43 // all templates.
44 //#include "generic/refineable_quad_mesh.template.cc"
45 
46 using namespace std;
47 
48 using namespace oomph;
49 
50 //==============================start_of_mesh======================
51 /// Refineable equivalent of the SimpleRectangularQuadMesh.
52 /// Refinement is performed by the QuadTree-based procedures
53 /// implemented in the RefineableQuadMesh base class.
54 //=================================================================
55 template<class ELEMENT>
57  public virtual SimpleRectangularQuadMesh<ELEMENT>,
58  public RefineableQuadMesh<ELEMENT>
59 {
60 
61 public:
62 
63  /// \short Pass number of elements in the horizontal
64  /// and vertical directions, and the corresponding dimensions.
65  /// Timestepper defaults to Static.
67  const unsigned &Ny,
68  const double &Lx, const double &Ly,
69  TimeStepper* time_stepper_pt=
70  &Mesh::Default_TimeStepper) :
71  SimpleRectangularQuadMesh<ELEMENT>(Nx,Ny,Lx,Ly,time_stepper_pt)
72  {
73  // Nodal positions etc. were created in constructor for
74  // SimpleRectangularQuadMesh<...> --> We only need to set up
75  // adaptivity information: Associate finite elements with their
76  // QuadTrees and plant them in a QuadTreeForest:
77  this->setup_quadtree_forest();
78 
79  } // end of constructor
80 
81 
82  /// Destructor: Empty
84 
85 }; // end of mesh
86 
87 
88 
89 //===== start_of_namespace=============================================
90 /// Namespace for exact solution for Poisson equation with "sharp step"
91 //=====================================================================
92 namespace TanhSolnForPoisson
93 {
94 
95  /// Parameter for steepness of "step"
96  double Alpha=5.0;
97 
98  /// Parameter for angle Phi of "step"
99  double TanPhi=0.0;
100 
101  /// Exact solution as a Vector
102  void get_exact_u(const Vector<double>& x, Vector<double>& u)
103  {
104  u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
105  }
106 
107  /// Source function required to make the solution above an exact solution
108  void get_source(const Vector<double>& x, double& source)
109  {
110  source = 2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
111  (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*
112  Alpha*Alpha*TanPhi*TanPhi+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))*Alpha*Alpha;
114  }
115 
116 } // end of namespace
117 
118 
119 
120 
121 
122 
123 
124 
125 //====== start_of_problem_class=======================================
126 /// 2D Poisson problem on rectangular domain, discretised with
127 /// refineable 2D QPoisson elements. The specific type of element is
128 /// specified via the template parameter.
129 //====================================================================
130 template<class ELEMENT>
131 class RefineablePoissonProblem : public Problem
132 {
133 
134 public:
135 
136  /// Constructor: Pass pointer to source function
137  RefineablePoissonProblem(PoissonEquations<2>::PoissonSourceFctPt
138  source_fct_pt);
139 
140  /// Destructor (empty)
142 
143  /// \short Update the problem specs before solve: Reset boundary conditions
144  /// to the values from the exact solution.
145  void actions_before_newton_solve();
146 
147  /// Update the problem after solve (empty)
149 
150  /// \short Doc the solution. DocInfo object stores flags/labels for where the
151  /// output gets written to
152  void doc_solution(DocInfo& doc_info);
153 
154  /// \short Overloaded version of the Problem's access function to
155  /// the mesh. Recasts the pointer to the base Mesh object to
156  /// the actual mesh type.
158  {
159  return dynamic_cast<SimpleRefineableRectangularQuadMesh<ELEMENT>*>(
160  Problem::mesh_pt());
161  }
162 
163 private:
164 
165  /// Pointer to source function
166  PoissonEquations<2>::PoissonSourceFctPt Source_fct_pt;
167 
168 }; // end of problem class
169 
170 
171 
172 
173 //=====start_of_constructor===============================================
174 /// Constructor for Poisson problem: Pass pointer to source function.
175 //========================================================================
176 template<class ELEMENT>
178  RefineablePoissonProblem(PoissonEquations<2>::PoissonSourceFctPt
179  source_fct_pt)
180  : Source_fct_pt(source_fct_pt)
181 {
182 
183  // Setup mesh
184 
185  // # of elements in x-direction
186  unsigned n_x=4;
187 
188  // # of elements in y-direction
189  unsigned n_y=4;
190 
191  // Domain length in x-direction
192  double l_x=1.0;
193 
194  // Domain length in y-direction
195  double l_y=2.0;
196 
197  // Build and assign mesh
198  Problem::mesh_pt() =
200 
201  // Create/set error estimator
202  mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
203 
204  // Set the boundary conditions for this problem: All nodes are
205  // free by default -- only need to pin the ones that have Dirichlet conditions
206  // here.
207  unsigned num_bound = mesh_pt()->nboundary();
208  for(unsigned ibound=0;ibound<num_bound;ibound++)
209  {
210  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
211  for (unsigned inod=0;inod<num_nod;inod++)
212  {
213  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
214  }
215  }
216 
217  // Complete the build of all elements so they are fully functional
218 
219  // Loop over the elements to set up element-specific
220  // things that cannot be handled by the (argument-free!) ELEMENT
221  // constructor: Pass pointer to source function
222  unsigned n_element = mesh_pt()->nelement();
223  for(unsigned i=0;i<n_element;i++)
224  {
225  // Upcast from GeneralsedElement to the present element
226  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
227 
228  //Set the source function pointer
229  el_pt->source_fct_pt() = Source_fct_pt;
230  }
231 
232  // Setup equation numbering scheme
233  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
234 
235 } // end of constructor
236 
237 
238 
239 
240 //========================================start_of_actions_before_newton_solve===
241 /// Update the problem specs before solve: (Re-)set boundary conditions
242 /// to the values from the exact solution.
243 //========================================================================
244 template<class ELEMENT>
246 {
247  // How many boundaries are there?
248  unsigned num_bound = mesh_pt()->nboundary();
249 
250  //Loop over the boundaries
251  for(unsigned ibound=0;ibound<num_bound;ibound++)
252  {
253  // How many nodes are there on this boundary?
254  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
255 
256  // Loop over the nodes on boundary
257  for (unsigned inod=0;inod<num_nod;inod++)
258  {
259  // Get pointer to node
260  Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
261 
262  // Extract nodal coordinates from node:
263  Vector<double> x(2);
264  x[0]=nod_pt->x(0);
265  x[1]=nod_pt->x(1);
266 
267  // Compute the value of the exact solution at the nodal point
268  Vector<double> u(1);
270 
271  // Assign the value to the one (and only) nodal value at this node
272  nod_pt->set_value(0,u[0]);
273  }
274  }
275 } // end of actions before solve
276 
277 
278 
279 //===============start_of_doc=============================================
280 /// Doc the solution: doc_info contains labels/output directory etc.
281 //========================================================================
282 template<class ELEMENT>
284 {
285 
286  ofstream some_file;
287  char filename[100];
288 
289  // Number of plot points: npts x npts
290  unsigned npts=5;
291 
292  // Output solution
293  //-----------------
294  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
295  doc_info.number());
296  some_file.open(filename);
297  mesh_pt()->output(some_file,npts);
298  some_file.close();
299 
300 
301  // Output exact solution
302  //----------------------
303  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
304  doc_info.number());
305  some_file.open(filename);
306  mesh_pt()->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
307  some_file.close();
308 
309  // Doc error and return of the square of the L2 error
310  //---------------------------------------------------
311  double error,norm;
312  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
313  doc_info.number());
314  some_file.open(filename);
315  mesh_pt()->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
316  error,norm);
317  some_file.close();
318 
319  // Doc L2 error and norm of solution
320  cout << "\nNorm of error : " << sqrt(error) << std::endl;
321  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
322 
323 } // end of doc
324 
325 
326 
327 
328 
329 
330 //===== start_of_main=====================================================
331 /// Driver code for 2D Poisson problem
332 //========================================================================
333 int main()
334 {
335 
336  //Set up the problem
337  //------------------
338 
339  // Create the problem with 2D nine-node refineable elements from the
340  // RefineableQuadPoissonElement family. Pass pointer to source function.
343 
344  // Create label for output
345  //------------------------
346  DocInfo doc_info;
347 
348  // Set output directory
349  doc_info.set_directory("RESLT");
350 
351  // Step number
352  doc_info.number()=0;
353 
354  // Check if we're ready to go:
355  //----------------------------
356  cout << "\n\n\nProblem self-test ";
357  if (problem.self_test()==0)
358  {
359  cout << "passed: Problem can be solved." << std::endl;
360  }
361  else
362  {
363  throw OomphLibError("Self test failed",
364  OOMPH_CURRENT_FUNCTION,
365  OOMPH_EXCEPTION_LOCATION);
366  }
367 
368 
369  // Set the orientation of the "step" to 45 degrees
371 
372  // Choose a large value for the steepness of the "step"
374 
375  // Solve the problem, performing up to 4 adaptive refinements
376  problem.newton_solve(4);
377 
378  //Output the solution
379  problem.doc_solution(doc_info);
380 
381 } //end of main
382 
383 
384 
385 
386 
387 
388 
389 
390 
SimpleRefineableRectangularQuadMesh< ELEMENT > * mesh_pt()
Overloaded version of the Problem's access function to the mesh. Recasts the pointer to the base Mesh...
int main()
Driver code for 2D Poisson problem.
void actions_after_newton_solve()
Update the problem after solve (empty)
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...
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...
double Alpha
Parameter for steepness of "step".
~RefineablePoissonProblem()
Destructor (empty)
double TanPhi
Parameter for angle Phi of "step".
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the exact solutio...
void get_source(const Vector< double > &x, double &source)
Source function required to make the solution above an exact solution.
RefineablePoissonProblem(PoissonEquations< 2 >::PoissonSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function.