← Back to branch summary

~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« back to all changes in this revision

Viewing changes to nfem/geo_element/reference_element.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/geo_element/reference_element.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
/// 
 
20
/// =========================================================================
 
21
 
 
22
#include "reference_element.h"
 
23
#include "rheolef/point.h"
 
24
 
 
25
 
 
26
namespace rheolef {
 
27
 
 
28
namespace prism {
 
29
#include "prism.icc"
 
30
} // namespace prism
 
31
 
 
32
namespace hexahedron {
 
33
#include "hexahedron.icc"
 
34
} // namespace hexahedron
 
35
 
 
36
#include "reference_element_declare.cc"
 
37
 
 
38
#ifdef TO_CLEAN
 
39
// static data members:
 
40
const   reference_element::dof_family_type
 
41
        reference_element::Lagrange       = 0,
 
42
        reference_element::Hermite        = 1,
 
43
        reference_element::dof_family_max = 2;
 
44
#endif // TO_CLEAN
 
45
 
 
46
// ==============================================================================
 
47
// generic
 
48
// ==============================================================================
 
49
void
 
50
reference_element::set_name (char name)
 
51
{
 
52
  _x = variant(name);
 
53
  check_macro (_x != max_variant, "undefined reference element `" << name << "'");
 
54
}
 
55
Float
 
56
measure (reference_element hat_K) {
 
57
    return hat_K_measure [hat_K.variant()];
 
58
}
 
59
reference_element::variant_type
 
60
reference_element::variant (char name)
 
61
{
 
62
  for (size_type i = 0; i < max_variant; i++) {
 
63
    if (_name [i] == name) return variant_type(i);
 
64
  }
 
65
  return variant_type(max_variant);
 
66
}
 
67
reference_element::variant_type
 
68
reference_element::variant (size_type n_vertex, size_type dim)
 
69
{
 
70
  size_type variant = 0;
 
71
  for (; variant < max_variant; variant++) {
 
72
    if (_dimension[variant] == dim && _n_vertex[variant] == n_vertex) {
 
73
      return variant_type(variant);
 
74
    }
 
75
  }
 
76
  error_macro ("undefined "<<dim<<"d reference element with "<<n_vertex << " vertices");
 
77
}
 
78
reference_element::size_type
 
79
reference_element::n_sub_edge (variant_type variant)
 
80
{
 
81
  switch (variant) {
 
82
    case reference_element::p: return 0;
 
83
    case reference_element::e: return 0;
 
84
    case reference_element::t: return 3;
 
85
    case reference_element::q: return 4;
 
86
    case reference_element::T: return 6;
 
87
    case reference_element::P: return 9;
 
88
    case reference_element::H: return 12;
 
89
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
 
90
  }
 
91
}
 
92
reference_element::size_type
 
93
reference_element::n_sub_face (variant_type variant)
 
94
{
 
95
  switch (variant) {
 
96
    case reference_element::p: return 0;
 
97
    case reference_element::e: return 0;
 
98
    case reference_element::t: return 0;
 
99
    case reference_element::q: return 0;
 
100
    case reference_element::T: return 4;
 
101
    case reference_element::P: return 5;
 
102
    case reference_element::H: return 6;
 
103
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
 
104
  }
 
105
}
 
106
reference_element::size_type
 
107
reference_element::n_subgeo (variant_type variant, size_type subgeo_dim)
 
108
{
 
109
#define _RHEOLEF_reference_element_case(VARIANT)                         \
 
110
    case   reference_element::VARIANT:                                \
 
111
      return reference_element_##VARIANT::n_subgeo (subgeo_dim);
 
112
 
 
113
  switch (variant) {
 
114
        _RHEOLEF_reference_element_case(p)
 
115
        _RHEOLEF_reference_element_case(e)
 
116
        _RHEOLEF_reference_element_case(t)
 
117
        _RHEOLEF_reference_element_case(q)
 
118
        _RHEOLEF_reference_element_case(T)
 
119
        _RHEOLEF_reference_element_case(P)
 
120
        _RHEOLEF_reference_element_case(H)
 
121
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
 
122
  }
 
123
#undef  _RHEOLEF_reference_element_case
 
124
}
 
125
reference_element::size_type
 
126
reference_element::subgeo_n_node (
 
127
    variant_type variant,
 
128
    size_type    order, 
 
129
    size_type    subgeo_dim,
 
130
    size_type    loc_isid)
 
131
{
 
132
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
 
133
    case   reference_element::VARIANT:                                  \
 
134
      return reference_element_##VARIANT::subgeo_n_node (order, subgeo_dim, loc_isid);
 
135
 
 
136
  switch (variant) {
 
137
    _RHEOLEF_geo_element_auto_case(p)
 
138
    _RHEOLEF_geo_element_auto_case(e)
 
139
    _RHEOLEF_geo_element_auto_case(t)
 
140
    _RHEOLEF_geo_element_auto_case(q)
 
141
    _RHEOLEF_geo_element_auto_case(T)
 
142
    _RHEOLEF_geo_element_auto_case(P)
 
143
    _RHEOLEF_geo_element_auto_case(H)
 
144
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
 
145
  }
 
146
#undef _RHEOLEF_geo_element_auto_case
 
147
}
 
148
reference_element::size_type
 
149
reference_element::subgeo_local_node (
 
150
    variant_type variant,
 
151
    size_type    order, 
 
152
    size_type    subgeo_dim,
 
153
    size_type    loc_isid,
 
154
    size_type    loc_jsidnod)
 
155
{
 
156
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
 
157
    case   reference_element::VARIANT:                                  \
 
158
      return reference_element_##VARIANT::subgeo_local_node (order, subgeo_dim, loc_isid, loc_jsidnod);
 
159
 
 
160
  switch (variant) {
 
161
    _RHEOLEF_geo_element_auto_case(p)
 
162
    _RHEOLEF_geo_element_auto_case(e)
 
163
    _RHEOLEF_geo_element_auto_case(t)
 
164
    _RHEOLEF_geo_element_auto_case(q)
 
165
    _RHEOLEF_geo_element_auto_case(T)
 
166
    _RHEOLEF_geo_element_auto_case(P)
 
167
    _RHEOLEF_geo_element_auto_case(H)
 
168
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
 
169
  }
 
170
#undef _RHEOLEF_geo_element_auto_case
 
171
}
 
172
reference_element::size_type
 
173
reference_element::first_inod_by_variant (
 
174
    variant_type variant,
 
175
    size_type    order, 
 
176
    variant_type subgeo_variant)
 
177
{
 
178
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
 
179
    case   reference_element::VARIANT:                                  \
 
180
      return reference_element_##VARIANT::first_inod_by_variant (order, subgeo_variant);
 
181
 
 
182
  switch (variant) {
 
183
    _RHEOLEF_geo_element_auto_case(p)
 
184
    _RHEOLEF_geo_element_auto_case(e)
 
185
    _RHEOLEF_geo_element_auto_case(t)
 
186
    _RHEOLEF_geo_element_auto_case(q)
 
187
    _RHEOLEF_geo_element_auto_case(T)
 
188
    _RHEOLEF_geo_element_auto_case(P)
 
189
    _RHEOLEF_geo_element_auto_case(H)
 
190
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
 
191
  }
 
192
#undef _RHEOLEF_geo_element_auto_case
 
193
}
 
194
reference_element::size_type
 
195
reference_element::n_node (variant_type variant, size_type order)
 
196
{
 
197
  switch (variant) {
 
198
    case reference_element::p: return 1;
 
199
    case reference_element::e: return order+1;
 
200
    case reference_element::t: return (order+1)*(order+2)/2;
 
201
    case reference_element::q: return (order+1)*(order+1);
 
202
    case reference_element::T: return (order+1)*(order+2)*(order+3)/6;
 
203
    case reference_element::P: return (order+1)*(order+1)*(order+2)/2;
 
204
    case reference_element::H: return (order+1)*(order+1)*(order+1);
 
205
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
 
206
  }
 
207
}
 
208
void
 
209
reference_element::init_local_nnode_by_variant (
 
210
    size_type                                               order, 
 
211
    boost::array<size_type,reference_element::max_variant>& sz)
 
212
{
 
213
  sz [reference_element::p] = 1;
 
214
  sz [reference_element::e] = order-1;
 
215
  sz [reference_element::t] = (order-1)*(order-2)/2;
 
216
  sz [reference_element::q] = (order-1)*(order-1);
 
217
  sz [reference_element::T] = (order-1)*(order-2)*(order-3)/6;
 
218
  sz [reference_element::P] = (order-1)*(order-1)*(order-2)/2;
 
219
  sz [reference_element::H] = (order-1)*(order-1)*(order-1);
 
220
}
 
221
// ==============================================================================
 
222
// point
 
223
// ==============================================================================
 
224
reference_element_p::size_type
 
225
reference_element_p::n_subgeo (size_type side_dim)
 
226
{
 
227
  return (side_dim == 0) ? 1 : 0;
 
228
}
 
229
reference_element_p::size_type
 
230
reference_element_p::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
 
231
{
 
232
  return (side_dim == 0) ? 1 : 0;
 
233
}
 
234
reference_element_p::size_type
 
235
reference_element_p::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
 
236
{
 
237
  return 0;
 
238
}
 
239
// edge 0d-lattice for high order elements, i <= 0
 
240
// convert to local inod geo_element number for rheolef
 
241
reference_element_p::size_type
 
242
reference_element_p::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
243
{
 
244
  return 0;
 
245
}
 
246
reference_element_p::size_type
 
247
reference_element_p::first_inod_by_variant (
 
248
    size_type    order, 
 
249
    size_type    subgeo_variant)
 
250
{
 
251
  if (subgeo_variant == reference_element::p) return 0;
 
252
  return 1;
 
253
}
 
254
// ==============================================================================
 
255
// edge
 
256
// ==============================================================================
 
257
reference_element_e::size_type
 
258
reference_element_e::n_subgeo (size_type side_dim)
 
259
{
 
260
  switch (side_dim) {
 
261
    case 0:  return 2;
 
262
    case 1:  return 1;
 
263
    default: return 0;
 
264
  }
 
265
}
 
266
reference_element_e::size_type
 
267
reference_element_e::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
 
268
{
 
269
  switch (side_dim) {
 
270
    case 0:  return 1;
 
271
    case 1:  return order+1;
 
272
    default: return 0;
 
273
  }
 
274
}
 
275
reference_element_e::size_type
 
276
reference_element_e::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
 
277
{
 
278
  switch (side_dim) {
 
279
    case 0:  return loc_isid;
 
280
    case 1:  return loc_jsidnod;
 
281
    default: return 0;
 
282
  }
 
283
}
 
284
// edge 1d-lattice (i) for high order elements, i <= order
 
285
// convert to local inod geo_element number for rheolef
 
286
reference_element_e::size_type
 
287
reference_element_e::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
288
{
 
289
  size_type i = ilat[0];
 
290
  if (i == 0) return 0;
 
291
  if (i == order) return 1;
 
292
  return 2 + (i - 1);
 
293
}
 
294
reference_element_e::size_type
 
295
reference_element_e::first_inod_by_variant (
 
296
    size_type    order, 
 
297
    size_type    subgeo_variant)
 
298
{
 
299
  if (subgeo_variant == reference_element::p) return 0;
 
300
  if (subgeo_variant == reference_element::e) return 2;
 
301
  return order+1;
 
302
}
 
303
// ==============================================================================
 
304
// triangle
 
305
// ==============================================================================
 
306
reference_element_t::size_type
 
307
reference_element_t::n_subgeo (size_type side_dim)
 
308
{
 
309
  switch (side_dim) {
 
310
    case 0:  return 3;
 
311
    case 1:  return 3;
 
312
    case 2:  return 1;
 
313
    default: return 0;
 
314
  }
 
315
}
 
316
reference_element_t::size_type
 
317
reference_element_t::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
 
318
{
 
319
  switch (side_dim) {
 
320
    case 0:  return 1;
 
321
    case 1:  return order+1;
 
322
    case 2:  return ((order+1)*(order+2))/2;
 
323
    default: return 0;
 
324
  }
 
325
}
 
326
reference_element_t::size_type
 
327
reference_element_t::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
 
328
{
 
329
  switch (side_dim) {
 
330
    case 0:  return loc_isid;
 
331
    case 1:  if (loc_jsidnod < 2) return (loc_isid + loc_jsidnod) % 3;         // edge-node is a vertex
 
332
    else                          return (order-1)*loc_isid + loc_jsidnod + 1; // edge-node is edge-internal
 
333
    case 2:  return loc_jsidnod;
 
334
    default: return 0;
 
335
  }
 
336
}
 
337
// triangle lattice (i,j) for high order elements, i+j <= order
 
338
// convert to local inod geo_element number for rheolef
 
339
reference_element_t::size_type
 
340
reference_element_t::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
341
{
 
342
  size_type i = ilat[0];
 
343
  size_type j = ilat[1];
 
344
  if (j == 0) { // horizontal edge [0,1]
 
345
    if (i == 0) return 0;
 
346
    if (i == order) return 1;
 
347
    return 3 + 0*(order-1) + (i - 1);
 
348
  }
 
349
  if (i + j == order) { // oblique edge ]1,2]
 
350
    if (j == order) return 2;
 
351
    return 3 + 1*(order-1) + (j - 1);
 
352
  }
 
353
  size_type j1 = order - j;
 
354
  if (i == 0) { // vertical edge ]2,0[
 
355
    return 3 + 2*(order-1) + (j1 - 1);
 
356
  }
 
357
  // internal face ]0,1[x]0,2[: i, then j
 
358
  size_type n_face_node = (order-1)*(order-2)/2;
 
359
  size_type ij = (n_face_node - (j1-1)*j1/2) + (i - 1);
 
360
  return 3 + 3*(order-1) + ij;
 
361
}
 
362
reference_element_t::size_type
 
363
reference_element_t::first_inod_by_variant (
 
364
    size_type    order, 
 
365
    size_type    subgeo_variant)
 
366
{
 
367
  if (subgeo_variant == reference_element::p) return 0;
 
368
  if (subgeo_variant == reference_element::e) return 3;
 
369
  if (subgeo_variant == reference_element::t) return 3 + 3*(order-1);
 
370
  return (order+1)*(order+2)/2;
 
371
}
 
372
// ==============================================================================
 
373
// quadrangle
 
374
// ==============================================================================
 
375
reference_element_q::size_type
 
376
reference_element_q::n_subgeo (size_type side_dim)
 
377
{
 
378
  switch (side_dim) {
 
379
    case 0:  return 4;
 
380
    case 1:  return 4;
 
381
    case 2:  return 1;
 
382
    default: return 0;
 
383
  }
 
384
}
 
385
reference_element_q::size_type
 
386
reference_element_q::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
 
387
{
 
388
  switch (side_dim) {
 
389
    case 0:  return 1;
 
390
    case 1:  return order+1;
 
391
    case 2:  return (order+1)*(order+1);
 
392
    default: return 0;
 
393
  }
 
394
}
 
395
reference_element_q::size_type
 
396
reference_element_q::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
 
397
{
 
398
  switch (side_dim) {
 
399
    case 0:  return loc_isid;
 
400
    case 1:  if (loc_jsidnod < 2) return (loc_isid + loc_jsidnod) % 4;         // edge-node is a vertex
 
401
    else                          return (order-1)*loc_isid + loc_jsidnod + 2; // edge-node is edge-internal
 
402
    case 2:  return loc_jsidnod;
 
403
    default: return 0;
 
404
  }
 
405
}
 
406
// quadrangle lattice (i,j) for high order elements, 0 <= i,j <= order
 
407
// convert to local inod geo_element number for rheolef
 
408
reference_element_q::size_type
 
409
reference_element_q::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
410
{
 
411
  size_type i = ilat[0];
 
412
  size_type j = ilat[1];
 
413
  if (j == 0) { // horizontal edge [0,1]
 
414
    if (i == 0) return 0;
 
415
    if (i == order) return 1;
 
416
    return 4 + 0*(order-1) + (i - 1);
 
417
  }
 
418
  if (i == order) { // vertical edge ]1,2]
 
419
    if (j == order) return 2;
 
420
    return 4 + 1*(order-1) + (j - 1);
 
421
  }
 
422
  size_type i1 = order - i;
 
423
  if (j == order) { // horizontal edge ]2,3]
 
424
    if (i == 0) return 3;
 
425
    return 4 + 2*(order-1) + (i1 - 1);
 
426
  }
 
427
  size_type j1 = order - j;
 
428
  if (i == 0) { // vertical edge ]3,0[
 
429
    return 4 + 3*(order-1) + (j1 - 1);
 
430
  }
 
431
  // internal face ]0,1[x]0,3[: i, then j
 
432
  size_type ij = (order-1)*(j-1) + (i-1);
 
433
  size_type n_bdry_node = 4 + 4*(order-1);
 
434
  return n_bdry_node + ij;
 
435
}
 
436
reference_element_q::size_type
 
437
reference_element_q::first_inod_by_variant (
 
438
    size_type    order, 
 
439
    size_type    subgeo_variant)
 
440
{
 
441
  if (subgeo_variant == reference_element::p) return 0;
 
442
  if (subgeo_variant == reference_element::e) return 4;
 
443
  if (subgeo_variant == reference_element::t) return 4 + 4*(order-1);
 
444
  if (subgeo_variant == reference_element::q) return 4 + 4*(order-1);
 
445
  return (order+1)*(order+1);
 
446
}
 
447
// ==============================================================================
 
448
// tetra
 
449
// ==============================================================================
 
450
reference_element_T::size_type
 
451
reference_element_T::n_subgeo (size_type side_dim)
 
452
{
 
453
  switch (side_dim) {
 
454
    case 0:  return 4;
 
455
    case 1:  return 6;
 
456
    case 2:  return 4;
 
457
    case 3:  return 1;
 
458
    default: return 0;
 
459
  }
 
460
}
 
461
reference_element_T::size_type
 
462
reference_element_T::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid) 
 
