Toggle navigation
Documentation
Big picture
The finite element method
The data structure
Not-so-quick guide
Optimisation
Order of action functions
Example codes and tutorials
List of example codes and tutorials
Meshing
Solvers
MPI parallel processing
Post-processing/visualisation
Other
Change log
Creating documentation
Coding conventions
Index
FAQ
Installation
Installation guide
Copyright
About
People
Contact/Get involved
Publications
Acknowledgements
Picture show
Go
src
meshes
simple_rectangular_tri_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_SIMPLE_RECTANGULAR_TRIMESH_TEMPLATE_CC
27
#define OOMPH_SIMPLE_RECTANGULAR_TRIMESH_TEMPLATE_CC
28
29
30
#include <algorithm>
31
32
#include "../generic/Telements.h"
33
34
// Simple 2D triangle mesh class
35
#include "
simple_rectangular_tri_mesh.template.h
"
36
37
38
namespace
oomph
39
{
40
//====================================================================
41
/// Constructor for simple 2D triangular mesh class:
42
///
43
/// n_x : number of elements in the x direction
44
///
45
/// n_y : number of elements in the y direction
46
///
47
/// l_x : length in the x direction
48
///
49
/// l_y : length in the y direction
50
///
51
/// Ordering of elements: 'lower left' to 'lower right' then 'upwards'
52
//====================================================================
53
template
<
class
ELEMENT>
54
SimpleRectangularTriMesh<ELEMENT>::SimpleRectangularTriMesh
(
55
const
unsigned
&
n_x
,
56
const
unsigned
&
n_y
,
57
const
double
&
l_x
,
58
const
double
& l_y,
59
TimeStepper
*
time_stepper_pt
)
60
: Nx(
n_x
), Ny(
n_y
), Lx(
l_x
), Ly(l_y)
61
{
62
using namespace
MathematicalConstants
;
63
64
// Mesh can only be built with 2D Telements.
65
MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
66
67
// Set number of boundaries
68
this->
set_nboundary
(4);
69
70
// Allocate the store for the elements
71
unsigned
n_element
= (
Nx
) * (
Ny
)*2;
72
Element_pt
.resize(
n_element
, 0);
73
74
// Create first element
75
Element_pt
[0] =
new
ELEMENT
;
76
77
// Currently this mesh only works for 3 and 6 noded triangles
78
if
((
finite_element_pt
(0)->
nnode_1d
() != 2) &&
79
(
finite_element_pt
(0)->nnode_1d() != 3) &&
80
(
finite_element_pt
(0)->nnode_1d() != 4))
81
{
82
throw
OomphLibError
(
83
"Currently this mesh only works for 3, 6 & 10-noded triangles"
,
84
OOMPH_CURRENT_FUNCTION
,
85
OOMPH_EXCEPTION_LOCATION
);
86
}
87
88
unsigned
n_node
= 0;
89
// Unless nnode_1d returned as !=2 default no extra nodes
90
unsigned
additional_nnode
= 0;
91
92
// Allocate storage if no extra nodes
93
if
(
finite_element_pt
(0)->
nnode_1d
() == 2)
94
{
95
n_node
= (
Nx
+ 1) * (
Ny
+ 1);
96
}
97
98
if
(
finite_element_pt
(0)->nnode_1d() == 3)
99
{
100
additional_nnode
= 3;
101
// Allocate storage for nodes (including extra)
102
n_node
= (2 *
Nx
+ 1) * (2 *
Ny
+ 1);
103
}
104
105
if
(
finite_element_pt
(0)->nnode_1d() == 4)
106
{
107
additional_nnode
= 7;
108
// Allocate storage for notes (including extra)
109
n_node
= (3 *
Nx
+ 1) * (3 *
Ny
+ 1);
110
}
111
112
Node_pt
.resize(
n_node
);
113
114
// Set up geometrical data
115
//------------------------
116
unsigned
long
node_count
= 0;
117
unsigned
long
element_count
= 0;
118
double
xinit
= 0.0,
yinit
= 0.0;
119
// Set the values of the increments
120
double
xstep
=
Lx
/ (
Nx
);
121
double
ystep
=
Ly
/ (
Ny
);
122
123
124
// FIRST NODE (lower left corner)
125
//-------------------------------
126
127
// Set the corner node
128
129
// Allocate memory for the corner node of the first element
130
//(which is not created loop as all subsequent bottom corners already exist)
131
Node_pt
[
node_count
] =
132
finite_element_pt
(0)->construct_node(1,
time_stepper_pt
);
133
134
// Set the pointer from the element
135
finite_element_pt
(0)->node_pt(1) =
Node_pt
[
node_count
];
136
137
// Don't increment node_count yet as we still need its containing box
138
//(position and boundaries for first node are allocated in the main loop)
139
140
// CREATE THE ELEMENTS (each has 3 local nodes)
141
//--------------------------------------------------
142
// Elements are created two at a time, the first being lower triangular
143
// and the second upper triangular, so that they form a containing box.
144
// Local nodes are numbered anti-clockwise with node_pt(1) being the
145
// right-angle corner of the triangle.
146
// The global node number, node_count, is considered to define
147
// a containing box, with node_count defined as the node number
148
// of the bottom left corner of the box.
149
for
(
unsigned
j
= 0;
j
<
Ny
+ 1; ++
j
)
150
{
151
for
(
unsigned
i
= 0;
i
<
Nx
+ 1; ++
i
)
152
{
153
// CURRENT BOX
154
//(nodes on RHS or top edge of domain do not define a box)
155
if
(
j
<
Ny
&&
i
<
Nx
)
156
{
157
// Create two elements
158
//------------------------------
159
if
(
element_count
!= 0)
// 0th element already exists
160
{
161
Element_pt
[
element_count
] =
new
ELEMENT
;
162
}
163
164
Element_pt
[
element_count
+ 1] =
new
ELEMENT
;
165
166
167
// Allocate memory for nodes in the current box
168
//--------------------------------------------
169
// If node on LHS of domain, allocate the top left corner node
170
if
(
i
== 0)
171
{
172
Node_pt
[
node_count
+
Nx
+ 1] =
173
finite_element_pt
(
element_count
)
174
->construct_node(0,
time_stepper_pt
);
175
}
176
177
// If on bottom row, allocate the bottom right corner node
178
if
(
j
== 0)
179
{
180
Node_pt
[
node_count
+ 1] =
finite_element_pt
(
element_count
)
181
->construct_node(2,
time_stepper_pt
);
182
}
183
184
// Always need to allocate top right corner node
185
Node_pt
[
node_count
+
Nx
+ 2] =
finite_element_pt
(
element_count
+ 1)
186
->construct_node(1,
time_stepper_pt
);
187
188
// Bottom left corner node is already allocated
189
190
191
// Set the pointers from the elements
192
//----------------------------------
193
finite_element_pt
(
element_count
)->node_pt(0) =
194
Node_pt
[
node_count
+
Nx
+ 1];
195
finite_element_pt
(
element_count
)->node_pt(1) =
Node_pt
[
node_count
];
196
finite_element_pt
(
element_count
)->node_pt(2) =
197
Node_pt
[
node_count
+ 1];
198
finite_element_pt
(
element_count
+ 1)->node_pt(0) =
199
Node_pt
[
node_count
+ 1];
200
finite_element_pt
(
element_count
+ 1)->node_pt(1) =
201
Node_pt
[
node_count
+
Nx
+ 2];
202
finite_element_pt
(
element_count
+ 1)->node_pt(2) =
203
Node_pt
[
node_count
+
Nx
+ 1];
204
205
element_count
+= 2;
206
}
207
208
// CURRENT (GLOBAL) NODE
209
// Set the position
210
Node_pt
[
node_count
]->x(0) =
xinit
+
i
*
xstep
;
211
Node_pt
[
node_count
]->x(1) =
yinit
+
j
*
ystep
;
212
213
// Add node to any relevant boundaries
214
if
(
j
== 0)
215
{
216
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
217
add_boundary_node
(0,
Node_pt
[
node_count
]);
218
}
219
if
(
j
==
Ny
)
220
{
221
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
222
add_boundary_node
(2,
Node_pt
[
node_count
]);
223
}
224
if
(
i
== 0)
225
{
226
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
227
add_boundary_node
(3,
Node_pt
[
node_count
]);
228
}
229
if
(
i
==
Nx
)
230
{
231
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
232
add_boundary_node
(1,
Node_pt
[
node_count
]);
233
}
234
235
// Increment counter
236
node_count
++;
237
}
238
}
239
240
241
if
(
additional_nnode
== 3)
242
{
243
// Reset element counter to allow looping over existing elements
244
// to add extra nodes.
245
// Note that the node_count is not reset so no entries are overwritten
246
element_count
= 0;
247
for
(
unsigned
j
= 0;
j
<
Ny
+ 1; ++
j
)
248
{
249
// Note: i counter reduced by 1 since i axis runs through middle of
250
// elements on x-axis
251
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
252
{
253
// The additional nodes will be added in stages for ease of
254
// application. Note: local numbering follows the anti-clockwise
255
// pattern, node 3 halfway between nodes 0-1, 4 between 1-2 and the
256
// 5th local node between 2-0.
257
258
// Restricted to stop creation of additional nodes outside the mesh
259
if
(
j
<
Ny
)
260
{
261
// If on the bottom row middle-node at bottom needs allocating
262
if
(
j
== 0)
263
{
264
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
265
->construct_node(4,
time_stepper_pt
);
266
}
267
268
// Due to directions of iteration node at middle of top box edge
269
// (in next element) always needs allocating
270
Node_pt
[
node_count
+
Nx
] =
finite_element_pt
(
element_count
+ 1)
271
->construct_node(4,
time_stepper_pt
);
272
273
// Set pointers to the top/bottom middle nodes
274
// Bottom node
275
finite_element_pt
(
element_count
)->node_pt(4) =
Node_pt
[
node_count
];
276
// Top node
277
finite_element_pt
(
element_count
+ 1)->node_pt(4) =
278
Node_pt
[
node_count
+
Nx
];
279
280
// Increase the element counter to add top/bottom nodes to
281
// next pair of element on next pass
282
element_count
+= 2;
283
}
// End extra top/bottom node construction
284
285
// Set position of current (Global) Node
286
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 0.5) *
xstep
;
287
Node_pt
[
node_count
]->x(1) =
yinit
+
j
*
ystep
;
288
289
// Add node to any applicable boundaries (node 4's can only be top
290
// or bottom boundaries)
291
if
(
j
== 0)
292
{
293
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
294
add_boundary_node
(0,
Node_pt
[
node_count
]);
295
}
296
if
(
j
==
Ny
)
297
{
298
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
299
add_boundary_node
(2,
Node_pt
[
node_count
]);
300
}
301
302
// Update node_count
303
node_count
++;
304
}
305
}
306
307
// Next stage of additional node implementation involes the middle left
308
// and right nodes (local number 3 on each triangle)
309
310
// Element counter reset for second loop over existing elements
311
element_count
= 0;
312
// Note: j counter reduced by 1 since j axis runs through middle of
313
// elements on y-axis
314
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
315
{
316
for
(
unsigned
i
= 0;
i
<
Nx
+ 1; ++
i
)
317
{
318
if
(
j
<
Ny
&&
i
<
Nx
)
319
{
320
// If on the left hand boundary the node at middle of the left
321
// side needs allocating
322
if
(
i
== 0)
323
{
324
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
325
->construct_node(3,
time_stepper_pt
);
326
}
327
328
// The mid node on the right hand side always needs constructing
329
// within the bounds of the mesh
330
Node_pt
[
node_count
+ 1] =
finite_element_pt
(
element_count
+ 1)
331
->construct_node(3,
time_stepper_pt
);
332
333
// Set pointers from the elements to new nodes
334
finite_element_pt
(
element_count
)->node_pt(3) =
Node_pt
[
node_count
];
335
finite_element_pt
(
element_count
+ 1)->node_pt(3) =
336
Node_pt
[
node_count
+ 1];
337
338
// Increase element_count by 2
339
element_count
+= 2;
340
}
// End extra left/right node construction
341
342
// Set position of current (Global) Node
343
Node_pt
[
node_count
]->x(0) =
xinit
+
i
*
xstep
;
344
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 0.5) *
ystep
;
345
346
// Add node to any applicable boundaries again - only be left/right
347
if
(
i
== 0)
348
{
349
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
350
add_boundary_node
(3,
Node_pt
[
node_count
]);
351
}
352
if
(
i
==
Nx
)
353
{
354
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
355
add_boundary_node
(1,
Node_pt
[
node_count
]);
356
}
357
358
// Update node_count
359
node_count
++;
360
}
361
}
362
363
// Final stage of inserting extra nodes - inclusion of the local
364
// number 5 (middle of hypot. edge)
365
366
element_count
= 0;
367
// Note: both i,j counters reduced by 1 since j axis runs through middle
368
// of elements in both x,y
369
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
370
{
371
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
372
{
373
// The middle node always needs constructing
374
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
375
->construct_node(5,
time_stepper_pt
);
376
377
// Set pointers from the elements to new nodes
378
finite_element_pt
(
element_count
)->node_pt(5) =
Node_pt
[
node_count
];
379
finite_element_pt
(
element_count
+ 1)->node_pt(5) =
380
Node_pt
[
node_count
];
381
382
// Increase element_count by 2
383
element_count
+= 2;
384
// End extra left/right node construction
385
386
// Set position of current (Global) Node
387
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 0.5) *
xstep
;
388
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 0.5) *
ystep
;
389
390
// All nodes are internal in this structure so no boundaries
391
// applicable
392
393
// Update node_count
394
node_count
++;
395
}
396
}
397
}
398
// End of extra nodes for 6 noded trianglur elements
399
400
401
if
(
additional_nnode
== 7)
402
{
403
// Reset element counter to allow looping over existing elements
404
// to add extra nodes.
405
// Note that the node_count is not reset so no entries are overwritten
406
element_count
= 0;
407
for
(
unsigned
j
= 0;
j
<
Ny
+ 1; ++
j
)
408
{
409
// Note: i counter reduced by 1 for implementation of lower element
410
// node 5 and upper element node 6's.
411
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
412
{
413
// The additional nodes will be added in stages for ease of
414
// application. Note: local numbering follows the anti-clockwise
415
// pattern, nodes 3 and 4 equidistant between nodes 0-1, 5 and 6
416
// between 1-2, 7 and 8 between 2-0 and last node 9 located
417
// (internally) at the centroid of the triangle.
418
419
// Restricted to stop creation of additional nodes outside the mesh
420
if
(
j
<
Ny
)
421
{
422
// If on the bottom row middle-left-node at bottom needs allocating
423
if
(
j
== 0)
424
{
425
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
426
->construct_node(5,
time_stepper_pt
);
427
}
428
429
// Due to directions of iteration node at middle left of top box
430
// edge (in next element) always needs allocating
431
Node_pt
[
node_count
+
Nx
] =
finite_element_pt
(
element_count
+ 1)
432
->construct_node(6,
time_stepper_pt
);
433
434
// Set pointers to the top/bottom middle nodes
435
// Bottom node
436
finite_element_pt
(
element_count
)->node_pt(5) =
Node_pt
[
node_count
];
437
// Top node
438
finite_element_pt
(
element_count
+ 1)->node_pt(6) =
439
Node_pt
[
node_count
+
Nx
];
440
441
// Increase the element counter to add top/bottom nodes to
442
// next pair of element on next pass
443
element_count
+= 2;
444
}
// End extra top/bottom node construction
445
446
// Set position of current (Global) Node
447
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 1.0 / 3.0) *
xstep
;
448
Node_pt
[
node_count
]->x(1) =
yinit
+
j
*
ystep
;
449
450
// Add node to any applicable boundaries (exactly as previous)
451
if
(
j
== 0)
452
{
453
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
454
add_boundary_node
(0,
Node_pt
[
node_count
]);
455
}
456
if
(
j
==
Ny
)
457
{
458
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
459
add_boundary_node
(2,
Node_pt
[
node_count
]);
460
}
461
462
// Update node_count
463
node_count
++;
464
}
465
}
466
467
468
// Next stage: bottom-middle-right node (node 6 in lower tri.el.)
469
// and top-middle-right node (node 5 in upper tri.el.)
470
471
element_count
= 0;
472
for
(
unsigned
j
= 0;
j
<
Ny
+ 1; ++
j
)
473
{
474
// Note: i counter as for above 5/6's
475
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
476
{
477
// Restricted to stop creation of additional nodes outside the mesh
478
if
(
j
<
Ny
)
479
{
480
// If on the bottom row middle-right-node at bottom needs allocating
481
if
(
j
== 0)
482
{
483
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
484
->construct_node(6,
time_stepper_pt
);
485
}
486
487
// Due to directions of iteration node at middle left of top box
488
// edge (in next element) always needs allocating
489
Node_pt
[
node_count
+
Nx
] =
finite_element_pt
(
element_count
+ 1)
490
->construct_node(5,
time_stepper_pt
);
491
492
// Set pointers to the top/bottom middle nodes
493
// Bottom node
494
finite_element_pt
(
element_count
)->node_pt(6) =
Node_pt
[
node_count
];
495
// Top node
496
finite_element_pt
(
element_count
+ 1)->node_pt(5) =
497
Node_pt
[
node_count
+
Nx
];
498
499
// Increase the element counter to add top/bottom nodes to
500
// next pair of element on next pass
501
element_count
+= 2;
502
}
// End extra top/bottom node construction
503
504
// Set position of current (Global) Node
505
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 2.0 / 3.0) *
xstep
;
506
Node_pt
[
node_count
]->x(1) =
yinit
+
j
*
ystep
;
507
508
// Add node to any applicable boundaries (exactly as previous)
509
if
(
j
== 0)
510
{
511
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
512
add_boundary_node
(0,
Node_pt
[
node_count
]);
513
}
514
if
(
j
==
Ny
)
515
{
516
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
517
add_boundary_node
(2,
Node_pt
[
node_count
]);
518
}
519
520
// Update node_count
521
node_count
++;
522
}
523
}
524
525
526
// Next stage of additional node implementation involes node 4 on lower
527
// tri. el. and node 3 on upper tri. el.
528
529
// Element counter reset for next loop over existing elements
530
element_count
= 0;
531
// Note: j counter reduced by 1 similar to others above
532
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
533
{
534
for
(
unsigned
i
= 0;
i
<
Nx
+ 1; ++
i
)
535
{
536
if
(
j
<
Ny
&&
i
<
Nx
)
537
{
538
// If on the left hand boundary the lower middle node of the left
539
// side needs allocating
540
if
(
i
== 0)
541
{
542
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
543
->construct_node(4,
time_stepper_pt
);
544
}
545
546
// The mid node on the right hand side always needs constructing
547
// within the bounds of the mesh
548
Node_pt
[
node_count
+ 1] =
finite_element_pt
(
element_count
+ 1)
549
->construct_node(3,
time_stepper_pt
);
550
551
// Set pointers from the elements to new nodes
552
finite_element_pt
(
element_count
)->node_pt(4) =
Node_pt
[
node_count
];
553
finite_element_pt
(
element_count
+ 1)->node_pt(3) =
554
Node_pt
[
node_count
+ 1];
555
556
// Increase element_count by 2
557
element_count
+= 2;
558
}
// End extra left/right node construction
559
560
// Set position of current (Global) Node
561
Node_pt
[
node_count
]->x(0) =
xinit
+
i
*
xstep
;
562
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 1.0 / 3.0) *
ystep
;
563
564
// Add node to any applicable boundaries again - only be left/right
565
if
(
i
== 0)
566
{
567
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
568
add_boundary_node
(3,
Node_pt
[
node_count
]);
569
}
570
if
(
i
==
Nx
)
571
{
572
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
573
add_boundary_node
(1,
Node_pt
[
node_count
]);
574
}
575
576
// Update node_count
577
node_count
++;
578
}
579
}
580
581
582
// Next stage of additional node implementation involes node 3 on lower
583
// tri. el. and node 4 on upper tri. el.
584
585
// Element counter reset for next loop over existing elements
586
element_count
= 0;
587
// Note: j counter reduced by 1 similar to others above
588
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
589
{
590
for
(
unsigned
i
= 0;
i
<
Nx
+ 1; ++
i
)
591
{
592
if
(
j
<
Ny
&&
i
<
Nx
)
593
{
594
// If on the left hand boundary the lower middle node of the left
595
// side needs allocating
596
if
(
i
== 0)
597
{
598
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
599
->construct_node(3,
time_stepper_pt
);
600
}
601
602
// The mid node on the right hand side always needs constructing
603
// within the bounds of the mesh
604
Node_pt
[
node_count
+ 1] =
finite_element_pt
(
element_count
+ 1)
605
->construct_node(4,
time_stepper_pt
);
606
607
// Set pointers from the elements to new nodes
608
finite_element_pt
(
element_count
)->node_pt(3) =
Node_pt
[
node_count
];
609
finite_element_pt
(
element_count
+ 1)->node_pt(4) =
610
Node_pt
[
node_count
+ 1];
611
612
// Increase element_count by 2
613
element_count
+= 2;
614
}
// End extra left/right node construction
615
616
// Set position of current (Global) Node
617
Node_pt
[
node_count
]->x(0) =
xinit
+
i
*
xstep
;
618
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 2.0 / 3.0) *
ystep
;
619
620
// Add node to any applicable boundaries again - only be left/right
621
if
(
i
== 0)
622
{
623
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
624
add_boundary_node
(3,
Node_pt
[
node_count
]);
625
}
626
if
(
i
==
Nx
)
627
{
628
this->
convert_to_boundary_node
(
Node_pt
[
node_count
]);
629
add_boundary_node
(1,
Node_pt
[
node_count
]);
630
}
631
632
// Update node_count
633
node_count
++;
634
}
635
}
636
637
638
// Next is the lower tri. el. totally internal node with local number 9
639
element_count
= 0;
640
// Note: both i,j counters reduced by 1
641
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
642
{
643
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
644
{
645
// The middle node always needs constructing
646
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
647
->construct_node(9,
time_stepper_pt
);
648
649
// Set pointers from the elements to new nodes
650
finite_element_pt
(
element_count
)->node_pt(9) =
Node_pt
[
node_count
];
651
652
// Increase element_count by 2
653
element_count
+= 2;
654
655
// Set position of current (Global) Node
656
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 1.0 / 3.0) *
xstep
;
657
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 1.0 / 3.0) *
ystep
;
658
659
// All nodes are internal in this structure so no boundaries
660
// applicable
661
662
// Update node_count
663
node_count
++;
664
}
665
}
666
667
668
// Next is the bottom hyp. node -
669
// lower tri. el. number 7, upper tri. el. number 8
670
element_count
= 0;
671
// Note: both i,j counters reduced by 1
672
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
673
{
674
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
675
{
676
// The node always needs constructing
677
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
678
->construct_node(7,
time_stepper_pt
);
679
680
// Set pointers from the elements to new nodes
681
finite_element_pt
(
element_count
)->node_pt(7) =
Node_pt
[
node_count
];
682
finite_element_pt
(
element_count
+ 1)->node_pt(8) =
683
Node_pt
[
node_count
];
684
685
// Increase element_count by 2
686
element_count
+= 2;
687
688
// Set position of current (Global) Node
689
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 2.0 / 3.0) *
xstep
;
690
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 1.0 / 3.0) *
ystep
;
691
692
// All nodes are internal in this structure so no boundaries
693
// applicable
694
695
// Update node_count
696
node_count
++;
697
}
698
}
699
700
701
// Next is the upper hyp. node -
702
// lower tri. el. number 8, upper tri. el. number 7
703
element_count
= 0;
704
// Note: both i,j counters reduced by 1
705
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
706
{
707
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
708
{
709
// The node always needs constructing
710
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
)
711
->construct_node(8,
time_stepper_pt
);
712
713
// Set pointers from the elements to new nodes
714
finite_element_pt
(
element_count
)->node_pt(8) =
Node_pt
[
node_count
];
715
finite_element_pt
(
element_count
+ 1)->node_pt(7) =
716
Node_pt
[
node_count
];
717
718
// Increase element_count by 2
719
element_count
+= 2;
720
721
// Set position of current (Global) Node
722
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 1.0 / 3.0) *
xstep
;
723
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 2.0 / 3.0) *
ystep
;
724
725
// All nodes are internal in this structure so no boundaries
726
// applicable
727
728
// Update node_count
729
node_count
++;
730
}
731
}
732
733
734
// Next is the upper tri. el. totally internal node with local number 9
735
element_count
= 0;
736
// Note: both i,j counters reduced by 1
737
for
(
unsigned
j
= 0;
j
<
Ny
; ++
j
)
738
{
739
for
(
unsigned
i
= 0;
i
<
Nx
; ++
i
)
740
{
741
// The middle node always needs constructing
742
Node_pt
[
node_count
] =
finite_element_pt
(
element_count
+ 1)
743
->construct_node(9,
time_stepper_pt
);
744
745
// Set pointers from the elements to new nodes
746
finite_element_pt
(
element_count
+ 1)->node_pt(9) =
747
Node_pt
[
node_count
];
748
749
// Increase element_count by 2
750
element_count
+= 2;
751
752
// Set position of current (Global) Node
753
Node_pt
[
node_count
]->x(0) =
xinit
+
double
(
i
+ 2.0 / 3.0) *
xstep
;
754
Node_pt
[
node_count
]->x(1) =
yinit
+
double
(
j
+ 2.0 / 3.0) *
ystep
;
755
756
// All nodes are internal in this structure so no boundaries
757
// applicable
758
759
// Update node_count
760
node_count
++;
761
}
762
}
763
}
764
}
765
766
}
// namespace oomph
767
#endif
oomph::SimpleRectangularTriMesh::SimpleRectangularTriMesh
SimpleRectangularTriMesh(const unsigned &n_x, const unsigned &n_y, const double &l_x, const double &l_y, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor n_x : number of elements in the x direction; n_y : number of elements in the y direction;...
Definition
simple_rectangular_tri_mesh.template.cc:54
oomph::SimpleRectangularTriMesh::Lx
double Lx
Length of mesh in x-direction.
Definition
simple_rectangular_tri_mesh.template.h:84
oomph::SimpleRectangularTriMesh::Ny
unsigned Ny
Number of elements in y directions.
Definition
simple_rectangular_tri_mesh.template.h:81
oomph::SimpleRectangularTriMesh::Nx
unsigned Nx
Number of elements in x direction.
Definition
simple_rectangular_tri_mesh.template.h:78
oomph::SimpleRectangularTriMesh::Ly
double Ly
Length of mesh in y-direction.
Definition
simple_rectangular_tri_mesh.template.h:87
oomph::TetgenMesh
Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen/.
Definition
tetgen_mesh.template.h:52
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition
annular_domain.h:35
simple_rectangular_tri_mesh.template.h