Skip to content
Snippets Groups Projects
Select Git revision
  • gmsh_4_13_1
  • master default protected
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • gmsh_4_14_0
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
40 results

adapt_mesh.py

Blame
  • adapt_mesh.py 4.79 KiB
    import numpy as np
    import sys
    import gmsh
    
    
    def triangle_max_edge(x):
        a = np.sum((x[:, 0, :] - x[:, 1, :])**2, 1)**0.5
        b = np.sum((x[:, 0, :] - x[:, 2, :])**2, 1)**0.5
        c = np.sum((x[:, 1, :] - x[:, 2, :])**2, 1)**0.5
        return np.maximum(a, np.maximum(b, c))
    
    
    class Mesh:
        def __init__(self):
            self.vtags, vxyz, _ = gmsh.model.mesh.getNodes()
            self.vxyz = vxyz.reshape((-1, 3))
            vmap = dict({j: i for i, j in enumerate(self.vtags)})
            self.triangles_tags, evtags = gmsh.model.mesh.getElementsByType(2)
            evid = np.array([vmap[j] for j in evtags])
            self.triangles = evid.reshape((self.triangles_tags.shape[-1], -1))
    
    
    def my_function(xyz):
        a = 6 * (np.hypot(xyz[..., 0] - .5, xyz[..., 1] - .5) - .2)
        f = np.real(np.arctanh(a + 0j))
        return f
    
    
    def compute_interpolation_error(nodes, triangles, f):
        uvw, weights = gmsh.model.mesh.getIntegrationPoints(2, "Gauss2")
        jac, det, pt = gmsh.model.mesh.getJacobians(2, uvw)
        numcomp, sf, _ = gmsh.model.mesh.getBasisFunctions(2, uvw, "Lagrange")
        sf = sf.reshape((weights.shape[0], -1))
        qx = pt.reshape((triangles.shape[0], -1, 3))
        det = np.abs(det.reshape((triangles.shape[0], -1)))
        f_vert = f(nodes)
        f_fem = np.dot(f_vert[triangles], sf)
        err_tri = np.sum((f_fem - f(qx))**2 * det * weights, 1)
        return f_vert, np.sqrt(err_tri)
    
    def compute_size_field(nodes, triangles, err, N):
        x = nodes[triangles]
        a = 2.
        d = 2.
        fact = (a**((2. + a) /
                    (1. + a)) + a**(1. / (1. + a))) * np.sum(err**(2. / (1. + a)))
        ri = err**(2. /
                   (2. *
                    (1 + a))) * a**(1. /
                                    (d *
                                     (1. + a))) * ((1. + a) * N / fact)**(1. / d)
        return triangle_max_edge(x) / ri
    
    
    print("Usage: adapt_mesh [initial lc] [target #elements] [dump files]")
    
    lc = 0.02
    N = 10000
    dumpfiles = False
    gui = True
    
    argv = sys.argv
    if '-nopopup' in sys.argv:
        gui = False
        argv.remove('-nopopup')
    
    if len(argv) > 1: lc = float(sys.argv[1])
    if len(argv) > 2: N = int(sys.argv[2])
    if len(argv) > 3: dumpfiles = int(sys.argv[3])
    
    gmsh.initialize()
    
    # create a geometrical gmsh.model
    gmsh.model.add("square")
    square = gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
    gmsh.model.occ.synchronize()
    
    # create intial uniform mesh
    pnts = gmsh.model.getBoundary([(2, square)], True, True, True)
    gmsh.model.mesh.setSize(pnts, lc)
    #gmsh.option.setNumber('Mesh.Algorithm', 6) # Frontal
    gmsh.model.mesh.generate(2)
    if dumpfiles: gmsh.write("mesh.msh")
    mesh = Mesh()
    
    # compute and visualize the interpolation error
    f_nod, err_ele = compute_interpolation_error(mesh.vxyz, mesh.triangles,
                                                 my_function)
    f_view = gmsh.view.add("nodal function")
    gmsh.view.addModelData(f_view, 0, "square", "NodeData", mesh.vtags,
                           f_nod[:, None])
    if dumpfiles: gmsh.view.write(f_view, "f.pos")
    err_view = gmsh.view.add("element-wise error")
    gmsh.view.addModelData(err_view, 0, "square", "ElementData",
                           mesh.triangles_tags, err_ele[:, None])
    if dumpfiles: gmsh.view.write(err_view, "err.pos")
    
    # compute and visualize the remeshing size field
    sf_ele = compute_size_field(mesh.vxyz, mesh.triangles, err_ele, N)
    sf_view = gmsh.view.add("mesh size field")
    gmsh.view.addModelData(sf_view, 0, "square", "ElementData",
                           mesh.triangles_tags, sf_ele[:, None])
    gmsh.plugin.setNumber("Smooth", "View", gmsh.view.getIndex(sf_view))
    gmsh.plugin.run("Smooth")
    if dumpfiles: gmsh.view.write(sf_view, "sf.pos")
    
    # create a new gmsh.model (to remesh the original gmsh.model in-place, the size field
    # should be created as a list-based view)
    gmsh.model.add("square2")
    gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
    gmsh.model.occ.synchronize()
    
    # mesh the new gmsh.model using the size field
    bg_field = gmsh.model.mesh.field.add("PostView")
    gmsh.model.mesh.field.setNumber(bg_field, "ViewTag", sf_view)
    gmsh.model.mesh.field.setAsBackgroundMesh(bg_field)
    #gmsh.option.setNumber('Mesh.Algorithm', 2) # Delaunay
    gmsh.model.mesh.generate(2)
    if dumpfiles: gmsh.write("mesh2.msh")
    mesh2 = Mesh()
    
    # compute and visualize the interpolation error on the adapted mesh
    f2_nod, err2_ele = compute_interpolation_error(mesh2.vxyz, mesh2.triangles,
                                                   my_function)
    f2_view = gmsh.view.add("nodal function on adapted mesh")
    gmsh.view.addModelData(f2_view, 0, "square2", "NodeData", mesh2.vtags,
                           f2_nod[:, None])
    if dumpfiles: gmsh.view.write(f2_view, "f2.pos")
    err2_view = gmsh.view.add("element-wise error on adapated mesh")
    gmsh.view.addModelData(err2_view, 0, "square2", "ElementData",
                           mesh2.triangles_tags, err2_ele[:, None])
    if dumpfiles: gmsh.view.write(err2_view, "err2.pos")
    
    # show everything in the gui
    if gui:
        gmsh.fltk.run()
    
    gmsh.finalize()