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