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// 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 #ifndef OOMPH_BACKWARD_STEP_MESH_TEMPLATE_CC
27 #define OOMPH_BACKWARD_STEP_MESH_TEMPLATE_CC
28 
29 #include <set>
30 #include <map>
32 
33 namespace oomph
34 {
35  //=================================================================
36  /// Actual build function. Pass overall number of elements in the
37  /// horizontal and vertical directions, nx and ny, and the corresponding
38  /// dimensions, lx and ly. nx_cut_out and ny_cut_out elements
39  /// are cut out from the lower right corner to create the
40  /// (reversed) backward step geometry. Timestepper defaults
41  /// to Steady.
42  //=================================================================
43  template<class ELEMENT>
45  const unsigned& ny,
46  const unsigned& nx_cut_out,
47  const unsigned& ny_cut_out,
48  const double& lx,
49  const double& ly)
50  {
51  // Mesh can only be built with 2D Qelements.
52  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
53 
54  // By default nobody's claiming any nodes
55  std::map<Node*, bool> keep;
56 
57  // Get elements outside lower right block
58  Vector<FiniteElement*> el_retained_pt;
59  Vector<FiniteElement*> el_killed_pt;
60  for (unsigned i = 0; i < nx; i++)
61  {
62  for (unsigned j = 0; j < ny; j++)
63  {
64  FiniteElement* el_pt = this->finite_element_pt(i + nx * j);
65  if ((i > (nx_cut_out - 1)) && (j < ny_cut_out))
66  {
67  el_killed_pt.push_back(el_pt);
68  }
69  else
70  {
71  el_retained_pt.push_back(el_pt);
72  unsigned nnod_el = el_pt->nnode();
73  for (unsigned jj = 0; jj < nnod_el; jj++)
74  {
75  keep[el_pt->node_pt(jj)] = true;
76  }
77  }
78  }
79  }
80 
81 
82  // By default nobody's claiming the nodes; also store old
83  // boundary ids
84  std::map<Node*, std::set<unsigned>> boundaries;
85  unsigned nnod = this->nnode();
86  for (unsigned j = 0; j < nnod; j++)
87  {
88  std::set<unsigned>* aux_pt = 0;
89  this->node_pt(j)->get_boundaries_pt(aux_pt);
90  if (aux_pt != 0)
91  {
92  boundaries[this->node_pt(j)] = (*aux_pt);
93  }
94  }
95 
96  // Remove information about boundary nodes
97  this->remove_boundary_nodes();
98 
99  // Reset number of boundaries
100  this->set_nboundary(6);
101 
102  // Kill superfluous nodes
103  Vector<Node*> node_backup_pt(this->Node_pt);
104  this->Node_pt.clear();
105  for (unsigned j = 0; j < nnod; j++)
106  {
107  Node* nod_pt = node_backup_pt[j];
108  if (keep[nod_pt])
109  {
110  this->Node_pt.push_back(nod_pt);
111  std::set<unsigned> aux = boundaries[nod_pt];
112  for (std::set<unsigned>::iterator it = boundaries[nod_pt].begin();
113  it != boundaries[nod_pt].end();
114  it++)
115  {
116  unsigned b = (*it);
117  if (b > 0) b += 2;
118  this->add_boundary_node(b, nod_pt);
119  }
120  }
121  else
122  {
123  delete node_backup_pt[j];
124  }
125  }
126 
127  // Add nodes to new boundary 1
128  unsigned i = nx_cut_out - 1;
129  for (unsigned j = 0; j < ny_cut_out; j++)
130  {
131  FiniteElement* el_pt = this->finite_element_pt(i + nx * j);
132  unsigned nnod_1d = el_pt->nnode_1d();
133  for (unsigned jj = 0; jj < nnod_1d; jj++)
134  {
135  unsigned jnod = (nnod_1d - 1) + jj * nnod_1d;
136  if (!(el_pt->node_pt(jnod)->is_on_boundary()))
137  {
138  Node* nod_pt = el_pt->node_pt(jnod);
139  this->convert_to_boundary_node(nod_pt);
140  }
141  this->add_boundary_node(1, el_pt->node_pt(jnod));
142  }
143  }
144 
145  // Add nodes to new boundary 2
146  unsigned j = ny_cut_out;
147  for (unsigned i = nx_cut_out; i < nx; i++)
148  {
149  FiniteElement* el_pt = this->finite_element_pt(i + nx * j);
150  unsigned nnod_1d = el_pt->nnode_1d();
151  for (unsigned jj = 0; jj < nnod_1d; jj++)
152  {
153  unsigned jnod = jj;
154  if (!(el_pt->node_pt(jnod)->is_on_boundary()))
155  {
156  Node* nod_pt = el_pt->node_pt(jnod);
157  this->convert_to_boundary_node(nod_pt);
158  }
159  this->add_boundary_node(2, el_pt->node_pt(jnod));
160  }
161  }
162 
163 
164  // Kill redundant elements
165  this->Element_pt.clear();
166  unsigned n_retained = el_retained_pt.size();
167  for (unsigned e = 0; e < n_retained; e++)
168  {
169  this->Element_pt.push_back(el_retained_pt[e]);
170  }
171  unsigned n_killed = el_killed_pt.size();
172  for (unsigned e = 0; e < n_killed; e++)
173  {
174  delete el_killed_pt[e];
175  }
176 
177  // Re-setup lookup scheme that establishes which elements are located
178  // on the mesh boundaries
179  this->setup_boundary_element_info();
180  }
181 
182 } // namespace oomph
183 
184 #endif
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.
////////////////////////////////////////////////////////////////////// //////////////////////////////...