Dear NGSolve developers;

I am trying to understand how NGSolve solves the problems with boundary conditions defined on surface that is symmetric about a coordinate plane. When solving

[tex]-\Delta u -k^2 u = 0 \quad {\rm{\in}} \quad \Omega \

\frac{\partial u}{\partial n} = 1.0 \quad {\rm{on}} \quad \Gamma_1\

\frac{\partial u}{\partial n} -\beta u = 0 \quad {\rm{on}} \quad \Gamma_2[/tex]

with \Gamma_1 (symmetric_srf) and \Gamma_2 (outer) nonintersecting, \Omega a sphere of given dimensions, and \Gamma_1 is symmetric about one of the coordinate planes. Taking advantage of the symmetry of \Gamma_1 when the problem is divided in half, I see the solution is scaled by a factor of 2 compared to the original problem. This is a bit surprising because in the original problem the pressure field is symmetric, i.e., du/dn = 0 on the symmetry plane. Hence I should have been able to take any subset of the full problem, and apply either Dirichlet or Neumann boundary conditions, taken from the full problem, at the boundary of the subset, and be able to compute the solution correctly, using the boundary conditions and sources internal to the subset. Note that when the full model is sliced into half to take advantage of the symmetry, i don’t specify any condition on the newly created surface (thus a natural bc is imposed there which i believe is already from the subset of the original problem). Could you comment on what i could be doing wrong here?

Below the relevant snippet

```
mesh = Mesh(ngmesh)
fes = H1(mesh, complex=True, order=2)
u = fes.TrialFunction()
v = fes.TestFunction()
g = 1;
beta = 1j*waveno
a = BilinearForm (fes, symmetric=True)
a += SymbolicBFI (grad(u)*grad(v) )
a += SymbolicBFI (-waveno*waveno*u*v)
a += SymbolicBFI(beta*u*v, definedon=mesh.Boundaries("outer"))
f = LinearForm (fes)
f += v * g * ds("symmetric_srf")
gfu = GridFunction (fes)
a.Assemble()
f.Assemble()
```

At first glance I do not see any problems with that. Could you post a minimal example?

Best

Hi Christopher,

MWE:

I tried to recreate the pressure doubling issue with a different geometry -

half and full model for sphere (R1=1) in sphere (R2=2) where the interior sphere surface is excited with neumann bc. I tried both an impedance and PML condition. I don’t see pressure doubling in the MWE which is reassuring.

Real Geometry where pressure doubling happens:

Unfortunately, I can’t post the original geometry STEP file here but the

GMSH mesh for the actual geometry is attached - realGeomHalfmodel.txt.

The attached code.py is practically what i used for the real geometry problem (except the bit where python calls gmsh and meshes using .geo file) where i get the pressure doubling. So i am guessing i must be doing something wrong when i create the mesh via the .geo file. Plane 6 is symmetrical about the vertical plane passing through it and is excited with neumann bc. In the full model, plane 6 would have been a full circle and the model would have contained the elements from the other half.

Attachment: code_2020-05-20-2.py

Attachment: fullmodel_2020-05-20-2.geo

Attachment: makemesh_2020-05-20-2.py

Attachment: halfmodel_2020-05-20-2.geo

Attachment: realGeomHalfmodel.txt

If you can send me the step file at clackner@cerbsim.com ? The attached mesh seems to have no surface elements.

Is there a specific reason you use gmsh for the meshing? If you use netgen you can have higher order curved boundaries which would be better for absorbing bc on a sphere…

My first guess is that something in the mesh is not correct that may cause the factor.

Best

Your observation is correct:

If you apply the source at the symmetry plane of the full model, the weak form leads to the interface condition

[du/dn] = g

i.e.

du/dn_left + du/dn_right = g,

i.e. by symmetry

du/dn_left = du/dn_right = g/2

while for the half problem the source term is

du/dn = g

in physical words:

if you apply the source term for the full problem, half of the energy is going to the left, and half to the right.

if you do the half problem, the full energy is going to this side.

Joachim

The reason for using gmsh is ease with which geo files can be edited and also ability to assign various physical tags to the surfaces/lines/volumes etc.

[quote=“joachim” post=2711]Your observation is correct:

If you apply the source at the symmetry plane of the full model, the weak form leads to the interface condition

[du/dn] = g

i.e.

du/dn_left + du/dn_right = g,

i.e. by symmetry

du/dn_left = du/dn_right = g/2

while for the half problem the source term is

du/dn = g

in physical words:

if you apply the source term for the full problem, half of the energy is going to the left, and half to the right.

if you do the half problem, the full energy is going to this side.

s

Joachim[/quote]

I am not sure if I understand your explanation: are you effectively saying then only half the source magnitude must be applied when doing load integrals in case of half problem? Perhaps it would work for a point source but not for distributed loading on the surface. In this case, either for a full or a half problem, the Neumann data g is integrated on a specific area. Wouldn’t ngsolve do the integral contributing to the load vector f as:

[tex]

f_i = \int_{\Gamma_1} g \phi_i ds

[/tex]

where \phi_i are the basis functions. So integral for the full problem will be performed on a full circle whereas that for the half problem is performed only on the half. So although the magnitude of the source term is the same in both cases, the half model only “sees” half its contribution due to the affected area being half. Thus no scaling should be required when dealing with half model (again assuming that the dp/dn=0 exists in the original problem and nothing needs to be said on the symmetrical plane) Does it make sense?