463
{
 
464
  switch (side_dim) {
 
465
    case 0:  return 1;
 
466
    case 1:  return order+1;
 
467
    case 2:  return (order+1)*(order+2)/2;
 
468
    case 3:  return ((order+1)*(order+2)*(order+3))/6;
 
469
    default: return 0;
 
470
  }
 
471
}
 
472
reference_element_T::size_type
 
473
reference_element_T::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod) 
 
474
{
 
475
  switch (side_dim) {
 
476
    case 0:  return loc_isid;
 
477
    case 1:  if (loc_jsidnod < 2) { // edge-node is a vertex
 
478
               if (loc_isid < 3) return (loc_isid + loc_jsidnod) % 3;
 
479
               else              return loc_jsidnod == 0 ? loc_isid - 3 : 3;
 
480
             } else {               // edge-node is internal to the edge
 
481
               return loc_isid*(order-1) + loc_jsidnod + 2;
 
482
             }
 
483
    case 2:  {
 
484
               if (loc_jsidnod < 3) { // face-node is a vertex
 
485
                 if (loc_isid == 3)    return loc_jsidnod + 1;
 
486
                 if (loc_jsidnod == 0) return 0;
 
487
                 if (loc_jsidnod == 2) return loc_isid + 1;
 
488
                 return ((loc_isid + 1) % 3) + 1;
 
489
               }
 
490
               // face-node is edge-internal(P2) or face-internal(P3) ; volume-internal(P4) impossible
 
491
               size_type last_edge_node_iloc = 3 + 3*(order-1);
 
492
               if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
 
493
#ifdef TO_CLEAN
 
494
                 extern const size_type geo_element_T_fac2edg_idx    [4][3];
 
495
                 extern const int       geo_element_T_fac2edg_orient [4][3];
 
496
#endif // TO_CLEAN
 
497
                 size_type order1 = order - 1; // avoid div by zero compiler error
 
498
                 size_type loc_jedg = (loc_jsidnod-3) / order1;
 
499
                 size_type loc_kedg = (loc_jsidnod-3) % order1;
 
500
                 if (geo_element_T_fac2edg_orient [loc_isid][loc_jedg] < 0) {
 
501
                        loc_kedg = order - loc_kedg - 2;
 
502
                 }
 
503
                 return (order-1)*geo_element_T_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg + 4;
 
504
               }
 
505
               // face-node is face-internal(P3)
 
506
               size_type ij_loc = (loc_jsidnod - last_edge_node_iloc);
 
507
               return 4 + 6*(order-1) + loc_isid*(order-1)*(order-2)/2 + ij_loc;
 
508
             }
 
509
    case 3:  return loc_jsidnod;
 
510
    default: return 0;
 
511
  }
 
512
}
 
513
// tetrahedron lattice (i,j,k) for high order elements, i+j+k <= order
 
514
// convert to local inod geo_element number for rheolef
 
515
reference_element_T::size_type
 
516
reference_element_T::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
517
{
 
518
  size_type i = ilat[0];
 
519
  size_type j = ilat[1];
 
520
  size_type k = ilat[2];
 
521
  size_type n_tot_edge_node = 6*(order-1);
 
522
  size_type n_face_node = (order-1)*(order-2)/2;
 
523
  if (k == 0) {
 
524
    // horizontal face 021 = [0,2]x[0,1]
 
525
    if (j == 0) { // left edge [0,1]
 
526
      if (i == 0) return 0;
 
527
      if (i == order) return 1;
 
528
      return 4 + 0*(order-1) + (i - 1);
 
529
    }
 
530
    if (i + j == order) { // oblique edge ]1,2[
 
531
      if (j == order) return 2;
 
532
      return 4 + 1*(order-1) + (j - 1);
 
533
    }
 
534
    size_type j1 = order - j;
 
535
    if (i == 0) { // right edge ]0,2[
 
536
      return 4 + 2*(order-1) + (j1 - 1);
 
537
    }
 
538
    // internal face 021 = ]0,2[x]0,1[: j, then i
 
539
    size_type i1 = order - i;
 
540
    size_type ji = (n_face_node - i1*(i1-1)/2) + (j - 1);
 
541
    return 4 + n_tot_edge_node + 0*n_face_node + ji;
 
542
  }
 
543
  if (i == 0) {
 
544
    // back face 032 =]0,3]x]0,2[
 
545
    if (j == 0) { // vertical edge ]0,3]
 
546
      if (k == order) return 3;
 
547
      return 4 + 3*(order-1) + (k - 1);
 
548
    }
 
549
    if (j + k == order) { // oblique edge ]2,3[
 
550
      return 4 + 5*(order-1) + (k - 1);
 
551
    }
 
552
    // internal face 032 =]0,3]x]0,2[, k, then j
 
553
    size_type j1 = order - j;
 
554
    size_type kj = (n_face_node - j1*(j1-1)/2) + (k - 1);
 
555
    return 4 + n_tot_edge_node + 1*n_face_node + kj;
 
556
  }
 
557
  if (j == 0) {
 
558
    // left face 013 = ]0,1[x]0,3[
 
559
    if (i + k == order) { // oblique edge ]1,3[
 
560
      return 4 + 4*(order-1) + (k - 1);
 
561
    }
 
562
    // internal face 013 = ]0,1[x]0,3[: i, then k
 
563
    size_type k1 = order - k;
 
564
    size_type ik = (n_face_node - k1*(k1-1)/2) + (i - 1);
 
565
    return 4 + n_tot_edge_node + 2*n_face_node + ik;
 
566
  }
 
567
  if (i + j + k == order) {
 
568
    // internal to oblique face 123 = ]1,2[x]1,3], j, then k
 
569
    size_type k1 = order - k;
 
570
    size_type jk = (n_face_node - k1*(k1-1)/2) + (j - 1);
 
571
    return 4 + n_tot_edge_node + 3*n_face_node + jk;
 
572
  }
 
573
  // internal to volume: i, then j, then k
 
574
  size_type n_tot_face_node = 4*n_face_node;
 
575
  size_type n_tot_node = (order+1)*(order+2)*(order+3)/6;
 
576
  size_type n_volume_node = (order-1)*(order-2)*(order-3)/6;
 
577
  size_type k1 = order - k;
 
578
  size_type j1 = order - j - k;
 
579
  size_type n_level_k_node = (k1-1)*(k1-2)/2;
 
580
  size_type ijk = (n_volume_node  - k1*(k1-1)*(k1-2)/6) 
 
581
                + (n_level_k_node - j1*(j1-1)/2)
 
582
                + (i - 1);
 
583
  return 4 + n_tot_edge_node + n_tot_face_node + ijk;
 
584
}
 
