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