diff --git a/pycalculix/feamodel.py b/pycalculix/feamodel.py index 85adf66..c09b06e 100644 --- a/pycalculix/feamodel.py +++ b/pycalculix/feamodel.py @@ -160,6 +160,22 @@ def set_ediv(self, items, ediv): items = self.get_items(items) for line in items: line.set_ediv(ediv) + + 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) def set_units(self, dist_unit='m', cfswitch=False): """Sets the units that will be displayed when plotting. @@ -1424,7 +1440,7 @@ def __read_inp(self, fname): L = line.split(',') L = [int(a.strip()) for a in L] enum = L[0] - nlist = [N.idget(a) for a in L[1:]] + nlist = [Dict_NodeIDs[a] for a in L[1:]] e = mesh.Element(enum, etype, nlist) faces = e.faces E.append(e) @@ -1477,7 +1493,7 @@ def __read_inp(self, fname): sets[set_type][set_name] = [] mode = 'set' f.close() - + # loop through sets and remove empty sets # store sets to delete todel = [] @@ -1569,30 +1585,78 @@ def __read_inp(self, fname): print('Nodes: %i' % len(N)) print('Done reading Calculix/Abaqus .inp file') - def mesh(self, fineness=1.0, mesher='gmsh'): + def mesh(self, size=1.0, meshmode='fineness', mesher='gmsh'): """Meshes all parts. Args: - fineness (float): 0.0001 - 1.0, how fine the mesh is. - - - Low numbers are very fine, higher numbers are coarser. + 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(fineness) + self.__mesh_gmsh(size, meshmode) elif mesher == 'cgx': - self.__mesh_cgx(fineness) + self.__mesh_cgx(size) - def __mesh_gmsh(self, fineness): + def __mesh_gmsh(self, size, meshmode): """Meshes all parts using the Gmsh mesher. Args: - fineness (float): 0.0001 - 1.0, how fine the mesh is. - - Low numbers are very fine, higher numbers are coarser. + + 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 + + """ geo = [] ids = {} @@ -1602,6 +1666,15 @@ def __mesh_gmsh(self, fineness): # 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 @@ -1622,7 +1695,7 @@ def __mesh_gmsh(self, fineness): geo.append(txtline) # set division if we have it - if line.ediv != None: + if line.ediv != None and meshmode=='fineness': ndiv = line.ediv+1 esize = line.length()/line.ediv if self.eshape == 'quad': @@ -1690,8 +1763,8 @@ def __mesh_gmsh(self, fineness): geo.append(txtline) # set the meshing options - geo.append('Mesh.CharacteristicLengthFactor = ' - +str(fineness)+'; //mesh fineness') + if meshmode == 'fineness': + geo.append('Mesh.CharacteristicLengthFactor = '+str(size)+'; //mesh fineness') geo.append('Mesh.RecombinationAlgorithm = 1; //blossom') if self.eshape == 'quad': @@ -1740,13 +1813,26 @@ def __mesh_gmsh(self, fineness): # read in the calculix mesh self.__read_inp(self.fname+'.inp') - def __mesh_cgx(self, fineness): + def __mesh_cgx(self, size, meshmode): """Meshes all parts using the Calculix cgx mesher. Args: - fineness (float): 0.0001 - 1.0, how fine the mesh is. - - Low numbers are very fine, higher numbers are coarser. + + 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 """ fbd = [] comps = [] @@ -1770,8 +1856,8 @@ def __mesh_cgx(self, fineness): cgx_elements['quad2plstrain'] = 'qu8e' cgx_elements['quad1plstrain'] = 'qu4e' - num = 1.0/fineness - emult = int(round(num)) # this converts fineness to mesh multiplier + num = 1.0/size + emult = int(round(num)) # this converts size to mesh multiplier # write all points for point in self.points: diff --git a/pycalculix/geometry.py b/pycalculix/geometry.py index 51e6ed5..749d6f5 100644 --- a/pycalculix/geometry.py +++ b/pycalculix/geometry.py @@ -67,6 +67,7 @@ def __init__(self, x, y, z=0): self.lines = [] self.arc_center = False self.on_part = True + self.esize = None base_classes.Idobj.__init__(self) def get_name(self): @@ -258,6 +259,10 @@ 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 + + def set_esize(self, size): + self.esize = size + class Line(base_classes.Idobj): @@ -352,6 +357,15 @@ def set_ediv(self, ediv): ediv (int): number of required elements on this line """ self.ediv = ediv + + 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) def signed_copy(self, sign): """Returns a SignLine copy of this Line with the passed sign.""" @@ -675,7 +689,11 @@ def set_pt(self, index, point): def set_ediv(self, ediv): """Applies the element divisions onto the parent line.""" self.line.set_ediv(ediv) - + + def set_esize(self, esize): + """Applies the element size onto the parent line.""" + self.line.set_esize(esize) + def set_lineloop(self, lineloop): """Sets the parent LineLoop""" # this is needed to cascade set ediv up to FEA model and down onto @@ -869,7 +887,7 @@ def save_to_points(self): """This stores this line in the line's child points.""" for point in self.allpoints: point.save_line(self) - + def set_ediv(self, ediv): """Sets the number of element divisions on the arc when meshing. @@ -877,7 +895,17 @@ def set_ediv(self, ediv): ediv (int): number of required elements on this arc """ self.ediv = ediv + + 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) + def signed_copy(self, sign): """Returns a SignArc instance of this Arc with the passed sign.""" return SignArc(self, sign) @@ -1256,6 +1284,10 @@ def set_pt(self, ind, point): def set_ediv(self, ediv): """Apply the element divisions onto the parent line.""" self.line.set_ediv(ediv) + + def set_esize(self, esize): + """Applies the element size onto the parent line.""" + self.line.set_esize(esize) def set_lineloop(self, lineloop): """Sets the parent LineLoop""" @@ -1657,7 +1689,7 @@ class Area(base_classes.Idobj): 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 - point (list): a list of all points making the area, excludes arc centers + 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 @@ -1983,6 +2015,16 @@ def set_etype(self, etype): """ self.etype = etype self.set_child_ccxtypes() + + 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) def update(self, line_list): """Updates the area's external lines to the passed list. diff --git a/pycalculix/material.py b/pycalculix/material.py index 7542117..b39aad9 100644 --- a/pycalculix/material.py +++ b/pycalculix/material.py @@ -10,16 +10,22 @@ class Material(base_classes.Idobj): name (str): a unique matl name Attributes: - id (int): material id number - density (float): density in mass/volume units - pratio (float): poisson's ratio (unitless) - youngs (float): young's modulus in Force/area = stress units - 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 + - 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 + -- Thermal expansion is in strain per temperature - dict['data'] = zip of (temp, thermal_expansion) - dict['tzero'] = the temperature zero point @@ -29,33 +35,48 @@ 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 = {} - def set_mech_props(self, density, youngs, pratio): + 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 + - 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 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) + - conductivity (float): Power / (distance-temp) + - spec_heat (float): specific heat (energy/(mass-temp) """ self.conductivity = conductivity self.spec_heat = spec_heat @@ -64,10 +85,10 @@ 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 + - 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 + - temps (list): list of temperatures + - tzero (float): temperature zero point """ self.thermal_exp = {} isnum = (isinstance(alphas, float) or isinstance(alphas, int)) @@ -82,9 +103,13 @@ def ccx(self): """Returns a list of text strings for ccx defining the material.""" res = [] res.append('*MATERIAL,NAME='+self.name) - if self.youngs != None: + + 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)) diff --git a/pycalculix/problem.py b/pycalculix/problem.py index 6f743df..2c1abaf 100644 --- a/pycalculix/problem.py +++ b/pycalculix/problem.py @@ -13,12 +13,12 @@ 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 + - 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: diff --git a/pycalculix/results_file.py b/pycalculix/results_file.py index d7a62e2..06db0e5 100644 --- a/pycalculix/results_file.py +++ b/pycalculix/results_file.py @@ -8,6 +8,7 @@ import matplotlib.pyplot as plt import matplotlib.colors as colors import matplotlib.cm as cmx +import numpy as np from . import base_classes # needed for RESFIELDS from . import mesh @@ -424,6 +425,182 @@ def set_time(self, time): else: print('Time %f is not in the loaded times. Valid times are:') print(self.steps) + + + 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' + + Kargs: + 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) + + 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):