Source code for pyNastran.bdf.mesh_utils.delete_bad_elements

# pylint: disable=C0103
"""
defines:
 - model = delete_bad_shells(model, max_theta=175., max_skew=70., max_aspect_ratio=100.,
                              max_taper_ratio=4.0)
 - eids_to_delete = get_bad_shells(model, xyz_cid0, nid_map, max_theta=175., max_skew=70.,
                                   max_aspect_ratio=100., max_taper_ratio=4.0)

"""
from typing import List
import numpy as np
from cpylog import SimpleLogger
from pyNastran.bdf.bdf import BDF

SIDE_MAP = {}
SIDE_MAP['CHEXA'] = {
    1 : [4, 3, 2, 1],
    2 : [1, 2, 6, 5],
    3 : [2, 3, 7, 6],
    4 : [3, 4, 8, 7],
    5 : [4, 1, 5, 8],
    6 : [5, 6, 7, 8],
}
PIOVER2 = np.pi / 2.
PIOVER3 = np.pi / 3.


[docs]def delete_bad_shells(model: BDF, min_theta: float=0.1, max_theta: float=175., max_skew: float=70., max_aspect_ratio: float=100., max_taper_ratio: float=4.0, max_warping: float=90.): """ Removes bad CQUAD4/CTRIA3 elements Parameters ---------- model : BDF () this should be equivalenced min_theta : float; default=0.1 the maximum interior angle (degrees) max_theta : float; default=175. the maximum interior angle (degrees) max_skew : float; default=70. the maximum skew angle (degrees) max_aspect_ratio : float; default=100. the max aspect ratio taper_ratio : float; default=4.0 the taper ratio; applies to CQUAD4s only max_warping: float: default=20.0 the maximum warp angle (degrees) """ xyz_cid0 = model.get_xyz_in_coord(cid=0, fdtype='float32') nid_map = get_node_map(model) eids_to_delete = get_bad_shells(model, xyz_cid0, nid_map, min_theta=min_theta, max_theta=max_theta, max_skew=max_skew, max_aspect_ratio=max_aspect_ratio, max_taper_ratio=max_taper_ratio, max_warping=max_warping) for eid in eids_to_delete: del model.elements[eid] model.log.info('deleted %s bad CTRIA3/CQUAD4s' % len(eids_to_delete)) model.validate() return model
[docs]def get_node_map(model): """gets an nid->inid mapper""" nid_map = {} for i, nid in enumerate(sorted(model.nodes.keys())): nid_map[nid] = i return nid_map
[docs]def get_bad_shells(model: BDF, xyz_cid0, nid_map, min_theta: float=0.1, max_theta: float=175., max_skew: float=70., max_aspect_ratio: float=100., max_taper_ratio: float=4.0, max_warping: float=60.0) -> List[int]: """ Get the bad shell elements Parameters ---------- model : BDF() the model object xyz_cid0 : (N, 3) float ndarray the xyz coordinates in cid=0 nid_map : dict[nid] : index nid : int the node id index : int the index of the node id in xyz_cid0 min_theta : float; default=0.1 the maximum interior angle (degrees) max_theta : float; default=175. the maximum interior angle (degrees) max_skew : float; default=70. the maximum skew angle (degrees) max_aspect_ratio : float; default=100. the max aspect ratio taper_ratio : float; default=4.0 the taper ratio; applies to CQUAD4s only max_warping: float: default=60.0 the maximum warp angle (degrees) Returns ------- eids_failed : List[int] element ids that fail the criteria shells with a edge length=0.0 are automatically added """ log = model.log min_theta_quad = min_theta min_theta_tri = min_theta min_theta_quad = np.radians(min_theta_quad) min_theta_tri = np.radians(min_theta_tri) max_theta = np.radians(max_theta) max_skew = np.radians(max_skew) max_warping = np.radians(max_warping) eids_failed = [] for eid, element in sorted(model.elements.items()): if element.type == 'CQUAD4': node_ids = element.node_ids #pid = element.Pid() #for nid in node_ids: #if nid is not None: #nid_to_pid_map[nid].append(pid) #self.eid_to_nid_map[eid] = node_ids n1, n2, n3, n4 = [nid_map[nid] for nid in node_ids] p1 = xyz_cid0[n1, :] p2 = xyz_cid0[n2, :] p3 = xyz_cid0[n3, :] p4 = xyz_cid0[n4, :] if _is_bad_quad(eid, p1, p2, p3, p4, log, max_aspect_ratio, min_theta_quad, max_theta, max_skew, max_taper_ratio, max_warping): eids_failed.append(eid) elif element.type == 'CTRIA3': node_ids = element.node_ids #pid = element.Pid() #self.eid_to_nid_map[eid] = node_ids #for nid in node_ids: #if nid is not None: #nid_to_pid_map[nid].append(pid) n1, n2, n3 = [nid_map[nid] for nid in node_ids] p1 = xyz_cid0[n1, :] p2 = xyz_cid0[n2, :] p3 = xyz_cid0[n3, :] if _is_bad_tri(eid, p1, p2, p3, log, max_aspect_ratio, min_theta_tri, max_theta, max_skew): eids_failed.append(eid) return eids_failed
[docs]def _is_bad_quad(eid: int, p1, p2, p3, p4, log: SimpleLogger, max_aspect_ratio: float, min_theta_quad: float, max_theta: float, max_skew: float, max_taper_ratio: float, max_warping: float) -> bool: """identifies if a CQUAD4 has poor quality""" is_bad_quad = True v21 = p2 - p1 v32 = p3 - p2 v43 = p4 - p3 v14 = p1 - p4 #aspect_ratio = max(p12, p23, p34, p14) / max(p12, p23, p34, p14) lengths = np.linalg.norm([v21, v32, v43, v14], axis=1) #assert len(lengths) == 3, lengths length_min = lengths.min() if length_min == 0.0: log.debug('eid=%s failed length_min check; length_min=%s' % (eid, length_min)) return is_bad_quad # ------------------------------------------------------------------------------- aspect_ratio = lengths.max() / length_min if aspect_ratio > max_aspect_ratio: log.debug('eid=%s failed aspect_ratio check; AR=%.2f' % (eid, aspect_ratio)) return is_bad_quad # ------------------------------------------------------------------------------- p12 = (p1 + p2) / 2. p23 = (p2 + p3) / 2. p34 = (p3 + p4) / 2. p14 = (p4 + p1) / 2. normal = np.cross(p3 - p1, p4 - p2) # v42 x v31 # e3 # 4-------3 # | | # |e4 | e2 # 1-------2 # e1 e13 = p34 - p12 e42 = p23 - p14 norm_e13 = np.linalg.norm(e13) norm_e42 = np.linalg.norm(e42) cos_skew1 = (e13 @ e42) / (norm_e13 * norm_e42) cos_skew2 = (e13 @ -e42) / (norm_e13 * norm_e42) skew = np.pi / 2. - np.abs(np.arccos(np.clip([cos_skew1, cos_skew2], -1., 1.))).min() if skew > max_skew: log.debug('eid=%s failed max_skew check; skew=%.2f' % (eid, np.degrees(skew))) return is_bad_quad # ------------------------------------------------------------------------------- area1 = 0.5 * np.linalg.norm(np.cross(-v14, v21)) # v41 x v21 area2 = 0.5 * np.linalg.norm(np.cross(-v21, v32)) # v12 x v32 area3 = 0.5 * np.linalg.norm(np.cross(v43, v32)) # v43 x v32 area4 = 0.5 * np.linalg.norm(np.cross(v14, -v43)) # v14 x v34 aavg = (area1 + area2 + area3 + area4) / 4. taper_ratioi = (abs(area1 - aavg) + abs(area2 - aavg) + abs(area3 - aavg) + abs(area4 - aavg)) / aavg if taper_ratioi > max_taper_ratio: log.debug('eid=%s failed taper_ratio check; taper=%.2f' % (eid, taper_ratioi)) return is_bad_quad # ------------------------------------------------------------------------------- #if 0: #areai = 0.5 * np.linalg.norm(normal) ## still kind of in development ## ## the ratio of the ideal area to the actual area ## this is an hourglass check #areas = [ #np.linalg.norm(np.cross(-v14, v21)), # v41 x v21 #np.linalg.norm(np.cross(v32, -v21)), # v32 x v12 #np.linalg.norm(np.cross(v43, -v32)), # v43 x v23 #np.linalg.norm(np.cross(v14, v43)), # v14 x v43 #] #area_ratioi1 = areai / min(areas) #area_ratioi2 = max(areas) / areai #area_ratioi = max(area_ratioi1, area_ratioi2) # ------------------------------------------------------------------------------- # ixj = k # dot the local normal with the normal vector # then take the norm of that to determine the angle relative to the normal # then take the sign of that to see if we're pointing roughly towards the normal # np.sign(np.linalg.norm(np.dot( # a x b = ab sin(theta) # a x b / ab = sin(theta) # sin(theta) < 0. -> normal is flipped n2 = np.sign(np.cross(v21, v32) @ normal) n3 = np.sign(np.cross(v32, v43) @ normal) n4 = np.sign(np.cross(v43, v14) @ normal) n1 = np.sign(np.cross(v14, v21) @ normal) n = np.array([n1, n2, n3, n4]) theta_additional = np.where(n < 0, 2*np.pi, 0.) norm_v21 = np.linalg.norm(v21) norm_v32 = np.linalg.norm(v32) norm_v43 = np.linalg.norm(v43) norm_v14 = np.linalg.norm(v14) cos_theta1 = (v21 @ -v14) / (norm_v21 * norm_v14) cos_theta2 = (v32 @ -v21) / (norm_v32 * norm_v21) cos_theta3 = (v43 @ -v32) / (norm_v43 * norm_v32) cos_theta4 = (v14 @ -v43) / (norm_v14 * norm_v43) interior_angle = np.arccos(np.clip( [cos_theta1, cos_theta2, cos_theta3, cos_theta4], -1., 1.)) theta = n * interior_angle + theta_additional theta_mini = theta.min() theta_maxi = theta.max() if theta_mini < min_theta_quad: log.debug('eid=%s failed min_theta check; theta=%.2f' % ( eid, np.degrees(theta_mini))) return is_bad_quad if theta_maxi > max_theta: log.debug('eid=%s failed max_theta check; theta=%.2f' % ( eid, np.degrees(theta_maxi))) return is_bad_quad # ------------------------------------------------------------------------------- # warping v31 = p3 - p1 v42 = p4 - p2 v41 = -v14 n123 = np.cross(v21, v31) n134 = np.cross(v31, v41) #v1 o v2 = v1 * v2 cos(theta) cos_warp1 = (n123 @ n134) / (np.linalg.norm(n123) * np.linalg.norm(n134)) # split the quad in the order direction and take the maximum of the two splits # 4---3 # | \ | # | \| # 1---2 n124 = np.cross(v21, v41) n234 = np.cross(v32, v42) cos_warp2 = (n124 @ n234) / (np.linalg.norm(n124) * np.linalg.norm(n234)) max_warpi = np.abs(np.arccos( np.clip([cos_warp1, cos_warp2], -1., 1.))).max() if max_warpi > max_warping: log.debug('eid=%s failed max_warping check; theta=%.2f' % ( eid, np.degrees(max_warpi))) #print('eid=%s theta_min=%-5.2f theta_max=%-5.2f ' #'skew=%-5.2f AR=%-5.2f taper_ratioi=%.2f' % ( #eid, #np.degrees(theta_mini), np.degrees(theta_maxi), #np.degrees(skew), aspect_ratio, taper_ratioi)) is_bad_quad = False return is_bad_quad
[docs]def _is_bad_tri(eid: int, p1, p2, p3, log: SimpleLogger, max_aspect_ratio: float, min_theta_tri: float, max_theta: float, max_skew: float) -> bool: """identifies if a CTRIA3 has poor quality""" is_bad_tri = True v21 = p2 - p1 v32 = p3 - p2 v13 = p1 - p3 lengths = np.linalg.norm([v21, v32, v13], axis=1) length_min = lengths.min() if length_min == 0.0: log.debug('eid=%s failed length_min check; length_min=%s' % ( eid, length_min)) return is_bad_tri #assert len(lengths) == 3, lengths aspect_ratio = lengths.max() / length_min if aspect_ratio > max_aspect_ratio: log.debug('eid=%s failed aspect_ratio check; AR=%s' % (eid, aspect_ratio)) return is_bad_tri cos_theta1 = (v21 @ -v13) / (np.linalg.norm(v21) * np.linalg.norm(v13)) cos_theta2 = (v32 @ -v21) / (np.linalg.norm(v32) * np.linalg.norm(v21)) cos_theta3 = (v13 @ -v32) / (np.linalg.norm(v13) * np.linalg.norm(v32)) theta = np.arccos(np.clip( [cos_theta1, cos_theta2, cos_theta3], -1., 1.)) theta_mini = theta.min() theta_maxi = theta.max() if theta_mini < min_theta_tri: log.debug('eid=%s failed min_theta check; theta=%s' % ( eid, np.degrees(theta_mini))) return is_bad_tri if theta_maxi > max_theta: log.debug('eid=%s failed max_theta check; theta=%s' % ( eid, np.degrees(theta_maxi))) return is_bad_tri # 3 # / \ # e3/ \ e2 # / /\ # / / \ # 1---/----2 # e1 e1 = (p1 + p2) / 2. e2 = (p2 + p3) / 2. e3 = (p3 + p1) / 2. e21 = e2 - e1 e31 = e3 - e1 e32 = e3 - e2 norm_e21 = np.linalg.norm(e21) norm_e31 = np.linalg.norm(e31) norm_e32 = np.linalg.norm(e32) e3_p2 = e3 - p2 e2_p1 = e2 - p1 e1_p3 = e1 - p3 norm_e3_p2 = np.linalg.norm(e3_p2) norm_e2_p1 = np.linalg.norm(e2_p1) norm_e1_p3 = np.linalg.norm(e1_p3) cos_skew1 = (e2_p1 @ e31) / (norm_e2_p1 * norm_e31) cos_skew2 = (e2_p1 @ -e31) / (norm_e2_p1 * norm_e31) cos_skew3 = (e3_p2 @ e21) / (norm_e3_p2 * norm_e21) cos_skew4 = (e3_p2 @ -e21) / (norm_e3_p2 * norm_e21) cos_skew5 = (e1_p3 @ e32) / (norm_e1_p3 * norm_e32) cos_skew6 = (e1_p3 @ -e32) / (norm_e1_p3 * norm_e32) skew = np.pi / 2. - np.abs(np.arccos( np.clip([cos_skew1, cos_skew2, cos_skew3, cos_skew4, cos_skew5, cos_skew6], -1., 1.) )).min() if skew > max_skew: log.debug('eid=%s failed max_skew check; skew=%s' % (eid, np.degrees(skew))) return is_bad_tri is_bad_tri = False # warping doesn't happen to CTRIA3s #print('eid=%s theta_min=%-5.2f theta_max=%-5.2f skew=%-5.2f AR=%-5.2f' % ( #eid, #np.degrees(theta_mini), np.degrees(theta_maxi), #np.degrees(skew), aspect_ratio)) return is_bad_tri
[docs]def element_quality(model, nids=None, xyz_cid0=None, nid_map=None): """ Gets various measures of element quality Parameters ---------- model : BDF() a cross-referenced model nids : (nnodes, ) int ndarray; default=None the nodes of the model in sorted order includes GRID, SPOINT, & EPOINTs xyz_cid0 : (nnodes, 3) float ndarray; default=None the associated global xyz locations nid_map : Dict[nid]->index; default=None a mapper dictionary Returns ------- quality : Dict[name] : (nelements, ) float ndarray Various quality metrics names : min_interior_angle, max_interior_angle, dideal_theta, max_skew_angle, max_aspect_ratio, area_ratio, taper_ratio, min_edge_length values : The result is ``np.nan`` if element type does not define the parameter. For example, CELAS1 doesn't have an aspect ratio. Notes ----- - pulled from nastran_io.py """ if nids is None or xyz_cid0 is None: out = model.get_displacement_index_xyz_cp_cd( fdtype='float64', idtype='int32', sort_ids=True) unused_icd_transform, icp_transform, xyz_cp, nid_cp_cd = out nids = nid_cp_cd[:, 0] xyz_cid0 = model.transform_xyzcp_to_xyz_cid( xyz_cp, nids, icp_transform, cid=0, in_place=False) if nid_map is None: nid_map = {} for i, nid in enumerate(nids): nid_map[nid] = i all_nids = nids del nids # these normals point inwards # 4 # / | \ # / | \ # 3-------2 # \ | / # \ | / # 1 _ctetra_faces = ( (0, 1, 2), # (1, 2, 3), (0, 3, 1), # (1, 4, 2), (0, 3, 2), # (1, 3, 4), (1, 3, 2), # (2, 4, 3), ) # these normals point inwards # # # # # /4-----3 # / / # / 5 / # / \ / # / \ / # 1---------2 _cpyram_faces = ( (0, 1, 2, 3), # (1, 2, 3, 4), (1, 4, 2), # (2, 5, 3), (2, 4, 3), # (3, 5, 4), (0, 3, 4), # (1, 4, 5), (0, 4, 1), # (1, 5, 2), ) # these normals point inwards # /6 # / | \ # / | \ # 3\ | \ # | \ /4-----5 # | \/ / # | / \ / # | / \ / # | / \ / # 1---------2 _cpenta_faces = ( (0, 2, 1), # (1, 3, 2), (3, 4, 5), # (4, 5, 6), (0, 1, 4, 3), # (1, 2, 5, 4), # bottom (1, 2, 5, 4), # (2, 3, 6, 5), # right (0, 3, 5, 2), # (1, 4, 6, 3), # left ) # these normals point inwards # 8----7 # /| /| # / | / | # / 5-/--6 # 4-----3 / # | / | / # | / | / # 1-----2 _chexa_faces = ( (4, 5, 6, 7), # (5, 6, 7, 8), (0, 3, 2, 1), # (1, 4, 3, 2), (1, 2, 6, 5), # (2, 3, 7, 6), (2, 3, 7, 6), # (3, 4, 8, 7), (0, 4, 7, 3), # (1, 5, 8, 4), (0, 6, 5, 4), # (1, 7, 6, 5), ) # quality nelements = len(model.elements) min_interior_angle = np.zeros(nelements, 'float32') max_interior_angle = np.zeros(nelements, 'float32') dideal_theta = np.zeros(nelements, 'float32') max_skew_angle = np.zeros(nelements, 'float32') max_warp_angle = np.zeros(nelements, 'float32') max_aspect_ratio = np.zeros(nelements, 'float32') #area = np.zeros(nelements, 'float32') area_ratio = np.zeros(nelements, 'float32') taper_ratio = np.zeros(nelements, 'float32') min_edge_length = np.zeros(nelements, 'float32') #normals = np.full((nelements, 3), np.nan, 'float32') #nids_list = [] ieid = 0 for unused_eid, elem in sorted(model.elements.items()): if ieid % 5000 == 0 and ieid > 0: print(' map_elements = %i' % ieid) etype = elem.type nids = None inids = None dideal_thetai = np.nan min_thetai = np.nan max_thetai = np.nan #max_thetai = np.nan max_skew = np.nan max_warp = np.nan aspect_ratio = np.nan #areai = np.nan area_ratioi = np.nan taper_ratioi = np.nan min_edge_lengthi = np.nan #normali = np.nan if etype in ['CTRIA3', 'CTRIAR', 'CTRAX3', 'CPLSTN3']: nids = elem.nodes inids = np.searchsorted(all_nids, nids) p1, p2, p3 = xyz_cid0[inids, :] out = tri_quality(p1, p2, p3) (areai, max_skew, aspect_ratio, min_thetai, max_thetai, dideal_thetai, min_edge_lengthi) = out #normali = np.cross(p1 - p2, p1 - p3) elif etype in ['CQUAD4', 'CQUADR', 'CPLSTN4', 'CQUADX4']: nids = elem.nodes inids = np.searchsorted(all_nids, nids) p1, p2, p3, p4 = xyz_cid0[inids, :] out = quad_quality(elem, p1, p2, p3, p4) (areai, taper_ratioi, area_ratioi, max_skew, aspect_ratio, min_thetai, max_thetai, dideal_thetai, min_edge_lengthi, max_warp) = out elif etype == 'CTRIA6': nids = elem.nodes if None in nids: inids = np.searchsorted(all_nids, nids[:3]) nids = nids[:3] p1, p2, p3 = xyz_cid0[inids, :] else: inids = np.searchsorted(all_nids, nids) p1, p2, p3, p4, unused_p5, unused_p6 = xyz_cid0[inids, :] out = tri_quality(p1, p2, p3) (areai, max_skew, aspect_ratio, min_thetai, max_thetai, dideal_thetai, min_edge_lengthi) = out elif etype == 'CQUAD8': nids = elem.nodes if None in nids: inids = np.searchsorted(all_nids, nids[:4]) nids = nids[:4] p1, p2, p3, p4 = xyz_cid0[inids, :] else: inids = np.searchsorted(all_nids, nids) p1, p2, p3, p4 = xyz_cid0[inids[:4], :] out = quad_quality(elem, p1, p2, p3, p4) (areai, taper_ratioi, area_ratioi, max_skew, aspect_ratio, min_thetai, max_thetai, dideal_thetai, min_edge_lengthi, max_warp) = out #normali = np.cross(p1 - p3, p2 - p4) elif etype == 'CSHEAR': nids = elem.nodes inids = np.searchsorted(all_nids, nids) p1, p2, p3, p4 = xyz_cid0[inids, :] out = quad_quality(elem, p1, p2, p3, p4) (areai, taper_ratioi, area_ratioi, max_skew, aspect_ratio, min_thetai, max_thetai, dideal_thetai, min_edge_lengthi, max_warp) = out elif etype == 'CTETRA': nids = elem.nodes if None in nids: nids = nids[:4] inids = np.searchsorted(all_nids, nids) min_thetai, max_thetai, dideal_thetai, min_edge_lengthi = get_min_max_theta( _ctetra_faces, nids, nid_map, xyz_cid0) elif etype == 'CHEXA': nids = elem.nodes if None in nids: nids = nids[:8] inids = np.searchsorted(all_nids, nids) min_thetai, max_thetai, dideal_thetai, min_edge_lengthi = get_min_max_theta( _chexa_faces, nids, nid_map, xyz_cid0) elif etype == 'CPENTA': nids = elem.nodes if None in nids: nids = nids[:6] inids = np.searchsorted(all_nids, nids) min_thetai, max_thetai, dideal_thetai, min_edge_lengthi = get_min_max_theta( _cpenta_faces, nids, nid_map, xyz_cid0) elif etype == 'CPYRAM': # TODO: assuming 5 nids = elem.nodes if None in nids: nids = nids[:5] inids = np.searchsorted(all_nids, nids) min_thetai, max_thetai, dideal_thetai, min_edge_lengthi = get_min_max_theta( _cpyram_faces, nids, nid_map, xyz_cid0) elif etype in ['CELAS2', 'CELAS4', 'CDAMP4']: # these can have empty nodes and have no property # CELAS1: 1/2 GRID/SPOINT and pid # CELAS2: 1/2 GRID/SPOINT, k, ge, and s # CELAS3: 1/2 SPOINT and pid # CELAS4: 1/2 SPOINT and k continue #nids = elem.nodes #assert nids[0] != nids[1] #if None in nids: #assert nids[0] is not None, nids #assert nids[1] is None, nids #nids = [nids[0]] #else: #nids = elem.nodes #assert nids[0] != nids[1] #inids = np.searchsorted(all_nids, nids) elif etype in ['CBUSH', 'CBUSH1D', 'CBUSH2D', 'CELAS1', 'CELAS3', 'CDAMP1', 'CDAMP2', 'CDAMP3', 'CDAMP5', 'CFAST', 'CGAP', 'CVISC']: #nids = elem.nodes #assert nids[0] != nids[1] #assert None not in nids, 'nids=%s\n%s' % (nids, elem) #inids = np.searchsorted(all_nids, nids) continue elif etype in ['CBAR', 'CBEAM']: nids = elem.nodes inids = np.searchsorted(all_nids, nids) p1, p2 = xyz_cid0[inids, :] min_edge_lengthi = np.linalg.norm(p2 - p1) elif etype in ['CROD', 'CTUBE']: nids = elem.nodes inids = np.searchsorted(all_nids, nids) p1, p2 = xyz_cid0[inids, :] min_edge_lengthi = np.linalg.norm(p2 - p1) #nnodes = 2 #dim = 1 elif etype == 'CONROD': nids = elem.nodes inids = np.searchsorted(all_nids, nids) p1, p2 = xyz_cid0[inids, :] min_edge_lengthi = np.linalg.norm(p2 - p1) #------------------------------ # rare #elif etype == 'CIHEX1': #nids = elem.nodes #pid = elem.pid #cell_type = cell_type_hexa8 #inids = np.searchsorted(all_nids, nids) #min_thetai, max_thetai, dideal_thetai, min_edge_lengthi = get_min_max_theta( #_chexa_faces, nids, nid_map, xyz_cid0) #nnodes = 8 #dim = 3 elif etype == 'CHBDYE': continue #self.eid_map[eid] = ieid #eid_solid = elem.eid2 #side = elem.side #element_solid = model.elements[eid_solid] #mapped_inids = SIDE_MAP[element_solid.type][side] #side_inids = [nid - 1 for nid in mapped_inids] #nodes = element_solid.node_ids ##nnodes = len(side_inids) #nids = [nodes[inid] for inid in side_inids] #inids = np.searchsorted(all_nids, nids) #if len(side_inids) == 4: #pass #else: #msg = 'element_solid:\n%s' % (str(element_solid)) #msg += 'mapped_inids = %s\n' % mapped_inids #msg += 'side_inids = %s\n' % side_inids #msg += 'nodes = %s\n' % nodes ##msg += 'side_nodes = %s\n' % side_nodes #raise NotImplementedError(msg) else: #raise NotImplementedError(elem) nelements -= 1 continue #nids_list.append(nnodes) #nids_list.extend(inids) #normals[ieid] = normali #eids_array[ieid] = eid #pids_array[ieid] = pid #dim_array[ieid] = dim #cell_types_array[ieid] = cell_type #cell_offsets_array[ieid] = cell_offset # I assume the problem is here #cell_offset += nnodes + 1 #eid_map[eid] = ieid min_interior_angle[ieid] = min_thetai max_interior_angle[ieid] = max_thetai dideal_theta[ieid] = dideal_thetai max_skew_angle[ieid] = max_skew max_warp_angle[ieid] = max_warp max_aspect_ratio[ieid] = aspect_ratio #area[ieid] = areai area_ratio[ieid] = area_ratioi taper_ratio[ieid] = taper_ratioi min_edge_length[ieid] = min_edge_lengthi ieid += 1 quality = { 'min_interior_angle' : min_interior_angle, 'max_interior_angle' : max_interior_angle, 'dideal_theta' : dideal_theta, 'max_skew_angle' : max_skew_angle, 'max_warp_angle' : max_warp_angle, 'max_aspect_ratio' : max_aspect_ratio, #'area' : area, 'area_ratio' : area_ratio, 'taper_ratio' : taper_ratio, 'min_edge_length' : min_edge_length, } return quality
[docs]def tri_quality(p1, p2, p3): """gets the quality metrics for a tri""" e1 = (p1 + p2) / 2. e2 = (p2 + p3) / 2. e3 = (p3 + p1) / 2. # 3 # / \ # e3/ \ e2 # / /\ # / / \ # 1---/----2 # e1 e21 = e2 - e1 e31 = e3 - e1 e32 = e3 - e2 e3_p2 = e3 - p2 e2_p1 = e2 - p1 e1_p3 = e1 - p3 v21 = p2 - p1 v32 = p3 - p2 v13 = p1 - p3 length21 = np.linalg.norm(v21) length32 = np.linalg.norm(v32) length13 = np.linalg.norm(v13) min_edge_length = min(length21, length32, length13) area = 0.5 * np.linalg.norm(np.cross(v21, v13)) ne31 = np.linalg.norm(e31) ne21 = np.linalg.norm(e21) ne32 = np.linalg.norm(e32) ne2_p1 = np.linalg.norm(e2_p1) ne3_p2 = np.linalg.norm(e3_p2) ne1_p3 = np.linalg.norm(e1_p3) cos_skew1 = (e2_p1 @ e31) / (ne2_p1 * ne31) cos_skew2 = (e2_p1 @ -e31) / (ne2_p1 * ne31) cos_skew3 = (e3_p2 @ e21) / (ne3_p2 * ne21) cos_skew4 = (e3_p2 @ -e21) / (ne3_p2 * ne21) cos_skew5 = (e1_p3 @ e32) / (ne1_p3 * ne32) cos_skew6 = (e1_p3 @ -e32) / (ne1_p3 * ne32) max_skew = np.pi / 2. - np.abs(np.arccos(np.clip([ cos_skew1, cos_skew2, cos_skew3, cos_skew4, cos_skew5, cos_skew6], -1., 1.))).min() lengths = np.linalg.norm([v21, v32, v13], axis=1) #assert len(lengths) == 3, lengths length_min = lengths.min() if length_min == 0.0: aspect_ratio = np.nan # assume length_min = length21 = nan, so: # cos_theta1 = nan # thetas = [nan, b, c] # min_theta = max_theta = dideal_theta = nan min_theta = np.nan max_theta = np.nan dideal_theta = np.nan else: aspect_ratio = lengths.max() / length_min cos_theta1 = (v21 @ -v13) / (length21 * length13) cos_theta2 = (v32 @ -v21) / (length32 * length21) cos_theta3 = (v13 @ -v32) / (length13 * length32) thetas = np.arccos(np.clip([cos_theta1, cos_theta2, cos_theta3], -1., 1.)) min_theta = thetas.min() max_theta = thetas.max() dideal_theta = max(max_theta - PIOVER3, PIOVER3 - min_theta) #theta_deg = np.degrees(np.arccos(max_cos_theta)) #if theta_deg < 60.: #print('p1=%s' % xyz_cid0[p1, :]) #print('p2=%s' % xyz_cid0[p2, :]) #print('p3=%s' % xyz_cid0[p3, :]) #print('theta1=%s' % np.degrees(np.arccos(cos_theta1))) #print('theta2=%s' % np.degrees(np.arccos(cos_theta2))) #print('theta3=%s' % np.degrees(np.arccos(cos_theta3))) #print('max_theta=%s' % theta_deg) #asdf return area, max_skew, aspect_ratio, min_theta, max_theta, dideal_theta, min_edge_length
[docs]def quad_quality(element, p1, p2, p3, p4): """gets the quality metrics for a quad""" v21 = p2 - p1 v32 = p3 - p2 v43 = p4 - p3 v14 = p1 - p4 length21 = np.linalg.norm(v21) length32 = np.linalg.norm(v32) length43 = np.linalg.norm(v43) length14 = np.linalg.norm(v14) min_edge_length = min(length21, length32, length43, length14) p12 = (p1 + p2) / 2. p23 = (p2 + p3) / 2. p34 = (p3 + p4) / 2. p14 = (p4 + p1) / 2. v31 = p3 - p1 v42 = p4 - p2 normal = np.cross(v31, v42) area = 0.5 * np.linalg.norm(normal) # still kind of in development # # the ratio of the ideal area to the actual area # this is an hourglass check areas = [ np.linalg.norm(np.cross(-v14, v21)), # v41 x v21 np.linalg.norm(np.cross(v32, -v21)), # v32 x v12 np.linalg.norm(np.cross(v43, -v32)), # v43 x v23 np.linalg.norm(np.cross(v14, v43)), # v14 x v43 ] # # for: # area=1; area1=0.5 -> area_ratioi1=2.0; area_ratio=2.0 # area=1; area1=2.0 -> area_ratioi2=2.0; area_ratio=2.0 min_area = min(areas) if min_area == 0.: print('nan min area; area=%g areas=%s:\n%s' % (area, areas, element)) #nodes = element.nodes #print(' n_%i = %s' % (nodes[0], p1)) #print(' n_%i = %s' % (nodes[1], p2)) #print(' n_%i = %s' % (nodes[2], p3)) #print(' n_%i = %s' % (nodes[3], p4)) area_ratio = np.nan #raise RuntimeError('bad quad...') else: area_ratioi1 = area / min_area area_ratioi2 = max(areas) / area area_ratio = max(area_ratioi1, area_ratioi2) area1 = 0.5 * areas[0] # v41 x v21 area2 = 0.5 * np.linalg.norm(np.cross(-v21, v32)) # v12 x v32 area3 = 0.5 * np.linalg.norm(np.cross(v43, v32)) # v43 x v32 area4 = 0.5 * np.linalg.norm(np.cross(v14, -v43)) # v14 x v34 aavg = (area1 + area2 + area3 + area4) / 4. taper_ratio = (abs(area1 - aavg) + abs(area2 - aavg) + abs(area3 - aavg) + abs(area4 - aavg)) / aavg # e3 # 4-------3 # | | # |e4 | e2 # 1-------2 # e1 e13 = p34 - p12 e42 = p23 - p14 ne42 = np.linalg.norm(e42) ne13 = np.linalg.norm(e13) cos_skew1 = (e13 @ e42) / (ne13 * ne42) cos_skew2 = (e13 @ -e42) / (ne13 * ne42) max_skew = np.pi / 2. - np.abs(np.arccos( np.clip([cos_skew1, cos_skew2], -1., 1.))).min() #aspect_ratio = max(p12, p23, p34, p14) / max(p12, p23, p34, p14) lengths = np.linalg.norm([v21, v32, v43, v14], axis=1) #assert len(lengths) == 3, lengths aspect_ratio = lengths.max() / lengths.min() cos_theta1 = (v21 @ -v14) / (length21 * length14) cos_theta2 = (v32 @ -v21) / (length32 * length21) cos_theta3 = (v43 @ -v32) / (length43 * length32) cos_theta4 = (v14 @ -v43) / (length14 * length43) #max_thetai = np.arccos([cos_theta1, cos_theta2, cos_theta3, cos_theta4]).max() # dot the local normal with the normal vector # then take the norm of that to determine the angle relative to the normal # then take the sign of that to see if we're pointing roughly towards the normal # # np.sign(np.linalg.norm(np.dot( # a x b = ab sin(theta) # a x b / ab = sin(theta) # sin(theta) < 0. -> normal is flipped normal2 = np.sign(np.cross(v21, v32) @ normal) normal3 = np.sign(np.cross(v32, v43) @ normal) normal4 = np.sign(np.cross(v43, v14) @ normal) normal1 = np.sign(np.cross(v14, v21) @ normal) n = np.array([normal1, normal2, normal3, normal4]) theta_additional = np.where(n < 0, 2*np.pi, 0.) theta = n * np.arccos(np.clip( [cos_theta1, cos_theta2, cos_theta3, cos_theta4], -1., 1.)) + theta_additional min_theta = theta.min() max_theta = theta.max() dideal_theta = max(max_theta - PIOVER2, PIOVER2 - min_theta) #print('theta_max = ', theta_max) # warp angle # split the quad and find the normals of each triangl # find the angle between the two triangles # # 4---3 # | / | # |/ | # 1---2 # v41 = -v14 n123 = np.cross(v21, v31) n134 = np.cross(v31, v41) #v1 o v2 = v1 * v2 cos(theta) cos_warp1 = (n123 @ n134) / (np.linalg.norm(n123) * np.linalg.norm(n134)) # split the quad in the order direction and take the maximum of the two splits # 4---3 # | \ | # | \| # 1---2 n124 = np.cross(v21, v41) n234 = np.cross(v32, v42) cos_warp2 = (n124 @ n234) / (np.linalg.norm(n124) * np.linalg.norm(n234)) max_warp = np.abs(np.arccos( np.clip([cos_warp1, cos_warp2], -1., 1.))).max() out = (area, taper_ratio, area_ratio, max_skew, aspect_ratio, min_theta, max_theta, dideal_theta, min_edge_length, max_warp) return out
[docs]def get_min_max_theta(faces, all_node_ids, nid_map, xyz_cid0): """get the min/max thetas for CTETRA, CPENTA, CHEXA, CPYRAM""" cos_thetas = [] ideal_theta = [] #print('faces =', faces) #assert len(faces) > 0, 'faces=%s nids=%s' % (faces, all_node_ids) for face in faces: if len(face) == 3: node_ids = all_node_ids[face[0]], all_node_ids[face[1]], all_node_ids[face[2]] n1, n2, n3 = [nid_map[nid] for nid in node_ids[:3]] v21 = xyz_cid0[n2, :] - xyz_cid0[n1, :] v32 = xyz_cid0[n3, :] - xyz_cid0[n2, :] v13 = xyz_cid0[n1, :] - xyz_cid0[n3, :] length21 = np.linalg.norm(v21) length32 = np.linalg.norm(v32) length13 = np.linalg.norm(v13) min_edge_length = min(length21, length32, length13) cos_theta1 = (v21 @ -v13) / (length21 * length13) cos_theta2 = (v32 @ -v21) / (length32 * length21) cos_theta3 = (v13 @ -v32) / (length13 * length32) cos_thetas.extend([cos_theta1, cos_theta2, cos_theta3]) ideal_theta.extend([PIOVER3, PIOVER3, PIOVER3]) elif len(face) == 4: try: node_ids = (all_node_ids[face[0]], all_node_ids[face[1]], all_node_ids[face[2]], all_node_ids[face[3]]) except KeyError: print(face) print(node_ids) raise n1, n2, n3, n4 = [nid_map[nid] for nid in node_ids[:4]] v21 = xyz_cid0[n2, :] - xyz_cid0[n1, :] v32 = xyz_cid0[n3, :] - xyz_cid0[n2, :] v43 = xyz_cid0[n4, :] - xyz_cid0[n3, :] v14 = xyz_cid0[n1, :] - xyz_cid0[n4, :] length21 = np.linalg.norm(v21) length32 = np.linalg.norm(v32) length43 = np.linalg.norm(v43) length14 = np.linalg.norm(v14) min_edge_length = min(length21, length32, length43, length14) cos_theta1 = (v21 @ -v14) / (length21 * length14) cos_theta2 = (v32 @ -v21) / (length32 * length21) cos_theta3 = (v43 @ -v32) / (length43 * length32) cos_theta4 = (v14 @ -v43) / (length14 * length43) cos_thetas.extend([cos_theta1, cos_theta2, cos_theta3, cos_theta4]) ideal_theta.extend([PIOVER2, PIOVER2, PIOVER2, PIOVER2]) else: raise NotImplementedError(face) thetas = np.arccos(cos_thetas) ideal_theta = np.array(ideal_theta) ideal_thetai = max((thetas - ideal_theta).max(), (ideal_theta - thetas).min()) min_thetai = thetas.min() max_thetai = thetas.max() return min_thetai, max_thetai, ideal_thetai, min_edge_length