Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

NGSolve coding issue

  • Dori
  • New Member
  • New Member
More
1 year 3 weeks ago - 1 year 2 weeks 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 
Last edit: 1 year 2 weeks ago by Dori.
More
1 year 2 weeks 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
 
  • Dori
  • New Member
  • New Member
More
1 year 1 week 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
 
Time to create page: 0.105 seconds