585
reference_element_T::size_type
 
586
reference_element_T::first_inod_by_variant (
 
587
    size_type    order, 
 
588
    size_type    subgeo_variant)
 
589
{
 
590
  if (subgeo_variant == reference_element::p) return 0;
 
591
  if (subgeo_variant == reference_element::e) return 4;
 
592
  if (subgeo_variant == reference_element::t) return 4 + 6*(order-1);
 
593
  if (subgeo_variant == reference_element::q) return 4 + 6*(order-1) + 4*((order-1)*(order-2)/2);
 
594
  if (subgeo_variant == reference_element::T) return 4 + 6*(order-1) + 4*((order-1)*(order-2)/2);
 
595
  return (order+1)*(order+2)*(order+3)/6;
 
596
}
 
597
reference_element_T::size_type
 
598
reference_element_T::face2edge (size_type loc_iface, size_type loc_iface_jedg)
 
599
{
 
600
  return geo_element_T_fac2edg_idx [loc_iface][loc_iface_jedg];
 
601
}
 
602
int
 
603
reference_element_T::face2edge_orient (size_type loc_iface, size_type loc_iface_jedg)
 
604
{
 
605
  return geo_element_T_fac2edg_orient [loc_iface][loc_iface_jedg];
 
606
}
 
607
// ==============================================================================
 
