################################################################################ # # This script is intended to be used with all quad8 models using internal solver # Analysys type must be Static 2D. # it estimates the error on each element from the difference between nodal and element Von Mises stress # and selects elements exceeding a tolerance to be refined # tolerance can be chosen as Absolute value in MPa or as a % of maximum Von Mises stress # # The script solves the model and saves it with the names user defines plus iteration number. # After running this script, elements to be refined will be in selection named "Refine" # In order to avoid generation of high valence nodes on refinement, # the scrips adds to named selection "Vertices" vertex nodes that are near to 2 or more # vertices to be refined from diferent elements, # then it Refinex3 to get the new refined model with all quad8 elements. # The new mesh may be smoothed if needed. # # WARNING: # If model has curved boundaries, x3 refinement may create new boundary nodes not exactly # in the right curved surface that may need adjust using Menu Mesh Tools > Project Onto Surface # ################################################################################ import math import sys import os import clr clr.AddReference('System.Windows.Forms') clr.AddReference('System.Drawing') from System.Windows.Forms import * def sort_reverse(lst, indexes, reverse=True): return [val for (_, val) in sorted(zip(indexes, lst), key=lambda x: \ x[0], reverse=reverse)] def is_float(string): try: float(string) return True except ValueError: return False def inputfloat(mensaje): entrada="" while not is_float(entrada): entrada=mw.input(mensaje) return entrada assert mw.version() >= 19 , "Version should be 19 or later" def elements_of_node(n): # this funtion returns the list of elements of a node return Elem_of_node[n] def vertices_of_element(e): # this funtion returns the list of vertices of an element vert=[] for v in [0, 2, 4, 6]: vert.append(mw.nodes(e)[v]) return vert def sides_of_element(e): # this funtion returns the list of sdes nodes of an element vert=[] for v in [0, 2, 4, 6]: vert.append(mw.nodes(e)[v]) return vert def neigh_vertices(n): # this funtion returns the list of vertices neighboring a node (if node is a vertex) neigh_index = [[2, 6], [0, 4], [2, 6], [0, 4]] # neighboring vertices indices depending of node index nv=[] for e in elements_of_node(n): nodos=list(mw.nodes(e)) ind=nodos.index(n) if ind % 2 == 0: # if n is a vertex, it adds neighbors for i in neigh_index[ind/2]: if nodos[i] not in nv: nv.append(nodos[i]) return nv ################################################################################ # End of auxiliary definitions # Script start archivo=mw.file_name() # current file name dlg = SaveFileDialog(); dlg.FileName = archivo; dlg.DefaultExt = "liml"; dlg.Filter = "Mecway files (*.liml)|*.liml"; dlg.Title = "Save this refinement "; if dlg.ShowDialog() != DialogResult.OK: # if Cancel was clicked mw.message("Canceled") sys.exit() archivonuevo = dlg.FileName # generic new file name RoA=mw.input("Input tolerance type\n\nA=Absolut\nR=Relativ\nOther Cancel ") # asks for tolerance type if RoA.upper() == "A" : TolA = float(inputfloat("Input Absolut error tolerancia (MPa):\n 0 to cancel")) # asks for Absolute tolerance if TolA == 0 : mw.message("Canceled") sys.exit() elif RoA.upper() == "R" : porc = float(inputfloat("Input Relative error tolerancia (%):\n 0 to cancel")) # asks for Relative tolerance if RolA == 0 : mw.message("Canceled") sys.exit() else: mw.message("Canceled") sys.exit() Smooth=mw.input("Smooth after refine? \n\nY=Yes\nN=No\nOther Cancel ") # asks for tolerance type if Smooth.upper() == "Y" : SmoothFlag = True elif Smooth.upper() == "N" : SmoothFlag = False else: mw.message("Canceled") sys.exit() elementos= mw.all_elements() # list of all elements for e in elementos: if mw.element_shape(e) != "quad8": # verify if all elements are quad8 mw.message("All elements must be quad8.\n Check element "+str(e)) sys.exit() Max=abs(int(mw.input("iterationations limit ?"))) ################################################################################ # Begins refinement iterationations iteration=1 Stop=False while not Stop and iteration<=Max: allElemNodes=[] for e in mw.all_elements(): allElemNodes.extend(mw.nodes(e)) # creates a list of all element nodes. notused=[] for n in mw.all_nodes(): if n not in allElemNodes: # checks if the node is in any element notused.append(n) if len(notused)>0 : mw.delete_nodes(notused) # deletes unused nodes mw.display() Elem_of_node = [[] for i in range(len(mw.all_nodes())+1)] # list of [list of elements of each node] must be recalculated for each new mesh for e in mw.all_elements(): for n in mw.nodes(e): Elem_of_node[n].append(e) salida=archivonuevo[0:len(archivonuevo)-5]+str(iteration) # new file name without extension. elementos= mw.all_elements() # list of all elements ############################################################################ # solves the model and check if solved OK mw.solve() mw.display() if len(solution.all_nodes()) == 0 : mw.message("Model solution failed") sys.exit() # Exits iteration if solution fails elementos= mw.all_elements() # list of all elements nodos = list(mw.all_nodes()) # list of global nodes without unused nodes allElemNodes=[] for e in mw.all_elements(): allElemNodes.extend(mw.nodes(e)) # creates a list of all element nodes. vertices=[] sides=[] for e in elementos: NofE=list(mw.nodes(e)) for i in [0, 2, 4, 6]: vertices.append(NofE[i]) for i in [1, 3, 5, 7]: sides.append(NofE[i]) for n in nodos: if n in vertices and n in sides: mw.message("Node "+str(n)+" is vertex and side") sys.exit() EPN = [0] * (len(nodos)+1) for n in nodos: EPN[n]=allElemNodes.count(n) # counts how many elements share the node if "Boundary" in mw.named_selections(): mw.delete_named_selection("Boundary") # deletes old "Boundary" selection if exists BoundaryVertices = mw.new_node_selection("Boundary") boundary=[] for s in sides: if EPN[s]==1: e=elements_of_node(s)[0] NofE=list(mw.nodes(e)) i=NofE.index(s) boundary.append(NofE[i-1]) mw.add_to_named_selection("Boundary", NofE[i-1]) if i<7: boundary.append(NofE[i+1]) mw.add_to_named_selection("Boundary", NofE[i+1]) else: boundary.append(NofE[0]) mw.add_to_named_selection("Boundary", NofE[0]) VMn = [] for n in solution.all_nodes(): VMn.append(solution.node_value("vonmises", n)) maxVM=max(VMn)*1e-6 # Maximum node VM stress in MPa EofN = [0] * len(nodos) # crea una lista de ceros del tamano de nodos ErrEst2n = [0] * len(nodos) ErrEstn = [0] * len(nodos) for e in elementos: NofE = mw.nodes(e) ErrEst2e=0 for n in range(0, len(NofE)): ErrNod=abs(solution.node_value("vonmises", NofE[n])-solution.element_node_value("vonmises", e, n+1))*1e-6 ErrEst2e=ErrEst2e+ErrNod**2 #suma los errores al cuadrado de cada nodo del elemento EofN[NofE[n]-1] = e # elemento al que pertenece el nodo # es el ultimo elemento si pertenece a varios ############################################################################################# # a continuacion se calcula para cada elemento el estimador de error en cada uno de sus nodos # por diferencia entre las tensiones elementales y nodales. # los cuadrados de esas diferencias se van promediando para cada nodo global para crear la # variable de campo de error # y a su vez se va calculando para cada elemento la RMS de esas diferencias en sus nodos compartidos. # Luego para los nodos no compartidos (que tienen diferencia nula) se define el estimador de error como # la RMS de las diferencias en los nodos compartidos del elemento que tiene el nodo no compartido. ############################################################################################# err_RMS2Tot=0 vol=0 # zeroes elements volumes sum err_RMSe = [] # creates a list of RMS error estimate in the nodes od each element for e in elementos: cuenta=0 Dif2e = 0 NofE = mw.nodes(e) vole=mw.volume(e) # elements volume vol=vol+vole # sums elements volume for n in range(0, len(NofE)): # for each node in the element if EPN[NofE[n]]>=1 : # if node is shared by several elements err2=(5*(solution.node_value("vonmises", NofE[n])-solution.element_node_value("vonmises", e, n+1))*1e-6)**2 # 5 times (experimental factor) differences between VM nodal and elemental stresses in MPa. Dif2e=Dif2e+err2 # sums squared differences in each node of the element ErrEst2n[NofE[n]-1]=ErrEst2n[NofE[n]-1]+err2/EPN[NofE[n]] # sums squared differences in global node cuenta=cuenta+1 # counts how many element nodes are shared with other elements err_RMSe.append(math.sqrt(Dif2e/cuenta)) # RMS of differences in element shared nodes added to list of element error estimates. err_RMS2Tot=err_RMS2Tot+vole*Dif2e/cuenta # mean squared differences in element shared nodes are volume integrated for n in nodos: if EPN[n]==1: ErrEst2n[n-1]=err_RMSe[EofN[n-1]-1]**2 # if a global node is not shared between several elements, the error estimates in node is equal to the error estimate of the element ErrEstn[n-1] = math.sqrt(ErrEst2n[n-1]) solution.set_variable("ErrorEst_MPa", ErrEstn) # creates Error estimator variable to be shown solution.set_variable("ErrorEst2_MPa2", ErrEst2n) # creates Error estimator sqared variable to be integrated if RoA.upper() == "A" : tol = TolA else: tol = maxVM*porc/100.0 if "Refine" in mw.named_selections(): mw.delete_named_selection("Refine") # deletes old "Refine" selection if exists Refi = mw.new_element_selection("Refine") # creates new "Refine" element selection Arefinar=[] # creates list of elements to be refined for e in elementos: if err_RMSe[e-1] >= tol: # select elements to be refined mw.add_to_named_selection(Refi, e) # adds element to the selection Arefinar.append(e) # adds element to the list## Refinar = len(Arefinar) # quantity of elements to be refined RMSErr = math.sqrt(err_RMS2Tot/vol) # whole model RMS error ############################################################################ # Output text results ############################################################################ f = open(salida+".txt","w") f.write("File "+salida+".txt\n") f.write(str(len(elementos))+" elements and "+str(len(nodos))+" nodes\n") f.write(str(Refinar)+" elements to be refined\n") f.write("Maximum node VM Stress = "+str(maxVM)+" MPa\n") f.write("Maximum element Error = "+str(max(err_RMSe))+" MPa\n") f.write("Error Tolerance = "+str(tol)+" MPa\n") f.write("Total Volume = "+str(vol)+" m3\n") f.write("Error RMS = "+str(RMSErr)+" MPa\n") f.write("\n") f.write("Ord. Elem Error\n") Elemorden=sort_reverse(elementos, err_RMSe) for ord in range(1, Refinar+1): e = Elemorden[ord-1] f.write(str(ord).rjust(4)+str(e).rjust(4)+" "+"{:13.8f}".format(err_RMSe[e-1])+"\n") f.close() os.startfile(salida+".txt") # opens saved text output file if "Vertices" in mw.named_selections(): mw.delete_named_selection("Vertices") # deletes old Vertices selection if exists mw.new_node_selection("Vertices") # creates new Vertices selection if "Vertvecin" in mw.named_selections(): mw.delete_named_selection("Vertvecin") # deletes old Vertvecin selection if exists mw.new_node_selection("Vertvecin") # creates new Vertvecin selection if Refinar == 0 : # Exit iterationation if there are no element to refine mw.save_as(salida+".liml") # saves Mecway model mw.message("No more refinement needed") Stop=True else: Stop=False # Now it improves list of vertices to be refined adding vertices that have 2 or more neighboring vertices (from different elements) to be refined vertref=[] # creates list of vertices to be refined for e in mw.named_selection("Refine"): # for each element to be refined for v in vertices_of_element(e): # for each of its vertex if v not in vertref: # but not in the list of vetices to be refined vertref.append(v) # appends vertex to list of vetices to be refined mw.add_to_named_selection("Vertices", v) # and also to selection of vertex to be refined vertvecin=[] # creates list of vertices neighboring to elements to be refined for e in mw.named_selection("Refine"): # for every element to be refined for v in vertices_of_element(e): # for each of its vertices for nb in neigh_vertices(v): # for each of its neighboring vertices if nb not in vertref and nb not in vertvecin: # if it is not to be refined and not in the list of vertices neighboring to elements to be refined vertvecin.append(nb) # appends nighbor vertex to list of vertices neighboring to elements to be refined mw.add_to_named_selection("Vertvecin", nb) # and also to selection Vertvecin for v in vertvecin: # for each vertex neighboring to elements to be refined vecinosref=[] # creates list of its neighboring vertices to be refined for vc in neigh_vertices(v): # if any neighbor of v if vc in vertref: # is to be refined if vc not in vecinosref: # and it is not allready in the list vecinosref.append(vc) # it is added to the list if len(vecinosref) > 2 : # if vertex v has more than 2 neighbors to be refined, vertref.append(v) # adds v to be refined mw.add_to_named_selection("Vertices", v) if len(vecinosref) == 2 : # if vertex v has 2 neighbors to be refined, arefinar=True for ele in elements_of_node(vecinosref[0]): # for any element of the first neighbor if vecinosref[1] in mw.nodes(ele): # if the other neighbor is in the same element arefinar=False # v will not be refined if arefinar: # but if the 2 neighbors belong to different element, vertref.append(v) # v will be refined mw.add_to_named_selection("Vertices", v) ############################################################################ # Cuando un elemento que no es de Refine tiene vertices a refinar, si los demas vertices son frontera, ponerlos a refinar ############################################################################ for e in elementos: if e not in Refi: # para los elementos no en Refi VoE=vertices_of_element(e) arefinar=True for v in VoE: if v not in vertref and v not in boundary: # si tiene algun vertice no a refinar y no frontera arefinar=False # no se debe refinar if arefinar: for v in VoE: if v not in vertref: vertref.append(v) mw.add_to_named_selection("Vertices", v) mw.save_as(salida+".liml") # saves Mecway model mw.refine_x3(vertref) # refines x3 iteration= iteration+1 if SmoothFlag: mw.smooth(mw.all_elements()) # smoothes mesh mw.smooth(mw.all_elements()) # smoothes mesh if iteration>Max: mw.message("Maximum iterations reached") salida=archivonuevo[0:len(archivonuevo)-5]+str(iteration) mw.save_as(salida+".liml") # saves last refined Mecway model