Interface domain is not working.

Main Forum Old-Forum - do not post here
1

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:

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:

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:

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:
[color=red]
NgException: FacetBilinearFormIntegrator::CalcFacetMatrix for inner facets not implemented!in Assemble BilinearForm ‘biform_from_py’
[/color]

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

2

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

3

Dear Joachim,

Thank you for the prompt reply.

No, I am not using u.Others().

I am creating these trace functions:

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,
https://ngsolve.org/media/kunena/attachments/1211/AcousticFSI.pdfI 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

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

Attachment: AcousticFSI.pdf

4

can you attach a complete example producing the exception ?

5

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 interface integrals - Kunena

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

https://ngsolve.org/media/kunena/attachments/1211/AcousticFSI.py

Thank you again.

Attachment: AcousticFSI_2019-12-30.py

6

You don’t use the pressure in the solid domain. Add a definedon, and your code runs through:

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 exception comes from the “skeleton=True” complex integrator (which you should not need).

Joachim

7

Dear Joachim,

Thank you for the answer. Now it looks like the interface domain is working. The numerical result is not correct yet, I may have made a mistake in the formulation.
I am checking the formulation again and I hope to figure it out soon.
Thank you for your help. If I don’t find the error, I send a new message.

Best Regards and happy new year,

Alan