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-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#ifndef OOMPH_BACKWARD_STEP_MESH_TEMPLATE_CC
27#define OOMPH_BACKWARD_STEP_MESH_TEMPLATE_CC
28
29#include <set>
30#include <map>
32
33namespace 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
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
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.
A general Finite Element class.
Definition: elements.h:1313
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:2218
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...