# Source code for sfepy.discrete.fem.linearizer

```"""
Linearization of higher order solutions for the purposes of visualization.
"""
from __future__ import absolute_import
import numpy as nm

from sfepy.linalg import dot_sequences
from sfepy.discrete.fem.refine import refine_reference
from six.moves import range

[docs]def get_eval_dofs(dofs, dof_conn, ps, ori=None):
"""
Get default function for evaluating field DOFs given a list of elements and
reference element coordinates.
"""
def _eval(iels, rx):
edofs = dofs[dof_conn[iels]]

if ori is not None:
eori = ori[iels]

else:
eori = None

bf = ps.eval_base(rx, ori=eori, force_axis=True)[...,0,:]
rvals = dot_sequences(bf, edofs)

return rvals

return _eval

[docs]def get_eval_coors(coors, conn, ps):
"""
Get default function for evaluating physical coordinates given a list of
elements and reference element coordinates.
"""
def _eval(iels, rx):
ecoors = coors[conn[iels]]
aux = ecoors.transpose((0, 2, 1))

bf = ps.eval_base(rx).squeeze()
phys_coors = nm.dot(aux, bf.T).transpose((0, 2, 1))
return phys_coors

return _eval

[docs]def create_output(eval_dofs, eval_coors, n_el, ps, min_level=0, max_level=2,
eps=1e-4):
"""
Create mesh with linear elements that approximates DOFs returned by
`eval_dofs()` corresponding to a higher order approximation with a relative
precision given by `eps`. The DOFs are evaluated in physical coordinates
returned by `eval_coors()`.
"""

def _get_msd(iels, rx, ree):
rvals = eval_dofs(iels, rx)
rng = rvals.max() - rvals.min()
n_components = rvals.shape[-1]

msd = 0.0
for ic in range(n_components):
rval = rvals[..., ic]

sd = rval[:, ree]
# ~ max. second derivative.
msd += nm.abs(sd[..., 0] + sd[..., 2]
- 2.0 * sd[..., 1]).max(axis=-1)

msd /= n_components

return msd, rng

rx0 = ps.geometry.coors

rc0 = ps.geometry.conn[None, :]
rx, rc, ree = refine_reference(ps.geometry, 1)

factor = rc.shape[0] / rc0.shape[0]

iels = nm.arange(n_el)
msd, rng = _get_msd(iels, rx, ree)
eps_r = rng * eps
flag = msd > eps_r

iels0 = flag0 = None

coors = []
conns = []
vdofs = []

inod = 0
for level in range(max_level + 1):
if level < min_level:
flag.fill(True) # Force refinement everywhere.

elif level == max_level:
# Last level - take everything.
flag.fill(False)

# Deal with finished elements.
if flag0 is not None:
ii = nm.searchsorted(iels0, iels)
expand_flag0 = flag0[ii].repeat(factor, axis=1)

else:
expand_flag0 = nm.ones_like(flag)

ie, ir = nm.where((flag == False) & (expand_flag0 == True))
if len(ie):
uie, iies = nm.unique(ie, return_inverse=True)

# Each (sub-)element has own coordinates - no shared vertices.
xes = eval_coors(iels[uie], rx0)
des = eval_dofs(iels[uie], rx0)

# Vectorize (how??) or use cython?
cc = []
vd = []
for ii, iie in enumerate(iies):
ce = rc0[ir[ii]]

xe = xes[iie]
cc.append(xe[ce])

de = des[iie]
vd.append(de[ce])

cc = nm.vstack(cc)
vd = nm.vstack(vd)

nc = cc.shape[0]
np = rc0.shape[1]
conn = nm.arange(nc, dtype=nm.int32).reshape((nc // np, np))

coors.append(cc)
conns.append(conn + inod)
vdofs.append(vd)

inod += nc

if not flag.any():
break

iels0 = iels
flag0 = flag

# Deal with elements to refine.
if level < max_level:
eflag = flag.sum(axis=1, dtype=bool)
iels = iels[eflag]

rc0 = rc
rx0 = rx
rx, rc, ree = refine_reference(ps.geometry, level + 2)

msd, rng = _get_msd(iels, rx, ree)
eps_r = rng * eps
flag = msd > eps_r

all_coors = nm.concatenate(coors, axis=0)
conn = nm.concatenate(conns, axis=0)
all_vdofs = nm.concatenate(vdofs, axis=0)

mat_ids = nm.zeros(conn.shape[0], dtype=nm.int32)

return level, all_coors, conn, all_vdofs, mat_ids
```