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