eighth_sphere_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// 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 3D poisson problem
31 
32 //Generic routines
33 #include "generic.h"
34 
35 // The Poisson equations
36 #include "poisson.h"
37 
38 // The mesh
39 #include "meshes/eighth_sphere_mesh.h"
40 
41 using namespace std;
42 
43 using namespace oomph;
44 
45 //=============start_of_namespace=====================================
46 /// Namespace for exact solution for Poisson equation with sharp step
47 //====================================================================
48 namespace TanhSolnForPoisson
49 {
50 
51  /// Parameter for steepness of step
52  double Alpha=1;
53 
54  ///\short Orientation (non-normalised x-component of unit vector in direction
55  /// of step plane)
56  double N_x=-1.0;
57 
58  ///\short Orientation (non-normalised y-component of unit vector in direction
59  /// of step plane)
60  double N_y=-1.0;
61 
62  ///\short Orientation (non-normalised z-component of unit vector in direction
63  /// of step plane)
64  double N_z=1.0;
65 
66 
67  ///\short Orientation (x-coordinate of step plane)
68  double X_0=0.0;
69 
70  ///\short Orientation (y-coordinate of step plane)
71  double Y_0=0.0;
72 
73  ///\short Orientation (z-coordinate of step plane)
74  double Z_0=0.0;
75 
76 
77  // Exact solution as a Vector
78  void get_exact_u(const Vector<double>& x, Vector<double>& u)
79  {
80  u[0] = tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-Y_0)*
81  N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*
82  N_z/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)));
83  }
84 
85  /// Exact solution as a scalar
86  void get_exact_u(const Vector<double>& x, double& u)
87  {
88  u = tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-Y_0)*
89  N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*
90  N_z/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)));
91  }
92 
93 
94  /// Source function to make it an exact solution
95  void get_source(const Vector<double>& x, double& source)
96  {
97 
98  double s1,s2,s3,s4;
99 
100  s1 = -2.0*tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-
101 Y_0)*N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*N_z/sqrt(N_x*N_x+N_y*N_y+N_z*
102 N_z)))*(1.0-pow(tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-
103 Y_0)*N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*N_z/sqrt(N_x*N_x+N_y*N_y+N_z*
104 N_z))),2.0))*Alpha*Alpha*N_x*N_x/(N_x*N_x+N_y*N_y+N_z*N_z);
105  s3 = -2.0*tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-
106 Y_0)*N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*N_z/sqrt(N_x*N_x+N_y*N_y+N_z*
107 N_z)))*(1.0-pow(tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-
108 Y_0)*N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*N_z/sqrt(N_x*N_x+N_y*N_y+N_z*
109 N_z))),2.0))*Alpha*Alpha*N_y*N_y/(N_x*N_x+N_y*N_y+N_z*N_z);
110  s4 = -2.0*tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-
111 Y_0)*N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*N_z/sqrt(N_x*N_x+N_y*N_y+N_z*
112 N_z)))*(1.0-pow(tanh(Alpha*((x[0]-X_0)*N_x/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[1]-
113 Y_0)*N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*N_z/sqrt(N_x*N_x+N_y*N_y+N_z*
114 N_z))),2.0))*Alpha*Alpha*N_z*N_z/(N_x*N_x+N_y*N_y+N_z*N_z);
115  s2 = s3+s4;
116  source = s1+s2;
117  }
118 
119 
120 } // end of namespace
121 
122 
123 
124 
125 ////////////////////////////////////////////////////////////////////////
126 ////////////////////////////////////////////////////////////////////////
127 ////////////////////////////////////////////////////////////////////////
128 
129 
130 
131 //=======start_of_class_definition====================================
132 /// Poisson problem in refineable eighth of a sphere mesh.
133 //====================================================================
134 template<class ELEMENT>
135 class EighthSpherePoissonProblem : public Problem
136 {
137 
138 public:
139 
140  /// Constructor: Pass pointer to source function
142  PoissonEquations<3>::PoissonSourceFctPt source_fct_pt);
143 
144  /// Destructor: Empty
146 
147  /// \short Overload generic access function by one that returns
148  /// a pointer to the specific mesh
149  RefineableEighthSphereMesh<ELEMENT>* mesh_pt()
150  {
151  return dynamic_cast<RefineableEighthSphereMesh<ELEMENT>*>(Problem::mesh_pt());
152  }
153 
154  /// Update the problem specs after solve (empty)
156 
157  /// \short Update the problem specs before solve:
158  /// Set Dirchlet boundary conditions from exact solution.
160  {
161  //Loop over the boundaries
162  unsigned num_bound = mesh_pt()->nboundary();
163  for(unsigned ibound=0;ibound<num_bound;ibound++)
164  {
165  // Loop over the nodes on boundary
166  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
167  for (unsigned inod=0;inod<num_nod;inod++)
168  {
169  Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
170  double u;
171  Vector<double> x(3);
172  x[0]=nod_pt->x(0);
173  x[1]=nod_pt->x(1);
174  x[2]=nod_pt->x(2);
176  nod_pt->set_value(0,u);
177  }
178  }
179  }
180 
181  /// Doc the solution
182  void doc_solution(DocInfo& doc_info);
183 
184 private:
185 
186  /// Pointer to source function
187  PoissonEquations<3>::PoissonSourceFctPt Source_fct_pt;
188 
189 }; // end of class definition
190 
191 
192 
193 
194 
195 //====================start_of_constructor================================
196 /// Constructor for Poisson problem on eighth of a sphere mesh
197 //========================================================================
198 template<class ELEMENT>
200  PoissonEquations<3>::PoissonSourceFctPt source_fct_pt) :
201  Source_fct_pt(source_fct_pt)
202 {
203 
204  // Setup parameters for exact tanh solution
205  // Steepness of step
207 
208  /// Create mesh for sphere of radius 5
209  double radius=5.0;
210  Problem::mesh_pt() = new RefineableEighthSphereMesh<ELEMENT>(radius);
211 
212  // Set error estimator
213  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
214  mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
215 
216  // Adjust error targets for adaptive refinement
217  if (CommandLineArgs::Argc>1)
218  {
219  // Validation: Relax tolerance to get nonuniform refinement during
220  // first step
221  mesh_pt()->max_permitted_error()=0.7;
222  mesh_pt()->min_permitted_error()=0.5;
223  }
224  else
225  {
226  mesh_pt()->max_permitted_error()=0.01;
227  mesh_pt()->min_permitted_error()=0.001;
228  } // end adjustment
229 
230  //Doc the mesh boundaries
231  ofstream some_file;
232  some_file.open("boundaries.dat");
233  mesh_pt()->output_boundaries(some_file);
234  some_file.close();
235 
236  // Set the boundary conditions for this problem: All nodes are
237  // free by default -- just pin the ones that have Dirichlet conditions
238  // here (all the nodes on the boundary)
239  unsigned num_bound = mesh_pt()->nboundary();
240  for(unsigned ibound=0;ibound<num_bound;ibound++)
241  {
242  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
243  for (unsigned inod=0;inod<num_nod;inod++)
244  {
245  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
246  }
247  } // end of pinning
248 
249 
250  //Find number of elements in mesh
251  unsigned n_element = mesh_pt()->nelement();
252 
253  // Loop over the elements to set up element-specific
254  // things that cannot be handled by constructor
255  for(unsigned i=0;i<n_element;i++)
256  {
257  // Upcast from FiniteElement to the present element
258  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
259 
260  //Set the source function pointer
261  el_pt->source_fct_pt() = Source_fct_pt;
262  }
263 
264  // Setup equation numbering
265  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
266 
267 } // end of constructor
268 
269 
270 
271 //========================start_of_doc====================================
272 /// Doc the solution
273 //========================================================================
274 template<class ELEMENT>
276 {
277  ofstream some_file;
278  char filename[100];
279 
280  // Number of plot points
281  unsigned npts;
282  npts=5;
283 
284 
285  // Output solution
286  //-----------------
287  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
288  doc_info.number());
289  some_file.open(filename);
290  mesh_pt()->output(some_file,npts);
291  some_file.close();
292 
293 
294  // Output exact solution
295  //----------------------
296  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
297  doc_info.number());
298  some_file.open(filename);
299  mesh_pt()->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
300  some_file.close();
301 
302 
303  // Doc error
304  //----------
305  double error,norm;
306  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
307  doc_info.number());
308  some_file.open(filename);
309  mesh_pt()->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
310  error,norm);
311  some_file.close();
312  cout << "error: " << sqrt(error) << std::endl;
313  cout << "norm : " << sqrt(norm) << std::endl << std::endl;
314 
315 } // end of doc
316 
317 
318 ////////////////////////////////////////////////////////////////////////
319 ////////////////////////////////////////////////////////////////////////
320 ////////////////////////////////////////////////////////////////////////
321 
322 
323 
324 //=========start_of_main=============================================
325 /// Driver for 3D Poisson problem in eighth of a sphere. Solution
326 /// has a sharp step. If there are
327 /// any command line arguments, we regard this as a validation run
328 /// and perform only a single adaptation.
329 //===================================================================
330 int main(int argc, char *argv[])
331 {
332 
333  // Store command line arguments
334  CommandLineArgs::setup(argc,argv);
335 
336  // Set up the problem with 27-node brick elements, pass pointer to
337  // source function
340 
341  // Setup labels for output
342  DocInfo doc_info;
343 
344  // 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  cout << "Self test: " << problem.self_test() << std::endl;
352 
353  // Solve the problem
354  problem.newton_solve();
355 
356  //Output solution
357  problem.doc_solution(doc_info);
358 
359  //Increment counter for solutions
360  doc_info.number()++;
361 
362  // Now do (up to) three rounds of fully automatic adapation in response to
363  // error estimate
364  unsigned max_solve;
365  if (CommandLineArgs::Argc>1)
366  {
367  // Validation run: Just one adaptation
368  max_solve=1;
369  cout << "Only doing one adaptation for validation" << std::endl;
370  }
371  else
372  {
373  // Up to four adaptations
374  max_solve=4;
375  }
376 
377  for (unsigned isolve=0;isolve<max_solve;isolve++)
378  {
379  // Adapt problem/mesh
380  problem.adapt();
381 
382  // Re-solve the problem if the adaptation has changed anything
383  if ((problem.mesh_pt()->nrefined() !=0)||
384  (problem.mesh_pt()->nunrefined()!=0))
385  {
386  problem.newton_solve();
387  }
388  else
389  {
390  cout << "Mesh wasn't adapted --> we'll stop here" << std::endl;
391  break;
392  }
393 
394  //Output solution
395  problem.doc_solution(doc_info);
396 
397  //Increment counter for solutions
398  doc_info.number()++;
399  }
400 
401 // pause("done");
402 
403 } // end of main
404 
405 
406 
407 
408 
409 
410 
411 
412 
double N_x
Orientation (non-normalised x-component of unit vector in direction of step plane) ...
double X_0
Orientation (x-coordinate of step plane)
Poisson problem in refineable eighth of a sphere mesh.
void actions_before_newton_solve()
Update the problem specs before solve: Set Dirchlet boundary conditions from exact solution...
PoissonEquations< 3 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
void get_exact_u(const Vector< double > &x, Vector< double > &u)
double Z_0
Orientation (z-coordinate of step plane)
~EighthSpherePoissonProblem()
Destructor: Empty.
double Alpha
Parameter for steepness of step.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
EighthSpherePoissonProblem(PoissonEquations< 3 >::PoissonSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function.
int main(int argc, char *argv[])
double Y_0
Orientation (y-coordinate of step plane)
double N_y
Orientation (non-normalised y-component of unit vector in direction of step plane) ...
void get_exact_u(const Vector< double > &x, double &u)
Exact solution as a scalar.
void doc_solution(DocInfo &doc_info)
Doc the solution.
double N_z
Orientation (non-normalised z-component of unit vector in direction of step plane) ...
RefineableEighthSphereMesh< ELEMENT > * mesh_pt()
Overload generic access function by one that returns a pointer to the specific mesh.
void get_source(const Vector< double > &x, double &source)
Source function to make it an exact solution.