608
// prism
 
609
// ==============================================================================
 
610
reference_element_P::size_type
 
611
reference_element_P::n_subgeo (size_type side_dim)
 
612
{
 
613
  switch (side_dim) {
 
614
    case 0:  return 6;
 
615
    case 1:  return 9;
 
616
    case 2:  return 5;
 
617
    case 3:  return 1;
 
618
    default: return 0;
 
619
  }
 
620
}
 
621
reference_element_P::size_type
 
622
reference_element_P::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
 
623
{
 
624
  switch (side_dim) {
 
625
    case 0:  return 1;
 
626
    case 1:  return order+1;
 
627
    case 2:  return (loc_isid < 2) ? (order+1)*(order+2)/2 : (order+1)*(order+1);
 
628
    case 3:  return ((order+1)*(order+1)*(order+2))/2;
 
629
    default: return 0;
 
630
  }
 
631
}
 
632
reference_element_P::size_type
 
633
reference_element_P::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
 
634
{
 
635
  switch (side_dim) {
 
636
    case 0:  return loc_isid;
 
637
    case 1:  if (loc_jsidnod < 2) { // edge-node is a vertex
 
638
               return prism::edge [loc_isid] [loc_jsidnod];
 
639
             } else {               // edge-node is internal to the edge
 
640
               return 6 + loc_isid*(order-1) + (loc_jsidnod - 2);
 
641
             }
 
642
    case 2:  {
 
643
               size_type n_vert_on_side = (loc_isid < 2) ? 3 : 4;
 
644
               if (loc_jsidnod < n_vert_on_side) { // face-node is a vertex
 
645
                return prism::face [loc_isid] [loc_jsidnod];
 
646
               }
 
647
               // face-node is edge-internal(P2) or face-internal(P3) ; volume-internal is impossible
 
648
               size_type last_edge_node_iloc = n_vert_on_side + n_vert_on_side*(order-1);
 
649
               if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
 
650
#ifdef TO_CLEAN
 
651
                 extern const size_type geo_element_P_fac2edg_idx    [5][4];
 
652
                 extern const int       geo_element_P_fac2edg_orient [5][4];
 
653
#endif // TO_CLEAN
 
654
                 size_type order1 = order - 1; // avoid div by zero compiler error
 
655
                 size_type loc_jedg = (loc_jsidnod - n_vert_on_side) / order1;
 
656
                 size_type loc_kedg = (loc_jsidnod - n_vert_on_side) % order1;
 
657
                 if (geo_element_P_fac2edg_orient [loc_isid][loc_jedg] < 0) {
 
658
                   loc_kedg = order - loc_kedg - 2;
 
659
                 }
 
660
                 return 6 + (order-1)*geo_element_P_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg;
 
661
               }
 
662
               // face-node is face-internal(P3)
 
663
               size_type ij_loc = (loc_jsidnod - last_edge_node_iloc);
 
664
               size_type shift_prev_faces;
 
665
               if (loc_isid < 2) {
 
666
                 shift_prev_faces = loc_isid*(order-1)*(order-2)/2;
 
667
               } else {
 
668
                 shift_prev_faces = 2*(order-1)*(order-2)/2 + (loc_isid-2)*(order-1)*(order-1);
 
669
               }
 
670
               return 6 + 9*(order-1) + shift_prev_faces + ij_loc;
 
671
             }
 
672
    case 3:  return loc_jsidnod;
 
673
    default: return 0;
 
674
  }
 
