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
About
People
Contact/Get involved
Publications
Acknowledgements
Copyright
Picture show
Go
src
generic
simple_cubic_scaffold_tet_mesh.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-2025 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
#include "
Telements.h
"
27
#include "
simple_cubic_scaffold_tet_mesh.h
"
28
29
30
namespace
oomph
31
{
32
//===========================================================
33
/// Scaffold mesh for cubic tet mesh.
34
//===========================================================
35
SimpleCubicScaffoldTetMesh::SimpleCubicScaffoldTetMesh
(
36
const
unsigned
&
n_x
,
37
const
unsigned
&
n_y
,
38
const
unsigned
&
n_z
,
39
const
double
&
l_x
,
40
const
double
& l_y,
41
const
double
&
l_z
,
42
TimeStepper
* time_stepper_pt)
43
{
44
// Storage for the pointers to the nodes
45
RankThreeTensor<Node*>
vertex_node_pt(
n_x
+ 1,
n_y
+ 1,
n_z
+ 1);
46
RankThreeTensor<Node*>
left_face_node_pt
(
n_x
+ 1,
n_y
+ 1,
n_z
+ 1);
47
RankThreeTensor<Node*>
front_face_node_pt
(
n_x
+ 1,
n_y
+ 1,
n_z
+ 1);
48
RankThreeTensor<Node*>
down_face_node_pt
(
n_x
+ 1,
n_y
+ 1,
n_z
+ 1);
49
RankThreeTensor<Node*>
central_node_pt
(
n_x
,
n_y
,
n_z
);
50
51
// "Lower left" corner
52
double
x_0
= 0.0;
53
double
y_0
= 0.0;
54
double
z_0
= 0.0;
55
56
// Increments
57
double
dx
=
l_x
/
double
(
n_x
);
58
double
dy
= l_y /
double
(
n_y
);
59
double
dz
=
l_z
/
double
(
n_z
);
60
61
// Total number of nodes
62
unsigned
nnod
= (
n_x
+ 1) * (
n_y
+ 1) * (
n_z
+ 1) +
// vertex nodes
63
n_x
*
n_y
*
n_z
+
// central nodes
64
(
n_x
*
n_y
) * (
n_z
+ 1) +
// faces
65
(
n_x
*
n_z
) * (
n_y
+ 1) + (
n_y
*
n_z
) * (
n_x
+ 1);
66
Node_pt
.reserve(
nnod
);
67
68
// Total number of elements
69
unsigned
nelem
= 24 *
n_x
*
n_y
*
n_z
;
70
Element_pt
.reserve(
nelem
);
71
72
// Six boundaries
73
set_nboundary
(6);
74
75
76
// Generate the vertex nodes of all cells
77
//=======================================
78
Node
*
node_pt
;
79
for
(
unsigned
k
= 0;
k
<
n_z
+ 1;
k
++)
80
{
81
for
(
unsigned
j
= 0;
j
<
n_y
+ 1;
j
++)
82
{
83
for
(
unsigned
i
= 0;
i
<
n_x
+ 1;
i
++)
84
{
85
// Need to work out whether to create a boundary node or not
86
if
(((
i
<
n_x
+ 1) && (
j
<
n_y
+ 1) && (
k
<
n_z
+ 1)) &&
87
((
i
== 0) || (
i
==
n_x
) || (
j
== 0) || (
j
==
n_y
) || (
k
== 0) ||
88
(
k
==
n_z
)))
89
{
90
node_pt
=
new
BoundaryNode<Node>
(3, 1, 0);
91
}
92
else
93
{
94
node_pt
=
new
Node
(3, 1, 0);
95
}
96
97
vertex_node_pt(
i
,
j
,
k
) =
node_pt
;
98
99
// Coordinates
100
node_pt
->
x
(0) =
x_0
+
double
(
i
) *
dx
;
101
node_pt
->
x
(1) =
y_0
+
double
(
j
) *
dy
;
102
node_pt
->
x
(2) =
z_0
+
double
(
k
) *
dz
;
103
104
// Real node?
105
if
((
i
<
n_x
+ 1) && (
j
<
n_y
+ 1) && (
k
<
n_z
+ 1))
106
{
107
// Add node
108
Node_pt
.push_back(
node_pt
);
109
110
// Left boundary is 0
111
if
(
i
== 0)
112
{
113
add_boundary_node
(0,
node_pt
);
114
}
115
// Right boundary is 1
116
else
if
(
i
==
n_x
)
117
{
118
add_boundary_node
(1,
node_pt
);
119
}
120
121
// Front boundary is 2
122
if
(
j
== 0)
123
{
124
add_boundary_node
(2,
node_pt
);
125
}
126
// Back boundary is 3
127
else
if
(
j
==
n_y
)
128
{
129
add_boundary_node
(3,
node_pt
);
130
}
131
132
// Down boundary is 4
133
if
(
k
== 0)
134
{
135
add_boundary_node
(4,
node_pt
);
136
}
137
// Up boundary is 5
138
else
if
(
k
==
n_z
)
139
{
140
add_boundary_node
(5,
node_pt
);
141
}
142
}
143
}
144
}
145
}
146
147
148
// Generate the face nodes of all cells
149
//=====================================
150
for
(
unsigned
k
= 0;
k
<
n_z
+ 1;
k
++)
151
{
152
for
(
unsigned
j
= 0;
j
<
n_y
+ 1;
j
++)
153
{
154
for
(
unsigned
i
= 0;
i
<
n_x
+ 1;
i
++)
155
{
156
// Node on front face of cell
157
//---------------------------
158
// Need to work out whether to create boundary node or not
159
if
(((
i
<
n_x
) && (
k
<
n_z
)) && ((
j
== 0) || (
j
==
n_y
)))
160
{
161
node_pt
=
new
BoundaryNode<Node>
(3, 1, 0);
162
}
163
else
164
{
165
node_pt
=
new
Node
(3, 1, 0);
166
}
167
front_face_node_pt
(
i
,
j
,
k
) =
node_pt
;
168
169
170
// Coordinates
171
node_pt
->
x
(0) =
x_0
+ 0.5 *
dx
+
double
(
i
) *
dx
;
172
node_pt
->
x
(1) =
y_0
+
double
(
j
) *
dy
;
173
node_pt
->
x
(2) =
z_0
+ 0.5 *
dz
+
double
(
k
) *
dz
;
174
175
// Real node?
176
if
((
i
<
n_x
) && (
k
<
n_z
))
177
{
178
// Add node to mesh
179
Node_pt
.push_back(
node_pt
);
180
181
// Front boundary is 2
182
if
(
j
== 0)
183
{
184
add_boundary_node
(2,
node_pt
);
185
}
186
// Back boundary is 3
187
else
if
(
j
==
n_y
)
188
{
189
add_boundary_node
(3,
node_pt
);
190
}
191
}
192
193
// Node on left face of cell
194
//--------------------------
195
// Work out whether to construct boundary node or not
196
if
(((
j
<
n_y
) && (
k
<
n_z
)) && ((
i
== 0) || (
i
==
n_x
)))
197
{
198
node_pt
=
new
BoundaryNode<Node>
(3, 1, 0);
199
}
200
else
201
{
202
node_pt
=
new
Node
(3, 1, 0);
203
}
204
left_face_node_pt
(
i
,
j
,
k
) =
node_pt
;
205
206
// Coordinates
207
node_pt
->
x
(0) =
x_0
+
double
(
i
) *
dx
;
208
node_pt
->
x
(1) =
y_0
+ 0.5 *
dy
+
double
(
j
) *
dy
;
209
node_pt
->
x
(2) =
z_0
+ 0.5 *
dz
+
double
(
k
) *
dz
;
210
211
212
// Real node?
213
if
((
j
<
n_y
) && (
k
<
n_z
))
214
{
215
// Add node to mesh
216
Node_pt
.push_back(
node_pt
);
217
218
// Left boundary is 0
219
if
(
i
== 0)
220
{
221
add_boundary_node
(0,
node_pt
);
222
}
223
// Right boundary is 1
224
else
if
(
i
==
n_x
)
225
{
226
add_boundary_node
(1,
node_pt
);
227
}
228
}
229
230
// Node on down face of cell
231
//--------------------------
232
if
(((
i
<
n_x
) && (
j
<
n_y
)) && ((
k
== 0) || (
k
==
n_z
)))
233
{
234
node_pt
=
new
BoundaryNode<Node>
(3, 1, 0);
235
}
236
else
237
{
238
node_pt
=
new
Node
(3, 1, 0);
239
}
240
down_face_node_pt
(
i
,
j
,
k
) =
node_pt
;
241
242
// Coordinates
243
node_pt
->
x
(0) =
x_0
+ 0.5 *
dx
+
double
(
i
) *
dx
;
244
node_pt
->
x
(1) =
y_0
+ 0.5 *
dy
+
double
(
j
) *
dy
;
245
node_pt
->
x
(2) =
z_0
+
double
(
k
) *
dz
;
246
247
248
// Real node?
249
if
((
i
<
n_x
) && (
j
<
n_y
))
250
{
251
// Add node to mesh
252
Node_pt
.push_back(
node_pt
);
253
254
// Down boundary is 4
255
if
(
k
== 0)
256
{
257
add_boundary_node
(4,
node_pt
);
258
}
259
// Up boundary is 5
260
else
if
(
k
==
n_z
)
261
{
262
add_boundary_node
(5,
node_pt
);
263
}
264
}
265
}
266
}
267
}
268
269
270
// Central nodes for all cells
271
for
(
unsigned
k
= 0;
k
<
n_z
;
k
++)
272
{
273
for
(
unsigned
j
= 0;
j
<
n_y
;
j
++)
274
{
275
for
(
unsigned
i
= 0;
i
<
n_x
;
i
++)
276
{
277
node_pt
=
new
Node
(3, 1, 0);
278
central_node_pt
(
i
,
j
,
k
) =
node_pt
;
279
Node_pt
.push_back(
node_pt
);
280
node_pt
->
x
(0) =
x_0
+ 0.5 *
dx
+
double
(
i
) *
dx
;
281
node_pt
->
x
(1) =
y_0
+ 0.5 *
dy
+
double
(
j
) *
dy
;
282
node_pt
->
x
(2) =
z_0
+ 0.5 *
dz
+
double
(
k
) *
dz
;
283
}
284
}
285
}
286
287
288
// Loop over blocks and create elements
289
TElement<3, 2>
*
el_pt
= 0;
290
for
(
unsigned
k
= 0;
k
<
n_z
;
k
++)
291
{
292
for
(
unsigned
j
= 0;
j
<
n_y
;
j
++)
293
{
294
for
(
unsigned
i
= 0;
i
<
n_x
;
i
++)
295
{
296
// FRONT FACE
297
//===========
298
299
// Left element on front face
300
//---------------------------
301
el_pt
=
new
TElement<3, 2>
;
302
303
// LFD
304
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
,
k
);
305
306
// Central front face
307
el_pt
->
node_pt
(1) =
front_face_node_pt
(
i
,
j
,
k
);
308
309
// LFU
310
el_pt
->
node_pt
(2) = vertex_node_pt(
i
,
j
,
k
+ 1);
311
312
// central node
313
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
314
315
Element_pt
.push_back(
el_pt
);
316
317
318
// Down element on front face
319
//---------------------------
320
el_pt
=
new
TElement<3, 2>
;
321
322
// LFD
323
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
,
k
);
324
325
// RFD
326
el_pt
->
node_pt
(1) = vertex_node_pt(
i
+ 1,
j
,
k
);
327
328
// Central front face
329
el_pt
->
node_pt
(2) =
front_face_node_pt
(
i
,
j
,
k
);
330
331
// Central node
332
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
333
334
Element_pt
.push_back(
el_pt
);
335
336
337
// Right element on front face
338
//---------------------------
339
el_pt
=
new
TElement<3, 2>
;
340
341
// RFD
342
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
,
k
);
343
344
// RFU
345
el_pt
->
node_pt
(1) = vertex_node_pt(
i
+ 1,
j
,
k
+ 1);
346
347
// Central front face
348
el_pt
->
node_pt
(2) =
front_face_node_pt
(
i
,
j
,
k
);
349
350
// Central node
351
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
352
353
Element_pt
.push_back(
el_pt
);
354
355
356
// Up element on front face
357
//---------------------------
358
el_pt
=
new
TElement<3, 2>
;
359
360
// RFU
361
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
,
k
+ 1);
362
363
// LFU
364
el_pt
->
node_pt
(1) = vertex_node_pt(
i
,
j
,
k
+ 1);
365
366
// Central front face
367
el_pt
->
node_pt
(2) =
front_face_node_pt
(
i
,
j
,
k
);
368
369
// Central node
370
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
371
372
Element_pt
.push_back(
el_pt
);
373
374
375
// RIGHT FACE
376
//===========
377
378
// Front element on right face
379
//---------------------------
380
el_pt
=
new
TElement<3, 2>
;
381
382
// RFD
383
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
,
k
);
384
385
// Central right face
386
el_pt
->
node_pt
(1) =
left_face_node_pt
(
i
+ 1,
j
,
k
);
387
388
// RFU
389
el_pt
->
node_pt
(2) = vertex_node_pt(
i
+ 1,
j
,
k
+ 1);
390
391
// Central node
392
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
393
394
Element_pt
.push_back(
el_pt
);
395
396
397
// Down element on right face
398
//---------------------------
399
el_pt
=
new
TElement<3, 2>
;
400
401
// RFD
402
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
,
k
);
403
404
// RBD
405
el_pt
->
node_pt
(1) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
);
406
407
// Central right face
408
el_pt
->
node_pt
(2) =
left_face_node_pt
(
i
+ 1,
j
,
k
);
409
410
// central node
411
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
412
413
Element_pt
.push_back(
el_pt
);
414
415
416
// Back element on right face
417
//---------------------------
418
el_pt
=
new
TElement<3, 2>
;
419
420
// RBD
421
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
);
422
423
// RBU
424
el_pt
->
node_pt
(1) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
+ 1);
425
426
// Central right face
427
el_pt
->
node_pt
(2) =
left_face_node_pt
(
i
+ 1,
j
,
k
);
428
429
// central node
430
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
431
432
Element_pt
.push_back(
el_pt
);
433
434
435
// Up element on right face
436
//---------------------------
437
el_pt
=
new
TElement<3, 2>
;
438
439
// RBU
440
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
+ 1);
441
442
// RFU
443
el_pt
->
node_pt
(1) = vertex_node_pt(
i
+ 1,
j
,
k
+ 1);
444
445
// Central right face
446
el_pt
->
node_pt
(2) =
left_face_node_pt
(
i
+ 1,
j
,
k
);
447
448
// Central node
449
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
450
451
Element_pt
.push_back(
el_pt
);
452
453
454
// UP FACE
455
//===========
456
457
// Front element on up face
458
//---------------------------
459
el_pt
=
new
TElement<3, 2>
;
460
461
// LFU
462
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
,
k
+ 1);
463
464
// RFU
465
el_pt
->
node_pt
(1) = vertex_node_pt(
i
+ 1,
j
,
k
+ 1);
466
467
// Central up face
468
el_pt
->
node_pt
(2) =
down_face_node_pt
(
i
,
j
,
k
+ 1);
469
470
// Central node
471
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
472
473
Element_pt
.push_back(
el_pt
);
474
475
476
// Front element on up face
477
//---------------------------
478
el_pt
=
new
TElement<3, 2>
;
479
480
// RFU
481
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
,
k
+ 1);
482
483
// RBU
484
el_pt
->
node_pt
(1) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
+ 1);
485
486
// Central up face
487
el_pt
->
node_pt
(2) =
down_face_node_pt
(
i
,
j
,
k
+ 1);
488
489
// central node
490
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
491
492
Element_pt
.push_back(
el_pt
);
493
494
495
// Back element on up face
496
//---------------------------
497
el_pt
=
new
TElement<3, 2>
;
498
499
// RBU
500
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
+ 1);
501
502
// LBU
503
el_pt
->
node_pt
(1) = vertex_node_pt(
i
,
j
+ 1,
k
+ 1);
504
505
// Central up face
506
el_pt
->
node_pt
(2) =
down_face_node_pt
(
i
,
j
,
k
+ 1);
507
508
// central node
509
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
510
511
Element_pt
.push_back(
el_pt
);
512
513
514
// Left element on up face
515
//---------------------------
516
el_pt
=
new
TElement<3, 2>
;
517
518
519
// LBU
520
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
+ 1,
k
+ 1);
521
522
// RFU
523
el_pt
->
node_pt
(1) = vertex_node_pt(
i
,
j
,
k
+ 1);
524
525
// Central up face
526
el_pt
->
node_pt
(2) =
down_face_node_pt
(
i
,
j
,
k
+ 1);
527
528
// Central node
529
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
530
531
Element_pt
.push_back(
el_pt
);
532
533
534
// DOWN FACE
535
//===========
536
537
// Front element on down face
538
//---------------------------
539
el_pt
=
new
TElement<3, 2>
;
540
541
// LFD
542
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
,
k
);
543
544
// RFD
545
el_pt
->
node_pt
(2) = vertex_node_pt(
i
+ 1,
j
,
k
);
546
547
// Central down face
548
el_pt
->
node_pt
(1) =
down_face_node_pt
(
i
,
j
,
k
);
549
550
// Central node
551
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
552
553
Element_pt
.push_back(
el_pt
);
554
555
556
// Front element on down face
557
//---------------------------
558
el_pt
=
new
TElement<3, 2>
;
559
560
// RFD
561
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
,
k
);
562
563
// RBD
564
el_pt
->
node_pt
(2) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
);
565
566
// Central down face
567
el_pt
->
node_pt
(1) =
down_face_node_pt
(
i
,
j
,
k
);
568
569
// central node
570
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
571
572
Element_pt
.push_back(
el_pt
);
573
574
575
// Back element on down face
576
//---------------------------
577
el_pt
=
new
TElement<3, 2>
;
578
579
// RBD
580
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
);
581
582
// LBD
583
el_pt
->
node_pt
(2) = vertex_node_pt(
i
,
j
+ 1,
k
);
584
585
// Central down face
586
el_pt
->
node_pt
(1) =
down_face_node_pt
(
i
,
j
,
k
);
587
588
// central node
589
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
590
591
Element_pt
.push_back(
el_pt
);
592
593
594
// Left element on down face
595
//---------------------------
596
el_pt
=
new
TElement<3, 2>
;
597
598
599
// LBD
600
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
+ 1,
k
);
601
602
// RFD
603
el_pt
->
node_pt
(2) = vertex_node_pt(
i
,
j
,
k
);
604
605
// Central down face
606
el_pt
->
node_pt
(1) =
down_face_node_pt
(
i
,
j
,
k
);
607
608
// Central node
609
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
610
611
Element_pt
.push_back(
el_pt
);
612
613
614
// BACK FACE
615
//===========
616
617
// Left element on back face
618
//---------------------------
619
el_pt
=
new
TElement<3, 2>
;
620
621
// LBD
622
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
+ 1,
k
);
623
624
// Central back face
625
el_pt
->
node_pt
(2) =
front_face_node_pt
(
i
,
j
+ 1,
k
);
626
627
// LBU
628
el_pt
->
node_pt
(1) = vertex_node_pt(
i
,
j
+ 1,
k
+ 1);
629
630
// central node
631
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
632
633
Element_pt
.push_back(
el_pt
);
634
635
636
// Down element on back face
637
//---------------------------
638
el_pt
=
new
TElement<3, 2>
;
639
640
// LBD
641
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
+ 1,
k
);
642
643
// RBD
644
el_pt
->
node_pt
(2) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
);
645
646
// Central back face
647
el_pt
->
node_pt
(1) =
front_face_node_pt
(
i
,
j
+ 1,
k
);
648
649
// Central node
650
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
651
652
Element_pt
.push_back(
el_pt
);
653
654
655
// Right element on back face
656
//---------------------------
657
el_pt
=
new
TElement<3, 2>
;
658
659
// RBD
660
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
);
661
662
// RBU
663
el_pt
->
node_pt
(2) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
+ 1);
664
665
// Central back face
666
el_pt
->
node_pt
(1) =
front_face_node_pt
(
i
,
j
+ 1,
k
);
667
668
// Central node
669
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
670
671
Element_pt
.push_back(
el_pt
);
672
673
674
// Up element on back face
675
//---------------------------
676
el_pt
=
new
TElement<3, 2>
;
677
678
// RBU
679
el_pt
->
node_pt
(0) = vertex_node_pt(
i
+ 1,
j
+ 1,
k
+ 1);
680
681
// LBU
682
el_pt
->
node_pt
(2) = vertex_node_pt(
i
,
j
+ 1,
k
+ 1);
683
684
// Central back face
685
el_pt
->
node_pt
(1) =
front_face_node_pt
(
i
,
j
+ 1,
k
);
686
687
// Central node
688
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
689
690
Element_pt
.push_back(
el_pt
);
691
692
693
// LEFT FACE
694
//===========
695
696
// Front element on left face
697
//---------------------------
698
el_pt
=
new
TElement<3, 2>
;
699
700
// LFD
701
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
,
k
);
702
703
// Central left face
704
el_pt
->
node_pt
(2) =
left_face_node_pt
(
i
,
j
,
k
);
705
706
// LFU
707
el_pt
->
node_pt
(1) = vertex_node_pt(
i
,
j
,
k
+ 1);
708
709
// Central node
710
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
711
712
Element_pt
.push_back(
el_pt
);
713
714
715
// Down element on left face
716
//---------------------------
717
el_pt
=
new
TElement<3, 2>
;
718
719
// LFD
720
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
,
k
);
721
722
// LBD
723
el_pt
->
node_pt
(2) = vertex_node_pt(
i
,
j
+ 1,
k
);
724
725
// Central left face
726
el_pt
->
node_pt
(1) =
left_face_node_pt
(
i
,
j
,
k
);
727
728
// central node
729
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
730
731
Element_pt
.push_back(
el_pt
);
732
733
734
// Back element on left face
735
//---------------------------
736
el_pt
=
new
TElement<3, 2>
;
737
738
// LBD
739
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
+ 1,
k
);
740
741
// LBU
742
el_pt
->
node_pt
(2) = vertex_node_pt(
i
,
j
+ 1,
k
+ 1);
743
744
// Central left face
745
el_pt
->
node_pt
(1) =
left_face_node_pt
(
i
,
j
,
k
);
746
747
// central node
748
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
749
750
Element_pt
.push_back(
el_pt
);
751
752
753
// Up element on left face
754
//---------------------------
755
el_pt
=
new
TElement<3, 2>
;
756
757
// LBU
758
el_pt
->
node_pt
(0) = vertex_node_pt(
i
,
j
+ 1,
k
+ 1);
759
760
// LFU
761
el_pt
->
node_pt
(2) = vertex_node_pt(
i
,
j
,
k
+ 1);
762
763
// Central left face
764
el_pt
->
node_pt
(1) =
left_face_node_pt
(
i
,
j
,
k
);
765
766
// Central node
767
el_pt
->
node_pt
(3) =
central_node_pt
(
i
,
j
,
k
);
768
769
Element_pt
.push_back(
el_pt
);
770
}
771
}
772
}
773
774
if
(
Node_pt
.
size
() !=
nnod
)
775
{
776
std::ostringstream
error_stream
;
777
error_stream
<<
"Some internal error in the constructor\n"
778
<<
"Actual number of nodes : "
<<
Node_pt
.
size
()
779
<<
"\ndoesn't match the prediction : "
<<
nnod
780
<< std::endl;
781
782
throw
OomphLibError
(
783
error_stream
.str(),
OOMPH_CURRENT_FUNCTION
,
OOMPH_EXCEPTION_LOCATION
);
784
}
785
786
787
if
(
Element_pt
.
size
() !=
nelem
)
788
{
789
std::ostringstream
error_stream
;
790
error_stream
<<
"Some internal error in the constructor\n"
791
<<
"Actual number of elements : "
<<
Element_pt
.
size
()
792
<<
"\ndoesn't match the prediction : "
<<
nelem
793
<< std::endl;
794
795
throw
OomphLibError
(
796
error_stream
.str(),
OOMPH_CURRENT_FUNCTION
,
OOMPH_EXCEPTION_LOCATION
);
797
}
798
799
// // Deliberately switch nodes to invert elements -- test for self_test()
800
// nelem=nelement();
801
// for (unsigned e=0;e<nelem;e++)
802
// {
803
// Node* node_pt=finite_element_pt(e)->node_pt(0);
804
// finite_element_pt(e)->node_pt(0)=
805
// finite_element_pt(e)->node_pt(1);
806
// finite_element_pt(e)->node_pt(1)=node_pt;
807
// }
808
}
809
810
811
}
// namespace oomph
Telements.h
i
cstr elem_len * i
Definition
cfortran.h:603
oomph::FiniteElement::size
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition
elements.cc:4320
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition
elements.h:2179
oomph::Mesh::add_boundary_node
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition
mesh.cc:243
oomph::Mesh::Node_pt
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition
mesh.h:183
oomph::Mesh::set_nboundary
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition
mesh.h:509
oomph::Mesh::node_pt
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition
mesh.h:440
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition
mesh.h:186
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition
nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition
nodes.h:1060
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition
oomph_definitions.h:229
oomph::SimpleCubicScaffoldTetMesh::SimpleCubicScaffoldTetMesh
SimpleCubicScaffoldTetMesh(const unsigned &n_x, const unsigned &n_y, const unsigned &n_z, const double &l_x, const double &l_y, const double &l_z, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements and dimensions of cube.
Definition
simple_cubic_scaffold_tet_mesh.cc:35
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition
Tadvection_diffusion_reaction_elements.h:66
oomph::TimeStepper
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition
timesteppers.h:231
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition
advection_diffusion_elements.cc:30
simple_cubic_scaffold_tet_mesh.h