mesh_from_tetgen_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 function for a simple test poisson problem using a mesh
31 // generated from an input file generated by the tetrahedra mesh generator.
32 // the files .node, .ele and .face should be specified in this order
33 
34 //Generic routines
35 #include "generic.h"
36 
37 // Poisson equations
38 #include "poisson.h"
39 
40 // Get the mesh
41 #include "meshes/tetgen_mesh.h"
42 
43 using namespace std;
44 
45 using namespace oomph;
46 
47 //=============start_of_namespace=====================================
48 /// Namespace for exact solution for Poisson equation with sharp step
49 //====================================================================
50 namespace TanhSolnForPoisson
51 {
52 
53  /// Parameter for steepness of step
54  double Alpha=3;
55 
56  ///\short Orientation (non-normalised x-component of unit vector in direction
57  /// of step plane)
58  double N_x=-1.0;
59 
60  ///\short Orientation (non-normalised y-component of unit vector in direction
61  /// of step plane)
62  double N_y=-1.0;
63 
64  ///\short Orientation (non-normalised z-component of unit vector in direction
65  /// of step plane)
66  double N_z=1.0;
67 
68 
69  ///\short Orientation (x-coordinate of step plane)
70  double X_0=0.0;
71 
72  ///\short Orientation (y-coordinate of step plane)
73  double Y_0=0.0;
74 
75  ///\short Orientation (z-coordinate of step plane)
76  double Z_0=0.0;
77 
78 
79  // Exact solution as a Vector
80  void get_exact_u(const Vector<double>& x, Vector<double>& u)
81  {
82  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)*
83  N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*
84  N_z/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)));
85  }
86 
87  /// Exact solution as a scalar
88  void get_exact_u(const Vector<double>& x, double& u)
89  {
90  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)*
91  N_y/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)+(x[2]-Z_0)*
92  N_z/sqrt(N_x*N_x+N_y*N_y+N_z*N_z)));
93  }
94 
95 
96  /// Source function to make it an exact solution
97  void get_source(const Vector<double>& x, double& source)
98  {
99 
100  double s1,s2,s3,s4;
101 
102  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]-
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)))*(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]-
105  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*
106  N_z))),2.0))*Alpha*Alpha*N_x*N_x/(N_x*N_x+N_y*N_y+N_z*N_z);
107  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]-
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)))*(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]-
110  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*
111  N_z))),2.0))*Alpha*Alpha*N_y*N_y/(N_x*N_x+N_y*N_y+N_z*N_z);
112  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]-
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)))*(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]-
115  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*
116  N_z))),2.0))*Alpha*Alpha*N_z*N_z/(N_x*N_x+N_y*N_y+N_z*N_z);
117  s2 = s3+s4;
118  source = s1+s2;
119  }
120 
121 
122 } // end of namespace
123 
124 
125 
126 
127 
128 //====================================================================
129 /// Micky mouse Poisson problem.
130 //====================================================================
131 template<class ELEMENT>
132 class PoissonProblem : public Problem
133 {
134 
135 public:
136 
137 
138  /// Constructor
139  PoissonProblem(PoissonEquations<3>::PoissonSourceFctPt source_fct_pt,
140  const string& node_file_name,
141  const string& element_file_name,
142  const string& face_file_name);
143 
144  /// Destructor (empty)
146 
147 
148  /// Update the problem specs before solve: (Re)set boundary conditions
150  {
151  //Loop over the boundaries
152  unsigned num_bound = mesh_pt()->nboundary();
153  for(unsigned ibound=0;ibound<num_bound;ibound++)
154  {
155  // Loop over the nodes on boundary
156  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
157  for (unsigned inod=0;inod<num_nod;inod++)
158  {
159  Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
160  double u;
161  Vector<double> x(3);
162  x[0]=nod_pt->x(0);
163  x[1]=nod_pt->x(1);
164  x[2]=nod_pt->x(2);
166  nod_pt->set_value(0,u);
167  }
168  }
169  }
170 
171  /// Update the problem specs before solve (empty)
173 
174 
175  //Access function for the specific mesh
177  {
178  return dynamic_cast<TetgenMesh<ELEMENT>*>(Problem::mesh_pt());
179  }
180 
181  /// Doc the solution
182  void doc_solution(const unsigned& nplot, DocInfo& doc_info);
183 
184 private:
185 
186  /// Pointer to source function
187  PoissonEquations<3>::PoissonSourceFctPt Source_fct_pt;
188 
189 };
190 
191 
192 
193 //========================================================================
194 /// Constructor for Poisson problem
195 //========================================================================
196 template<class ELEMENT>
198  PoissonProblem(PoissonEquations<3>::PoissonSourceFctPt source_fct_pt,
199  const string& node_file_name,
200  const string& element_file_name,
201  const string& face_file_name)
202  : Source_fct_pt(source_fct_pt)
203 {
204 
205  // Setup parameters for exact tanh solution
206 
207  // Steepness of step
209 
210 
211  //Create mesh
212  Problem::mesh_pt() = new TetgenMesh<ELEMENT>(node_file_name,
213  element_file_name,
214  face_file_name);
215 
216  // Set the boundary conditions for this problem: All nodes are
217  // free by default -- just pin the ones that have Dirichlet conditions
218  // here.
219  unsigned num_bound = mesh_pt()->nboundary();
220  for(unsigned ibound=0;ibound<num_bound;ibound++)
221  {
222  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
223  for (unsigned inod=0;inod<num_nod;inod++)
224  {
225  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
226  }
227  }
228 
229  // Complete the build of all elements so they are fully functional
230 
231  //Find number of elements in mesh
232  unsigned n_element = mesh_pt()->nelement();
233 
234  // Loop over the elements to set up element-specific
235  // things that cannot be handled by constructor
236  for(unsigned i=0;i<n_element;i++)
237  {
238  // Upcast from GeneralElement to the present element
239  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
240 
241  //Set the source function pointer
242  el_pt->source_fct_pt() = Source_fct_pt;
243  }
244 
245  // Setup equation numbering scheme
246  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
247 
248 }
249 
250 
251 
252 //========================================================================
253 /// Doc the solution
254 //========================================================================
255 template<class ELEMENT>
256 void PoissonProblem<ELEMENT>::doc_solution(const unsigned& nplot,
257  DocInfo& doc_info)
258 {
259 
260  ofstream some_file;
261  char filename[100];
262 
263 
264  // Doc local node numbering
265  //-------------------------
266  sprintf(filename,"%s/node_numbering%i.dat",doc_info.directory().c_str(),
267  doc_info.number());
268  some_file.open(filename);
269  FiniteElement* el_pt=mesh_pt()->finite_element_pt(0);
270  unsigned nnode=el_pt->nnode();
271  unsigned ndim=el_pt->node_pt(0)->ndim();
272  for (unsigned j=0;j<nnode;j++)
273  {
274  for (unsigned i=0;i<ndim;i++)
275  {
276  some_file << el_pt->node_pt(j)->x(i) << " " ;
277  }
278  some_file << j << std::endl;
279  }
280  some_file.close();
281 
282  // Output boundaries
283  //------------------
284  sprintf(filename,"%s/boundaries%i.dat",doc_info.directory().c_str(),
285  doc_info.number());
286  some_file.open(filename);
287  mesh_pt()->output_boundaries(some_file);
288  some_file.close();
289 
290 
291  // Output solution
292  //----------------
293  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
294  doc_info.number());
295  some_file.open(filename);
296  mesh_pt()->output(some_file,nplot);
297  some_file.close();
298 
299 
300  // Output exact solution
301  //----------------------
302  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
303  doc_info.number());
304  some_file.open(filename);
305  mesh_pt()->output_fct(some_file,nplot,TanhSolnForPoisson::get_exact_u);
306  some_file.close();
307 
308  // Doc error
309  //----------
310  double error,norm;
311  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
312  doc_info.number());
313  some_file.open(filename);
314  mesh_pt()->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
315  error,norm);
316  some_file.close();
317  cout << "error: " << sqrt(error) << std::endl;
318  cout << "norm : " << sqrt(norm) << std::endl << std::endl;
319 
320 } // end of doc
321 
322 
323 
324 
325 //========================================================================
326 /// Demonstrate how to solve Poisson problem
327 //========================================================================
328 int main(int argc, char* argv[])
329 {
330 
331 // Store command line arguments
332  CommandLineArgs::setup(argc,argv);
333 
334  // Check number of command line arguments: Need exactly three.
335  if (argc!=4)
336  {
337  std::string error_message =
338  "Wrong number of command line arguments.\n";
339  error_message +=
340  "Must specify the following file names \n";
341  error_message +=
342  "filename.node then filename.ele then filename.face\n";
343 
344  throw OomphLibError(error_message,
345  OOMPH_CURRENT_FUNCTION,
346  OOMPH_EXCEPTION_LOCATION);
347  }
348 
349  // Convert arguments to strings that specify the input file names
350  string node_file_name(argv[1]);
351  string element_file_name(argv[2]);
352  string face_file_name(argv[3]);
353 
354 
355  // Label for output
356  DocInfo doc_info;
357 
358  // Output directory
359  doc_info.set_directory("RESLT");
360 
361  // Number of output points per edge
362  unsigned nplot=2;
363 
364  // Do the problem with linear elements
365  //------------------------------------
366  {
369  node_file_name,element_file_name,face_file_name);
370 
371  // Solve the problem
372  problem.newton_solve();
373 
374  //Output solution with 2 points per edge
375  nplot=2;
376  problem.doc_solution(nplot,doc_info);
377 
378  //Increment counter for solutions
379  doc_info.number()++;
380 
381  //Output solution with 5 points per edge
382  nplot=5;
383  problem.doc_solution(nplot,doc_info);
384 
385  //Increment counter for solutions
386  doc_info.number()++;
387 
388 
389  }
390 
391 
392 
393  // Do the problem with quadratic elements
394  //---------------------------------------
395  {
398  node_file_name,element_file_name,face_file_name);
399 
400  // Solve the problem
401  problem.newton_solve();
402 
403  //Output solution with 2 points per edge
404  nplot=2;
405  problem.doc_solution(nplot,doc_info);
406 
407  //Increment counter for solutions
408  doc_info.number()++;
409 
410  //Output solution with 5 points per edge
411  nplot=5;
412  problem.doc_solution(nplot,doc_info);
413 
414  //Increment counter for solutions
415  doc_info.number()++;
416 
417  }
418 
419 
420 }
421 
422 
423 
double N_x
Orientation (non-normalised x-component of unit vector in direction of step plane) ...
TetgenMesh< ELEMENT > * mesh_pt()
double X_0
Orientation (x-coordinate of step plane)
~PoissonProblem()
Destructor (empty)
void actions_after_newton_solve()
Update the problem specs before solve (empty)
void get_exact_u(const Vector< double > &x, Vector< double > &u)
double Z_0
Orientation (z-coordinate of step plane)
double Alpha
Parameter for steepness of step.
Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen//.
void actions_before_newton_solve()
Update the problem specs before solve: (Re)set boundary conditions.
PoissonEquations< 3 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
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.
double N_z
Orientation (non-normalised z-component of unit vector in direction of step plane) ...
void doc_solution(const unsigned &nplot, DocInfo &doc_info)
Doc the solution.
Micky mouse Poisson problem.
PoissonProblem(PoissonEquations< 3 >::PoissonSourceFctPt source_fct_pt, const string &node_file_name, const string &element_file_name, const string &face_file_name)
Constructor.
void get_source(const Vector< double > &x, double &source)
Source function to make it an exact solution.
int main(int argc, char *argv[])
Demonstrate how to solve Poisson problem.