1import sys,petsc4py 2petsc4py.init(sys.argv) 3from petsc4py import PETSc 4import numpy as np 5 6 7def sensor(x, y): 8 """ 9 Classic hyperbolic sensor function for testing 10 multi-scale anisotropic mesh adaptation: 11 12 f:[-1, 1]² → R, 13 f(x, y) = sin(50xy)/100 if |xy| > 2π/50 else sin(50xy) 14 15 (mapped to have domain [0,1]² in this case). 16 """ 17 xy = (2*x - 1)*(2*y - 1) 18 ret = np.sin(50*xy) 19 if np.abs(xy) > 2*np.pi/50: 20 ret *= 0.01 21 return ret 22 23 24# Set metric parameters 25h_min = 1.0e-10 # Minimum tolerated metric magnitude ~ cell size 26h_max = 1.0e-01 # Maximum tolerated metric magnitude ~ cell size 27a_max = 1.0e+05 # Maximum tolerated anisotropy 28targetComplexity = 10000.0 # Analogous to number of vertices in adapted mesh 29p = 1.0 # Lᵖ normalization order 30 31# Create a uniform mesh 32OptDB = PETSc.Options() 33dim = OptDB.getInt('dim', 2) 34numEdges = 10 35simplex = True 36plex = PETSc.DMPlex().createBoxMesh([numEdges]*dim, simplex=simplex) 37plex.distribute() 38plex.view() 39viewer = PETSc.Viewer().createVTK("anisotropic_mesh_0.vtk", "w") 40viewer(plex) 41 42# Do four mesh adaptation iterations 43for i in range(1, 5): 44 vStart, vEnd = plex.getDepthStratum(0) 45 46 # Create a P1 sensor function 47 comm = plex.getComm() 48 fe = PETSc.FE().createLagrange(dim, 1, simplex, 1, -1, comm=comm) 49 plex.setField(0, fe) 50 plex.createDS() 51 f = plex.createLocalVector() 52 csec = plex.getCoordinateSection() 53 coords = plex.getCoordinatesLocal() 54 pf = f.getArray() 55 pcoords = coords.getArray() 56 for v in range(vStart, vEnd): 57 off = csec.getOffset(v) 58 x = pcoords[off] 59 y = pcoords[off+1] 60 pf[off//dim] = sensor(x, y) 61 f.setName("Sensor") 62 viewer = PETSc.Viewer().createVTK(f"sensor_{i}.vtk", "w") 63 viewer(f) 64 65 # Recover the gradient of the sensor function 66 dmGrad = plex.clone() 67 fe = PETSc.FE().createLagrange(dim, dim, simplex, 1, -1, comm=comm) 68 dmGrad.setField(0, fe) 69 dmGrad.createDS() 70 g = dmGrad.createLocalVector() 71 plex.computeGradientClementInterpolant(f, g) 72 g.setName("Gradient") 73 viewer = PETSc.Viewer().createVTK(f"gradient_{i}.vtk", "w") 74 viewer(g) 75 76 # Recover the Hessian of the sensor function 77 dmHess = plex.clone() 78 H = dmHess.metricCreate() 79 dmGrad.computeGradientClementInterpolant(g, H) 80 H.setName("Hessian") 81 viewer = PETSc.Viewer().createVTK(f"hessian_{i}.vtk", "w") 82 viewer(H) 83 84 # Obtain a metric by Lᵖ normalization 85 dmHess.metricSetMinimumMagnitude(h_min) 86 dmHess.metricSetMaximumMagnitude(h_max) 87 dmHess.metricSetMaximumAnisotropy(a_max) 88 dmHess.metricSetNormalizationOrder(p) 89 dmHess.metricSetTargetComplexity(targetComplexity) 90 metric = dmHess.metricCreate() 91 det = dmHess.metricDeterminantCreate() 92 dmHess.metricNormalize(H, metric, det) 93 metric.setName("Metric") 94 viewer = PETSc.Viewer().createVTK(f"metric_{i}.vtk", "w") 95 viewer(metric) 96 97 # Call adapt routine - boundary label None by default 98 plex = plex.adaptMetric(metric) 99 plex.distribute() 100 plex.view() 101 102 # Write to VTK file 103 viewer = PETSc.Viewer().createVTK(f"anisotropic_mesh_{i}.vtk", "w") 104 viewer(plex) 105