simple_spine_channel.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 for 2-D Channel with changing width, using Taylor Hood
27// and Crouzeix Raviart elements.
28
29// Generic oomph-lib header
30#include "generic.h"
31
32// Navier Stokes headers
33#include "navier_stokes.h"
34
35// The mesh
36#include "meshes/simple_rectangular_quadmesh.h"
37
38using namespace std;
39
40using namespace oomph;
41
42/// ////////////////////////////////////////////////////////////////////
43/// ////////////////////////////////////////////////////////////////////
44
45
46
47//=========================================================================
48/// Geometric object representing a wavy wall, parametrised by
49/// \f[ x = \zeta \f]
50/// \f[ y = A (1 - \cos(2\pi\zeta/L) )\f]
51//=========================================================================
52class WavyWall : public GeomObject
53{
54
55public:
56
57 /// Constructor: Pass wavelength and amplitude
58 WavyWall(const double& l, const double& amplitude)
59 : GeomObject(1,2), L(l), A(amplitude) {}
60
61 /// Position vector to wavy wall.
62 void position(const Vector<double>& zeta, Vector<double>& r) const
63 {
64 r[0]=zeta[0];
65 r[1]=A*(1.0-cos(2.0*MathematicalConstants::Pi*zeta[0]/L));
66 }
67
68protected:
69
70 /// Wavelength
71 double L;
72
73 /// Amplitude
74 double A;
75
76
77};
78
79
80/// ////////////////////////////////////////////////////////////////////
81/// ////////////////////////////////////////////////////////////////////
82
83
84//======================================================================
85/// Simple spine mesh class derived from standard 2D mesh.
86/// Vertical lines of nodes are located on spines.
87//======================================================================
88template <class ELEMENT>
89class SimpleSpineMesh : public SimpleRectangularQuadMesh<ELEMENT >,
90 public SpineMesh
91{
92
93public:
94
95 /// Constructor: Pass number of elements in x-direction,
96 /// number of elements in y-direction, length in x direction, initial
97 /// height of mesh, and pointer to timestepper (defaults to
98 /// Steady timestepper)
99 SimpleSpineMesh(const unsigned &nx,
100 const unsigned &ny,
101 const double &lx,
102 const double &h,
103 GeomObject* substrate_pt,
104 TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper);
105
106
107 /// Function pointer to function, h(x), that may be used
108 /// to specify the "height" of the domain, by assigning the function
109 /// values to the spine heights.
110 typedef double (*HeightFctPt)(const double& x);
111
112 /// Access function: Pointer to height function
114 {
115 return Height_fct_pt;
116 }
117
118
119 /// Height function -- this is called by update_spine_heights()
120 /// when spine heights are assigned.
121 double height_fct(const double& x)
122 {
123 // Resolve function pointer if non-NULL
124 if (Height_fct_pt==0)
125 {
126 return Default_height;
127 }
128 else
129 {
130 return Height_fct_pt(x);
131 }
132 }
133
134
135 /// Update the spine heights according to the function specified
136 /// by height_fct_pt().
138 {
139 unsigned n_spine=nspine();
140 for (unsigned i=0;i<n_spine;i++)
141 {
142 // Zero-th geometric parameter of the spine is the x-coordinate
143 // of its basis
144 double x=spine_pt(i)->geom_parameter(0);
145
146 // Set spine height
147 spine_pt(i)->height()=height_fct(x);
148 }
149 }
150
151
152 /// General node update function implements pure virtual function
153 /// defined in SpineMesh base class and performs specific node update
154 /// actions: Nodes are located along vertical "spines" that emanate
155 /// from the "substrate" (the lower wall) specified by a
156 /// GeomObject.
157 virtual void spine_node_update(SpineNode* spine_node_pt)
158 {
159
160 //Get fraction along the spine
161 double w = spine_node_pt->fraction();
162
163 // Get height of spine
164 double h = spine_node_pt->spine_pt()->height();
165
166 // Get x-position of spine origin (as vector)
167 Vector<double> x_spine_origin(1);
168 x_spine_origin[0]=spine_node_pt->spine_pt()->geom_parameter(0);
169
170 // Get position vector to substrate (lower wall) for this spine
171 Vector<double> spine_base(2);
172 spine_node_pt->spine_pt()->geom_object_pt(0)->
173 position(x_spine_origin,spine_base);
174
175 //Update the nodal position
176 spine_node_pt->x(0) = spine_base[0];
177 spine_node_pt->x(1) = spine_base[1]+w*h;
178 }
179
180
181private:
182
183 /// Default height
185
186 /// Pointer to height function
188
189 /// Pointer to GeomObject that specifies the "substrate" (the lower wall)
190 GeomObject* Substrate_pt;
191
192};
193
194
195
196//===========================================================================
197/// Constructor for spine 2D mesh: Pass number of elements in x-direction,
198/// number of elements in y-direction, axial length and height of layer,
199/// pointer to geometric object that specifies the substrate (the lower wall)
200/// and pointer to timestepper (defaults to Static timestepper).
201///
202/// The mesh contains a layer of spine-ified fluid elements (of type ELEMENT;
203/// e.g SpineElement<QCrouzeixRaviartElement<2>)
204//===========================================================================
205template<class ELEMENT>
207 const unsigned &ny,
208 const double &lx,
209 const double &h,
210 GeomObject* substrate_pt,
211 TimeStepper* time_stepper_pt) :
212 SimpleRectangularQuadMesh<ELEMENT>(nx,ny,lx,h,time_stepper_pt),
213 Default_height(h), Height_fct_pt(0), Substrate_pt(substrate_pt)
214
215{
216
217 // We've already called the constructor for the SimpleRectangularQuadMesh
218 // which sets up nodal positons/topology etc. Now attach the spine
219 // information
220
221
222 // Allocate memory for the spines and fractions along spines
223
224 // Read out number of linear points in the element
225 unsigned np = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
226
227 // Number of spines
228 unsigned nspine = (np-1)*nx+1;
229 Spine_pt.reserve(nspine);
230
231 // Set up x-increments between spine origins
232 double x_lo=0.0;
233 double dx=lx/double(nx);
234
235 //FIRST SPINE
236 //-----------
237
238 //Element 0
239 //Node 0
240 //Assign the new spine with specified height
241 Spine* new_spine_pt=new Spine(h);
242 Spine_pt.push_back(new_spine_pt);
243
244 // Get pointer to node
245 SpineNode* nod_pt=element_node_pt(0,0);
246
247 //Pass the pointer to the spine to the node
248 nod_pt->spine_pt() = new_spine_pt;
249
250 //Set the node's fraction along the spine
251 nod_pt->fraction() = 0.0;
252
253 // Set the pointer to the mesh that implements the update fct for this node
254 nod_pt->spine_mesh_pt() = this;
255
256 // Set update fct id (not really needed here...)
257 nod_pt->node_update_fct_id()=0;
258
259
260 // Create vector that stores additional geometric parameters that
261 // are required during the execution of the remesh function:
262 // For our remesh function, we need the x-coordinate of the
263 // spine origin.
264
265 // Vector with enough storage for one geometic parameter
266 Vector<double> parameters(1);
267
268 // Store the x-coordinate in it
269 parameters[0]=x_lo;
270
271 // Pass the parameter to the spine
272 new_spine_pt->set_geom_parameter(parameters);
273
274
275
276 // The remesh function involving this spine only involves
277 // a single geometric object: The substrate
278 Vector<GeomObject*> geom_object_pt(1);
279 geom_object_pt[0] = Substrate_pt;
280
281 // Pass geom object(s) to spine
282 new_spine_pt->set_geom_object_pt(geom_object_pt);
283
284
285 //Loop vertically along the first spine
286
287 //Loop over the elements
288 for(unsigned i=0;i<ny;i++)
289 {
290 //Loop over the vertical nodes, apart from the first
291 for(unsigned l1=1;l1<np;l1++)
292 {
293 // Get pointer to node
294 SpineNode* nod_pt=element_node_pt(i*nx,l1*np);
295
296 //Pass the pointer to the spine to the node
297 nod_pt->spine_pt() = new_spine_pt;
298
299 //Set the fraction
300 nod_pt->fraction()=(double(i)+double(l1)/double(np-1))/double(ny);
301
302 // Pointer to the mesh that implements the update fct
303 nod_pt->spine_mesh_pt() = this;
304
305 // Set update fct id (not really needed here)
306 nod_pt->node_update_fct_id()=0;
307
308 }
309 } // end loop over elements
310
311
312 //LOOP OVER OTHER SPINES
313 //----------------------
314
315 //Now loop over the elements horizontally
316 for(unsigned j=0;j<nx;j++)
317 {
318 //Loop over the nodes in the elements horizontally, ignoring
319 //the first column
320 unsigned npmax=np;
321 for(unsigned l2=1;l2<npmax;l2++)
322 {
323 //Assign the new spine with specified height
324 new_spine_pt=new Spine(h);
325
326 // Add to collection of spines
327 Spine_pt.push_back(new_spine_pt);
328
329 // Get the node
330 SpineNode* nod_pt=element_node_pt(j,l2);
331
332 //Pass the pointer to the spine to the node
333 nod_pt->spine_pt() = new_spine_pt;
334
335 //Set the fraction
336 nod_pt->fraction() = 0.0;
337
338 // Pointer to the mesh that implements the update fct
339 nod_pt->spine_mesh_pt() = this;
340
341 // Set update fct id (not really needed here)
342 nod_pt->node_update_fct_id()=0;
343
344 // Create vector that stores additional geometric parameters that
345 // are required during the execution of the remesh function:
346 // For our remesh function, we need the x-coordinate of the
347 // spine origin.
348
349 // Vector with enough storage for one geometic parameter
350 Vector<double> parameters(1);
351
352 // Store the x-coordinate in it
353 parameters[0]=x_lo+double(j)*dx + double(l2)*dx/double(np-1);
354
355 // Pass the parameter to the spine
356 new_spine_pt->set_geom_parameter(parameters);
357
358
359 // The remesh function involving this spine only involves
360 // a single geometric object: The substrate
361 Vector<GeomObject*> geom_object_pt(1);
362 geom_object_pt[0] = Substrate_pt;
363
364 // Pass geom object(s) to spine
365 new_spine_pt->set_geom_object_pt(geom_object_pt);
366
367 //Loop vertically along the spine
368 //Loop over the elements
369 for(unsigned i=0;i<ny;i++)
370 {
371 //Loop over the vertical nodes, apart from the first
372 for(unsigned l1=1;l1<np;l1++)
373 {
374 // Get the node
375 SpineNode* nod_pt=element_node_pt(i*nx+j,l1*np+l2);
376
377 //Set the pointer to the spine
378 nod_pt->spine_pt() = new_spine_pt;
379
380 //Set the fraction
381 nod_pt->fraction()=(double(i)+double(l1)/double(np-1))/double(ny);
382
383 // Pointer to the mesh that implements the update fct
384 nod_pt->spine_mesh_pt() = this;
385
386 // Set update fct id (not really needed here)
387 nod_pt->node_update_fct_id()=0;
388 }
389 }
390 }
391 }
392
393} // end of constructor
394
395
396
397
398
399/// ////////////////////////////////////////////////////////////////////
400/// ////////////////////////////////////////////////////////////////////
401
402
403
404//==start_of_namespace===================================================
405/// Namespace for physical parameters
406//=======================================================================
408{
409
410 /// Reynolds number
411 double Re=100;
412
413 /// Start of indented region
414 double X_indent_start=0.5;
415
416 /// Length of indented region
417 double L=1.0;
418
419 /// Total length of domain
420 double L_total=4.0;
421
422 /// Undeformed height of domain
423 double H=1.0;
424
425 /// Amplitude of indentation
426 double A=-0.6;
427
428 /// Height of domain
429 double height(const double& x)
430 {
431 if ((x>X_indent_start)&&(x<(X_indent_start+L)))
432 {
433 return H + A*sin(MathematicalConstants::Pi*(x-X_indent_start)/L);
434 }
435 else
436 {
437 return H;
438 }
439 }
440
441} // end_of_namespace
442
443
444
445
446
447/// ////////////////////////////////////////////////////////////////////
448/// ////////////////////////////////////////////////////////////////////
449/// ////////////////////////////////////////////////////////////////////
450
451
452
453//==start_of_problem_class============================================
454/// Channel flow through a non-uniform channel whose geometry is defined
455/// by a spine mesh.
456//====================================================================
457template<class ELEMENT>
458class ChannelSpineFlowProblem : public Problem
459{
460
461public:
462
463 /// Constructor
465
466 /// Destructor: (empty)
468
469 /// Update the problem specs before solve. Update the nodal
470 /// positions
472 {
473 // Update the mesh
474 mesh_pt()->node_update();
475
476 } // end_of_actions_before_newton_solve
477
478 /// Update the after solve (empty)
480
481 /// Doc the solution
482 void doc_solution(DocInfo& doc_info);
483
484 /// Overload access to mesh
486 {
487 return dynamic_cast<SimpleSpineMesh<ELEMENT>*>(Problem::mesh_pt());
488 }
489
490
491private:
492
493 /// Width of channel
494 double Ly;
495
496
497}; // end_of_problem_class
498
499
500
501//==start_of_constructor==================================================
502/// Constructor for ChannelSpineFlow problem
503//========================================================================
504template<class ELEMENT>
506{
507
508 // Setup mesh
509
510 // Domain length in x-direction
512
513 // Domain length in y-direction
515
516 // # of elements in x-direction
517 unsigned Nx=40;
518
519 // # of elements in y-direction
520 unsigned Ny=10;
521
522 // Substrate (lower wall): A wavy wall
523 GeomObject* substrate_pt=
525
526 // Build and assign mesh -- pass pointer to geometric object
527 // that represents the sinusoidal bump on the upper wall
528 Problem::mesh_pt() = new SimpleSpineMesh<ELEMENT>(Nx,Ny,Lx,Ly,substrate_pt);
529
530 // Set function pointer for height function
531 mesh_pt()->height_fct_pt()=&Global_Physical_Variables::height;
532
533 // Update spine heights according to the function just specified
534 mesh_pt()->update_spine_heights();
535
536 // Update nodal positions
537 mesh_pt()->node_update();
538
539 // Pin all spine heights
540 unsigned nspine=mesh_pt()->nspine();
541 for (unsigned i=0;i<nspine;i++)
542 {
543 mesh_pt()->spine_pt(i)->spine_height_pt()->pin(0);
544 }
545
546
547 // Set the boundary conditions for this problem: All nodes are
548 // free by default -- just pin the ones that have Dirichlet conditions
549 // here: All boundaries are Dirichlet boundaries, except on boundary 1
550 unsigned num_bound = mesh_pt()->nboundary();
551 for(unsigned ibound=0;ibound<num_bound;ibound++)
552 {
553 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
554 for (unsigned inod=0;inod<num_nod;inod++)
555 {
556 if (ibound!=1)
557 {
558 // Loop over values (u and v velocities)
559 for (unsigned i=0;i<2;i++)
560 {
561 mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
562 }
563 }
564 else
565 {
566 // Parallel outflow ==> no-slip
567 mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
568 }
569 }
570 } // end loop over boundaries
571
572
573 // No slip on stationary upper and lower walls (boundaries 0 and 2)
574 // and parallel outflow (boundary 1)
575 for (unsigned ibound=0;ibound<num_bound-1;ibound++)
576 {
577 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
578 for (unsigned inod=0;inod<num_nod;inod++)
579 {
580 if (ibound!=1)
581 {
582 for (unsigned i=0;i<2;i++)
583 {
584 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
585 }
586 }
587 else
588 {
589 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
590 }
591 }
592 }
593
594
595 // Setup parabolic inflow along boundary 3:
596 unsigned ibound=3;
597 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
598 for (unsigned inod=0;inod<num_nod;inod++)
599 {
600 double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
601 // Parallel, parabolic inflow
602 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,y*(Ly-y));
603 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
604 }
605
606
607 // Loop over the elements to set up element-specific
608 // things that cannot be handled by constructor: Pass
609 // pointer to Reynolds number
610 unsigned n_element = mesh_pt()->nelement();
611 for(unsigned e=0;e<n_element;e++)
612 {
613 // Upcast from GeneralisedElement to the present element
614 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
615 //Set the Reynolds number
616 el_pt->re_pt() = &Global_Physical_Variables::Re;
617 }
618
619 // Setup equation numbering scheme
620 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
621
622} // end_of_constructor
623
624
625
626
627//==start_of_doc_solution=================================================
628/// Doc the solution
629//========================================================================
630template<class ELEMENT>
632{
633
634 ofstream some_file;
635 char filename[100];
636
637 // Number of plot points
638 unsigned npts=5;
639
640 // Output solution
641 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
642 doc_info.number());
643 some_file.open(filename);
644 mesh_pt()->output(some_file,npts);
645 some_file.close();
646
647} // end_of_doc_solution
648
649
650
651//==start_of_main======================================================
652/// Driver for channel flow problem with spine mesh.
653//=====================================================================
654int main()
655{
656
657 // Set output directory
658 DocInfo doc_info;
659 doc_info.set_directory("RESLT");
660 doc_info.number()=0;
661
662 // Solve problem with Taylor Hood elements
663 //---------------------------------------
664 {
665 //Build problem
667
668 // Solve the problem
669 problem.newton_solve();
670
671 //Output solution
672 problem.doc_solution(doc_info);
673
674 // Step number
675 doc_info.number()++;
676
677 } // end of Taylor Hood elements
678
679
680 // Solve problem with Crouzeix Raviart elements
681 //--------------------------------------------
682 {
683 // Build problem
685 problem;
686
687 // Solve the problem with automatic adaptation
688 problem.newton_solve();
689
690 //Output solution
691 problem.doc_solution(doc_info);
692 // Step number
693 doc_info.number()++;
694
695 } // end of Crouzeix Raviart elements
696
697
698
699} // end_of_main
700
701
702
703
704// //====================================================================
705// /// Create the files to illustrate the sparse node update operation
706// //====================================================================
707// void doc_sparse_node_update()
708// {
709
710// // Set output directory
711// DocInfo doc_info;
712// doc_info.set_directory("RESLT");
713
714// // Setup mesh
715
716// // # of elements in x-direction
717// unsigned Nx=5;
718
719// // # of elements in y-direction
720// unsigned Ny=5;
721
722// // Domain length in x-direction
723// double Lx=5.0;
724
725// // Domain length in y-direction
726// double Ly=1.0;
727
728// // Build and assign mesh
729// SimpleSpineMesh<SpineElement<QTaylorHoodElement<2> > >* mesh_pt =
730// new SimpleSpineMesh<SpineElement<QTaylorHoodElement<2> > >(Nx,Ny,Lx,Ly);
731
732// // Update *all* nodal positions
733// mesh_pt->node_update();
734
735// unsigned count=0;
736// ofstream some_file;
737// char filename[100];
738
739// // Number of plot points
740// unsigned npts=5;
741
742// // Output solution
743// sprintf(filename,"%s/mesh_update%i.dat",doc_info.directory().c_str(),
744// count);
745// count++;
746
747// // Loop over spines
748// unsigned n_node = mesh_pt->nnode();
749// for (unsigned inode=0;inode<n_node;inode++)
750// {
751// SpineNode* node_pt=dynamic_cast<SpineNode*>(mesh_pt->node_pt(inode));
752// node_pt->node_update();
753// // Output solution
754// some_file.open(filename);
755// sprintf(filename,"%s/mesh_update%i.dat",doc_info.directory().c_str(),
756// count);
757// count++;
758// mesh_pt->output(some_file,npts);
759// some_file.close();
760// }
761// }
762
763
764
765
766
767
768
769
770
771
772
773
774
775
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void doc_solution(DocInfo &doc_info)
Doc the solution.
void actions_after_newton_solve()
Update the after solve (empty)
double Ly
Width of channel.
SimpleSpineMesh< ELEMENT > * mesh_pt()
Overload access to mesh.
void actions_before_newton_solve()
Update the problem specs before solve. Update the nodal positions.
~ChannelSpineFlowProblem()
Destructor: (empty)
////////////////////////////////////////////////////////////////////
void update_spine_heights()
Update the spine heights according to the function specified by height_fct_pt().
virtual void spine_node_update(SpineNode *spine_node_pt)
General node update function implements pure virtual function defined in SpineMesh base class and per...
double(* HeightFctPt)(const double &x)
Function pointer to function, h(x), that may be used to specify the "height" of the domain,...
SimpleSpineMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &h, GeomObject *substrate_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction,...
double Default_height
Default height.
GeomObject * Substrate_pt
Pointer to GeomObject that specifies the "substrate" (the lower wall)
double height_fct(const double &x)
Height function – this is called by update_spine_heights() when spine heights are assigned.
HeightFctPt Height_fct_pt
Pointer to height function.
HeightFctPt & height_fct_pt()
Access function: Pointer to height function.
////////////////////////////////////////////////////////////////////
double A
Amplitude.
double L
Wavelength.
WavyWall(const double &l, const double &amplitude)
Constructor: Pass wavelength and amplitude.
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector to wavy wall.
////////////////////////////////////////////////////////////////////
double L
Length of indented region.
double A
Amplitude of indentation.
double X_indent_start
Start of indented region.
double L_total
Total length of domain.
double height(const double &x)
Height of domain.
double H
Undeformed height of domain.
int main()
Driver for channel flow problem with spine mesh.