675
}
 
676
reference_element_P::size_type
 
677
reference_element_P::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
678
{
 
679
  size_type i = ilat[0];
 
680
  size_type j = ilat[1];
 
681
  size_type k = ilat[2];
 
682
  size_type i1 = order - i;
 
683
  size_type j1 = order - j;
 
684
  size_type k1 = order - k;
 
685
  size_type n_node_edge = order-1;
 
686
  size_type n_tot_edge_node = 9*(order-1);
 
687
  size_type n_node_face_tri = (order-1)*(order-2)/2;
 
688
  size_type n_node_face_qua = (order-1)*(order-1);
 
689
  if (k == 0) {
 
690
    // horizontal face 021 = [0,1]x[0,2]
 
691
    if (j == 0) { // left edge [0,1]
 
692
      if (i == 0) return 0;
 
693
      if (i == order) return 1;
 
694
      return 6 + 0*(order-1) + (i - 1);
 
695
    }
 
696
    if (i + j == order) { // oblique edge ]1,2[
 
697
      if (j == order) return 2;
 
698
      return 6 + 1*(order-1) + (j - 1);
 
699
    }
 
700
    if (i == 0) { // right edge ]0,2[
 
701
      return 6 + 2*(order-1) + (j1 - 1);
 
702
    }
 
703
    // internal face 021 = ]0,2[x]0,1[: j, then i
 
704
    size_type ji = (n_node_face_tri - i1*(i1-1)/2) + (j - 1);
 
705
    return 6 + n_tot_edge_node + 0*n_node_face_tri + ji;
 
706
  }
 
707
  if (k == order) {
 
708
    // horizontal face 354 = [3,4]x[3,5]
 
709
    if (j == 0) { // left edge [3,4]
 
710
      if (i == 0) return 3;
 
711
      if (i == order) return 4;
 
712
      return 6 + 6*(order-1) + (i - 1);
 
713
    }
 
714
    if (i + j == order) { // oblique edge ]4,5[
 
715
      if (j == order) return 5;
 
716
      return 6 + 7*(order-1) + (j - 1);
 
717
    }
 
718
    if (i == 0) { // right edge ]3,5[
 
719
      return 6 + 8*(order-1) + (j1 - 1);
 
720
    }
 
721
    // internal face 354 = ]3,4[x]3,5[: i, then j
 
722
    size_type ij = (n_node_face_tri - j1*(j1-1)/2) + (i - 1);
 
723
    return 6 + n_tot_edge_node + 1*n_node_face_tri + ij;
 
724
  }
 
725
  if (j == 0) {
 
726
    // left face 0143 = [0,1]x]0,3[
 
727
    if (i == 0) {
 
728
      // vertical edge ]0,3[
 
729
      return 6 + 3*(order-1) + (k - 1);
 
730
    }
 
731
    if (i == order) {
 
732
      // vertical edge ]1,4[
 
733
      return 6 + 4*(order-1) + (k - 1);
 
734
    }
 
735
    // internal face 01243 = ]0,1[x]0,3[: i then k
 
736
    size_type ik = n_node_edge*(k-1) + (i-1);
 
737
    return 6 + n_tot_edge_node + 2*n_node_face_tri + 0*n_node_face_qua + ik;
 
738
  }
 
739
  if (i == 0) {
 
740
    // back face 0352 = ]0,2]x]0,3[
 
741
    if (j == order) {
 
742
      // vertical edge ]2,5[
 
743
      return 6 + 5*(order-1) + (k - 1);
 
744
    }
 
745
    // internal face 0352 = ]0,2[x]0,3[: k then j
 
746
    size_type kj = n_node_edge*(j-1) + (k-1);
 
747
    return 6 + n_tot_edge_node + 2*n_node_face_tri + 2*n_node_face_qua + kj;
 
748
  }
 
749
  if (i + j == order) {
 
750
    // internal face 1254 = ]1,2[x]1,4[: j then k
 
751
    size_type jk = n_node_edge*(k-1) + (j-1);
 
752
    return 6 + n_tot_edge_node + 2*n_node_face_tri + 1*n_node_face_qua + jk;
 
753
  }
 
754
  // internal volume ]0,1[x]0,2[x][0,3[: i, then j, then k
 
755
  size_type n_tot_face_node = 2*n_node_face_tri + 3*n_node_face_qua;
 
756
  size_type ij = (n_node_face_tri - (j1-1)*j1/2) + (i - 1);
 
757
  size_type ijk = n_node_face_tri*(k-1) + ij;
 
758
  return 6 + n_tot_edge_node + n_tot_face_node + ijk;
 
759
}
 
