backward_step_mesh.template.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 #ifndef OOMPH_BACKWARD_STEP_MESH_TEMPLATE_CC
31 #define OOMPH_BACKWARD_STEP_MESH_TEMPLATE_CC
32 
33 #include<set>
34 #include<map>
36 
37 namespace oomph
38 {
39 
40 
41 
42 
43 //=================================================================
44  /// \short Actual build function. Pass overall number of elements in the
45  /// horizontal and vertical directions, nx and ny, and the corresponding
46  /// dimensions, lx and ly. nx_cut_out and ny_cut_out elements
47  /// are cut out from the lower right corner to create the
48  /// (reversed) backward step geometry. Timestepper defaults
49  /// to Steady.
50 //=================================================================
51  template<class ELEMENT>
53  const unsigned &nx,
54  const unsigned &ny,
55  const unsigned& nx_cut_out,
56  const unsigned& ny_cut_out,
57  const double &lx,
58  const double &ly)
59  {
60 
61  // Mesh can only be built with 2D Qelements.
62  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
63 
64  // By default nobody's claiming any nodes
65  std::map<Node*,bool> keep;
66 
67  // Get elements outside lower right block
68  Vector<FiniteElement*> el_retained_pt;
69  Vector<FiniteElement*> el_killed_pt;
70  for (unsigned i=0;i<nx;i++)
71  {
72  for (unsigned j=0;j<ny;j++)
73  {
74  FiniteElement* el_pt=this->finite_element_pt(i+nx*j);
75  if ((i>(nx_cut_out-1))&&(j<ny_cut_out))
76  {
77  el_killed_pt.push_back(el_pt);
78  }
79  else
80  {
81  el_retained_pt.push_back(el_pt);
82  unsigned nnod_el=el_pt->nnode();
83  for (unsigned jj=0;jj<nnod_el;jj++)
84  {
85  keep[el_pt->node_pt(jj)]=true;
86  }
87  }
88  }
89  }
90 
91 
92  // By default nobody's claiming the nodes; also store old
93  // boundary ids
94  std::map<Node*,std::set<unsigned> > boundaries;
95  unsigned nnod=this->nnode();
96  for (unsigned j=0;j<nnod;j++)
97  {
98  std::set<unsigned>* aux_pt=0;
99  this->node_pt(j)->get_boundaries_pt(aux_pt);
100  if (aux_pt!=0)
101  {
102  boundaries[this->node_pt(j)]=(*aux_pt);
103  }
104  }
105 
106  // Remove information about boundary nodes
107  this->remove_boundary_nodes();
108 
109  // Reset number of boundaries
110  this->set_nboundary(6);
111 
112  // Kill superfluous nodes
113  Vector<Node*> node_backup_pt(this->Node_pt);
114  this->Node_pt.clear();
115  for (unsigned j=0;j<nnod;j++)
116  {
117  Node* nod_pt=node_backup_pt[j];
118  if (keep[nod_pt])
119  {
120  this->Node_pt.push_back(nod_pt);
121  std::set<unsigned> aux=boundaries[nod_pt];
122  for (std::set<unsigned>::iterator it=boundaries[nod_pt].begin();
123  it!=boundaries[nod_pt].end();it++)
124  {
125  unsigned b=(*it);
126  if (b>0) b+=2;
127  this->add_boundary_node(b,nod_pt);
128  }
129  }
130  else
131  {
132  delete node_backup_pt[j];
133  }
134  }
135 
136  // Add nodes to new boundary 1
137  unsigned i=nx_cut_out-1;
138  for (unsigned j=0;j<ny_cut_out;j++)
139  {
140  FiniteElement* el_pt=this->finite_element_pt(i+nx*j);
141  unsigned nnod_1d=el_pt->nnode_1d();
142  for (unsigned jj=0;jj<nnod_1d;jj++)
143  {
144  unsigned jnod=(nnod_1d-1)+jj*nnod_1d;
145  if (!(el_pt->node_pt(jnod)->is_on_boundary()))
146  {
147  Node* nod_pt=el_pt->node_pt(jnod);
148  this->convert_to_boundary_node(nod_pt);
149  }
150  this->add_boundary_node(1,el_pt->node_pt(jnod));
151  }
152  }
153 
154  // Add nodes to new boundary 2
155  unsigned j=ny_cut_out;
156  for (unsigned i=nx_cut_out;i<nx;i++)
157  {
158  FiniteElement* el_pt=this->finite_element_pt(i+nx*j);
159  unsigned nnod_1d=el_pt->nnode_1d();
160  for (unsigned jj=0;jj<nnod_1d;jj++)
161  {
162  unsigned jnod=jj;
163  if (!(el_pt->node_pt(jnod)->is_on_boundary()))
164  {
165  Node* nod_pt=el_pt->node_pt(jnod);
166  this->convert_to_boundary_node(nod_pt);
167  }
168  this->add_boundary_node(2,el_pt->node_pt(jnod));
169  }
170  }
171 
172 
173 
174  // Kill redundant elements
175  this->Element_pt.clear();
176  unsigned n_retained=el_retained_pt.size();
177  for (unsigned e=0;e<n_retained;e++)
178  {
179  this->Element_pt.push_back(el_retained_pt[e]);
180  }
181  unsigned n_killed=el_killed_pt.size();
182  for (unsigned e=0;e<n_killed;e++)
183  {
184  delete el_killed_pt[e];
185  }
186 
187  // Re-setup lookup scheme that establishes which elements are located
188  // on the mesh boundaries
189  this->setup_boundary_element_info();
190  }
191 
192 }
193 
194 #endif
cstr elem_len * i
Definition: cfortran.h:607
A general Finite Element class.
Definition: elements.h:1271
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary, so return false. This will be overloaded by BoundaryNodes.
Definition: nodes.h:1293
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
e
Definition: cfortran.h:575
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2139
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void build_mesh(const unsigned &nx, const unsigned &ny, const unsigned &nx_cut_out, const unsigned &ny_cut_out, const double &lx, const double &ly)
Actual build function.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134