- Thank you received: 0
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.
Interface domain is not working.
4 years 4 months ago #2245
by Amad
Interface domain is not working. was created by Amad
Hello everyone,
I am trying to solve an Acoustic FSI problem with the HDG method, but an error has occurred. The following lines are the (global) code I wrote.
Creating geometry:
Spaces generated:
bilinear form:
Then, I got the following error:
NgException: FacetBilinearFormIntegrator::CalcFacetMatrix<Complex> for inner facets not implemented!in Assemble BilinearForm 'biform_from_py'
1- Assuming I wrote the correct formulation, why did I get this error?
2- Am I skipping any steps? Do I need to include anything else?
Finally, I tried to solve the same problem with Galerkin FEM and got the same error.
Please let me know what I did wrong.
Thanks in advance and have a nice new year!
Best Regards,
Alan
I am trying to solve an Acoustic FSI problem with the HDG method, but an error has occurred. The following lines are the (global) code I wrote.
Creating geometry:
Code:
geometry = SplineGeometry()
# point coordinates ...
pnts = [ (-0.5,-0.5), (0,-0.5), (0,0.5), (-0.5,0.5), (0.5,-0.5), (0.5,0.5)]
pnums = [geometry.AppendPoint(*p) for p in pnts]
# start-point, end-point, boundary-condition, domain on left side, domain on right side:
geometry.Append(["line", pnums[0], pnums[1]], bc="fluid", leftdomain=1, rightdomain=0)
geometry.Append(["line", pnums[1], pnums[2]], bc="interface", leftdomain=1, rightdomain=2)
geometry.Append(["line", pnums[2], pnums[3]], bc="fluid", leftdomain=1, rightdomain=0)
geometry.Append(["line", pnums[3], pnums[0]], bc="fluid", leftdomain=1, rightdomain=0)
geometry.Append(["line", pnums[1], pnums[4]], bc="solid", leftdomain=2, rightdomain=0)
geometry.Append(["line", pnums[4], pnums[5]], bc="solid", leftdomain=2, rightdomain=0)
geometry.Append(["line", pnums[5], pnums[2]], bc="solid", leftdomain=2, rightdomain=0)
geometry.SetMaterial(1,"fluid")
geometry.SetMaterial(2,"solid")
mesh = Mesh(geometry.GenerateMesh(maxh=0.1))
Spaces generated:
Code:
U_s = VectorL2(mesh, order=order, complex=True, definedon="solid")
F_s = FacetFESpace(mesh, order=order, complex=True, definedon="solid")
U_f = VectorL2(mesh, order=order, complex=True, definedon="fluid")
F_f = FacetFESpace(mesh, order=order, complex=True, definedon="fluid")
Q = L2(mesh, order=order, complex=True) # Pressure space
fes = FESpace([Q, U_f, U_s, F_f, F_f, F_s, F_s], highest_order_dc=True)
(p_f, u_f, u_s, uhatx_f, uhaty_f, uhatx_s, uhaty_s), \
(q_f, v_f, v_s, vhatx_f, vhaty_f, vhatx_s, vhaty_s) = fes.TnT()
bilinear form:
Code:
a = BilinearForm(fes, condense=condense)
# fluid domain
dS = dx(element_boundary=True, definedon=mesh.Materials("fluid"))
a += (.........) * dx("fluid") # volume elements
a += (.........) * dS # boundary elements
a += (.........) * ds(definedon=mesh.Boundaries("fluid")) # fluid boundary domain
# solid domain
dS = dx(element_boundary=True, definedon=mesh.Materials("solid"))
a += (.........) * dx("solid") # volume elements
a += (.........) * dS # boundary elements
a += (.........) * ds(definedon=mesh.Boundaries("solid")) # solid boundary domain
# solid domain
dS = dx(element_boundary=True, definedon=mesh.Materials("interface"))
a += (.........) * dx("solid") # volume elements
a += (.........) * dS # boundary elements
a += (.........) * ds(definedon=mesh.Boundaries("solid")) # solid boundary domain
# interface fluid-solid domain
a += (.........) * ds(definedon=mesh.Boundaries("interface"))
# source
f = LinearForm(fes)
f += (.........) * ds(definedon=mesh.Boundaries("fluid")) # fluid
f += (.........) * ds(definedon=mesh.Boundaries("solid")) # solid
a.Assemble()
f.Assemble()
Then, I got the following error:
NgException: FacetBilinearFormIntegrator::CalcFacetMatrix<Complex> for inner facets not implemented!in Assemble BilinearForm 'biform_from_py'
1- Assuming I wrote the correct formulation, why did I get this error?
2- Am I skipping any steps? Do I need to include anything else?
Finally, I tried to solve the same problem with Galerkin FEM and got the same error.
Please let me know what I did wrong.
Thanks in advance and have a nice new year!
Best Regards,
Alan
4 years 3 months ago #2246
by joachim
Replied by joachim on topic Interface domain is not working.
do your integrals involve u.Other() terms ?
You should not need them for HDG, but they are certainly the cause for the appearance of the FacetBilinearFormIntegrator (which complains about the complex form).
Best, Joachim
You should not need them for HDG, but they are certainly the cause for the appearance of the FacetBilinearFormIntegrator (which complains about the complex form).
Best, Joachim
4 years 3 months ago #2248
by Amad
Replied by Amad on topic Interface domain is not working.
Dear Joachim,
Thank you for the prompt reply.
No, I am not using u.Others().
I am creating these trace functions:
After that, I am using these functions to integrate into the interface, for example. I have attached a pdf of the equation I am trying to code using NGSolve. For example, in the interface domain, I used equation (32) in (30) and got the equation as follows
Is this correct? Am I missing something?
Thank you again for the help.
Best Regards,
Alan
Thank you for the prompt reply.
No, I am not using u.Others().
I am creating these trace functions:
Code:
uhat_s = CoefficientFunction( (uhatx_s, uhaty_s) )
vhat_s = CoefficientFunction( (vhatx_s, vhaty_s) )
uhat_tr_s = CoefficientFunction( (uhatx_s.Trace(), uhaty_s.Trace()) )
vhat_tr_s = CoefficientFunction( (vhatx_s.Trace(), vhaty_s.Trace()) )
uhat_f = CoefficientFunction( (uhatx_f, uhaty_f) )
vhat_f = CoefficientFunction( (vhatx_f, vhaty_f) )
uhat_tr_f = CoefficientFunction( (uhatx_f.Trace(), uhaty_f.Trace()) )
vhat_tr_f = CoefficientFunction( (vhatx_f.Trace(), vhaty_f.Trace()) )
After that, I am using these functions to integrate into the interface, for example. I have attached a pdf of the equation I am trying to code using NGSolve. For example, in the interface domain, I used equation (32) in (30) and got the equation as follows
Code:
interfaceTerm = p_f*(v_f - vhat_tr_f)*n + q_f*(u_f - rho*omega_s**2*uhat_tr_s)*n + \
beta_f*(u_f - rho*omega_s**2*uhat_tr_s)*n*(v_f - vhat_tr_f)*n
interface_int = interfaceTerm *ds(definedon=mesh.Boundaries("interface"))
a += interface_int
Is this correct? Am I missing something?
Thank you again for the help.
Best Regards,
Alan
Attachments:
4 years 3 months ago #2249
by joachim
Replied by joachim on topic Interface domain is not working.
can you attach a complete example producing the exception ?
4 years 3 months ago #2251
by Amad
Replied by Amad on topic Interface domain is not working.
Dear Joachim,
I have attached code that is not working. The code works very well if we choose only one domain, i.e., only Helmholtz or elastic wave problem. But when trying to solve the coupled problem, no. So I think I'm doing something wrong on the interface.
Also, I just realised that the exception occurs when I'm trying to identify the interface that was developed in this post ngsolve.org/forum/ngspy-forum/114-interface-integrals#491
Thank you again.
I have attached code that is not working. The code works very well if we choose only one domain, i.e., only Helmholtz or elastic wave problem. But when trying to solve the coupled problem, no. So I think I'm doing something wrong on the interface.
Also, I just realised that the exception occurs when I'm trying to identify the interface that was developed in this post ngsolve.org/forum/ngspy-forum/114-interface-integrals#491
Code:
def GetInterfaceIndicator(mesh):
dom_indicator = GridFunction(L2(mesh,order=0))
dom_indicator.Set(CoefficientFunction([1,2]))
ffes = FacetFESpace(mesh,order=0)
dom_ind_average = GridFunction(ffes)
mf = BilinearForm(ffes,symmetric=True)
mf += SymbolicBFI(ffes.TrialFunction()*ffes.TestFunction(),element_boundary=True)
mf.Assemble()
ff = LinearForm(ffes)
ff += SymbolicLFI(dom_indicator*ffes.TestFunction(),element_boundary=True)
ff.Assemble()
dom_ind_average.vec.data = mf.mat.Inverse() * ff.vec
interface_indicator = BitArray(ffes.ndof)
for i in range(len(dom_ind_average.vec)):
if (dom_ind_average.vec[i] > 1.75) or (dom_ind_average.vec[i] < 1.25):
interface_indicator[i] = False
else:
interface_indicator[i] = True
print(sum(interface_indicator))
return interface_indicator
interface_indicator = GetInterfaceIndicator(mesh)
interfacebfi = SymbolicBFI (interfaceTerm, VOL, skeleton=True)
interfacebfi.SetDefinedOnElements(interface_indicator)
a += interfacebfi
Thank you again.
Attachments:
4 years 3 months ago #2252
by joachim
Replied by joachim on topic Interface domain is not working.
You don't use the pressure in the solid domain. Add a definedon, and your code runs through:
Otherwise, you get error messages like
ZGETRF::info = 1
ZGETRI::info = 1
which stem from Lapack complaining about singular matrices in local condensation.
Yes, as you have observed, the FacetBilinearFormIntegrator::CalcFacetMatrix<Complex> exception comes from the "skeleton=True" complex integrator (which you should not need).
Joachim
Code:
Q = L2(mesh, order=order, complex=True , definedon="fluid")
Otherwise, you get error messages like
ZGETRF::info = 1
ZGETRI::info = 1
which stem from Lapack complaining about singular matrices in local condensation.
Yes, as you have observed, the FacetBilinearFormIntegrator::CalcFacetMatrix<Complex> exception comes from the "skeleton=True" complex integrator (which you should not need).
Joachim
Time to create page: 0.160 seconds
