Last active
May 5, 2021 04:07
-
-
Save BinWang0213/3b32419ec7da69778349d4841f0795be to your computer and use it in GitHub Desktop.
ngsolve HDG solver
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
from ngsolve import * | |
from netgen.geom2d import SplineGeometry | |
from netgen.read_gmsh import ReadGmsh | |
from ngsolve.webgui import Draw | |
#Gmsh read support | |
#gmesh = ReadGmsh('gmsh_msh.msh') | |
#mesh = Mesh( gmesh ); | |
geo = SplineGeometry() | |
geo.AddRectangle((0, 0), (2, 0.41), bcs = ("default", "outlet", "default", "inlet")) | |
geo.AddCircle((0.2, 0.2), r=0.05, leftdomain = 0, rightdomain = 1, bc = "default") | |
mesh = Mesh( geo.GenerateMesh(maxh = 0.08) ); | |
#order = 3 | |
#mesh.Curve(order) | |
#Draw(mesh) | |
print("ElementBoundary=",mesh.GetBoundaries()) | |
print("NumEles=",mesh.ne) | |
print("NumEdges=",mesh.nedge) | |
print("NumVertex=",mesh.nv) | |
#PkPk-1 scheme | |
#scheme = 'PkP0' | |
scheme = 'PkPk-1' | |
order=2 | |
V1 = HDiv ( mesh, order = order, dirichlet = "default") | |
V2 = FESpace ( "vectorfacet", mesh, order = order, dirichlet = "default" ) | |
Q = L2(mesh, order = order-1, lowest_order_wb = True) | |
V = FESpace ([V1,V2,Q]) | |
#PkP0 scheme with the minimum DOF | |
if(scheme == 'PkP0'): | |
V1 = HDiv ( mesh, order = order, dirichlet = "default", hodivfree = True ) | |
V2 = FESpace ( "vectorfacet", mesh, order = order, dirichlet = "default" ) | |
Q = L2(mesh, order = 0, lowest_order_wb = True) | |
V = FESpace ([V1,V2,Q]) | |
#lowest_order_wb = following code | |
#hdiv dofs + TangentialFacet dofs | |
# 2d - version | |
#offset = (order +1) * (order + 2) + 3 * (order +1) | |
#for el in mesh.Elements(): | |
# dofs = V.GetDofNrs(el) | |
# V.SetCouplingType(dofs[offset] , COUPLING_TYPE.INTERFACE_DOF ) | |
#Change pressure as interface dof, we can not eliminate pressure out of the linear system | |
#only the local dofs from HDiv will be eliminated | |
#V.SetCouplingType(V.Range(2), COUPLING_TYPE.INTERFACE_DOF) | |
V.FinalizeUpdate() | |
u, uhat, p = V.TrialFunction() | |
v, vhat, q = V.TestFunction() | |
nu = 0.001 | |
alpha = 4 | |
condense=True | |
gradu = CoefficientFunction ( (grad(u),), dims=(mesh.dim,mesh.dim) ) | |
gradv = CoefficientFunction ( (grad(v),), dims=(mesh.dim,mesh.dim) ) | |
n = specialcf.normal(mesh.dim) | |
h = specialcf.mesh_size | |
def tang(vec): | |
return vec - (vec*n)*n | |
a = BilinearForm ( V, symmetric=True, condense=condense) | |
a += SymbolicBFI ( nu*InnerProduct ( gradu, gradv) ) | |
a += SymbolicBFI ( nu*InnerProduct ( gradu.trans * n, tang(vhat-v) ), element_boundary=True ) | |
a += SymbolicBFI ( nu*InnerProduct ( gradv.trans * n, tang(uhat-u) ), element_boundary=True ) | |
a += SymbolicBFI ( nu*alpha*order*order/h * InnerProduct ( tang(vhat-v), tang(uhat-u) ), element_boundary=True ) | |
a += SymbolicBFI ( -div(u)*q -div(v)*p ) | |
a.Assemble() | |
f = LinearForm ( V ) | |
#Pressure Neumann BCs | |
p_in = CoefficientFunction( 1.0 ) | |
f += - InnerProduct(p_in*n,v.Trace())*ds(definedon=mesh.Boundaries("inlet")) | |
p_out = CoefficientFunction( 0.2 ) | |
f += - InnerProduct(p_out*n,v.Trace())*ds(definedon=mesh.Boundaries("outlet")) | |
f.Assemble() | |
#Check static condensation | |
ndof_total = V.ndof | |
ndof_local=0 | |
for i in range(V.ndof): | |
if(V.CouplingType(i)==COUPLING_TYPE.LOCAL_DOF): | |
ndof_local+=1 | |
print("V1.ndof =", V1.ndof,V.Range(0)) | |
print("V2.ndof =", V2.ndof,V.Range(1)) | |
print("Q.ndof =", Q.ndof,V.Range(2)) | |
print(f"Local DOFs: {ndof_local}, Global DOFs: {ndof_total-ndof_local}, " | |
f"CompressRatio: {100-ndof_local/ndof_total*100:.2f}%") | |
u = GridFunction(V) | |
#we don't need to explict setup no slip BC here | |
#noslip = CoefficientFunction( (0,0) ) | |
#gfu.components[0].Set(noslip, definedon=mesh.Boundaries("default")) | |
r = f.vec.CreateVector() | |
#Dirichlet BC setup | |
r.data = f.vec - a.mat * u.vec | |
#Condense RHS | |
if condense: | |
r.data += a.harmonic_extension_trans * r | |
#Solve a condensed system | |
u.vec.data += a.mat.Inverse(V.FreeDofs(condense)) * r | |
#Recover the full solution | |
if condense: | |
u.vec.data += a.harmonic_extension * u.vec | |
u.vec.data += a.inner_solve * r | |
vel = u.components[0] | |
pres = u.components[2] | |
#gives a list of the integral on each part of the bnd of u*n | |
flux_bd = Integrate(vel * n, mesh, BND, region_wise = True) | |
Markers = list(set(mesh.GetBoundaries())) | |
flux = {m:0.0 for m in Markers} | |
for i,bd in enumerate(mesh.GetBoundaries()): | |
flux[bd] += flux_bd[i] | |
print(flux) | |
Draw( Norm(vel), mesh, "absvel(hdiv)") | |
Draw( Norm(pres), mesh, "pressure") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#Cell-wise interpolation (continue from hdg.py) | |
vel = u.components[0] | |
ele=ElementId(VOL, 118) | |
#Get element pts | |
Pts=[] | |
for Vi in mesh[ele].vertices: | |
Pts+=[np.array(mesh[Vi].point)] | |
Pts=np.array(Pts) | |
#Get the centroid of cell 118 | |
cntr= np.mean(Pts,axis=0) #physical coordinate | |
cntr_ngs = mesh(*cntr) | |
cntr_local = cntr_ngs.pnt #local coordinate | |
#----------NGS Velocity interpolation ------ | |
vel = u.components[0] | |
print('Ngs Interp=', vel(cntr_ngs) ) | |
#----------Manully Velocity interpolation ------ | |
#Get shape function of BDMk | |
vel = u.components[0] | |
FE=vel.space.GetFE(ele) | |
shapeFunc=FE.CalcShape(cntr_local[0],cntr_local[1]).NumPy() | |
#Get solution coeff for BDMk | |
uIdx=np.array(vel.space.GetDofNrs(ele)) | |
print(uIdx) | |
coeff = u.vec.FV().NumPy()[uIdx] | |
#Get piola trans (not sure how to get this in ngsolve) | |
Jac=np.array( [[ Pts[1,0]-Pts[0,0],Pts[2,0]-Pts[0,0] ], #[(x2-x1, x3-x1), | |
[ Pts[1,1]-Pts[0,1],Pts[2,1]-Pts[0,1]] ] ) # (y2-y1, y3-y1)] | |
trans=Jac/np.linalg.det(Jac) | |
print(Jac) | |
#Interpolate velocity by shapefunc*coeff*trans_mat | |
shapeFunc_trans = np.dot(trans,shapeFunc.T) | |
vel_interp = np.dot(shapeFunc_trans, coeff) | |
print('Manual Interp=',vel_interp) | |
#Check if two interpolation has the same value | |
assert np.all( np.isclose(vel_interp, vel(cntr_ngs)) ) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
TODO list:
The codes runs smoothly, but the results are not correct.
In fenics we can do
flux_bd = Integrate(vel * n, mesh, BND, region_wise = True)
Not sure how to set a piecewise material coefficient. Namely, we need to set material properties for each cell.
In fenics, we can use DG0 and set value for each cell below.
I can not get the same velocity from my manually velocity interpoaltion
Ngsolve seems to use different trans matrix other than piola transform
Installing ngsolve on cluster now.. Will try this later on