Source code for sfepy.discrete.fem.facets

"""
Helper functions related to mesh facets and Lagrange FE approximation.

Line: ori - iter:

0 - iter0
1 - iter1

Triangle: ori - iter:

0 - iter21
1 - iter12
3 - iter02
4 - iter20
6 - iter10
7 - iter01

Possible couples:

1, 4, 7 <-> 0, 3, 6

Square: ori - iter:

 0 - iter10x01y
 7 - iter10y01x
11 - iter01y01x
30 - iter01x10y
33 - iter10x10y
52 - iter01y10x
56 - iter10y10x
63 - iter01x01y

Possible couples:

7, 33, 52, 63 <-> 0, 11, 30, 56

_quad_ori_groups:

i < j < k < l

all faces are permuted to

l --- k
|     |
|     |
i --- j

ijkl

which is the same as

l --- j
|     |
|     |
i --- k

ikjl

k --- l
|     |
|     |
i --- j

ijlk

- start at one vertex and go around clock-wise or anticlock-wise
-> 8 groups of 3
-> same face nodes order in
ijkl (63), ikjl (59), ijlk (31)
ilkj (11), iklj (15), iljk (43)
jkli ( 7), jlki ( 3), kjli ( 6)
kjil (56), jkil (57), ljik (48)
lijk (52), likj (20), kijl (60)
lkji ( 0), ljki ( 4), klji ( 1)
klij (33), lkij (32), jlik (41)
jilk (30), kilj (22), jikl (62)
"""
from __future__ import print_function
from __future__ import absolute_import
import numpy as nm
import six
from six.moves import range

_quad_ori_groups = {
    0 : 0,
    1 : 0,
    3 : 7,
    4 : 0,
    6 : 7,
    7 : 7,
    11 : 11,
    15 : 11,
    20 : 52,
    22 : 30,
    30 : 30,
    31 : 63,
    32 : 33,
    33 : 33,
    41 : 33,
    43 : 11,
    48 : 56,
    52 : 52,
    56 : 56,
    57 : 56,
    59 : 63,
    60 : 52,
    62 : 30,
    63 : 63,
}

