Skip to content
Snippets Groups Projects
Commit 40e62086 authored by Tscherpel, Tabea's avatar Tscherpel, Tabea
Browse files

Upload New File

parent 957bc4f0
No related branches found
No related tags found
No related merge requests found
#########################
# author: T.Tscherpel
# created 21/06/2023
#########################
from dolfin import *
import numpy as np
from mesh_routines import *
from mshr import *
import matplotlib.pyplot as plt
import tikzplotlib
## define 3D meshes
dim = 3 # space dimension
s = 20.0 # scaling factor
mesh1 = meshcube(2,s) # mesh of cube domain
mesh2 = meshdiamond(50,s) # mesh of diamond shaped domain with 50 slices
refpt = [.1111111,.33333333,0.33333333]
## plot initial meshes
plt.figure(1)
plot(mesh1)
plt.savefig("cube_mesh.png")
tikzplotlib.save("cube_mesh.tex")
plt.figure(2)
plot(mesh2)
plt.savefig("diamond_mesh.png")
tikzplotlib.save("diamond_mesh.tex")
## uniform refinement
# with maximal 5 full refinements and at most 10^6 vertices in the final mesh
vlist1, ebarlist1 = ebar(mesh1,0,6,10e+6,refpt)
vlist2, ebarlist2 = ebar(mesh2,0,6,10e+6,refpt)
plt.figure(3)
plt.semilogx(vlist1,ebarlist1,'-o')
plt.semilogx(vlist2,ebarlist2,'-s')
plt.semilogx(vlist2,14*np.ones(len(vlist1)),'--k')
plt.ylim([7, 15])
plt.legend(['$\mathcal{T}_0 = \mathcal{T}_{F}$',\
'$\mathcal{T}_0 = \mathcal{T}_{D,50}$'],\
loc="upper left")
plt.xlabel('V')
plt.ylabel('$\overline{e}$')
plt.grid()
tikzplotlib.save("ebar_unif.tex")
plt.savefig("ebar_unif.png")
## adaptive refinement (58 times)
vlist3, ebarlist3 = ebar(mesh1,1,58,10e+6,refpt)
vlist4, ebarlist4 = ebar(mesh2,1,58,10e+6,refpt)
plt.figure(4)
plt.semilogx(vlist3,ebarlist3,'-o')
plt.semilogx(vlist4,ebarlist4,'-s')
plt.semilogx(vlist4,14*np.ones(len(vlist4)),'--k')
plt.ylim([7, 15])
plt.legend(['$\mathcal{T}_0 = \mathcal{T}_{F}$',\
'$\mathcal{T}_0 = \mathcal{T}_{D,50}$'],\
loc="upper left")
plt.xlabel('V')
plt.ylabel('$\overline{e}$')
plt.grid()
tikzplotlib.save("ebar_adapt.tex")
plt.savefig("ebar_adapt.png")
plt.show()
####################
##### routines #####
####################
def meshcube(iR,s):
# input: iR initial refinement, s scaling factor leading to cube [0,s]^3
# output: regular mesh of unit cube
mesh = UnitCubeMesh(iR,iR,iR)
x = mesh.coordinates()
x[:, :] *= s
mesh.bounding_box_tree().build(mesh)
return(mesh)
####################
def meshdiamond(m,s):
# input: m number of tetrahedra on the spine, s scaling factor for the whole mesh
#
Pm = np.array([0.0, 0.0, -1.0])
Pp = np.array([0.0, 0.0, 1.0])
P0 = np.array([0.0, 0.0, 0.0])
al = 2.*pi/m
Plist = np.zeros((m,3))
for j in range(0,m):
a = al*j
Plist[j,:] = [cos(a),sin(a),0.0]
editor = MeshEditor()
mesh = Mesh()
editor.open(mesh, 'tetrahedron',3 , 3)
editor.init_vertices(m+3) # number of vertices
editor.init_cells(2*m) # number of cells
editor.add_vertex(0, P0)
editor.add_vertex(1, Pm)
editor.add_vertex(2, Pp)
for j in range(0,m):
editor.add_vertex(j+3, Plist[j,:])
for j in range(0,m-1):
editor.add_cell(2*j,np.array([0, j+3, j+4, 1], dtype=np.uintp))
editor.add_cell(2*j+1,np.array([0, j+3, j+4, 2], dtype=np.uintp))
editor.add_cell(2*(m-1),np.array([0, m+2, 3, 1], dtype=np.uintp))
editor.add_cell(2*m-1,np.array([0, m+2, 3, 2], dtype=np.uintp))
editor.close()
# scale mesh
x = mesh.coordinates()
x[:, :] *= s
mesh.bounding_box_tree().build(mesh)
return(mesh)
####################
def meshquant(mesh):
# input: mesh
# output: v number of vertices, e number of edges
mesh.init()
v = mesh.num_vertices()
e = mesh.num_edges()
return(v,e)
####################
def ebar(mesh,ind,maxref,vmax,refpt):
# input: mesh, ind indicator for uniform (0) or adaptive (1) refinement
# refmax maximal number of refinements,
# vmax maximal number of vertices
# refpt point for adaptive refinement
# output: vlist list of vertices, ebarlist list of ebar values
vlist = []
elist = []
v,e = meshquant(mesh)
vlist.append(v)
elist.append(e)
j = 0
while j < maxref and v < vmax:
if ind == 0:
mesh = refine(mesh) # uniform refinement (Carey-Plaza)
if ind == 1:
mesh = adaptiveref(mesh,refpt) # uniform refinement (Carey-Plaza)
v,e = meshquant(mesh)
vlist.append(v)
elist.append(e)
j = j + 1
ebarlist = np.divide(np.multiply(2,elist),vlist)
return(vlist,ebarlist)
####################
def adaptiveref(mesh,refpt):
# input: mesh, refpt point towards which adaptive refinement is performed
# output: refined mesh
DG = FunctionSpace(mesh,'DG',0)
dg = Function(DG)
nrCells = dg.vector()[:].size
dg.vector()[:] = range(0,nrCells)
val = np.int(dg((refpt)))
dg.vector()[:] = np.zeros(nrCells)
cmark = MeshFunction("bool", mesh, mesh.topology().dim())
cmark.set_all(False)
cell = Cell(mesh,val)
cmark[cell] = True
mesh = refine(mesh,cmark)
return(mesh)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment