# *************************************************************************** # * Copyright (c) 2017 Bernd Hahnebach * # * * # * This file is part of the FreeCAD CAx development system. * # * * # * 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__ = "FreeCAD FEM import tools" __author__ = "Bernd Hahnebach" __url__ = "http://www.freecadweb.org" ## @package importToolsFem # \ingroup FEM # \brief FreeCAD FEM import tools import FreeCAD from FreeCAD import Console def get_FemMeshObjectMeshGroups( fem_mesh_obj ): """ Get mesh groups from mesh. """ # this method is not really needed. It is used in Fenics mesh only. # there was an exception handling if there was no Group property, but # any FemMesh should have the Group property # if not it would be a bug SMESH return fem_mesh_obj.FemMesh.Groups def get_FemMeshObjectOrder( fem_mesh_obj ): """ Gets element order. Element order counting based on number of nodes on edges. Edge with 2 nodes -> linear elements, Edge with 3 nodes -> quadratic elements, and so on. No edges in mesh -> not determined. (Is this possible? Seems to be a very degenerate case.) If there are edges with different number of nodes appearing, return list of orders. """ presumable_order = None edges = fem_mesh_obj.FemMesh.Edges if edges != (): edges_length_set = list({len(fem_mesh_obj.FemMesh.getElementNodes(e)) for e in edges}) # only need set to eliminate double entries if len(edges_length_set) == 1: presumable_order = edges_length_set[0] - 1 else: presumable_order = [el - 1 for el in edges_length_set] else: Console.PrintMessage( "Found no edges in mesh: Element order determination does not work without them.\n" ) return presumable_order def get_FemMeshObjectDimension( fem_mesh_obj ): """ Count all entities in an abstract sense, to distinguish which dimension the mesh is (i.e. linemesh, facemesh, volumemesh) """ dim = None if fem_mesh_obj.FemMesh.Nodes != (): dim = 0 if fem_mesh_obj.FemMesh.Edges != (): dim = 1 if fem_mesh_obj.FemMesh.Faces != (): dim = 2 if fem_mesh_obj.FemMesh.Volumes != (): dim = 3 return dim def get_FemMeshObjectElementTypes( fem_mesh_obj, remove_zero_element_entries=True ): """ Spit out all elements in the mesh with their appropriate dimension. """ FreeCAD_element_names_dims = { "Node": 0, "Edge": 1, "Hexa": 3, "Polygon": 2, "Polyhedron": 3, "Prism": 3, "Pyramid": 3, "Quadrangle": 2, "Tetra": 3, "Triangle": 2} eval_dict = locals() # to access local variables from eval elements_list_with_zero = [( eval("fem_mesh_obj.FemMesh." + s + "Count", eval_dict), s, d ) for (s, d) in FreeCAD_element_names_dims.items()] # ugly but necessary if remove_zero_element_entries: elements_list = [(num, s, d) for (num, s, d) in elements_list_with_zero if num > 0] else: elements_list = elements_list_with_zero return elements_list def get_MaxDimElementFromList( elem_list ): """ Gets element with the maximal dimension in the mesh to determine cells. """ elem_list.sort(key=lambda t: t[2]) return elem_list[-1] def make_femmesh( mesh_data ): """ makes an FreeCAD FEM Mesh object from FEM Mesh data """ import Fem mesh = Fem.FemMesh() m = mesh_data if ("Nodes" in m) and (len(m["Nodes"]) > 0): FreeCAD.Console.PrintLog("Found: nodes\n") if ( ("Seg2Elem" in m) or ("Seg3Elem" in m) or ("Tria3Elem" in m) or ("Tria6Elem" in m) or ("Quad4Elem" in m) or ("Quad8Elem" in m) or ("Tetra4Elem" in m) or ("Tetra10Elem" in m) or ("Penta6Elem" in m) or ("Penta15Elem" in m) or ("Hexa8Elem" in m) or ("Hexa20Elem" in m) ): nds = m["Nodes"] FreeCAD.Console.PrintLog("Found: elements\n") for i in nds: n = nds[i] mesh.addNode(n[0], n[1], n[2], i) elms_hexa8 = m["Hexa8Elem"] for i in elms_hexa8: e = elms_hexa8[i] mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i) elms_penta6 = m["Penta6Elem"] for i in elms_penta6: e = elms_penta6[i] mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5]], i) elms_tetra4 = m["Tetra4Elem"] for i in elms_tetra4: e = elms_tetra4[i] mesh.addVolume([e[0], e[1], e[2], e[3]], i) elms_tetra10 = m["Tetra10Elem"] for i in elms_tetra10: e = elms_tetra10[i] mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9]], i) elms_penta15 = m["Penta15Elem"] for i in elms_penta15: e = elms_penta15[i] mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9], e[10], e[11], e[12], e[13], e[14]], i) elms_hexa20 = m["Hexa20Elem"] for i in elms_hexa20: e = elms_hexa20[i] mesh.addVolume([ e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9], e[10], e[11], e[12], e[13], e[14], e[15], e[16], e[17], e[18], e[19] ], i) elms_tria3 = m["Tria3Elem"] for i in elms_tria3: e = elms_tria3[i] mesh.addFace([e[0], e[1], e[2]], i) elms_tria6 = m["Tria6Elem"] for i in elms_tria6: e = elms_tria6[i] mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5]], i) elms_quad4 = m["Quad4Elem"] for i in elms_quad4: e = elms_quad4[i] mesh.addFace([e[0], e[1], e[2], e[3]], i) elms_quad8 = m["Quad8Elem"] for i in elms_quad8: e = elms_quad8[i] mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i) elms_seg2 = m["Seg2Elem"] for i in elms_seg2: e = elms_seg2[i] mesh.addEdge([e[0], e[1]], i) elms_seg3 = m["Seg3Elem"] for i in elms_seg3: e = elms_seg3[i] mesh.addEdge([e[0], e[1], e[2]], i) Console.PrintLog( "imported mesh: {} nodes, {} HEXA8, {} PENTA6, {} TETRA4, {} TETRA10, {} PENTA15\n" .format( len(nds), len(elms_hexa8), len(elms_penta6), len(elms_tetra4), len(elms_tetra10), len(elms_penta15) ) ) Console.PrintLog( "imported mesh: {} " "HEXA20, {} TRIA3, {} TRIA6, {} QUAD4, {} QUAD8, {} SEG2, {} SEG3\n" .format( len(elms_hexa20), len(elms_tria3), len(elms_tria6), len(elms_quad4), len(elms_quad8), len(elms_seg2), len(elms_seg3) ) ) else: Console.PrintError("No Elements found!\n") else: Console.PrintError("No Nodes found!\n") return mesh def make_dict_from_femmesh( femmesh ): """ Converts FemMesh into dictionary structure which can immediately used from importToolsFem.make_femmesh(mesh_data) to create a valid FEM mesh. """ # this dict can be easily saved and reloaded by yaml # see importYamlJasonMesh for a implementation mesh_data = {} seg2 = [] seg3 = [] tri3 = [] tri6 = [] quad4 = [] quad8 = [] tet4 = [] tet10 = [] hex8 = [] hex20 = [] pent6 = [] pent15 = [] # associations for lengths of tuples to different # edge, face, and volume elements len_to_edge = {2: seg2, 3: seg3} len_to_face = {3: tri3, 6: tri6, 4: quad4, 8: quad8} len_to_volume = { 4: tet4, 10: tet10, 8: hex8, 20: hex20, 6: pent6, 15: pent15 } # analyze edges for e in femmesh.Edges: t = femmesh.getElementNodes(e) len_to_edge[len(t)].append((e, t)) # analyze faces for f in femmesh.Faces: t = femmesh.getElementNodes(f) len_to_face[len(t)].append((f, t)) # analyze volumes for v in femmesh.Volumes: t = femmesh.getElementNodes(v) len_to_volume[len(t)].append((v, t)) mesh_data = { "Nodes": dict([(k, (v.x, v.y, v.z)) for (k, v) in femmesh.Nodes.items()]), "Seg2Elem": dict(seg2), "Seg3Elem": dict(seg3), "Tria3Elem": dict(tri3), "Tria6Elem": dict(tri6), "Quad4Elem": dict(quad4), "Quad8Elem": dict(quad8), "Tetra4Elem": dict(tet4), "Tetra10Elem": dict(tet10), "Hexa8Elem": dict(hex8), "Hexa20Elem": dict(hex20), "Penta6Elem": dict(pent6), "Penta15Elem": dict(pent15), "Groups": dict([( group_num, ( femmesh.getGroupName(group_num), femmesh.getGroupElements(group_num) ) ) for group_num in femmesh.Groups]) } # no pyr5, pyr13? # no groups? return mesh_data def fill_femresult_mechanical( res_obj, result_set ): """ fills a FreeCAD FEM mechanical result object with result data """ if "number" in result_set: eigenmode_number = result_set["number"] else: eigenmode_number = 0 if "time" in result_set: step_time = result_set["time"] step_time = round(step_time, 2) # if disp exists, fill res_obj.NodeNumbers and # res_obj.DisplacementVectors as well as stress and strain # furthermore the eigenmode number if "disp" in result_set: disp = result_set["disp"] res_obj.DisplacementVectors = list(map((lambda x: x), disp.values())) res_obj.NodeNumbers = list(disp.keys()) # fill res_obj.NodeStressXX etc if they exist in result_set # list values are just added # Should we check if the key in stress and strain dict # is the same as the number in NodeNumbers? if "stress" in result_set: stress = result_set["stress"] Sxx = [] Syy = [] Szz = [] Sxy = [] Sxz = [] Syz = [] # values_S .. stress_tensor .. (Sxx, Syy, Szz, Sxy, Sxz, Syz) for i, values_S in enumerate(stress.values()): Sxx.append(values_S[0]) Syy.append(values_S[1]) Szz.append(values_S[2]) Sxy.append(values_S[3]) Sxz.append(values_S[4]) Syz.append(values_S[5]) res_obj.NodeStressXX = Sxx res_obj.NodeStressYY = Syy res_obj.NodeStressZZ = Szz res_obj.NodeStressXY = Sxy res_obj.NodeStressXZ = Sxz res_obj.NodeStressYZ = Syz # fill res_obj.NodeStrainXX etc if they exist in result_set if "strain" in result_set: strain = result_set["strain"] Exx = [] Eyy = [] Ezz = [] Exy = [] Exz = [] Eyz = [] # values_E .. straintuple .. (Exx, Eyy, Ezz, Exy, Exz, Eyz) for i, values_E in enumerate(strain.values()): Exx.append(values_E[0]) Eyy.append(values_E[1]) Ezz.append(values_E[2]) Exy.append(values_E[3]) Exz.append(values_E[4]) Eyz.append(values_E[5]) res_obj.NodeStrainXX = Exx res_obj.NodeStrainYY = Eyy res_obj.NodeStrainZZ = Ezz res_obj.NodeStrainXY = Exy res_obj.NodeStrainXZ = Exz res_obj.NodeStrainYZ = Eyz # fill Equivalent Plastic strain if they exist if "peeq" in result_set: Peeq = result_set["peeq"] if len(Peeq) > 0: if len(Peeq.values()) != len(disp.values()): # how is this possible? An example is needed! Console.PrintError("PEEQ seams to have exptra nodes.\n") Pe = [] Pe_extra_nodes = list(Peeq.values()) nodes = len(disp.values()) for i in range(nodes): Pe_value = Pe_extra_nodes[i] Pe.append(Pe_value) res_obj.Peeq = Pe else: res_obj.Peeq = list(Peeq.values()) # fill eigenmode number if they exist if eigenmode_number > 0: res_obj.Eigenmode = eigenmode_number # fill res_obj.Temperature if they exist # TODO, check if it is possible to have Temperature without disp # we would need to set NodeNumbers than if "temp" in result_set: Temperature = result_set["temp"] if len(Temperature) > 0: if len(Temperature.values()) != len(disp.values()): Temp = [] Temp_extra_nodes = list(Temperature.values()) nodes = len(disp.values()) for i in range(nodes): # how is this possible? An example is needed! Console.PrintError("Temperature seams to have exptra nodes.\n") Temp_value = Temp_extra_nodes[i] Temp.append(Temp_value) res_obj.Temperature = list(map((lambda x: x), Temp)) else: res_obj.Temperature = list(map((lambda x: x), Temperature.values())) res_obj.Time = step_time # fill res_obj.MassFlow if "mflow" in result_set: MassFlow = result_set["mflow"] if len(MassFlow) > 0: res_obj.MassFlowRate = list(map((lambda x: x), MassFlow.values())) res_obj.Time = step_time # disp does not exist, res_obj.NodeNumbers needs to be set res_obj.NodeNumbers = list(MassFlow.keys()) # fill res_obj.NetworkPressure, disp does not exist, see MassFlow if "npressure" in result_set: NetworkPressure = result_set["npressure"] if len(NetworkPressure) > 0: res_obj.NetworkPressure = list(map((lambda x: x), NetworkPressure.values())) res_obj.Time = step_time return res_obj