FENICSx¶
Docker¶
Jak se dostat k FENICSxu přes Docker a mít přístup k aktuálnímu adresáři:
docker run -ti --name <název kontejneru> -v $(pwd):/root/shared dolfinx/dolfinx
Zastavení a návrat z kontejneru:
exit
Výpis všech kontejnerů:
docker ps -a
Spuštění a zastavení existujících kontejnerů:
docker start <název nebo id kontejneru>
docker stop <název nebo id kontejneru>
Vymazání zastavených kontenjenrů:
docker rm <název nebo id kontejneru 1> <název nebo id kontejneru 2> ...
Vstup do bashe spuštěného kontejneru, bash se otevře v adresáři ./root:
docker exec -ti <název nebo id kontejneru> bash
Kopírování souborů z lokálního adresáře do adresáře ./root kontejneru.
docker cp ./complex.py <název nebo id kontejneru>:./root
Kopírování souborů z kontejneru z adresáře ./root do lokálního adresáře:
docker cp <název nebo id kontejneru>:./root ./
Realný a komplexní řešič FENICSx¶
FENICSx se může přepínat pro reálný a komplexní obor úloh. K přepnutí mezi oběma řešiči se provede následujícími příkazy přímo v Dockeru:
source dolfinx-real-mode
source dolfinx-complex-mode
FENICSx 0.10 - Homogenní trhlina v piezoelektrickém ortotropním materiálu¶
from dolfinx import mesh,fem,plot,io,default_scalar_type,default_real_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import ufl,basix
import numpy as np
#-----[ materialy ]-----
C_11=np.zeros((2,1))
C_12=np.zeros((2,1))
C_23=np.zeros((2,1))
C_22=np.zeros((2,1))
C_44=np.zeros((2,1))
C_66=np.zeros((2,1))
e_11=np.zeros((2,1))
e_12=np.zeros((2,1))
e_26=np.zeros((2,1))
omega_11=np.zeros((2,1))
omega_22=np.zeros((2,1))
#-----[ ]-----
#-----[ PZT-5H ]-----
#-----[ ]-----
#-----[ elasticke konstanty >Pa< ]-----
C_11=11.7e10
C_12=5.30e10
C_23=5.50e10
C_22=12.6e10
C_44=3.53e10
C_66=(C_22-C_23)/2
#-----[ piezoelectric coefficients >C/m2< ]-----
e_12=-6.50
e_11=23.30
e_26=17.00
#-----[ dielectric permittivity >10e-9 C/Vm< ]-----
omega_11=13.0e-9
omega_22=15.1e-9
#-----[ ]-----
#-----[ PZT-4 ]-----
#-----[ ]-----
#-----[cislo materialu]-----
#m1 = 1
#-----[elasticke konstanty >Pa<]-----
#C_11[m1,0]=11.3e10
#C_12[m1,0]=7.43e10
#C_23[m1,0]=7.78e10
#C_22[m1,0]=13.9e10
#C_44[m1,0]=2.56e10
#C_66[m1,0]=(C_22[m1,0]-C_23[m1,0])/2
#-----[piezoelectric coefficients >C/m2<]-----
#e_12[m1,0]=-6.98
#e_11[m1,0]=13.84
#e_26[m1,0]=13.44
#-----[dielectric parmittivity >10e-9 C/Vm<]-----
#omega_11[m1,0]=5.47e-9
#omega_22[m1,0]=6.00e-9
#-----[ ]-----
#-----[ BaTiO3 ]-----
#-----[ ]-----
#-----[cislo materialu]-----
#m1 = 1
#-----[elasticke konstanty >Pa<]-----
#C_11[m1,0]=14.6e10
#C_12[m1,0]=6.60e10
#C_23[m1,0]=6.60e10
#C_22[m1,0]=15.0e10
#C_44[m1,0]=4.40e10
#C_66[m1,0]=(C_22[m1,0]-C_23[m1,0])/2
#-----[piezoelectric coefficients >C/m2<]-----
#e_12[m1,0]=-4.35
#e_11[m1,0]=17.50
#e_26[m1,0]=11.40
#-----[dielectric parmittivity >10e-9 C/Vm<]-----
#omega_11[m1,0]=11.2e-9
#omega_22[m1,0]=9.87e-9
#-----[ ]-----
#-----[ PZT-7A ]-----
#-----[ ]-----
#-----[cislo materialu]-----
#m2=1
#-----[elasticke konstanty >Pa<]-----
#C_11[m2,0]=13.1e10
#C_12[m2,0]=7.42e10
#C_23[m2,0]=7.62e10
#C_22[m2,0]=14.8e10
#C_44[m2,0]=2.54e10
#C_66[m2,0]=(C_22[m2,0]-C_23[m2,0])/2
#-----[piezoelectric coefficients >C/m2<]-----
#e_12[m2,0]=-2.10
#e_11[m2,0]=9.50
#e_26[m2,0]=9.70
#-----[dielectric permittivity >10e-9 C/Vm<]-----
#omega_11[m2,0]=7.35e-9
#omega_22[m2,0]=8.11e-9
#-----[ matice voigtovy notace ]-----
C=np.zeros((6,6))
C[0,0]=C_11
C[0,1]=C_12
C[0,2]=C_12
C[1,1]=C_22
C[1,2]=C_23
C[2,2]=C_22
C[3,3]=C_66
C[4,4]=C_44
C[5,5]=C_44
#-----[ podminky symetrie ]------
C[1,0]=C[0,1]
C[2,0]=C[0,2]
C[2,1]=C[1,2]
C[3,0]=C[0,3]
C[3,1]=C[1,3]
C[3,2]=C[2,3]
C[4,0]=C[0,4]
C[4,1]=C[1,4]
C[4,2]=C[2,4]
C[4,3]=C[3,4]
C[5,0]=C[0,5]
C[5,1]=C[1,5]
C[5,2]=C[2,5]
C[5,3]=C[3,5]
C[5,4]=C[4,5]
e_=np.zeros((3,6))
e_[0,0]=e_11
e_[0,1]=e_12
e_[0,2]=e_12
e_[2,4]=e_26
e_[1,5]=e_26
omega_=np.zeros((3,3))
omega_[0,0]=omega_11
omega_[1,1]=omega_22
omega_[2,2]=omega_22
#-----[ ]-----
#-----[ GMSH ]-----
#-----[ ]-----
#-----[ vytvoreni site pomoci gmsh ]-----
#-----[ princip je stejny jako ]-----
#-----[ ve skriptu GMSH ]-----
import gmsh
gmsh.initialize()
gmsh.model.add("t1")
lc=1e-2
gmsh.model.geo.addPoint( .0, .0, .0, lc, 1)
gmsh.model.geo.addPoint( .2, .0, .0, lc, 2)
gmsh.model.geo.addPoint( .2, .2, .0, lc, 3)
gmsh.model.geo.addPoint(-.2, .2, .0, lc, 4)
gmsh.model.geo.addPoint(-.2, .0, .0, lc, 5)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 5, 4)
gmsh.model.geo.addLine(5, 1, 5)
gmsh.model.geo.addCurveLoop([1, 2, 3, 4, 5], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [1], 1)
gmsh.model.addPhysicalGroup(1, [2], 2)
gmsh.model.addPhysicalGroup(1, [3], 3)
gmsh.model.addPhysicalGroup(1, [4], 4)
gmsh.model.addPhysicalGroup(1, [5], 5)
gmsh.model.addPhysicalGroup(2, [1], 1)
#-----[ generace site a souboru *.geo_unrolled a *msh ]-----
#-----[ soubory nejsou nutne k nacteni do FENICSXu ]-----
gmsh.model.mesh.generate(2)
gmsh.write("t1.geo_unrolled")
gmsh.write("t1.msh")
#-----[ nacteni gmsh site do FENICSxu ]-----
#-----[ a ukonceni GMSH ]-----
gmsh_model_rank=0
mesh_comm=MPI.COMM_WORLD
mesh_data=io.gmsh.model_to_mesh(gmsh.model,mesh_comm,gmsh_model_rank,gdim=2)
gmsh.finalize()
#-----[ ]-----
#-----[ FENICSx ]-----
#-----[ ]-----
#-----[ matice poddajnosti a tuhosti ]-----
C=ufl.as_matrix([[C[0,0],C[0,1],0.],[C[1,0],C[1,1],0.],[0.,0.,C[5,5]]])
print('matice tuhosti materilu:')
print('C=')
print(C)
e=ufl.as_matrix([[e_[0,0],e_[0,1],e_[0,5]],[e_[1,0],e_[1,1],e_[1,5]]])
print('piezoelektricke koeficienty materialu:')
print('e=')
print(e)
omega=ufl.as_matrix([[omega_[0,0],omega_[0,1]],[omega_[1,0],omega_[1,1]]])
print('permitivity dialektrik:')
print('omega=')
print(omega)
#-----[ definice velicin pruznosti ]-----
def eps(u):
Du=ufl.grad(u)
return 0.5*(Du+Du.T)
def E(phi):
Ephi=-ufl.grad(phi)
return Ephi
def strain2voigt(e):
return ufl.as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s):
return ufl.as_tensor([[s[0],s[2]],[s[2],s[1]]])
def sigma(u,phi):
su=voigt2stress(ufl.dot(C,strain2voigt(eps(u))))
sphi=voigt2stress(ufl.dot(e.T,E(phi)))
return su-sphi
def D(u,phi):
Du=ufl.dot(e,strain2voigt(eps(u)))
Dphi=ufl.dot(omega,E(phi))
return Du+Dphi
#-----[ objekt konecnoprvkove oblasti >domain< ]-----
#-----[ a hran prvku >ft< ]-----
domain=mesh_data.mesh
ft=mesh_data.facet_tags
#-----[ definice plosneho integralu "dx" pres celou oblast >domain< ]-----
#-----[ a definice krivkoveho integralu "ds" po hranach prvku >ft< ]-----
#-----[ s fyzikalnim GMSH indexem >3< ]-----
dx=ufl.Measure("dx",domain=domain)
ds3=ufl.Measure("ds",subdomain_data=ft,subdomain_id=3)
#-----[ definice objemoveho a plosneho zatizeni jako ]-----
#-----[ konstantich vektoru pres celou oblast >domain<, ]-----
#-----[ ktere musi byt typu >default_scalar_type< ]-----
T=fem.Constant(domain,default_scalar_type((0,1.e6)))
f=fem.Constant(domain,default_scalar_type((0,0)))
#-----[ definice protoru funkci >V< ]-----
#-----[ dimenze ulohy je >domain.geometry.dim=2< ]-----
#-----[ posuvy >u< jsou vektorovy >shape=(2,)< ]-----
#-----[ trojuhelnikovy >domain.topology.cell_name() -> "trinagle"< ]-----
#-----[ prvek >V_el< s >"Lagrange"< polynomy radu >2< ]-----
#-----[ el.potencial >phi< je skalarni ]-----
#-----[ trojuhelnikovy >domain.topology.cell_name() -> "triangle"< ]-----
#-----[ prvek >W_el< s >"Lagrange"< polynomy radu >1< ]-----
#-----[ oba prvky jsou namixovany na prvek >VW_el< ]-----
#-----[ a z nej vytvoren prostor bazovych funkci >V< ]-----
dim=domain.geometry.dim
V_el=basix.ufl.element("Lagrange",domain.topology.cell_name(),2,shape=(dim,))
W_el=basix.ufl.element("Lagrange",domain.topology.cell_name(),1)
VW_el=basix.ufl.mixed_element([V_el,W_el])
V=fem.functionspace(domain,VW_el)
#-----[ podprostory funkci ]-----
#-----[ >V.sub(0)< je podprostor vektoru posuvu >u< prostoru >V< ]-----
#-----[ >funkce collapse()< vytvori z podprostoru >V.sub(0)< novy ]-----
#-----[ prostor >V0< s odpovidajicimi dofs >V0_dofs< ]-----
#-----[ >V.sub(0).sub(0)< a >V.sub(0).sub(1)< jsou podprostory ]-----
#-----[ slozek >u_x< a >u_y< vektoroveho podprostoru >V.sub(0)< ]-----
#-----[ funkce >collapse()< opet vytvori z techto podprostoru nove ]-----
#-----[ prostory >V0x< a >V0y< slozek vektoru posuvu >u< spolu s ]-----
#-----[ odpovidajicimi dofs >V0x_dofs< a >V0y_dofs<, podobne podprostor ]-----
#-----[ >V.sub(1)< odpovidajici el.potencialu je preveden na prostor ]-----
#-----[ >V1< se stupni volnosti >V1_dofs< ]-----
V0,V0_dofs=V.sub(0).collapse()
V0x,V0x_dofs=V.sub(0).sub(0).collapse()
V0y,V0y_dofs=V.sub(0).sub(1).collapse()
V1,V1_dofs=V.sub(1).collapse()
#-----[ hledana a testovaci funkce ]-----
(u_trial,phi_trial)=ufl.TrialFunctions(V)
(u_test,phi_test)=ufl.TestFunctions(V)
#-----[ >facets_dim< je dimenze hran prvku ]-----
facets_dim=domain.topology.dim-1
#-----[ dirichletova okrajova podminka >bc_displ_0< ]-----
#-----[ posuvu >u_x=0< a >u_y=0< v bode [0,0] ]-----
#-----[ >pin_zero(x)< je funkce definujici bod o souradnici [0,0] ]-----
#-----[ funkce >locate_dofs_geometrical()< vytahne "geometricky" ]-----
#-----[ dofs z podprostoru a prostoru posuvu >V.sub(0)< a >V0< na zaklade ]-----
#-----[ funkce >pin_zero(x)<, >fixed_displ_0< je konecnoprvkova funkce jako ]-----
#-----[ interpolace nulovych hodnot posuvu (dofs) v uzlovych bodech konecnoprvkove ]-----
#-----[ site, nulove hodnoty posuvu v techto bodech jsou dane funkci >displ_0(x)< ]-----
def pin_zero(x):
return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],0.0))
def displ_0(x):
return np.vstack((np.full(x.shape[1],0.*x[0]+0.*x[1]),
np.full(x.shape[1],0.*x[0]+0.*x[1])))
dof_V=fem.locate_dofs_geometrical(V=(V.sub(0),V0),marker=pin_zero)
fixed_displ_0=fem.Function(V0)
fixed_displ_0.interpolate(displ_0)
bc_displ_0=fem.dirichletbc(value=fixed_displ_0,dofs=dof_V,V=V.sub(0))
#-----[ dirichletova okrajova podminka >bc_displ_1< ]-----
#-----[ posuvu >u_y=0< podel hranice >1<, tj. >x[0]>0< a >x[1]=0< ]-----
#-----[ hrany prvku podel teto hranice se mouhou urcit bud geometricky kombinaci ]-----
#-----[ funkce >pin_1(x)< a funkce >locate_entities_boundary< nebo ]-----
#-----[ pomoci physical tags z gmsh pomoci funkce >find(physical_tag)<, ]-----
#-----[ funkce >locate_dofs_topological()< vytahne "topologicky" ]-----
#-----[ dofs z podprostoru a prostoru posuvu >V.sub(0).sub(1)< a >V0y< na zaklade ]-----
#-----[ cisel hran prvku >boundary_facets_1<, ]-----
#-----[ >fixed_displ_y< je konecnoprvkova funkce na prostoru >V0y< jako ]-----
#-----[ interpolace nulovych hodnot >displ_y_0(x)< posuvu >u_y< (dofs) v uzlovych bodech ]-----
#-----[ hran prvku na hranici >1<, konecnoprvkove ]-----
#-----[ site, nulove hodnoty posuvu v techto bodech jsou dane funkci >displ_y_0(x)< ]-----
#bud
def pin_1(x):
return np.logical_and(x[0]>=0,np.isclose(x[1],0.0))
boundary_facets_1=mesh.locate_entities_boundary(msh=domain,dim=facets_dim,marker=pin_1)
#nebo
boundary_facets_1=mesh_data.facet_tags.find(1)
def displ_y_0(x):
return np.array([0.0*x[0],])
dofs_V_1=fem.locate_dofs_topological(V=(V.sub(0).sub(1),V0y),entity_dim=facets_dim,entities=boundary_facets_1)
fixed_displ_y=fem.Function(V0y)
fixed_displ_y.interpolate(displ_y_0)
bc_displ_1=fem.dirichletbc(value=fixed_displ_y,dofs=dofs_V_1,V=V.sub(0).sub(1))
#-----[ dirichletova okrajova podminka >bc_potenc< ]-----
#-----[ el.potencialu >phi=0< podel hranice >1<, tj. >x[0]>0< a >x[1]=0< ]-----
#-----[ hrany prvku podel teto hranice se mouhou urcit bud geometricky kombinaci ]-----
#-----[ funkce pin_1(x) a funkce >locate_entities_boundary< nebo ]-----
#-----[ pomoci physical tags z gmsh pomoci funkce >find(physical_tag)<, ]-----
#-----[ stejne jako u dirichletovy okrajove podminku >bc_displ_1<, ]-----
#-----[ funkce >locate_dofs_topological()< vytahne "topologicky" ]-----
#-----[ dofs z podprostoru a prostoru posuvu >V.sub(1)< a >V1< na zaklade ]-----
#-----[ cisel hran prvku >boundary_facets_1<, ]-----
#-----[ fixed_potenc je konecnoprvkova funkce na prostoru >V1< jako ]-----
#-----[ interpolace nulovych hodnot potenc_0(x) el.potencialu >phi< (dofs) ]-----
#-----[ v uzlovych bodech hran prvku na hranici >1<, konecnoprvkove ]-----
#-----[ site, nulove hodnoty posuvu v techto bodech jsou dane funkci potenc_0(x) ]-----
def potenc_0(x):
return np.array([0.0*x[0],])
dofs_W=fem.locate_dofs_topological(V=(V.sub(1),V1),entity_dim=facets_dim,entities=boundary_facets_1)
fixed_potenc=fem.Function(V1)
fixed_potenc.interpolate(potenc_0)
bc_potenc=fem.dirichletbc(value=fixed_potenc,dofs=dofs_W,V=V.sub(1))
#-----[ sestaveni dirichletovych okrajovych podminek do pole ]-----
bcs=[bc_displ_0,bc_displ_1,bc_potenc]
#-----[ definice variacni ulohy - leva strana variacni ulohy ]-----
a1=-ufl.inner(sigma(u_trial,phi_trial),eps(u_test))*ufl.dx
a2=ufl.inner(D(u_trial,phi_trial),E(phi_test))*ufl.dx
a=a1+a2
#-----[ definice variacni ulohy - prava strana variacnich uloh ]-----
l1=ufl.dot(f,u_test)*ufl.dx-ufl.dot(T,u_test)*ds3
l=l1
#----[ reseni variacni ulohy ]-----
solution=LinearProblem(a=a,L=l,bcs=bcs,
petsc_options_prefix="piezo_",
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh=solution.solve()
u,phi=uh.sub(0).collapse(),uh.sub(1).collapse()
#-----[ ulozeni pro ParaView ]-----
with io.VTKFile(MPI.COMM_WORLD,"piezo solution/u.pvd","w") as vtk:
vtk.write_function(u,0.0) # 0.0 = time step
with io.VTKFile(MPI.COMM_WORLD,"piezo solution/phi.pvd","w") as vtk:
vtk.write_function(phi,0.0) # 0.0 = time step
#-----[ ]-----
#-----[ konec ]-----
#-----[ ]-----