fish_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-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 solution of 2D Poisson equation in fish-shaped domain with
27 // adaptivity
28 
29 
30 // Generic oomph-lib headers
31 #include "generic.h"
32 
33 // The Poisson equations
34 #include "poisson.h"
35 
36 // The fish mesh
37 #include "meshes/fish_mesh.h"
38 
39 using namespace std;
40 
41 using namespace oomph;
42 
43 
44 //============ start_of_namespace=====================================
45 /// Namespace for const source term in Poisson equation
46 //====================================================================
48 {
49 
50  /// Strength of source function: default value -1.0
51  double Strength=-1.0;
52 
53 /// Const source function
54  void source_function(const Vector<double>& x, double& source)
55  {
56  source = Strength;
57  }
58 
59 } // end of namespace
60 
61 
62 
63 
64 //======start_of_problem_class========================================
65 /// Refineable Poisson problem in fish-shaped domain.
66 /// Template parameter identifies the element type.
67 //====================================================================
68 template<class ELEMENT>
69 class RefineableFishPoissonProblem : public Problem
70 {
71 
72 public:
73 
74  /// Constructor
76 
77  /// Destructor: Empty
79 
80  /// Update the problem specs after solve (empty)
82 
83  /// Update the problem specs before solve (empty)
85 
86  /// Overloaded version of the problem's access function to
87  /// the mesh. Recasts the pointer to the base Mesh object to
88  /// the actual mesh type.
89  RefineableFishMesh<ELEMENT>* mesh_pt()
90  {
91  return dynamic_cast<RefineableFishMesh<ELEMENT>*>(Problem::mesh_pt());
92  }
93 
94  /// Doc the solution. Output directory and labels are specified
95  /// by DocInfo object
96  void doc_solution(DocInfo& doc_info);
97 
98 }; // end of problem class
99 
100 
101 
102 
103 
104 //===========start_of_constructor=========================================
105 /// Constructor for adaptive Poisson problem in fish-shaped
106 /// domain.
107 //========================================================================
108 template<class ELEMENT>
110 {
111 
112  // Build fish mesh -- this is a coarse base mesh consisting
113  // of four elements. We'll refine/adapt the mesh later.
114  Problem::mesh_pt()=new RefineableFishMesh<ELEMENT>;
115 
116  // Create/set error estimator
117  mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
118 
119  // Set the boundary conditions for this problem: All nodes are
120  // free by default -- just pin the ones that have Dirichlet conditions
121  // here. Since the boundary values are never changed, we set
122  // them here rather than in actions_before_newton_solve().
123  unsigned n_bound = mesh_pt()->nboundary();
124  for(unsigned i=0;i<n_bound;i++)
125  {
126  unsigned n_node = mesh_pt()->nboundary_node(i);
127  for (unsigned n=0;n<n_node;n++)
128  {
129  // Pin the single scalar value at this node
130  mesh_pt()->boundary_node_pt(i,n)->pin(0);
131 
132  // Assign the homogenous boundary condition for the one and only
133  // nodal value
134  mesh_pt()->boundary_node_pt(i,n)->set_value(0,0.0);
135  }
136  }
137 
138  // Loop over elements and set pointers to source function
139  unsigned n_element = mesh_pt()->nelement();
140  for(unsigned e=0;e<n_element;e++)
141  {
142  // Upcast from FiniteElement to the present element
143  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
144 
145  //Set the source function pointer
146  el_pt->source_fct_pt() = &ConstSourceForPoisson::source_function;
147  }
148 
149  // Setup the equation numbering scheme
150  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
151 
152 } // end of constructor
153 
154 
155 
156 
157 //=======start_of_doc=====================================================
158 /// Doc the solution in tecplot format.
159 //========================================================================
160 template<class ELEMENT>
162 {
163 
164  ofstream some_file;
165  char filename[100];
166 
167  // Number of plot points in each coordinate direction.
168  unsigned npts;
169  npts=5;
170 
171  // Output solution
172  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
173  doc_info.number());
174  some_file.open(filename);
175  mesh_pt()->output(some_file,npts);
176  some_file.close();
177 
178 } // end of doc
179 
180 
181 
182 
183 
184 
185 //=====================start_of_incremental===============================
186 /// Demonstrate how to solve 2D Poisson problem in
187 /// fish-shaped domain with mesh adaptation. First we solve on the original
188 /// coarse mesh. Next we do a few uniform refinement steps and re-solve.
189 /// Finally, we enter into an automatic adapation loop.
190 //========================================================================
192 {
193 
194  //Set up the problem with nine-node refineable Poisson elements
196 
197  // Setup labels for output
198  //------------------------
199  DocInfo doc_info;
200 
201  // Set output directory
202  doc_info.set_directory("RESLT_incremental");
203 
204  // Step number
205  doc_info.number()=0;
206 
207 
208  // Doc (default) refinement targets
209  //----------------------------------
210  problem.mesh_pt()->doc_adaptivity_targets(cout);
211 
212  // Solve/doc the problem on the initial, very coarse mesh
213  //-------------------------------------------------------
214 
215  // Solve the problem
216  problem.newton_solve();
217 
218  //Output solution
219  problem.doc_solution(doc_info);
220 
221  //Increment counter for solutions
222  doc_info.number()++;
223 
224 
225  // Do three rounds of uniform mesh refinement and re-solve
226  //--------------------------------------------------------
227  problem.refine_uniformly();
228  problem.refine_uniformly();
229  problem.refine_uniformly();
230 
231  // Solve the problem
232  problem.newton_solve();
233 
234  //Output solution
235  problem.doc_solution(doc_info);
236 
237  //Increment counter for solutions
238  doc_info.number()++;
239 
240 
241  // Now do (up to) four rounds of fully automatic adapation in response to
242  //-----------------------------------------------------------------------
243  // error estimate
244  //---------------
245  unsigned max_solve=4;
246  for (unsigned isolve=0;isolve<max_solve;isolve++)
247  {
248  // Adapt problem/mesh
249  problem.adapt();
250 
251  // Re-solve the problem if the adaptation has changed anything
252  if ((problem.mesh_pt()->nrefined() !=0)||
253  (problem.mesh_pt()->nunrefined()!=0))
254  {
255  problem.newton_solve();
256  }
257  else
258  {
259  cout << "Mesh wasn't adapted --> we'll stop here" << std::endl;
260  break;
261  }
262 
263  //Output solution
264  problem.doc_solution(doc_info);
265 
266  //Increment counter for solutions
267  doc_info.number()++;
268  }
269 
270 
271 } // end of incremental
272 
273 
274 
275 
276 
277 
278 //================================start_black_box=========================
279 /// Demonstrate how to solve 2D Poisson problem in
280 /// fish-shaped domain with fully automatic mesh adaptation
281 //========================================================================
283 {
284 
285  //Set up the problem with nine-node refineable Poisson elements
287 
288  // Setup labels for output
289  //------------------------
290  DocInfo doc_info;
291 
292  // Set output directory
293  doc_info.set_directory("RESLT_fully_automatic");
294 
295  // Step number
296  doc_info.number()=0;
297 
298 
299  // Doc (default) refinement targets
300  //----------------------------------
301  problem.mesh_pt()->doc_adaptivity_targets(cout);
302 
303 
304  // Solve/doc the problem with fully automatic adaptation
305  //------------------------------------------------------
306 
307  // Maximum number of adaptations:
308  unsigned max_adapt=5;
309 
310  // Solve the problem; perform up to specified number of adaptations.
311  problem.newton_solve(max_adapt);
312 
313  //Output solution
314  problem.doc_solution(doc_info);
315 
316 } // end black box
317 
318 
319 
320 
321 //=================start_of_main==========================================
322 /// Demonstrate how to solve 2D Poisson problem in
323 /// fish-shaped domain with mesh adaptation.
324 //========================================================================
325 int main()
326 {
327  // Solve with adaptation, docing the intermediate steps
329 
330  // Solve directly, with fully automatic adaptation
332 
333 } // end of main
334 
335 
336 
337 
338 // //Do few unrefinements and keep doc-ing
339 // for (unsigned i=0;i<5;i++)
340 // {
341 // problem.unrefine_uniformly();
342 
343 // //Output solution
344 // problem.doc_solution(doc_info);
345 
346 // //Increment counter for solutions
347 // doc_info.number()++;
348 // }
Refineable Poisson problem in fish-shaped domain. Template parameter identifies the element type.
Definition: fish_poisson.cc:70
virtual ~RefineableFishPoissonProblem()
Destructor: Empty.
Definition: fish_poisson.cc:78
void actions_before_newton_solve()
Update the problem specs before solve (empty)
Definition: fish_poisson.cc:84
RefineableFishPoissonProblem()
Constructor.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
Definition: fish_poisson.cc:81
RefineableFishMesh< ELEMENT > * mesh_pt()
Overloaded version of the problem's access function to the mesh. Recasts the pointer to the base Mesh...
Definition: fish_poisson.cc:89
void doc_solution(DocInfo &doc_info)
Doc the solution. Output directory and labels are specified by DocInfo object.
void solve_with_fully_automatic_adaptation()
Demonstrate how to solve 2D Poisson problem in fish-shaped domain with fully automatic mesh adaptatio...
void solve_with_incremental_adaptation()
Demonstrate how to solve 2D Poisson problem in fish-shaped domain with mesh adaptation....
int main()
Demonstrate how to solve 2D Poisson problem in fish-shaped domain with mesh adaptation.
Namespace for const source term in Poisson equation.
Definition: fish_poisson.cc:48
double Strength
Strength of source function: default value -1.0.
Definition: fish_poisson.cc:51
void source_function(const Vector< double > &x, double &source)
Const source function.
Definition: fish_poisson.cc:54