[docs] def build_orientation_map(n_fp): """ The keys are binary masks of the lexicographical ordering of facet vertices. A bit i set to one means `v[i] < v[i+1]`. The values are `[original_order, permutation]`, where `permutation` can be used to sort facet vertices lexicographically. Hence `permuted_facet = facet[permutation]`. """ from sfepy.linalg import permutations indices = list(range(n_fp)) cmps = [(i1, i2) for i2 in indices for i1 in indices[:i2]] powers = [2**ii for ii in range(len(cmps))] ori_map = {} for indx in permutations(indices): key = 0 sign = 1 for ip, power in enumerate(powers): i1, i2 = cmps[ip] less = (indx[i1] < indx[i2]) key += power * less if not less: sign *= -1 isort = nm.argsort(indx) ori_map[key] = [indx, isort] return ori_map, cmps, powers
[docs] def iter0(num): for ir in range(num - 1, -1, -1): yield ir
[docs] def iter1(num): for ir in range(num): yield ir
ori_line_to_iter = { 0 : iter0, 1 : iter1, }
[docs] def make_line_matrix(order): if (order < 2): return nm.zeros((0, 0), dtype=nm.int32) oo = order - 1 mtx = nm.arange(oo, dtype=nm.int32) return mtx
[docs] def iter01(num): for ir in range(num - 1, -1, -1): for ic in range(ir + 1): yield ir, ic
[docs] def iter10(num): for ir in range(num - 1, -1, -1): for ic in range(ir, -1, -1): yield ir, ic
[docs] def iter02(num): for ic in range(num): for ir in range(num - 1, ic - 1, -1): yield ir, ic
[docs] def iter20(num): for ic in range(num): for ir in range(ic, num): yield ir, ic
[docs] def iter12(num): for idiag in range(num): irs, ics = nm.diag_indices(num - idiag) for ii in range(irs.shape[0] - 1, -1, -1): yield irs[ii] + idiag, ics[ii]
[docs] def iter21(num): for idiag in range(num): irs, ics = nm.diag_indices(num - idiag) for ii in range(irs.shape[0]): yield irs[ii] + idiag, ics[ii]
ori_triangle_to_iter = { 0 : iter21, 1 : iter12, 3 : iter02, 4 : iter20, 6 : iter10, 7 : iter01, }
[docs] def make_triangle_matrix(order): if (order < 3): return nm.zeros((0, 0), dtype=nm.int32) oo = order - 2 mtx = nm.zeros((oo, oo), dtype=nm.int32) for ii, (ir, ic) in enumerate(iter01(oo)): mtx[ir, ic] = ii return mtx
[docs] def iter01x01y(num): for ir in range(num): for ic in range(num): yield ir, ic
[docs] def iter01y01x(num): for ir, ic in iter01x01y(num): yield ic, ir
[docs] def iter10x01y(num): for ir in range(num - 1, -1, -1): for ic in range(num): yield ir, ic
[docs] def iter10y01x(num): for ir, ic in iter10x01y(num): yield ic, ir
[docs] def iter01x10y(num): for ir in range(num): for ic in range(num - 1, -1, -1): yield ir, ic
[docs] def iter01y10x(num): for ir, ic in iter01x10y(num): yield ic, ir
[docs] def iter10x10y(num): for ir in range(num - 1, -1, -1): for ic in range(num - 1, -1, -1): yield ir, ic
[docs] def iter10y10x(num): for ir, ic in iter10x10y(num): yield ic, ir
ori_square_to_iter = { 0 : iter10x01y, 7 : iter10y01x, 11 : iter01y01x, 30 : iter01x10y, 33 : iter10x10y, 52 : iter01y10x, 56 : iter10y10x, 63 : iter01x01y, }
[docs] def make_square_matrix(order): if (order < 2): return nm.zeros((0, 0), dtype=nm.int32) oo = order - 1 mtx = nm.arange(oo * oo, dtype=nm.int32) mtx.shape = (oo, oo) return mtx
[docs] def get_facet_dof_permutations(n_fp, order): """ Prepare DOF permutation vector for each possible facet orientation. """ from sfepy.base.base import dict_to_array if n_fp == 2: mtx = make_line_matrix(order) ori_map = ori_line_to_iter fo = order - 1 elif n_fp == 3: mtx = make_triangle_matrix(order) ori_map = ori_triangle_to_iter fo = order - 2 elif n_fp == 4: mtx = make_square_matrix(order) ori_map = {} for key, val in six.iteritems(_quad_ori_groups): ori_map[key] = ori_square_to_iter[val] fo = order - 1 else: raise ValueError('unsupported number of facet points! (%d)' % n_fp) dof_perms = {} for key, itfun in six.iteritems(ori_map): dof_perms[key] = [mtx[ii] for ii in itfun(fo)] dof_perms = dict_to_array(dof_perms) return dof_perms
if __name__ == '__main__': order = 5 mtx = make_triangle_matrix(order) print(mtx) oo = order - 2 print([mtx[ir, ic] for ir, ic in ori_triangle_to_iter[0](oo)]) print([mtx[ir, ic] for ir, ic in ori_triangle_to_iter[1](oo)]) print([mtx[ir, ic] for ir, ic in ori_triangle_to_iter[3](oo)]) print([mtx[ir, ic] for ir, ic in ori_triangle_to_iter[4](oo)]) print([mtx[ir, ic] for ir, ic in ori_triangle_to_iter[6](oo)]) print([mtx[ir, ic] for ir, ic in ori_triangle_to_iter[7](oo)]) order = 4 mtx = make_square_matrix(order) print(mtx) oo = order - 1 print([mtx[ir, ic] for ir, ic in ori_square_to_iter[0](oo)]) print([mtx[ir, ic] for ir, ic in ori_square_to_iter[7](oo)]) print([mtx[ir, ic] for ir, ic in ori_square_to_iter[11](oo)]) print([mtx[ir, ic] for ir, ic in ori_square_to_iter[30](oo)]) print([mtx[ir, ic] for ir, ic in ori_square_to_iter[33](oo)]) print([mtx[ir, ic] for ir, ic in ori_square_to_iter[52](oo)]) print([mtx[ir, ic] for ir, ic in ori_square_to_iter[56](oo)]) print([mtx[ir, ic] for ir, ic in ori_square_to_iter[63](oo)])