760
reference_element_P::size_type
 
761
reference_element_P::first_inod_by_variant (
 
762
    size_type    order, 
 
763
    size_type    subgeo_variant)
 
764
{
 
765
  if (subgeo_variant == reference_element::p) return 0;
 
766
  if (subgeo_variant == reference_element::e) return 6;
 
767
  if (subgeo_variant == reference_element::t) return 6 + 9*(order-1);
 
768
  if (subgeo_variant == reference_element::q) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2);
 
769
  if (subgeo_variant == reference_element::T) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2) + 3*(order-1)*(order-1);
 
770
  if (subgeo_variant == reference_element::P) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2) + 3*(order-1)*(order-1);
 
771
  return ((order+1)*(order+1)*(order+2))/2;
 
772
}
 
773
// ==============================================================================
 
774
// hexa
 
775
// ==============================================================================
 
776
 
 
777
reference_element_H::size_type
 
778
reference_element_H::n_subgeo (size_type side_dim)
 
779
{
 
780
  switch (side_dim) {
 
781
    case 0:  return 8;
 
782
    case 1:  return 12;
 
783
    case 2:  return 6;
 
784
    case 3:  return 1;
 
785
    default: return 0;
 
786
  }
 
787
}
 
788
reference_element_H::size_type
 
789
reference_element_H::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
 
790
{
 
791
  switch (side_dim) {
 
792
    case 0:  return 1;
 
793
    case 1:  return order+1;
 
794
    case 2:  return (order+1)*(order+1);
 
795
    case 3:  return (order+1)*(order+1)*(order+1);
 
796
    default: return 0;
 
797
  }
 
798
}
 
799
reference_element_H::size_type
 
800
reference_element_H::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
 
801
{
 
802
  switch (side_dim) {
 
803
    case 0:  return loc_isid;
 
804
    case 1:  if (loc_jsidnod < 2) { // edge-node is a vertex
 
805
                    if (loc_isid < 4) return (loc_isid + loc_jsidnod) % 4;
 
806
               else if (loc_isid < 8) return loc_isid - 4 + loc_jsidnod*4;
 
807
               else                   return loc_jsidnod == 0 ? loc_isid - 4 : (loc_isid - 7)%4 + 4;
 
808
             } else {               // edge-node is internal to the edge
 
809
               return loc_isid*(order-1) + loc_jsidnod + 6; // TODO : BUG? +8 a la place de +6 ??
 
810
             }
 
811
    case 2:  {
 
812
               if (loc_jsidnod < 4) { // face-node is a vertex
 
813
                      return hexahedron::face [loc_isid][loc_jsidnod];
 
814
               } 
 
815
               // face-node is edge-internal(P2) or face-internal(P2) ; volume-internal(P2) impossible
 
816
               size_type last_edge_node_iloc = 4 + 4*(order-1);
 
817
               if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
 
818
#ifdef TO_CLEAN
 
819
                 extern const size_t geo_element_H_fac2edg_idx    [6][4];
 
820
                 extern const int    geo_element_H_fac2edg_orient [6][4];
 
821
#endif // TO_CLEAN
 
822
                 size_type order1 = order - 1; // avoid div by zero compiler error
 
823
                 size_type loc_jedg = (loc_jsidnod-4) / order1;
 
824
                 size_type loc_kedg = (loc_jsidnod-4) % order1;
 
825
                 if (geo_element_H_fac2edg_orient [loc_isid][loc_jedg] < 0) {
 
826
                   loc_kedg = order - loc_kedg - 2;
 
827
                 }
 
828
                 return (order-1)*geo_element_H_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg + 8;
 
829
               }
 
830
               // face-node is face-internal(P2)
 
831
               return 8 + 12*(order-1) + (order-1)*(order-1)*loc_isid + (loc_jsidnod - last_edge_node_iloc);
 
832
             }
 
833
    case 3:  return loc_jsidnod;
 
834
    default: return 0;
 
835
  }
 
836
}
 
837
// edge 0d-lattice for high order elements, i <= 0
 
838
// convert to local inod geo_element number for rheolef
 
839
typedef reference_element_H::size_type size_type;
 
840
static
 
841
reference_element_H::size_type
 
