Mesh ve 2-D

generovani geometrie

kruh, vyšší hustota ve středu [Point(3)]

    cl__1 = 0.1;
    Point(1) = {-0.5, 0, 0, cl__1};
    Point(2) = {0, 0.5, 0, cl__1};
    Point(3) = {0, 0, 0, 0.02};
    Point(4) = {0, -0.5, 0, cl__1};
    Point(5) = {0.5, 0, 0, cl__1};
    Circle(1) = {1, 3, 2};
    Circle(2) = {2, 3, 5};
    Circle(3) = {5, 3, 4};
    Circle(4) = {4, 3, 1};
    Line(5) = {1, 3};
    Line(6) = {3, 5};
    Line Loop(8) = {1, 2, -6, -5};
    Plane Surface(8) = {8};
    Line Loop(10) = {4, 5, 6, 3};
    Plane Surface(10) = {10};

generování skriptově pomocí gmsh

geometry='''
cl__1 = %.4f;
Point(1) = {-0.5, 0, 0, cl__1};
Point(2) = {0, 0.5, 0, cl__1};
Point(3) = {0, 0, 0, %.4f};
Point(6) = {%.3f, %.3f, 0, %.4f};
Point(4) = {0, -0.5, 0, cl__1};
Point(5) = {0.5, 0, 0, cl__1};
Circle(1) = {1, 3, 2};
Circle(2) = {2, 3, 5};
Circle(3) = {5, 3, 4};
Circle(4) = {4, 3, 1};
Line(5) = {1, 6};
Line(6) = {6, 5};
Line Loop(8) = {1, 2, -6, -5};
Plane Surface(8) = {8};
Line Loop(10) = {4, 5, 6, 3};
Plane Surface(10) = {10};
'''
xc,yc=0.3,0
cl1=0.2
cl2=0.05
clint=(0.5*cl2+xc*cl1)/(0.5+xc)

open("work/excent3.geo","w").write(geometry%(cl1,clint,xc,yc,cl2))

450

import subprocess as sp
#wdir="C:/Users/Admin/Dropbox/Notebook/Temp/"
wdir="/home/jupyter/note/FEM/"
#gmsh_run='C:/Program Files/Gmsh/gmsh.exe'
gmsh_run="/usr/bin/gmsh"
pars=[gmsh_run,'work/excent3.geo','-2','-o','work/excent3.msh']
#pars+=['-clmin',str(0.005)]
sp.run(pars,cwd=wdir)

CompletedProcess(args=['/usr/bin/gmsh', 'work/excent3.geo', '-2', '-o', 'work/excent3.msh'], returncode=0)

soubor .msh vytvořen v textovém formátu = lze načíst ručně

import numpy as np
import sys
#sys.path.append("/home/limu/Code/Monty")
import fem_laplace as fl
%matplotlib inline
from matplotlib import pyplot as pl

raději ale použijeme udržovanou knihovnu

import meshio

remesh=meshio.read("work/excent.msh")
pts=remesh.points
pts=remesh.points[:,:2]
pl.plot(pts[:,0],pts[:,1],'.')

[<matplotlib.lines.Line2D at 0x7f2c35e979b0>]

len(remesh.cells_dict['triangle'])

590

def pplot(cor,mshi=0.2):
    tri=pts[cor]
    cen=tri.mean(0)
    tri=tri*(1-mshi)+cen*mshi
    for i in range(3):
        pl.plot([tri[i][0],tri[(i+1)%3][0]], [tri[i][1],tri[(i+1)%3][1]],'r')
#jenom nekktere elementy
for c in remesh.cells_dict['triangle'][::3]:
    pplot(c)
np.min(remesh.cells_dict['triangle'])

0

dist2=((remesh.points[:,0]-xc)**2+(remesh.points[:,1]-yc)**2)
centid=dist2.argmin()
centid

208

pts[centid],xc,yc

(array([0.31859877, 0.01728179]), 0.3, 0)

#body sousedici se zdrojovym bodem
faces=remesh.cells_dict['triangle']
centfaces=np.where(np.sum(faces==centid,axis=1))[0]
faces[centfaces]

array([[158, 208, 295],
       [ 60, 295, 208],
       [ 60, 208, 254],
       [158, 304, 208],
       [208, 304, 254]], dtype=int32)

pts.shape

(313, 2)

surf=[]
for cor in remesh.cells_dict['triangle']:
    fac=pts[cor].copy()
    #fac[:,2]=1
    fac=fac[1:]-fac[0]
    surf.append(abs(np.linalg.det(fac)/2))
pl.hist(surf,20);

#rychlejsi metoda (narocnejsi na pamet)
allcor=pts[remesh.cells_dict['triangle']]
#allcor[:,:,2]=1
lencor=allcor[:,1:]-allcor[:,0].reshape(len(allcor),1,2)
max(abs(surf2-surf))

0.0

ind=[[0,1],[1,2],[2,0]]
leng=[]
for cor in remesh.cells_dict['triangle']:
    fac=pts[cor].copy()
    leng.append(np.sqrt([((fac[d][1]-fac[d][0])**2).sum() for d in ind]))
leng=np.array(leng)
leng.shape

(590, 3)

zz=np.r_[1,-1].reshape(2,1)
leng2=np.sqrt([((allcor[:,d]*zz).sum(1)**2).sum(1) for d in ind]).T

pl.hist(leng.flat,20);
np.where(leng>0.1)

(array([ 10,  11,  12,  13,  24,  25,  58,  59, 219, 219, 220, 220, 221,
        222, 235, 235, 236, 237, 238, 241, 245, 245, 246, 247, 248, 273,
        274, 279, 279, 279, 286, 287, 297, 297, 303, 304, 307, 308, 319,
        359, 378, 416, 417, 509]),
 array([1, 1, 2, 0, 1, 2, 2, 0, 1, 2, 1, 2, 0, 2, 1, 2, 1, 2, 0, 1, 0, 1,
        2, 2, 0, 0, 2, 0, 1, 2, 0, 2, 0, 1, 0, 2, 2, 0, 2, 0, 2, 2, 0, 2]))

#pomer nejdelsi a nejkratsi strany trojuhelnika (hledame prilis "spicate" elementy) 
rat=leng.min(axis=1)/leng.max(axis=1)
pl.hist(rat);

#umime secist plochu?
sum(surf),np.pi*0.5**2

(0.7803612880645129, 0.7853981633974483)

#body na okraji (totez co v uvodu)
dist=np.sqrt(pts[:,0]**2+pts[:,1]**2)
amax=dist.max()
edgsel=dist>amax-0.01
edgcor=pts[edgsel]
pl.plot(edgcor[:,0],edgcor[:,1],'.')
sum(edgsel)

64

#sousedni body

nint=len(pts)-sum(edgsel) #pocet internich bodu
edgind=np.r_[:len(edgsel)][edgsel]
nearcor=[cor for cor in remesh.cells_dict['triangle'] if np.any([a in edgind for a in cor])]

nearind=set(np.array(nearcor).ravel())-set(edgind)
#nearind=[i for i in boundind if not i in edgind]
nearind=list(nearind)
len(nearind)

74

edgcor=pts[edgsel]
pl.plot(edgcor[:,0],edgcor[:,1],'.')
nearcor=pts[nearind]
pl.plot(nearcor[:,0],nearcor[:,1],'.')

[<matplotlib.lines.Line2D at 0x24da397e0f0>]

ukol