From 1bd16797cdf4469f1d7dcf51da2df4148ddedb6e Mon Sep 17 00:00:00 2001 From: Manu Date: Tue, 12 Sep 2017 11:38:49 +0900 Subject: [PATCH] improved and fixed gmsh import --- libsrc/include/nginterface_v2_impl.hpp | 2 +- libsrc/meshing/python_mesh.cpp | 7 + python/read_gmsh.py | 304 +++++++++++++++++++++---- 3 files changed, 265 insertions(+), 48 deletions(-) diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index d31fb55d..6c92708b 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -81,7 +81,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const if (mesh->GetDimension() == 3) ret.mat = mesh->GetCD2NamePtr(el.edgenr-1); else - ret.mat = nullptr; + ret.mat = mesh->GetMaterialPtr(el.si); } ret.points.num = el.GetNP(); diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 8d03e358..f280ddac 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -213,6 +213,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) (*instance)[i] = py::extract(vertices[i])(); instance->SetIndex(index); } + else if (py::len(vertices) == 5) + { + new (instance) Element(PYRAMID); + for (int i=0; i < 5; i++) + (*instance)[i] = py::extract(vertices[i])(); + instance->SetIndex(index); + } else if (py::len(vertices) == 6) { new (instance) Element(PRISM); diff --git a/python/read_gmsh.py b/python/read_gmsh.py index f9012868..1afc0e52 100644 --- a/python/read_gmsh.py +++ b/python/read_gmsh.py @@ -1,73 +1,283 @@ from netgen.meshing import * +import re +import shutil +import os -def ReadGmsh (filename): - f = open(filename, 'r') - mesh = Mesh(dim=3) - - pointmap = { } - facedescriptormap = { } - - +def import_gmsh(filename, extension, dim): + """ filename: name of the filename + extension: geo or msh + dim: dimension of the mesh + """ + if extension == "geo": + filename_tmp = check_geofile(filename) + if dim == 1: + os.system("gmsh -1 -saveall -format msh -o ./" + filename + ".msh ./" + filename_tmp + ".geo") + elif dim == 2: + os.system("gmsh -2 -saveall -format msh -o ./" + filename + ".msh ./" + filename_tmp + ".geo") + elif dim == 3: + os.system("gmsh -3 -saveall -format msh -o ./" + filename + ".msh ./" + filename_tmp + ".geo") + + return ReadGmsh(filename) + + +def check_geofile(filename): + mod = False + with open(filename + ".geo", 'r') as f: + with open(filename + "_tmp.geo", 'w') as fw: + line = "start" + while line: + line = f.readline() + if line[0:8] == "Physical": + print('WARNING: Physical groups detected - Physical names are not taken into account.') + print('The physical group is discarded.') + fw.write("//" + line) + mod = True + else: + fw.write(line) + + if mod: + filename_tmp = filename + "_tmp" + else: + filename_tmp = filename + os.remove(filename + "_tmp.geo") + + return filename_tmp + + +def complete_geofile(filename): + line_entity = [] + surf_entity = [] + vol_entity = [] + line_phys_gr = [] + surf_phys_gr = [] + vol_phys_gr = [] + + num_lines = sum(1 for line in open(filename + ".geo", 'r')) + 1 + + with open(filename + ".geo", 'r') as f: + for i in range(num_lines): + line = f.readline() + if not line: + line = [" "] + + if line[0:2] == '//': + pass + + elif line[0:5] == 'Line(': + line_entity += {line.split()[0][5:-1]} + + elif line[0:14] == 'Plane Surface(': + surf_entity += {line.split()[1][8:-1]} + + elif line[0:7] == 'Volume(': + vol_entity += {line.split()[0][7:-1]} + + elif line[0:13] == "Physical Line": + start = line.index('{') + end = line.index('}') + number = re.sub(r"\s", "", line[start + 1:end]).split(sep=",") + line_phys_gr += number + + elif line[0:16] == "Physical Surface": + start = line.index('{') + end = line.index('}') + number = re.sub(r"\s", "", line[start + 1:end]).split(sep=",") + surf_phys_gr += number + + elif line[0:15] == "Physical Volume": + start = line.index('{') + end = line.index('}') + number = re.sub(r"\s", "", line[start + 1:end]).split(sep=",") + vol_phys_gr += number + + line_diff = list(set(line_entity) - set(line_phys_gr)) + surf_diff = list(set(surf_entity) - set(surf_phys_gr)) + vol_diff = list(set(vol_entity) - set(vol_phys_gr)) + filename_geo = filename + if line_diff or surf_diff or vol_diff: + shutil.copy2(filename + '.geo', filename + '_tmp.geo') + filename_geo = filename + "_tmp" + if line_diff: + fid = open(filename + ".geo", 'a') + fid.write('//+\n// Added by NGSolve\n') + fid.write('Physical Line("Remaining_line") = {' + re.sub(r'[\']', '', str(line_diff))[1:-1] + '}; \n') + fid.close() + + if surf_diff: + fid = open(filename + ".geo", 'a') + fid.write('//+\n// Added by NGSolve\n') + fid.write('Physical Surface("Remaining_surf") = {' + re.sub(r'[\']', '', str(surf_diff))[1:-1] + '}; \n') + fid.close() + + if vol_diff: + fid = open(filename + ".geo", 'a') + fid.write('//+\n// Added by NGSolve\n') + fid.write('Physical Volume("Remaining_vol") = {' + re.sub(r'[\']', '', str(vol_diff))[1:-1] + '}; \n') + fid.close() + + return filename_geo + + +def ReadGmsh(filename): + meshdim = 1 + with open(filename + ".msh", 'r') as f: + while f.readline().split()[0] != "$Elements": + pass + nelem = int(f.readline()) + for i in range(nelem): + line = f.readline() + eltype = int(line.split()[1]) + if eltype > 1 and eltype != 15: + meshdim = 2 + if eltype > 3 and eltype != 15: + meshdim = 3 + break + + f = open(filename + ".msh", 'r') + mesh = Mesh(dim=meshdim) + + pointmap = {} + facedescriptormap = {} + namemap = {} + materialmap = {} + bbcmap = {} + while True: line = f.readline() - if line == "": break - - if line.split()[0]=="$Nodes": + if line == "": + break + + if line.split()[0] == "$PhysicalNames": + print('WARNING: Physical groups detected - Be sure to define them for every geometrical entity.') + numnames = int(f.readline()) + for i in range(numnames): + f.readline + line = f.readline() + namemap[int(line.split()[1])] = line.split()[2][1:-1] + + if line.split()[0] == "$Nodes": num = int(f.readline().split()[0]) - print ("reading",num,"nodes") for i in range(num): line = f.readline() - nodenum,x,y,z = line.split()[0:4] - pnum = mesh.Add (MeshPoint(Pnt(float(x),float(y),float(z)))) + nodenum, x, y, z = line.split()[0:4] + pnum = mesh.Add(MeshPoint(Pnt(float(x), float(y), float(z)))) pointmap[int(nodenum)] = pnum - if line.split()[0]=="$Elements": + if line.split()[0] == "$Elements": num = int(f.readline().split()[0]) - print ("reading",num,"elements") - + for i in range(num): line = f.readline().split() elmnum = int(line[0]) elmtype = int(line[1]) numtags = int(line[2]) - tag1 = int(line[3]) - - if elmtype == 2: # 3-node trig + # the first tag is the physical group nr, the second tag is the group nr of the dim + tags = [int(line[3 + k]) for k in range(numtags)] + + if elmtype == 1: # 2-node line + num_nodes = 2 + elif elmtype == 2: # 3-node trig num_nodes = 3 - elif elmtype == 3: # 4-node quad + elif elmtype == 3: # 4-node quad num_nodes = 4 - elif elmtype == 4: # 4-node tet + elif elmtype == 4: # 4-node tet num_nodes = 4 + elif elmtype == 5: # 8-node hex + num_nodes = 8 + elif elmtype == 6: # 6-node prism + num_nodes = 6 + elif elmtype == 7: # 5-node pyramid + num_nodes = 5 + elif elmtype == 15: # 1-node point + num_nodes = 1 else: - print ("element type", elmtype, "not implemented") - - nodenums = line[3+numtags:3+numtags+num_nodes] - nodenums2 = [ pointmap[int(nn)] for nn in nodenums ] + raise Exception("element type", elmtype, "not implemented") + nodenums = line[3 + numtags:3 + numtags + num_nodes] + nodenums2 = [pointmap[int(nn)] for nn in nodenums] - if elmtype in [2,3]: # 2d elements + if elmtype == 1: + if meshdim == 3: + if tags[1] in bbcmap: + index = bbcmap[tags[1]] + else: + index = len(bbcmap) + if len(namemap): + mesh.SetCD2Name(index + 1, namemap[tags[0]]) + else: + mesh.SetCD2Name(index+1, "line" + str(tags[1])) + bbcmap[tags[1]] = index - # generate facedescriptor for boundary conditon number, - # first tag is used as bc-number - # element index maps into facedescriptor array - if tag1 in facedescriptormap.keys(): - fdindex = facedescriptormap[tag1] + elif meshdim == 2: + if tags[1] in facedescriptormap.keys(): + index = facedescriptormap[tags[1]] + else: + index = len(facedescriptormap) + 1 + fd = FaceDescriptor(bc=index) + if len(namemap): + fd.bcname = namemap[tags[0]] + else: + fd.bcname = 'line' + str(tags[1]) + mesh.SetBCName(index - 1, fd.bcname) + mesh.Add(fd) + facedescriptormap[tags[1]] = index else: - fd = FaceDescriptor(bc=tag1) - fdindex = mesh.Add(fd) - facedescriptormap[tag1] = fdindex - - mesh.Add (Element2D(fdindex, nodenums2)) - - if elmtype == 4: # 4-node tet - nodenums2 = [ pointmap[int(nn)] for nn in nodenums ] - mesh.Add (Element3D(tag1, [nodenums2[0],nodenums2[1],nodenums2[3],nodenums2[2]])) - print (nodenums2) - # print (mesh) - # print (pointmap) + if tags[1] in materialmap: + index = materialmap[tags[1]] + else: + index = len(materialmap) + 1 + if len(namemap): + mesh.SetMaterial(index, namemap[tags[0]]) + else: + mesh.SetMaterial(index, "line" + str(tags[1])) + materialmap[tags[1]] = index + + mesh.Add(Element1D(index=index, vertices=nodenums2)) + + if elmtype in [2, 3]: # 2d elements + if meshdim == 3: + if tags[1] in facedescriptormap.keys(): + index = facedescriptormap[tags[1]] + else: + index = len(facedescriptormap) + 1 + fd = FaceDescriptor(bc=index) + if len(namemap): + fd.bcname = namemap[tags[0]] + else: + fd.bcname = "surf" + str(tags[1]) + mesh.SetBCName(index - 1, fd.bcname) + mesh.Add(fd) + facedescriptormap[tags[1]] = index + else: + if tags[1] in materialmap: + index = materialmap[tags[1]] + else: + index = len(materialmap) + 1 + if len(namemap): + mesh.SetMaterial(index, namemap[tags[0]]) + else: + mesh.SetMaterial(index, "surf" + str(tags[1])) + materialmap[tags[1]] = index + + mesh.Add(Element2D(index, nodenums2)) + + if elmtype in [4, 5, 6, 7]: # volume elements + if tags[1] in materialmap: + index = materialmap[tags[1]] + else: + index = len(materialmap) + 1 + if len(namemap): + mesh.SetMaterial(index, namemap[tags[0]]) + else: + mesh.SetMaterial(index, "vol" + str(tags[1])) + materialmap[tags[1]] = index + + nodenums2 = [pointmap[int(nn)] for nn in nodenums] + + if elmtype == 4: + mesh.Add(Element3D(index, [nodenums2[0], nodenums2[1], nodenums2[3], nodenums2[2]])) + elif elmtype in [5, 6, 7]: + mesh.Add(Element3D(index, nodenums2)) return mesh - -# ReadGMSH ("test.gmsh") -