jorgensd / jorgensd/MultiMeshShapeOpt_code
MultiMesh in solving eigenvalue problem of Stokes operator
- 主要言語
- Python
- スター
- 3
- フォーク
- 1
- PR マージ指標
- 30日以内にマージされた PR はありません
説明
I am following your code to use multimesh to solve eigenvalue problem of the Stokes operator.
By selecting test function from a rawer mesh, I want to weaken the divergence-free condition.
When ``mesh_u`` and ``mesh_div_free`` are the same mesh, the eigenfunction should be completely divergence-free.
But the code below failed to apply the divergence-free condition to solution ``u``.
Could you give me any hint where the code is wrong?
```
from fenics import *
mesh_u=UnitSquareMesh(32, 32, diagonal="right")
mesh_div_free=UnitSquareMesh(16, 16, diagonal="right")
mesh_R=UnitSquareMesh(1, 1, diagonal="right")
multimesh = MultiMesh()
multimesh.add(mesh_u)
multimesh.add(mesh_div_free)
multimesh.add(mesh_R)
multimesh.build()
degree=3
CG = VectorElement("CG", triangle, degree,2)
DG = FiniteElement("DG", triangle, degree-1)
Real = FiniteElement("Real", triangle, 0)
MixedV = MultiMeshFunctionSpace(multimesh, MixedElement([CG, DG, Real]))
V = MultiMeshSubSpace(MixedV, 0)
(u, p, alpha) = TrialFunctions(MixedV)
(v, q, beta) = TestFunctions(MixedV)
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx + alpha*q*dx + beta*p*dx
b = inner(u,v)*dx
f = Constant((0.0, 0.0))
L = inner(f,v)*dx
mvc = MeshValueCollection("size_t", mesh_u, 1)
mf = cpp.mesh.MeshFunctionSizet(mesh_u, mvc)
noslip = Constant((0.0, 0.0))
marker=100
bc0=MultiMeshDirichletBC(V, noslip, mf, marker, 0)
A = assemble_multimesh(a)
B = assemble_multimesh(b)
bc0.apply(A)
bc0.zero(B)
# downcast to PETSc matrices
MatA = as_backend_type(A)
MatB = as_backend_type(B)
import scipy.io as myio
import scipy.sparse as sp
row_, col_, val_ = MatA.mat().getValuesCSR()
mat_A = sp.csr_matrix((val_,col_,row_))
row_,col_,val_ = MatB.mat().getValuesCSR()
mat_B = sp.csr_matrix((val_,col_,row_))
#Another approach to remove boundary value DOF
bd_int = inner(u,v)*ds
bd_int = assemble_multimesh(bd_int)
Mat_bd_int = as_backend_type(bd_int)
row_,col_,val_ = Mat_bd_int.mat().getValuesCSR()
Mat_bd_int = sp.csr_matrix((val_,col_,row_))
myio.savemat("./matrix/matrix.mat",{'A':mat_A,'B':mat_B,'BD_Int':Mat_bd_int})
```
(I exported the matrix to mat file and then solve it in Octave or MATLAB.)
コントリビューションガイド
このリポジトリのコントリビューションガイドは索引されていません
評価
この issue はまだ評価されていません。