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 DalcinOptDB = PETSc.Options() 7*55a74a43SLisandro Dalcin 8*55a74a43SLisandro Dalcindim = OptDB.getInt('dim', 2) 9*55a74a43SLisandro Dalcinplex = PETSc.DMPlex().createBoxMesh([4]*dim, simplex=True) 10*55a74a43SLisandro Dalcinplex.distribute() 11*55a74a43SLisandro Dalcinplex.view() 12*55a74a43SLisandro Dalcin 13*55a74a43SLisandro Dalcin# Create two metric tensor fields corresponding to uniform mesh sizes of 0.1 and 0.2 14*55a74a43SLisandro Dalcinmetric1 = plex.metricCreateUniform(100.0) 15*55a74a43SLisandro Dalcinmetric2 = plex.metricCreateUniform(25.0) 16*55a74a43SLisandro Dalcin 17*55a74a43SLisandro Dalcin# The metrics can be combined using intersection, the result of which corresponds to 18*55a74a43SLisandro Dalcin# the maximum ellipsoid at each point 19*55a74a43SLisandro Dalcinmetric = plex.metricCreate() 20*55a74a43SLisandro Dalcinplex.metricIntersection2(metric1, metric2, metric) 21*55a74a43SLisandro Dalcinmetric1.axpy(-1, metric) 22*55a74a43SLisandro Dalcinassert np.isclose(metric1.norm(), 0.0) 23*55a74a43SLisandro Dalcin 24*55a74a43SLisandro Dalcin# Call adapt routine - boundary label None by default 25*55a74a43SLisandro Dalcinnewplex = plex.adaptMetric(metric) 26*55a74a43SLisandro Dalcinnewplex.view() 27*55a74a43SLisandro Dalcin 28*55a74a43SLisandro Dalcin# Write to VTK file 29*55a74a43SLisandro Dalcinviewer = PETSc.Viewer().createVTK('base_mesh.vtk', 'w') 30*55a74a43SLisandro Dalcinviewer(plex) 31*55a74a43SLisandro Dalcinviewer = PETSc.Viewer().createVTK('isotropic_mesh.vtk', 'w') 32*55a74a43SLisandro Dalcinviewer(newplex) 33