hermite_element_quad_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// the mesh assembly functions for the hermite element quad mesh
27#ifndef OOMPH_HERMITE_QUAD_MESH_TEMPLATE_CC
28#define OOMPH_HERMITE_QUAD_MESH_TEMPLATE_CC
29
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36// oomph-lib includes
38
39
40namespace oomph
41{
42 //=============================================================================
43 /// sets the generalised position of the node (i.e. - x_i, dx_i/ds_0,
44 /// dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the
45 /// node from which its position can be determined.
46 //=============================================================================
47 template<class ELEMENT>
49 const unsigned& node_num_x, const unsigned& node_num_y, Node* node_pt)
50 {
51 // get the generalised macro element position of the node
52 DenseMatrix<double> m_gen(4, 2);
53 generalised_macro_element_position_of_node(node_num_x, node_num_y, m_gen);
54
55 // copy the macro element coordinates
56 Vector<double> node_macro_position(2);
57 node_macro_position[0] = m_gen(0, 0);
58 node_macro_position[1] = m_gen(0, 1);
59
60 // get the global position of the nodes
61 Vector<double> x_node(2);
62 Domain_pt->macro_element_pt(0)->macro_map(node_macro_position, x_node);
63
64 // get the jacobian of the mapping from macro to eulerian coordinates
65 DenseMatrix<double> jacobian_MX(2, 2);
66 Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian(
67 node_macro_position, jacobian_MX);
68
69
70 // get the jacobian2 of the mapping from macro to eulerian coordinates
71 DenseMatrix<double> jacobian2_MX(3, 2);
72 Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian2(
73 node_macro_position, jacobian2_MX);
74
75 // set x_0
76 node_pt->x_gen(0, 0) = x_node[0];
77 // set x_1
78 node_pt->x_gen(0, 1) = x_node[1];
79
80 // set dxi/dsj for i = 0,1 and j = 0,1
81 // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
82 for (unsigned i = 0; i < 2; i++)
83 {
84 for (unsigned j = 0; j < 2; j++)
85 {
86 node_pt->x_gen(j + 1, i) = m_gen(j + 1, 0) * jacobian_MX(0, i) +
87 m_gen(j + 1, 1) * jacobian_MX(1, i);
88 }
89 }
90
91 // set d2xi/ds0ds1 for i = 0,1
92 // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0
93 // + d2m1/ds0ds1*dxi/dm1
94 // + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
95 // + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0
96 for (unsigned i = 0; i < 2; i++)
97 {
98 node_pt->x_gen(3, i) =
99 m_gen(3, 0) * jacobian_MX(0, i) + m_gen(3, 1) * jacobian_MX(1, i) +
100 m_gen(1, 0) * (m_gen(2, 0) * jacobian2_MX(0, i) +
101 m_gen(2, 1) * jacobian2_MX(2, i)) +
102 m_gen(1, 1) *
103 (m_gen(2, 0) * jacobian2_MX(2, i) + m_gen(2, 1) * jacobian2_MX(1, i));
104 }
105 }
106
107
108 //=============================================================================
109 /// sets the generalised position of the node (i.e. - x_i, dx_i/ds_0,
110 /// dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the
111 /// node from which its position can be determined. Also sets coordinates
112 /// on boundary vector for the node to be the generalised position of the node
113 /// in macro element coordinates
114 //=============================================================================
115 template<class ELEMENT>
117 const unsigned& node_num_x,
118 const unsigned& node_num_y,
119 BoundaryNode<Node>* node_pt)
120 {
121 // get the generalised macro element position of the node
122 DenseMatrix<double> m_gen(4, 2, 0.0);
123 generalised_macro_element_position_of_node(node_num_x, node_num_y, m_gen);
124
125 // copy the macro element coordinates
126 Vector<double> node_macro_position(2);
127 node_macro_position[0] = m_gen(0, 0);
128 node_macro_position[1] = m_gen(0, 1);
129
130 // get the global position of the nodes
131 Vector<double> x_node(2);
132 Domain_pt->macro_element_pt(0)->macro_map(node_macro_position, x_node);
133
134 // get the jacobian of the mapping from macro to eulerian coordinates
135 DenseMatrix<double> jacobian_MX(2, 2);
136 Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian(
137 node_macro_position, jacobian_MX);
138
139 // get the jacobian2 of the mapping from macro to eulerian coordinates
140 DenseMatrix<double> jacobian2_MX(3, 2);
141 Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian2(
142 node_macro_position, jacobian2_MX);
143
144 // set x_0
145 node_pt->x_gen(0, 0) = x_node[0];
146 // set x_1
147 node_pt->x_gen(0, 1) = x_node[1];
148
149 // set dxi/dsj for i = 0,1 and j = 0,1
150 // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
151 for (unsigned i = 0; i < 2; i++)
152 {
153 for (unsigned j = 0; j < 2; j++)
154 {
155 node_pt->x_gen(j + 1, i) = m_gen(j + 1, 0) * jacobian_MX(0, i) +
156 m_gen(j + 1, 1) * jacobian_MX(1, i);
157 }
158 }
159
160 // set d2xi/ds0ds1 for i = 0,1
161 // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0
162 // + d2m1/ds0ds1*dxi/dm1
163 // + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
164 // + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0
165 for (unsigned i = 0; i < 2; i++)
166 {
167 node_pt->x_gen(3, i) =
168 m_gen(3, 0) * jacobian_MX(0, i) + m_gen(3, 1) * jacobian_MX(1, i) +
169 m_gen(1, 0) * (m_gen(2, 0) * jacobian2_MX(0, i) +
170 m_gen(2, 1) * jacobian2_MX(2, i)) +
171 m_gen(1, 1) *
172 (m_gen(2, 0) * jacobian2_MX(2, i) + m_gen(2, 1) * jacobian2_MX(1, i));
173 }
174
175
176 // pass generalised macro element position to vector to pass to boundary
177 // node
178 for (unsigned b = 0; b < 4; b++)
179 {
180 if (node_pt->is_on_boundary(b))
181 {
182 Vector<double> boundary_macro_position(2, 0);
183 boundary_macro_position[0] = m_gen(0, b % 2);
184 boundary_macro_position[1] = m_gen(1 + b % 2, b % 2);
185 node_pt->set_coordinates_on_boundary(b, boundary_macro_position);
186 }
187 }
188 }
189
190
191 //=============================================================================
192 /// computes the generalised position of the node at position
193 /// (node_num_x, node_num_y) in the macro element coordinate scheme.
194 /// index 0 of m_gen : 0 - m_i
195 /// 1 - dm_i/ds_0
196 /// 2 - dm_i/ds_1
197 /// 3 - d2m_i/ds_0ds_1 (where i is index 1 of m_gen)
198 //=============================================================================
199 template<class ELEMENT>
201 const unsigned& node_num_x,
202 const unsigned& node_num_y,
203 DenseMatrix<double>& m_gen)
204 {
205 // obtain position of node
206 Vector<double> s_macro_N(2);
207 macro_coordinate_position(node_num_x, node_num_y, s_macro_N);
208
209 // set position of node in macro element coordinates
210 m_gen(0, 0) = s_macro_N[0];
211 m_gen(0, 1) = s_macro_N[1];
212
213 // if on left hand edge
214 if (node_num_x == 0)
215 {
216 // first dm0ds0 and dm1ds0
217
218 // get position of right node
219 Vector<double> s_macro_R(2);
220 macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
221
222 // compute dm0ds0
223 m_gen(1, 0) = 0.5 * (s_macro_R[0] - s_macro_N[0]);
224
225 // compute dm1ds0
226 m_gen(1, 1) = 0.5 * (s_macro_R[1] - s_macro_N[1]);
227
228 // now d2m0/ds0ds1 and d2m1/ds0ds1
229
230 // if lower left hand corner
231 if (node_num_y == 0)
232 {
233 // get position of upper right node
234 Vector<double> s_macro_UR(2);
235 macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
236
237 // get position of upper node
238 Vector<double> s_macro_U(2);
239 macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
240
241 // compute d2m0/ds0ds1
242 m_gen(3, 0) =
243 0.5 * (0.5 * (s_macro_UR[0] - s_macro_U[0]) - m_gen(1, 0));
244
245 // compute d2m1/ds0ds1
246 m_gen(3, 1) =
247 0.5 * (0.5 * (s_macro_UR[1] - s_macro_U[1]) - m_gen(1, 1));
248 }
249
250 // else if upper left hand corner
251 else if (node_num_y == Nelement[1])
252 {
253 // get position of lower right node
254 Vector<double> s_macro_DR(2);
255 macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
256
257 // get position of lower node
258 Vector<double> s_macro_D(2);
259 macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
260
261 // compute d2m0/ds0ds1
262 m_gen(3, 0) =
263 0.5 * (m_gen(1, 0) - 0.5 * (s_macro_DR[0] - s_macro_D[0]));
264
265 // compute d2m0/ds0ds1
266 m_gen(3, 1) =
267 0.5 * (m_gen(1, 1) - 0.5 * (s_macro_DR[1] - s_macro_D[1]));
268 }
269
270 // else left hand edge (not corner)
271 else
272 {
273 // get position of lower right node
274 Vector<double> s_macro_DR(2);
275 macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
276
277 // get position of lower node
278 Vector<double> s_macro_D(2);
279 macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
280
281 // get position of upper right node
282 Vector<double> s_macro_UR(2);
283 macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
284
285 // get position of upper node
286 Vector<double> s_macro_U(2);
287 macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
288
289 // compute d2m0/ds0ds1
290 m_gen(3, 0) =
291 0.5 * std::min(m_gen(1, 0) - 0.5 * (s_macro_DR[0] - s_macro_D[0]),
292 0.5 * (s_macro_UR[0] - s_macro_U[0]) - m_gen(1, 0));
293
294 // compute d2m1/ds0ds1
295 m_gen(3, 1) =
296 0.5 * std::min(m_gen(1, 1) - 0.5 * (s_macro_DR[1] - s_macro_D[1]),
297 0.5 * (s_macro_UR[1] - s_macro_U[1]) - m_gen(1, 1));
298 }
299 }
300
301 // if on right hand edge
302 else if (node_num_x == Nelement[0])
303 {
304 // first dm0ds0 and dm1ds0
305
306 // get position of left node
307 Vector<double> s_macro_L(2);
308 macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
309
310 // compute dm0ds0
311 m_gen(1, 0) = 0.5 * (s_macro_N[0] - s_macro_L[0]);
312
313 // compute dm1ds0
314 m_gen(1, 1) = 0.5 * (s_macro_N[1] - s_macro_L[1]);
315
316 // now d2m0/ds0ds1 and d2m1/ds0ds1
317
318 // if lower right hand corner
319 if (node_num_y == 0)
320 {
321 // get position of upper left node
322 Vector<double> s_macro_UL(2);
323 macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
324
325 // get position of upper node
326 Vector<double> s_macro_U(2);
327 macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
328
329 // compute d2m0/ds0ds1
330 m_gen(3, 0) =
331 0.5 * (0.5 * (s_macro_U[0] - s_macro_UL[0]) - m_gen(1, 0));
332
333 // compute d2m1/ds0ds1
334 m_gen(3, 1) =
335 0.5 * (0.5 * (s_macro_U[1] - s_macro_UL[1]) - m_gen(1, 1));
336 }
337
338 // else if upper right hand corner
339 else if (node_num_y == Nelement[1])
340 {
341 // get position of lower left node
342 Vector<double> s_macro_DL(2);
343 macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
344
345 // get position of lower node
346 Vector<double> s_macro_D(2);
347 macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
348
349 // compute d2m0/ds0ds1
350 m_gen(3, 0) =
351 0.5 * (m_gen(1, 0) - 0.5 * (s_macro_D[0] - s_macro_DL[0]));
352
353 // compute d2m0/ds0ds1
354 m_gen(3, 1) =
355 0.5 * (m_gen(1, 1) - 0.5 * (s_macro_D[1] - s_macro_DL[1]));
356 }
357
358 // else left hand edge (not corner)
359 else
360 {
361 // get position of lower left node
362 Vector<double> s_macro_DL(2);
363 macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
364
365 // get position of lower node
366 Vector<double> s_macro_D(2);
367 macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
368
369 // get position of upper left node
370 Vector<double> s_macro_UL(2);
371 macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
372
373 // get position of upper node
374 Vector<double> s_macro_U(2);
375 macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
376
377 // compute d2m0/ds0ds1
378 m_gen(3, 0) =
379 0.5 * std::min(m_gen(1, 0) - 0.5 * (s_macro_D[0] - s_macro_DL[0]),
380 0.5 * (s_macro_U[0] - s_macro_UL[0]) - m_gen(1, 0));
381
382 // compute d2m0/ds0ds1
383 m_gen(3, 1) =
384 0.5 * std::min(m_gen(1, 1) - 0.5 * (s_macro_D[1] - s_macro_DL[1]),
385 0.5 * (s_macro_U[1] - s_macro_UL[1]) - m_gen(1, 1));
386 }
387 }
388 //
389 else
390 {
391 // get position of left node
392 Vector<double> s_macro_L(2);
393 macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
394
395 // get position of right node
396 Vector<double> s_macro_R(2);
397 macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
398
399 // compute dm0/ds0
400 m_gen(1, 0) = 0.5 * std::min(s_macro_N[0] - s_macro_L[0],
401 s_macro_R[0] - s_macro_N[0]);
402
403 // compute dm1/ds0
404 Vector<double> node_space(2);
405 node_space[0] = s_macro_N[1] - s_macro_L[1];
406 node_space[1] = s_macro_R[1] - s_macro_N[1];
407 // set nodal dof
408 if (node_space[0] > 0 && node_space[1] > 0)
409 m_gen(1, 1) = 0.5 * std::min(node_space[0], node_space[1]);
410 else if (node_space[0] < 0 && node_space[1] < 0)
411 m_gen(1, 1) = 0.5 * std::max(node_space[0], node_space[1]);
412 else
413 m_gen(1, 1) = 0;
414 }
415
416 // if node in lower row
417 if (node_num_y == 0)
418 {
419 // now dm0ds1 and dm1ds1
420
421 // get position of upper node
422 Vector<double> s_macro_U(2);
423 macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
424
425 // compute dm0ds1
426 m_gen(2, 0) = 0.5 * (s_macro_U[0] - s_macro_N[0]);
427
428 // compute dm1ds1
429 m_gen(2, 1) = 0.5 * (s_macro_U[1] - s_macro_N[1]);
430
431 // and if not corner node
432 if (node_num_x > 0 && node_num_x < Nelement[0])
433 {
434 // now dm0/ds0ds1 and dm1/ds0ds1
435
436 // get position of upper left node
437 Vector<double> s_macro_UL(2);
438 macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
439
440 // get position of left node
441 Vector<double> s_macro_L(2);
442 macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
443
444 // get position of right node
445 Vector<double> s_macro_R(2);
446 macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
447
448 // get position of upper right node
449 Vector<double> s_macro_UR(2);
450 macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
451
452 // set dm0/ds0ds1
453 m_gen(3, 0) =
454 0.5 * std::min(m_gen(2, 0) - 0.5 * (s_macro_UL[0] - s_macro_L[0]),
455 0.5 * (s_macro_UR[0] - s_macro_R[0]) - m_gen(2, 0));
456
457 // set dm1/ds0ds1
458 m_gen(3, 1) =
459 0.5 * std::min(m_gen(2, 1) - 0.5 * (s_macro_UL[1] - s_macro_L[1]),
460 0.5 * (s_macro_UR[1] - s_macro_R[1]) - m_gen(2, 1));
461 }
462 }
463
464 // else if node in upper row
465 else if (node_num_y == Nelement[1])
466 {
467 // now dm0ds1 and dm1ds1
468
469 // get position of lower node
470 Vector<double> s_macro_D(2);
471 macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
472
473 // compute dm0ds1
474 m_gen(2, 0) = 0.5 * (s_macro_N[0] - s_macro_D[0]);
475
476 // compute dm1ds1
477 m_gen(2, 1) = 0.5 * (s_macro_N[1] - s_macro_D[1]);
478
479 // and if not corner node
480 if (node_num_x > 0 && node_num_x < Nelement[0])
481 {
482 // now dm0/ds0ds1 and dm1/ds0ds1
483
484 // get position of lower left node
485 Vector<double> s_macro_DL(2);
486 macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
487
488 // get position of left node
489 Vector<double> s_macro_L(2);
490 macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
491
492 // get position of right node
493 Vector<double> s_macro_R(2);
494 macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
495
496 // get position of lower right node
497 Vector<double> s_macro_DR(2);
498 macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
499
500 // set dm0/ds0ds1
501 m_gen(3, 0) =
502 0.5 * std::min(m_gen(2, 0) - 0.5 * (s_macro_L[0] - s_macro_DL[0]),
503 0.5 * (s_macro_R[0] - s_macro_DR[0]) - m_gen(2, 0));
504
505 // set dm1/ds0ds1
506 m_gen(3, 1) =
507 0.5 * std::min(m_gen(2, 1) - 0.5 * (s_macro_L[1] - s_macro_DL[1]),
508 0.5 * (s_macro_R[1] - s_macro_DR[1]) - m_gen(2, 1));
509 }
510 }
511 else
512 {
513 // get position of upper node
514 Vector<double> s_macro_U(2);
515 macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
516
517 // get position of lower node
518 Vector<double> s_macro_D(2);
519 macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
520
521 // compute dm0/ds1
522 Vector<double> node_space(2);
523 node_space[0] = s_macro_N[0] - s_macro_D[0];
524 node_space[1] = s_macro_U[0] - s_macro_N[0];
525 // set nodal coordinate
526 if (node_space[0] > 0 && node_space[1] > 0)
527 m_gen(2, 0) = 0.5 * std::min(node_space[0], node_space[1]);
528 else if (node_space[0] < 0 && node_space[1] < 0)
529 m_gen(2, 0) = 0.5 * std::max(node_space[0], node_space[1]);
530 else
531 m_gen(2, 0) = 0;
532
533 // compute dm1/ds1
534 m_gen(2, 1) = 0.5 * std::min(s_macro_N[1] - s_macro_D[1],
535 s_macro_U[1] - s_macro_N[1]);
536
537 // for interior nodes
538 if (node_num_x > 0 && node_num_x < Nelement[0])
539 {
540 // get position of left node
541 Vector<double> s_macro_L(2);
542 macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
543
544 // get position of upper left node
545 Vector<double> s_macro_UL(2);
546 macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
547
548 // get position of upper right node
549 Vector<double> s_macro_UR(2);
550 macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
551
552 // get position of right node
553 Vector<double> s_macro_R(2);
554 macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
555
556 // get position of lower right node
557 Vector<double> s_macro_DR(2);
558 macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
559
560 // get position of lower left node
561 Vector<double> s_macro_DL(2);
562 macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
563
564 Vector<double> node_space(2);
565 // comute dm0/ds0ds1 wrt to node above and below
566 node_space[0] =
567 m_gen(1, 0) - 0.5 * std::min(s_macro_D[0] - s_macro_DL[0],
568 s_macro_DR[0] - s_macro_D[0]);
569 node_space[1] = 0.5 * std::min(s_macro_U[0] - s_macro_UL[0],
570 s_macro_UR[0] - s_macro_U[0]) -
571 m_gen(1, 0);
572 // set nodal dof
573 if (node_space[0] > 0 && node_space[1] > 0)
574 m_gen(3, 0) = 0.5 * std::min(node_space[0], node_space[1]);
575 else if (node_space[0] < 0 && node_space[1] < 0)
576 m_gen(3, 0) = 0.5 * std::max(node_space[0], node_space[1]);
577 else
578 m_gen(3, 0) = 0;
579
580 // comute dm1/ds0ds1 wrt node left and right
581 node_space[0] =
582 m_gen(2, 1) - 0.5 * std::min(s_macro_L[0] - s_macro_DL[0],
583 s_macro_UL[0] - s_macro_L[0]);
584 node_space[1] = 0.5 * std::min(s_macro_R[0] - s_macro_DR[0],
585 s_macro_UR[0] - s_macro_R[0]) -
586 m_gen(2, 1);
587 // set nodal dof
588 if (node_space[0] > 0 && node_space[1] > 0)
589 m_gen(3, 1) = 0.5 * std::min(node_space[0], node_space[1]);
590 else if (node_space[0] < 0 && node_space[1] < 0)
591 m_gen(3, 1) = 0.5 * std::max(node_space[0], node_space[1]);
592 else
593 m_gen(3, 1) = 0;
594 }
595 }
596 }
597
598
599 //=============================================================================
600 /// Generic mesh construction function to build the mesh
601 //=============================================================================
602 template<class ELEMENT>
604 {
605 // Mesh can only be built with 2D QHermiteElements.
606 MeshChecker::assert_geometric_element<QHermiteElementBase, ELEMENT>(2);
607
608 // Set the number of boundaries
609 set_nboundary(4);
610
611 // Allocate the store for the elements
612 Element_pt.resize(Nelement[0] * Nelement[1]);
613
614 // Allocate the memory for the first element
615 Element_pt[0] = new ELEMENT;
616 finite_element_pt(0)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
617
618 // Can now allocate the store for the nodes
619 Node_pt.resize((1 + Nelement[0]) * (1 + Nelement[1]));
620
621 // Set up geometrical data
622 unsigned long node_count = 0;
623
624 // Now assign the topology
625 // Boundaries are numbered 0 1 2 3 from the bottom proceeding anticlockwise
626 // Pinned value are denoted by an integer value 1
627 // Thus if a node is on two boundaries, ORing the values of the
628 // boundary conditions will give the most restrictive case (pinning)
629
630
631 // we first create the lowest row of elements
632 // ##########################################
633 // ##########################################
634
635
636 // FIRST ELEMENT - lower left hand corner
637 // **************************************
638
639
640 // LOWER LEFT HAND NODE
641
642 // Set the corner node
643 // Allocate memory for the node
644 Node_pt[node_count] =
645 finite_element_pt(0)->construct_boundary_node(0, time_stepper_pt);
646
647 // Push the node back onto boundaries
648 add_boundary_node(0, Node_pt[node_count]);
649 add_boundary_node(3, Node_pt[node_count]);
650
651 // set the position of the boundary node
652 set_position_of_boundary_node(
653 0, 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
654
655 // Increment the node number
656 node_count++;
657
658 // LOWER RIGHT HAND NODE
659
660 // Allocate memory for the node
661 Node_pt[node_count] =
662 finite_element_pt(0)->construct_boundary_node(1, time_stepper_pt);
663
664 // Push the node back onto boundaries
665 add_boundary_node(0, Node_pt[node_count]);
666
667 // If we only have one column then the RHS node is on the right boundary
668 if (Nelement[0] == 1)
669 {
670 add_boundary_node(1, Node_pt[node_count]);
671 }
672
673 // set the position of the node
674 set_position_of_boundary_node(
675 1, 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
676
677 // Increment the node number
678 node_count++;
679
680 // UPPER LEFT HAND NODE
681
682 // Allocate memory for the node
683 Node_pt[node_count] =
684 finite_element_pt(0)->construct_boundary_node(2, time_stepper_pt);
685
686 // Push the node back onto boundaries
687 add_boundary_node(3, Node_pt[node_count]);
688
689 // If we only have one row, then the top node is on the top boundary
690 if (Nelement[1] == 1)
691 {
692 add_boundary_node(2, Node_pt[node_count]);
693 }
694
695 // set the position of the boundary node
696 set_position_of_boundary_node(
697 0, 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
698
699 // Increment the node number
700 node_count++;
701
702 // UPPER RIGHT NODE
703
704 // Allocate the memory for the node
705 Node_pt[node_count] =
706 finite_element_pt(0)->construct_node(3, time_stepper_pt);
707
708 // If we only have one column then the RHS node is on the right boundary
709 if (Nelement[0] == 1)
710 {
711 add_boundary_node(1, Node_pt[node_count]);
712 }
713 // If we only have one row, then the top node is on the top boundary
714 if (Nelement[1] == 1)
715 {
716 add_boundary_node(2, Node_pt[node_count]);
717 }
718
719 // if the node is a boundary node
720 if (Nelement[0] == 1 || Nelement[1] == 1)
721 {
722 // set the position of the boundary node
723 set_position_of_boundary_node(
724 1, 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
725 }
726 else
727 {
728 // set the position of the node
729 set_position_of_node(1, 1, Node_pt[node_count]);
730 }
731
732 // Increment the node number
733 node_count++;
734
735 // END OF FIRST ELEMENT
736
737
738 // CENTRE OF FIRST ROW OF ELEMENTS
739 // *******************************
740
741
742 // Now loop over the first row of elements, apart from final element
743 for (unsigned j = 1; j < (Nelement[0] - 1); j++)
744 {
745 // Allocate memory for new element
746 Element_pt[j] = new ELEMENT;
747 finite_element_pt(j)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
748
749 // LOWER LEFT NODE
750
751 // lower left hand node column of nodes is same as lower right hand node
752 // in neighbouring element
753 finite_element_pt(j)->node_pt(0) = finite_element_pt(j - 1)->node_pt(1);
754
755 // LOWER RIGHT NODE
756
757 // Allocate memory for the nodes
758 Node_pt[node_count] =
759 finite_element_pt(j)->construct_boundary_node(1, time_stepper_pt);
760
761 // Push the node back onto boundaries
762 add_boundary_node(0, Node_pt[node_count]);
763
764 // set the position of the boundary node
765 set_position_of_boundary_node(
766 j + 1, 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
767
768 // Increment the node number
769 node_count++;
770
771 // UPPER LEFT NODE
772
773 // First column of nodes is same as neighbouring element
774 finite_element_pt(j)->node_pt(2) = finite_element_pt(j - 1)->node_pt(3);
775
776 // UPPER RIGHT NODE
777
778 // Allocate memory for the nodes
779 Node_pt[node_count] =
780 finite_element_pt(j)->construct_node(3, time_stepper_pt);
781
782 // If we only have one row, then the top node is on the top boundary
783 if (Nelement[0] == 1)
784 {
785 add_boundary_node(2, Node_pt[node_count]);
786 }
787
788 // set the position of the (boundary) node
789 if (Nelement[0] == 1)
790 set_position_of_boundary_node(
791 j + 1, 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
792 else
793 set_position_of_node(j + 1, 1, Node_pt[node_count]);
794
795 // Increment the node number
796 node_count++;
797 }
798
799 // FINAL ELEMENT IN FIRST ROW (lower right corner element)
800 // **************************
801
802
803 // Only allocate if there is more than one element in the row
804 if (Nelement[0] > 1)
805 {
806 // Allocate memory for element
807 Element_pt[Nelement[0] - 1] = new ELEMENT;
808 finite_element_pt(Nelement[0] - 1)
809 ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
810
811 // LOWER LEFT NODE
812
813 // First column of nodes is same as neighbouring element
814 finite_element_pt(Nelement[0] - 1)->node_pt(0) =
815 finite_element_pt(Nelement[0] - 2)->node_pt(1);
816
817 // LOWER RIGHT NODE
818
819 // If we have an Xperiodic mesh then the final node is the first node in
820 // the row
821 if (Xperiodic == true)
822 {
823 // Note that this is periodic in the x-direction
824 finite_element_pt(Nelement[0] - 1)->node_pt(1) =
825 finite_element_pt(0)->node_pt(0);
826 }
827 // Otherwise it's a new final node
828 else
829 {
830 // Allocate memory for the node
831 Node_pt[node_count] = finite_element_pt(Nelement[0] - 1)
832 ->construct_boundary_node(1, time_stepper_pt);
833 }
834
835 // Push the node back onto boundaries
836 add_boundary_node(0, Node_pt[node_count]);
837 add_boundary_node(1, Node_pt[node_count]);
838
839 // set the position of the boundary node
840 set_position_of_boundary_node(
841 Nelement[0], 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
842
843 // if not periodic mesh
844 if (Xperiodic == false)
845 {
846 node_count++;
847 }
848
849 // UPPER LEFT NODE
850
851 // same as upper right node in element to left
852 finite_element_pt(Nelement[0] - 1)->node_pt(2) =
853 finite_element_pt(Nelement[0] - 2)->node_pt(3);
854
855 // If we only have one row, then the top node is on the top boundary
856 if (Nelement[1] == 1)
857 {
858 add_boundary_node(2, Node_pt[node_count]);
859 }
860
861 // set the position of the (boundary) node
862 if (Nelement[1] == 1)
863 set_position_of_boundary_node(
864 Nelement[0] - 1,
865 1,
866 dynamic_cast<BoundaryNode<Node>*>(
867 finite_element_pt(Nelement[0] - 2)->node_pt(3)));
868
869 // UPPER RIGHT NODE
870
871 // If we have an Xperiodic mesh then the nodes in the final column are
872 // those in the first column
873 if (Xperiodic == true)
874 {
875 // Note that these are periodic in the x-direction
876 finite_element_pt(Nelement[0] - 1)->node_pt(3) =
877 finite_element_pt(0)->node_pt(2);
878 }
879 // Otherwise we need new nodes for the final column
880 else
881 {
882 // Allocate memory for node
883 Node_pt[node_count] = finite_element_pt(Nelement[0] - 1)
884 ->construct_boundary_node(3, time_stepper_pt);
885 }
886
887 // If we only have one row, then the top node is on the top boundary
888 if (Nelement[1] == 1)
889 {
890 add_boundary_node(2, Node_pt[node_count]);
891 }
892
893 // boundary 1 node
894 add_boundary_node(1, Node_pt[node_count]);
895
896 // set the position of the boundary node
897 set_position_of_boundary_node(
898 Nelement[0], 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
899
900 // Increment the node number
901 if (Xperiodic == false)
902 {
903 // Increment the node number
904 node_count++;
905 }
906 }
907
908 // now create all remaining central rows
909 // #####################################
910 // #####################################
911
912
913 // Loop over remaining element rows in the mesh
914 for (unsigned i = 1; i < (Nelement[1] - 1); i++)
915 {
916 // set the first element in the row
917 // ********************************
918
919
920 // Allocate memory for element
921 Element_pt[Nelement[0] * i] = new ELEMENT;
922 finite_element_pt(Nelement[0] * i)
923 ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
924
925 // The first row of nodes is copied from the element below
926 for (unsigned l = 0; l < 2; l++)
927 {
928 finite_element_pt(Nelement[0] * i)->node_pt(l) =
929 finite_element_pt(Nelement[0] * (i - 1))->node_pt(2 + l);
930 }
931
932 // UPPER LEFT HAND NODE
933
934 // Allocate memory for node
935 Node_pt[node_count] = finite_element_pt(Nelement[0] * i)
936 ->construct_boundary_node(2, time_stepper_pt);
937
938 // Push the node back onto boundaries
939 add_boundary_node(3, Node_pt[node_count]);
940
941 // set the position of the boundary node
942 set_position_of_boundary_node(
943 0, i + 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
944
945 // Increment the node number
946 node_count++;
947
948 // UPPER RIGHT HAND NODE
949
950 // If we only have one column, the node could be on the boundary
951 if (Nelement[0] == 1)
952 {
953 Node_pt[node_count] = finite_element_pt(Nelement[0] * i)
954 ->construct_boundary_node(3, time_stepper_pt);
955 }
956 else
957 {
958 Node_pt[node_count] = finite_element_pt(Nelement[0] * i)
959 ->construct_node(3, time_stepper_pt);
960 }
961
962 // If we only have one column then the RHS node is on the
963 // right boundary
964 if (Nelement[0] == 1)
965 {
966 add_boundary_node(1, Node_pt[node_count]);
967 }
968
969 // set the position of the (boundary) node
970 if (Nelement[0] == 1)
971 set_position_of_boundary_node(
972 1, i + 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
973 else
974 set_position_of_node(1, i + 1, Node_pt[node_count]);
975
976 // Increment the node number
977 node_count++;
978
979
980 // loop over central elements in row
981 // *********************************
982
983 // Now loop over the rest of the elements in the row, apart from the last
984 for (unsigned j = 1; j < (Nelement[0] - 1); j++)
985 {
986 // Allocate memory for new element
987 Element_pt[Nelement[0] * i + j] = new ELEMENT;
988 finite_element_pt(Nelement[0] * i + j)
989 ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
990
991 // LOWER LEFT AND RIGHT NODES
992
993 // The first row is copied from the lower element
994 for (unsigned l = 0; l < 2; l++)
995 {
996 finite_element_pt(Nelement[0] * i + j)->node_pt(l) =
997 finite_element_pt(Nelement[0] * (i - 1) + j)->node_pt(2 + l);
998 }
999
1000 // UPPER LEFT NODE
1001
1002 // First column is same as neighbouring element
1003 finite_element_pt(Nelement[0] * i + j)->node_pt(2) =
1004 finite_element_pt(Nelement[0] * i + (j - 1))->node_pt(3);
1005
1006 // UPPER RIGHT NODE
1007
1008 // Allocate memory for the nodes
1009 Node_pt[node_count] = finite_element_pt(Nelement[0] * i + j)
1010 ->construct_node(3, time_stepper_pt);
1011
1012 // set position of node
1013 set_position_of_node(j + 1, i + 1, Node_pt[node_count]);
1014
1015 // Increment the node number
1016 node_count++;
1017 }
1018
1019
1020 // final element in row
1021 // ********************
1022
1023
1024 // Only if there is more than one column
1025 if (Nelement[0] > 1)
1026 {
1027 // Allocate memory for element
1028 Element_pt[Nelement[0] * i + Nelement[0] - 1] = new ELEMENT;
1029 finite_element_pt(Nelement[0] * i + Nelement[0] - 1)
1030 ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1031
1032 // LOWER LEFT AND RIGHT NODES
1033
1034 // The first row is copied from the lower element
1035 for (unsigned l = 0; l < 2; l++)
1036 {
1037 finite_element_pt(Nelement[0] * i + Nelement[0] - 1)->node_pt(l) =
1038 finite_element_pt(Nelement[0] * (i - 1) + Nelement[0] - 1)
1039 ->node_pt(2 + l);
1040 }
1041
1042 // UPPER LEFT NODE
1043
1044 // First node is same as neighbouring element
1045 finite_element_pt(Nelement[0] * i + Nelement[0] - 1)->node_pt(2) =
1046 finite_element_pt(Nelement[0] * i + Nelement[0] - 2)->node_pt(3);
1047
1048 // UPPER RIGHT NODE
1049
1050 // If we have an Xperiodic mesh then this node is the same
1051 // as the first node
1052 if (Xperiodic == true)
1053 {
1054 // Allocate memory for node, periodic in x-direction
1055 finite_element_pt(Nelement[0] * i + Nelement[0] - 1)->node_pt(3) =
1056 finite_element_pt(Nelement[0] * i)->node_pt(2);
1057 }
1058 // Otherwise allocate memory for a new node
1059 else
1060 {
1061 Node_pt[node_count] =
1062 finite_element_pt(Nelement[0] * i + Nelement[0] - 1)
1063 ->construct_boundary_node(3, time_stepper_pt);
1064 }
1065
1066 // Push the node back onto boundaries
1067 add_boundary_node(1, Node_pt[node_count]);
1068
1069 // set position of boundary node
1070 set_position_of_boundary_node(
1071 Nelement[0],
1072 i + 1,
1073 dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1074
1075 // Increment the node number
1076 node_count++;
1077 }
1078 }
1079
1080
1081 // final element row
1082 // #################
1083 // #################
1084
1085
1086 // only if there is more than one row
1087 if (Nelement[1] > 1)
1088 {
1089 // first element in upper row (upper left corner)
1090 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1091
1092 // Allocate memory for element
1093 Element_pt[Nelement[0] * (Nelement[1] - 1)] = new ELEMENT;
1094 finite_element_pt(Nelement[0] * (Nelement[1] - 1))
1095 ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1096
1097 // LOWER LEFT AND LOWER RIGHT NODES
1098
1099 // The first row of nodes is copied from the element below
1100 for (unsigned l = 0; l < 2; l++)
1101 {
1102 finite_element_pt(Nelement[0] * (Nelement[1] - 1))->node_pt(l) =
1103 finite_element_pt(Nelement[0] * (Nelement[1] - 2))->node_pt(2 + l);
1104 }
1105
1106 // UPPER LEFT NODE
1107
1108 // Allocate memory for node
1109 Node_pt[node_count] = finite_element_pt(Nelement[0] * (Nelement[1] - 1))
1110 ->construct_boundary_node(2, time_stepper_pt);
1111
1112 // Push the node back onto boundaries
1113 add_boundary_node(2, Node_pt[node_count]);
1114 add_boundary_node(3, Node_pt[node_count]);
1115
1116 // set position of boundary node
1117 set_position_of_boundary_node(
1118 0, Nelement[1], dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1119
1120 // Increment the node number
1121 node_count++;
1122
1123 // UPPER RIGHT NODE
1124
1125 // Allocate memory for the node
1126 Node_pt[node_count] = finite_element_pt(Nelement[0] * (Nelement[1] - 1))
1127 ->construct_boundary_node(3, time_stepper_pt);
1128
1129 // Push the node back onto boundaries
1130 add_boundary_node(2, Node_pt[node_count]);
1131
1132 // set position of boundary node
1133 set_position_of_boundary_node(
1134 1, Nelement[1], dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1135
1136 // Increment the node number
1137 node_count++;
1138
1139
1140 // loop over central elements of upper element row
1141 // ***********************************************
1142
1143
1144 // Now loop over the rest of the elements in the row, apart from the last
1145 for (unsigned j = 1; j < (Nelement[0] - 1); j++)
1146 {
1147 // Allocate memory for element
1148 Element_pt[Nelement[0] * (Nelement[1] - 1) + j] = new ELEMENT;
1149 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)
1150 ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1151
1152 // LOWER LEFT AND LOWER RIGHT NODES
1153
1154 // The first row is copied from the lower element
1155 for (unsigned l = 0; l < 2; l++)
1156 {
1157 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)->node_pt(l) =
1158 finite_element_pt(Nelement[0] * (Nelement[1] - 2) + j)
1159 ->node_pt(2 + l);
1160 }
1161
1162 // UPPER LEFT NODE
1163
1164 // First column is same as neighbouring element
1165 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)->node_pt(2) =
1166 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + (j - 1))
1167 ->node_pt(3);
1168
1169 // UPPER RIGHT NODE
1170
1171 // Allocate memory for node
1172 Node_pt[node_count] =
1173 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)
1174 ->construct_boundary_node(3, time_stepper_pt);
1175
1176 // Push the node back onto boundaries
1177 add_boundary_node(2, Node_pt[node_count]);
1178
1179 // set position of boundary node
1180 set_position_of_boundary_node(
1181 j + 1,
1182 Nelement[1],
1183 dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1184
1185 // Increment the node number
1186 node_count++;
1187
1188 } // End of loop over central elements in row
1189
1190
1191 // final element in row (upper right corner)
1192 // *****************************************
1193
1194 // Only if there is more than one column
1195 if (Nelement[0] > 1)
1196 {
1197 // Allocate memory for element
1198 Element_pt[Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1] =
1199 new ELEMENT;
1200 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1201 ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1202
1203 // LOWER LEFT AND LOWER RIGHT NODES
1204
1205 // The first row is copied from the lower element
1206 for (unsigned l = 0; l < 2; l++)
1207 {
1208 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1209 ->node_pt(l) =
1210 finite_element_pt(Nelement[0] * (Nelement[1] - 2) + Nelement[0] - 1)
1211 ->node_pt(2 + l);
1212 }
1213
1214 // UPPER LEFT NODE
1215
1216 // First column is same as neighbouring element
1217 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1218 ->node_pt(2) =
1219 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 2)
1220 ->node_pt(3);
1221
1222 // UPPER RIGHT NODE
1223
1224 // If we have an Xperiodic mesh, the node must be copied from
1225 // the first column
1226 if (Xperiodic == true)
1227 {
1228 // Allocate memory for node, periodic in x-direction
1229 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1230 ->node_pt(3) =
1231 finite_element_pt(Nelement[0] * (Nelement[1] - 1))->node_pt(2);
1232 }
1233 // Otherwise, allocate new memory for node
1234 else
1235 {
1236 Node_pt[node_count] =
1237 finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1238 ->construct_boundary_node(3, time_stepper_pt);
1239 }
1240
1241 // Push the node back onto boundaries
1242 add_boundary_node(1, Node_pt[node_count]);
1243 add_boundary_node(2, Node_pt[node_count]);
1244
1245 // set position of boundary node
1246 set_position_of_boundary_node(
1247 Nelement[0],
1248 Nelement[1],
1249 dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1250
1251 // Increment the node number
1252 node_count++;
1253 }
1254 }
1255
1256 // Setup boundary element lookup schemes
1257 setup_boundary_element_info();
1258 }
1259
1260
1261 //=============================================================================
1262 /// Setup lookup schemes which establish which elements are located
1263 /// next to which boundaries (Doc to outfile if it's open). Specific version
1264 /// for HermiteQuadMesh to ensure that the order of the elements in
1265 /// Boundary_element_pt matches the actual order along the boundary. This is
1266 /// required when hijacking the BiharmonicElement to apply the
1267 /// BiharmonicFluidBoundaryElement in
1268 /// BiharmonicFluidProblem::impose_traction_free_edge(...)
1269 //================================================================
1270 template<class ELEMENT>
1272 std::ostream& outfile)
1273 {
1274 bool doc = false;
1275 if (outfile) doc = true;
1276
1277 // Number of boundaries
1278 unsigned nbound = nboundary();
1279
1280 // Wipe/allocate storage for arrays
1281 Boundary_element_pt.clear();
1282 Face_index_at_boundary.clear();
1283 Boundary_element_pt.resize(nbound);
1284 Face_index_at_boundary.resize(nbound);
1285
1286 // Temporary vector of sets of pointers to elements on the boundaries:
1287 Vector<Vector<FiniteElement*>> vector_of_boundary_element_pt;
1288 vector_of_boundary_element_pt.resize(nbound);
1289
1290 // Matrix map for working out the fixed local coord for elements on boundary
1292
1293
1294 // Loop over elements
1295 //-------------------
1296 unsigned nel = nelement();
1297 for (unsigned e = 0; e < nel; e++)
1298 {
1299 // Get pointer to element
1300 FiniteElement* fe_pt = finite_element_pt(e);
1301
1302 if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
1303
1304 // Only include 2D elements! Some meshes contain interface elements too.
1305 if (fe_pt->dim() == 2)
1306 {
1307 // Loop over the element's nodes and find out which boundaries they're
1308 // on
1309 // ----------------------------------------------------------------------
1310 unsigned nnode_1d = fe_pt->nnode_1d();
1311
1312 // Loop over nodes in order
1313 for (unsigned i0 = 0; i0 < nnode_1d; i0++)
1314 {
1315 for (unsigned i1 = 0; i1 < nnode_1d; i1++)
1316 {
1317 // Local node number
1318 unsigned j = i0 + i1 * nnode_1d;
1319
1320 // Get pointer to vector of boundaries that this
1321 // node lives on
1322 std::set<unsigned>* boundaries_pt = 0;
1323 fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
1324
1325 // If the node lives on some boundaries....
1326 if (boundaries_pt != 0)
1327 {
1328 // Loop over boundaries
1329 // unsigned nbound=(*boundaries_pt).size();
1330 for (std::set<unsigned>::iterator it = boundaries_pt->begin();
1331 it != boundaries_pt->end();
1332 ++it)
1333 {
1334 // Add pointer to finite element to set for the appropriate
1335 // boundary -- storage in set makes sure we don't count elements
1336 // multiple times
1337 unsigned temp_size = vector_of_boundary_element_pt[*it].size();
1338 bool temp_flag = true;
1339 for (unsigned temp_i = 0; temp_i < temp_size; temp_i++)
1340 {
1341 if (vector_of_boundary_element_pt[*it][temp_i] == fe_pt)
1342 temp_flag = false;
1343 }
1344 if (temp_flag)
1345 vector_of_boundary_element_pt[*it].push_back(fe_pt);
1346
1347 // For the current element/boundary combination, create
1348 // a vector that stores an indicator which element boundaries
1349 // the node is located (boundary_identifier=-/+1 for nodes
1350 // on the left/right boundary; boundary_identifier=-/+2 for
1351 // nodes on the lower/upper boundary. We determine these indices
1352 // for all corner nodes of the element and add them to a vector
1353 // to a vector. This allows us to decide which face of the
1354 // element coincides with the boundary since the (quad!) element
1355 // must have exactly two corner nodes on the boundary.
1356 if (boundary_identifier(*it, fe_pt) == 0)
1357 {
1358 boundary_identifier(*it, fe_pt) = new Vector<int>;
1359 }
1360
1361 // Are we at a corner node?
1362 if (((i0 == 0) || (i0 == nnode_1d - 1)) &&
1363 ((i1 == 0) || (i1 == nnode_1d - 1)))
1364 {
1365 // Create index to represent position relative to s_0
1366 (*boundary_identifier(*it, fe_pt))
1367 .push_back(1 * (2 * i0 / (nnode_1d - 1) - 1));
1368
1369 // Create index to represent position relative to s_1
1370 (*boundary_identifier(*it, fe_pt))
1371 .push_back(2 * (2 * i1 / (nnode_1d - 1) - 1));
1372 }
1373 }
1374 }
1375 // else
1376 // {
1377 // oomph_info << "...does not live on any boundaries " <<
1378 // std::endl;
1379 // }
1380 }
1381 }
1382 }
1383 // else
1384 //{
1385 // oomph_info << "Element " << e << " does not qualify" << std::endl;
1386 //}
1387 }
1388
1389
1390 // Now copy everything across into permanent arrays
1391 //-------------------------------------------------
1392
1393 // Note: vector_of_boundary_element_pt contains all elements
1394 // that have (at least) one corner node on a boundary -- can't copy
1395 // them across into Boundary_element_pt
1396 // yet because some of them might have only one node on the
1397 // the boundary in which case they don't qualify as
1398 // boundary elements!
1399
1400 // Loop over boundaries
1401 //---------------------
1402 for (unsigned i = 0; i < nbound; i++)
1403 {
1404 // Number of elements on this boundary (currently stored in a set)
1405 unsigned nel = vector_of_boundary_element_pt[i].size();
1406
1407 // Allocate storage for the coordinate identifiers
1408 Face_index_at_boundary[i].resize(nel);
1409
1410 // Loop over elements that have at least one corner node on this boundary
1411 //-----------------------------------------------------------------------
1412 unsigned e_count = 0;
1414 for (IT it = vector_of_boundary_element_pt[i].begin();
1415 it != vector_of_boundary_element_pt[i].end();
1416 it++)
1417 {
1418 // Recover pointer to element
1419 FiniteElement* fe_pt = *it;
1420
1421 // Initialise count for boundary identiers (-2,-1,1,2)
1422 std::map<int, unsigned> count;
1423
1424 // Loop over coordinates
1425 for (unsigned ii = 0; ii < 2; ii++)
1426 {
1427 // Loop over upper/lower end of coordinates
1428 for (int sign = -1; sign < 3; sign += 2)
1429 {
1430 count[(ii + 1) * sign] = 0;
1431 }
1432 }
1433
1434 // Loop over boundary indicators for this element/boundary
1435 unsigned n_indicators = (*boundary_identifier(i, fe_pt)).size();
1436 for (unsigned j = 0; j < n_indicators; j++)
1437 {
1438 count[(*boundary_identifier(i, fe_pt))[j]]++;
1439 }
1440 delete boundary_identifier(i, fe_pt);
1441
1442 // Determine the correct boundary indicator by checking that it
1443 // occurs twice (since two corner nodes of the element's boundary
1444 // need to be located on the domain boundary
1445 int indicator = -10;
1446
1447 // Check that we're finding exactly one boundary indicator
1448 bool found = false;
1449
1450 // Loop over coordinates
1451 for (unsigned ii = 0; ii < 2; ii++)
1452 {
1453 // Loop over upper/lower end of coordinates
1454 for (int sign = -1; sign < 3; sign += 2)
1455 {
1456 if (count[(ii + 1) * sign] == 2)
1457 {
1458 // Check that we haven't found multiple boundaries
1459 if (found)
1460 {
1461 std::string error_message =
1462 "Trouble: Multiple boundary identifiers!\n";
1463 error_message += "Elements should only have at most 2 corner ";
1464 error_message += "nodes on any one boundary.\n";
1465
1466 throw OomphLibError(error_message,
1467 OOMPH_CURRENT_FUNCTION,
1468 OOMPH_EXCEPTION_LOCATION);
1469 }
1470 found = true;
1471 indicator = (ii + 1) * sign;
1472 }
1473 }
1474 }
1475
1476 // Element has exactly two corner nodes on boundary
1477 if (found)
1478 {
1479 // Add to permanent storage
1480 Boundary_element_pt[i].push_back(fe_pt);
1481
1482 // Now convert boundary indicator into information required
1483 // for FaceElements
1484 switch (indicator)
1485 {
1486 case -2:
1487
1488 // s_1 is fixed at -1.0:
1489 Face_index_at_boundary[i][e_count] = -2;
1490 break;
1491
1492 case -1:
1493
1494 // s_0 is fixed at -1.0:
1495 Face_index_at_boundary[i][e_count] = -1;
1496 break;
1497
1498
1499 case 1:
1500
1501 // s_0 is fixed at 1.0:
1502 Face_index_at_boundary[i][e_count] = 1;
1503 break;
1504
1505 case 2:
1506
1507 // s_1 is fixed at 1.0:
1508 Face_index_at_boundary[i][e_count] = 2;
1509 break;
1510
1511 default:
1512
1513 throw OomphLibError("Never get here",
1514 OOMPH_CURRENT_FUNCTION,
1515 OOMPH_EXCEPTION_LOCATION);
1516 }
1517
1518 // Increment counter
1519 e_count++;
1520 }
1521 }
1522 }
1523
1524
1525 // Doc?
1526 //-----
1527 if (doc)
1528 {
1529 // Loop over boundaries
1530 for (unsigned i = 0; i < nbound; i++)
1531 {
1532 unsigned nel = Boundary_element_pt[i].size();
1533 outfile << "Boundary: " << i << " is adjacent to " << nel << " elements"
1534 << std::endl;
1535
1536 // Loop over elements on given boundary
1537 for (unsigned e = 0; e < nel; e++)
1538 {
1539 FiniteElement* fe_pt = Boundary_element_pt[i][e];
1540 outfile << "Boundary element:" << fe_pt
1541 << " Face index of element along boundary is "
1542 << Face_index_at_boundary[i][e] << std::endl;
1543 }
1544 }
1545 }
1546
1547
1548 // Lookup scheme has now been setup yet
1549 Lookup_for_elements_next_boundary_is_setup = true;
1550 }
1551} // namespace oomph
1552
1553#endif
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh....
Definition: nodes.h:2242
void set_coordinates_on_boundary(const unsigned &b, const Vector< double > &boundary_zeta)
Set the vector of coordinates on mesh boundary b Final overload.
Definition: nodes.h:2574
bool is_on_boundary() const
Test whether the node lies on a boundary Final overload.
Definition: nodes.h:2532
A general Finite Element class.
Definition: elements.h:1313
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
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
void generalised_macro_element_position_of_node(const unsigned &node_num_x, const unsigned &node_num_y, DenseMatrix< double > &m_gen)
computes the generalised position of the node at position (node_num_x, node_num_y) in the macro eleme...
virtual void setup_boundary_element_info()
Setup lookup schemes which establish whic elements are located next to mesh's boundaries (wrapper to ...
void set_position_of_boundary_node(const unsigned &node_num_x, const unsigned &node_num_y, BoundaryNode< Node > *node_pt)
sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, dx_i/ds_1 & d2x_i/ds_0ds_1 for i = ...
void set_position_of_node(const unsigned &node_num_x, const unsigned &node_num_y, Node *node_pt)
sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, dx_i/ds_1 & d2x_i/ds_0ds_1 for i = ...
virtual void build_mesh(TimeStepper *time_stepper_pt)
Generic mesh construction function to build the mesh.
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1365
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
Definition: nodes.h:1126
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...