mesh_from_geompack_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-2022 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 function for a simple test poisson problem using a mesh
27// generated from an input file generated by the quadrilateral mesh generator
28// Geompack.
29// files .mh2 and .cs2 should be specified in this order
30
31//Generic routines
32#include "generic.h"
33
34// Poisson equations
35#include "poisson.h"
36
37// The mesh
38#include "meshes/geompack_mesh.h"
39
40using namespace std;
41
42
43using namespace oomph;
44
45//====================================================================
46/// Namespace for exact solution for Poisson equation with sharp step
47//====================================================================
49{
50
51 /// Parameter for steepness of step
52 double Alpha;
53
54 /// Parameter for angle of step
55 double Beta;
56
57
58 /// Exact solution as a Vector
59 void get_exact_u(const Vector<double>& x, Vector<double>& u)
60 {
61 u[0]=tanh(1.0-Alpha*(Beta*x[0]-x[1]));
62 }
63
64
65 /// Exact solution as a scalar
66 void get_exact_u(const Vector<double>& x, double& u)
67 {
68 u=tanh(1.0-Alpha*(Beta*x[0]-x[1]));
69 }
70
71
72 /// Source function to make it an exact solution
73 void get_source(const Vector<double>& x, double& source)
74 {
75 source = 2.0*tanh(-1.0+Alpha*(Beta*x[0]-x[1]))*
76 (1.0-pow(tanh(-1.0+Alpha*(Beta*x[0]-x[1])),2.0))*
77 Alpha*Alpha*Beta*Beta+2.0*tanh(-1.0+Alpha*(Beta*x[0]-x[1]))*
78 (1.0-pow(tanh(-1.0+Alpha*(Beta*x[0]-x[1])),2.0))*Alpha*Alpha;
79 }
80
81}
82
83
84
85
86
87
88
89//====================================================================
90/// Micky mouse Poisson problem.
91//====================================================================
92template<class ELEMENT>
93class PoissonProblem : public Problem
94{
95
96public:
97
98
99 /// Constructor
100 PoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt,
101 const string& mesh_file_name,
102 const string& curve_file_name);
103
104 /// Destructor (empty)
106
107 /// Update the problem specs before solve: (Re)set boundary conditions
109 {
110 //Loop over the boundaries
111 unsigned num_bound = mesh_pt()->nboundary();
112 for(unsigned ibound=0;ibound<num_bound;ibound++)
113 {
114 // Loop over the nodes on boundary
115 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
116 for (unsigned inod=0;inod<num_nod;inod++)
117 {
118 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
119 double u;
120 Vector<double> x(2);
121 x[0]=nod_pt->x(0);
122 x[1]=nod_pt->x(1);
124 nod_pt->set_value(0,u);
125 }
126 }
127 }
128
129 /// Update the problem specs before solve (empty)
131 {}
132
133
134 //Access function for the specific mesh
136 {
137 return dynamic_cast<GeompackQuadMesh<ELEMENT>*>(Problem::mesh_pt());
138 }
139
140 /// Doc the solution
141 void doc_solution(DocInfo& doc_info);
142
143private:
144
145 /// Pointer to source function
146 PoissonEquations<2>::PoissonSourceFctPt Source_fct_pt;
147
148};
149
150
151
152//========================================================================
153/// Constructor for Poisson problem
154//========================================================================
155template<class ELEMENT>
157 PoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt,
158 const string& mesh_file_name,
159 const string& curve_file_name)
160 : Source_fct_pt(source_fct_pt)
161{
162
163 // Setup parameters for exact tanh solution
164
165 // Steepness of step
167
168 // Orientation of step
170
171 //Create mesh
172 Problem::mesh_pt() = new GeompackQuadMesh<ELEMENT>(mesh_file_name,
173 curve_file_name);
174
175 // Set the boundary conditions for this problem: All nodes are
176 // free by default -- just pin the ones that have Dirichlet conditions
177 // here.
178 unsigned num_bound = mesh_pt()->nboundary();
179 for(unsigned ibound=0;ibound<num_bound;ibound++)
180 {
181 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
182 for (unsigned inod=0;inod<num_nod;inod++)
183 {
184 mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
185 }
186 }
187
188 // Complete the build of all elements so they are fully functional
189
190 //Find number of elements in mesh
191 unsigned n_element = mesh_pt()->nelement();
192
193 // Loop over the elements to set up element-specific
194 // things that cannot be handled by constructor
195 for(unsigned i=0;i<n_element;i++)
196 {
197 // Upcast from GeneralElement to the present element
198 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
199
200 //Set the source function pointer
201 el_pt->source_fct_pt() = Source_fct_pt;
202 }
203
204 // Setup equation numbering scheme
205 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
206
207}
208
209
210
211//========================================================================
212/// Doc the solution
213//========================================================================
214template<class ELEMENT>
216{
217
218 ofstream some_file;
219 char filename[100];
220
221 // Number of plot points
222 unsigned npts;
223 npts=5;
224
225
226
227 // Output boundaries
228 //------------------
229 sprintf(filename,"%s/boundaries.dat",doc_info.directory().c_str());
230 some_file.open(filename);
231 mesh_pt()->output_boundaries(some_file);
232 some_file.close();
233
234
235 // Output solution
236 //----------------
237 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
238 doc_info.number());
239 some_file.open(filename);
240 mesh_pt()->output(some_file,npts);
241 some_file.close();
242
243
244 // Output exact solution
245 //----------------------
246 sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
247 doc_info.number());
248 some_file.open(filename);
249 mesh_pt()->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
250 some_file.close();
251
252
253// Doc error
254 //----------
255 double error,norm;
256 sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
257 doc_info.number());
258 some_file.open(filename);
259 mesh_pt()->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
260 error,norm);
261 some_file.close();
262 cout << "error: " << sqrt(error) << std::endl;
263 cout << "norm : " << sqrt(norm) << std::endl << std::endl;
264
265
266 }
267
268
269
270
271
272
273
274/// /////////////////////////////////////////////////////////////////////
275/// /////////////////////////////////////////////////////////////////////
276/// /////////////////////////////////////////////////////////////////////
277
278
279
280//========================================================================
281/// Demonstrate how to solve Poisson problem
282//========================================================================
283int main(int argc, char* argv[])
284{
285// Store command line arguments
286 CommandLineArgs::setup(argc,argv);
287
288 // Check number of command line arguments: Need exactly two.
289 if (argc!=3)
290 {
291 std::string error_message =
292 "Wrong number of command line arguments.\n";
293 error_message +=
294 "Must specify the following file names \n";
295 error_message +=
296 "filename.mh2 then filename.cs2\n";
297
298 throw OomphLibError(error_message,
299 OOMPH_CURRENT_FUNCTION,
300 OOMPH_EXCEPTION_LOCATION);
301 }
302
303 // Convert arguments to strings that specify the input file names
304 string mesh_file_name(argv[1]);
305 string curve_file_name(argv[2]);
306
307 //Set up the problem
310 mesh_file_name,curve_file_name);
311
312 /// Label for output
313 DocInfo doc_info;
314
315 // Output directory
316 doc_info.set_directory("RESLT");
317
318 // Solve the problem
319 problem.newton_solve();
320
321 //Output solution
322 problem.doc_solution(doc_info);
323
324}
325
326
327
Micky mouse Poisson problem.
PoissonEquations< 2 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
void actions_before_newton_solve()
Update the problem specs before solve: (Re)set boundary conditions.
GeompackQuadMesh< ELEMENT > * mesh_pt()
PoissonProblem(PoissonEquations< 2 >::PoissonSourceFctPt source_fct_pt, const string &mesh_file_name, const string &curve_file_name)
Constructor.
void actions_after_newton_solve()
Update the problem specs before solve (empty)
void doc_solution(DocInfo &doc_info)
Doc the solution.
~PoissonProblem()
Destructor (empty)
Quadrilateral mesh generator; Uses input from Geompack++. See: http://members.shaw....
int main(int argc, char *argv[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Namespace for exact solution for Poisson equation with sharp step.
double Beta
Parameter for angle of step.
void get_source(const Vector< double > &x, double &source)
Source function to make it 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.
////////////////////////////////////////////////////////////////////// //////////////////////////////...