842
H_ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
843
{
 
844
  size_type i = ilat[0];
 
845
  size_type j = ilat[1];
 
846
  size_type k = ilat[2];
 
847
  size_type i1 = order - i;
 
848
  size_type j1 = order - j;
 
849
  size_type k1 = order - k;
 
850
  size_type n_node_edge = order-1;
 
851
  size_type n_node_face = (order-1)*(order-1);
 
852
  if (k == 0) { // bottom face [0,3]x[0,1]
 
853
    if (j == 0) { // bottom-left edge [0,1]
 
854
      if (i == 0) return 0;
 
855
      if (i == order) return 1;
 
856
      return 8 + 0*n_node_edge + (i - 1);
 
857
    }
 
858
    if (j == order) { // bottom-right edge [3,2]
 
859
      if (i == 0) return 3;
 
860
      if (i == order) return 2;
 
861
      return 8 + 2*n_node_edge + (i1 - 1);
 
862
    }
 
863
    if (i == 0) { // bottom-back edge ]3,0[
 
864
      return 8 + 3*n_node_edge + (j1 - 1);
 
865
    }
 
866
    if (i == order) { // bottom-front edge ]1,2[
 
867
      return 8 + 1*n_node_edge + (j - 1);
 
868
    }
 
869
    // internal face ]0,3[x]0,1[: j then i (TODO: not checked since Pk(H) only with k <= 2 yet
 
870
    size_type ji = n_node_edge*(i-1) + (j-1);
 
871
    return 8 + 12*n_node_edge + 0*n_node_face + ji;
 
872
  }
 
873
  if (k == order) { // top face [4,5]x[4,7]
 
874
    if (j == 0) { // top-left edge [4,5]
 
875
      if (i == 0) return 4;
 
876
      if (i == order) return 5;
 
877
      return 8 + 8*n_node_edge + (i - 1);
 
878
    }
 
879
    if (j == order) { // top-right edge [7,6]
 
880
      if (i == 0) return 7;
 
881
      if (i == order) return 6;
 
882
      return 8 + 10*n_node_edge + (i1 - 1);
 
883
    }
 
884
    if (i == 0) { // top-back edge ]7,4[
 
885
      return 8 + 11*n_node_edge + (j1 - 1);
 
886
    }
 
887
    if (i == order) { // top-front edge ]5,6[
 
888
      return 8 + 9*n_node_edge + (j - 1);
 
889
    }
 
890
    // internal face ]4,5[x]4,7[: i then j
 
891
    size_type ij = n_node_edge*(j-1) + (i-1);
 
892
    return 8 + 12*n_node_edge + 3*n_node_face + ij;
 
893
  }
 
894
  if (j == 0) { // left face ]0,1[x[0,4]
 
895
    if (i == 0) { // left-back edge [0,4]
 
896
      return 8 + 4*n_node_edge + (k - 1);
 
897
    }
 
898
    if (i == order) { // left-front edge [1,5]
 
899
      return 8 + 5*n_node_edge + (k - 1);
 
900
    }
 
901
    // internal face ]0,1[x]0,4[: i then k
 
902
    size_type ik = n_node_edge*(k-1) + (i-1);
 
903
    return 8 + 12*n_node_edge + 2*n_node_face + ik;
 
904
  }
 
905
  if (j == order) { // right face ]2,3[x[2,6]
 
906
    if (i == 0) { // right-back edge [3,7]
 
907
      return 8 + 7*n_node_edge + (k - 1);
 
908
    }
 
909
    if (i == order) { // right-front edge [2,6]
 
910
      return 8 + 6*n_node_edge + (k - 1);
 
911
    }
 
912
    // internal face ]2,3[x]2,6[: i1 then k
 
913
    size_type i1k = n_node_edge*(k-1) + (i1-1);
 
914
    return 8 + 12*n_node_edge + 5*n_node_face + i1k;
 
915
  }
 
916
  if (i == 0) { // internal back face ]0,4[x[0,3]: k then j
 
917
    size_type kj = n_node_edge*(j-1) + (k-1);
 
918
    return 8 + 12*n_node_edge + 1*n_node_face + kj;
 
919
  }
 
920
  if (i == order) { // internal front face ]1,2[x[1,5]: j then k
 
921
    size_type jk = n_node_edge*(k-1) + (j-1);
 
922
    return 8 + 12*n_node_edge + 4*n_node_face + jk;
 
923
  }
 
924
  // internal volume: i then j then k
 
925
  size_type n_node_bdry = 8 + 12*n_node_edge + 6*n_node_face;
 
926
  size_type ijk = (order-1)*((order-1)*(k-1) + (j-1)) + (i-1);
 
927
  return n_node_bdry + ijk;
 
928
}
 
929
reference_element_H::size_type
 
930
reference_element_H::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
 
931
{
 
932
  size_type loc_inod = H_ilat2loc_inod (order, ilat);
 
933
  return H_ilat2loc_inod (order, ilat);
 
934
}
 
935
reference_element_H::size_type
 
936
reference_element_H::first_inod_by_variant (
 
937
    size_type    order, 
 
938
    size_type    subgeo_variant)
 
939
{
 
940
  if (subgeo_variant == reference_element::p) return 0;
 
941
  if (subgeo_variant == reference_element::e) return 8;
 
942
  if (subgeo_variant == reference_element::t) return 8 + 12*(order-1);
 
943
  if (subgeo_variant == reference_element::q) return 8 + 12*(order-1);
 
944
  if (subgeo_variant == reference_element::T) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
 
945
  if (subgeo_variant == reference_element::P) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
 
946
  if (subgeo_variant == reference_element::H) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
 
947
  return (order+1)*(order+1)*(order+1);
 
948
}
 
949
 
 
950
} // namespace rheolef