# *************************************************************************** # * * # * Copyright (c) 2016 Frantisek Loeffelmann * # * * # * This program is free software; you can redistribute it and/or modify * # * it under the terms of the GNU Lesser General Public License (LGPL) * # * as published by the Free Software Foundation; either version 2 of * # * the License, or (at your option) any later version. * # * for detail see the LICENCE text file. * # * * # * This program is distributed in the hope that it will be useful, * # * but WITHOUT ANY WARRANTY; without even the implied warranty of * # * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * # * GNU Library General Public License for more details. * # * * # * You should have received a copy of the GNU Library General Public * # * License along with this program; if not, write to the Free Software * # * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * # * USA * # * * # *************************************************************************** __title__ = "FemMesh to Mesh converter" __author__ = "Frantisek Loeffelmann, Ulrich Brammer, Bernd Hahnebach" __url__ = "https://www.freecad.org" ## @package FwmMesh2Mesh # \ingroup FEM import time import FreeCAD # import Mesh """ from os.path import join the_file = join(FreeCAD.getResourceDir(), "examples", "FemCalculixCantilever3D.FCStd") fc_file = FreeCAD.openDocument(the_file) fem_mesh = fc_file.getObject("Box_Mesh").FemMesh # do not remove the _ result = fc_file.getObject("CCX_Results") scale = 1 # displacement scale factor from femmesh import femmesh2mesh out_mesh = femmesh2mesh.femmesh_2_mesh(fem_mesh, result, scale) import Mesh Mesh.show(Mesh.Mesh(out_mesh)) """ # These dictionaries list the nodes, that define faces of an element. # The key is the face number, used internally by FreeCAD. # The list contains the nodes in the element for each face. tetFaces = { 1: [0, 1, 2], 2: [0, 3, 1], 3: [1, 3, 2], 4: [2, 3, 0]} pentaFaces = { 1: [0, 1, 2], 2: [3, 5, 4], 3: [0, 3, 4, 1], 4: [1, 4, 5, 2], 5: [0, 2, 5, 3]} hexaFaces = { # hexa8 or hexa20 (ignoring mid-nodes) 1: [0, 1, 2, 3], 2: [4, 7, 6, 5], 3: [0, 4, 5, 1], 4: [1, 5, 6, 2], 5: [2, 6, 7, 3], 6: [3, 7, 4, 0]} pyraFaces = { # pyra5 or pyra13 (ignoring mid-nodes) 1: [0, 1, 2, 3], 2: [0, 4, 1], 3: [1, 4, 2], 4: [2, 4, 3], 5: [3, 4, 0]} face_dicts = { 4: tetFaces, 5: pyraFaces, 6: pentaFaces, 8: hexaFaces, 10: tetFaces, 13: pyraFaces, 15: pentaFaces, 20: hexaFaces} def femmesh_2_mesh(myFemMesh, myResults=None, myDispScale=1): shiftBits = 20 # allows a million nodes, needs to be higher for more nodes in a FemMesh # This code generates a dict and a faceCode for each face of all elements # All faceCodes are than sorted. start_time = time.process_time() faceCodeList = [] faceCodeDict = {} if myFemMesh.VolumeCount > 0: for ele in myFemMesh.Volumes: element_nodes = myFemMesh.getElementNodes(ele) # print("element_node: ", element_nodes) faceDef = face_dicts[len(element_nodes)] for key in faceDef: nodeList = [] codeList = [] faceCode = 0 shifter = 0 for nodeIdx in faceDef[key]: nodeList.append(element_nodes[nodeIdx]) codeList.append(element_nodes[nodeIdx]) codeList.sort() for node in codeList: faceCode += (node << shifter) # x << n: x shifted left by n bits = Multiplication shifter += shiftBits # print("codeList: ", codeList) faceCodeDict[faceCode] = nodeList faceCodeList.append(faceCode) elif myFemMesh.FaceCount > 0: for ele in myFemMesh.Faces: element_nodes = myFemMesh.getElementNodes(ele) # print("element_node: ", element_nodes) faceDef = {1: [0, 1, 2]} for key in faceDef: nodeList = [] codeList = [] faceCode = 0 shifter = 0 for nodeIdx in faceDef[key]: nodeList.append(element_nodes[nodeIdx]) codeList.append(element_nodes[nodeIdx]) codeList.sort() for node in codeList: faceCode += (node << shifter) # x << n: x shifted left by n bits = Multiplication shifter += shiftBits # print("codeList: ", codeList) faceCodeDict[faceCode] = nodeList faceCodeList.append(faceCode) faceCodeList.sort() allFaces = len(faceCodeList) actFaceIdx = 0 singleFaces = [] # Here we search for faces, which do not have a counterpart. # These are the faces on the surface of the mesh. while actFaceIdx < allFaces: if actFaceIdx < (allFaces - 1): if faceCodeList[actFaceIdx] == faceCodeList[actFaceIdx + 1]: actFaceIdx += 2 else: # print("found a single Face: ", faceCodeList[actFaceIdx]) singleFaces.append(faceCodeList[actFaceIdx]) actFaceIdx += 1 else: FreeCAD.Console.PrintMessage("Found a last Face: {}\n".format(faceCodeList[actFaceIdx])) singleFaces.append(faceCodeList[actFaceIdx]) actFaceIdx += 1 output_mesh = [] if myResults: FreeCAD.Console.PrintMessage("{}\n".format(myResults.Name)) for myFace in singleFaces: face_nodes = faceCodeDict[myFace] dispVec0 = myResults.DisplacementVectors[myResults.NodeNumbers.index(face_nodes[0])] dispVec1 = myResults.DisplacementVectors[myResults.NodeNumbers.index(face_nodes[1])] dispVec2 = myResults.DisplacementVectors[myResults.NodeNumbers.index(face_nodes[2])] triangle = [myFemMesh.getNodeById(face_nodes[0]) + dispVec0 * myDispScale, myFemMesh.getNodeById(face_nodes[1]) + dispVec1 * myDispScale, myFemMesh.getNodeById(face_nodes[2]) + dispVec2 * myDispScale] output_mesh.extend(triangle) # print("my triangle: ", triangle) if len(face_nodes) == 4: dispVec3 = myResults.DisplacementVectors[myResults.NodeNumbers.index(face_nodes[3])] triangle = [myFemMesh.getNodeById(face_nodes[2]) + dispVec2 * myDispScale, myFemMesh.getNodeById(face_nodes[3]) + dispVec3 * myDispScale, myFemMesh.getNodeById(face_nodes[0]) + dispVec0 * myDispScale] output_mesh.extend(triangle) # print("my 2. triangle: ", triangle) else: for myFace in singleFaces: face_nodes = faceCodeDict[myFace] triangle = [myFemMesh.getNodeById(face_nodes[0]), myFemMesh.getNodeById(face_nodes[1]), myFemMesh.getNodeById(face_nodes[2])] output_mesh.extend(triangle) # print("my triangle: ", triangle) if len(face_nodes) == 4: triangle = [myFemMesh.getNodeById(face_nodes[2]), myFemMesh.getNodeById(face_nodes[3]), myFemMesh.getNodeById(face_nodes[0])] output_mesh.extend(triangle) # print("my 2. triangle: ", triangle) end_time = time.process_time() FreeCAD.Console.PrintMessage( "Mesh by surface search method: {}\n".format(end_time - start_time) ) return output_mesh