+"""This module stores base classes and functions.
+
+Attributes:
+ RESFILEDS (dict): stores results fields under results type
+
+ key (str):
+ - 'displ': displacement
+ - 'stress': stress
+ - 'strain': strain
+ - 'force': force
+ value (list):
+ - 'displ' = 'ux','uy','uz','utot'
+ - 'stress' = 'Sx','Sy','Sz','Sxy','Syz','Szx','Seqv','S1','S2','S3'
+ - 'strain' = 'ex','ey','ez','exy','eyz','ezx','eeqv','e1','e2','e3'
+ - 'force' = 'fx','fy','fz'
+ FIELDTYPE (dict): the inverse dict of RESFIELDS
+
+ For example: key 'ux' --> value: 'displ'
+"""
+from math import ceil
+from . import environment
+
+RESFIELDS = {}
+RESFIELDS['displ'] = 'ux,uy,uz,utot'.split(',')
+RESFIELDS['stress'] = 'Sx,Sy,Sz,Sxy,Syz,Szx,Seqv,S1,S2,S3'.split(',')
+RESFIELDS['strain'] = 'ex,ey,ez,exy,eyz,ezx,eeqv,e1,e2,e3'.split(',')
+RESFIELDS['force'] = 'fx,fy,fz'.split(',')
+
+#FIELDTYPE is a dict that inverts the RESFIELDS dictionary mapping.
+FIELDTYPE = {}
+for (k, v) in RESFIELDS.items():
+ for vi in v:
+ FIELDTYPE[vi] = k
+
+[docs]class Idobj(object):
+ """Makes an object that stores an id number.
+
+ This is a base class for nodes, lines, areas etc.
+
+ Attributes:
+ id (int): the unique id number for the item
+ """
+
+ def __init__(self):
+ self.id = -1
+
+
+
+ def __hash__(self):
+ """Sets the hash of the item to the id number.
+
+ This allows one to make sets of these items.
+ """
+ return self.id
+
+[docs]def plot_finish(plt, fname, display):
+ """Display and or save plot."""
+ if fname != '':
+ # save the image
+ fname += '.png'
+ if environment.DPI != None:
+ plt.savefig(fname, dpi=environment.DPI, bbox_inches='tight')
+ else:
+ plt.savefig(fname, bbox_inches='tight')
+ print('File %s was saved.' % fname)
+
+ if display:
+ plt.tight_layout()
+ plt.show()
+
+ # remove all figures
+ plt.close()
+
+[docs]def plot_set_bounds(plt, axials, radials):
+ """Sets the axial and radial bounds of the shown plot."""
+ vert = max(radials) - min(radials)
+ horiz = max(axials) - min(axials)
+ vadder = (vert)/5
+ hadder = (horiz)/5
+ (vmax, vmin) = (max(radials)+vadder, min(radials)-vadder)
+ (hmax, hmin) = (max(axials)+hadder, min(axials)-hadder)
+ plt.xlim(hmin, hmax)
+ plt.ylim(vmin, vmax)
+
+[docs]class Itemlist(list):
+ """Makes a custom list used to store lists of non-mesh items.
+
+ non-mesh items = Point, Line, Arc, Area, Part, Components, Loads
+ This class allows us to automatically assign id numbers to items in the
+ list. Thes id numbers are needed when we send geometry out for meshing.
+ All instances of this class start as empty lists.
+ """
+
+ def __init__(self):
+ super().__init__() # these lists start empty
+
+[docs] def get_ids(self):
+ """Returns a list of ids. Loops through all items in self."""
+ return [a.id for a in self]
+
+[docs] def get_next_id(self):
+ """Returns the next unused id."""
+ ids = self.get_ids()
+ minid = 0
+ if len(ids) == 0:
+ return minid # list is empty so return the minid number
+ else:
+ ids = sorted(ids)
+ maxid = ids[-1]
+ unused = list(set(list(range(minid, maxid+2))) - set(ids))
+ return unused[0] # return first available id number
+
+[docs] def append(self, item):
+ """Adds an item to the list and sets the item's id.
+
+ Args:
+ item (Point or Line or Arc or Area or Part): item to add to the list
+ """
+ idnum = self.get_next_id()
+ item.set_id(idnum)
+ super().append(item)
+ return item
+
+
+[docs]class Meshlist(list):
+ """Makes a custom list used to store lists of nodes and elements.
+
+ All instances of this class start as empty lists.
+ """
+
+ def __init__(self):
+ list.__init__([])
+
+[docs] def get_minid(self):
+ """Returns the min id number in the list."""
+ ids = [x.id for x in self]
+ return min(ids)
+
+[docs] def get_maxid(self):
+ """Returns the max id number in the list."""
+ ids = [x.id for x in self]
+ return max(ids)
+
+[docs] def idget(self, idnum):
+ """Returns an item with the passed id number.
+
+ Args:
+ idnum (int): passed id number, we want the item that has this
+
+ Returns:
+ item or None: returns item if found, None if not found
+ """
+ for item in self:
+ if item.id == idnum:
+ return item
+ #print ('ID %i was not found!' % (idnum))
+ return None
+
+[docs] def set_minid(self, val):
+ """Sets the min id to the passed val.
+
+ Args:
+ val (int): New min id number in the list
+ """
+ print('Re-indexing elements to min element number of: %i' % val)
+ print('Old min:%i max:%i' % (self.get_minid(), self.get_maxid()))
+ minid = self.get_minid()
+ offset = minid - val
+ for item in self:
+ item.id -= offset
+ print('New min:%i max:%i' % (self.get_minid(), self.get_maxid()))
+
+
+[docs]def chunk_list(inlist, size):
+ """Returns a list of lists where each list <= size length.
+
+ Splits inlist into list of lists where each child list len <= size.
+
+ Args:
+ inlist (list): list that we want cut into smaller lists
+ size (int): max length of small lists returned
+
+ Returns:
+ res (list of lists): list of list where each child list <= size length
+ """
+ res = []
+ numlists = ceil(len(inlist)/size)
+ for ind in range(numlists):
+ res.append(inlist[ind*size:(ind+1)*size])
+ return res
+
+[docs]def listify(items):
+ """Returns a list of items. If items is an object it puts it in a list.
+
+ If a list is passed, it returns the list without changing it.
+ If an object is passed, it returns a list with the object in it.
+
+ Args:
+ items (object or list): item or list to convert to a list
+ """
+ if not isinstance(items, list):
+ items = [items]
+ return items
+
+"""This module stores the CadImporter class, which is used to load CAD parts."""
+
+import collections
+import math
+import os
+# needed to prevent dxfgrabber from crashing on import
+os.environ['DXFGRABBER_CYTHON'] = 'OFF'
+import dxfgrabber # needed for dxf files
+import subprocess # needed to run gmsh to make geo files
+
+from . import geometry
+from . import partmodule
+
+[docs]class CadImporter(object):
+ """Makes an object which can import cad parts.
+
+ Args:
+ feamodel (FeaModel): model that we want to import parts into
+ layer (int): layer to import, all entities will be flattened
+ to one plane. -1 mean all layers, any other number is a specific layer
+ swapxy (bool): True rotates the part from x axial to y axial
+ scale (str): string telling the unit conversion scalar
+
+ * Any option of 'fromunit-tounit' using the below units
+ * mm, m, in, ft
+ * Examples 'mm-m' 'm-in' 'ft-mm' 'in-ft'
+ * Default value is '' and does not apply a scale factor
+
+ Attributes:
+ __fea (FeaModel): model
+ __fname (str): file we want to import
+ __layer (int): layer to import
+ __swapxy (bool): If true, swap part from xy to yx orientation
+ """
+
+ def __init__(self, feamodel, fname='', layer=-1, swapxy=False, scale=''):
+ self.__fea = feamodel
+ self.__fname = fname
+ self.__layer = layer
+ self.__swapxy = swapxy
+ self.__scale = scale
+
+[docs] def load(self):
+ """Loads the self.__fname cad file
+
+ Returns:
+ list: list of Part
+ """
+ if self.__fname == '':
+ print('You must pass in a file name to load!')
+ return []
+
+ ext = os.path.splitext(self.__fname)[1]
+ first_pt = None
+ if len(self.__fea.points) > 0:
+ first_pt = self.__fea.points[0]
+ if ext == '.dxf':
+ parts = self.__load_dxf()
+ elif ext in ['.brep', '.brp', '.iges', '.igs', '.step', '.stp']:
+ self.__make_geo()
+ parts = self.__load_geo()
+ last_pt = None
+ if first_pt != None:
+ if len(self.__fea.points) > 2:
+ last_pt = self.__fea.points[-1]
+ if self.__scale != '':
+ # call scale
+ pass
+ return parts
+
+ def __fix_tuple(self, xy_tup):
+ """Adjusts the point to be in the right plane (yx)"""
+ if self.__swapxy:
+ return xy_tup[::-1]
+ return xy_tup
+
+ @staticmethod
+ def __find_make_pt(xy_tup, points_dict):
+ """
+ Returns a point if it exists within geometry.ACC
+ or makes, stores and returns the point if it doesn't exist
+ """
+ point = points_dict.get(xy_tup)
+ if point is not None:
+ return point
+ xy_point = geometry.Point(xy_tup[0], xy_tup[1])
+ for realpoint in points_dict.values():
+ dist = (xy_point - realpoint).length()
+ if dist < geometry.ACC:
+ return realpoint
+ points_dict[xy_tup] = xy_point
+ return xy_point
+
+ def __get_pts_lines(self, lines, arcs):
+ """Returns a set of points, and a list of Lines and Arcs
+
+ Args:
+ lines (dxfgrabber LINE list): dxf lines
+ arcs (dxfgrabber ARC list): dxf arcs
+
+ Returns:
+ list: [list of points, list of Line and Arc]
+ """
+ # store unique points
+ points_dict = {}
+ all_lines = []
+ for ind, line in enumerate(lines):
+ tup = self.__fix_tuple((line.start[0], line.start[1]))
+ start = self.__find_make_pt(tup, points_dict)
+ tup = self.__fix_tuple((line.end[0], line.end[1]))
+ end = self.__find_make_pt(tup, points_dict)
+ line = geometry.Line(start, end)
+ all_lines.append(line)
+ for ind, arc in enumerate(arcs):
+ # dxfgrabber arcs are stored ccw when looking at xy plane
+ # x horizontal
+ # y vertical
+ tup = self.__fix_tuple((arc.center[0], arc.center[1]))
+ center = self.__find_make_pt(tup, points_dict)
+ sign = -1
+ if self.__swapxy:
+ sign = 1
+ startangle = arc.start_angle*sign
+ endangle = arc.end_angle*sign
+ angle = endangle - startangle
+ if arc.end_angle < arc.start_angle:
+ angle = angle + 360*sign
+ """
+ print('---------------------------------------')
+ print('| ARC')
+ print('center: %s' % center)
+ print('startangle: %f' % startangle)
+ print('endangle: %f' % endangle)
+ print('traversed_angle: %f' % angle)
+ """
+ start_vect = geometry.Point(0, arc.radius)
+ if self.__swapxy == False:
+ start_vect = geometry.Point(arc.radius, 0)
+ start_vect.rot_ccw_deg(arc.start_angle*sign)
+ end_vect = geometry.Point(0, arc.radius)
+ if self.__swapxy == False:
+ end_vect = geometry.Point(arc.radius, 0)
+ end_vect.rot_ccw_deg(arc.end_angle*sign)
+ start = center + start_vect
+ start_tup = (start.x, start.y)
+ end = center + end_vect
+ end_tup = (end.x, end.y)
+ start = self.__find_make_pt(start_tup, points_dict)
+ end = self.__find_make_pt(end_tup, points_dict)
+ rvect = start - center
+ if abs(angle) <= 90:
+ arc = geometry.Arc(start, end, center)
+ all_lines.append(arc)
+ print('1 arc made')
+ continue
+ #print(' %s' % arc)
+ pieces = math.ceil(abs(angle)/90)
+ print('%i arcs being made' % pieces)
+ points = [start, end]
+ # 2 pieces need 3 points, we have start + end already --> 1 pt
+ inserts = pieces + 1 - 2
+ piece_ang = angle/pieces
+ #print('piece_ang = %f' % piece_ang)
+ while inserts > 0:
+ rvect.rot_ccw_deg(piece_ang)
+ point = center + rvect
+ tup = (point.x, point.y)
+ point = self.__find_make_pt(tup, points_dict)
+ points.insert(-1, point)
+ inserts = inserts - 1
+ for ind in range(len(points)-1):
+ #print(' %s' % arc)
+ arc = geometry.Arc(points[ind], points[ind+1], center)
+ all_lines.append(arc)
+ for line in all_lines:
+ line.save_to_points()
+ return [list(points_dict.values()), all_lines]
+
+ def __make_geo(self):
+ """Makes a gmsh geo file given a step, iges, or brep input"""
+ # gmsh freecad_part.iges -o out_iges.geo -0
+ fname_list = self.__fname.split('.')
+ geo_file = fname_list[0]+'.geo'
+ runstr = "%s %s -o %s -0" % (environment.GMSH, self.__fname, geo_file)
+ print(runstr)
+ subprocess.call(runstr, shell=True)
+ print('Wrote file: %s' % geo_file)
+
+ def __load_geo(self):
+ """Loads in a gmsh geo file and returns a list of parts
+
+ Returns:
+ list: list of Part
+ """
+ pass
+ # process any splines? and turn them into arcs
+ # http://www.mathopenref.com/constcirclecenter.html
+ # find max dist between points
+ # double it
+ # select two segments
+ # draw normal lines
+ # find intersections, that is the center
+
+ @staticmethod
+ def __dangling_points(all_points):
+ return [point for point in all_points
+ if len(point.lines) == 1 and not point.arc_center]
+
+ def __load_dxf(self):
+ """Loads in a dxf file and returns a list of parts
+
+ Returns:
+ list: list of Part
+ """
+ print('Loading file: %s' % self.__fname)
+ dwg = dxfgrabber.readfile(self.__fname)
+ lines = [item for item in dwg.entities if item.dxftype == 'LINE']
+ arcs = [item for item in dwg.entities if item.dxftype == 'ARC']
+ if self.__layer > -1:
+ lines = [item for item in lines if item.layer == self.__layer]
+ arcs = [item for item in arcs if item.layer == self.__layer]
+ print('File read.')
+ print('Loaded %i lines' % len(lines))
+ print('Loaded %i arcs' % len(arcs))
+ print('Loaded %i line segments, lines or arcs' %
+ (len(lines)+len(arcs)))
+ # get all points and Line and Arc using pycalculix entities
+ print('Converting to pycalculix lines arcs and points ...')
+ all_points, all_lines = self.__get_pts_lines(lines, arcs)
+ print('Loaded %i line segments, lines or arcs' % len(all_lines))
+ print('Loaded %i points' % len(all_points))
+ # for point in all_points:
+ # print('%s %s' % (point, point.lines))
+ # for line in all_lines:
+ # print('%s %s' % (line, line.points))
+
+ # remove all lines that are not part of areas
+ dangling_points = self.__dangling_points(all_points)
+ pruned_geometry = bool(dangling_points)
+ while dangling_points:
+ for point in dangling_points:
+ all_points.remove(point)
+ print('Removed point= %s' % point)
+ dangling_line = list(point.lines)[0]
+ point.unset_line(dangling_line)
+ if dangling_line in all_lines:
+ all_lines.remove(dangling_line)
+ print('Removed line= %s' % dangling_line)
+ dangling_points = self.__dangling_points(all_points)
+ if pruned_geometry:
+ print('Remaining line segments: %i' % len(all_lines))
+ print('Remaining points: %i' % len(all_points))
+
+ # make line all_loops now
+ all_loops = []
+ line = all_lines[0]
+ this_loop = geometry.LineLoop()
+ while len(all_lines) > 0:
+ this_loop.append(line)
+ all_lines.remove(line)
+ if this_loop.closed == True:
+ all_loops.append(this_loop)
+ this_loop = geometry.LineLoop()
+ if all_lines:
+ line = all_lines[0]
+ continue
+ point = line.pt(1)
+ other_lines = point.lines - set([line])
+ if len(other_lines) > 1:
+ # note: one could exclude connected segment nodes
+ # make disconnected line all_loops, then have another
+ # loop to connect thos disconnected line all_loops
+ print('One point was connected to > 2 lines.')
+ print('Only import simple part all_loops, or surfaces.')
+ raise Exception('Import geometry is too complex')
+ next_line = list(other_lines)[0]
+ if line.pt(1) != next_line.pt(0):
+ next_line.reverse()
+ line = next_line
+
+ # find exterior loops
+ exterior_loops = []
+ for ind, loop in enumerate(all_loops):
+ other_loops = all_loops[ind+1:]
+ other_loops.extend(exterior_loops)
+ is_exterior = True
+ for other_loop in other_loops:
+ if loop.inside(other_loop):
+ is_exterior = False
+ break
+ if is_exterior:
+ # exterior must be clockwise
+ if loop.ccw:
+ loop.reverse()
+ exterior_loops.append(loop)
+ # remove the found part exterior loops from all_loops
+ for exterior_loop in exterior_loops:
+ all_loops.remove(exterior_loop)
+ # each part in parts is a list of line all_loops
+ # [exterior, hole1, hole2]
+ parts = [[exterior_loop] for exterior_loop in exterior_loops]
+ # now place the child hole loops after the part exterior loop
+ for part_loops in parts:
+ exterior_loop = part_loops[0]
+ # find child holes
+ for hole_loop in all_loops:
+ if hole_loop.inside(exterior_loop):
+ hole_loop.hole = True
+ # holes must be ccw
+ if not hole_loop.ccw:
+ hole_loop.reverse()
+ part_loops.append(hole_loop)
+ # remove child holes from loop list
+ for hole_loop in part_loops[1:]:
+ all_loops.remove(hole_loop)
+
+ # make parts
+ parts_list = []
+ for part_loops in parts:
+ this_part = partmodule.Part(self.__fea)
+ for ind, loop in enumerate(part_loops):
+ is_hole = loop.hole
+ start = loop[0].pt(0)
+ this_part.goto(start.x, start.y, is_hole)
+ for item in loop:
+ if isinstance(item, geometry.Line):
+ end = item.pt(1)
+ this_part.draw_line_to(end.x, end.y)
+ elif isinstance(item, geometry.Arc):
+ end = item.pt(1)
+ center = item.actr
+ this_part.draw_arc(end.x, end.y, center.x, center.y)
+ parts_list.append(this_part)
+ print('Parts created: %i' % len(parts_list))
+ return parts_list
+
+"""This module defines the Component class, which is used to make components.
+
+Components are collections that loads are applied to.
+"""
+from . import base_classes
+from . import geometry
+
+[docs]class Component(base_classes.Idobj):
+ """Makes a component to store loads, boundary conditions, matl, thickness.
+
+ A component must store a list of items where the type is identical for
+ all items.
+
+ Args:
+ item_list (list): list of items to store
+ ctype (str): type of component. Options:
+
+ - 'nodes': all element nodes (corner and mid-side)
+ - 'n1': element corner nodes
+ - 'faces': faces
+ - 'elements': elements
+
+ cname (str): component name
+ """
+
+ def __init__(self, item_list, ctype, cname=''):
+ self.items = item_list
+ self.ctype = ctype
+ self.name = cname+'_'+ctype
+ base_classes.Idobj.__init__(self)
+
+ def __eq__(self, other):
+ """ Returns bool telling if self == other."""
+ if isinstance(other, Component):
+ if self.items == other.items and self.ctype == other.ctype:
+ return True
+ return False
+ else:
+ return False
+
+ def __hash__(self):
+ """Returns a hash of self instance."""
+ return self.id
+
+[docs] def set_name(self, name):
+ """Set the component name to the passed name.
+
+ Args:
+ name (str): passed name, component name is changed to this
+ """
+ self.name = name
+
+[docs] def get_name(self):
+ """Returns the component name created with its id number."""
+ return 'C'+str(self.id)
+
+[docs] def get_children(self):
+ """Returns a list the child ctype items of the current component."""
+ res = []
+ for i in self.items:
+ children = getattr(i, self.ctype)
+ if isinstance(children, list):
+ res += children
+ else:
+ res += [children]
+ return res
+
+[docs] def ccx(self):
+ """Writes a component for Calculix ccx solver.
+
+ Returns a list strings, where each item is a line to write
+ to a Calculix .inp text file.
+
+ Note:
+ Only node and element components should use this.
+ """
+ res = []
+ items_per_line = 6
+ firstline = ''
+ if self.ctype == 'nodes':
+ # node component, displacement
+ firstline = '*NSET,NSET='+self.name
+ elif self.ctype == 'elements':
+ # element component, materials
+ firstline = '*ELSET,ELSET='+self.name
+ elif self.ctype == 'faces':
+ # face component, contact surface
+ firstline = '*SURFACE,NAME=%s,TYPE=ELEMENT' % self.name
+ res.append(firstline)
+ items = self.get_children()
+ if self.ctype == 'faces':
+ for face in items:
+ enum = face.element.id
+ fnum = face.id
+ res.append("%i,S%i" % (enum, fnum))
+ elif self.ctype in ['nodes', 'elements']:
+ grouped_items = base_classes.chunk_list(items, items_per_line)
+ for group in grouped_items:
+ item_ids = [str(x.id) for x in group]
+ line = ', '.join(item_ids)
+ if group != grouped_items[-1]:
+ line += ','
+ res.append(line)
+ return res
+
+[docs] def write_cgx(self):
+ """Writes a component for Calculix cgx prerocessor.
+
+ Returns a list strings, where each item is a line to write
+ to a Calculix .fbd text file.
+
+ Note:
+ Only line and point node components should use this.
+ """
+ res = []
+
+ # set the type of component
+ parent = 'l'
+ if (isinstance(self.items[0], geometry.Line) or
+ isinstance(self.items[0], geometry.Arc)):
+ parent = 'l'
+ elif isinstance(self.items[0], geometry.Point):
+ parent = 'p'
+
+ # pull the list of entities
+ alist = ' '.join([a.get_name() for a in self.items])
+
+ if self.ctype == 'n1':
+ # write component of line element first order nodes
+ res.append('seta '+self.name+' '+parent+' '+alist)
+ elif self.ctype == 'nodes':
+ res.append('seta '+self.name+' '+parent+' '+alist)
+ res.append('comp '+self.name+' do')
+ # write component of all line element nodes
+ elif self.ctype == 'f':
+ # write component of line element faces
+ res.append('seta '+self.name+' '+parent+' '+alist)
+ res.append('comp '+self.name+' do')
+ res.append('comp '+self.name+' do')
+ return res
+
+[docs] def write_gmsh(self):
+ """Writes a component for gmsh mesher.
+
+ Returns a list strings, where each item is a line to write
+ to a Calculix .geo text file.
+
+ Note:
+ Only line and point node components should use this.
+ """
+ res = []
+
+ # set the type of component
+ parent = 'l'
+ if (isinstance(self.items[0], geometry.Line) or
+ isinstance(self.items[0], geometry.Arc)):
+ parent = 'l'
+ elif isinstance(self.items[0], geometry.Point):
+ parent = 'p'
+
+ # pull the list of entities
+ alist = ','.join([str(a.id) for a in self.items])
+
+ line = ''
+ if parent == 'l':
+ line = "Physical Line('"+self.name+"') = {" + alist + '};'
+ elif parent == 'p':
+ line = "Physical Point('"+self.name+"') = {" + alist + '};'
+ res.append(line)
+ return res
+
+"""This module stores connector classes, like Contact.
+"""
+from . import base_classes
+
+[docs]class SurfaceInteraction(base_classes.Idobj):
+ """Makes a surface interaction object.
+
+ Args:
+ int_type (str): interaction type
+
+ * 'EXPONENTIAL'
+ * 'LINEAR'
+
+ *args: following arguments
+
+ * int_type = 'EXPONENTIAL'
+ * c0, then p0 must be passed
+
+ * int_type = 'LINEAR'
+ * k must be passed
+
+ """
+
+ def __init__(self, int_type, *args):
+ self.int_type = int_type
+ self.c0 = None
+ self.p0 = None
+ self.k = None
+ if self.int_type == 'EXPONENTIAL':
+ self.c0 = args[0]
+ self.p0 = args[1]
+ elif self.int_type == 'LINEAR':
+ self.k = args[0]
+ base_classes.Idobj.__init__(self)
+
+ @property
+ def name(self):
+ """SurfaceInteraction name."""
+ return 'SI%i' % self.id
+
+[docs] def ccx(self):
+ """Writes the surface interaction for the ccx solver."""
+ res = ["*SURFACE INTERACTION,NAME=%s" % self.name]
+ res.append("*SURFACE BEHAVIOR,PRESSURE-OVERCLOSURE=%s" % self.int_type)
+ if self.int_type == 'EXPONENTIAL':
+ res.append("%f,%f" % (self.c0, self.p0))
+ elif self.int_type == 'LINEAR':
+ res.append("%e" % self.k)
+ #res.append("%e,1." % self.k)
+ return res
+
+[docs]class Contact(base_classes.Idobj):
+ """Makes a contact which will be between lines which have faces on them.
+
+ Args:
+ master_comp (Component): component of master element faces
+ slave_comp (Component): component of slave element faces
+ surf_int (SurfaceInteraction): object which stores closure behavior
+ surf_to_surf (bool): if True surface to surface is used, if False node
+ to surface is used
+ """
+
+ def __init__(self, master_comp, slave_comp, surf_int, surf_to_surf=True):
+ self.master_comp = master_comp
+ self.slave_comp = slave_comp
+ self.surf_int = surf_int
+ self.surf_to_surf = surf_to_surf
+ base_classes.Idobj.__init__(self)
+
+ @property
+ def name(self):
+ """Contact name."""
+ return 'CONT%i' % self.id
+
+[docs] def ccx(self):
+ """Writes the contact pair for the ccx solver."""
+ res = []
+ end = ''
+ int_name = self.surf_int.name
+ m_name = self.master_comp.name
+ s_name = self.slave_comp.name
+ if self.surf_to_surf:
+ end = ',TYPE=SURFACE TO SURFACE'
+ line = "*CONTACT PAIR,INTERACTION=%s%s" % (int_name, end)
+ res.append(line)
+ line = "%s,%s" % (s_name, m_name)
+ res.append(line)
+ return res
+
+"""This module sets the dpi and the paths to gmsh, ccx, and cgx.
+
+Attributes:
+ DPI (None or float): if high dpi windows 8 monitor value is set, otherwise
+ value is None.
+ CCX (str): path to Calculix ccx, the sovler.
+ CGX (str): path to Calculix cgx, the preprocessor/postprocessor/mesher.
+ GMSH (str): path to the Gmsh mesher.
+"""
+
+import sys # needed to check if 32 or 64 bit interpreter
+import os # used to delete files written by cgx
+import platform # need this to check for win and do the below dpi fix
+import ctypes # needed to see if it's 32-bit and to get correct win version
+
+# http://stackoverflow.com/questions/19128219/detect-windows-8-1-in-python/22325767#22325767
+[docs]class OSVERSIONINFOEXW(ctypes.Structure):
+ """Returns object w/ attributes that will identify window os version"""
+ _fields_ = [('dwOSVersionInfoSize', ctypes.c_ulong),
+ ('dwMajorVersion', ctypes.c_ulong),
+ ('dwMinorVersion', ctypes.c_ulong),
+ ('dwBuildNumber', ctypes.c_ulong),
+ ('dwPlatformId', ctypes.c_ulong),
+ ('szCSDVersion', ctypes.c_wchar*128),
+ ('wServicePackMajor', ctypes.c_ushort),
+ ('wServicePackMinor', ctypes.c_ushort),
+ ('wSuiteMask', ctypes.c_ushort),
+ ('wProductType', ctypes.c_byte),
+ ('wReserved', ctypes.c_byte)]
+
+[docs]def get_version():
+ """Get's the OS version. Returns a float of OS_MAJOR.OS_MINOR
+ """
+ os_version = OSVERSIONINFOEXW()
+ os_version.dwOSVersionInfoSize = ctypes.sizeof(os_version)
+ retcode = ctypes.windll.Ntdll.RtlGetVersion(ctypes.byref(os_version))
+ if retcode != 0:
+ print("Failed to get OS version")
+ return -1.0
+
+ major = str(os_version.dwMajorVersion)
+ minor = str(os_version.dwMinorVersion)
+ return float(major+'.'+minor)
+
+[docs]def get_dpi():
+ """Returns an int of current DPI for windows or None."""
+ if platform.system() == 'Windows':
+ version = get_version()
+ if version >= 6.2:
+ # Windows 8.0 is 6.2.9200
+ # Only run on windows 8 or higher to check high dpi monitors
+
+ user32 = ctypes.windll.user32
+ w_curr = user32.GetSystemMetrics(0)
+ user32.SetProcessDPIAware()
+ w_phys = user32.GetSystemMetrics(0)
+ curr_dpi = round(w_phys*96/w_curr, 0)
+
+ from pylab import rcParams
+ rcParams['figure.dpi'] = curr_dpi
+ return curr_dpi
+ else:
+ return None
+
+[docs]def get_paths():
+ """Returns a list of paths to: [ccx, cgx, gmsh]."""
+ ccx = 'ccx'
+ cgx = 'cgx'
+ gmsh = 'gmsh'
+ return [ccx, cgx, gmsh]
+
+DPI = get_dpi()
+CCX, CGX, GMSH = get_paths()
+
+""" This module stores the FeaModel class.
+
+This class is the highest level object in a pycalculix program.
+It stores all parts, loads, constraints, mesh, problem, and results_file
+objects.
+"""
+
+import matplotlib.pyplot as plt
+from matplotlib.collections import PatchCollection # element plotting
+import matplotlib.colors as colors
+import matplotlib.cm as cmx
+import subprocess # used to launch meshers cgx and gmsh
+import os # used to delete files written by cgx
+from numpy.core.function_base import linspace # need to make contours
+
+from . import environment
+from . import base_classes
+from . import geometry
+from . import material
+from . import components
+from . import connectors
+from . import loads
+from . import mesh
+from . import partmodule
+from . import problem
+from . import selector
+
+# element colors, 0-1, 0=black, 1=whate
+ECOLOR = '.4'
+FCOLOR = '.9'
+CMAP = 'jet' # for results
+GEOM_CMAP = 'Pastel1' # color map for parts or areas
+
+[docs]class FeaModel(object):
+ """Makes a FeaModel instance.
+
+ Parts, area, lines, arcs, points, nodes, elements, faces, models,
+ components, and loads are stored in this object.
+
+ Args:
+ model_name (str): project name for this FeaModel, this is a file prefix
+ ccx (None or str): path to Calculix ccx solver, pass this when you want
+ to overwite the default program location.
+ None means the default envionment.CCX is used.
+ cgx (None or str): path to Calculix cgx mesher, pass this when you want
+ to overwite the default program location.
+ None means the default envionment.CGX is used.
+ gmsh (None or str): path to gmsh mesher, pass this when you want
+ to overwite the default program location.
+ None means the default envionment.GMSH is used.
+
+ Attributes:
+ fname (str): FeaModel project file name prefix
+ points (Itemlist): list of all Point
+ lines (Itemlist): list of all Line and Arc
+ signlines (Itemlist): list of all SignLine and SignArc
+ lineloops (Itemlist): list of all LineLoop, each contains SignLine SignArc
+ areas (Itemlist):list of all Area, each contains LineLoop(s)
+ parts (Itemlist): list of all Part
+ matls (Itemlist): list of all materials
+ components (Itemlist): list of all components
+ loads (dict): a dictionary of loads
+
+ | Key (float): the load time point
+ | loads[time] = list of loads for that time step
+ | See method set_time to change the current time.
+ | Time = 0.0 stores constant loads, such as:
+ | material, thickness
+
+ contacts (Itemlist): list of contacts
+ surfints (Itemlist): list of surface interactions
+ problems (Itemlist): list of problems
+ nodes (Meshlist): list of all mesh nodes
+ eshape (str): element shape
+
+ - 'quad': quadrilateral elements (Default)
+ - 'tri': triangle elements
+ eorder (int): element order, 1 or 2
+
+ - 1: elements meshed with corner nodes only
+ - 2: (Default) elements meshed with corner and midside nodes
+ elements (Meshlist): list of all mesh elements
+ faces (list): list of all element faces, includes non-exterior ones
+ view (Selector): currently selected items
+
+ Attributes:
+ parts: list of selected parts
+ areas: list of selected areas
+ lines: list of selected signed lines or arcs
+ points: list of selected points
+ elements: list of selected elements
+ faces: list of selected faces
+ nodes: list of selected nodes
+ time (float): current model time value, defaults to 1.0
+ units (dict): a dict to store units for different fields
+
+ Keys:
+ - 'displ': displacement or location
+ - 'force': force
+ - 'stress': stress
+ - 'temp': temperature
+ - 'density': density (mass/volume)
+ - 'time': time
+
+ Returns:
+ Text string describing the units for the given key field.
+
+ See the set_units method.
+
+ For example when units have been set to metric, the below
+ values are returned.
+
+ - 'dist' --> 'm'
+ - 'density' --> 'kg/(m^3)'
+ - 'time' --> 's'
+ """
+ def __init__(self, model_name, ccx=None, cgx=None, gmsh=None):
+ self.fname = model_name
+ self.points = base_classes.Itemlist()
+ self.lines = base_classes.Itemlist()
+ self.signlines = base_classes.Itemlist()
+ self.lineloops = base_classes.Itemlist()
+ self.areas = base_classes.Itemlist()
+ self.parts = base_classes.Itemlist()
+ self.matls = base_classes.Itemlist()
+ self.components = base_classes.Itemlist()
+ self.loads = {} # store loads by time
+ self.contacts = base_classes.Itemlist()
+ self.surfints = base_classes.Itemlist()
+ self.problems = base_classes.Itemlist()
+ self.nodes = base_classes.Meshlist()
+ self.elements = base_classes.Meshlist()
+ self.faces = []
+ self.view = selector.Selector(self)
+ self.time = 1.0 # 0 is model set-up, 1 is first step, etc
+ self.eshape = 'quad'
+ self.eorder = 2
+ self.units = {}
+
+ # fix the paths to the needed programs, ccx, cgx, and gmsh
+ if ccx != None:
+ environment.CCX = ccx
+ if cgx != None:
+ environment.CGX = cgx
+ if gmsh != None:
+ environment.GMSH = gmsh
+
+[docs] def set_ediv(self, items, ediv):
+ """Sets the number of elements on the passed line.
+
+ Args:
+ items (str or SignLine or SignArc or list): lines to set ediv on
+
+ - str: 'L0'
+ - list of str ['L0', 'L1']
+ - list of SignLine or SignArc part.bottom or part.hole[-1]
+
+ ediv (int): number of elements to mesh on the line
+ """
+ items = self.get_items(items)
+ for line in items:
+ line.set_ediv(ediv)
+
+[docs] def set_esize(self, items, esize):
+ """Sets the element size on the passed line.
+
+ Args:
+ items (str or SignLine or SignArc or list): lines or points to set esize on
+
+ - str: 'L0'
+ - list of str ['L0', 'L1', 'P3']
+ - list of SignLine or SignArc part.bottom or part.hole[-1]
+
+ esize (float): size of the mesh elements on the line
+ """
+ items = self.get_items(items)
+ for item in items:
+ item.set_esize(esize)
+
+[docs] def set_units(self, dist_unit='m', cfswitch=False):
+ """Sets the units that will be displayed when plotting.
+
+ Picks a unit set based on the passed distance unit.
+ That unit set is printed to the console when set.
+ Defaults to MKS units (meter-kilogram-second)
+
+ ======== ===== ====== =========== =============== ====
+ Distance Force Stress Temperature Density Time
+ ======== ===== ====== =========== =============== ====
+ 'm' 'N' 'Pa' 'K' 'kg/(m^3)' 's'
+ 'mm' 'N' 'MPa' 'K' 'tonne/(mm^3)' 's'
+ 'in' 'lbf' 'psi' 'R' 'slinch/(in^3)' 's'
+ 'ft' 'lbf' 'psf' 'R' 'slug/(ft^3)' 's'
+ ======== ===== ====== =========== =============== ====
+
+ See get_units method or returning text strings based on unit types.
+
+ Args:
+ dist_unit (str): string of distance unit. Options:
+
+ - 'm': meter
+ - 'mm': milimeter
+ - 'in': inch
+ - 'ft': foot
+
+ cfswitch (bool): Celsius/Fahrenheit temperature switch.
+ Default is False.
+
+ If True, this switches temperature from K-->C or from R-->F
+ Default keeps units as shown in above table.
+ """
+ keys = ['dist', 'force', 'stress', 'temp', 'density', 'time']
+
+ m_newton = ['m', 'N', 'Pa', 'K', 'kg/(m^3)', 's']
+ mm_newton = ['mm', 'N', 'MPa', 'K', 'tonne/(mm^3)', 's']
+ in_lbf = ['in', 'lbf', 'psi', 'R', 'slinch/(in^3)', 's']
+ ft_lbf = ['ft', 'lbf', 'psf', 'R', 'slug/(ft^3)', 's']
+ unit_systems = [m_newton, mm_newton, in_lbf, ft_lbf]
+ vals = [usys for usys in unit_systems if usys[0] == dist_unit][0]
+
+ if cfswitch:
+ # switches K-->C and R-->F
+ newunit = {'K':'C', 'R':'F'}
+ vals['temp'] = newunit[vals['temp']]
+
+ # set values
+ adict = dict(zip(keys, vals))
+ adict['displ'] = adict['dist']
+ print('================================================')
+ print('Units have been set to %s_%s' % (adict['dist'], adict['force']))
+ for key in adict:
+ print('For %s use %s' % (key, adict[key]))
+ print('================================================')
+ self.units = adict
+
+[docs] def get_units(self, *args):
+ """Returns units for the passed arguments. Accepts and returns a list
+ of strings.
+
+ Options for inputs:
+ - 'dist'
+ - 'displ' (same units as 'dist')
+ - 'force'
+ - 'stress'
+ - 'temp'
+ - 'density'
+ - 'time'
+ """
+ res = []
+ for arg in args:
+ mystr = ''
+ if arg in self.units:
+ # check if the requested item is in the units dict
+ mystr = ' ('+ self.units[arg] + ')'
+ elif arg in base_classes.FIELDTYPE:
+ # return units on results dict item, ex: stress, strain etc
+ ftype = base_classes.FIELDTYPE[arg]
+ if ftype in self.units:
+ mystr = ' ('+ self.units[ftype] + ')'
+ res.append(mystr)
+ return res
+
+[docs] def set_time(self, time):
+ """Sets the time in the FeaModel (preprocessor).
+
+ This time is used when setting loads.
+
+ Args:
+ time (float): the time to set
+ """
+ self.time = time
+
+[docs] def get_item(self, item):
+ """Returns an item given a string identifying the item.
+
+ Args:
+ item (str): 'A0' or 'P0' or 'L0' etc.
+ """
+ if item[0] == 'P':
+ # get point
+ items = self.points
+ num = int(item[1:])
+ res = [a for a in items if a.id == num]
+ return res[0]
+ elif item[0] == 'L' or item[1] == 'L':
+ # get line
+ items = self.signlines
+ num = int(item[1:])
+ items = [a for a in items if a.get_name() == item]
+ return items[0]
+ elif item[0] == 'A':
+ # get area
+ items = self.areas
+ num = int(item[1:])
+ items = [a for a in items if a.id == num]
+ return items[0]
+ elif item[0] == 'E':
+ # get element
+ items = self.elements
+ items = [a for a in items if a.get_name() == item]
+ return items[0]
+ elif item[0] == 'N':
+ # get node
+ items = self.nodes
+ items = [a for a in items if a.get_name() == item]
+ return items[0]
+ else:
+ print('Unknown item! Please pass the name of a point, line or area!')
+
+[docs] def get_items(self, items):
+ """Returns a list of correctly typed items.
+
+ Input can be a single string or item, or a list of strings identifying
+ items.
+ """
+ res_items = base_classes.listify(items)
+ for ind, item in enumerate(res_items):
+ if isinstance(item, str):
+ res_items[ind] = self.get_item(item)
+ return res_items
+
+[docs] def make_matl(self, name):
+ """Makes and returns a new material.
+
+ Args:
+ name (str): the material's name
+ """
+ mat = material.Material(name)
+ self.matls.append(mat)
+ return mat
+
+[docs] def make_part(self):
+ """Makes and returns a new part."""
+ #p = part.Part(self)
+ #self.parts.append(p)
+ return partmodule.Part(self)
+
+[docs] def make_problem(self, problem_type='struct', parts='all'):
+ """Makes and returns a new problem, which can be solved.
+
+ Args:
+ problem_type (str): problem type, options:
+ 'struct': structural
+ parts (str Part or list of Part): Parts the model will analyze.
+
+ Options:
+ - 'all': add all parts. This is the default.
+ - Part: add the single part
+ - list of Part: add these parts
+ """
+ if parts == 'all':
+ parts = self.parts
+ prob = problem.Problem(self, problem_type, parts)
+ return prob
+
+[docs] def print_summary(self):
+ """Prints a summary of the items in the database.
+ """
+ items = ['parts', 'areas', 'signlines', 'points', 'elements', 'faces',
+ 'nodes']
+ names = ['parts', 'areas', 'lines', 'points', 'elements', 'faces',
+ 'nodes']
+
+ spacer = '----------------------------------'
+ print(spacer)
+ print("Items in the database:")
+ for (item_str, name) in zip(items, names):
+ items = getattr(self, item_str)
+ print(' %s: %i' % (name, len(items)))
+ print(spacer)
+
+[docs] def plot_nodes(self, fname='', display=True, title='Nodes', nnum=False):
+ """Plots the selected nodes.
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown
+ title (str): the plot's title
+ nnum (bool): if True node numbers are plotted
+ """
+ nodes = self.view.nodes
+ if len(nodes) > 0:
+ # plotting elements
+ fig = plt.figure()
+ axis = fig.add_subplot(111)
+
+ # plot nodes, this is quicker than individual plotting
+ axs = [node.y for node in nodes]
+ rads = [node.x for node in nodes]
+ axis.scatter(axs, rads, s=7, color='black')
+ if nnum:
+ for node in nodes:
+ node.label(axis)
+
+ # set units
+ [d_unit] = self.get_units('dist')
+ plt.title(title)
+ plt.xlabel('axial, y'+d_unit)
+ plt.ylabel('radial, x'+d_unit)
+ plt.axis('scaled')
+
+ # extract max and min for plot window
+ radials = [n.x for n in nodes]
+ axials = [n.y for n in nodes]
+
+ # finish pot
+ base_classes.plot_set_bounds(plt, axials, radials)
+ base_classes.plot_finish(plt, fname, display)
+
+ else:
+ # no elements exist or no elemnts are selected
+ res = ''
+ if len(self.nodes) == 0:
+ res = 'No nodes exist! Try meshing your parts!'
+ else:
+ res = 'No nodes are selected! Select some!'
+ print(res)
+
+[docs] def plot_elements(self, fname='', display=True, title='Elements',
+ enum=False, nshow=False, nnum=False):
+ """Plots the selected elements.
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown
+ title (str): the plot's title
+ enum (bool): if True element numbers are plotted
+ nshow (bool): True=plot nodes, False=don't plot them
+ nnum (bool): if True node numbers are plotted
+ """
+ nodes = self.view.nodes
+ elements = self.view.elements
+ if len(elements) > 0:
+ # plotting elements
+ fig = plt.figure()
+ axis = fig.add_subplot(111)
+ polys = []
+ for element in elements:
+ poly = element.get_poly()
+ polys.append(poly)
+ coll = PatchCollection(polys, facecolors=FCOLOR, edgecolors=ECOLOR)
+ axis.add_collection(coll)
+
+ # plot element numbers
+ if enum:
+ for element in elements:
+ element.label(axis)
+
+ # plot nodes, this is quicker than individual plotting
+ if nshow:
+ axs = [node.y for node in nodes]
+ rads = [node.x for node in nodes]
+ axis.scatter(axs, rads, s=7, color='black')
+ if nnum:
+ for node in nodes:
+ node.label(axis)
+
+ # set units
+ [d_unit] = self.get_units('dist')
+ plt.title(title)
+ plt.xlabel('axial, y'+d_unit)
+ plt.ylabel('radial, x'+d_unit)
+ plt.axis('scaled')
+
+ # extract nodes for elements if no nodes are selected
+ if len(nodes) == 0:
+ tmp = set()
+ for element in elements:
+ tmp.update(element.nodes)
+ nodes = list(tmp)
+ # extract max and min for plot window
+ radials = [n.x for n in nodes]
+ axials = [n.y for n in nodes]
+
+ # finish pot
+ base_classes.plot_set_bounds(plt, axials, radials)
+ base_classes.plot_finish(plt, fname, display)
+
+ else:
+ # no elements exist or no elemnts are selected
+ res = ''
+ if len(self.elements) == 0:
+ res = 'No elements exist! Try meshing your parts!'
+ else:
+ res = 'No elements are selected! Select some!'
+ print(res)
+
+[docs] def plot_pressures(self, fname='', display=True):
+ """Plots the load step pressures on the passed part(s) or selected
+ elements.
+
+ This is an element plot, with arrows showing pressure magnitude and
+ directions.
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown. Default = True
+ """
+ elements = self.view.elements
+ faces = self.view.faces
+ nodes = self.view.nodes
+
+ # plot all elements and store length, length determins min pressure arrow
+ if len(elements) > 0:
+ # plotting elements
+ fig = plt.figure()
+ axis = fig.add_subplot(111, aspect='equal')
+
+ # plot polys and calculate avg face length
+ polys = []
+ face_len = []
+ for element in elements:
+ face_len.append(element.face[1].length())
+ poly = element.get_poly()
+ polys.append(poly)
+ coll = PatchCollection(polys, edgecolors=ECOLOR, facecolors=FCOLOR)
+ axis.add_collection(coll)
+ face_len = sum(face_len)/len(face_len)
+
+ # store pressures we'll want to plot: list of [face, pval]
+ plist = []
+ for load in self.loads[self.time]:
+ if load.ltype in ['press', 'press_fluid']:
+ loadlist = load.get_list()
+ for [face, pval] in loadlist:
+ if face in faces:
+ plist.append([face, pval])
+ pressures = [abs(pval) for [face, pval] in plist]
+
+ # set max and min pressure bounds
+ pmin = min(pressures)
+ pmax = max(pressures)
+ # extract nodes for elements if no nodes are selected
+ if len(nodes) == 0:
+ tmp = set()
+ for element in elements:
+ tmp.update(element.nodes)
+ nodes = list(tmp)
+ radials = [p.x for p in nodes]
+ axials = [p.y for p in nodes]
+ # arrow length = arrow_min + (pval-pmin)*mult
+ arrow_min = face_len
+ mult = 0
+ if pmax != pmin:
+ adelta = max(axials) - min(axials)
+ rdelta = max(radials) - min(radials)
+ delta = max(adelta, rdelta)
+ dist = delta*0.2
+ mult = dist/(pmax - pmin)
+
+ # make tick list for later plot, and color map
+ cbar_val = None
+ tick_list = []
+ cmap = None
+ if pmax != pmin:
+ # we have a range of values we're plotting
+ tick_list = linspace(pmin, pmax, 8)
+ cmap = plt.get_cmap(CMAP)
+ else:
+ cbar_val = pmin
+ pmax = pmin + 1.0
+ pmin = pmin - 1.0
+ tick_list = [pmin, cbar_val, pmax] # default 3 values to plot one solid color
+ cmap = colors.ListedColormap(['b', 'b']) # default to plot one val
+
+ # set color contours for arrows
+ cnorm = colors.Normalize(vmin=pmin, vmax=pmax)
+ scalarmap = cmx.ScalarMappable(norm=cnorm, cmap=cmap)
+ scalarmap.set_array([])
+
+ # make arrows, store axials + radials
+ for [face, pval] in plist:
+ [face_point, unit] = face.get_mnorm()
+ arrow_length = arrow_min + (abs(pval) - pmin)*mult
+ other_point = face_point + unit*arrow_length
+ radials.append(other_point.x)
+ axials.append(other_point.y)
+
+ # compression
+ tail = other_point
+ delta = face_point - other_point
+ if pval < 0:
+ # tension
+ tail = face_point
+ delta = other_point - face_point
+
+ headwidth = face_len*0.2
+ headlength = face_len*0.3
+ colorval = scalarmap.to_rgba(abs(pval))
+ plt.arrow(tail.y, tail.x, delta.y, delta.x,
+ color=colorval,
+ head_width=headwidth, head_length=headlength,
+ length_includes_head=True)
+
+ # set the horizontal and vertical axes
+ base_classes.plot_set_bounds(plt, axials, radials)
+
+ # set units and titles
+ [d_unit, p_unit, t_unit] = self.get_units('dist', 'stress', 'time')
+ tstr = 'Pressures %s\nTime=%f%s' % (p_unit, self.time, t_unit)
+ plt.title(tstr)
+ plt.xlabel('axial, y'+d_unit)
+ plt.ylabel('radial, x'+d_unit)
+
+ # set the colorbar
+ cbar = plt.colorbar(scalarmap, orientation='vertical',
+ ticks=tick_list)
+ if cbar_val != None:
+ cbar.ax.set_yticklabels(['', str(cbar_val), ''])
+ base_classes.plot_finish(plt, fname, display)
+
+ else:
+ res = ''
+ if len(self.elements) == 0:
+ res = 'No elements exist! Try meshing your parts!'
+ else:
+ res = 'No elements are selected! Select some!'
+ print(res)
+
+[docs] def plot_constraints(self, fname='', display=True):
+ """Plots the constraints on the passed part(s) or selected items.
+
+ This is an element and node plot, with arrows showing displacement
+ constraints.
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown. Default = True
+ """
+ elements = self.view.elements
+ nodes = self.view.nodes
+
+ # plot all elements
+ if len(elements) > 0:
+ # plotting elements
+ fig = plt.figure()
+ axis = fig.add_subplot(111, aspect='equal')
+
+ # plot polys and calculate avg face length
+ polys = []
+ face_len = []
+ for element in elements:
+ face_len.append(element.face[1].length())
+ poly = element.get_poly()
+ polys.append(poly)
+ coll = PatchCollection(polys, edgecolors=ECOLOR, facecolors=FCOLOR)
+ axis.add_collection(coll)
+ face_len = sum(face_len)/len(face_len)
+
+ # store displacements we'll plot: list of [node, dict ux,uy,uz]
+ ulist = []
+ vals = []
+ for load in self.loads[self.time]:
+ if load.ltype in ['ux', 'uy', 'uz']:
+ alist = load.get_list()
+ for nodelist in alist:
+ node = nodelist[0]
+ if node in nodes:
+ ulist += [nodelist]
+ vals += list(nodelist[1].values())
+
+ # check min and max bounds
+ pmin = min(vals)
+ pmax = max(vals)
+
+ # make tick list for later plot, and color map
+ cbar_val = None
+ tick_list = []
+ cmap = None
+ if pmax != pmin:
+ # we have a range of values we're plotting
+ tick_list = linspace(pmin, pmax, 8)
+ cmap = plt.get_cmap(CMAP)
+ else:
+ cbar_val = pmin
+ pmax = pmin + 1.0
+ pmin = pmin - 1.0
+ tick_list = [pmin, cbar_val, pmax] # default 3 values to plot one solid color
+ cmap = colors.ListedColormap(['b', 'b']) # default to plot one val
+
+ # set color contours for arrows
+ cnorm = colors.Normalize(vmin=pmin, vmax=pmax)
+ scalarmap = cmx.ScalarMappable(norm=cnorm, cmap=cmap)
+ scalarmap.set_array([])
+
+ # make arrows for displacements
+ # extract nodes for elements if no nodes are selected
+ if len(nodes) == 0:
+ tmp = set()
+ for element in elements:
+ tmp.update(element.nodes)
+ nodes = list(tmp)
+ radials = [p.x for p in nodes]
+ axials = [p.y for p in nodes]
+ pvect = {'ux':geometry.Point(1, 0, 0),
+ 'uy':geometry.Point(0, 1, 0),
+ 'uz':geometry.Point(0, 0, 1)}
+ for [node, udict] in ulist:
+ for (key, val) in udict.items():
+ headw = face_len*0.4
+ headl = headw
+ point = geometry.Point(node.x, node.y, node.z)
+ unit = pvect[key]
+
+ # nonzero displacement, draw from node
+ tail = point
+ head = tail + unit*val
+ store_point = head
+ if val == 0:
+ # zero displ, draw to node
+ head = point
+ tail = head - unit*face_len
+ store_point = tail
+ delta = head - tail
+ radials.append(store_point.x)
+ axials.append(store_point.y)
+
+ colorVal = scalarmap.to_rgba(val)
+ plt.arrow(tail.y, tail.x, delta.y, delta.x,
+ color=colorVal, head_width=headw,
+ head_length=headl, length_includes_head=True)
+
+ # set the horizontal and vertical axes
+ base_classes.plot_set_bounds(plt, axials, radials)
+
+ # set units + titles
+ [d_unit, t_unit] = self.get_units('dist', 'time')
+ tstr = 'Constraints %s\nTime=%f%s' % (d_unit, self.time, t_unit)
+ plt.title(tstr)
+ plt.xlabel('axial, y'+d_unit)
+ plt.ylabel('radial, x'+d_unit)
+
+ # set the colorbar
+ cbar = plt.colorbar(scalarmap, orientation='vertical',
+ ticks=tick_list)
+ if cbar_val != None:
+ cbar.ax.set_yticklabels(['', str(cbar_val), ''])
+ base_classes.plot_finish(plt, fname, display)
+
+ else:
+ res = ''
+ if len(self.elements) == 0:
+ res = 'No elements exist! Try meshing your parts!'
+ else:
+ res = 'No elements are selected! Select some!'
+ print(res)
+
+[docs] def plot_multiple(self, fname='', display=True, title='',
+ styledict={'labels':['points', 'lines', 'areas', 'parts'],
+ 'items':['points, lines, areas']}):
+ """Plots items of type styledict['items'], labels styledict['labels']
+
+ Only the items currently selected will be plotted.
+ """
+ # http://stackoverflow.com/questions/470690/how-to-automatically-generate-n-distinct-colors/4382138#4382138
+ # http://stackoverflow.com/questions/2328339/how-to-generate-n-different-colors-for-any-natural-number-n
+ # start plotting
+ fig = plt.figure()
+ axis = fig.add_subplot(111)
+
+ # check if we're plotting parts or areas and pick a color map
+ colormaker, ids = None, None
+ color_plot = False
+ if 'parts' in styledict['items']:
+ color_plot = True
+ ids = [item.id for item in self.view.parts]
+ if 'areas' in styledict['items']:
+ styledict['items'].remove('areas')
+ if 'areas' in styledict['items']:
+ color_plot = True
+ ids = [item.id for item in self.view.areas]
+ if color_plot:
+ cmap = plt.get_cmap(GEOM_CMAP)
+ norm = colors.Normalize(vmin=min(ids), vmax=max(ids))
+ colormaker = cmx.ScalarMappable(norm=norm, cmap=cmap)
+
+ # plot the items
+ for item_type in ['points', 'lines', 'areas', 'parts']:
+ plot_on = item_type in styledict['items']
+ label_on = item_type in styledict['labels']
+ items = getattr(self.view, item_type)
+ if plot_on and label_on:
+ for item in items:
+ if color_plot and item_type in ['areas', 'parts']:
+ color = colormaker.to_rgba(item.id)
+ item.plot(axis, True, color)
+ else:
+ item.plot(axis, True)
+ elif plot_on:
+ for item in items:
+ if color_plot and item_type in ['areas', 'parts']:
+ color = colormaker.to_rgba(item.id)
+ item.plot(axis, False, color)
+ else:
+ item.plot(axis, False)
+ elif label_on:
+ for item in items:
+ item.label(axis)
+
+ # set the horizontal and vertical axes
+ points = self.view.points
+ radials = [point.x for point in points]
+ axials = [point.y for point in points]
+ base_classes.plot_set_bounds(plt, axials, radials)
+
+ # set units
+ [d_unit] = self.get_units('dist')
+
+ # show plot
+ plt.title(title)
+ plt.xlabel('axial, y'+d_unit)
+ plt.ylabel('radial, x'+d_unit)
+ axis.set_aspect('equal')
+ base_classes.plot_finish(plt, fname, display)
+
+[docs] def plot_parts(self, fname='', display=True, title='Parts', label=True):
+ """Plots selected parts
+
+ Defaults to displaying labels on: parts and filling all part areas.
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown. Default = True
+ title (str): the plot's title
+ label (bool): if True all Part objects will be labeled
+ """
+ item = 'parts'
+ styledict = {'items':[], 'labels':[]}
+ if label:
+ styledict['items'].append(item)
+ styledict['labels'].append(item)
+ else:
+ styledict['items'].append(item)
+
+ self.plot_multiple(fname, display, title, styledict)
+
+[docs] def plot_areas(self, fname='', display=True, title='Areas', label=True):
+ """Plots selected areas
+
+ Defaults to displaying labels on: areas and filling all part areas.
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown. Default = True
+ title (str): the plot's title
+ label (bool): if True all areas will be labeled
+ """
+ item = 'areas'
+ styledict = {'items':[], 'labels':[]}
+ if label:
+ styledict['items'].append(item)
+ styledict['labels'].append(item)
+ else:
+ styledict['items'].append(item)
+
+ self.plot_multiple(fname, display, title, styledict)
+
+[docs] def plot_lines(self, fname='', display=True, title='Lines', label=True):
+ """Plots selected lines and arcs
+
+ Defaults to displaying labels on: lines and arcs
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown. Default = True
+ title (str): the plot's title
+ label (bool): if True all lines and arcs will be labeled
+ """
+ item = 'lines'
+ styledict = {'items':[], 'labels':[]}
+ if label:
+ styledict['items'].append(item)
+ styledict['labels'].append(item)
+ else:
+ styledict['items'].append(item)
+
+ self.plot_multiple(fname, display, title, styledict)
+
+[docs] def plot_points(self, fname='', display=True, title='Points', label=True):
+ """Plots selected points
+
+ Defaults to displaying labels on: points
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown. Default = True
+ title (str): the plot's title
+ label (bool): if True all points will be labeled
+ """
+ item = 'points'
+ styledict = {'items':[], 'labels':[]}
+ if label:
+ styledict['items'].append(item)
+ styledict['labels'].append(item)
+ else:
+ styledict['items'].append(item)
+
+ self.plot_multiple(fname, display, title, styledict)
+
+[docs] def plot_geometry(self, fname='', display=True, title='Geometry', pnum=True,
+ lnum=True, anum=True, afill=True):
+ """Plots selected geometry items: areas, lines, points.
+
+ Defaults to displaying labels on: areas, lines, and points, and filling
+ in the areas.
+
+ Args:
+ fname (str): png image file prefix, if given, image will be saved
+ display (bool): if True, interactive plot will be shown, if False
+ plot will not be shown. Default = True
+ title (str): the plot's title
+ pnum (bool): if True all Point objects will be labeled
+ lnum (bool): if True all Line and Arc objects will be labeled
+ anum (bool): if True all Area objects will be labeled
+ afill (bool): if True all Area objects will be filled
+ """
+ styledict = {'items':[], 'labels':[]}
+ if pnum:
+ styledict['items'].append('points')
+ styledict['labels'].append('points')
+ if lnum:
+ styledict['items'].append('lines')
+ styledict['labels'].append('lines')
+ if anum:
+ styledict['labels'].append('areas')
+ if afill:
+ styledict['items'].append('areas')
+ elif anum == False and afill:
+ styledict['items'].append('areas')
+
+ self.plot_multiple(fname, display, title, styledict)
+
+ @staticmethod
+ def __get_cname(items):
+ """Returns a component name prefix, for labeling lists of items.
+
+ Args:
+ items (list): a list of the items we'll be making a component of
+ """
+ cname = ''
+ if len(items) == 1:
+ if items[0] == 'all':
+ cname = 'all'
+ else:
+ cname = items[0].get_name()
+ else:
+ cname = items[0].get_name()+'-'+items[-1].get_name()
+ return cname
+
+ def __get_make_comp(self, comp):
+ """Stores component if it doesn't exist, returns it if it does.
+
+ Args:
+ comp (Component): component we want to store/get
+ """
+ if comp not in self.components:
+ comp = self.components.append(comp)
+ else:
+ ind = self.components.index(comp)
+ comp = self.components[ind]
+ return comp
+
+ def __get_make_surfint(self, surfint):
+ """Stores surfac interaction if it doesn't exist, returns it if it does.
+
+ Args:
+ surfint (connectors.SurfaceInteraction): item to get or make
+ """
+ items = [item for item in self.surfints if item.int_type == surfint.int_type]
+ if surfint.int_type == 'LINEAR':
+ for item in items:
+ if item.k == surfint.k:
+ return item
+ elif surfint.type == 'EXPONENTIAL':
+ for item in items:
+ if item.c0 == surfint.c0 and item.p0 == surfint.p0:
+ return item
+ surfint = self.surfints.append(surfint)
+ return surfint
+
+[docs] def register(self, item):
+ """Adds an item to the feamodel.
+
+ Item is added to the correct list and its id is updated.
+ Allowed Items:
+
+ * Point
+ * Line
+ * Arc
+ * SignLine
+ * SignArc
+ * LineLoop
+ * Area
+ * Part
+ """
+ class_name = item.__class__.__name__
+ class_to_list = {'Point': 'points',
+ 'Line': 'lines',
+ 'Arc': 'lines',
+ 'SignLine': 'signlines',
+ 'SignArc': 'signlines',
+ 'LineLoop': 'lineloops',
+ 'Area': 'areas',
+ 'Part': 'parts'}
+ if class_name in class_to_list:
+ list_name = class_to_list[class_name]
+ getattr(self, list_name).append(item)
+ else:
+ print('ERROR: the item you passed must be a geometry item!')
+ print('You passed a %s' % class_name)
+ message = ("Allowed items: Point, Line, Arc, SignLine, SignArc, "
+ "LineLoop, Area, Part")
+ print(message)
+
+ def __add_load(self, load, time):
+ """Adds a load to the FeaModel.
+
+ Args:
+ load (Load or Load_Linear): th eload to add
+ time (float): the load's time point. The first time point is 1.0
+ """
+ if time in self.loads:
+ self.loads[time].append(load)
+ else:
+ self.loads[time] = [load]
+ return load
+
+[docs] def scale(self, unitstr, point_first=None, point_last=None):
+ """Scales the points from [point_first, ..., point_last] using unitstr
+
+ If point_first and point_last are both None, all points are scaled.
+
+ Args:
+ unitstr (str): string scalar 'fromunit-tounit' using the below units:
+
+ * mm, m, in, ft
+ * Examples 'mm-m' 'm-in' 'ft-mm' 'in-ft'
+ * Default value is '' and does not apply a scale factor
+
+ point_first (Point or None or str): the first point to scale,
+ string point names may be passed
+ point_last (Point or None or str): the last point to scale,
+ string point names may be passed
+ """
+ # convert to points if passing in string point names 'P0'
+ if isinstance(point_first, str):
+ point_first = self.get_item(point_first)
+ if isinstance(point_last, str):
+ point_last = self.get_item(point_last)
+
+ ind_first = 0
+ if point_first != None:
+ ind_first = self.points.index(point_first)
+ ind_last = -1
+ if point_last != None:
+ ind_last = self.points.index(point_last)
+
+ unit_size = {'m':1.0, 'mm':.001, 'ft':.3048, 'in':0.0254}
+ units = unitstr.split('-')
+ from_unit, to_unit = units[0], units[1]
+ scalar = unit_size[from_unit]/unit_size[to_unit]
+ # scale all points in set
+ # update all dependent arcs
+ # update all dependent line loops
+ # update all dependent areas
+ # update all dependent parts
+
+[docs] def set_gravity(self, grav, items):
+ """Sets gravity on the elements in items.
+
+ Assumes gravity acts in the -x direction with magnitude grav.
+
+ Args:
+ grav (float): gravity acceleration, MUST BE POSTIVE
+ items (Area or Part or list): items gravity acts on
+
+ - Area: gravity acts on elements in this area
+ - Part: gravity acts on elements in this part
+ - list of Part or Area: gravity acts on their elements
+ """
+ items = self.get_items(items)
+ ctype = 'elements'
+ cname = self.__get_cname(items)
+
+ # make compoenet
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # make load
+ ltype = 'gravity'
+ time = self.time
+ load = loads.ConstLoad(ltype, comp, grav)
+ self.__add_load(load, time)
+ return load
+
+[docs] def set_rpm(self, rpm, items):
+ """Sets rpm rotation load on the items.
+
+ Args:
+ rpm (float): rotation per minute
+ items (Area or Part or list): items rotation acts on
+
+ - Area: rotation acts on elements in this area
+ - Part: rotation acts on elements in this part
+ - list of Part or Area: rotation acts on their elements
+ """
+ # applies rpm to the items
+ items = self.get_items(items)
+ ctype = 'elements'
+ cname = self.__get_cname(items)
+
+ # make compoenet
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # make load
+ ltype = 'rpm'
+ time = self.time
+ load = loads.ConstLoad(ltype, comp, rpm)
+ self.__add_load(load, time)
+ return load
+
+[docs] def set_radps(self, radps, items):
+ """Sets radians per second rotation load on the items.
+
+ Args:
+ radps (float): radians per second
+ items (Area or Part or list): items rotation acts on
+
+ - Area: rotation acts on elements in this area
+ - Part: rotation acts on elements in this part
+ - list of Part or Area: rotation acts on their elements
+ """
+ items = self.get_items(items)
+ ctype = 'elements'
+ cname = self.__get_cname(items)
+
+ # make compoenet
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # make load
+ ltype = 'radps'
+ time = self.time
+ load = loads.ConstLoad(ltype, comp, radps)
+ self.__add_load(load, time)
+ return load
+
+[docs] def set_fluid_press(self, items, rho, g, xo, po):
+ """This sets a fluid presure load on items.
+
+ This fluid pressure is dependednt on the x axis.
+ g must be positive.
+
+ - P = f(x)
+ - P(xo) = po
+ - P(x) = po + rho*g*(xo - x)
+
+ Args:
+ items (str or Line or Arc or list): items to set pressure on
+
+ - str: string name of Line or Arc item, for example 'L0'
+ - Line or Arc: set presssure on this
+ - list or Line or Arc: set pressure on these
+ rho (float): fluid density in mass/volume
+ g (+float): gravity in dist/(t^2), MUST BE POSITIVE
+ xo (float): sea level height, MUST BE POSITIVE
+ po (float): sea level pressure, MUST BE POSITIVE
+ """
+ items = self.get_items(items)
+ ctype = 'faces'
+ cname = self.__get_cname(items)
+
+ # make compoenet
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # make load
+ ltype = 'press_fluid'
+ mult = rho*g
+ load = loads.LinearLoad(ltype, comp, po, mult, xo)
+ self.__add_load(load, self.time)
+ return load
+
+[docs] def set_load(self, ltype, items, lval, ldir=None):
+ """Sets a pressure or force load on line(s).
+
+ The force load is divided by the number of nodes and applied on
+ each node.
+
+ The pressure load is applied to the child faces under the line.
+ Positive is compression, negative is tension.
+
+ Args:
+ ltype (str): 'press' or 'force'
+ items (str or Line or Arc or list): items to set load on
+
+ - str: string name of Line or Arc item, for example 'L0'
+ - Line or Arc: set load on this
+ - list or Line or Arc: set load on these
+ lval (float): load value.
+
+ - For ltype = 'press' this is in force/area units
+ - For ltype = 'force' this is in force units
+ ldir (None or str): load direction. Defaults to None
+
+ - str: when ltype='load', we need to set ldir to 'x' or 'y'
+ - None: when ltype='press'
+ """
+ items = self.get_items(items)
+ ctype = 'nodes'
+ if ltype == 'press':
+ ctype = 'faces'
+ cname = self.__get_cname(items)
+
+ # make compoenet
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # make load
+ if ltype == 'force':
+ ltype = 'f'+ldir # for example fx
+ load = loads.ConstLoad(ltype, comp, lval)
+ self.__add_load(load, self.time)
+ return load
+
+[docs] def set_constr(self, ltype, items, axis, val=0.0):
+ """Sets a displacement constraint on the passed item(s).
+
+ Args:
+ ltype (str): 'fix' or 'displ'
+
+ - 'fix': val arg should not be passed
+ - 'displ': val arg must be passed
+ items (str, SignLine, SignArc, Point or list): item(s) to apply the
+ constraint on
+
+ - str: this string must be a line or point name
+ - SignLine or SignArc: the constraint will be applied
+ to this item
+ - Point: the constraint will be applied on this item
+ - list of SignLine or SignArc: the constraint will be appplied
+ on these
+ - list of Point: the constraint will be appplied on these
+ - list of str names: pass line or point names, constr
+ applied on them
+
+ axis (str): load axis, 'x' or 'y'
+ val (float): displacement value, defaults to 0.0, only needs to be
+ set when ltype = 'displ'
+ """
+ items = self.get_items(items)
+ cname = self.__get_cname(items)
+ ctype = 'nodes'
+
+ # make compoenet
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # make load
+ ltype = 'u'+axis # for example ux
+ load = loads.ConstLoad(ltype, comp, val)
+ self.__add_load(load, self.time)
+ return load
+
+[docs] def set_contact_linear(self, master_lines, slave_lines, kval, many_si=False):
+ """Sets contact between master and slave lines.
+
+ Slave lines are on the more compliant or more curved object.
+
+ Args:
+ master_lines (list): list of SignLine or SignArc
+ slave_lines (list): list of SignLine or SignArc
+ kval (float): stiffness, 5 to 50 times the youngs modulus
+ of the touching matl
+ many_si (bool): True, make unique surface interaction for every contact
+ False, use existing surface interaction if we can
+ """
+
+ master_items = self.get_items(master_lines)
+ master_cname = self.__get_cname(master_items)
+ ctype = 'faces'
+ master_comp = components.Component(master_items, ctype, master_cname)
+ master_comp = self.__get_make_comp(master_comp)
+
+ slave_items = self.get_items(slave_lines)
+ slave_cname = self.__get_cname(slave_items)
+ ctype = 'faces'
+ slave_comp = components.Component(slave_items, ctype, slave_cname)
+ slave_comp = self.__get_make_comp(slave_comp)
+
+ surf_int = connectors.SurfaceInteraction('LINEAR', kval)
+ if many_si:
+ surf_int = self.surfints.append(surf_int)
+ else:
+ surf_int = self.__get_make_surfint(surf_int)
+
+ cont = connectors.Contact(master_comp, slave_comp, surf_int, True)
+ self.contacts.append(cont)
+
+
+[docs] def set_eshape(self, eshape='quad', eorder=2):
+ """Sets the element shape and order to use when meshing the model.
+
+ Args:
+ eshape (str): element shape
+
+ - 'quad': quadrilaterials
+ - 'tri': triangles
+ eorder (int): element order, default=2
+
+ - 1: corner nodes only (3 and 4 noded elements)
+ - 2: corder nodes and midside nodes (6 and 8 noded elements)
+ """
+ self.eshape = eshape # quad or tri
+ self.eorder = eorder # 1 or 2
+
+[docs] def set_etype(self, etype, items, thick=None):
+ """Sets the element type, and thickness on areas or parts.
+
+ Args:
+ etype (str): element type
+
+ - 'plstress': plane stress
+ - 'plstrain': plane strain
+ - 'axisym': axisymmetric
+ items (str or Area or Part or list): set element type on these
+
+ - str: string name of Area or Part item. Example: 'A0', 'PART0'
+ - Area: set element type on the elements in this area
+ - Part: set element type on the elements in this part
+ - list of Area or Part: set element type on their elements
+ thick (float or None): element thickness
+
+ - None: default, used for axisymmetric element type
+ - float: thickness value to use for plane stress or plane
+ strain elements
+ """
+ items = self.get_items(items)
+
+ # set the element types on the areas, this is used to fix
+ # elements when importing them from the inp file
+ for item in items:
+ if isinstance(item, geometry.Area):
+ item.set_etype(etype)
+ if isinstance(item, partmodule.Part):
+ for area in item.areas:
+ area.set_etype(etype)
+
+ # set a thickness component if needed
+ if etype != 'axisym' and thick != None:
+
+ # make component for element thickness
+ cname = self.__get_cname(items)
+ ctype = 'nodes'
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # add load
+ ltype = 'nodal_thickness'
+ time = 0.0
+ load = loads.ConstLoad(ltype, comp, thick)
+ self.__add_load(load, time)
+ return load
+ return None
+
+[docs] def set_matl(self, matl, items):
+ """Sets the material on Part or Area items.
+
+ Args:
+ matl (Material): material to assign to items
+ items (Part or Area or list): items that will have matl assigned
+
+ - Part: assign matl to this
+ - Area: assign matl to this
+ - list of Area or Part: assign matl to these
+ """
+ items = self.get_items(items)
+ cname = self.__get_cname(items)
+
+ # store the material if it's not already
+ if matl not in self.matls:
+ matl = self.matls.append(matl)
+
+ # make component
+ ctype = 'elements'
+ comp = components.Component(items, ctype, cname)
+ comp = self.__get_make_comp(comp)
+
+ # make load
+ ltype = 'matl'
+ time = 0.0
+ load = loads.ConstLoad(ltype, comp, matl)
+ self.__add_load(load, time)
+ return load
+
+ def __read_inp(self, fname):
+ """Reads in the mesh from a Calculix inp file.
+
+ All nodes, elements, and faces are read. Child nodes, elements and
+ faces are assigned to their geometry parents (points, lines, arcs,
+ areas, parts.
+
+ Args:
+ fname (str): file name to read, must include '.inp' extension
+ """
+ f = open(fname, 'r')
+ mode = None
+ set_name = None
+ set_type = None
+
+ items = [] # holder for nodes or elements in nsets or esets
+ N = base_classes.Meshlist() # store nodes
+ E = base_classes.Meshlist() # store elements, allows renumbering before putting int model
+ F = [] # store faces
+ sets = {'E':{}, 'N':{}} # store sets
+ etype = ''
+ Dict_NodeIDs={}
+ Dict_ElemIDs={}
+
+ # read in input file
+ for line in f:
+ if line[0] != '*':
+ if mode == 'nmake':
+ L = line.split(',')
+ L = [a.strip() for a in L]
+ (nnum, x, y, z) = (int(L[0]), float(L[1]), float(L[2]), float(L[3]))
+ node = mesh.Node(nnum, x, y, z)
+ N.append(node)
+ Dict_NodeIDs[nnum]=node
+ elif mode == 'emake':
+ L = line.split(',')
+ L = [int(a.strip()) for a in L]
+ enum = L[0]
+ nlist = [Dict_NodeIDs[a] for a in L[1:]]
+ e = mesh.Element(enum, etype, nlist)
+ faces = e.faces
+ E.append(e)
+ Dict_ElemIDs[enum]=e
+ F += faces
+ sets[set_type][set_name].append(e)
+ elif mode == 'set':
+ L = line.split(',')
+ L = [a.strip() for a in L]
+ L = [int(a) for a in L if a != '']
+ items = []
+ if set_type == 'E':
+ items = [Dict_ElemIDs[a] for a in L if a in Dict_ElemIDs.keys()]
+ elif set_type == 'N':
+ items = [Dict_NodeIDs[a] for a in L]
+ if items == [None]*len(items):
+ pass # the elements were not found
+ else:
+ sets[set_type][set_name] += items
+
+ # mode setting
+ if '*Node' in line or '*NODE' in line:
+ mode = 'nmake'
+ elif '*Element' in line or '*ELEMENT' in line:
+ L = line.split(',') # split it based on commas
+ e = L[1].split('=')
+ etype = e[1]
+
+ # exclude T elements made in gmsh
+ if etype[0] != 'T':
+ e = L[2].split('=')
+ set_name = e[1].strip()
+ set_type = 'E'
+ sets[set_type][set_name] = []
+ mode = 'emake'
+ else:
+ mode = None
+ elif '*ELSET' in line:
+ L = line.split(',')
+ e = L[1].split('=')
+ set_name = e[1].strip()
+ set_type = 'E'
+ sets[set_type][set_name] = []
+ mode = 'set'
+ elif '*NSET' in line:
+ L = line.split(',')
+ e = L[1].split('=')
+ set_name = e[1].strip()
+ set_type = 'N'
+ sets[set_type][set_name] = []
+ mode = 'set'
+ f.close()
+
+ # loop through sets and remove empty sets
+ # store sets to delete
+ todel = []
+ for (set_type, set_dict) in sets.items():
+ for (set_name, item_list) in set_dict.items():
+ if item_list == []:
+ todel.append({'set_type':set_type, 'set_name':set_name})
+ # delete the empty sets
+ for adict in todel:
+ (set_type, set_name) = (adict['set_type'], adict['set_name'])
+ del sets[set_type][set_name]
+ #print('Empty set type:%s name:%s deleted' % (set_type, set_name))
+
+ # this resets the min element to number 1
+ if E.get_minid() > 1:
+ E.set_minid(1)
+
+ #-----------------------------------
+ # Node and element assignment back onto parts, areas, lines, points
+ #-----------------------------------
+ # assign elements + nodes to feamodel
+ self.elements = E
+ self.faces = F
+
+ # remove arc center nodes from imported node set
+ # those nodes have no elements under them
+ torem = []
+ for node in N:
+ if len(node.elements) == 0:
+ torem.append(node)
+ for node in torem:
+ N.remove(node)
+ # assign nodes to feamodel
+ self.nodes = N
+
+ for part in self.parts:
+ # assign part element and node sets
+ pname = part.get_name()
+ part.elements = sets['E'][pname]
+ part.nodes = sets['N'][pname]
+
+ # assign all nodes and elements to areas, fix the element types
+ for area in part.areas:
+ aname = area.get_name()
+ area.elements = sets['E'][aname]
+ area.nodes = sets['N'][aname]
+ area.set_child_ccxtypes() #paint element types on elements
+
+ # assign the child nodes to points
+ pts = part.points
+ for point in pts:
+ ndist = []
+ for node in part.nodes:
+ p_tmp = geometry.Point(node.x, node.y)
+ dist = point - p_tmp
+ dist = dist.length()
+ ndist.append({'dist':dist, 'node':node})
+ # sort the list by dist, sorts low to high
+ ndist = sorted(ndist, key=lambda k: k['dist'])
+ point.nodes = [ndist[0]['node']]
+ #print('Point %s = node %s' % (pt, pt.nodes))
+
+ # assign the nodes and n1 and faces to lines
+ slines = part.signlines
+ for sline in slines:
+ lname = sline.get_name()
+ area = sline.lineloop.parent
+ nodes = sets['N'][sline.line.get_name()]
+ n1 = [n for n in nodes if n.order == 1]
+ sline.nodes = nodes
+ sline.n1 = n1
+ # make a set of all faces that contain the line's node
+ allfaces = set()
+ for node in nodes:
+ allfaces.update(node.faces)
+ # make a set of all faces on the line only
+ linefaces = set()
+ for face in allfaces:
+ if set(face.nodes).issubset(set(nodes)):
+ linefaces.add(face)
+ # and only the ones whose element centers are in the area
+ faces = set()
+ for face in linefaces:
+ is_lineface = area.contains_point(face.element.center)
+ if is_lineface:
+ faces.add(face)
+ sline.faces = list(faces)
+ print('Elements: %i' % len(E))
+ print('Nodes: %i' % len(N))
+ print('Done reading Calculix/Abaqus .inp file')
+
+[docs] def mesh(self, size=1.0, meshmode='fineness', mesher='gmsh'):
+ """Meshes all parts.
+
+ Args:
+ size (float):
+
+ - if meshmode == 'fineness' (default):
+ - mesh size is adapted to geometry size
+ - set size = 0.0001 - 1.0, to define how fine the mesh is.
+ - Low numbers are very fine, higher numbers are coarser.
+
+ - if meshmode == 'esize':
+ - element size is kept constant
+ - choose it depending on geometry size
+ - it should be reduced e.g. at arcs with small radius, by calling line.esize function
+
+ meshmode (str):
+
+ - 'fineness': adapt mesh size to geometry
+ - 'esize': keep explicitly defined element size
+
+ meshmode is changed to 'esize' is used if esize property is set to points or lines
+
+ mesher (str): the mesher to use
+
+ - 'gmsh': mesh with Gmsh, this is reccomended, it allows holes
+ - 'cgx': mesh with Calculix cgx, it doesn't allow holes
+ """
+
+
+ #check if element size is set to points and change meshmode if necessary
+ for pt in self.points:
+ if pt.esize != None:
+ if meshmode=='fineness': print('meshmode is changed to esize, because elementsize was defined on points!')
+ meshmode = 'esize'
+
+
+ #if meshmode esize is chosen: ediv's on lines and arcs are transformed to element sizes on start and end point
+ if meshmode == 'esize':
+ for line in self.lines:
+ if line.ediv != None:
+ line.pt(0).set_esize(line.length()/line.ediv)
+ line.pt(1).set_esize(line.length()/line.ediv)
+
+
+ if mesher == 'gmsh':
+ self.__mesh_gmsh(size, meshmode)
+ elif mesher == 'cgx':
+ self.__mesh_cgx(size)
+
+ def __mesh_gmsh(self, size, meshmode, timeout=20):
+ """Meshes all parts using the Gmsh mesher.
+
+ Args:
+
+ size (float):
+
+ - if meshmode == 'fineness' (default):
+ - mesh size is adapted to geometry size
+ - set size = 0.0001 - 1.0, to define how fine the mesh is.
+ - Low numbers are very fine, higher numbers are coarser.
+
+ - if meshmode == 'esize':
+ - element size is kept constant
+ - choose it depending on geometry size
+
+ meshmode (str):
+
+ - 'fineness': adapt mesh size to geometry
+ - 'esize': keep explicitly defined element size
+
+ timeout (int): time in seconds before the process throws a
+ subprocess.TimeoutExpired
+
+ """
+ geo = []
+ ids = {}
+ ids['line'] = {}
+ ids['plane_surface'] = {}
+
+ # write all points
+ for pt in self.points:
+ txtline = 'Point(%i) = {%f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0)
+
+ if meshmode == 'esize':
+ #add element size to points
+ if pt.esize == None:
+ txtline = txtline.replace('}', ', %f}' % (size if self.eshape=='tri' else size*2.))
+ #txtline = txtline.replace('}', ', %f}' % (size))
+ else:
+ txtline = txtline.replace('}', ', %f}' % (pt.esize if self.eshape=='tri' else pt.esize*2.))
+ #txtline = txtline.replace('}', ', %f}' % (pt.esize))
+ geo.append(txtline)
+
+ # start storing an index number
+ ind = self.points[-1].id + 1
+
+ # write all lines
+ for line in self.lines:
+ lnum = line.id
+ ids['line'][lnum] = ind
+ pt1 = line.pt(0).id
+ pt2 = line.pt(1).id
+ txtline = ''
+ if isinstance(line, geometry.Arc):
+ ctr = line.actr.id
+ txtline = 'Circle(%i) = {%i, %i, %i};' % (ind, pt1, ctr, pt2)
+ else:
+ txtline = 'Line(%i) = {%i,%i};' % (ind, pt1, pt2)
+ geo.append(txtline)
+
+ # set division if we have it
+ if line.ediv != None and meshmode=='fineness':
+ ndiv = line.ediv+1
+ esize = line.length()/line.ediv
+ if self.eshape == 'quad':
+ ndiv = line.ediv/2+1
+ esize = esize*2
+ # this is needed because quad recombine
+ # splits 1 element into 2
+ txtline = 'Transfinite Line{%i} = %i;' % (ind, ndiv)
+ print('LINE ELEMENT SIZE: %f, MAKES %i ELEMENTS'
+ % (line.length()/line.ediv, line.ediv))
+ geo.append(txtline)
+ geo.append('Characteristic Length {%i,%i} = %f;'
+ % (pt1, pt2, esize))
+ ind += 1
+
+ # write all areas
+ for area in self.areas:
+ if area.closed:
+ aid = area.id
+ aname = area.get_name()
+ loop_ids = []
+ loops = [area.exlines] + area.holes
+ for loop in loops:
+ txtline = 'Line Loop(%i) = ' % (ind)
+ loop_ids.append(str(ind))
+ line_ids = []
+ for sline in loop:
+ lid = ids['line'][sline.line.id]
+ prefix = ''
+ if sline.sign == -1:
+ prefix = '-'
+ line_ids.append('%s%i' % (prefix, lid))
+ txtline = txtline + '{'+ ','.join(line_ids)+'};'
+ geo.append(txtline)
+ ind += 1
+ loop_ids = ','.join(loop_ids)
+ geo.append('Plane Surface(%i) = {%s};' % (ind, loop_ids))
+ ids['plane_surface'][aid] = ind
+ geo.append("Physical Surface('%s') = {%i};" % (aname, ind))
+ ind += 1
+
+ # write part area components
+ for part in self.parts:
+ # make components for each part
+ txtline = "Physical Surface('%s') = " % (part.get_name())
+ area_ids = []
+ for area in part.areas:
+ if area.closed:
+ aid = ids['plane_surface'][area.id]
+ area_ids.append(str(aid))
+ txtline = txtline + '{' + ','.join(area_ids) + '};'
+ geo.append(txtline)
+
+ # write all line componenets so we can get nodes out
+ for line in self.lines:
+ lid = ids['line'][line.id]
+ txtline = "Physical Line('%s') = {%i};" % (line.get_name(), lid)
+ geo.append(txtline)
+
+ # write node componenets
+ # ERROR: node list is not produced by gmsh
+ for point in self.points:
+ pname = point.get_name()
+ txtline = "Physical Point('%s') = {%i};" % (pname, point.id)
+ geo.append(txtline)
+
+ # set the meshing options
+ if meshmode == 'fineness':
+ geo.append('Mesh.CharacteristicLengthFactor = '+str(size)+'; //mesh fineness')
+ geo.append('Mesh.RecombinationAlgorithm = 1; //blossom')
+
+ if self.eshape == 'quad':
+ geo.append('Mesh.RecombineAll = 1; //turns on quads')
+ geo.append('Mesh.SubdivisionAlgorithm = 1; // quadrangles only')
+ #geo.append('Mesh.RecombinationAlgorithm = 1; //turns on blossom needed for quad')
+
+ order = self.eorder
+ geo.append('Mesh.CharacteristicLengthExtendFromBoundary = 1;')
+ geo.append('Mesh.CharacteristicLengthMin = 0;')
+ geo.append('Mesh.CharacteristicLengthMax = 1e+022;')
+ # use this so small circles are meshed finely
+ # this is broken in Gmsh 4.3.0
+ # geo.append('Mesh.CharacteristicLengthFromCurvature = 1;')
+ geo.append('Mesh.CharacteristicLengthFromPoints = 1;')
+ # geo.append('Mesh.Algorithm = 2; //delauny') #okay for quads
+ geo.append('Mesh.Algorithm = 8; //delquad = delauny for quads')
+ geo.append('Mesh.ElementOrder = '
+ +str(order)
+ +'; //linear or second set here')
+ if order == 2:
+ geo.append('Mesh.SecondOrderIncomplete=1; //no face node w/ 2nd order')
+ geo.append('Mesh.SaveGroupsOfNodes = 1; // save node groups')
+
+ # write geo file to the local directory
+ fname = self.fname+'.geo'
+ fout = self.fname+'.inp'
+ outfile = open(fname, 'w')
+ for line in geo:
+ #print (line)
+ outfile.write(line+'\n')
+ outfile.close()
+ print('File: %s was written' % fname)
+
+ # run file in bg mode, -2 is 2d mesh, makes required inp file
+ runstr = "%s %s -2 -o %s" % (environment.GMSH, fname, fout)
+ print(runstr)
+ subprocess.check_call(runstr, timeout=timeout, shell=True)
+ print('File: %s was written' % fout)
+ print('Meshing done!')
+
+ # write gmsh msh file, for manual checking only
+ # not required by pycalculix
+ runstr = "%s %s -2 -o %s" % (environment.GMSH, fname,
+ self.fname+'.msh')
+ subprocess.check_call(runstr, timeout=timeout, shell=True)
+ print('File: %s.msh was written' % self.fname)
+
+ # read in the calculix mesh
+ self.__read_inp(self.fname+'.inp')
+
+ def __mesh_cgx(self, size, meshmode, timeout=20):
+ """Meshes all parts using the Calculix cgx mesher.
+
+ Args:
+
+ size (float):
+
+ - if meshmode == 'fineness' (default):
+ - mesh size is adapted to geometry size
+ - set size = 0.0001 - 1.0, to define how fine the mesh is.
+ - Low numbers are very fine, higher numbers are coarser.
+
+ - if meshmode == 'esize': NOT TESTED WITH CGX
+ - element size is kept constant
+ - choose it depending on geometry size
+
+ meshmode (str):
+
+ - 'fineness': adapt mesh size to geometry
+ - 'esize': keep explicitly defined element size NOT TESTED WITH CGX
+
+ timeout (int): time in seconds before the process throws a
+ subprocess.TimeoutExpired
+ """
+ fbd = []
+ comps = []
+ cfiles = []
+
+ # Calculix CGX elements
+ # axisymmetric
+ cgx_elements = {}
+ cgx_elements['tri2axisym'] = 'tr6c'
+ cgx_elements['tri1axisym'] = 'tr3c'
+ cgx_elements['quad2axisym'] = 'qu8c'
+ cgx_elements['quad1axisym'] = 'qu4c'
+ # plane stress
+ cgx_elements['tri2plstress'] = 'tr6s'
+ cgx_elements['tri1plstress'] = 'tr3s'
+ cgx_elements['quad2plstress'] = 'qu8s'
+ cgx_elements['quad1plstress'] = 'qu4s'
+ # plane strain
+ cgx_elements['tri2plstrain'] = 'tr6e'
+ cgx_elements['tri1plstrain'] = 'tr3e'
+ cgx_elements['quad2plstrain'] = 'qu8e'
+ cgx_elements['quad1plstrain'] = 'qu4e'
+
+ num = 1.0/size
+ emult = int(round(num)) # this converts size to mesh multiplier
+
+ # write all points
+ for point in self.points:
+ linestr = 'pnt %s %f %f %f' % (point.get_name(), point.x, point.y, 0.0)
+ fbd.append(linestr)
+ # gmsh can't make node componenets so don't do it in cgx
+ #L = 'seta %s p %s' % (pt.get_name(), pt.get_name())
+ #comps.append(L)
+
+ # write all lines
+ for line in self.lines:
+ lname = line.get_name()
+ pt1 = line.pt(0).get_name()
+ pt2 = line.pt(1).get_name()
+ linestr = ''
+ if isinstance(line, geometry.Arc):
+ # line is arc
+ pctr = line.actr.get_name()
+ linestr = 'line %s %s %s %s' % (lname, pt1, pt2, pctr)
+ else:
+ # straight line
+ linestr = 'line %s %s %s' % (lname, pt1, pt2)
+ # set division if we have it
+ if line.ediv != None:
+ ndiv = self.eorder*line.ediv
+ linestr += ' '+str(int(ndiv))
+ fbd.append(linestr)
+ linestr = 'seta %s l %s' % (lname, lname)
+ comps.append(linestr)
+ cfiles.append(lname)
+
+ # write all areas
+ for area in self.areas:
+ if area.closed:
+ linestr = 'gsur '+area.get_name()+' + BLEND '
+ line_ids = []
+ for line in area.signlines:
+ lname = '+ '+line.get_name()
+ if line.sign == -1:
+ lname = '- '+line.get_name()[1:]
+ line_ids.append(lname)
+ linestr = linestr + ' '.join(line_ids)
+ fbd.append(linestr)
+ # add area component, nodes + elements
+ astr = 'seta %s s %s' % (area.get_name(), area.get_name())
+ fbd.append(astr)
+ cfiles.append(area.get_name())
+
+ # write part area components
+ for apart in self.parts:
+ # make components for each part
+ # seta P0 s s0 s1
+ line = 'seta %s s ' % apart.get_name()
+ cfiles.append(apart.get_name())
+ area_ids = []
+ for area in apart.areas:
+ if area.closed:
+ area_ids.append(area.get_name())
+ line = line + ' '.join(area_ids)
+ fbd.append(line)
+
+ # mesh all areas
+ for area in apart.areas:
+ aname = area.get_name()
+ estr = self.eshape+str(self.eorder)+area.etype
+ cgx_etype = cgx_elements[estr]
+ fbd.append('elty %s %s' % (aname, cgx_etype))
+ fbd.append('div all mult %i' % emult)
+ fbd.append('mesh all')
+
+ # save mesh file
+ fbd.append('send all abq')
+
+ # add line and area components
+ fbd += comps
+
+ # save component node and element sets
+ for comp in cfiles:
+ # select nodes under
+ fbd.append('comp %s do' % (comp,))
+ fbd.append('send %s abq names' % (comp,))
+
+ # this orients the view correctly
+ # this is the same as switching to the z+ orientation
+ # y us axial, x is radial
+ fbd.append('rot z')
+ fbd.append('rot c -90')
+ fbd.append('plot e all')
+
+ # write fbd file to the local directory
+ fname = self.fname+'.fbd'
+ f = open(fname, 'w')
+ for line in fbd:
+ #print (line)
+ f.write(line+'\n')
+ f.close()
+ print('File: %s was written' % fname)
+
+ # run file in bg mode
+ runstr = "%s -bg %s" % (environment.CGX, fname)
+ p = subprocess.check_call(runstr, timeout=timeout, shell=True)
+ print('Meshing done!')
+
+ # assemble the output files into a ccx input file
+ inp = []
+ files = ['all.msh']
+ files += [f+'.nam' for f in cfiles]
+ for fname in files:
+ infile = open(fname, 'r')
+ for line in infile:
+ # cgx adds E and N prfixes on sets after =, get rid of these
+ if '=' in line and fname != 'all.msh':
+ L = line.split('=')
+ line = L[0] + '=' + L[1][1:]
+ inp.append(line.strip())
+ else:
+ inp.append(line.strip())
+ infile.close()
+
+ # delete file
+ os.remove(fname)
+
+ # write out inp file
+ fname = self.fname+'.inp'
+ outfile = open(fname, 'w')
+ for line in inp:
+ #print (line)
+ outfile.write(line+'\n')
+ outfile.close()
+ print('File: %s was written' % fname)
+
+ # read in the calculix mesh
+ self.__read_inp(fname)
+
+"""This module stores geometry classes, which are used to make parts.
+"""
+from math import atan2, pi, cos, sin, radians, degrees
+from matplotlib.patches import Arc as AArc
+from matplotlib.patches import PathPatch
+from matplotlib.path import Path
+from numpy.linalg import det
+from numpy import array
+from numpy.linalg import solve as linsolve
+
+from . import base_classes
+
+#accuracy for small numbers in math below
+ACC = .00001 # original
+
+LDIVS = {}
+
+[docs]def get_text_hv(angle):
+ """Returns (ha, va) text horizontal and vertical alignment for line label.
+
+ This makes the line label text inside its area, assuming areas closed CW.
+
+ Args:
+ angle (float): line normal vector in degrees, between 180 and -180
+ """
+ horiz, vert = ('', '')
+ if abs(angle) < 67.5:
+ horiz = 'right'
+ elif 67.5 <= abs(angle) < 112.5:
+ horiz = 'center'
+ elif 112.5 <= abs(angle) <= 180:
+ horiz = 'left'
+ if abs(angle) < 22.5 or abs(angle) > 157.5:
+ vert = 'center'
+ else:
+ if angle > 0:
+ vert = 'top'
+ else:
+ vert = 'bottom'
+ return (horiz, vert)
+
+
+[docs]class Point(base_classes.Idobj):
+ """Makes a point or vector.
+
+ Args:
+ x (float): vertical coordinate of point
+ y (float): horizontal coordinate of point
+ z (float): in-page coordinate of point, defaults to zero
+
+ Attributes:
+ x (float): vertical coordinate of point
+ y (float): horizontal coordinate of point
+ z (float): in-page coordinate of point, defaults to zero
+ nodes (Node): a list of nodes in a mesh at this point.
+ List length is 1.
+ nodes (list): list of nodes under point
+ lines (set): set of lines which use this point
+ arc_center (bool): True if this point is an arc center
+ """
+
+ def __init__(self, x, y, z=0):
+ self.x = x
+ self.y = y
+ self.z = z
+ self.nodes = []
+ self.lines = set()
+ self.arc_center = False
+ self.on_part = True
+ self.esize = None
+ base_classes.Idobj.__init__(self)
+
+
+
+ def __eq__(self, other):
+ """Checks x and y equality to other point.
+
+ Args:
+ other (Point): passed point
+ """
+ if isinstance(other, Point):
+ if self.x == other.x and self.y == other.y:
+ return True
+ return False
+ else:
+ return False
+
+ def __add__(self, other):
+ """Returns new point = self + other.
+
+ Args:
+ other (Point): passed point
+
+ Returns:
+ Point = self + other (new point, does not modify self)
+ """
+ return Point(self.x+other.x, self.y+other.y, self.z+other.z)
+
+ def __sub__(self, other):
+ """Returns new point = self - other.
+
+ Args:
+ other (Point): passed point
+
+ Returns:
+ Point = self - other (new point, does not modify self)
+ """
+ return Point(self.x-other.x, self.y-other.y, self.z-other.z)
+
+ def __hash__(self):
+ """Set the hash of the item to the id number.
+
+ This allows one to make sets of these items.
+ """
+ return self.id
+
+ def __mul__(self, factor):
+ """Returns new point = self*factor.
+
+ Args:
+ factor (float): factor to multiply the Point by
+
+ Returns:
+ Point = self*factor (new point, does not modify self)
+ """
+ return Point(self.x*factor, self.y*factor, self.z*factor)
+
+ def __truediv__(self, other):
+ """Returns new point = self/other.
+
+ Args:
+ other (float or Point): item we're dividing by
+
+ Returns:
+ Point or factor. Does not modify self.
+ """
+ if isinstance(other, float):
+ return Point(self.x/other, self.y/other, self.z/other)
+ elif isinstance(other, Point):
+ res = []
+ if other.x != 0:
+ res.append(self.x/other.x)
+ if other.y != 0:
+ res.append(self.y/other.y)
+ if other.z != 0:
+ res.append(self.z/other.z)
+ factor = sum(res)/len(res)
+ return factor
+
+[docs] def set_line(self, line):
+ """Saves the line or arc in the point."""
+ if line not in self.lines:
+ self.lines.add(line)
+ if self.arc_center == True:
+ if len(self.lines) > 1:
+ for line in self.lines:
+ if isinstance(line, Line):
+ self.arc_center = False
+ break
+
+[docs] def unset_line(self, line):
+ """Removes the line or arc in the point."""
+ if line not in self.lines:
+ return
+ self.lines.remove(line)
+
+[docs] def length(self):
+ """Returns the length of this point or vector."""
+ res = (self.x**2 + self.y**2)**0.5
+ return res
+
+[docs] def label(self, axis):
+ """Labels the point on a Matplotlib axis.
+
+ Args:
+ axis (Matplotlib Axis): Matplotlib Axis
+ """
+ axis.annotate(self.get_name(), (self.y, self.x))
+
+[docs] def plot(self, axis, label=True):
+ """Plots the point on the passed matplotlib axis.
+
+ Args:
+ axis (Matplotlib axis): plate to plot the point
+ pnum (bool): True turns on point labeling
+ """
+ axis.scatter(self.y, self.x)
+ if label:
+ self.label(axis)
+
+[docs] def make_unit(self):
+ """Modifies self vector to make it a unit vector."""
+ lenval = self.length()
+ self.x = self.x*1/lenval
+ self.y = self.y*1/lenval
+
+[docs] def ang_rad(self):
+ """Returns angle in radians, assume point is vector from 0,0."""
+ vert = self.x
+ horiz = self.y
+ rad_angle = atan2(vert, horiz)
+ return rad_angle
+
+[docs] def ang_deg(self):
+ """Returns angle in degrees, assume point is vector from 0,0."""
+ rad_angle = self.ang_rad()
+ deg = rad_angle * 180.0 / pi
+ return deg
+
+[docs] def axrad(self):
+ """Returns a tuple of (ax, rad) for plotting."""
+ return (self.y, self.x)
+
+[docs] def radax(self):
+ """Returns a tuple of (rad, ax) for checking if inside."""
+ return (self.x, self.y)
+
+[docs] def rot_ccw_deg(self, ang):
+ """Rotates the current vector by ccw degrees about 0,0. Returns self.
+
+ Args:
+ ang (float): angle or rotation in degrees, counter-clockwise (ccw)
+ """
+ ang = radians(ang)
+ axial = self.y*cos(ang) - self.x*sin(ang)
+ rad = self.y*sin(ang) + self.x*cos(ang)
+ axial = round(axial, 5)
+ rad = round(rad, 5)
+ (self.x, self.y) = (rad, axial)
+ return self
+
+[docs] def ang_bet_rad(self, other):
+ """Returns angle between self and other vector in radians.
+
+ Assumes vectors start at 0,0
+
+ Args:
+ other (Point): passed vector
+
+ Returns:
+ angle (float): radians between self vector and other vector
+ """
+ avect = Point(self.y*other.x - self.x*other.y, self.y*other.y + self.x*other.x)
+ ang = avect.ang_rad()
+ return ang
+
+[docs] def ang_bet_deg(self, other):
+ """Returns angle between self and other vector in degrees.
+
+ Assumes vectors start at 0,0
+
+ Args:
+ other (Point): passed vector
+
+ Returns:
+ angle (float): degrees between self vector and other vector
+ """
+ avect = Point(self.y*other.x - self.x*other.y, self.y*other.y + self.x*other.x)
+ ang = avect.ang_deg()
+ return ang
+
+ def __str__(self):
+ """Returns string listing object type, id number, and coordinates"""
+ val = 'Point %s, (x, y)=(%.3f, %.3f)' % (self.get_name(), self.x, self.y)
+ return val
+
+
+
+
+
+[docs]class Line(base_classes.Idobj):
+ """Stores a line from points start_pt to end_pt.
+
+ Args:
+ start_pt (Point): first point
+ end_pt (Point): second point
+
+ Attributes:
+ points (list of Point) : list of the line's points [start_pt, end_pt]
+ ediv (None or float): number of elements on the line
+ midpt (Point): a mid-point between start_pt and end_pt
+ nodes (list of Node): a list of meshed nodes on the line
+ edge (bool): True if line is a non-shared edge in an area
+ signlines (list): list of SignLine that use this Line
+ """
+
+ def __init__(self, start_pt, end_pt):
+ self.points = [start_pt, end_pt]
+ self.midpt = self.mid()
+ self.ediv = None
+ self.nodes = []
+ self.edge = True
+ self.signlines = []
+ base_classes.Idobj.__init__(self)
+
+
+
+ @property
+ def allpoints(self):
+ """Returns all line defining points."""
+ return self.points
+
+ @property
+ def areas(self):
+ """Returns all areas that signlines use."""
+ areas = set()
+ for sline in self.signlines:
+ areas.add(sline.parent)
+ return list(areas)
+
+ def __eq__(self, other):
+ """Returns bool telling if self == other.
+
+ Returns True if the lines are the same type, and contain the same
+ list of points.
+ """
+ if isinstance(other, self.__class__):
+ mypts = self.points
+ otherpts = other.points
+ same = True
+ for (mypt, otherpt) in zip(mypts, otherpts):
+ if mypt != otherpt:
+ same = False
+ break
+ return same
+ else:
+ return False
+
+ def __hash__(self):
+ """Set the hash of the item to the id number.
+
+ This allows one to make sets of these items.
+ """
+ return self.id
+
+[docs] def add_signline(self, signline):
+ """Adds a signline to self.signlines."""
+ if signline not in self.signlines:
+ self.signlines.append(signline)
+
+
+
+[docs] def length(self):
+ """Returns the length of the line."""
+ return (self.pt(0)-self.pt(1)).length()
+
+[docs] def save_to_points(self):
+ """This stores this line in the line's child points."""
+ for point in self.allpoints:
+ point.set_line(self)
+
+[docs] def set_ediv(self, ediv):
+ """Sets the number of element divisions on the line when meshing.
+
+ Args:
+ ediv (int): number of required elements on this line
+ """
+ self.ediv = ediv
+
+[docs] def set_esize(self, esize):
+ """Sets the size of mesh elements on the line when meshing.
+
+ Args:
+ esize (float): size of mesh elements on this line
+ """
+ for i in range(2):
+ self.pt(i).set_esize(esize)
+
+[docs] def signed_copy(self, sign):
+ """Returns a SignLine copy of this Line with the passed sign."""
+ return SignLine(self, sign)
+
+[docs] def pt(self, ind):
+ """Returns the start or end point of the line.
+
+ Args:
+ ind (int): index of the point we want, 0=start, 1=end
+
+ Returns:
+ Point: requested Point
+ """
+ return self.points[ind]
+
+[docs] def set_pt(self, ind, point):
+ """Update the line's point at ind to the new point.
+
+ Args:
+ ind (int): index of the point we're updating, 0=start, 1=end
+ point (Point): new point we are assigning to Line
+ """
+ # remove this line from the old point
+ self.points[ind].lines.remove(self)
+ # update the point
+ self.points[ind] = point
+ self.midpt = self.mid()
+ self.save_to_points()
+ other_point = self.pt(int(not ind))
+ if point == other_point:
+ # we are making a zero length line
+ other_point.lines.remove(self)
+ if self.id != -1:
+ if len(self.signlines) > 0:
+ area = self.signlines[0].lineloop.parent
+ area.part.fea.lines.remove(self)
+ for sline in self.signlines:
+ sline.lineloop.remove(sline)
+
+[docs] def mid(self):
+ """Calculates and returns a midpoint between start and end points."""
+ point = (self.pt(0) + self.pt(1))*0.5
+ return point
+
+[docs] def touches(self, other):
+ """Checks if a line is connected to a passed line.
+
+ Args:
+ other (Line or Arc): passed other line
+
+ Returns:
+ True: beginning or end point connect to other line
+ False: lines are not connected
+ """
+ [start_pt, end_pt] = [self.pt(0), self.pt(1)]
+ if start_pt in other.points or end_pt in other.points:
+ return True
+ else:
+ return False
+
+[docs] def offset(self, dist):
+ """Returns a new line offset from the dirrent line by value dist.
+
+ Right hand rule is followed assuming x is the origina line.
+ Z points out of page.
+ Y is the direction of line offset.
+
+ Args:
+ dist (float): distance to offset the new line by
+
+ Returns:
+ tmpline (Line): new Line which is offset from self by dist
+ """
+ lvect = self.get_perp_vec()
+ lvect.make_unit()
+ # scale it by the new amount
+ lvect = lvect*dist
+ # add the vector onto the defining points
+ start_pt = self.pt(0) + lvect
+ end_pt = self.pt(1) + lvect
+ tmpline = Line(start_pt, end_pt)
+ return tmpline
+
+[docs] def get_abc(self):
+ """ Returns a list of abc terms for a line definition.
+
+ ax + by + c = 0
+
+ Returns:
+ [a, b, c]: list of the above line terms
+ """
+ delta = self.pt(1)-self.pt(0)
+ dx_ = delta.x
+ dy_ = delta.y
+ (a__, b__, c__) = (0.0, 0.0, 0.0)
+ if dy_ == 0:
+ # vertical radial line
+ # y = 5 --> 0 = -y + 5 --> 0 = 0x -y + 5
+ (b__, c__) = (-1.0, self.pt(0).y)
+ elif dx_ == 0:
+ # horizontal axial line
+ # x = 5 --> 0 = -x + 5 --> 0 = -x +0y + 5
+ (a__, c__) = (-1.0, self.pt(0).x)
+ else:
+ # y = ax + c --> 0 = ax - y + c
+ slope = dy_/dx_
+ offset = self.pt(0).y - slope*self.pt(0).x
+ (a__, b__, c__) = (slope, -1.0, offset)
+ return [a__, b__, c__]
+
+[docs] def get_perp_vec(self, point=None):
+ """ Returns vector perpendicular to current line.
+
+ Vector is created by rotating the current line ccw 90 degrees
+ around the start point.
+
+ Args:
+ point (Point): defaults to None, unused for Line, used for Arc
+
+ Returns:
+ lvect (Point): vector that is perpendiculat to the current line
+ """
+ lvect = self.pt(1) - self.pt(0)
+ lvect.rot_ccw_deg(90)
+ return lvect
+
+[docs] def get_tan_vec(self, point):
+ """ Returns vector tangent to the current line.
+
+ Vector points from the passed point to a unit location further away.
+
+ Args:
+ point (Point): start point
+
+ Returns:
+ lvect (Point): vector that is tangent to the current line
+ """
+ end = self.pt(0)
+ if point == end:
+ end = self.pt(1)
+ lvect = end - point
+ lvect.make_unit()
+ return lvect
+
+[docs] def arc_tang_intersection(self, point, mag):
+ """Returns an intersection point on this line of an arc at center point
+
+ Args:
+ point (Point): arc center point
+ mag (float): passed radius distance from the arc center to the line
+
+ Returns:
+ newpt (Point or None):
+ Point: the intersection point
+ None: returns None if no intersection point exists
+ """
+ vect = self.get_perp_vec()
+ vect.make_unit()
+ pt2 = point + vect*(-2*mag)
+
+ tmpline = Line(point, pt2)
+ newpt = tmpline.intersects(self)
+ if type(newpt) == type(None):
+ print('Intersection failed!')
+ return newpt
+
+[docs] def coincident(self, point):
+ """Checks to see if point is on the line.
+
+ Args:
+ point (Point): input point to check
+
+ Returns:
+ bool: True if point on line or False if not
+ """
+ # check if point is almost equal to end points
+ dist = point - self.pt(0), point - self.pt(1)
+ dist = [item.length() for item in dist]
+ if min(dist) < ACC:
+ return True
+ # check if point is on the lines
+ [a__, b__, c__] = self.get_abc()
+ remainder = abs(a__*point.x + b__*point.y + c__)
+ # print(' coincident: %s %s' % (self.get_name(), point))
+ # print(' remainder: %.6f' % remainder)
+ if remainder < ACC:
+ # point is on the line equation, but is it between end pts
+ xvals = [apoint.x for apoint in self.points]
+ yvals = [apoint.y for apoint in self.points]
+ xin = min(xvals)-ACC <= point.x <= max(xvals)+ACC
+ yin = min(yvals)-ACC <= point.y <= max(yvals)+ACC
+ if xin and yin:
+ # we're on the line
+ return True
+ else:
+ # we're off of the line
+ # print(' outside of bounds!')
+ return False
+ else:
+ # point is not on line equation
+ return False
+
+[docs] def intersects(self, other):
+ """Checks if other line intersects this line.
+
+ Args:
+ other (Line, Arc, Signline, SignArc): other line to check
+
+ Returns:
+ Point or None:
+ Point: intersection point
+ None: None if lines don't intersect
+ """
+ # other is cutline when checking
+ if isinstance(other, Line) or isinstance(other, SignLine):
+ terms_1 = self.get_abc()
+ terms_2 = other.get_abc()
+ p_1, p_2 = self.pt(0), self.pt(1)
+ p_3, p_4 = other.pt(0), other.pt(1)
+ dist = [p_1 - p_3, p_1 - p_4, p_2 - p_3, p_2 - p_4]
+ vals = [p_1, p_1, p_2, p_2]
+ dist = [item.length() for item in dist]
+ # check each of self's point to see if it == other's points
+ for (distance, val) in zip(dist, vals):
+ if distance < ACC:
+ return Point(val.x, val.y)
+ if terms_1[:2] == terms_2[:2]:
+ # lines are parallel and will not have an intersection point
+ #print('Line-Line intersection = None: same slope')
+ return None
+ denom = (p_1.x-p_2.x)*(p_3.y-p_4.y) - (p_1.y-p_2.y)*(p_3.x-p_4.x)
+ if abs(denom) < ACC:
+ # lines are nearly identical slope
+ #print('Line-Line intersection = None because of denom')
+ return None
+ else:
+ # intersection calculation
+ a_array = array([[terms_1[0], terms_1[1]],
+ [terms_2[0], terms_2[1]]])
+ b_array = array([-terms_1[2], -terms_2[2]])
+ result = linsolve(a_array, b_array)
+ point = Point(result[0], result[1])
+ # check against existing points on self
+ dist = [point - p_1, point - p_2]
+ vals = [p_1, p_2]
+ dist = [item.length() for item in dist]
+ for (distance, loc) in zip(dist, vals):
+ if distance < ACC:
+ point = Point(loc.x, loc.y)
+ break
+ on_self = self.coincident(point)
+ on_other = other.coincident(point)
+ if on_self and on_other:
+ """
+ s_name = self.get_name()
+ o_name = other.get_name()
+ print('X Checking intersection of %s and %s' % (s_name, o_name))
+ print(' %s' % point)
+ print(' (on_self, on_other): (%s, %s)' % (on_self, on_other))
+ """
+ return point
+ else:
+ """
+ print('Line-Line intersection = None point not on lines')
+ print(' %s' % point)
+ print(' on_self: %s' % on_self)
+ print(' on_other: %s' % on_other)
+ """
+ return None
+ elif isinstance(other, Arc) or isinstance(other, SignArc):
+ # arc line intersection
+ return other.intersects(self)
+
+ def __str__(self):
+ """Returns string listing object type, name, and points"""
+ start = self.pt(0)
+ end = self.pt(1)
+ val = 'Line %s, start: %s end: %s' % (self.get_name(), start, end)
+ return val
+
+[docs]class SignLine(Line, base_classes.Idobj):
+ """Makes a signed line.
+
+ Attributes:
+ line (Line): parent line
+ sign (int): 1 = positive, or -1 = negative
+ lineloop (None or LineLoop): parent LineLoop, default is None
+ faces (list): list of element faces
+ n1 (list): list of first order nodes on line
+ nodes (list): nodes on the line
+ """
+ def __init__(self, parent_line, sign):
+ self.line = parent_line
+ self.sign = sign
+ self.lineloop = None
+ self.faces = []
+ self.n1 = []
+ base_classes.Idobj.__init__(self)
+
+[docs] def pt(self, index):
+ """Returns the point at the start or end of the line.
+
+ Args:
+ index (int): 0 = start point, 1 = end point
+ """
+ if self.sign == -1:
+ return self.line.points[not index]
+ else:
+ return self.line.points[index]
+
+[docs] def set_pt(self, index, point):
+ """Sets the start or end point of the line to a new point.
+
+ Args:
+ index (int): 0 for start point, 1 for end point
+ point (Point): the new point we are updating to
+ """
+ self.line.set_pt(index, point)
+
+[docs] def set_ediv(self, ediv):
+ """Applies the element divisions onto the parent line."""
+ self.line.set_ediv(ediv)
+
+[docs] def set_esize(self, esize):
+ """Applies the element size onto the parent line."""
+ self.line.set_esize(esize)
+
+[docs] def set_lineloop(self, lineloop):
+ """Sets the parent LineLoop"""
+ # this is needed to cascade set ediv up to FEA model and down onto
+ # other lines with this ediv
+ self.lineloop = lineloop
+
+[docs] def get_name(self):
+ """Returns the SignLine name."""
+ if self.sign == -1:
+ return '-L'+str(self.line.id)
+ else:
+ return 'L'+str(self.line.id)
+
+[docs] def label(self, ax):
+ """Labels the line on the matplotlib ax axis.
+
+ Args:
+ ax (Matplotlib axis): Matplotlib axis
+ """
+ # find angle of perpendiclar line so we can set text alignment
+ lvect = self.get_perp_vec()
+ ang = lvect.ang_deg()
+ horiz, vert = get_text_hv(ang)
+ axial = self.midpt.y
+ radial = self.midpt.x
+ ax.text(axial, radial, self.get_name(), ha=horiz, va=vert)
+
+[docs] def plot(self, ax, label=True):
+ """Draws and labels the line onto the passed Matplotlib axis.
+
+ Args:
+ ax (Matlotlib axis): Matplotlib axis
+ label (bool): if True label the line
+ """
+ lax = [pt.y for pt in self.points]
+ lrad = [pt.x for pt in self.points]
+ # draw line
+ ax.plot(lax, lrad)
+ if label:
+ self.label(ax)
+
+[docs] def signed_copy(self, sign):
+ """Returns a SignLine copy of this line."""
+ realsign = self.sign*sign
+ return SignLine(self.line, realsign)
+
+ @property
+ def nodes(self):
+ """Gets the line nodes."""
+ return self.line.nodes
+
+ @nodes.setter
+ def nodes(self, nlist):
+ """Sets the line nodes."""
+ self.line.nodes = nlist
+
+ @property
+ def midpt(self):
+ """Returns the mid point."""
+ return self.line.midpt
+
+ @property
+ def edge(self):
+ """Returns a bool saying if this is an edge line."""
+ return self.line.edge
+
+ @edge.setter
+ def edge(self, isedge):
+ """Sets the line nodes."""
+ self.line.edge = isedge
+
+ @property
+ def points(self):
+ """Returns a list of the line's points [pstart, pend]"""
+ pstart = self.pt(0)
+ pend = self.pt(1)
+ return [pstart, pend]
+
+
+
+
+[docs]class Arc(base_classes.Idobj):
+ """Makes an arc from points start_pt to end_pt about center actr.
+
+ Arcs should be 90 degrees maximum.
+
+ Arguments:
+ start_pt (Point): first point
+ end_pt (Point): second point
+ actr (Point): arc center
+
+ Attributes:
+ points (list of Point) : list of the line's points [start_pt, end_pt]
+ allpoints (list of Point): [start_pt, end_pt, actr]
+ actr (Point): arc center
+ radius (float): radius of the arc
+ concavity (str):
+
+ - 'concave': concave arc
+ - 'convex': convex arc
+
+ midpt (Point): a mid-point between p1 and p2
+ ediv (None or int): number of elements on arc
+ nodes (list of Node): a list of meshed nodes on the line
+ edge (bool): True if arc is a non-shared edge in an area
+ signlines (list): list of SignArc that use this Arc
+ """
+
+ def __init__(self, start_pt, end_pt, actr):
+ self.points = [start_pt, end_pt]
+ self.actr = actr
+ self.actr.arc_center = True
+ self.radius = (start_pt - actr).length()
+ self.concavity = self.get_concavity()
+ self.midpt = self.mid()
+ self.ediv = None
+ self.nodes = []
+ self.edge = True
+ self.signlines = []
+ base_classes.Idobj.__init__(self)
+
+[docs] def reverse(self):
+ """Reverses self."""
+ self.points.reverse()
+ self.concavity = self.get_concavity()
+
+ @property
+ def allpoints(self):
+ """Returns all points including arc center"""
+ return self.points + [self.actr]
+
+ @property
+ def areas(self):
+ """Returns all areas that signlines use."""
+ areas = set()
+ for sline in self.signlines:
+ areas.add(sline.parent)
+ return list(areas)
+
+ def __eq__(self, other):
+ """Returns bool telling if self == other
+
+ Returns True if the lines are the same type, and contain the same
+ list of points.
+ """
+ if isinstance(other, self.__class__):
+ # check to see if the end points are the same
+ mypts = self.points
+ otherpts = other.points
+ same = True
+ for (mypt, otherpt) in zip(mypts, otherpts):
+ if mypt != otherpt:
+ same = False
+ break
+ # flip to false if arc centers are different
+ myactr = self.actr
+ otheractr = other.actr
+ if myactr != otheractr:
+ same = False
+ return same
+ else:
+ return False
+
+ def __hash__(self):
+ """Set the hash of the item to the id number.
+
+ This allows one to make sets of these items.
+ """
+ return self.id
+
+[docs] def add_signline(self, signline):
+ """Adds a signline to self.signlines."""
+ if signline not in self.signlines:
+ self.signlines.append(signline)
+
+
+
+[docs] def length(self):
+ """Returns the length of the arc."""
+ a = self.pt(0)-self.actr
+ b = self.pt(1)-self.actr
+ ang = abs(a.ang_bet_rad(b))
+ res = self.radius*ang
+ return res
+
+[docs] def save_to_points(self):
+ """This stores this line in the line's child points."""
+ for point in self.allpoints:
+ point.set_line(self)
+
+[docs] def set_ediv(self, ediv):
+ """Sets the number of element divisions on the arc when meshing.
+
+ Args:
+ ediv (int): number of required elements on this arc
+ """
+ self.ediv = ediv
+
+[docs] def set_esize(self, esize):
+ """Sets the size of mesh elements on the arc when meshing.
+
+ Args:
+ esize (float): size of mesh elements on this arc
+ """
+
+ for i in range(2):
+ self.pt(i).set_esize(esize)
+
+[docs] def signed_copy(self, sign):
+ """Returns a SignArc instance of this Arc with the passed sign."""
+ return SignArc(self, sign)
+
+[docs] def pt(self, ind):
+ """Returns the start or end point of the arc.
+
+ Args:
+ ind (int): index of the point we want, 0=start, 1=end
+
+ Returns:
+ Point: requested Point
+ """
+ return self.points[ind]
+
+[docs] def set_pt(self, ind, point):
+ """Update the arc's point at ind to the new point.
+
+ Args:
+ ind (int): index of the point we're updating, 0=start, 1=end
+ point (Point): new point we are assigning to Line
+ """
+ # remove this line from the old point
+ self.points[ind].lines.remove(self)
+ # set the new point
+ self.points[ind] = point
+ self.midpt = self.mid()
+ self.save_to_points()
+ other_point = self.pt(int(not ind))
+ if point == other_point:
+ other_point.lines.remove(self)
+ for sline in self.signlines:
+ sline.lineloop.remove(sline)
+ if self.id != -1:
+ if len(self.signlines) > 0:
+ area = self.signlines[0].lineloop.parent
+ area.part.fea.lines.remove(self)
+
+[docs] def mid(self):
+ """Caclulates and returns the midpoint of the arc."""
+ pt = self.get_pt_at(0.5)
+ return pt
+
+[docs] def touches(self, other):
+ """Checks if this arc is connected to a passed line.
+
+ Args:
+ other (Line or Arc): passed other line
+
+ Returns:
+ True: beginning or end point connect to other line
+ False: lines are not connected
+ """
+ [p1, p2] = [self.pt(0), self.pt(1)]
+ if p1 in other.points or p2 in other.points:
+ return True
+ else:
+ return False
+
+[docs] @staticmethod
+ def sgn(val):
+ """Returns sign of the passed value.
+
+ Helper function for self.interects code.
+
+ Args:
+ val (float): passed value
+
+ Returns:
+ 1 or -1: sign of the passed value, assumes sign of zero = 1
+ """
+ if val < 0:
+ return -1
+ else:
+ return 1
+
+[docs] def get_ang(self):
+ """ Returns angle from beginning to end of arc from arc ctr in degrees
+
+ Notes:
+ Answer is between 180 and -180
+ Positive is concave
+ Negative is convex
+ """
+ a = self.pt(0)-self.actr
+ b = self.pt(1)-self.actr
+ ang = a.ang_bet_deg(b)
+ return ang
+
+[docs] def get_ang_rad(self):
+ """ Returns angle from beginning to end of arc from arc ctr in rad
+
+ Notes:
+ Answer is between pi and -pi
+ Positive is concave
+ Negative is convex
+ """
+ a = self.pt(0)-self.actr
+ b = self.pt(1)-self.actr
+ ang = a.ang_bet_rad(b)
+ return ang
+
+[docs] def get_concavity(self, clockwise=True):
+ """Returns concavity of this arc, 'concave' or 'convex'."""
+ ang = self.get_ang()
+ res = ''
+ if clockwise:
+ res = 'concave'
+ if ang < 0:
+ res = 'convex'
+ else:
+ res = 'convex'
+ if ang < 0:
+ res = 'concave'
+ return res
+
+[docs] def get_pt_at(self, nondim):
+ """Returns point on arc between start and end point at nondim dist.
+
+ Args:
+ nondim (float): distance between start and end point, 0=start, 1=end
+
+ Returns:
+ res (Point): requested point on arc
+ """
+ a = self.pt(0)-self.actr
+
+ # neg is convex, pos is concave
+ ang = self.get_ang()
+
+ ang_delta = nondim*ang
+ a.rot_ccw_deg(ang_delta)
+
+ res = self.actr + a
+ return res
+
+[docs] def get_perp_vec(self, pt):
+ """ Returns vector perpendicular to current arc.
+
+ Right hand rule is used:
+ X+ is the direction of the arc from start to end.
+ Z+ is out of the page
+ Y+ is the direction of the perpendicular vector
+
+ Args:
+ pt (Point): point on arc where we want the perpendicular vector
+
+ Returns:
+ resv (Point): perpendicular vector
+ """
+ resv = None
+ if self.concavity == 'convex':
+ resv = pt - self.actr
+ elif self.concavity == 'concave':
+ resv = self.actr - pt
+ return resv
+
+[docs] def get_tan_vec(self, point):
+ """ Returns vector tangent to the current arc.
+
+ Vector points from the passed point to a unit location further away.
+
+ Args:
+ point (Point): start point
+
+ Returns:
+ resv (Point): vector that is tangent to the current line
+ """
+ if self.concavity == 'concave':
+ resv = self.actr - point
+ elif self.concavity == 'convex':
+ resv = point - self.actr
+ if point == self.pt(0):
+ resv.rot_ccw_deg(-90)
+ else:
+ resv.rot_ccw_deg(90)
+ resv.make_unit()
+ return resv
+
+[docs] def get_verts_codes(self, plot=True):
+ """Returns a list of [verts, codes] vertices for plotting."""
+ # http://www.spaceroots.org/documents/ellipse/node19.html
+ # http://www.tinaja.com/glib/bezarc1.pdf
+ # http://stackoverflow.com/questions/734076/geometrical-arc-to-bezier-curve
+
+ # find the angle of rotation to the middle of the arc
+ ang_offset = self.get_pt_at(0.5) - self.actr
+ ang_offset = ang_offset.ang_deg()
+
+ # store the half angle as theta
+ theta = radians(abs(self.get_ang()))*0.5
+
+ ax1 = (4 - cos(theta))/3
+ rad1 = ((1 - cos(theta))*(cos(theta) - 3))/(3*sin(theta))
+ pt1 = Point(rad1, ax1)*self.radius
+ pt2 = Point(-rad1, ax1)*self.radius
+
+ # fix the angles
+ pt1.rot_ccw_deg(ang_offset)
+ pt2.rot_ccw_deg(ang_offset)
+
+ # fix the location
+ pt1 += self.actr
+ pt2 += self.actr
+
+ # need to change the order if the arc is concave or convex
+ verts = [pt1, pt2, self.pt(1)]
+ if self.concavity == 'convex':
+ verts = [pt2, pt1, self.pt(1)]
+
+ if plot:
+ verts = [point.axrad() for point in verts]
+ else:
+ verts = [point.radax() for point in verts]
+ codes = [Path.CURVE4, Path.CURVE4, Path.CURVE4]
+ return [verts, codes]
+
+[docs] def coincident(self, pt):
+ """Checks to see if pt is on the arc.
+
+ Args:
+ pt (Point): input point to check
+
+ Returns:
+ bool: True if pt on arc or False if not
+ """
+ # (x - xc)^2 + (y - yc)^2 - r^2 = 0
+ dist = pt - self.pt(0), pt - self.pt(1)
+ dist = [item.length() for item in dist]
+ if min(dist) < ACC:
+ # True if point almost equal to end points
+ return True
+ remainder = abs((pt.x - self.actr.x)**2 + (pt.y - self.actr.y)**2
+ - self.radius**2)
+ if remainder < ACC:
+ # point is on the circle, but is it in the arc?
+ # check if the angle traversed to the point is within the ang
+ # traversed by the arc
+ a1 = self.pt(0)-self.actr
+ a2 = pt-self.actr
+ ang_pt = a1.ang_bet_deg(a2)
+
+ # this is the angle traversed byt the arc
+ ang = self.get_ang()
+ minval = min(0, ang)
+ maxval = max(0, ang)
+ if minval <= ang_pt <= maxval:
+ # point is on our arc
+ return True
+ else:
+ # point on circle but outside our arc
+ return False
+ else:
+ # point not on arc
+ return False
+
+[docs] def intersects(self, other):
+ """Checks if other line intersects this arc.
+
+ Args:
+ other (Line or Arc): other line to check
+
+ Returns:
+ Point or None:
+ Point: intersection point
+ None: None if lines don't intersect
+ """
+ if isinstance(other, Line) or isinstance(other, SignLine):
+ #print('|C Checking arc-line intersection')
+ #print(' %s' % self)
+ # arc-line intersection
+ # reset the math to be about the circle centroid
+ # http://mathworld.wolfram.com/Circle-LineIntersection.html
+ pt1 = other.pt(0) - self.actr
+ pt2 = other.pt(1) - self.actr
+ # for us x and y are reversed because of y axial
+ (x1, y1, x2, y2) = (pt1.y, pt1.x, pt2.y, pt2.x)
+ dx = x2 - x1
+ dy = y2 - y1
+ rad = self.radius
+ dr = (dx**2 + dy**2)**(0.5)
+ D = x1*y2 - x2*y1
+ dist = [self.pt(0) - other.pt(0), self.pt(0) - other.pt(1),
+ self.pt(1) - other.pt(0), self.pt(1) - other.pt(1)]
+ vals = [self.pt(0), self.pt(0), self.pt(1), self.pt(1)]
+ dist = [item.length() for item in dist]
+ for (distance, val) in zip(dist, vals):
+ if distance < ACC:
+ return Point(val.x, val.y)
+ discr = 1.0*(rad**2)*(dr**2) - D**2
+ if discr < 0:
+ # no intersection
+ #print('none returned because of discriminant')
+ return None
+ else:
+ # tangent or intersection
+ rsq = rad**2
+ drsq = dr**2
+ x1 = (D*dy + self.sgn(dy)*dx*(rsq*drsq-D**2)**(0.5))/drsq
+ y1 = (-D*dx + abs(dy)*(rsq*(dr**2)-D**2)**(0.5))/drsq
+ x2 = (D*dy - self.sgn(dy)*dx*(rsq*drsq-D**2)**(0.5))/drsq
+ y2 = (-D*dx - abs(dy)*(rsq*drsq-D**2)**(0.5))/drsq
+
+ # convert result points back to global coordinates
+ pt1 = Point(y1, x1) + self.actr
+ pt2 = Point(y2, x2) + self.actr
+ res = []
+ # check that the resultant points are on line and arcs
+ if self.coincident(pt1) and other.coincident(pt1):
+ for apoint in self.points:
+ dist = pt1 - apoint
+ dist = dist.length()
+ if dist < ACC:
+ pt1 = Point(apoint.x, apoint.y)
+ break
+ res.append(pt1)
+ if pt1 != pt2 and self.coincident(pt2) and other.coincident(pt2):
+ for apoint in self.points:
+ dist = pt2 - apoint
+ dist = dist.length()
+ if dist < ACC:
+ pt2 = Point(apoint.x, apoint.y)
+ break
+ res.append(pt2)
+ if len(res) == 1:
+ return res[0]
+ elif len(res) > 1:
+ return res
+ elif len(res) == 0:
+ #print('none returned because no point from math')
+ return None
+ elif isinstance(other, Arc):
+ # arc-arc intersection
+ # I have not yet written the code here
+ return None
+
+ def __str__(self):
+ """Returns string listing object type, name, and points"""
+ name = self.get_name()
+ start, end = self.pt(0), self.pt(1)
+ actr = self.actr
+ val = 'Arc %s, start: %s end: %s center: %s' % (name, start, end, actr)
+ return val
+
+[docs]class SignArc(Arc, base_classes.Idobj):
+ """Makes a signed arc.
+
+ Attributes:
+ line (Arc): parent arc
+ sign (int): 1 = positive, or -1 = negative
+ lineloop (None or LineLoop): parent LineLoop, default is None
+ nodes (list): nodes on the arc
+ faces (list): list of element faces
+ n1 (list): list of first order nodes on line
+ """
+
+ def __init__(self, parent_line, sign):
+ self.line = parent_line
+ self.sign = sign
+ self.lineloop = None
+ self.faces = []
+ self.n1 = []
+ base_classes.Idobj.__init__(self)
+
+[docs] def pt(self, ind):
+ """Returns the point at the start or end of the line."""
+ if self.sign == -1:
+ return self.line.points[not ind]
+ else:
+ return self.line.points[ind]
+
+[docs] def set_pt(self, ind, point):
+ """Sets the passed point to at the index location."""
+ self.line.set_pt(ind, point)
+
+[docs] def set_ediv(self, ediv):
+ """Apply the element divisions onto the parent line."""
+ self.line.set_ediv(ediv)
+
+[docs] def set_esize(self, esize):
+ """Applies the element size onto the parent line."""
+ self.line.set_esize(esize)
+
+[docs] def set_lineloop(self, lineloop):
+ """Sets the parent LineLoop"""
+ # this is needed to cascade set ediv up to FEA model and down onto
+ # other lines with this ediv
+ self.lineloop = lineloop
+
+[docs] def get_name(self):
+ """Returns the line name."""
+ if self.sign == -1:
+ return '-L'+str(self.line.id)
+ else:
+ return 'L'+str(self.line.id)
+
+[docs] def label(self, ax):
+ """Labels the arc on the matplotlib ax axis.
+
+ Args:
+ ax (Matplotlib axis): Matplotlib axis
+ """
+ lvect = self.get_perp_vec(self.midpt)
+ ang = lvect.ang_deg()
+ horiz, vert = get_text_hv(ang)
+ (axial, radial) = (self.midpt.y, self.midpt.x)
+ ax.text(axial, radial, self.get_name(), ha=horiz, va=vert)
+
+[docs] def plot(self, ax, label=True):
+ """Draws and labels the arc onto the passed Matplotlib axis.
+
+ Args:
+ ax: Matplotlib axis
+ label (bool): if True, label the arc
+ """
+ ctr = self.actr
+ vect1 = None
+ sign = 1
+ if self.concavity == 'concave':
+ vect1 = self.pt(0)-ctr
+ else:
+ sign = -1
+ vect1 = self.pt(1)-ctr
+ rad = self.radius
+ ang1 = (vect1.ang_deg())
+ ang2 = ang1 + sign*self.get_ang()
+
+ # matplotlib assumes ccw arc drawing, calculix assumes cw drawing
+ a = AArc(xy=[ctr.y, ctr.x], width=2*rad, height=2*rad, angle=0,
+ theta1=ang1, theta2=ang2)
+ ax.add_artist(a)
+ if label:
+ self.label(ax)
+
+[docs] def signed_copy(self, sign):
+ """Returns a SignArc copy of this arc."""
+ realsign = self.sign*sign
+ return SignArc(self.line, realsign)
+
+ @property
+ def nodes(self):
+ """Gets the line nodes."""
+ return self.line.nodes
+
+ @nodes.setter
+ def nodes(self, nlist):
+ """Sets the line nodes."""
+ self.line.nodes = nlist
+
+ @property
+ def actr(self):
+ """Return the arc center."""
+ return self.line.actr
+
+ @property
+ def radius(self):
+ """Return the arc radius."""
+ return self.line.radius
+
+ @property
+ def midpt(self):
+ """Return the mid point."""
+ return self.line.midpt
+
+ @property
+ def edge(self):
+ """Returns a bool saying if this is an edge line."""
+ return self.line.edge
+
+ @edge.setter
+ def edge(self, isedge):
+ """Sets the line nodes."""
+ self.line.edge = isedge
+
+ @property
+ def points(self):
+ """Return a list of the arc's points [pstart, pend]"""
+ pstart = self.pt(0)
+ pend = self.pt(1)
+ return [pstart, pend]
+
+ @property
+ def concavity(self):
+ """Return the arc's concavity 'concave' or 'convex'"""
+ if self.sign == 1:
+ return self.line.concavity
+ else:
+ if self.line.concavity == 'concave':
+ return 'convex'
+ else:
+ return 'concave'
+
+
+
+
+[docs]class LineLoop(list):
+ """Makes a LineLoop, which stores multiple Line/Arc or SignLine/SignArc.
+
+ Args:
+ items (list): Line or Arc or SignLine or SignArc
+ hole_bool (bool): True if hole, False otherwise
+ parent (Area or None): parent Area, defaults to None
+
+ Attributes:
+ closed (bool): True if closed
+ area (float): the loop's area, can be pos or neg. Pos = cw, neg = ccw
+ hole (bool): True if loop is a hole
+ ccw (bool): True if loop is counter-clockwise
+ center (Point): a point at the loop centroid
+ parent (Area or None): the loop's parent Area
+ """
+
+ def __init__(self, items=[], hole_bool=False, parent=None):
+ super().__init__(items)
+ self.id = -1
+ for item in self:
+ if isinstance(item, SignLine) or isinstance(item, SignArc):
+ item.lineloop = self
+ if item not in item.line.signlines:
+ item.line.add_signline(item)
+ self.hole = hole_bool
+ self.parent = parent
+
+
+
+ def __hash__(self):
+ """Sets the hash of the item to the id number.
+
+ This allows one to make sets of these items.
+ """
+ return self.id
+
+ @property
+ def closed(self):
+ """Returns True if Closed, False if not."""
+ if len(self) < 2:
+ return False
+ else:
+ for ind, curr_line in enumerate(self):
+ lower_ind = ind-1
+ prev_line = self[lower_ind]
+ if curr_line.pt(0) != prev_line.pt(1):
+ return False
+ return True
+
+ @property
+ def center(self):
+ """Returns None or a Point at the line loop's center."""
+ if self.closed == False:
+ return None
+ else:
+ cx_sum = 0
+ cy_sum = 0
+ a_sum = 0
+ arcs = []
+ for item in self:
+ if isinstance(item, Line) or isinstance(item, SignLine):
+ start, end = item.pt(0), item.pt(1)
+ det_val = det([[start.x, end.x], [start.y, end.y]])
+ cx_sum += (start.x + end.x)*det_val
+ cy_sum += (start.y + end.y)*det_val
+ a_sum += det_val
+ elif isinstance(item, Arc) or isinstance(item, SignArc):
+ # concave: draw to actr, adder is neg
+ # convex: draw to actr, adder is pos
+ tmp_list = [[item.pt(0), item.actr],
+ [item.actr, item.pt(1)]]
+ for [start, end] in tmp_list:
+ det_val = det([[start.x, end.x], [start.y, end.y]])
+ cx_sum += (start.x + end.x)*det_val
+ cy_sum += (start.y + end.y)*det_val
+ a_sum += det_val
+ arcs.append(item)
+ area_val = a_sum*0.5
+ val_list = []
+ if area_val != 0:
+ cx_val = (1/(6*area_val))*cx_sum
+ cy_val = (1/(6*area_val))*cy_sum
+ val_list.append([area_val, cx_val, cy_val])
+ for arc in arcs:
+ ang = arc.get_ang_rad()
+ area_val = -0.5*ang*(arc.radius**2)
+ half_ang = abs(ang)*0.5
+ vect = arc.midpt - arc.actr
+ vect.make_unit()
+ mag = 2*arc.radius*sin(half_ang)/(3*half_ang)
+ vect = vect*mag
+ center = arc.actr + vect
+ val_list.append([area_val, center.x, center.y])
+ a_sum = sum([aval for [aval, cx, cy] in val_list])
+ cxa_sum = sum([cx*aval for [aval, cx, cy] in val_list])
+ cya_sum = sum([cy*aval for [aval, cx, cy] in val_list])
+ cx_val = cxa_sum/a_sum
+ cy_val = cya_sum/a_sum
+ return Point(cx_val, cy_val)
+
+ @property
+ def ccw(self):
+ """Returns bool telling if the area is closed and clockwise."""
+ if self.closed == False:
+ return False
+ else:
+ area_val = self.area
+ if area_val < 0:
+ return True
+ else:
+ return False
+
+ @property
+ def area(self):
+ """Returns the loop's signed area.
+
+ Returns 0 if not closed.
+ Pos = cw, neg = ccw
+
+ Formula from: http://mathworld.wolfram.com/PolygonArea.html
+ """
+ res = 0
+ adder = 0
+ arcs = []
+ if self.closed == False:
+ return res
+ else:
+ for item in self:
+ if isinstance(item, Line) or isinstance(item, SignLine):
+ start, end = item.pt(0), item.pt(1)
+ det_matrix = [[start.x, end.x], [start.y, end.y]]
+ res += det(det_matrix)
+ elif isinstance(item, Arc) or isinstance(item, SignArc):
+ # concave: draw to actr, adder is neg
+ # convex: draw to actr, adder is pos
+ start, end = item.pt(0), item.actr
+ det_matrix = [[start.x, end.x], [start.y, end.y]]
+ res += det(det_matrix)
+ start, end = item.actr, item.pt(1)
+ det_matrix = [[start.x, end.x], [start.y, end.y]]
+ res += det(det_matrix)
+ arcs.append(item)
+ res = 0.5*res
+ # need to loop through arcs here, now that we know if cw or ccw
+ for arc in arcs:
+ ang = arc.get_ang_rad()
+ tmp_adder = -0.5*ang*(arc.radius**2)
+ adder += tmp_adder
+ res = res + adder
+ return res
+
+[docs] def reverse(self):
+ """Reverse the loop direction cw<-->ccw
+
+ This modifies the order of lines, and flips each line.
+ """
+ for line in self:
+ line.reverse()
+ super().reverse()
+
+[docs] def set_parent(self, parent):
+ """Sets the line loop's parent part."""
+ self.parent = parent
+
+[docs] def get_patch(self):
+ """Returns a patch of the LineLoop, or None of not closed."""
+ if self.closed:
+ codes, verts = [], []
+ # start poly
+ codes += [Path.MOVETO]
+ verts += [self[0].pt(0).radax()]
+ for line in self:
+ if isinstance(line, SignLine) or isinstance(line, Line):
+ verts.append(line.pt(1).radax())
+ codes.append(Path.LINETO)
+ elif isinstance(line, SignArc) or isinstance(line, Arc):
+ lverts, lcodes = line.get_verts_codes(plot=False)
+ verts += lverts
+ codes += lcodes
+ if line == self[-1]:
+ # close poly if at last line
+ verts.append((0, 0))
+ codes.append(Path.CLOSEPOLY)
+ path = Path(verts, codes)
+ patch = PathPatch(path, facecolor='yellow', linewidth=0.0)
+ return patch
+ return None
+
+[docs] def contains_point(self, point):
+ """Returns bool telling if pt is inside the LineLoop.
+
+ Args:
+ pt (Point): point we are checking
+
+ Returns:
+ contains (bool): True if in this area False if not
+ """
+ contains = False
+ if self.closed:
+ patch = self.get_patch()
+ contains = patch.contains_point(point.radax())
+ return contains
+
+[docs] def inside(self, other):
+ """Checks to see if self is inside other.
+
+ Only checks self's points.
+ """
+ if self.closed == False:
+ return False
+ if isinstance(other, LineLoop):
+ points = set()
+ for sline in self:
+ points.update(sline.points)
+ patch = other.get_patch()
+ for point in points:
+ contains = patch.contains_point(point.radax())
+ if contains == False:
+ return False
+ return True
+ else:
+ return False
+
+[docs] def append(self, item):
+ """Adds an item to this LineLoop
+
+ Args:
+ item (SingArc or SignLine): item to add to the list
+ """
+ if isinstance(item, SignLine) or isinstance(item, SignArc):
+ item.set_lineloop(self)
+ super().append(item)
+ return item
+
+[docs] def remove(self, item):
+ """Removes an SignLine or SignArc item from this LineLoop
+
+ Args:
+ item (SignArc or SignLine): item to remove
+ """
+ if isinstance(item, SignLine) or isinstance(item, SignArc):
+ item.set_lineloop(None)
+ if item not in item.line.signlines:
+ print('''===THIS SLINE NOT IN SLINE.LINE.SIGNLINES!===''')
+ print('Removing sline %s from loop in area %s'
+ % (item.get_name(), self.parent.get_name()))
+ print('SLINE.LINE %s' % item.line)
+ print('SLINE.LINE.SIGNLINES %r' % item.line.signlines)
+ super().remove(item)
+
+[docs] def insert(self, index, item):
+ """Inserts a SignLine or SignArc item into this LineLoop
+
+ Args:
+ item (SingArc or SignLine): item to insert
+ """
+ if isinstance(item, SignLine) or isinstance(item, SignArc):
+ item.set_lineloop(self)
+ super().insert(index, item)
+
+ def __str__(self):
+ """Prints the LineLoop definition."""
+ res = 'Number of items %i\n' % len(self)
+ for line in self:
+ res += '%s\n' % line
+ return res
+
+[docs]class Area(base_classes.Idobj):
+ """ Makes an area.
+
+ Area is closed in a clockwise direction.
+
+ Args:
+ part (Part): parent Part
+ line_list (list): list of Lines and Arcs that close the area
+
+ Attributes:
+ part (Part): parent part
+ closed (boolean): True if closed, False if not
+ exlines (LineLoop): LineLoop of signlines that define the area exterior
+ signlines (list): list of signed lines or arcs that define the area
+ lines (list): a list of all the lines that make the area, includes hole
+ lines
+ points (list): a list of all points making the area, excludes arc centers
+ allpoints (list): a list of all points, includes arc centers
+ holepoints (list): a list of hole points, excludes arc centers
+ matl (Matl): the material fo the area
+ etype (str): element type of area. Options are:
+ 'plstress': plane stress
+ 'plstrain': plane strain
+ 'axisym': axisymmetric
+ area (double): in-plane part area (ballpark number, arcs not calc right)
+ holes (list): list of LineLoop of signlines
+ center (Point or None): the area centroid point. None before closed.
+ nodes (list): child mesh nodes that are in the area
+ elements (list): child elements that are in the area
+ """
+
+ def __init__(self, part, line_list=[]):
+ base_classes.Idobj.__init__(self)
+ self.part = part # parent part reference
+ self.closed = False
+ self.matl = None
+ self.etype = None
+ self.area = 0.0
+ self.holes = []
+ self.center = None
+ self.nodes = []
+ self.elements = []
+ exloop = LineLoop(line_list, False, self)
+ self.part.fea.register(exloop)
+ self.exlines = exloop
+ if self.exlines.closed == True:
+ self.close()
+
+ def __hash__(self):
+ """Returns the item's id as its hash."""
+ return self.id
+
+ @property
+ def lines(self):
+ """Returns a list of all Line and Arc under this area."""
+ mylist = []
+ loops = [self.exlines] + self.holes
+ for loop in loops:
+ line_list = [line.line for line in loop]
+ mylist += line_list
+ return mylist
+
+ @property
+ def signlines(self):
+ """Returns a list of all SignLine and SignArc under this area."""
+ mylist = []
+ loops = [self.exlines] + self.holes
+ for loop in loops:
+ mylist += loop
+ return mylist
+
+ @property
+ def points(self):
+ """Returns a list of all points under this area; excludes arc centers."""
+ pts = set()
+ for line in self.lines:
+ pts.update(line.points)
+ return list(pts)
+
+ @property
+ def allpoints(self):
+ """Returns a list of all points under this area; includes arc centers."""
+ pts = set()
+ for line in self.lines:
+ pts.update(line.allpoints)
+ return list(pts)
+
+ @property
+ def holepoints(self):
+ """Returns a list of all points in this area's holes."""
+ pts = set()
+ for hole in self.holes:
+ for sline in hole:
+ pts.update(sline.points)
+ return list(pts)
+
+ def __str__(self):
+ """Returns a string listing the area's name."""
+ return self.get_name()
+
+
+
+[docs] def close(self):
+ """Closes the area."""
+ # sets closed=True
+ # calculates the area center
+ self.closed = True
+ if self.exlines.ccw == True:
+ self.exlines.reverse()
+ self.area, self.center = self.calc_area_center()
+
+[docs] def line_from_startpt(self, pt):
+ """Returns the signline in the area that start at the passed pt.
+
+ Args:
+ pt (Point): the start point of the line one wants
+
+ Returns:
+ match (None or SignLine or SignArc):
+ | Returns None if no line is found, otherwise the correct line or
+ | arc which starts with pt will be
+ | returned.
+ """
+ # returns the line that starts on the given point
+ matches = []
+ area_slines = self.signlines
+ for sline in area_slines:
+ if sline.pt(0) == pt:
+ matches.append(sline)
+ if len(matches) == 1:
+ return matches[0]
+ elif len(matches) > 1:
+ # we have a saw cut in the area and we're cutting on one of its pts
+ print('>1 slines starting with %s' % pt.get_name())
+ for match in matches:
+ slines = set(match.line.signlines)
+ if len(slines) == 2 and slines.issubset(area_slines):
+ # this is the saw cut one
+ pass
+ else:
+ # we want the non saw cut one
+ return match
+ return None
+
+[docs] def calc_area_center(self):
+ """Calculates and returns the area and centroid Point.
+
+ Returns:
+ list: [area, Point]
+ """
+ loops = [self.exlines] + self.holes
+ val_list = []
+ for loop in loops:
+ if loop.closed == True:
+ aval, cval = loop.area, loop.center
+ val_list.append([aval, cval])
+ a_sum = sum([aval[0] for aval in val_list])
+ cxa_sum = sum([center.x*aval for [aval, center] in val_list])
+ cya_sum = sum([center.y*aval for [aval, center] in val_list])
+ cx_val = cxa_sum/a_sum
+ cy_val = cya_sum/a_sum
+ center = Point(cx_val, cy_val)
+ return [a_sum, center]
+
+[docs] def get_maxlength(self):
+ """Returns max distance between points in an area."""
+ pts = self.points
+ maxlen = 0.0
+ # loop through points checking dist to next point
+ for ind, p1 in enumerate(pts[:-1]):
+ for p2 in pts[ind:]:
+ vect = p1 - p2
+ dist = vect.length()
+ if dist > maxlen:
+ maxlen = dist
+ return maxlen
+
+[docs] def rem_sline(self, sline):
+ """Removes a sline from this area."""
+ loops = [self.exlines] + self.holes
+ for loop in loops:
+ if sline in loop:
+ print('Line %s removed from %s'
+ % (sline.get_name(), self.get_name()))
+ loop.remove(sline)
+ if sline.id != -1:
+ line = sline.line
+ self.part.fea.signlines.remove(sline)
+ if len(line.signlines) == 0:
+ self.part.fea.lines.remove(line)
+
+[docs] def add_sline(self, signline):
+ """Adds signline to the area definition.
+
+ SignLine is added to the end of the lines list.
+ Area is closed if needed.
+ """
+ # this adds a line to the area
+ self.exlines.append(signline)
+ if self.exlines.closed:
+ self.close()
+
+[docs] def add_hole_sline(self, signline):
+ """Adds signline to the area hole definition.
+
+ Line is added to the end of the lines list.
+ Area is closed if needed.
+
+ Returns:
+ closed (bool): boolean telling if the hole was closed
+ """
+ make_new = (len(self.holes) == 0)
+ if make_new == False:
+ if self.holes[-1].closed == True:
+ make_new = True
+ if make_new:
+ hole = LineLoop([], hole_bool=True, parent=self)
+ hole = self.part.fea.lineloops.append(hole)
+ self.holes.append(hole)
+ self.holes[-1].append(signline)
+ # reverse the hole if it's closed in the cw direction
+ if self.holes[-1].closed == True:
+ if self.holes[-1].ccw == False:
+ self.holes[-1].reverse()
+ self.area, self.center = self.calc_area_center()
+ return self.holes[-1].closed
+
+[docs] def get_patch(self):
+ """Returns a patch of the area."""
+ if self.closed:
+ loops = [self.exlines] + self.holes
+ codes, verts = [], []
+ for loop in loops:
+ # start poly
+ codes += [Path.MOVETO]
+ verts += [loop[0].pt(0).axrad()]
+ for sline in loop:
+ if isinstance(sline, SignLine):
+ verts.append(sline.pt(1).axrad())
+ codes.append(Path.LINETO)
+ elif isinstance(sline, SignArc):
+ lverts, lcodes = sline.get_verts_codes()
+ verts += lverts
+ codes += lcodes
+ if sline == loop[-1]:
+ # close poly if at last line
+ verts.append((0, 0))
+ codes.append(Path.CLOSEPOLY)
+ path = Path(verts, codes, _interpolation_steps=2)
+ patch = PathPatch(path, linewidth=0.0)
+ return patch
+
+[docs] def label(self, axis):
+ """Labels the area on a Matplotlib axis
+
+ Args:
+ axis (Matplotlib Axis): Matplotlib Axis
+ """
+ axis.text(self.center.y, self.center.x, self.get_name(),
+ ha='center', va='center')
+
+[docs] def plot(self, axis, label=True, color='yellow'):
+ """Plots the area on a Matplotlib axis.
+
+ Args:
+ axis (Matplotlib axis): Matplotlib axis
+ label (bool): if True, label area
+ """
+ # http://stackoverflow.com/questions/8919719/how-to-plot-a-complex-polygon
+ if self.closed:
+ # plot area patch
+ patch = self.get_patch()
+ patch.set_color(color)
+ axis.add_patch(patch)
+ # label ares
+ if label:
+ self.label(axis)
+
+[docs] def contains_point(self, point):
+ """Returns bool telling if pt is inside the area.
+
+ Args:
+ pt (Point): point we are checking
+
+ Returns:
+ contains (bool): True if in this area False if not
+ """
+ patch = self.get_patch()
+ contains = patch.contains_point(point.axrad())
+ return contains
+
+[docs] def line_insert(self, lgiven, lnew, after=True):
+ """Inserts line lnew before or after the given line.
+
+ Args:
+ lgiven (SignLine or SignArc): the line we are inserting relative to
+ lnew (SignLine or SignArc): the new line we will insert
+ after (bool): if True insert lnew after lgiven, false insert before
+
+ Returns:
+ bool: True if succeeded, False if failed
+ """
+ if lgiven in self.signlines:
+ print('Inserting line %s into area %s next to %s' %
+ (lnew.get_name(), self.get_name(), lgiven.get_name()))
+ lineloop = lgiven.lineloop
+ type = 'exterior lines'
+ if lineloop.hole:
+ type = 'hole lines'
+ ind = lineloop.index(lgiven)
+ if after:
+ ind += 1
+ print(' Inserting line into %s outer lines' % self.get_name())
+ lineloop.insert(ind, lnew)
+ return True
+ else:
+ print('ERROR: Passed line was not in %s!' % self.get_name())
+ print('You must pass a line in this area!')
+ return False
+
+[docs] def set_child_ccxtypes(self):
+ """Assigns the child element ccx types."""
+ if self.etype != None and self.elements != []:
+ # if etype is set and elements have been assigned
+ for element in self.elements:
+ element.set_ccxtype(self.etype)
+
+[docs] def set_etype(self, etype):
+ """Sets the area etype.
+
+ Sets the ccxtype on elements if they have been assigned to this area.
+
+ Args:
+ etype (str): element type
+ 'plstress' plane stress
+ 'plstrain': plane strain
+ 'axisym': axisymmetric
+ """
+ self.etype = etype
+ self.set_child_ccxtypes()
+
+[docs] def set_esize(self, esize):
+ """Sets the size of mesh elements on the area when meshing.
+
+ Args:
+ esize (float): size of mesh elements on this area
+ """
+ #size if set if no mesh size on points is peresent
+ for p in self.points:
+ p.set_esize(esize)
+
+[docs] def update(self, line_list):
+ """Updates the area's external lines to the passed list.
+
+ If the area is cosed, self.close() will be called.
+
+ Args:
+ line_list (list): list of SignLine and SignArc items
+ """
+ to_rem = list(self.exlines)
+ for sline in to_rem:
+ self.exlines.remove(sline)
+ for sline in line_list:
+ self.exlines.append(sline)
+ if self.exlines.closed == True:
+ self.close()
+
+
+import glob
+import os
+import re
+import shutil
+import shlex
+import subprocess
+import sys
+from urllib.parse import urlparse
+from urllib.parse import ParseResult
+import zipfile
+
+from sys import platform as platform
+
+import requests
+
+WINDOWS_PLATFORMS = ['win32', 'win64']
+LINUX_PLATFORMS = ['linux', 'linux2']
+
+[docs]def printos(os, bits):
+ """Prints the operating system and bit size"""
+ print('Detected {} {} bit'.format(os, bits))
+
+[docs]def add():
+ """Installs the fea tools on windows/mac/linux"""
+ osname = None
+ is_64bit = sys.maxsize > 2**32
+ bitsize_dict = {True: 64, False: 32}
+ bitsize = bitsize_dict[is_64bit]
+ if platform in LINUX_PLATFORMS:
+ printos('Linux', bitsize)
+ ubuntu_add()
+ elif platform == "darwin":
+ printos('Mac OS X', bitsize)
+ mac_add()
+ elif platform in WINDOWS_PLATFORMS:
+ printos('Windows', bitsize)
+ windows_add(bitsize)
+ print('Done!')
+
+[docs]def remove():
+ """Removes the fea tools on windows/mac/linux"""
+ osname = None
+ is_64bit = sys.maxsize > 2**32
+ bitsize_dict = {True: 64, False: 32}
+ bitsize = bitsize_dict[is_64bit]
+ if platform in LINUX_PLATFORMS:
+ printos('Linux', bitsize)
+ ubuntu_remove()
+ elif platform == "darwin":
+ printos('Mac OS X', bitsize)
+ mac_remove()
+ elif platform in WINDOWS_PLATFORMS:
+ printos('Windows', bitsize)
+ windows_remove(bitsize)
+ print('Done!')
+
+[docs]def ubuntu_add():
+ """Adds ccx and gmsh programs on ubuntu, uses apt"""
+ gmsh_installed = shutil.which('gmsh')
+ if not gmsh_installed:
+ print('Installing gmsh')
+ command_line = "sudo apt-get install gmsh"
+ subprocess.check_call(command_line, shell=True)
+ else:
+ print('gmsh present')
+ ccx_installed = shutil.which('ccx')
+ if not ccx_installed:
+ print('Installing calculix (ccx)')
+ command_line = "sudo apt-get install calculix-ccx"
+ subprocess.check_call(command_line, shell=True)
+ else:
+ print('calculix (ccx) present')
+
+[docs]def ubuntu_remove():
+ """Removes ccx and gmsh programs on ubuntu, uses apt"""
+ ccx_installed = shutil.which('ccx')
+ if not ccx_installed:
+ print('calculix (ccx) is not on your system')
+ else:
+ print('Removing calculix (ccx)')
+ command_line = "sudo apt-get remove calculix-ccx"
+ subprocess.check_call(command_line, shell=True)
+ gmsh_installed = shutil.which('gmsh')
+ if not gmsh_installed:
+ print('gmsh is not on your system')
+ else:
+ print('Removing gmsh')
+ command_line = "sudo apt-get remove gmsh"
+ subprocess.check_call(command_line, shell=True)
+
+[docs]def mac_add():
+ """Adds ccx and gmsh programs on mac, uses brew"""
+ brew_installed = shutil.which('brew')
+ if not brew_installed:
+ print('Installing brew')
+ url = 'https://raw.githubusercontent.com/Homebrew/install/master/install'
+ command_line = "/usr/bin/ruby -e \"$(curl -fsSL %s)\"" % url
+ subprocess.check_call(command_line, shell=True)
+ else:
+ print('brew present')
+ gmsh_installed = shutil.which('gmsh')
+ if not gmsh_installed:
+ print('Installing gmsh')
+ folder_path = os.path.dirname(os.path.abspath(__file__))
+ dmginstall_path = os.path.join(folder_path, 'dmginstall.sh')
+ url = 'http://gmsh.info/bin/MacOSX/gmsh-3.0.5-MacOSX.dmg'
+ command_line = '%s %s' % (dmginstall_path, url)
+ print('command_line=%s' % command_line)
+ subprocess.check_call(command_line, shell=True)
+ gmsh_path = '/Applications/Gmsh.app/Contents/MacOS/gmsh'
+ command_line = "ln -s %s /usr/local/bin/gmsh" % gmsh_path
+ subprocess.check_call(command_line, shell=True)
+ else:
+ print('gmsh present')
+ ccx_installed = shutil.which('ccx')
+ if not ccx_installed:
+ mac_add_ccx()
+ else:
+ print('calculix (ccx) present')
+
+[docs]def mac_add_ccx():
+ print('Installing calculix (ccx)')
+ command_line = "brew install brewsci/science/calculix-ccx"
+ subprocess.check_call(command_line, shell=True)
+ ccx_path = find_brew_binary_location('calculix-ccx', 'ccx')
+ if not ccx_path:
+ raise Exception('Failed to find ccx binary')
+ # link the ccx command to the ccx binary
+ command_line = "ln -s %s /usr/local/bin/ccx" % ccx_path
+ subprocess.check_call(command_line, shell=True)
+ # check to see if we have the fortran that ccx needs and return if we do
+ gcc7_to_path = "/usr/local/opt/gcc/lib/gcc/7"
+
+ needed_fortran_path = "%s/libgfortran.4.dylib" % gcc7_to_path
+ if os.path.isfile(needed_fortran_path):
+ if os.path.islink(link_to_path):
+ print('Linked gcc7 fortran found and used')
+ else:
+ print('System gcc7 fortran found and used')
+ print('Finished installing calculix (ccx)')
+ return
+ # install gcc@7 with brew if we don't have it
+ paths = glob.glob("/usr/local/Cellar/gcc@7/*")
+ if not paths:
+ command_line = "brew install gcc@7"
+ subprocess.check_call(command_line, shell=True)
+ print('Installed gcc@7 (needed by calculix)')
+ paths = glob.glob("/usr/local/Cellar/gcc@7/*")
+ # link gcc@7 from_path to to_path
+ gcc7_from_path = "%s/lib/gcc/7" % paths[-1]
+ command_line = "ln -s %s %s" % (gcc7_from_path, gcc7_to_path)
+ subprocess.check_call(command_line, shell=True)
+ if not os.path.isfile(needed_fortran_path):
+ raise Exception("Install failed, libgfortran.4.dylib (a library that "
+ "ccx needs) was NOT FOUND")
+ print('Finished installing calculix (ccx)')
+
+
+[docs]def find_brew_binary_location(package_folder, search_string):
+ """Finds the location of a binary installed by homebrew"""
+ match_str = '/usr/local/Cellar/%s/**/*%s*' % (package_folder,
+ search_string)
+ paths = glob.glob(match_str, recursive=True)
+ for path in paths:
+ if os.access(path, os.X_OK):
+ return path
+ return None
+
+[docs]def mac_remove():
+ """Removes ccx and gmsh programs on mac, uses brew"""
+ brew_installed = shutil.which('brew')
+ if not brew_installed:
+ print('Installing brew')
+ url = 'https://raw.githubusercontent.com/Homebrew/install/master/install'
+ command_line = "/usr/bin/ruby -e \"$(curl -fsSL %s)\"" % url
+ subprocess.check_call(command_line, shell=True)
+ else:
+ print('brew present')
+ ccx_installed = shutil.which('ccx')
+ if not ccx_installed:
+ print('calculix (ccx) is not on your system')
+ else:
+ print('Removing calculix (ccx)')
+ # remove link to ccx
+ command_line = "rm /usr/local/bin/ccx"
+ subprocess.check_call(command_line, shell=True)
+ print('Calculix (ccx) symlink was removed')
+ # remove binary
+ command_line = "brew uninstall calculix-ccx"
+ subprocess.check_call(command_line, shell=True)
+ print('Calculix (ccx) binary was removed')
+ # remove gcc7 link if it exists
+ gcc7_sys_path = "/usr/local/opt/gcc/lib/gcc/7"
+ if os.path.islink(gcc7_sys_path):
+ command_line = "rm %s" % gcc7_sys_path
+ subprocess.check_call(command_line, shell=True)
+ print('gcc@7 symlink was removed (it was needed by calculix)')
+ # remove gcc@7 if it exists
+ gcc7_brew_path = "/usr/local/Cellar/gcc@7"
+ if os.path.isdir(gcc7_brew_path):
+ command_line = "brew uninstall gcc@7"
+ subprocess.check_call(command_line, shell=True)
+ print('gcc@7 was removed (it was needed by calculix)')
+
+ gmsh_installed = shutil.which('gmsh')
+ if not gmsh_installed:
+ print('gmsh is not on your system')
+ else:
+ print('Removing gmsh')
+ command_line = "rm /usr/local/bin/gmsh"
+ subprocess.check_call(command_line, shell=True)
+ command_line = "rm -rf /Applications/Gmsh.app"
+ subprocess.check_call(command_line, shell=True)
+
+
+[docs]def windows_add(bitsize):
+ """Adds ccx and gmsh programs on windows"""
+ gmsh_installed = shutil.which('gmsh')
+ if not gmsh_installed:
+ print('Installing gmsh')
+ url = 'http://gmsh.info/bin/Windows/'
+ win_add_gmsh(bitsize, url, 'gmsh', '3.0.5')
+ else:
+ print('gmsh present')
+
+ ccx_installed = shutil.which('ccx')
+ if not ccx_installed:
+ print('Installing calculix (ccx)')
+ url = 'https://sourceforge.net/projects/calculixforwin/files/03.2/'
+ win_add_ccx(bitsize, url, 'ccx')
+ else:
+ print('calculix (ccx) present')
+
+[docs]def windows_remove(bitsize):
+ """Removes ccx and gmh programs on windows"""
+ gmsh_installed = shutil.which('gmsh')
+ if not gmsh_installed:
+ print('gmsh is not on your system')
+ else:
+ print('Removing gmsh')
+ env_path = os.getenv('VIRTUAL_ENV', sys.exec_prefix)
+ scripts_folder = '%s\Scripts\\' % env_path
+ remove_like(scripts_folder, 'gmsh')
+ ccx_installed = shutil.which('ccx')
+ if not ccx_installed:
+ print('ccx is not on your system')
+ else:
+ print('Removing ccx')
+ env_path = os.getenv('VIRTUAL_ENV', sys.exec_prefix)
+ scripts_folder = '%s\Scripts\\' % env_path
+ ccx_folder = '%sccx' % scripts_folder
+ add_remove_dll_links(ccx_folder, scripts_folder, add=False)
+ remove_like(scripts_folder, 'ccx')
+
+[docs]def remove_like(search_path, name):
+ """
+ Searches for files or folders matching name in search_path and deletes them
+ """
+ match_str = '%s*%s*' % (search_path, name)
+ paths = glob.glob(match_str)
+ for path in paths:
+ print('Removing %s' % path)
+ if os.path.isfile(path):
+ os.unlink(path)
+ elif os.path.isdir(path):
+ shutil.rmtree(path)
+
+[docs]def win_add_gmsh(bitsize, binaries_url, program_name, version_str):
+ """Installs a program from an apache server file listing page"""
+ # Needs the user agent header to exist for it to send back the content
+ user_agent = ('Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:56.0) '
+ 'Gecko/20100101 Firefox/56.0')
+ headers = {'User-Agent': user_agent}
+
+ zipfile_regex = '.*%s-%s.*' % (program_name, version_str)
+ # Windows 10 64 bit, gmsh version: Gmsh 3.0.7-git-35176e26 hangs
+ zipfile_name = zipfile_by_bitsize(binaries_url, headers, zipfile_regex,
+ bitsize)
+
+ zipfile_url = binaries_url + zipfile_name
+ print('Downloading %s from %s' % (program_name, zipfile_url))
+ response = requests.get(zipfile_url, stream=True, headers=headers)
+ with open(zipfile_name, 'wb') as out_file:
+ for chunk in response.iter_content(chunk_size=1024):
+ if chunk: # filter out keep-alive new chunks
+ out_file.write(chunk)
+ print('Unzipping %s' % program_name)
+ zip_ref = zipfile.ZipFile(zipfile_name, 'r')
+ zipfile_folder_name = os.path.dirname(zip_ref.namelist()[0])
+ zip_ref.extractall(None)
+ zip_ref.close()
+ print('Removing %s zipfile' % program_name)
+ os.remove(zipfile_name)
+
+ env_path = os.getenv('VIRTUAL_ENV', sys.exec_prefix)
+ folder_to = '%s\Scripts\%s' % (env_path, zipfile_folder_name)
+ exe_loc = '%s\%s.exe' % (folder_to, program_name)
+ exe_link = '%s\Scripts\%s.exe' % (env_path, program_name)
+ paths = [exe_link, folder_to]
+ for path in paths:
+ if os.path.isfile(path):
+ os.unlink(path)
+ elif os.path.isdir(path):
+ shutil.rmtree(path)
+ print('Installing %s to %s' % (program_name, folder_to))
+ shutil.move(zipfile_folder_name, folder_to)
+ os.link(exe_loc, exe_link)
+
+[docs]def zipfile_by_bitsize(binaries_url, headers, zipfile_regex, bitsize):
+ """Returns the url linking to the correct zipfile"""
+ # this is used by ccx and gmsh
+ res = requests.get(binaries_url, headers=headers)
+ html = res.text
+ urls = re.findall(r'href=[\'"]?([^\'" >]+)', html)
+ pattern = re.compile(zipfile_regex)
+ urls = [url for url in urls if pattern.match(url)]
+ urls = urls[-2:]
+ url_choices = {32: urls[0], 64: urls[1]}
+ if 'win32' in urls[1] or 'Windows32' in urls[1]:
+ url_choices = {32: urls[1], 64: urls[0]}
+ return url_choices[bitsize]
+
+[docs]def href_from_link_text(url, headers, link_text):
+ """
+ Returns the url for a link with link_text description, if this
+ function fails, it raises an error with tells the users a link
+ which creates an issue on the pycalculix repo
+ """
+ response = requests.get(url, headers=headers)
+ html = response.text.replace("\"", "'")
+ link_text_pos = html.find(link_text)
+ if link_text_pos == -1:
+ issue_url = ('https://github.com/spacether/pycalculix/issues/new?title='
+ 'CCX%20zip%20download%20fails%20on%20windows&body=Please%2'
+ '0update%20the%20installer.py%20get_direct_url%20function.'
+ '%20It%20no%20longer%20works')
+ raise ValueError('Unable to download file because there was not a link '
+ 'with text=\'%s\' on the page at url=%s\nTo fix this, '
+ 'please click the \'Submit new issue\' button '
+ 'here:\n%s' % (link_text, url, issue_url))
+ href_pos = html[:link_text_pos].rfind('href')
+ first_char = html.find("'", href_pos)+1
+ last_char = html.find("'", first_char)
+ return html[first_char:last_char]
+
+[docs]def get_direct_url(url, headers):
+ """Gets the zip direct download link from the project download page"""
+ direct_download_url = href_from_link_text(url,
+ headers,
+ 'Problems Downloading')
+ parsed_download_url = urlparse(direct_download_url)
+ if parsed_download_url.scheme not in ['http', 'https']:
+ # url is relative, and is missing the scheme and netloc
+ parsed_parent_url = urlparse(url)
+ parsed_download_url = ParseResult(parsed_parent_url.scheme,
+ parsed_parent_url.netloc,
+ parsed_download_url.path,
+ parsed_download_url.params,
+ parsed_download_url.query,
+ parsed_download_url.fragment)
+ direct_download_url = parsed_download_url.geturl()
+ direct_download_url = href_from_link_text(direct_download_url,
+ headers,
+ 'direct link')
+ return direct_download_url
+
+[docs]def win_add_ccx(bitsize, binaries_url, program_name):
+ """Installs ccx on a windows computer"""
+ # Needs the user agent header to exist for it to send back the content
+ user_agent = ('Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:56.0) '
+ 'Gecko/20100101 Firefox/56.0')
+ headers = {'User-Agent': user_agent}
+
+ zipfile_regex = '.+CL.+win.+zip/download$'
+ zipfile_webpage_url = zipfile_by_bitsize(binaries_url, headers,
+ zipfile_regex, bitsize)
+
+ zipfile_name = zipfile_webpage_url.split('/')[-2]
+ zipfile_url = get_direct_url(zipfile_webpage_url, headers)
+
+ print('Downloading %s from %s' % (zipfile_name, zipfile_url))
+ response = requests.get(zipfile_url, stream=True, headers=headers)
+ with open(zipfile_name, 'wb') as out_file:
+ for chunk in response.iter_content(chunk_size=1024):
+ if chunk: # filter out keep-alive new chunks
+ out_file.write(chunk)
+ print('Unzipping %s' % program_name)
+ zip_ref = zipfile.ZipFile(zipfile_name, 'r')
+ zipfile_folder_name = zip_ref.namelist()[0].split('/')[0]
+ zipfile_ccx_path = zipfile_folder_name + '/bin/ccx'
+ for member in zip_ref.namelist():
+ if zipfile_ccx_path in member:
+ zip_ref.extract(member)
+ zip_ref.close()
+ print('Removing %s zipfile' % program_name)
+ os.remove(zipfile_name)
+
+ folder_from = '%s\\bin\\ccx' % zipfile_folder_name
+ env_path = os.getenv('VIRTUAL_ENV', sys.exec_prefix)
+ scripts_path = '%s\Scripts' % env_path
+ folder_to = '%s\%s' % (scripts_path, program_name)
+
+ binary_name = 'ccx212.exe'
+ exe_loc = '%s\%s' % (folder_to, binary_name)
+ exe_link = '%s\%s.exe' % (scripts_path, program_name)
+ paths = [exe_link, folder_to]
+ for path in paths:
+ if os.path.isfile(path):
+ os.unlink(path)
+ elif os.path.isdir(path):
+ shutil.rmtree(path)
+
+ print('Installing %s to %s' % (program_name, folder_to))
+ shutil.move(folder_from, folder_to)
+ shutil.rmtree(zipfile_folder_name)
+ add_remove_dll_links(folder_to, scripts_path, add=True)
+ os.link(exe_loc, exe_link)
+
+[docs]def add_remove_dll_links(folder_dlls, folder_dll_links, add=True):
+ """
+ Finds all dlls in folder_from and adds hard links to them in folder_to
+ Or removes the hard links if add=False
+ """
+ match_str = '%s\*.dll' % folder_dlls
+ dll_paths = glob.glob(match_str)
+ for dll_path in dll_paths:
+ dll_name = dll_path.split('\\')[-1]
+ dll_link_path = '%s\%s' % (folder_dll_links, dll_name)
+ if add:
+ os.link(dll_path, dll_link_path)
+ elif os.path.isfile(dll_link_path):
+ os.unlink(dll_link_path)
+
+"""This module stores load classes.
+"""
+from math import pi
+
+[docs]class ConstLoad(object):
+ """Makes a load. Many types are supported.
+
+ Possible load types:
+ - forces
+ - displacements
+ - thickness (on nodes)
+ - matl (on elements)
+ - pressure (or stress on faces)
+ - gravity (on elements)
+ - rotation about an axis (on elements)
+
+ Args:
+ ltype (string): load type:
+
+ - 'fx','fy','fz': force on each axis
+ - 'ux','uy','uz': displacement on each axis
+ - 'nodal_thickness': thickness on element nodes
+ - 'matl': matl on elements
+ - 'press': pressure, + is tension, - is compresion
+ - 'gravity': gravity in x axis direction, - goes towards y axis
+ - 'radps': radians per second rotation
+ - 'rpm': rotations per minute rotation
+
+ comp (Component): component that the load is applied to
+ val (double or Matl): value of the load
+ Matl is needed when setting material loads
+ otherwise a double value is used to describe the load
+ """
+
+ def __init__(self, ltype, comp, val):
+ self.ltype = ltype
+ self.comp = comp
+ self.val = val
+
+[docs] def get_list(self):
+ """Returns a list of [item, val] for plotting."""
+ res = []
+ if self.ltype == 'press':
+ # pressure
+ faces = self.comp.get_children()
+ for face in faces:
+ res.append([face, self.val])
+ elif self.ltype in ['ux', 'uy', 'uz']:
+ nodes = self.comp.get_children()
+ load_dict = {self.ltype:self.val}
+ for node in nodes:
+ res.append([node, load_dict])
+ return res
+
+[docs] def ccx(self):
+ """Writes a load for Calculix ccx solver.
+
+ Returns a list strings, where each item is a line to write
+ to a Calculix .inp text file.
+ """
+ res = []
+ cname = self.comp.name
+ if self.ltype == 'press':
+ # pressure
+ res.append('** Pressure on component: '+cname)
+ res.append('*DLOAD')
+ faces = self.comp.get_children()
+ for face in faces:
+ myline = '%i,P%i,%f' % (face.element.id, face.id, self.val)
+ res.append(myline)
+ elif self.ltype == 'gravity':
+ # gravity
+ res.append('** Gravity on component: '+cname)
+ res.append('*DLOAD')
+ res.append('%s,GRAV,%f,-1.,0.,0.' % (cname, self.val))
+ elif self.ltype[0] == 'u':
+ # displacement
+ axis = self.ltype[1]
+ anum = {'x':1, 'y':2, 'z':3}
+ axis = anum[axis]
+ res.append('*BOUNDARY')
+ res.append('%s,%i,%i,%f' % (cname, axis, axis, self.val))
+ elif self.ltype[0] == 'f':
+ # force
+ res.append('** Force on component: '+cname)
+ nodes = len(self.comp.items)
+ fval = self.val/nodes
+ axis = self.ltype[1]
+ anum = {'x':1, 'y':2, 'z':3}
+ axis = anum[axis]
+ res.append('*CLOAD')
+ res.append('%s,%i,%f' % (cname, axis, fval))
+ elif self.ltype == 'nodal_thickness':
+ res.append('*NODAL THICKNESS')
+ res.append('%s, %f' % (cname, self.val))
+ elif self.ltype == 'radps':
+ radsq = self.val**2 # square the speed in radians per sec
+ res.append('** Rotation force on component: '+cname)
+ res.append('** Speed = '+str(self.val)+' radians/sec')
+ res.append('*DLOAD')
+ res.append('%s,CENTRIF,%f,0.,0.,0.,0.,1.,0.' % (cname, radsq))
+ elif self.ltype == 'rpm':
+ rad = 2.0*pi*(self.val/60.0)
+ radsq = rad**2 # square the speed in radians per sec
+ res.append('** Rotation force on component: '+cname)
+ res.append('** Speed = '+str(self.val)+' rotations/minute')
+ res.append('*DLOAD')
+ res.append('%s,CENTRIF,%f,0.,0.,0.,0.,1.,0.' % (cname, radsq))
+ elif self.ltype == 'matl':
+ mat = self.val.name
+ res.append('*SOLID SECTION,ELSET=%s,MATERIAL=%s' %
+ (cname, mat))
+
+ return res
+
+
+[docs]class LinearLoad(object):
+ """Makes a load which varies depending on x, y, or z location.
+
+ Load Equation:
+ Load(x) = const + mult*(xo - x)
+
+ This load is used to set water pressure, where p = po + rho*g*(xo - x)
+ For water pressure this becomes:
+
+ - P(x) = po + (rho*g)*(xo-x)
+ - P(x) = po + mult*(xo-x)
+ - g is positive, xo > x
+
+ Args:
+ ltype (str): load type 'press_fluid' is the only valid option
+ comp (Component): component that the load acts on
+ const (float): constant term in load equation
+ mult (float): mult term in the load equation
+ xo (float): xo term in the load equation
+ axis (str): axis that the load depends on
+ 'x' or 'y' are valid options
+ """
+
+ def __init__(self, ltype, comp, const, mult, xo, axis='x'):
+ self.ltype = ltype
+ self.comp = comp
+ self.const = const
+ self.mult = mult
+ self.xo = xo
+ self.axis = axis
+
+[docs] def get_val(self, xval):
+ """Returns the load function value at the given location.
+
+ Args:
+ xval (float): the location where we want to find the load
+ """
+ res = self.const + (self.xo - xval)*self.mult
+ return res
+
+[docs] def get_list(self):
+ """Returns a list of [item, val] for plotting."""
+ res = []
+ faces = self.comp.get_children()
+ for face in faces:
+ xvals = [getattr(n, self.axis) for n in face.nodes]
+ xavg = sum(xvals)/len(xvals)
+ pval = self.get_val(xavg)
+ res.append([face, pval])
+ return res
+
+[docs] def ccx(self):
+ """Writes a load for Calculix ccx solver.
+
+ Returns a list strings, where each item is a line to write
+ to a Calculix .inp text file.
+ """
+ res = []
+ res.append('** Fluid pressure on component: '+self.comp.name)
+ faces = self.comp.get_children()
+ for face in faces:
+ res.append('*DLOAD')
+ xvals = [getattr(n, self.axis) for n in face.nodes]
+ xavg = sum(xvals)/len(xvals)
+ pval = self.get_val(xavg)
+ myline = '%i,P%i,%f' % (face.element.id, face.id, pval)
+ res.append(myline)
+ return res
+
+"""This module stores material classes.
+"""
+
+from . import base_classes
+
+[docs]class Material(base_classes.Idobj):
+ """Makes a linear elastic meterial.
+
+ Args:
+ name (str): a unique matl name
+
+ Attributes:
+ - id (int): material id number
+ - mechtype (str): defines mechanical material type: linear or non-linear material
+ - density (float): density in mass/volume units
+ - pratio (float): poisson's ratio (unitless)
+ - youngs (float): young's modulus in Force/area = stress units
+ - mechtpye (str): mechanical material type
+ options: 'linear' or 'nonlinear'
+ - exponent (float): exponent of Ramberg-Osgood stress-strain equation
+ - yield_stress (float): yield stress of material
+ - yield_offset (float): yield offset of Ramberg-Osgood stress-strain equation
+ - conductivity (float): thermal conductivity, Power / (distance-temp)
+ - spec_heat (float): specific heat (energy/(mass-temp)
+ - thermal_exp (dict): a dict storing temperature dependent thermal
+ -- expansion properties
+
+ -- Thermal expansion is in strain per temperature
+
+ - dict['data'] = zip of (temp, thermal_expansion)
+ - dict['tzero'] = the temperature zero point
+ """
+
+ def __init__(self, name):
+ self.name = name
+ base_classes.Idobj.__init__(self)
+ # mechanical
+ self.mechtype = None
+ self.density = None
+ self.youngs = None
+ self.pratio = None
+ self.exponent = None
+ self.yield_stress = None
+ self.yield_offset = None
+ # thermal
+ self.conductivity = None
+ self.spec_heat = None
+ # thermal growth
+ self.thermal_exp = {}
+
+[docs] def set_mech_props(self, density, youngs, pratio, mechtype='linear', exponent=0., yield_stress=0., yield_offset=0.):
+ """Sets the mechanical properties: density, youngs, poisson_ratio.
+
+ Args:
+ - density (float): density in mass/volume units
+ - pratio (float): poisson's ratio (unitless)
+ - youngs (float): young's modulus in Force/area = stress units
+
+ Kargs:
+ - mechtpye (str): mechanical material type
+ options: 'linear' or 'nonlinear'
+ - exponent (float): exponent of Ramberg-Osgood stress-strain equation
+ - yield_stress (float): yield stress of material
+ - yield_offset (float): yield offset of Ramberg-Osgood stress-strain equation
+ """
+ self.density = density
+ self.youngs = youngs
+ self.pratio = pratio
+ self.mechtype = mechtype
+ self.exponent = exponent
+ self.yield_stress = yield_stress
+ self.yield_offset = yield_offset
+
+[docs] def set_therm_props(self, conductivity, spec_heat):
+ """Sets the thermal properties: conductivity, specifc_heat.
+
+ Args:
+ - conductivity (float): Power / (distance-temp)
+ - spec_heat (float): specific heat (energy/(mass-temp)
+ """
+ self.conductivity = conductivity
+ self.spec_heat = spec_heat
+
+[docs] def set_therm_expan(self, alphas, temps=None, tzero=None):
+ """Sets the thermal expansion of the material.
+
+ Args:
+ - alphas (float or list): list of thermal expansion alphas
+ length must be the same as the passed temps
+ - temps (list): list of temperatures
+ - tzero (float): temperature zero point
+ """
+ self.thermal_exp = {}
+ isnum = (isinstance(alphas, float) or isinstance(alphas, int))
+ if isinstance(alphas, list) and isinstance(temps, list):
+ self.thermal_exp['data'] = zip(temps, alphas)
+ elif temps == None and isnum:
+ self.thermal_exp['data'] = [alphas]
+ if tzero != None:
+ self.thermal_exp['tzero'] = tzero
+
+[docs] def ccx(self):
+ """Returns a list of text strings for ccx defining the material."""
+ res = []
+ res.append('*MATERIAL,NAME='+self.name)
+
+ if self.mechtype == 'linear':
+ res.append('*ELASTIC')
+ res.append(str(self.youngs)+','+str(self.pratio))
+ elif self.mechtype == 'nonlinear':
+ res.append('*DEFORMATION PLASTICITY')
+ res.append(','.join([str(self.youngs),str(self.pratio),str(self.yield_stress),str(self.exponent),str(self.yield_offset)]))
+ if self.density != None:
+ res.append('*DENSITY')
+ res.append(str(self.density))
+ if self.conductivity != None:
+ res.append('*CONDUCTIVITY')
+ res.append(str(self.conductivity))
+ if self.spec_heat != None:
+ res.append('*SPECIFIC HEAT')
+ res.append(str(self.spec_heat))
+ if self.thermal_exp != {}:
+ if 'tzero' in self.thermal_exp:
+ res.append('*EXPANSION,ZERO='+str(self.thermal_exp['tzero']))
+ else:
+ res.append('*EXPANSION')
+ for pair in self.thermal_exp['data']:
+ if len(pair) == 1:
+ res.append(str(pair[0]))
+ else:
+ # temp, val
+ res.append('%f %f' % (pair[0], pair[1]))
+ return res
+
+"""This module stores classes that make a finite element analysis mesh.
+"""
+from matplotlib.patches import Polygon # needed for plotting elements
+from . import geometry
+
+[docs]class Element(object):
+ """Makes a mesh element.
+
+ Supports 4 noded quads, 8 noded quads, 3 noded tris, 6 noded tris.
+
+ Args:
+ enum (int): element id number
+ ccxtype (str): ccx element type string, 4 characters long
+
+ Description: element_shape + element_order + element_type
+
+ - element_shape = 'quad' or 'tri'
+ - 'quad' = quadrangle
+ - 'tri' = triangle
+ - element_order = '1' or '2'
+ - '1' = corner nodes only (3 or 4 noded element)
+ - '2' = corner nodes and midside nodes (6 or 8 noded element)
+ - element_type = 'plstress' or 'plstrain' or 'axisym'
+ - 'plstress' = plane stress
+ - 'plstrain' = plane strain
+ - 'axisym' = axisymmetric
+
+ axisymmetric
+ - 'CAX6' : 'tri2axisym'
+ - 'CAX3' : 'tri1axisym'
+ - 'CAX8' : 'quad2axisym'
+ - 'CAX4' : 'quad1axisym'
+ plane stress
+ - 'CPS6' : 'tri2plstress'
+ - 'CPS3' : 'tri1plstress'
+ - 'CPS8' : 'quad2plstress'
+ - 'CPS4' : 'quad1plstress'
+ plane strain
+ - 'CPE6' : 'tri2plstrain'
+ - 'CPE3' : 'tri1plstrain'
+ - 'CPE8' : 'quad2plstrain'
+ - 'CPE4' : 'quad1plstrain'
+
+ nlist (list): list of Node, the nodes that define the element
+
+ Attributes:
+ id (int): element id number
+ ccxtype (str): ccx element type string, 4 characters long
+
+ Description: element_shape + element_order + element_type
+
+ - element_shape = 'quad' or 'tri'
+ - 'quad' = quadrangle
+ - 'tri' = triangle
+ - element_order = '1' or '2'
+ - '1' = corner nodes only (3 or 4 noded element)
+ - '2' = corner nodes and midside nodes (6 or 8 noded element)
+ - element_type = 'plstress' or 'plstrain' or 'axisym'
+ - 'plstress' = plane stress
+ - 'plstrain' = plane strain
+ - 'axisym' = axisymmetric
+
+ axisym
+ - 'CAX6' : 'tri2axisym'
+ - 'CAX3' : 'tri1axisym'
+ - 'CAX8' : 'quad2axisym'
+ - 'CAX4' : 'quad1axisym'
+ plane stress
+ - 'CPS6' : 'tri2plstress'
+ - 'CPS3' : 'tri1plstress'
+ - 'CPS8' : 'quad2plstress'
+ - 'CPS4' : 'quad1plstress'
+ plane strain
+ - 'CPE6' : 'tri2plstrain'
+ - 'CPE3' : 'tri1plstrain'
+ - 'CPE8' : 'quad2plstrain'
+ - 'CPE4' : 'quad1plstrain'
+
+ node (dict): dictionary of nodes, keys are int >= 1
+ face (dict): dictionary of faces, keys are int >= 1
+ nodes (list): list of element nodes
+ faces (list): list of element faces
+ center (Point): the center point of the element
+ """
+
+ def __init__(self, enum, ccxtype, nlist):
+ self.id = enum
+ self.ccxtype = ccxtype
+ self.node = {}
+ self.face = {}
+
+ # calculate the number of faces
+ fnum = len(nlist)
+ if fnum > 4:
+ fnum = int(fnum/2)
+
+ # store nodes
+ for (ind, node) in enumerate(nlist):
+ if ind >= fnum:
+ node.set_order(2)
+ node.add_element(self)
+ self.node[ind+1] = node
+
+ # store faces
+ nodes = nlist[0:fnum]+[nlist[0]]
+ for ind in range(fnum):
+ node1 = nodes[ind]
+ node2 = nodes[ind+1]
+ face = Face(ind+1, node1, node2, self)
+ if len(nlist) > fnum:
+ # we have second order nodes
+ face.set_nmid(self.node[ind+1+fnum])
+ self.face[ind+1] = face
+
+ # set the element center
+ self.center = self.calc_center()
+
+ def __hash__(self):
+ """Returns the item's id as its hash."""
+ return self.id
+
+ @property
+ def faces(self):
+ """Returns list of element faces."""
+ return self.face.values()
+
+ @property
+ def nodes(self):
+ """Returns list of element nodes."""
+ return self.node.values()
+
+[docs] def get_poly(self):
+ """Returns a polygon of the element."""
+ nodes = [node for node in self.nodes if node.order == 1]
+ # make a list of [[ax, rad],[... ] for making the polygon
+ xycoords = [[node.y, node.x] for node in nodes]
+ polygon = Polygon(xycoords, closed=True)
+ return polygon
+
+[docs] def label(self, ax):
+ """Labels the element on the passed matplotlib axis."""
+ ax.text(self.center.y, self.center.x, self.get_name(),
+ ha='center', va='center')
+
+[docs] def calc_center(self):
+ """Returns the element center as a Point."""
+ pts = [node for node in self.nodes if node.order == 1]
+ axials = [point.y for point in pts]
+ radials = [point.x for point in pts]
+ axial = sum(axials)/len(axials)
+ radial = sum(radials)/len(radials)
+ return geometry.Point(radial, axial)
+
+[docs] def get_tris(self):
+ """Returns a list of triangles for plotting. Triangles are node ids.
+
+ CCX elements are closed in a CCW direction relative to xy normal
+ CCX elements are closed in a CW direction relative to yx mine
+ ---> I have to draw in a clockwise direction, otherwise get error
+ Triangles are closed in a counter clockwise (CCW) direction to yx mine
+ http://matplotlib.org/api/tri_api.html#matplotlib.tri.Triangulation
+ """
+ res = []
+ numnodes = len(self.nodes)
+ if numnodes == 3:
+ # triangle, first order = 1 triangle
+ res.append([self.node[1].id, self.node[3].id, self.node[2].id])
+ elif numnodes == 4:
+ # quad, first order = 2 triangles
+ res.append([self.node[1].id, self.node[3].id, self.node[2].id])
+ res.append([self.node[1].id, self.node[4].id, self.node[3].id])
+ elif numnodes == 6:
+ # tri, second order = 4 triangles
+ res.append([self.node[1].id, self.node[6].id, self.node[4].id])
+ res.append([self.node[4].id, self.node[5].id, self.node[2].id])
+ res.append([self.node[6].id, self.node[3].id, self.node[5].id])
+ res.append([self.node[6].id, self.node[5].id, self.node[4].id])
+ elif numnodes == 8:
+ # quad, second order = 6 triangles
+ res.append([self.node[1].id, self.node[8].id, self.node[5].id])
+ res.append([self.node[5].id, self.node[6].id, self.node[2].id])
+ res.append([self.node[6].id, self.node[7].id, self.node[3].id])
+ res.append([self.node[8].id, self.node[4].id, self.node[7].id])
+ res.append([self.node[8].id, self.node[7].id, self.node[6].id])
+ res.append([self.node[8].id, self.node[6].id, self.node[5].id])
+ return res
+
+[docs] def get_area(self):
+ """Returns the element area."""
+ asum = 0.0
+ for face in self.faces:
+ node1 = face.nodes[0]
+ node2 = face.nodes[1]
+ asum += node2.y*node1.x - node2.x*node1.y
+ asum = asum*0.5
+ return asum
+
+[docs] def set_ccxtype(self, etype):
+ """Sets the ccx element type given an etype string.
+
+ Args:
+ etype (string): element type string
+
+ - 'plstress' = plane stress
+ - 'plstrain' = plane strain
+ - 'axisym' = axisymmetric
+ """
+ fnum = len(self.faces)
+ shape = 'tri'
+ if fnum == 4:
+ shape = 'quad'
+ order = '2'
+ if len(self.nodes) == fnum:
+ order = '1'
+ estr = shape+order+etype
+
+ # Calculix CCX elements
+ # axisymmetric
+ _ccx_elements = {}
+ _ccx_elements['tri2axisym'] = 'CAX6'
+ _ccx_elements['tri1axisym'] = 'CAX3'
+ _ccx_elements['quad2axisym'] = 'CAX8'
+ _ccx_elements['quad1axisym'] = 'CAX4'
+ # plane stress
+ _ccx_elements['tri2plstress'] = 'CPS6'
+ _ccx_elements['tri1plstress'] = 'CPS3'
+ _ccx_elements['quad2plstress'] = 'CPS8'
+ _ccx_elements['quad1plstress'] = 'CPS4'
+ # plane strain
+ _ccx_elements['tri2plstrain'] = 'CPE6'
+ _ccx_elements['tri1plstrain'] = 'CPE3'
+ _ccx_elements['quad2plstrain'] = 'CPE8'
+ _ccx_elements['quad1plstrain'] = 'CPE4'
+ # set element ccxtype
+ self.ccxtype = _ccx_elements[estr]
+
+[docs] def ccx(self):
+ """Returns text string defining the element for a ccx input file."""
+ nids = [str(n.id) for n in self.node.values()]
+ val = str(self.id)+', '+', '.join(nids)
+ return val
+
+[docs] def get_name(self):
+ """Returns the element name based on id number."""
+ name = 'E%i' % (self.id)
+ return name
+
+ def __str__(self):
+ """Returns the element name."""
+ nids = [str(node.id) for node in self.nodes]
+ return '%s nodes: %s' % (self.get_name(), ','.join(nids))
+
+
+[docs]class Face(object):
+ """Makes an element face.
+
+ Args:
+ fnum (int): element face number
+ node1 (Node): first node
+ node2 (Node): second node
+ element (Element): parent element
+
+ Attributes:
+ id (int): element face number
+ nodes (list): list of nodes [node1, node2]
+ nmid (Node or None): midside node, defaults to None
+ element (Element): parent element
+ ext (bool): stores if face is external
+ """
+
+ def __init__(self, fnum, node1, node2, element):
+ self.id = fnum
+ self.nmid = None
+ self.ext = False
+ self.element = element
+ self.nodes = [node1, node2]
+ for node in self.nodes:
+ node.add_face(self)
+
+ def __hash__(self):
+ """Returns a hash = self.element.id*4, Lets us make sets of faces."""
+ # we assume the largest number of faces is 4, face id is 1-4, need 0-3
+ hval = self.element.id*4+(self.id-1)
+ return hval
+
+[docs] def length(self):
+ """Return the length of the face."""
+ point0 = geometry.Point(self.nodes[0].x, self.nodes[0].y)
+ point1 = geometry.Point(self.nodes[1].x, self.nodes[1].y)
+ vect = point1 - point0
+ return vect.length()
+
+[docs] def get_mnorm(self):
+ """Returns a list: [midpoint, normal_vector]
+
+ Returns:
+ midpoint (Point): face midpoint
+ norm_vect (Point): face normal vector, it is a unit vector
+ """
+ # return midpt, normal
+ point0 = geometry.Point(self.nodes[0].x, self.nodes[0].y)
+ point1 = geometry.Point(self.nodes[1].x, self.nodes[1].y)
+ midpt = point0 + point1
+ midpt = midpt*0.5
+ norm_vect = point1 - point0
+ norm_vect.rot_ccw_deg(90)
+ norm_vect.make_unit()
+ return [midpt, norm_vect]
+
+[docs] def set_ext(self):
+ """Sets the ext boolean flag to True. Face is external."""
+ self.ext = True
+
+[docs] def set_nmid(self, nmid):
+ """Sets the face midside node, nmid.
+
+ Args:
+ nmid (Point): passed node to set as the face midside node
+ """
+ self.nmid = nmid
+ nmid.add_face(self)
+ self.nodes = [self.nodes[0], self.nmid, self.nodes[-1]]
+
+ def __eq__(self, other):
+ """Compare this face to other face. True if node lists equal.
+
+ Args:
+ other (Face): the other face we are comparing
+
+ Returns:
+ bool: True if face node lists are equal, False otherwise.
+ """
+ if self.nodes == other.nodes:
+ return True
+ else:
+ return False
+
+ def __str__(self):
+ """Returns string listing object type, id number and nodes."""
+ node1 = self.nodes[0].id
+ node2 = self.nodes[1].id
+ fid = self.id
+ eid = self.element.id
+ tup = (fid, eid, node1, node2)
+ res = 'Face %i on element %i with nodes [%i,%i]' % tup
+ return res
+
+
+[docs]class Node(object):
+ """Makes a mesh node.
+
+ Args:
+ node_num (int): node number
+ x (float): x-coordinate
+ y (float): y-coordinate
+ z (float): z-coordinate
+
+ Attributes:
+ id (int): node number
+ x (float): x-coordinate
+ y (float): y-coordinate
+ z (float): z-coordinate
+ order (int): 1 or 2, 1 if it is a corner node, 2 if it is a midside node
+ elements (set): a set of elements that contain this node
+ faces (set): a set of faces that contain this node
+ """
+
+ def __init__(self, node_num, x, y, z):
+ self.id = node_num
+ self.x = x
+ self.y = y
+ self.z = z
+ self.order = 1
+ self.elements = set()
+ self.faces = set()
+
+[docs] def set_order(self, order):
+ """Sets the node order, must be 1 or 2.
+
+ Args:
+ order (int): node order, 1 is corner node, 2 is midside node
+ """
+ self.order = order
+
+[docs] def add_element(self, element):
+ """Adds an element to this node's set of elements.
+
+ Args:
+ element (Element): element to add to set self.elements
+ """
+ self.elements.add(element)
+
+[docs] def add_face(self, face):
+ """Adds a face to this node's set of faces.
+
+ Args:
+ face (Face): face to add to set self.faces
+ """
+ self.faces.add(face)
+
+[docs] def label(self, ax):
+ """Labels the node on the passed Matplotlib axis."""
+ ax.annotate(self.get_name(), (self.y, self.x))
+
+ def __hash__(self):
+ """Returns the node id as the hash."""
+ return self.id
+
+ def __eq__(self, other):
+ """Returns boolean of id equality to other Node.
+
+ Args:
+ other (Node): the node to compare to this node.
+
+ Returns:
+ boolean: True if ids are equal, False if they are not
+ """
+ if isinstance(other, Node):
+ return self.id == other.id
+ else:
+ return False
+
+[docs] def ccx(self):
+ """Returns a string defining the node for a ccx inp file."""
+ val = '%i, %f, %f, %f' % (self.id, self.x, self.y, self.z)
+ return val
+
+[docs] def get_name(self):
+ """Returns the string node name."""
+ name = 'N%i' % (self.id)
+ return name
+
+ def __str__(self):
+ """Returns string listing object type, id number and (x,y) coords."""
+ val = 'Node, id %i, (x,y)=(%f,%f)' % (self.id, self.x, self.y)
+ return val
+
+"""This module stores the Part class. It is used to make 2D parts.
+"""
+import numpy as np # needed for linspace on hole creation
+
+from . import base_classes
+from . import geometry #point, line, area
+
+[docs]class Part(base_classes.Idobj):
+ """This makes a part.
+
+ Args:
+ parent: parent FeaModel
+
+ Attributes:
+ __fea (FeaModel): parent FeaModel
+ points (list): list or part points, excludes arc centers
+ allpoints (list): list or part points, includes arc centers
+ lines (list): list of all Line and Arc that make the part
+ signlines (list): list of all SignLine and SignArc that make the part
+ __cursor (Point): location of a cursor drawing the part
+ __holemode (bool): if True, lines will be added to holes, otherwise,
+ they'll be added to areas
+ areas (list): list of Area that make up the part
+ left (list): list the parts leftmost lines, they must be vertical
+ right (list): list the parts rightmost lines, they must be vertical
+ top (list): list the parts top lines, they must be horizontal
+ bottom (list): list the parts bottom lines, they must be horizontal
+ center (Point): the area centroid of the part
+ nodes (list): list of part's nodes
+ elements (list): list of part's elements
+ """
+
+ def __init__(self, feamodel):
+ self.fea = feamodel
+ self.__cursor = geometry.Point(0, 0)
+ self.areas = [] # top area is the buffer
+ # make the buffer
+ area = self.fea.areas.append(geometry.Area(self, []))
+ self.areas.append(area)
+ base_classes.Idobj.__init__(self)
+ self.left = None
+ self.right = None
+ self.top = None
+ self.bottom = None
+ self.center = None
+ self.nodes = []
+ self.elements = []
+ self.__holemode = False
+ self.fea.parts.append(self)
+
+ def __hash__(self):
+ """Returns the item's id as its hash."""
+ return self.id
+
+ @property
+ def lines(self):
+ """Returns list of part lines."""
+ lines = set()
+ for area in self.areas:
+ lines.update(area.lines)
+ return list(lines)
+
+ @property
+ def signlines(self):
+ """Returns list of part signline and signarc."""
+ lines = set()
+ for area in self.areas:
+ lines.update(area.signlines)
+ return list(lines)
+
+ @property
+ def points(self):
+ """Returns list of part points, excludes arc centers."""
+ points = set()
+ lines = self.lines
+ for line in lines:
+ points.update(line.points)
+ return list(points)
+
+ @property
+ def allpoints(self):
+ """Returns list of part points, includes arc centers."""
+ points = set()
+ lines = self.lines
+ for line in lines:
+ points.update(line.allpoints)
+ return list(points)
+
+[docs] def get_item(self, item):
+ """"Returns the part's item(s) requested by the passed string.
+
+ Args:
+ item (str): string requesting item(s)
+
+ * Valid examples: 'P0', 'L0', 'left', 'A0'
+
+ Returns:
+ item(s) or None: If items are found they are returned
+
+ * If there is only one item it is returned
+ * If there are multiple items, they are returned as a list
+ * If no items are found None is returned
+ """
+ if item in ['left', 'right', 'top', 'bottom']:
+ items = getattr(self, item)
+ return items
+ elif item[0] == 'P':
+ # get point
+ items = self.points
+ num = int(item[1:])
+ res = [a for a in items if a.id == num]
+ return res[0]
+ elif item[0] == 'L':
+ # get line
+ items = self.signlines
+ num = int(item[1:])
+ items = [a for a in items if a.id == num]
+ return items[0]
+ elif item[0] == 'A':
+ # get area
+ items = self.areas
+ num = int(item[1:])
+ items = [a for a in items if a.id == num]
+ return items[0]
+ else:
+ print('Unknown item! Please pass the name of a point, line or area!')
+ return None
+
+[docs] def get_name(self):
+ """Returns the part name based on id number."""
+ return 'PART'+str(self.id)
+
+ def __set_side(self, side):
+ """Sets the part.side to a list of lines on that side of the part.
+
+ Used to set the part.left, part.right, part.top, part.bottom sides.
+
+ Args:
+ side (string): 'left', 'right', 'top','bottom'
+ """
+ # set index and axis, ind=0 is low side, ind=-1 is high side
+ inds = {'left':0, 'right':-1, 'top':-1, 'bottom':0}
+ axes = {'left':'y', 'right':'y', 'top':'x', 'bottom':'x'}
+ ind = inds[side]
+ axis = axes[side]
+
+ # loc = 'left', ind = 0, axis = 'y'
+ points = self.points
+ # sort the points low to high
+ points = sorted(points, key=lambda pt: getattr(pt, axis))
+ # store the target value
+ target_value = getattr(points[ind], axis)
+ res = []
+ lines = self.signlines
+ for sline in lines:
+ if isinstance(sline, geometry.SignLine):
+ pt_axis_vals = [getattr(pt, axis) for pt in sline.points]
+ pt_dist_vals = [abs(target_value - pt_axis_val) for pt_axis_val in pt_axis_vals]
+ if all([pt_dist_val < geometry.ACC for pt_dist_val in pt_dist_vals]):
+ # line is on the left side
+ res.append(sline)
+ setattr(self, side, res)
+
+[docs] def goto(self, x, y, holemode=False):
+ """Moves the part cursor to a location.
+
+ If that location has a point at it, use it.
+ If not, make a new point at that location.
+
+ Args:
+ x (float): x-coordinate of the point to go to
+ y (float): y-coordinate of the point to go to
+ holemode (bool): if True, we start drawing a hole here, otherwise
+ we start drawing an area
+
+ Returns:
+ self.__cursor (Point): returns the updated cursor point
+ """
+ [pnew, already_exists] = self.__make_get_pt(x, y)
+
+ if already_exists:
+ if self.areas[-1].closed == True:
+ # make a new area if the old area is already closed and we're
+ # going to an existing point
+ area = self.fea.areas.append(geometry.Area(self, []))
+ self.areas.append(area)
+
+ # return cursor
+ self.__cursor = pnew
+ self.__holemode = holemode
+ return self.__cursor
+
+ def __get_point(self, point):
+ """Returns point if found, None otherwise."""
+ points = self.allpoints
+ found_point = None
+ for apoint in points:
+ dist = point - apoint
+ dist = dist.length()
+ if dist < geometry.ACC:
+ # point already exists in part, use it
+ found_point = apoint
+ break
+ return found_point
+
+ def __make_get_pt(self, x, y):
+ """Gets a point if it exists, makes it if it doesn't. Returns the point.
+
+ Use this when you need a point made in the part, and you want to
+ use an extant point if one is available.
+
+ Args:
+ x (float): point x-coordinate
+ y (float): point y-coordinate
+
+ Returns:
+ list:
+ list[0]: Point
+ list[1]: boolean, True = the point already existed
+ """
+ thept = geometry.Point(x, y)
+ pfound = self.__get_point(thept)
+ pexists = True
+ if pfound == None:
+ pfound = thept
+ self.fea.register(pfound)
+ pexists = False
+ return [pfound, pexists]
+
+ def __calc_area_center(self):
+ """Calculates and returns the part and centroid Point.
+
+ Returns:
+ list: [area, Point]
+ """
+ val_list = []
+ for area in self.areas:
+ if area.closed == True:
+ aval, cval = area.area, area.center
+ val_list.append([aval, cval])
+ a_sum = sum([aval[0] for aval in val_list])
+ cxa_sum = sum([center.x*aval for [aval, center] in val_list])
+ cya_sum = sum([center.y*aval for [aval, center] in val_list])
+ cx_val = cxa_sum/a_sum
+ cy_val = cya_sum/a_sum
+ center = geometry.Point(cx_val, cy_val)
+ return [a_sum, center]
+
+ def __make_get_sline(self, lnew):
+ """Returns a signed line or arc, makes it if it needs to.
+
+ Args:
+ lnew (Line or Arc or SignLine or SignArc): Line or Arc to make
+
+ Returns:
+ list:
+ list[0]: SignLine or SignArc
+ list[1]: boolean, True = the line already existed
+ """
+ lpos = lnew.signed_copy(1)
+ lneg = lnew.signed_copy(-1)
+ # get part's signed lines
+ slines = self.signlines
+ lexists = False
+ for sline in slines:
+ if lpos == sline:
+ lexists = True
+ signline_new = sline
+ break
+ elif lneg == sline:
+ lexists = True
+ signline_new = sline.signed_copy(-1)
+ signline_new.edge = False
+ self.fea.register(signline_new)
+ break
+ else:
+ # fired when we haven't broken out of the loop, the line is new
+ if isinstance(lnew, geometry.SignLine):
+ lnew = geometry.Line(lnew.pt(0), lnew.pt(1))
+ self.fea.register(lnew)
+ lnew.save_to_points()
+ signline_new = lnew.signed_copy(1)
+ self.fea.register(signline_new)
+ signline_new.line.add_signline(signline_new)
+ return [signline_new, lexists]
+
+[docs] def draw_circle(self, center_x, center_y, radius, num_arcs=4):
+ """Draws a circle area and adds it to the part.
+
+ Args:
+ center_x (float): x-axis hole center
+ center_y (float): y-axis hole center
+ radius (float): hole radius
+ num_arcs (int): number of arcs to use, must be >= 3
+
+ Returns:
+ loop (geometry.LineLoop): a LineLoop list of SignArc
+ """
+ center = geometry.Point(center_x, center_y)
+ rvect = geometry.Point(0, radius)
+ start = center + rvect
+ self.goto(start.x, start.y)
+ angles = np.linspace(360/num_arcs,360,num_arcs, endpoint=True)
+ for ang in angles:
+ point = geometry.Point(0, radius).rot_ccw_deg(ang)
+ point = point + center
+ self.draw_arc(point.x, point.y, center.x, center.y)
+ loop = self.areas[-1].exlines
+ self.__update()
+ return loop
+
+[docs] def draw_hole(self, center_x, center_y, radius, num_arcs=4, filled=False):
+ """Makes a hole in the part.
+
+ Args:
+ center_x (float): x-axis hole center
+ center_y (float): y-axis hole center
+ radius (float): hole radius
+ num_arcs (int): number of arcs to use, must be >= 3
+ filled (bool): whether to fill the hole
+
+ * True: makes a new area in the part
+
+ Returns:
+ hole_lines (list or None): list of hole SignLine or SignArc
+
+ * Returns None if hole was not made.
+ """
+ center = geometry.Point(center_x, center_y)
+ area = self.__area_from_pt(center)
+ if area == None:
+ print("You can't make a hole here until there's an area here!")
+ return None
+ else:
+ # make points
+ rvect = geometry.Point(0, radius)
+ start = center + rvect
+ self.goto(start.x, start.y, holemode=True)
+ angles = np.linspace(360/num_arcs,360,num_arcs, endpoint=True)
+ for ang in angles:
+ point = geometry.Point(0, radius).rot_ccw_deg(ang)
+ point = point + center
+ self.draw_arc(point.x, point.y, center.x, center.y)
+ # make new area
+ if filled:
+ # reverse order, reverse sline directions, store in feamodel
+ slines = list(area.holes[-1])
+ slines.reverse()
+ slines = [sline.signed_copy(-1) for sline in slines]
+ slines = [self.__make_get_sline(sline)[0] for sline in slines]
+ anew = self.fea.areas.append(geometry.Area(self, slines))
+ self.areas.append(anew)
+ self.__update()
+ return area.holes[-1]
+
+[docs] def draw_arc_angle(self, degrees_ccw, center_x, center_y):
+ """Makes an arc and adds it to the part.
+
+ | Current point is the first arc point.
+ | degrees_ccw is the swept angle in degrees, counterclockwise
+ | (center_x, center_y) is the arc center
+
+ | Degrees: Traversed angle of arc must be < 180 degrees
+
+ Args:
+ degrees_ccw (float): arc swept angle in degrees, counterclockwise
+ center_x (float): arc center point x-coordinate
+ center_y (float): arc center point y-coordinate
+
+ Returns:
+ list: [arc, arc_start_point, arc_end_point]
+ """
+ center = geometry.Point(center_x, center_y)
+ radius_vector = self.__cursor - center
+ radius_vector.rot_ccw_deg(degrees_ccw)
+ end = center + radius_vector
+ return self.draw_arc(end.x, end.y, center_x, center_y)
+
+[docs] def draw_arc(self, end_x, end_y, center_x, center_y):
+ """Makes an arc and adds it to the part.
+
+ | Current point is the first arc point.
+ | (end_x, end_y) is the end point
+ | (center_x, center_y) is the arc center
+
+ | Degrees: Traversed angle of arc must be < 180 degrees
+ | Radians: Traversed angle of arc must be < Pi
+
+ Args:
+ end_x (float): arc end point x-coordinate
+ end_y (float): arc end point y-coordinate
+ center_x (float): arc center point x-coordinate
+ center_y (float): arc center point y-coordinate
+
+ Returns:
+ list: [arc, arc_start_point, arc_end_point]
+ """
+ pold = self.__cursor
+ # make arc center point
+ ctr = self.__make_get_pt(center_x, center_y)[0]
+ # make arc end point
+ self.__cursor = self.__make_get_pt(end_x, end_y)[0]
+
+ # make arc
+ arc = self.__make_get_sline(geometry.Arc(pold, self.__cursor, ctr))[0]
+
+ if self.__holemode:
+ area = self.__area_from_pt(self.__cursor)
+ if area != None:
+ closed = area.add_hole_sline(arc)
+ if closed:
+ self.__holemode = False
+ else:
+ print('You must have a closed area here before making a hole!')
+ else:
+ self.areas[-1].add_sline(arc)
+ return [arc, pold, self.__cursor]
+
+[docs] def draw_line_delta(self, delta_x, delta_y):
+ """Draws a line a relative distance, and adds it to the part.
+
+ Args:
+ delta_x (float): x-axis delta distance to draw the line
+ delta_y (float): y-axis delta distance to draw the line
+
+ Returns:
+ list: [line, point_start, point_end]
+ """
+ x = self.__cursor.x + delta_x
+ y = self.__cursor.y + delta_y
+ return self.draw_line_to(x, y)
+
+[docs] def draw_line_rad(self, dx_rad):
+ """Draws a line a relative radial distance, and adds it to the part.
+
+ Args:
+ dx_rad (float): x-axis delta distance to draw the line
+
+ Returns:
+ list: [line, point_start, point_end]
+ """
+ return self.draw_line_delta(dx_rad, 0.0)
+
+[docs] def draw_line_ax(self, dy_ax):
+ """Draws a line a relative axial distance, and adds it to the part.
+
+ Args:
+ dy_ax (float): y-axis delta distance to draw the line
+
+ Returns:
+ list: [line, point_start, point_end]
+ """
+ return self.draw_line_delta(0.0, dy_ax)
+
+[docs] def draw_line_to(self, x, y):
+ """Draws a line to the given location, and adds it to the part.
+
+ Args:
+ x (float): x-axis coordinate of the end point
+ y (float): y-axis coordinate of the end point
+
+ Returns:
+ list: [SignLine, point_start, point_end]
+ """
+ pold = self.__cursor
+ self.__cursor = self.__make_get_pt(x, y)[0]
+ sline = self.__make_get_sline(geometry.Line(pold, self.__cursor))[0]
+
+ if self.__holemode:
+ area = self.__area_from_pt(self.__cursor)
+ if area != None:
+ closed = area.add_hole_sline(sline)
+ if closed:
+ self.__holemode = False
+ self.__update()
+ else:
+ print('You must have a closed area here before making a hole!')
+ else:
+ # drawing in last area
+ self.areas[-1].add_sline(sline)
+ # check for closure of the area
+ if self.areas[-1].closed:
+ self.__update()
+ return [sline, pold, self.__cursor]
+
+ def __get_maxlength(self):
+ """Returns the max distance between points in the part."""
+ points = self.points
+ maxlen = 0.0
+ # loop through points checking dist to next point
+ for ind, point_1 in enumerate(points[:-1]):
+ for point_2 in points[ind:]:
+ vect = point_1 - point_2
+ dist = vect.length()
+ if dist > maxlen:
+ maxlen = dist
+ return maxlen
+
+ def __area_from_pt(self, point):
+ """Returns the area that the point is inside.
+
+ Args:
+ point (Point): the point we are asking about
+
+ Returns:
+ Area or None:
+ Area is the found area
+ None is returned if the point is not in one of this part's areas
+ """
+ for area in self.areas:
+ if area.contains_point(point):
+ return area
+ return None
+
+[docs] def fillet_lines(self, line1, line2, radius):
+ """Fillets the given lines in the part.
+
+ This inserts an arc in the part tangent to the two given lines.
+
+ Args:
+ line1 (SignLine): line that the arc starts on, arc is tangent
+ line2 (SignLine): line that the arc ends on, arc is tangent
+ radius (float): arc radius size
+
+ Returns:
+ list: [arc, start_point, end_point]
+ """
+ # check if the lines are touching
+ if not line1.line.touches(line2.line):
+ print('ERROR: Cannot fillet! Lines must touch!')
+ return
+
+ if line1.line.pt(1) == line2.line.pt(0):
+ first_line = line1
+ second_line = line2
+ elif line2.line.pt(1) == line1.line.pt(0):
+ first_line = line2
+ second_line = line1
+ else:
+ print('ERROR: Sign lines must both be going in CW or CCW '
+ 'direction. The two passed lines are going in '
+ 'different directions. Unable to fillet them.')
+ return
+
+ tmp = self.__cursor
+ # offset the lines, assuming area is being traced clockwise
+ # get the intersection point
+ magnitude = radius
+ l1_off = first_line.offset(magnitude)
+ l2_off = second_line.offset(magnitude)
+ ctrpt = l1_off.intersects(l2_off)
+ if ctrpt == None:
+ # flip the offset direction if lines don't intersect
+ magnitude = -radius
+ l1_off = first_line.offset(magnitude)
+ l2_off = second_line.offset(magnitude)
+ ctrpt = l1_off.intersects(l2_off)
+
+ # now we have an intersecting point
+ p1_new = first_line.arc_tang_intersection(ctrpt, magnitude)
+ p2_new = second_line.arc_tang_intersection(ctrpt, magnitude)
+ rempt = first_line.pt(1)
+
+ p1_new = self.__make_get_pt(p1_new.x, p1_new.y)[0]
+ ctrpt = self.__make_get_pt(ctrpt.x, ctrpt.y)[0]
+ p2_new = self.__make_get_pt(p2_new.x, p2_new.y)[0]
+
+ # make the new arc
+ arc = self.__make_get_sline(geometry.Arc(p1_new, p2_new, ctrpt))[0]
+
+ # put the arc in the right location in the area
+ area = first_line.lineloop.parent
+ area.line_insert(first_line, arc)
+ print('Arc inserted into area %i' % (area.id))
+
+ # edit the adjacent lines to replace the removed pt
+ first_line.set_pt(1, arc.pt(0))
+ second_line.set_pt(0, arc.pt(1))
+ # del old pt, store new points for the arc
+ self.fea.points.remove(rempt)
+
+ # reset the cursor to where it should be
+ self.__cursor = tmp
+ return [arc, arc.pt(0), arc.pt(1)]
+
+[docs] def fillet_all(self, radius):
+ """Fillets all external lines not within 10 degrees of tangency
+
+ Args:
+ radius (float): the fillet radius to use
+
+ Returns:
+ arcs (list): list of SignArc
+ """
+ pairs = []
+ for area in self.areas:
+ for ind, sline in enumerate(area.exlines):
+ prev_sline = area.exlines[ind-1]
+ this_point = sline.pt(0)
+ if len(this_point.lines) == 2:
+ # only fillet lines that are not shared by other areas
+ if (isinstance(sline, geometry.SignLine)
+ and isinstance(prev_sline, geometry.SignLine)):
+ # only fillet lines
+ perp1 = prev_sline.get_perp_vec(this_point)
+ perp2 = sline.get_perp_vec(this_point)
+ ang = perp1.ang_bet_deg(perp2)
+ is_tangent = (-10 <= ang <= 10)
+ if is_tangent == False:
+ pairs.append([prev_sline, sline])
+ arcs = []
+ for pair in pairs:
+ arc = self.fillet_lines(pair[0], pair[1], radius)[0]
+ arcs.append(arc)
+ return arcs
+
+[docs] def label(self, axis):
+ """Labels the part on a Matplotlib axis
+
+ Args:
+ axis (Matplotlib Axis): Matplotlib Axis
+ """
+ axis.text(self.center.y, self.center.x, self.get_name(),
+ ha='center', va='center')
+
+[docs] def plot(self, axis, label=True, color='yellow'):
+ """Plots the part on the passed Matplotlib axis.
+
+ Args:
+ axis (Matplotlib axis): the axis we will plot on
+ label (bool): True displays the part label
+ color (tuple): r,g,b,a matplotlib color tuple
+ """
+ patches = []
+ for area in self.areas:
+ if area.closed:
+ patches.append(area.get_patch())
+ for patch in patches:
+ patch.set_color(color)
+ axis.add_patch(patch)
+ # apply the label
+ if label:
+ self.label(axis)
+
+ def __cut_line(self, point, line):
+ """Cuts the passed line at the passed point.
+
+ The passed line is cut into two lines. All areas that included the
+ original line are updated.
+
+ Args:
+ line (Line or Arc): the line to cut, must be Line or Arc
+ point (Point): the location on the line that we will cut it
+
+ Returns:
+ list: [pnew, lnew]
+ pnew: the new point we created to cut the original line
+ lnew: the new line we created, the end half of the orignal line
+ """
+ pnew = self.__make_get_pt(point.x, point.y)[0]
+ if point.id != -1:
+ # if passed point already exists, use it
+ pnew = point
+ pend = line.pt(1)
+ line.set_pt(1, pnew) # shortens the line
+ new_prim = geometry.Line(pnew, pend)
+ if isinstance(line, geometry.Arc):
+ new_prim = geometry.Arc(pnew, pend, line.actr)
+ new_sline = self.__make_get_sline(new_prim)[0]
+ # insert the new line into existing areas
+ is_line = isinstance(line, geometry.Line) or isinstance(line, geometry.Arc)
+ is_sline = isinstance(line, geometry.SignLine) or isinstance(line, geometry.SignArc)
+ print('Cutting line (is_line, is_sline, signlines) (%s, %s, %i)' % (is_line, is_sline, len(line.signlines)))
+ slines = line.signlines
+ for sline in slines:
+ area = sline.lineloop.parent
+ if sline.sign == 1:
+ # cutting line in clockwise area, where line is pos
+ area.line_insert(sline, new_sline)
+ elif sline.sign == -1:
+ # cutting line in clockwise area, where line is neg
+ rev_sline = self.__make_get_sline(new_sline.signed_copy(-1))[0]
+ area.line_insert(sline, rev_sline, after=False)
+ return [pnew, new_sline]
+
+ def __cut_area(self, area, start_pt, end_pt):
+ """Cuts the part area from start_pt to end_pt."""
+ # select the line portions that define the areas
+ # list[:low] excludes low index
+ # list [high:] includes high index
+ # we want the line which start with the point
+ lpre_start = area.line_from_startpt(start_pt)
+ lpre_end = area.line_from_startpt(end_pt)
+ if lpre_start == None or lpre_end == None:
+ self.fea.plot_geometry()
+ print(area.exlines)
+ istart = area.exlines.index(lpre_start)
+ iend = area.exlines.index(lpre_end)
+ low = min(istart, iend)
+ high = max(istart, iend)
+
+ # lists of lines for areas
+ beg = area.exlines[:low]
+ mid = area.exlines[low:high]
+ end = area.exlines[high:]
+
+ # make cut line for [beg + cut + end] area
+ start_pt = mid[0].pt(0)
+ end_pt = mid[-1].pt(1)
+ fwd = geometry.Line(start_pt, end_pt)
+ rev = geometry.Line(end_pt, start_pt)
+
+ # update existing area
+ cline = self.__make_get_sline(fwd)[0]
+ alist_curr = beg + [cline] + end
+ area.update(alist_curr)
+
+ # make new area
+ cline_rev = self.__make_get_sline(rev)[0]
+ alist_other = mid + [cline_rev]
+ anew = geometry.Area(self, alist_other)
+ self.fea.register(anew)
+ self.areas.append(anew)
+
+ # fix holes
+ self.__store_holes()
+
+ def __merge_hole(self, area, start_pt, end_pt):
+ """Merges the hole into its area with a line between passed points."""
+ # line will be drawn from start point on exlines to end point on hole
+ hole_points = area.holepoints
+ if start_pt in hole_points:
+ tmp = start_pt
+ start_pt = end_pt
+ end_pt = tmp
+ lpre_start = area.line_from_startpt(start_pt)
+ hole_line = area.line_from_startpt(end_pt)
+ if lpre_start == None or hole_line == None:
+ self.fea.plot_geometry()
+ ind = area.exlines.index(lpre_start)
+
+ # store sections of the area
+ beg = area.exlines[:ind]
+ end = area.exlines[ind:]
+ thehole = None
+ mid = []
+ for hole in area.holes:
+ for sline in hole:
+ if sline == hole_line:
+ ind = hole.index(sline)
+ mid = hole[ind:] + hole[:ind]
+ thehole = hole
+ break
+ if mid != []:
+ break
+
+ fwd = geometry.Line(start_pt, end_pt)
+ fwd_sline = self.__make_get_sline(fwd)[0]
+ rev_sline = fwd_sline.signed_copy(-1)
+ self.fea.register(rev_sline)
+ rev_sline.line.add_signline(rev_sline)
+ alist_curr = beg + [fwd_sline] + mid + [rev_sline] + end
+ area.holes.remove(thehole)
+ area.update(alist_curr)
+
+ def __get_cut_line(self, cutline):
+ """Returns a cut line beginning and ending on the part."""
+ # find all intersections
+ lines = self.lines
+ points = set()
+
+ # add line intersections
+ for line in lines:
+ newpt = line.intersects(cutline)
+ if newpt != None:
+ points.add(newpt)
+
+ # loop through intersection points, storing distance
+ points = list(points)
+ for (ind, point) in enumerate(points):
+ dist = point - cutline.pt(0)
+ dist = dist.length()
+ pdict = {'dist': dist, 'point': point}
+ points[ind] = pdict
+
+ # sort the points by dist, lowest to highest, return first cut
+ points = sorted(points, key=lambda k: k['dist'])
+ start = points[0]['point']
+ end = points[1]['point']
+ new_cut = geometry.Line(start, end)
+ return new_cut
+
+ def __cut_with_line(self, cutline, debug):
+ """Cuts the part using the passed line.
+
+ Args:
+ cutline (Line): line to cut the area with
+ debug (list): bool for printing, bool for plotting after every cut
+ """
+ # find all intersections
+ lines = self.lines
+ points = set()
+
+ # add line intersections
+ for line in lines:
+ if debug[0]:
+ print('Checking X between %s and cutline' % line.get_name())
+ newpt = line.intersects(cutline)
+ if debug[0]:
+ print(' Intersection: %s' % newpt)
+ if newpt != None:
+ points.add(newpt)
+
+ # loop through intersection points, storing distance and lines to cut
+ points = list(points)
+ for (ind, point) in enumerate(points):
+ dist = point - cutline.pt(0)
+ dist = dist.length()
+ pdict = {'dist': dist}
+ realpt = self.__get_point(point)
+ # we only want to store lines to cut here
+ if realpt == None or realpt.arc_center == True:
+ # we could have an existing arc center on a line that needs to
+ # be cut
+ if realpt == None:
+ realpt = point
+ for line in lines:
+ point_on_line = line.coincident(realpt)
+ if point_on_line and point not in line.points:
+ pdict['line'] = line
+ break
+ pdict['point'] = realpt
+ points[ind] = pdict
+
+ # sort the points by dist, lowest to highest
+ points = sorted(points, key=lambda k: k['dist'])
+ if debug[0]:
+ print('==================================')
+ print('Points on the cutline!------------')
+ for pdict in points:
+ print(pdict['point'])
+ print(' dist %.3f' % pdict['dist'])
+ if 'line' in pdict:
+ print(' X cut line: %s' % pdict['line'])
+ print('==================================')
+
+ # loop through the points cutting areas
+ for ind in range(len(points)):
+ pdict = points[ind]
+ start_pt = pdict['point']
+ if 'line' in pdict:
+ # cut the line and point to the real new point
+ print('Cut through line %s' % pdict['line'].get_name())
+ pnew = self.__cut_line(start_pt, pdict['line'])[0]
+ points[ind]['point'] = pnew
+ start_pt = pnew
+ end_pt = None
+ pavg = None
+ area = None
+ if ind > 0:
+ # find the area we're working on
+ end_pt = points[ind-1]['point']
+ pavg = start_pt + end_pt
+ pavg = pavg*0.5
+ area = self.__area_from_pt(pavg)
+
+ if area == None:
+ # stop cutting if we are trying to cut through a holes
+ print('No area found at point avg, no cut made')
+ break
+ start_hole = start_pt in area.holepoints
+ end_hole = end_pt in area.holepoints
+ if start_hole and end_hole and area != None:
+ print('Trying to join holes, no cut made')
+ break
+ # stop cutting if we are trying to join holes
+ if end_hole == True or start_hole == True:
+ print('Merging hole in %s' % area.get_name())
+ self.__merge_hole(area, start_pt, end_pt)
+ else:
+ print('Cutting %s' % area.get_name())
+ self.__cut_area(area, start_pt, end_pt)
+ if debug[1]:
+ self.fea.plot_geometry()
+
+ def __store_holes(self):
+ """Puts all holes in their correct areas"""
+ holes = []
+ for area in self.areas:
+ holes += area.holes
+ for hole in holes:
+ hole_area = hole.parent
+ for area in self.areas:
+ is_inside = hole.inside(area.exlines)
+ if is_inside == True:
+ if area != hole_area:
+ # delete the hole from the old area, move it to the new
+ hole.set_parent(area)
+ hole_area.holes.remove(hole)
+ hole_area.close()
+ area.holes.append(hole)
+ area.close()
+ afrom, ato = hole_area.get_name(), area.get_name()
+ print('Hole moved from %s to %s' % (afrom, ato))
+
+ def __vect_to_line(self, point, cvect):
+ """Returns a cutting line at a given point and cutting vector.
+
+ Args:
+ point (Point): the location we are cutting from
+ cvect (Point): the vector direction of the cut from pt
+
+ Returns:
+ cutline (Line): cut line
+ """
+ cvect.make_unit()
+ vsize = self.__get_maxlength()
+ endpt = point + cvect*vsize
+ cutline = geometry.Line(point, endpt)
+ cutline = self.__get_cut_line(cutline)
+ return cutline
+
+ def __chunk_area(self, area, mode, exclude_convex, debug):
+ """Cuts the passed area into regular smaller areas.
+
+ The cgx mesher only accepts areas which are 3-5 sides
+ so one may need to call this before using that mesher.
+ Cuts are made perpendicular to tangent points or at
+ internal corners.
+ At internal corners two perpendicular cuts are made.
+
+ Args:
+ area (Area): the area to cut into smaller areas
+ mode (str): 'both', 'holes' or 'ext' chunks the area using the
+ points form this set. See part.chunk
+ exclude_convex (bool): If true exclude cutting convex tangent points
+ debug (list): bool for printing, bool for plotting after every cut
+ """
+ # store the cuts first, then cut after
+ cuts = [] # each item is a dict with a pt and vect in it
+ loops = []
+ cut_point_sets = []
+ if mode == 'holes':
+ loops = area.holes
+ elif mode == 'ext':
+ loops = [area.exlines]
+ elif mode == 'both':
+ loops = area.holes + [area.exlines]
+ for loop in loops:
+ for ind, line in enumerate(loop):
+ line_pre = loop[ind-1]
+ line_post = line
+ point = line_pre.pt(1)
+ perp1 = line_pre.get_perp_vec(point)
+ perp2 = line_post.get_perp_vec(point)
+ #tan1 = line_pre.get_tan_vec(point)
+ #tan2 = line_post.get_tan_vec(point)
+ # flip these vectors later to make them cut the area(s)
+ ang = perp1.ang_bet_deg(perp2)
+ cut = {}
+ make_cut = True
+ pre_arc = isinstance(line_pre, geometry.SignArc)
+ post_arc = isinstance(line_post, geometry.SignArc)
+ if pre_arc or post_arc:
+ if pre_arc and post_arc:
+ if (line_pre.concavity == 'convex'
+ and line_post.concavity == 'convex'
+ and exclude_convex == True):
+ make_cut = False
+ else:
+ # only one is an arc
+ arc = line_pre
+ if post_arc:
+ arc = line_post
+ if arc.concavity == 'convex' and exclude_convex == True:
+ make_cut = False
+ is_tangent = (-10 <= ang <= 10)
+ is_int_corner = (45 <= ang <= 135)
+ """
+ print('-------------------')
+ print('%s' % point)
+ print('Angle is %.3f' % ang)
+ print('Make cut %s' % make_cut)
+ print('is_tangent %s' % is_tangent)
+ print('is_int_corner %s' % is_int_corner)
+ """
+ if is_tangent:
+ if make_cut == True:
+ # tangent
+ cut = {'pt':point, 'vect':perp1*-1}
+ cut_line = self.__vect_to_line(cut['pt'], cut['vect'])
+ pset = set(cut_line.points)
+ if pset not in cut_point_sets:
+ cut_point_sets.append(pset)
+ cut['line'] = cut_line
+ cuts.append(cut)
+ elif is_int_corner:
+ # internal corner
+ cut = {'pt':point, 'vect':perp1*-1}
+ cut_line = self.__vect_to_line(cut['pt'], cut['vect'])
+ pset = set(cut_line.points)
+ if pset not in cut_point_sets:
+ cut_point_sets.append(pset)
+ cut['line'] = cut_line
+ cuts.append(cut)
+ cut = {'pt':point, 'vect':perp2*-1}
+ cut_line = self.__vect_to_line(cut['pt'], cut['vect'])
+ pset = set(cut_line.points)
+ if pset not in cut_point_sets:
+ cut_point_sets.append(pset)
+ cut['line'] = cut_line
+ cuts.append(cut)
+ elif ang < 0:
+ # external corner
+ # do not split these
+ pass
+
+ # do the cuts
+ for cut in cuts:
+ print('--------------------')
+ print('Cut point:', cut['pt'].get_name())
+ print('Cut line:', cut['line'])
+ self.__cut_with_line(cut['line'], debug)
+
+[docs] def chunk(self, mode='both', exclude_convex = True, debug=[0, 0]):
+ """Chunks all areas in the part.
+
+ Args:
+ mode (str): area chunking mode
+
+ - 'both': cuts areas using holes and exterior points
+ - 'holes': cut areas using holes points only
+ - 'ext': cut areas using exterior points only
+
+ exclude_convex (bool): If true exclude cutting convex tangent points
+ """
+ for area in self.areas:
+ if area.closed:
+ min_sides = 5
+ has_holes = len(area.holes) > 0
+ ext_gr = len(area.exlines) > min_sides
+ both_false = (has_holes == False and ext_gr == False)
+ if mode == 'holes' and has_holes:
+ self.__chunk_area(area, mode, exclude_convex, debug)
+ elif (mode == 'both'
+ and (has_holes or ext_gr or not exclude_convex)):
+ self.__chunk_area(area, mode, exclude_convex, debug)
+ elif mode == 'ext' and (ext_gr or not exclude_convex):
+ self.__chunk_area(area, mode, exclude_convex, debug)
+ else:
+ aname = area.get_name()
+ val = 'Area %s was not chunked because it had' % aname
+ adder = ''
+ if mode == 'both' and both_false:
+ adder = '<= %i lines and no holes.' % min_sides
+ elif has_holes == False and (mode in ['both', 'holes']):
+ adder = 'no holes.'
+ elif ext_gr == False and (mode in ['both', 'ext']):
+ adder = '<= %i lines.' % min_sides
+ print('%s %s' % (val, adder))
+
+ # store the left, right, top, and bottom lines
+ self.__update()
+
+ def __update(self):
+ """Updates the left, right, top, bottom sides and area and center"""
+ self.__set_side('left')
+ self.__set_side('right')
+ self.__set_side('top')
+ self.__set_side('bottom')
+ self.area, self.center = self.__calc_area_center()
+
+ def __str__(self):
+ """Returns string listing object type, id number and name."""
+ val = 'Part, id=%i name=%s' % (self.id, self.get_name())
+ return val
+
+"""This module stores the Problem class which is used to solve different types
+of analysis.
+"""
+
+import os # used to see if there is a results file
+import subprocess # used to launch ccx solver
+
+from . import environment
+from . import base_classes
+from . import results_file
+
+[docs]class Problem(base_classes.Idobj):
+ """Makes a problem which can be analyzed with Calculix ccx.
+
+ Args:
+ - feamodel (FeaModel): the parent FeaModel
+ - problem_type (str): model type, options:
+ -- 'struct': structural
+ - fname (str): file prefix for the problem .inp and results files
+ If value is '' it will default to the project name of the FeaModel
+
+ Attributes:
+ fea (FeaModel): parent FeaModel
+ __ptype (str): problem type, options:
+
+ - 'struct': structural
+ solved (bool): boolean storing whether or not the analysis was solved
+ fname (str): the problem database and results file prefix
+ rfile (Results_File):
+ Results_File is loaded in after the model has been solved
+ """
+ def __init__(self, feamodel, problem_type, fname=''):
+ self.fea = feamodel
+ self.__ptype = problem_type
+ self.solved = False
+ if fname == '':
+ fname = self.fea.fname
+ self.fname = fname
+ self.rfile = results_file.ResultsFile(self)
+ base_classes.Idobj.__init__(self)
+ self.fea.problems.append(self)
+
+ @staticmethod
+ def __get_ntxt(nodes):
+ """Returns list of strings defining all nodes.
+
+ Args:
+ nodes (list): list of all nodes
+ """
+ res = []
+ res.append('*NODE, NSET=nodes')
+ for node in nodes:
+ res.append(node.ccx())
+ print('INFO: %i nodes' % len(nodes))
+ return res
+
+ def __get_etxt(self, elements):
+ """Returns list of strings defining all elements.
+
+ Args:
+ elements (list): list of all elements
+ """
+ res = []
+ ccxtypes = set([e.ccxtype for e in elements])
+ eall_written = False
+ eall = []
+ for ccxtype in ccxtypes:
+ setname = ccxtype
+ if len(ccxtypes) == 1:
+ setname = 'EAll'
+ eall_written = True
+ res.append('*ELEMENT, TYPE=%s, ELSET=%s' % (ccxtype, setname))
+ elist = [e for e in elements if e.ccxtype == ccxtype]
+ print('INFO: %i %s elements' % (len(elist), ccxtype))
+ for element in elist:
+ res.append(element.ccx())
+ eall += elist
+ print('INFO: %i total elements' % len(eall))
+ if eall_written == False:
+ tmp = self.__get_eset('EALL', eall)
+ res += tmp
+ return res
+
+ @staticmethod
+ def __get_ctxt(components):
+ """Returns list of strings defining all components.
+
+ Args:
+ components (list): list of all components
+ """
+ res = []
+ for comp in components:
+ res += comp.ccx()
+ return res
+
+ @staticmethod
+ def __get_eset(name, elements):
+ """Returns list of strings defining components of elements.
+
+ Args:
+ name (str): component name
+ elements (list): list of component elements
+ """
+ res = []
+ items_per_line = 6
+ res.append('*ELSET,ELSET='+name)
+ grouped_els = base_classes.chunk_list(elements, items_per_line)
+ for group in grouped_els:
+ item_ids = [str(e.id) for e in group]
+ line = ', '.join(item_ids)
+ if group != grouped_els[-1]:
+ line += ','
+ res.append(line)
+ return res
+
+ @staticmethod
+ def __fix_line(line, fstr):
+ """Fixes the line and returns a fixed line.
+
+ Args:
+ line (str): line to fix
+ fstr (str): format string to use
+ """
+ items = fstr.split(',')
+ node_num_size = int(items[2][1:])
+ pre_len = 3 + node_num_size
+ # we are missing + prefix on some numbers, some are too long
+ start_str = line[:pre_len]
+ end_str = line[pre_len:]
+ fields = end_str.count('E')
+ low_ind = 0
+ values = []
+ while len(values) < fields:
+ high_ind = end_str.find('E', low_ind)+5
+ value = float(end_str[low_ind:high_ind])
+ values.append(value)
+ low_ind = high_ind
+ values = ['%12.5e' % val for val in values]
+ end_str = ''.join(values)
+ new_str = start_str + end_str + '\n'
+ return new_str
+
+ def __fix_frd(self):
+ """Fixes the frd file on win32 systems. Text formatting fixed
+
+ On win32 Calculix, results file formatting is not written correctly.
+ Nodal results and stresses are not written as fixed length fields.
+ """
+ if 'win32' in environment.CCX:
+ frd_file = self.fname+'.frd'
+ lines = []
+ try:
+ with open(frd_file, "r") as infile:
+ lines = infile.readlines()
+ numlines = len(lines)
+ ind = 0
+ fix = False
+ fstr = ''
+ while ind < numlines:
+ line = lines[ind]
+ if '1PSTEP' in line:
+ # we are in a results block
+ ind += 1
+ line = lines[ind]
+ format_ = int(line.split()[-1])
+ if format_ == 0:
+ fstr = "1X,I2,I5,6E12.5"
+ elif format_ == 1:
+ fstr = "1X,I2,I10,6E12.5"
+ ind += 1
+ line = lines[ind]
+ if 'DISPR' in line or 'FORCR' in line:
+ ind += 5
+ else:
+ ind += 7
+ fix = True
+ line = lines[ind]
+ if line[:3] == ' -3':
+ fix = False
+ if fix:
+ lines[ind] = self.__fix_line(line, fstr)
+ ind += 1
+ with open(frd_file, "w") as outfile:
+ outfile.writelines(lines)
+ print('File %s had its formatting fixed!' % frd_file)
+ except IOError:
+ print('Error reading .frd file!')
+
+[docs] def solve(self):
+ """Solves the model in Calculix ccx."""
+ inp = []
+
+ # store what results we'll be outputting for each type of analysis
+ out_el = {}
+ out_el['struct'] = 'E,S' # strain, stress
+ out_node = {}
+ out_node['struct'] = 'RF,U' # reaction forces, displacement
+
+ if self.__ptype == 'struct':
+ # store nodes, elements, components
+ box = {}
+ box['nodes'] = []
+ box['elements'] = []
+ box['components'] = set()
+
+ # store all loads in the parts in this model
+ load_dict = self.fea.loads
+
+ # store all nodes, elements, and part element sets
+ box['nodes'] = self.fea.view.nodes
+ box['elements'] = self.fea.view.elements
+
+ # store all node and element components
+ for time in load_dict:
+ for load in load_dict[time]:
+ if load.ltype not in ['press', 'press_fluid']:
+ box['components'].add(load.comp)
+ for contact in self.fea.contacts:
+ box['components'].add(contact.master_comp)
+ box['components'].add(contact.slave_comp)
+
+ box['nodes'] = self.__get_ntxt(box['nodes'])
+ box['elements'] = self.__get_etxt(box['elements'])
+ box['components'] = self.__get_ctxt(box['components'])
+
+ # add text definition for nodes, elelents, components
+ inp += box['nodes']+['']
+ inp += box['elements']+['']
+ inp += box['components']+['']
+
+ # read in all materials
+ for matl in self.fea.matls:
+ inp += matl.ccx()
+
+ # write all steps and loads
+ for time in sorted(load_dict.keys()):
+ if time == 0.0:
+ # this is for thicknesses and materials
+ for load in load_dict[time]:
+ inp += load.ccx()
+ for surf_interaction in self.fea.surfints:
+ inp += surf_interaction.ccx()
+ for contact in self.fea.contacts:
+ inp += contact.ccx()
+ else:
+ # only write times >= 1
+ inp.append('*STEP')
+ inp.append('*STATIC')
+
+ for load in load_dict[time]:
+ inp += load.ccx()
+
+ # make output frd file for cgx
+ inp.append('*EL FILE')
+ inp.append(out_el[self.__ptype])
+ inp.append('*NODE FILE')
+ inp.append(out_node[self.__ptype])
+
+ # make output dat file for integration point results
+ inp.append('*EL PRINT,ELSET=EALL')
+ inp.append('S')
+
+ # end step
+ inp.append('*END STEP')
+
+ # write CCX inp file to the local directory
+ fname = self.fname+'.inp'
+ with open(fname, 'w') as outfile:
+ for line in inp:
+ #print (line)
+ outfile.write(line+'\n')
+ print('File: %s was written' % (fname))
+
+ # run file
+ runstr = "%s %s" % (environment.CCX, self.fname)
+ print(runstr)
+ subprocess.check_call(runstr, shell=True)
+ print('Solving done!')
+
+ # select the probem's parts and load the results file
+ frd_file = self.fname+'.frd'
+ if os.path.isfile(frd_file):
+ self.solved = True
+ self.__fix_frd()
+ self.rfile.load()
+ else:
+ print('ERROR: results .frd file was not written!')
+
+"""This module stores the Results_File class."""
+
+import collections
+import math # used for metric number conversion
+import os #need to check if results file exists
+import re # used to get info from frd file
+import subprocess # used to check ccx version
+
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import matplotlib.cm as cmx
+import numpy as np
+from numpy.lib.polynomial import roots # need to find S1-S3
+from numpy.core.function_base import linspace # need to make contours
+
+
+from . import base_classes # needed for RESFIELDS
+from . import environment
+from . import mesh
+
+CMAP = 'jet'
+
+[docs]class ResultsFile(object):
+ """Makes a results file.
+
+ Args:
+ problem (Problem): problem that was solved
+
+ Attributes:
+ __problem (Problem): parent problem
+ __steps (list): a list of float time steps
+ __results (dict): a dict storing the results data
+ results[step]['node'][nnum][field] --> value
+ results[step]['element'][enum]['avg'][field] --> value
+ results[step]['element'][enum]['max'][field] --> value
+ results[step]['element'][enum]['min'][field] --> value
+ results[step]['element'][enum]['ipoints'][ipnum][field] --> value
+ field = 'ux' or 'Seqv' or 'ey' etc.
+ __time (float): current time we are looking at, defaults to -1
+ When a file is loaded in, the first time step is loaded.
+ """
+
+ def __init__(self, problem):
+ self.__problem = problem
+ self.__steps = [] # this stores a list of time steps
+ self.__results = {} # stores results, nested dicts
+ self.__time = -1
+ if self.__problem.solved:
+ self.load()
+
+ @property
+ def steps(self):
+ """Returns a list of loaded time steps.
+
+ Note: this is read only, you can not assign a value to it.
+ """
+ return self.__steps
+
+ @property
+ def time(self):
+ """Returns the current time (float) in the results file.
+
+ Note: this is read only, you can not assign a value to it.
+ """
+ return self.__time
+
+ @staticmethod
+ def __metric_num(number, sig_figs=3, sci=False):
+ """Returns string of number, with only 10**3 suffixes.
+
+ If 0 <= number < 1000 no suffix is added.
+ This is useful for quantities usually given in metric: stress, displacement.
+
+ Args:
+ number (float or int): the number we want converted to __metric_number
+ sig_figs (int): number of significant figures to use, right of decimal
+ sci (bool): True means use scientific formatting, False use metric
+ """
+ if sci:
+ format_str = "%.{}e".format(sig_figs)
+ my_str = format_str % number
+ else:
+ format_str = "%.{}f".format(sig_figs)
+ my_str = format_str % number
+ if number != 0:
+ # get the scientific exponent
+ exp = math.floor(math.log10(abs(number)))
+ metric_exp = exp - (exp % 3)
+ new_float = number/(10**metric_exp)
+ if metric_exp != 0:
+ format_str = "%.{}fe%i".format(sig_figs)
+ my_str = format_str % (new_float, metric_exp)
+ return my_str
+
+[docs] def load(self):
+ """Loads the results file with problem.fname prefix."""
+ self.__read_frd() # read nodal results
+ self.__read_dat() # read element integration pt results
+
+[docs] def nplot(self, field, fname='', display=True, levels=21, gradient=False,
+ gmult=1.0, max_val=None, min_val=None, title=''):
+ """Plots nodal results.
+
+ Args:
+ field (str): results item to plot, examples: 'ux', 'ey', 'Seqv'
+ fname (str): prefix of png file name, if writing an image
+ display (bool): True = interactively show the plot
+ levels (int): number of levels to use in the colorbar
+ gradient (bool): True = results plotted with gradient
+ False = results plotted with filled areas
+ gmult (int): geometric multiplier on displacement of nodes
+ displayed_node_loc = model_node_loc + gmult*node_displacement
+ max_val (float or None): max value in the colorbar
+
+ - None: max from selected data used
+ - float: use the passed float
+ min_val (float): min value in the colorbar
+
+ - None: min from selected data used
+ - float: use the passed float
+ title (str): third line in the plot title
+ """
+ # store the selected nodes and elements
+ sel = {}
+ sel['nodes'] = self.__problem.fea.view.nodes
+ sel['elements'] = self.__problem.fea.view.elements
+ sel['faces'] = self.__problem.fea.view.faces
+
+ # sort nodes low to high so index is correct
+ # we have index to id below so showing subsets works
+ sel['nodes'] = list(sel['nodes'])
+ sel['nodes'] = sorted(sel['nodes'], key=lambda k: k.id)
+
+ # store results at nodes
+ axials = []
+ radials = []
+ zvals = []
+ id_to_ind = {}
+ for node in sel['nodes']:
+ id_to_ind[node.id] = len(axials)
+ axi = node.y + gmult*self.__results[self.__time]['node'][node.id]['uy']
+ rad = node.x + gmult*self.__results[self.__time]['node'][node.id]['ux']
+ axials.append(axi)
+ radials.append(rad)
+ zvals.append(self.__results[self.__time]['node'][node.id][field])
+
+ # make a list of triangles, given by indices, looping anticlockwise
+ triangles = []
+ mylist = []
+ if len(sel['elements']) > 0:
+ mylist = sel['elements']
+ elif len(sel['faces']) > 0:
+ mylist = sel['faces']
+ for element in mylist:
+ tris = element.get_tris() # list of triangle nodes
+ for tri in tris:
+ for ind, nid in enumerate(tri):
+ tri[ind] = id_to_ind[nid] # convert id to index
+ triangles += tris
+
+ # check to see if selected nodes and elements are
+ # in the parent model's nodes and elements
+ fig = plt.figure()
+ ax_ = fig.add_subplot(111)
+
+ # need to set tick list here
+ vmin = min(zvals)
+ vmax = max(zvals)
+ stop_plot = False
+ if max_val != None and min_val == None:
+ if max_val < vmin:
+ stop_plot = True
+ print('Error:')
+ print(' Only max was passed but it is < the data min!')
+ print(' Pass a max_val that is > the data min of %f' % vmin)
+ else:
+ vmax = max_val
+ elif min_val != None and max_val == None:
+ if min_val > vmax:
+ stop_plot = True
+ print('Error:')
+ print(' Only min was passed but it is > the data max!')
+ print(' Pass a min_val that is < the data max of %f' % vmax)
+ else:
+ vmin = min_val
+ elif max_val != None and min_val != None:
+ if max_val < min_val:
+ stop_plot = True
+ print('Error:')
+ print(' Min and max passed, but min > max!')
+ print(' Pass a min_val that is < max_val')
+ else:
+ vmax = max_val
+ vmin = min_val
+ # exit if stop plot flag is on
+ if stop_plot:
+ return None
+
+ tick_list = [vmin]
+ if vmax != vmin:
+ # we have a range of values we're plotting
+ tick_list = linspace(vmin, vmax, levels+1)
+
+ # plot using a gradient(shaded) or levels
+ # code required for the colorbar, needs to go before plotting for colormap
+ cnorm = colors.Normalize(vmin=vmin, vmax=vmax)
+ cmap = colors.ListedColormap(['b', 'b']) # default to plot one val
+ if vmax != vmin:
+ # we have a range of values we're plotting
+ if gradient:
+ cmap = plt.get_cmap(CMAP)
+ else:
+ cmap = plt.get_cmap('jet', levels)
+ cmap.set_under('0.3', 0.8)
+ cmap.set_over('0.7', 0.8)
+ if gradient or len(tick_list) == 1:
+ # This one is shaded
+ plt.tripcolor(axials, radials, triangles, zvals, shading='gouraud',
+ cmap=cmap, norm=cnorm)
+ else:
+ # this one is not shaded
+ plt.tricontourf(axials, radials, triangles, zvals, levels=tick_list,
+ cmap=cmap, norm=cnorm, extend='both')
+
+ scalarmap = cmx.ScalarMappable(norm=cnorm, cmap=cmap)
+ scalarmap.set_array([])
+ cbar = plt.colorbar(scalarmap, orientation='vertical', ticks=tick_list)
+
+ scibool = False
+ if field[0] == 'e':
+ # strain plotting, use scientific numbering
+ scibool = True
+ met_max = self.__metric_num(max(zvals), sci=scibool)
+ met_min = self.__metric_num(min(zvals), sci=scibool)
+ label = 'Max: %s\nMin: %s' % (met_max, met_min)
+ tick_list = [self.__metric_num(tick, sci=scibool) for tick in tick_list]
+ cbar.ax.set_yticklabels(tick_list)
+ cbar.ax.set_xlabel(label, labelpad=10, x=0, ha='left')
+ cbar.ax.xaxis.set_label_position('top')
+
+ # set the horizontal and vertical axes
+ base_classes.plot_set_bounds(plt, axials, radials)
+
+ # set units
+ alist = self.__problem.fea.get_units(field, 'dist', 'time')
+ [f_unit, d_unit, t_unit] = alist
+
+ # set plot axes
+ plot_title = ('Node %s%s\nTime=%f%s' %
+ (field, f_unit, self.__time, t_unit))
+ if title != '':
+ plot_title += '\n%s' % title
+ plt.title(plot_title)
+ plt.xlabel('axial, y'+d_unit)
+ plt.ylabel('radial, x'+d_unit)
+ ax_.set_aspect('equal')
+ if gmult != 1:
+ ax_.xaxis.set_ticklabels([])
+ ax_.yaxis.set_ticklabels([])
+ base_classes.plot_finish(plt, fname, display)
+
+[docs] def eplot(self, field, fname='', display=True, levels=21,
+ gmult=1.0, mode='avg', max_val=None, min_val=None, title=''):
+ """Plots element results.
+
+ Args:
+ field (str): results item to plot. Only stresses supported.
+ Examples: 'Sx', 'Sxy', 'S1', 'Seqv' etc.
+ fname (str): prefix of png file name, if writing an image
+ display (bool): True = interactively show the plot
+ levels (int): number of levels to use in the colorbar
+ gmult (int): geometric multiplier on displacement of nodes
+ displayed_node_loc = model_node_loc + gmult*node_displacement
+ mode (str): the type of element result to plot
+
+ - 'avg': integration points averaged to avg element result
+ - 'max': max value of field in the integration points plotted
+ - 'min': min value of field in the integration points plotted
+ max_val (float or None): max value in the colorbar
+
+ - None: max from selected data used
+ - float: use the passed float
+ min_val (float): min value in the colorbar
+
+ - None: min from selected data used
+ - float: use the passed float
+ title (str): third line in the plot title
+ """
+ # store the selected nodes and elements
+ sel = {}
+ sel['nodes'] = self.__problem.fea.view.nodes
+ sel['elements'] = self.__problem.fea.view.elements
+ sel['faces'] = self.__problem.fea.view.faces
+
+ # sort nodes low to high so index is correct
+ # we have index to id below so showing subsets works
+ sel['nodes'] = list(sel['nodes'])
+ sel['nodes'] = sorted(sel['nodes'], key=lambda k: k.id)
+
+ # store results at nodes
+ axials = []
+ radials = []
+ zvals = []
+ id_to_ind = {}
+ for node in sel['nodes']:
+ id_to_ind[node.id] = len(axials)
+ axi = node.y + gmult*self.__results[self.__time]['node'][node.id]['uy']
+ rad = node.x + gmult*self.__results[self.__time]['node'][node.id]['ux']
+ axials.append(axi)
+ radials.append(rad)
+
+ # make a list of triangles, given by indices, looping anticlockwise
+ triangles = []
+ mylist = []
+ if len(sel['elements']) > 0:
+ mylist = sel['elements']
+ elif len(sel['faces']) > 0:
+ mylist = sel['faces']
+ for ele in mylist:
+ val = self.__results[self.__time]['element'][ele.id][mode][field]
+ tris = ele.get_tris() # list of triangle nodes defined by node id
+ for tri in tris:
+ zvals.append(val)
+ for ind, nid in enumerate(tri):
+ tri[ind] = id_to_ind[nid] # convert id to index
+ triangles += tris
+
+ # check to see if selected nodes and elements are
+ # in the parent model's nodes and elements
+
+ fig = plt.figure()
+ ax_ = fig.add_subplot(111)
+
+ # need to set tick list here
+ vmin = min(zvals)
+ vmax = max(zvals)
+ stop_plot = False
+ if max_val != None and min_val == None:
+ if max_val < vmin:
+ stop_plot = True
+ print('Error:')
+ print(' Only max was passed but it is < the data min!')
+ print(' Pass a max_val that is > the data min of %f' % vmin)
+ else:
+ vmax = max_val
+ elif min_val != None and max_val == None:
+ if min_val > vmax:
+ stop_plot = True
+ print('Error:')
+ print(' Only min was passed but it is > the data max!')
+ print(' Pass a min_val that is < the data max of %f' % vmax)
+ else:
+ vmin = min_val
+ elif max_val != None and min_val != None:
+ if max_val < min_val:
+ stop_plot = True
+ print('Error:')
+ print(' Min and max passed, but min > max!')
+ print(' Pass a min_val that is < max_val')
+ else:
+ vmax = max_val
+ vmin = min_val
+ # exit if stop plot flag is on
+ if stop_plot:
+ return None
+
+ tick_list = [vmin]
+ if vmax != vmin:
+ # we have a range of values we're plotting
+ tick_list = linspace(vmin, vmax, levels+1)
+
+ # code required for the colorbar, needs to go before plotting for cmap
+ cnorm = colors.Normalize(vmin=vmin, vmax=vmax)
+ cmap = colors.ListedColormap(['b', 'b']) # default to plot one val
+ if vmax != vmin:
+ # we have a range of values we're plotting
+ cmap = plt.get_cmap(CMAP, levels)
+ cmap.set_under('0.3', 0.8)
+ cmap.set_over('0.7', 0.8)
+
+ # plot using levels
+ plt.tripcolor(axials, radials, triangles, zvals,
+ shading='flat', cmap=cmap, norm=cnorm)
+
+ scalarmap = cmx.ScalarMappable(norm=cnorm, cmap=cmap)
+ scalarmap.set_array([])
+ cbar = plt.colorbar(scalarmap, orientation='vertical', ticks=tick_list)
+ scibool = False
+ if field[0] == 'e':
+ # strain plotting, use scientific numbering
+ scibool = True
+ met_max = self.__metric_num(max(zvals), sci=scibool)
+ met_min = self.__metric_num(min(zvals), sci=scibool)
+ label = 'Max: %s\nMin: %s' % (met_max, met_min)
+ tick_list = [self.__metric_num(tick, sci=scibool) for tick in tick_list]
+ cbar.ax.set_yticklabels(tick_list)
+ cbar.ax.set_xlabel(label, labelpad=10, x=0, ha='left')
+ cbar.ax.xaxis.set_label_position('top')
+
+ # set the horizontal and vertical axes
+ base_classes.plot_set_bounds(plt, axials, radials)
+
+ # set units
+ alist = self.__problem.fea.get_units(field, 'dist', 'time')
+ [f_unit, d_unit, t_unit] = alist
+
+ # set plot axes
+ plot_title = ('Element %s %s%s\nTime=%f%s' %
+ (mode, field, f_unit, self.__time, t_unit))
+ if title != '':
+ plot_title += '\n%s' % title
+ plt.title(plot_title)
+ plt.xlabel('axial, y'+d_unit)
+ plt.ylabel('radial, x'+d_unit)
+ ax_.set_aspect('equal')
+ if gmult != 1:
+ ax_.xaxis.set_ticklabels([])
+ ax_.yaxis.set_ticklabels([])
+ base_classes.plot_finish(plt, fname, display)
+
+[docs] def set_time(self, time):
+ """Sets the time point we're looking at in the results file.
+
+ Args:
+ time (float): time we are setting
+ """
+ if time in self.steps:
+ self.__time = time
+ print('Results file time set to: %f' % (self.__time))
+ else:
+ print('Time %f is not in the loaded times. Valid times are:')
+ print(self.steps)
+
+
+[docs] def plot_gradient(self, start_point, end_point, field, fname='', display=True, title='', max_val=None, min_val=None, curve_fitting=True, n_poly=3, n_subpoints=500, legend=True):
+ """Create diagram with data projected onto line on the undeformed geometry.
+
+ Args:
+ start_point [float, float]: starting point of line. [x, y]
+ end_point [float, float]: end point of line. Example: [x, y]
+ field (str): results item to plot, examples: 'ux', 'ey', 'Seqv'
+ fname (str): prefix of png file name, if writing an image
+ display (bool): True = interactively show the plot
+ title (str): third line in the plot title
+ max_val (float or None): max value in the y-axis
+ - None: max from selected data used
+ - float: use the passed float
+ min_val (float or None): min value in the y-axis
+ - None: min from selected data used
+ - float: use the passed float
+ curve_fitting (bool): True = a curve is fitted to the gradient
+ n_poly (int): numbers of polygons for fitting
+ n_subpoints (int): numbers of points the line is subdivided into
+ legend (bool): True = legend with fitted equation is shown
+ """
+
+ # store the selected nodes and elements
+ sel = {}
+ sel['nodes'] = self.__problem.fea.view.nodes
+
+ # sort nodes low to high so index is correct
+ # we have index to id below so showing subsets works
+ sel['nodes'] = list(sel['nodes'])
+ sel['nodes'] = sorted(sel['nodes'], key=lambda k: k.id)
+
+ # store results at nodes
+ node_position = np.zeros((len(sel['nodes']),2))
+ field_values = np.zeros(len(sel['nodes']))
+
+ for idx, node in enumerate(sel['nodes']):
+
+ node_position[idx] = [node.x, node.y]
+ field_values[idx] = self.__results[self.__time]['node'][node.id][field]
+
+
+ #create subpoints on line
+ subpoints = np.zeros((n_subpoints, 3)) #[x, y, line position]
+
+ subpoints[:,0] = np.linspace(start_point[0], end_point[0], n_subpoints)
+ subpoints[:,1] = np.linspace(start_point[1], end_point[1], n_subpoints)
+ subpoints[:,2] = np.arange(n_subpoints) / n_subpoints * np.sqrt(np.sum( (np.array(start_point) - np.array(end_point))**2))
+
+ #calculate weighted field value for every subpoint
+ wfield = np.zeros(n_subpoints)
+
+ for idx in range(n_subpoints):
+
+ #calculate inverse of distance from nodes to subpoints
+ dist = np.sqrt(np.sum((node_position-subpoints[idx,0:2])**2,axis=1))
+
+ #calculte weighted field value
+ #dist[dist < 1E-10] = 1E-10
+ #inv_dist = 1. / dist**3
+ #wfield[idx] = np.average(field_values, weights=inv_dist)
+
+ #use nearest value
+ wfield[idx] = field_values[min(range(len(dist)),key=dist.__getitem__)]
+
+
+ #plot diagram
+
+ fig = plt.figure(figsize=(10,6))
+ ax_ = fig.add_subplot(111)
+
+ plt.plot(subpoints[:,2], wfield, '-r', linewidth=2.5, label=field)
+
+ if curve_fitting==True:
+ #execute curve fitting if needed
+ poly = np.polyfit(subpoints[:,2], wfield, n_poly)
+
+ #string for equation of fitted function
+ funcstring = [str(np.round(poly[i]))+u'*x^'+str(np.arange(n_poly,0,-1)[i]) for i in range(n_poly)]
+ funcstring.append(str(np.round(poly[-1])))
+ funcstring = '+'.join(funcstring)
+
+ func = np.poly1d(poly)
+
+ plt.plot(subpoints[:,2], func(subpoints[:,2]), '--k', linewidth=1.5, label=funcstring)
+
+
+ # set units
+ alist = self.__problem.fea.get_units(field, 'dist', 'time')
+ [f_unit, d_unit, t_unit] = alist
+
+ # set plot axes
+ plot_title = ('Gradient %s%s\nTime=%f%s' %(field, f_unit, self.__time, t_unit))
+ if title != '':
+ plot_title += '\n%s' % title
+ plt.title(plot_title)
+ plt.xlabel('path position'+d_unit)
+ plt.ylabel(field + ' ' +f_unit)
+
+ #show legend if needed
+ if legend == True:
+ plt.legend()
+
+ #set limits on y-axis
+ if min_val!=None:
+ plt.gca().set_ylim(bottom=min_val)
+ if max_val!=None:
+ plt.gca().set_ylim(top=max_val)
+
+ plt.grid()
+ base_classes.plot_finish(plt, fname, display)
+
+[docs] def get_relative_gradient(self, start_point, end_point, field, n_poly=3, n_subpoints=500):
+ """Calculte relative stress gradient (gradient/start_value)
+
+ Args:
+ start_point [(float), (float)]: starting point of line. [x, y]
+ end_point [(float), (float)]: end point of line. Example: [x, y]
+ field (str): results item to plot, examples: 'ux', 'ey', 'Seqv'
+
+ Kargs:
+ n_poly (int): numbers of polygons for fitting, min=2
+ n_subpoints (int): numbers of points the line is subdivided into
+ """
+
+ # store the selected nodes and elements
+ sel = {}
+ sel['nodes'] = self.__problem.fea.view.nodes
+
+ # sort nodes low to high so index is correct
+ # we have index to id below so showing subsets works
+ sel['nodes'] = list(sel['nodes'])
+ sel['nodes'] = sorted(sel['nodes'], key=lambda k: k.id)
+
+ # store results at nodes
+ node_position = np.zeros((len(sel['nodes']),2))
+ field_values = np.zeros(len(sel['nodes']))
+
+ for idx, node in enumerate(sel['nodes']):
+
+ node_position[idx] = [node.x, node.y]
+ field_values[idx] = self.__results[self.__time]['node'][node.id][field]
+
+
+ #create subpoints on line
+ subpoints = np.zeros((n_subpoints, 3)) #[x, y, line position]
+
+ subpoints[:,0] = np.linspace(start_point[0], end_point[0], n_subpoints)
+ subpoints[:,1] = np.linspace(start_point[1], end_point[1], n_subpoints)
+ subpoints[:,2] = np.arange(n_subpoints) / n_subpoints * np.sqrt(np.sum( (np.array(start_point) - np.array(end_point))**2))
+
+ #calculate weighted field value for every subpoint
+ wfield = np.zeros(n_subpoints)
+
+ for idx in range(n_subpoints):
+
+ #calculate inverse of distance from nodes to subpoints
+ dist = np.sqrt(np.sum((node_position-subpoints[idx,0:2])**2,axis=1))
+
+ #use nearest value
+ wfield[idx] = field_values[min(range(len(dist)),key=dist.__getitem__)]
+
+
+ #curve fitting
+ poly = np.polyfit(subpoints[:,2], wfield, n_poly)
+
+ rel_grad = abs(poly[-2])/abs(poly[-1])
+
+ return rel_grad
+
+
+
+
+
+ @staticmethod
+ def __utot(vals):
+ """Returns the total displacement distance, given [dx,dy,dz].
+
+ Args:
+ vals (list): [dx, dy, dz] list of displacements in x, y, and z axes
+
+ Returns:
+ res (float): displacement
+ """
+ # computes sum of the squares
+ res = [a**2 for a in vals]
+ res = (sum(res))**0.5
+ return res
+
+ @staticmethod
+ def __seqv(vals):
+ """Returns the Von Mises stress, which will be stored as 'Seqv'.
+
+ Args:
+ vals (list): list of six stresses [s11,s22,s33,s12,s13,s23]
+
+ Returns:
+ res (float): Von Mises stress
+ """
+ [s11, s22, s33, s12, s13, s23] = vals
+ aval = s11 - s22
+ bval = s22 - s33
+ cval = s33 - s11
+ dval = s12**2 + s23**2 +s13**2
+ res = (0.5*(aval**2 + bval**2 + cval**2 +6*dval))**0.5
+ return res
+
+ @staticmethod
+ def __principals(vals):
+ """Returns principal stresses [S1,S2,S3].
+
+ Args:
+ vals (list): six stresses [s11,s22,s33,s12,s13,s23]
+
+ Returns:
+ res (list): principal stresses [S1,S2,S3] stresses are high-to-low
+ """
+ # calculates and returns principal stresses, S1, S2, S3
+ [s11, s22, s33, s12, s13, s23] = vals
+ aval = 1
+ bval = (s11 + s22 + s33)*-1.0
+ cval = (s11*s22 + s11*s33 + s22*s33 - s12**2 - s13**2 - s23**2)
+ dval = (s11*s22*s33 + 2*s12*s13*s23 - s11*(s23**2) - s22*(s13**2)
+ - s33*(s12**2))*-1.0
+ res = list(roots([aval, bval, cval, dval]))
+ res = sorted(res, reverse=True)
+ return res
+
+ def __get_data_dict(self, time, type_str):
+ """Returns the data dict at the correct time for element or node.
+
+ Args:
+ time (float): None or the time we want, if None use current time
+ type_str: 'element' or 'node'
+
+ Returns:
+ res (dict or None): dictionary with field values in it
+ None if the time was invalid
+ """
+ res = self.__results[self.__time][type_str]
+ if time != None:
+ if time not in self.steps:
+ print('Error: passed time is not in steps!')
+ print(' Pass a time in the steps:')
+ print(self.steps)
+ return None
+ else:
+ res = self.__results[time][type_str]
+ return res
+
+[docs] def get_nmax(self, field, time=None):
+ """Returns the max value of node results field in selected nodes.
+
+ Reports results for the current time.
+
+ Args:
+ field (str): results field, for example 'ux', 'ey', 'S1', 'Seqv'
+ time (None or float): the time to query
+
+ - None: uses the current time
+ - float: uses the passed float time
+ Returns:
+ res (float): max value
+ """
+ nodes = self.__problem.fea.view.nodes
+ node_ids = [node.id for node in nodes]
+ data_dict = self.__get_data_dict(time, 'node')
+ if data_dict == None:
+ return None
+ ndicts = [data_dict[nid] for nid in node_ids]
+ res = [ndict[field] for ndict in ndicts]
+ res = max(res)
+ return res
+
+[docs] def get_nmin(self, field, time=None):
+ """Returns the min value of node results field in selected nodes.
+
+ Reports results for the current time.
+
+ Args:
+ field (str): results field, for example 'ux', 'ey', 'S1', 'Seqv'
+ time (None or float): the time to query
+
+ - None: uses the current time
+ - float: uses the passed float time
+ Returns:
+ res (float): min value
+ """
+ nodes = self.__problem.fea.view.nodes
+ node_ids = [node.id for node in nodes]
+ data_dict = self.__get_data_dict(time, 'node')
+ if data_dict == None:
+ return None
+ ndicts = [data_dict[nid] for nid in node_ids]
+ res = [ndict[field] for ndict in ndicts]
+ res = min(res)
+ return res
+
+[docs] def get_nval(self, node, field, time=None):
+ """Returns the field result value under node.
+
+ Result will be returned whether or not passed node is selected.
+
+ Args:
+ node (str or Node): node we are asking about
+ field (str): the results item we want: 'ux', 'Sy', 'Seqv', 'fx'
+ time (None or float): the time to query
+
+ - None: uses the current time
+ - float: uses the passed float time
+ Returns:
+ res (float or None): float value if field exists, None otherwise
+ """
+ items = self.__problem.fea.get_item(node)
+ if len(items) == 1:
+ if isinstance(items[0], mesh.Node):
+ nnum = items[0].id
+ data_dict = self.__get_data_dict(time, 'node')
+ if data_dict == None:
+ return None
+ ndict = data_dict[nnum]
+ if field in ndict:
+ res = ndict[field]
+ return res
+ else:
+ print('Passed field is not in the results!')
+ return None
+ else:
+ print('You did not pass in a node!')
+ print('A single node or string node name must be passed in!')
+ return None
+ else:
+ print('A single node or string node name must be passed in!')
+ return None
+
+[docs] def get_fsum(self, item):
+ """Returns the force sum on nodes under a given point or line.
+
+ Reports results for the current time.
+
+ Args:
+ item (Point or SignLine): item that has reaction forces on its nodes
+
+ Returns:
+ list: [fx, fy, fz] reaction forces in each axis, force units
+ """
+ (fxx, fyy, fzz) = ([], [], [])
+ nodes = item.nodes
+ nodes = [n.id for n in nodes]
+ for node in nodes:
+ f_x = self.__results[self.__time]['node'][node]['fx']
+ f_y = self.__results[self.__time]['node'][node]['fy']
+ f_z = self.__results[self.__time]['node'][node]['fz']
+ if f_x != 0 or f_y != 0 or f_z != 0:
+ fxx.append(f_x)
+ fyy.append(f_y)
+ fzz.append(f_z)
+ fxx = sum(fxx)
+ fyy = sum(fyy)
+ fzz = sum(fzz)
+ return [fxx, fyy, fzz]
+
+[docs] def get_emax(self, field, time=None, mode='avg'):
+ """Returns the max results field value of selected elements at curent time.
+
+ Args:
+ field (str): results field, stresses supported 'S1', 'Sx', etc.
+ time (None or float): the time to query
+
+ - None: uses the current time
+ - float: uses the passed float time
+ mode (str): type of element result to give back
+
+ - 'max': for each element only use the max value of field over
+ all of its integration points
+ - 'min': for each element only use the min value of field over
+ all of its integration points
+ - 'avg': for each element only use an average of all integration
+ points in the eleemnt. Principal streses and Seqv are
+ calculated after averaging 6 stress components.
+ Returns:
+ res (float): max value
+ """
+ res = []
+ elements = self.__problem.fea.view.elements
+ data_dict = self.__get_data_dict(time, 'element')
+ if data_dict == None:
+ return None
+ for element in elements:
+ enum = element.id
+ edict = data_dict[enum][mode]
+ res.append(edict[field])
+ res = max(res)
+ return res
+
+[docs] def get_emin(self, field, time=None, mode='avg'):
+ """Returns the min results field value of selected elements at curent time.
+
+ Args:
+ field (str): results field, stresses supported 'S1', 'Sx', etc.
+ time (None or float): the time to query
+
+ - None: uses the current time
+ - float: uses the passed float time
+ mode (str): type of element result to give back
+
+ - 'max': for each element only use the max value of field over
+ all of its integration points
+ - 'min': for each element only use the min value of field over
+ all of its integration points
+ - 'avg': for each element only use an average of all integration
+ points in the eleemnt. Principal streses and Seqv are
+ calculated after averaging 6 stress components.
+ Returns:
+ res (float): min value
+ """
+ res = []
+ elements = self.__problem.fea.view.elements
+ data_dict = self.__get_data_dict(time, 'element')
+ if data_dict == None:
+ return None
+ for element in elements:
+ enum = element.id
+ edict = data_dict[enum][mode]
+ res.append(edict[field])
+ res = min(res)
+ return res
+
+[docs] def get_eval(self, element, field, time=None, mode='avg'):
+ """Returns the field result value under element.
+
+ Result will be returned whether or not passed element is selected.
+
+ Args:
+ element (str or Element): element we are asking about
+ field (str): the results item we want: 'Sy', 'Seqv'
+ mode (str): the type of element result to get
+
+ - 'avg': integration points averaged to avg element result
+ - 'max': max value of field in the integration points plotted
+ - 'min': min value of field in the integration points plotted
+ Returns:
+ res (float or None): float value if field exists, None otherwise
+ """
+ items = self.__problem.fea.get_item(element)
+ if len(items) == 1:
+ if isinstance(items[0], mesh.Element):
+ enum = items[0].id
+ data_dict = self.__get_data_dict(time, 'element')
+ if data_dict == None:
+ return None
+ edict = data_dict[enum][mode]
+ if field in edict:
+ res = edict[field]
+ return res
+ else:
+ print('Passed field is not in the results!')
+ return None
+ else:
+ print('You did not pass in a element!')
+ print('A single element or string element name must be given!')
+ return None
+ else:
+ print('A single element or string element name must be given!')
+ return None
+
+ @staticmethod
+ def __get_vals(fstr, line):
+ """Returns a list of typed items based on an input format string.
+
+ Args:
+ fst (str): C format string, commas separate fields
+ line (str): line string to parse
+
+ Returns:
+ res (list): list of typed items extracted from the line
+ """
+ res = []
+ fstr = fstr.split(',')
+ thestr = str(line)
+ for item in fstr:
+ if item[0] == "'":
+ # strip off the char quaotes
+ item = item[1:-1]
+ # this is a string entry, grab the val out of the line
+ ind = len(item)
+ fwd = thestr[:ind]
+ thestr = thestr[ind:]
+ res.append(fwd)
+ else:
+ # format is: 1X, A66, 5E12.5, I12
+ # 1X is number of spaces
+ (mult, ctype) = (1, None)
+ m_pat = re.compile(r'^\d+') # find multiplier
+ c_pat = re.compile(r'[XIEA]') # find character
+ if m_pat.findall(item) != []:
+ mult = int(m_pat.findall(item)[0])
+ ctype = c_pat.findall(item)[0]
+ if ctype == 'X':
+ # we are dealing with spaces, just reduce the line size
+ thestr = thestr[mult:]
+ elif ctype == 'A':
+ # character string only, add it to results
+ fwd = thestr[:mult].strip()
+ thestr = thestr[mult:]
+ res.append(fwd)
+ else:
+ # IE, split line into m pieces
+ w_pat = re.compile(r'[IE](\d+)') # find the num after char
+ width = int(w_pat.findall(item)[0])
+ while mult > 0:
+ # only add items if we have enough line to look at
+ if width <= len(thestr):
+ substr = thestr[:width]
+ thestr = thestr[width:]
+ substr = substr.strip() # remove space padding
+ if ctype == 'I':
+ substr = int(substr)
+ elif ctype == 'E':
+ substr = float(substr)
+ res.append(substr)
+ mult -= 1
+ return res
+
+ @staticmethod
+ def __get_first_dataline(infile):
+ """
+ Reads infile until a line with data is found, then returns it
+ A line that starts with ' -1' has data
+ """
+ while True:
+ line = infile.readline()
+ if line[:3] == ' -1':
+ return line
+
+ def __store_time(self, time):
+ """Stores the passed time in the results steps"""
+ if time not in self.__steps:
+ self.__steps.append(time)
+ if time not in self.__results:
+ new_dict = {'node': collections.defaultdict(dict),
+ 'element': collections.defaultdict(dict)}
+ self.__results[time] = new_dict
+
+ def __modearr_estrsresults(self, infile, line):
+ """Returns an array of line, mode, rfstr, time"""
+ words = line.strip().split()
+ # add time if not present
+ time = float(words[-1])
+ self.__store_time(time)
+ # set mode
+ rfstr = "I10,2X,I2,6E14.2"
+ mode = 'stress'
+ infile.readline()
+ line = infile.readline()
+ return [line, mode, rfstr, time]
+
+ def __modearr_nresults(self, infile):
+ """Returns an array of line, mode, rfstr, time"""
+ line = infile.readline()
+ fstr = "1X,' 100','C',6A1,E12.5,I12,20A1,I2,I5,10A1,I2"
+ tmp = self.__get_vals(fstr, line)
+ #[key, code, setname, value, numnod, text, ictype, numstp, analys, format_]
+ time, format_ = tmp[3], tmp[9]
+
+ # set results format to short, long or binary
+ # only short and long are parsed so far
+ if format_ == 0:
+ rfstr = "1X,I2,I5,6E12.5"
+ elif format_ == 1:
+ rfstr = "1X,I2,I10,6E12.5"
+ elif format_ == 2:
+ # binary
+ pass
+
+ # set the time
+ self.__store_time(time)
+
+ # get the name to determine if stress or displ
+ line = infile.readline()
+ fstr = "1X,I2,2X,8A1,2I5"
+ # [key, name, ncomps, irtype]
+ name = self.__get_vals(fstr, line)[1]
+
+ mode_by_name = {'DISP': 'displ',
+ 'STRESS': 'stress',
+ 'TOSTRAIN': 'strain',
+ 'FORC': 'force'}
+ mode = mode_by_name[name]
+ print('Reading '+mode+' storing: '+
+ ','.join(base_classes.RESFIELDS[mode]))
+
+ line = self.__get_first_dataline(infile)
+ return [line, mode, rfstr, time]
+
+ def _save_node_displ(self, line, rfstr, time, mode='displ'):
+ """Saves node displacement"""
+ node, ux_, uy_, uz_ = self.__get_vals(rfstr, line)[1:]
+ labs = base_classes.RESFIELDS[mode]
+ vals = [ux_, uy_, uz_]
+ utot = self.__utot(vals)
+ vals.append(utot)
+ adict = self.__results[time]['node'][node]
+ for (label, val) in zip(labs, vals):
+ adict[label] = val
+
+ def _save_node_stress(self, line, rfstr, time, mode='stress'):
+ """Saves node stress"""
+ tmp = self.__get_vals(rfstr, line)
+ # [key, node, sx, sy, sz, sxy, syz, szx]
+ node, sxx, syy, szz, sxy, syz, szx = tmp[1:]
+ labs = base_classes.RESFIELDS[mode]
+ vals = [sxx, syy, szz, sxy, syz, szx]
+ seqv = self.__seqv(vals)
+ s_1, s_2, s_3 = self.__principals(vals)
+ vals.append(seqv)
+ vals += [s_1, s_2, s_3]
+ adict = self.__results[time]['node'][node]
+ for (label, val) in zip(labs, vals):
+ adict[label] = val
+
+ def _save_node_strain(self, line, rfstr, time, mode='strain'):
+ """Saves node strain"""
+ tmp = self.__get_vals(rfstr, line)
+ # [key, node, ex, ey, ez, exy, eyz, ezx]
+ node, exx, eyy, ezz, exy, eyz, ezx = tmp[1:]
+ labs = base_classes.RESFIELDS[mode]
+ vals = [exx, eyy, ezz, exy, eyz, ezx]
+ eeqv = self.__seqv(vals)
+ e_1, e_2, e_3 = self.__principals(vals)
+ vals.append(eeqv)
+ vals += [e_1, e_2, e_3]
+ adict = self.__results[time]['node'][node]
+ for (label, val) in zip(labs, vals):
+ adict[label] = val
+
+ def _save_node_force(self, line, rfstr, time, mode='force'):
+ """Saves node force"""
+ # [key, node, fx, fy, fz]
+ node, f_x, f_y, f_z = self.__get_vals(rfstr, line)[1:]
+ labs = base_classes.RESFIELDS[mode]
+ vals = [f_x, f_y, f_z]
+ adict = self.__results[time]['node'][node]
+ for (label, val) in zip(labs, vals):
+ adict[label] = val
+
+ def _save_ele_stress(self, line, rfstr, time,
+ mode='stress'):
+ """Saves element integration point stresses"""
+ labels = ['Sx', 'Sy', 'Sz', 'Sxy', 'Sxz', 'Syz']
+ vals = self.__get_vals(rfstr, line)
+ # element_number, integration_pt_number
+ enum, ipnum = vals[0], vals[1]
+ stress_vals = vals[2:]
+
+ adict = {}
+ for (label, val) in zip(labels, stress_vals):
+ adict[label] = val
+ if enum not in self.__results[time]['element']:
+ start_val = {'ipoints': {},
+ 'avg': {},
+ 'min': {},
+ 'max': {}}
+ self.__results[time]['element'][enum] = start_val
+ # each line is an integration point result
+ self.__results[time]['element'][enum]['ipoints'][ipnum] = adict
+
+[docs] def check_ccx_version(self, timeout=1):
+ """Raises an exception of the calculix ccx version is too old"""
+ runstr = "%s -version" % (environment.CCX)
+ try:
+ output_str = subprocess.check_output(runstr,
+ timeout=timeout,
+ shell=True)
+ except (subprocess.CalledProcessError, subprocess.TimeoutExpired) as ex:
+ output_str = ex.output
+ output_str = str(output_str, 'utf-8')
+ matches = re.findall(r'\d+\.\d+', output_str)
+ version_number = matches[-1]
+ print('Using Calculix ccx version=%s '
+ '(trailing characters like the p in 2.8p are omitted)'
+ % version_number)
+ major_version, minor_version = [int(f) for f
+ in version_number.split('.')]
+ if major_version <= 2 and minor_version <= 8:
+ raise Exception('Your version of calculix ccx is too old! '
+ 'Please update it to version >=2.8 with '
+ 'the command:\npycalculix-add-feaprograms')
+ version = float(output_str.strip().split()[-1])
+ # extract version with regex
+
+ def __read_frd(self):
+ """
+ Reads a ccx frd results file which contains nodal results.
+ The file format is desribed in the cgx docs
+ """
+ fname = self.__problem.fname+'.frd'
+ if not os.path.isfile(fname):
+ print("Error: %s file not found" % fname)
+ return
+ # frd reading uses formatting from ccx 2.8 or higher so
+ # throw an exception if our version is too old
+ self.check_ccx_version(timeout=1)
+ infile = open(fname, 'r')
+ print('Loading nodal results from file: '+fname)
+ mode = None
+ time = 0.0
+ rfstr = ''
+ while True:
+ line = infile.readline()
+ if not line:
+ break
+
+ # set the results mode
+ if '1PSTEP' in line:
+ # we are in a results block
+ arr = self.__modearr_nresults(infile)
+ line, mode, rfstr, time = arr
+
+ # set mode to none if we hit the end of a resuls block
+ if line[:3] == ' -3':
+ mode = None
+ if not mode:
+ continue
+
+ node_data_saver = getattr(self, '_save_node_' + mode)
+ node_data_saver(line, rfstr, time)
+
+ infile.close()
+ print('The following times have been read:')
+ print(self.__steps)
+ print('Nodal results from file: %s have been read.' % fname)
+ self.set_time(self.__steps[0])
+
+ def __read_dat(self):
+ """
+ Reads ccx dat results file.
+ It has element integration point results.
+ """
+ fname = self.__problem.fname+'.dat'
+ if not os.path.isfile(fname):
+ print('Error: %s file not found' % fname)
+ return
+ infile = open(fname, 'r')
+ print('Loading element results from file: '+fname)
+ mode = None
+ rfstr = ''
+ time = 0.0
+ while True:
+ line = infile.readline()
+ if not line:
+ break
+
+ # check for stress, we skip down to the line data when
+ # we call __modearr_estrsresults
+ if 'stress' in line:
+ arr = self.__modearr_estrsresults(infile, line)
+ line, mode, rfstr, time = arr
+
+ # reset the read type if we hit a blank line
+ if line.strip() == '':
+ mode = None
+ if not mode:
+ continue
+
+ # store stress results
+ self._save_ele_stress(line, rfstr, time)
+
+ infile.close()
+
+ # loop over all element results, calculating avg element result
+ # by averaging integration point vals
+ for time in self.__steps:
+ for edict in self.__results[time]['element'].values():
+ ipoints = edict['ipoints'].values()
+ stress_types = ['Sx', 'Sy', 'Sz', 'Sxy', 'Sxz', 'Syz']
+ strslist_by_strstype = collections.defaultdict(list)
+ # set stress values in max, min, avg locations
+ # of non-summary stress components
+ for stress_type in stress_types:
+ stress_vals = [ipt[stress_type] for ipt in ipoints]
+ stress_avg = sum(stress_vals)/len(stress_vals)
+ stress_max = max(stress_vals)
+ stress_min = min(stress_vals)
+ edict['avg'][stress_type] = stress_avg
+ edict['max'][stress_type] = stress_max
+ edict['min'][stress_type] = stress_min
+ strslist_by_strstype[stress_type] = stress_vals
+ # for each element, calc Seqv, S1, S2, S3
+ # at each integration point
+ for ipt in ipoints:
+ stress_vals = [ipt[stress_type] for stress_type
+ in stress_types]
+ seqv = self.__seqv(stress_vals)
+ [s_1, s_2, s_3] = self.__principals(stress_vals)
+ strslist_by_strstype['Seqv'].append(seqv)
+ strslist_by_strstype['S1'].append(s_1)
+ strslist_by_strstype['S2'].append(s_2)
+ strslist_by_strstype['S3'].append(s_3)
+ # now at the element level, store the avg, max, and min
+ # of the Seqv and principal stresses we found at
+ # integration points
+ for stress_type in ['Seqv', 'S1', 'S2', 'S3']:
+ stress_vals = strslist_by_strstype[stress_type]
+ stress_avg = sum(stress_vals)/len(stress_vals)
+ stress_max = max(stress_vals)
+ stress_min = min(stress_vals)
+ edict['avg'][stress_type] = stress_avg
+ edict['max'][stress_type] = stress_max
+ edict['min'][stress_type] = stress_min
+
+
+ print('The following times have been read:')
+ print(self.__steps)
+ print('Element results from file: %s have been read.' % fname)
+
+"""This module stores the Selector object
+
+Selector is used in FeaModel to store the model's selected set.
+"""
+
+from . import partmodule
+from . import geometry
+from . import mesh
+from . import components
+
+[docs]class Selector(object):
+ """Makes a selector which stores the selected items in the FeaModel.
+
+ Args:
+ feamodel (FeaModel): the parent FeaModel object
+
+ Attributes:
+ __fea (FeaModel): the parent FeaModel
+ __all (bool): If everything is selected all=True
+ parts (set): currently selected parts
+ areas (set): currently selected areas
+ lines (set): currently selected signed lines or arcs
+ points (set): currently selected points
+ elements (set): currently selected elements
+ faces (set): currently selected faces
+ nodes (set): currently selected nodes
+ __parts (set): internal storage set of parts
+ __areas (set): internal storage set of areas
+ __slines (set): internal storage set of signed lines and arcs
+ __points (set): internal storage set of points
+ __elements (set): internal storage set of elements
+ __faces (set): internal storage set of faces
+ __nodes (set): internal storage set of nodes
+ """
+ def __init__(self, feamodel):
+ self.__fea = feamodel
+ self.__all = True
+ self.__parts = set()
+ self.__areas = set()
+ self.__slines = set()
+ self.__points = set()
+ self.__elements = set()
+ self.__faces = set()
+ self.__nodes = set()
+
+ @property
+ def parts(self):
+ """Returns selected parts."""
+ if self.__all:
+ return self.__fea.parts
+ else:
+ return self.__parts
+
+ @property
+ def areas(self):
+ """Returns selected areas."""
+ if self.__all:
+ return self.__fea.areas
+ else:
+ return self.__areas
+
+ @property
+ def lines(self):
+ """Returns selected signed lines or arcs."""
+ if self.__all:
+ return self.__fea.signlines
+ else:
+ return self.__slines
+
+ @property
+ def points(self):
+ """Returns selected points."""
+ if self.__all:
+ points = [pnt for pnt in self.__fea.points if pnt.arc_center == False]
+ return points
+ else:
+ return self.__points
+
+ @property
+ def elements(self):
+ """Returns selected elements."""
+ if self.__all:
+ return self.__fea.elements
+ else:
+ return self.__elements
+
+ @property
+ def faces(self):
+ """Returns selected faces."""
+ if self.__all:
+ return self.__fea.faces
+ else:
+ return self.__faces
+
+ @property
+ def nodes(self):
+ """Returns selected nodes."""
+ if self.__all:
+ return self.__fea.nodes
+ else:
+ return self.__nodes
+
+ @staticmethod
+ def __add_items(inclusive, items, item, parent_list, this_set):
+ """Adds an item to items if appropriate.
+
+ Helper function for select_above_all
+ """
+ if inclusive == False:
+ items.add(item)
+ else:
+ parent_set = set(parent_list)
+ if parent_set.issubset(this_set):
+ items.add(item)
+
+ def __add_select(self, items):
+ """Adds items to the selection set."""
+ for item in items:
+ if isinstance(item, partmodule.Part):
+ self.__parts.add(item)
+ elif isinstance(item, geometry.Area):
+ self.__areas.add(item)
+ elif (isinstance(item, geometry.SignLine)
+ or isinstance(item, geometry.SignArc)):
+ self.__slines.add(item)
+ elif isinstance(item, geometry.Point):
+ self.__points.add(item)
+ elif isinstance(item, mesh.Element):
+ self.__elements.add(item)
+ elif isinstance(item, mesh.Face):
+ self.__faces.add(item)
+ elif isinstance(item, mesh.Node):
+ self.__nodes.add(item)
+
+[docs] def allsel_under(self, sel_type, full=True, byfaces=False):
+ """Selects all progeny under current sel_type items.
+
+ Args:
+ sel_type (str): the item type to select under
+
+ - 'parts'
+ - 'areas'
+ - 'lines'
+ - 'points' (needs full=True)
+ - 'elements'
+ - 'faces'
+
+ full (bool): if False selects under one branch
+
+ - solid modeling branch: area->line->keypoints->points
+ - mesh branch: elements->faces->nodes
+
+ Cross branch connections:
+ - faces are under lines
+ - elements are under areas
+ - nodes are under lines and points
+
+ If True, selects all items connected to sub items:
+ - Parts:
+
+ - Areas
+ - Elements, element faces + element nodes
+ - Lines, line points
+
+ - Areas:
+
+ - Elements, element faces + element nodes
+ - Lines, line points
+
+ -Lines:
+
+ - Faces, nodes
+ - byfaces=True: elements above faces
+ - byfaces=False: elements above line nodes selected
+
+ -Points:
+
+ - Nodes
+ - byfaces=True: no elements selected
+ - byfaces=False: elements above point nodes selected
+
+ byfaces (bool): If true, the child elements will only be selected
+ if a parent face was selected. If false just an element node
+ needs to be selected to select the element.
+
+ -This is only applicable when sel_type == 'lines'
+ """
+ sets = ['parts', 'areas', 'lines', 'points', 'elements', 'faces']
+ if isinstance(sel_type, str):
+ if sel_type in sets:
+ this_set = getattr(self, sel_type)
+ for item in this_set:
+ if isinstance(item, partmodule.Part):
+ # areas
+ self.select_below_all('parts')
+ # slines
+ self.select_below_all('areas')
+ # points
+ self.select_below_all('lines')
+ # select related elements if full is on
+ if full:
+ elements = item.elements
+ self.select(elements, also=True)
+ # sel faces
+ self.select_below_all('elements')
+ # sel nodes
+ self.select_below_all('faces')
+
+ elif isinstance(item, geometry.Area):
+ # slines
+ self.select_below_all('areas')
+ # points
+ self.select_below_all('lines')
+ # select related elements if full is on
+ if full:
+ elements = item.elements
+ self.select(elements, also=True)
+ # sel faces
+ self.select_below_all('elements')
+ # sel nodes
+ self.select_below_all('faces')
+ elif (isinstance(item, geometry.SignLine)
+ or isinstance(item, geometry.SignArc)):
+ # sel points
+ self.select_below_all('lines')
+ # select related faces, nodes, and elements
+ if full:
+ # select faces
+ faces = item.faces
+ self.select(faces, also=True)
+ # sel nodes
+ self.select_below_all('faces')
+ # byfaces logic
+ if byfaces:
+ # sel elements above faces
+ self.select_above_all('faces', inclusive=False)
+ else:
+ # select elements connected to nodes
+ elements = set()
+ for node in self.nodes:
+ elements.update(node.elements)
+ self.select(elements, also=True)
+ elif isinstance(item, geometry.Point):
+ if full:
+ # sel nodes and elements
+ self.select(item.nodes, also=True)
+ if byfaces == False:
+ for node in item.nodes:
+ self.select(node.elements, also=True)
+ else:
+ print('The passed sel_type must be a valid type.')
+ else:
+ print('The passed sel_type must be a string!')
+
+[docs] def select_below(self):
+ """Selects children below currently selected items.
+
+ Only one set must be selected.
+ """
+ sets = ['parts', 'areas', 'lines', 'points', 'elements', 'faces',
+ 'nodes']
+ sel_sets = []
+ for set_name in sets:
+ this_set = getattr(self, set_name)
+ if len(this_set) > 0:
+ sel_sets.append(set_name)
+ # check if number of visible sets is one
+ if len(sel_sets) == 0 or len(sel_sets) > 1:
+ print('Only one set must be selected to use this function.')
+ elif len(sel_sets) == 1:
+ sel_type = sel_sets[0]
+ if sel_type == 'nodes':
+ print("One can't select anything below nodes!'")
+ else:
+ self.select_below_all(sel_type)
+
+[docs] def select_below_all(self, sel_type):
+ """Selects children below all currently selected items of sel_type.
+
+ Args:
+ sel_type (str): type to select items below
+
+ - 'parts': selects areas
+ - 'areas': selects lines (type SignLine and SignArc)
+ - 'lines': selects points
+ - 'elements': selects faces
+ - 'faces': selects nodes
+ """
+ items = set()
+ if sel_type == 'parts':
+ for part in self.parts:
+ for area in part.areas:
+ items.add(area)
+ elif sel_type == 'areas':
+ for area in self.areas:
+ items.update(area.signlines)
+ elif sel_type == 'lines':
+ for sline in self.lines:
+ items.update(sline.line.points)
+ elif sel_type == 'elements':
+ for element in self.elements:
+ items.update(element.faces)
+ elif sel_type == 'faces':
+ for face in self.faces:
+ items.update(face.nodes)
+ self.__add_select(items)
+
+[docs] def select_above(self, inclusive=False):
+ """Selects parents above currently selected items.
+
+ Only one set must be selected.
+ """
+ sets = ['parts', 'areas', 'lines', 'points', 'elements', 'faces',
+ 'nodes']
+ sel_sets = []
+ for set_name in sets:
+ this_set = getattr(self, set_name)
+ if len(this_set) > 0:
+ sel_sets.append(set_name)
+ # check if number of visible sets is one
+ if len(sel_sets) == 0 or len(sel_sets) > 1:
+ print('Only one set must be selected to use this function.')
+ elif len(sel_sets) == 1:
+ sel_type = sel_sets[0]
+ if sel_type == 'parts':
+ print("One can't select anything above parts!'")
+ else:
+ self.select_above_all(sel_type, inclusive)
+
+[docs] def select_above_all(self, sel_type, inclusive=False):
+ """Selects parents above all currently selected items of sel_type.
+
+ Args:
+ sel_type (str): type to select items above
+
+ - 'areas': selects parts
+ - 'lines': selects areas
+ - 'points': selects signlines or signarcs
+ - 'faces': selects elements
+ - 'nodes': selects faces
+
+ inclusive (bool): if True, all child entities must be selected for
+ the parent to be selected. If false all parents connected to
+ the child are selected. Exampes:
+
+ - 'areas' + True: all part areas must be selected to select part
+ - 'areas' + False: parts with selected areas in them are selected
+ - 'lines' + True: all lines must be selected to select area
+ - 'lines' + False: areas with selected lines in them are selected
+ """
+ items = set()
+ if sel_type == 'areas':
+ this_set = self.areas
+ for area in this_set:
+ item = area.parent
+ parent_list = area.parent.areas
+ self.__add_items(inclusive, items, item, parent_list, this_set)
+ elif sel_type == 'lines':
+ this_set = self.lines
+ for sline in this_set:
+ item = sline.parent
+ parent_list = sline.parent.signlines
+ self.__add_items(inclusive, items, item, parent_list, this_set)
+ elif sel_type == 'points':
+ this_set = set(self.points)
+ for point in this_set:
+ for parent_line in point.lines:
+ for sline in parent_line.signlines:
+ item = sline
+ parent_list = sline.points
+ self.__add_items(inclusive, items, item, parent_list,
+ this_set)
+ elif sel_type == 'faces':
+ this_set = set(self.faces)
+ for face in this_set:
+ item = face.element
+ parent_list = item.faces
+ self.__add_items(inclusive, items, item, parent_list, this_set)
+ elif sel_type == 'nodes':
+ this_set = set(self.nodes)
+ for node in this_set:
+ for face in node.faces:
+ item = face
+ parent_list = face.nodes
+ self.__add_items(inclusive, items, item, parent_list,
+ this_set)
+ self.__add_select(items)
+
+
+
+[docs] def select(self, items='', also=False):
+ """Selects an item or a list of items.
+
+ Agrs:
+ items (str or item or list): item(s) to select
+ also (bool): whether to add the items to the current selection set
+
+ - True: add the items to the current selection set
+ - False: only select the passed items
+ """
+ # empty the selection set if also == False
+ if items == '':
+ print('You must pass in an item or item(s)!')
+ else:
+ if also == False:
+ self.select_none()
+ # Turn on 'all' mode if items == None, otherwise, get the items
+ if items == 'all':
+ self.__all = True
+ else:
+ items = self.__fea.get_items(items)
+ for (ind, item) in enumerate(items):
+ # if a component was passed, get its child items instead
+ if isinstance(item, components.Component):
+ real_items = item.get_children()
+ items[ind] = real_items
+ self.__add_select(items)
+
+[docs] def select_all(self, sel_type='all', also=False):
+ """Selects all items of sel_type from the full feamodel.
+
+ All items currently selected items are removed from the selection set
+ before the all set is selected.
+
+ Args:
+ sel_type (str): the items to select
+
+ - 'all'
+ - 'parts'
+ - 'areas'
+ - 'lines'
+ - 'points'
+ - 'elements'
+ - 'faces'
+ - 'nodes'
+ also (bool): if False, empties the selection set before selecting
+ all items, if True adds new items to the current set
+ """
+ if also == False:
+ self.__all = False
+ self.select_none()
+ # now select all items of type
+ if sel_type == 'all':
+ self.select()
+ elif sel_type == 'parts':
+ self.__add_select(self.__fea.parts)
+ elif sel_type == 'areas':
+ self.__add_select(self.__fea.areas)
+ elif sel_type == 'lines':
+ self.__add_select(self.__fea.signlines)
+ elif sel_type == 'points':
+ points = [pnt for pnt in self.__fea.points if pnt.arc_center == False]
+ self.__add_select(points)
+ elif sel_type == 'elements':
+ self.__add_select(self.__fea.elements)
+ elif sel_type == 'faces':
+ self.__add_select(self.__fea.faces)
+ elif sel_type == 'nodes':
+ self.__add_select(self.__fea.nodes)
+
+[docs] def deselect(self, items):
+ """Removes the passed item or items from the selection set.
+ """
+ self.__all = False
+ for item in items:
+ if isinstance(item, partmodule.Part):
+ self.__parts.discard(item)
+ elif isinstance(item, geometry.Area):
+ self.__areas.discard(item)
+ elif (isinstance(item, geometry.SignLine)
+ or isinstance(item, geometry.SignArc)):
+ self.__slines.discard(item)
+ elif isinstance(item, geometry.Point):
+ self.__points.discard(item)
+ elif isinstance(item, mesh.Element):
+ self.__elements.discard(item)
+ elif isinstance(item, mesh.Face):
+ self.__faces.discard(item)
+ elif isinstance(item, mesh.Node):
+ self.__nodes.discard(item)
+
+[docs] def deselect_all(self, sel_type):
+ """Deselects all items of sel_type.
+
+ All items currently selected items are removed from the selection set
+ before the all set is selected.
+
+ Args:
+ sel_type (str): the items to select
+
+ - 'all'
+ - 'parts'
+ - 'areas'
+ - 'lines'
+ - 'points'
+ - 'elements'
+ - 'faces'
+ - 'nodes'
+ """
+ if self.__all == True:
+ self.__parts = self.parts
+ self.__areas = self.areas
+ self.__slines = self.lines
+ self.__points = self.points
+ self.__elements = self.elements
+ self.__faces = self.faces
+ self.__nodes = self.nodes
+ self.__all = False
+ if sel_type == 'all':
+ self.select_none()
+ elif sel_type == 'parts':
+ self.__parts = set()
+ elif sel_type == 'areas':
+ self.__areas = set()
+ elif sel_type == 'lines':
+ self.__slines = set()
+ elif sel_type == 'points':
+ self.__points = set()
+ elif sel_type == 'elements':
+ self.__elements = set()
+ elif sel_type == 'faces':
+ self.__faces = set()
+ elif sel_type == 'nodes':
+ self.__nodes = set()
+
+[docs] def select_none(self):
+ """Empties out the selection object.
+ """
+ self.__all = False
+ self.__parts = set()
+ self.__areas = set()
+ self.__slines = set()
+ self.__points = set()
+ self.__elements = set()
+ self.__faces = set()
+ self.__nodes = set()
+
+[docs] def print_summary(self):
+ """Prints a summary of the number so items selected in each sel_type.
+ """
+ sets = ['parts', 'areas', 'lines', 'points', 'elements', 'faces',
+ 'nodes']
+ spacer = '----------------------------------'
+ print(spacer)
+ print("View's currently selected items:")
+ for sel_type in sets:
+ items = getattr(self, sel_type)
+ print(' %s: %i selected' % (sel_type, len(items)))
+ print(spacer)
+
+