"""
Basic uniform mesh refinement functions.
"""
from __future__ import absolute_import
import numpy as nm
from sfepy.discrete.fem import Mesh
from six.moves import range
[docs]
def refine_1_2(mesh_in):
"""
Refines 1D mesh by cutting each element in half
"""
if not nm.all(mesh_in.coors[:-1] <= mesh_in.coors[1:]):
raise ValueError("1D Mesh for refinement must have sorted coors array")
cmesh = mesh_in.cmesh
c_centres = cmesh.get_centroids(cmesh.tdim)
new_coors = nm.zeros((2*mesh_in.n_nod - 1, 1))
new_coors[0::2] = mesh_in.coors
new_coors[1::2] = c_centres
n_nod = len(new_coors)
new_conn = nm.arange(n_nod, dtype=nm.int32).repeat(2)[1:-1].reshape((-1, 2))
new_mat_id = cmesh.cell_groups.repeat(2)
mesh = Mesh.from_data(mesh_in.name + '_r', new_coors, None, [new_conn],
[new_mat_id], mesh_in.descs)
return mesh
[docs]
def refine_2_3(mesh_in):
"""
Refines mesh out of triangles by cutting cutting each edge in half
and making 4 new finer triangles out of one coarser one.
"""
cmesh = mesh_in.cmesh
# Unique edge centres.
e_centres = cmesh.get_centroids(cmesh.tdim - 1)
# New coordinates after the original ones.
coors = nm.r_[mesh_in.coors, e_centres]
o1 = mesh_in.n_nod
cc = cmesh.get_conn(cmesh.tdim, cmesh.tdim - 1)
conn = mesh_in.get_conn('2_3')
n_el = conn.shape[0]
e_nodes = cc.indices.reshape((n_el, 3)) + o1
c = nm.c_[conn, e_nodes].T
new_conn = nm.vstack([c[0], c[3], c[5],
c[3], c[4], c[5],
c[1], c[4], c[3],
c[2], c[5], c[4]]).T
new_conn = new_conn.reshape((4 * n_el, 3))
new_mat_id = cmesh.cell_groups.repeat(4)
mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
[new_mat_id], mesh_in.descs )
return mesh
[docs]
def refine_2_4(mesh_in):
"""
Refines mesh out of quadrilaterals by cutting cutting each edge in
half and making 4 new finer quadrilaterals out of one coarser one.
"""
cmesh = mesh_in.cmesh
# Unique edge centres.
e_centres = cmesh.get_centroids(cmesh.tdim - 1)
# Unique element centres.
centres = cmesh.get_centroids(cmesh.tdim)
# New coordinates after the original ones.
coors = nm.r_[mesh_in.coors, e_centres, centres]
o1 = mesh_in.n_nod
o2 = o1 + e_centres.shape[0]
cc = cmesh.get_conn(cmesh.tdim, cmesh.tdim - 1)
conn = mesh_in.get_conn('2_4')
n_el = conn.shape[0]
e_nodes = cc.indices.reshape((n_el, 4)) + o1
nodes = nm.arange(n_el) + o2
c = nm.c_[conn, e_nodes, nodes].T
new_conn = nm.vstack([c[0], c[4], c[8], c[7],
c[1], c[5], c[8], c[4],
c[2], c[6], c[8], c[5],
c[3], c[7], c[8], c[6]]).T
new_conn = new_conn.reshape((4 * n_el, 4))
new_mat_id = cmesh.cell_groups.repeat(4)
mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
[new_mat_id], mesh_in.descs )
return mesh
[docs]
def refine_3_4(mesh_in):
"""
Refines tetrahedra by cutting each edge in half and making 8 new
finer tetrahedra out of one coarser one. Old nodal coordinates come
first in `coors`, then the new ones. The new tetrahedra are similar
to the old one, no degeneration is supposed to occur as at most 3
congruence classes of tetrahedra appear, even when re-applied
iteratively (provided that `conns` are not modified between two
applications - ordering of vertices in tetrahedra matters not only
for positivity of volumes).
References:
- Juergen Bey: Simplicial grid refinement: on Freudenthal s algorithm and
the optimal number of congruence classes, Numer.Math. 85 (2000),
no. 1, 1--29, or
- Juergen Bey: Tetrahedral grid refinement, Computing 55 (1995),
no. 4, 355--378, or
http://citeseer.ist.psu.edu/bey95tetrahedral.html
"""
cmesh = mesh_in.cmesh
# Unique edge centres.
e_centres = cmesh.get_centroids(cmesh.tdim - 2)
# New coordinates after the original ones.
coors = nm.r_[mesh_in.coors, e_centres]
o1 = mesh_in.n_nod
cc = cmesh.get_conn(cmesh.tdim, cmesh.tdim - 2)
conn = mesh_in.get_conn('3_4')
n_el = conn.shape[0]
e_nodes = cc.indices.reshape((n_el, 6)) + o1
c = nm.c_[conn, e_nodes].T
new_conn = nm.vstack([c[0], c[4], c[6], c[7],
c[4], c[1], c[5], c[8],
c[6], c[5], c[2], c[9],
c[7], c[8], c[9], c[3],
c[4], c[6], c[7], c[8],
c[4], c[6], c[8], c[5],
c[6], c[7], c[8], c[9],
c[6], c[5], c[9], c[8]]).T
new_conn = new_conn.reshape((8 * n_el, 4))
new_mat_id = cmesh.cell_groups.repeat(8)
mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
[new_mat_id], mesh_in.descs )
return mesh
[docs]
def refine_3_8(mesh_in):
"""
Refines hexahedral mesh by cutting cutting each edge in half and
making 8 new finer hexahedrons out of one coarser one.
"""
cmesh = mesh_in.cmesh
# Unique edge centres.
e_centres = cmesh.get_centroids(cmesh.tdim - 2)
# Unique face centres.
f_centres = cmesh.get_centroids(cmesh.tdim - 1)
# Unique element centres.
centres = cmesh.get_centroids(cmesh.tdim)
# New coordinates after the original ones.
coors = nm.r_[mesh_in.coors, e_centres, f_centres, centres]
o1 = mesh_in.n_nod
o2 = o1 + e_centres.shape[0]
o3 = o2 + f_centres.shape[0]
ecc = cmesh.get_conn(cmesh.tdim, cmesh.tdim - 2)
fcc = cmesh.get_conn(cmesh.tdim, cmesh.tdim - 1)
conn = mesh_in.get_conn('3_8')
n_el = conn.shape[0]
st = nm.vstack
e_nodes = ecc.indices.reshape((n_el, 12)) + o1
f_nodes = fcc.indices.reshape((n_el, 6)) + o2
nodes = nm.arange(n_el) + o3
c = nm.c_[conn, e_nodes, f_nodes, nodes].T
new_conn = st([c[0], c[8], c[20], c[11], c[16], c[22], c[26], c[21],
c[1], c[9], c[20], c[8], c[17], c[24], c[26], c[22],
c[2], c[10], c[20], c[9], c[18], c[25], c[26], c[24],
c[3], c[11], c[20], c[10], c[19], c[21], c[26], c[25],
c[4], c[15], c[23], c[12], c[16], c[21], c[26], c[22],
c[5], c[12], c[23], c[13], c[17], c[22], c[26], c[24],
c[6], c[13], c[23], c[14], c[18], c[24], c[26], c[25],
c[7], c[14], c[23], c[15], c[19], c[25], c[26], c[21]]).T
new_conn = new_conn.reshape((8 * n_el, 8))
new_mat_id = cmesh.cell_groups.repeat(8)
mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
[new_mat_id], mesh_in.descs )
return mesh
[docs]
def refine_reference(geometry, level):
"""
Refine reference element given by `geometry`.
Notes
-----
The error edges must be generated in the order of the connectivity
of the previous (lower) level.
"""
from sfepy.discrete.fem import FEDomain
from sfepy.discrete.fem.geometry_element import geometry_data
gcoors, gconn = geometry.coors, geometry.conn
if level == 0:
return gcoors, gconn, None
gd = geometry_data[geometry.name]
conn = nm.array([gd.conn], dtype=nm.int32)
mat_id = conn[:, 0].copy()
mat_id[:] = 0
mesh = Mesh.from_data('aux', gd.coors, None, [conn],
[mat_id], [geometry.name])
domain = FEDomain('aux', mesh)
for ii in range(level):
domain = domain.refine()
coors = domain.mesh.coors
conn = domain.get_conn()
n_el = conn.shape[0]
if geometry.name == '2_3':
aux_conn = conn.reshape((n_el // 4, 4, 3))
ir = [[0, 1, 2], [2, 2, 3], [3, 3, 0]]
ic = [[0, 0, 0], [0, 1, 0], [0, 1, 0]]
elif geometry.name == '2_4':
aux_conn = conn.reshape((n_el // 4, 4, 4))
ir = [[0, 0, 1], [1, 1, 2], [2, 2, 3], [3, 3, 0], [0, 0, 2], [3, 3, 1]]
ic = [[0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0], [1, 2, 1], [1, 2, 1]]
elif geometry.name == '3_4':
aux_conn = conn.reshape((n_el // 8, 8, 4))
ir = [[0, 0, 1], [1, 1, 2], [2, 0, 0], [3, 1, 1], [3, 2, 2], [3, 0, 0]]
ic = [[0, 1, 1], [1, 2, 2], [2, 2, 0], [3, 3, 1], [3, 3, 2], [3, 3, 0]]
elif geometry.name == '3_8':
aux_conn = conn.reshape((n_el // 8, 8, 8))
ir = [[0, 0, 1], [1, 1, 2], [2, 2, 3], [3, 0, 0], [0, 0, 2], [0, 0, 1],
[0, 0, 1], [1, 1, 2], [2, 2, 3], [3, 0, 0], [0, 0, 2], [0, 0, 1],
[4, 4, 5], [5, 5, 6], [6, 6, 7], [7, 4, 4], [4, 4, 6], [4, 4, 5],
[0, 0, 4], [1, 1, 5], [2, 2, 6], [3, 3, 7],
[0, 0, 4], [1, 1, 5], [2, 2, 6], [0, 0, 4],
[0, 0, 4]]
ic = [[0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 3, 0], [1, 2, 1], [3, 2, 1],
[4, 5, 4], [4, 5, 4], [4, 5, 4], [4, 7, 4], [5, 6, 5], [7, 6, 5],
[0, 3, 0], [0, 3, 0], [0, 3, 0], [0, 1, 0], [3, 2, 3], [1, 2, 3],
[0, 4, 0], [0, 4, 0], [0, 4, 0], [0, 4, 0],
[1, 5, 3], [1, 5, 3], [1, 5, 3], [3, 7, 1],
[2, 6, 2]]
else:
raise ValueError('unsupported geometry! (%s)' % geometry.name)
conn = nm.array(conn, dtype=nm.int32)
error_edges = aux_conn[:, ir, ic]
return coors, conn, error_edges