NGSolve coding issue
- Dori
- Topic Author
- New Member
Less
More
1 year 8 months ago - 1 year 8 months ago #4740
by Dori
NGSolve coding issue was created by Dori
Hello,
I am trying to code a problem that involves two domains connected by an interface. There is pressure in the domains and at the interface. Therefore, I have a 2D and a 1D space. I tried to develop these pieces of code, but I am getting an anumerical factorization error( The matrix is singular). Any help would be greatly appreciated.
V = VectorL2(mesh, order=order)
Qi = H1(mesh, order=order, orderinner=0, definedon=mesh.Boundaries("gammai"), dirichlet="gamma1|gamma2")
Q1 = L2(mesh,order=order,definedon=mesh.Materials("Domain1"))
Q2 = L2(mesh,order=order,definedon=mesh.Materials("Domain2"))
Qbar1 = FacetFESpace(mesh, order=order,definedon=mesh.Materials("Domain1"), dirichlet="gamma1")
Qbar2 = FacetFESpace(mesh, order=order, definedon=mesh.Materials("Domain2"), dirichlet="gamma2")
X=FESpace([V,Q1,Qbar1,Q2,Qbar2,Qi])
(u,p1,pbar1,p2,pbar2,p_i),(v,q1,qbar1,q2,qbar2,qi)=X.TnT()
dx1 = dx(definedon=mesh.Materials("Domain1"))
ds1 = dx(element_boundary = True, definedon=mesh.Materials("Domain1"))
dx2 = dx(definedon=mesh.Materials("Domain2"))
ds2 = dx(element_boundary = True, definedon=mesh.Materials("Domain2"))
dsi = ds(skeleton = True, definedon=mesh.Boundaries("gammai"))
a=BilinearForm(X,check_unused=False)
a += u*v*dx
a += p1*q1*dx1
a += p2*q2*dx2
a += pbar1*qbar1*ds1
a += pbar2*qbar2*ds2
a += p_i*qi*dsi
f=LinearForm(X)
gfu=GridFunction(X)
a.Assemble()
f.Assemble()
res = f.vec.CreateVector()
res.data = f.vec - a.mat * gfu.vec
gfu.vec.data=a.mat.Inverse(X.FreeDofs(),inverse="umfpack")*res
I am trying to code a problem that involves two domains connected by an interface. There is pressure in the domains and at the interface. Therefore, I have a 2D and a 1D space. I tried to develop these pieces of code, but I am getting an anumerical factorization error( The matrix is singular). Any help would be greatly appreciated.
V = VectorL2(mesh, order=order)
Qi = H1(mesh, order=order, orderinner=0, definedon=mesh.Boundaries("gammai"), dirichlet="gamma1|gamma2")
Q1 = L2(mesh,order=order,definedon=mesh.Materials("Domain1"))
Q2 = L2(mesh,order=order,definedon=mesh.Materials("Domain2"))
Qbar1 = FacetFESpace(mesh, order=order,definedon=mesh.Materials("Domain1"), dirichlet="gamma1")
Qbar2 = FacetFESpace(mesh, order=order, definedon=mesh.Materials("Domain2"), dirichlet="gamma2")
X=FESpace([V,Q1,Qbar1,Q2,Qbar2,Qi])
(u,p1,pbar1,p2,pbar2,p_i),(v,q1,qbar1,q2,qbar2,qi)=X.TnT()
dx1 = dx(definedon=mesh.Materials("Domain1"))
ds1 = dx(element_boundary = True, definedon=mesh.Materials("Domain1"))
dx2 = dx(definedon=mesh.Materials("Domain2"))
ds2 = dx(element_boundary = True, definedon=mesh.Materials("Domain2"))
dsi = ds(skeleton = True, definedon=mesh.Boundaries("gammai"))
a=BilinearForm(X,check_unused=False)
a += u*v*dx
a += p1*q1*dx1
a += p2*q2*dx2
a += pbar1*qbar1*ds1
a += pbar2*qbar2*ds2
a += p_i*qi*dsi
f=LinearForm(X)
gfu=GridFunction(X)
a.Assemble()
f.Assemble()
res = f.vec.CreateVector()
res.data = f.vec - a.mat * gfu.vec
gfu.vec.data=a.mat.Inverse(X.FreeDofs(),inverse="umfpack")*res
Last edit: 1 year 8 months ago by Dori.
1 year 8 months ago #4747
by joachim
Replied by joachim on topic NGSolve coding issue
Hi Dori,
if your mesh is not too fine, you can check the stiffness matrix by hand.
The pi*qi*dsi integral does not need the skeleton=True argument - it needs only the pressure at the interface.
Best,
Joachim
if your mesh is not too fine, you can check the stiffness matrix by hand.
The pi*qi*dsi integral does not need the skeleton=True argument - it needs only the pressure at the interface.
Best,
Joachim
- Dori
- Topic Author
- New Member
Less
More
1 year 8 months ago #4764
by Dori
Replied by Dori on topic NGSolve coding issue
I am grateful for your assistance. However, I am experiencing another problem. Within my model, there is an equation in the interface that should be calculated in a linear form, but the code does not seem to be executing this correctly. No matter what coefficient function I apply to this right-hand side, the output remains consistent.
Any assistance with resolving this issue would be greatly appreciated.
Best,
Dori
Any assistance with resolving this issue would be greatly appreciated.
Best,
Dori
Time to create page: 0.097 seconds