Problem in dealing with different geometries

[Sihao Wang]

Hi,

I have a problem for using bempp to solve green's integral equation. The code is attached. In my code, when the grid is a sphere, the result is decent. However, when I changed it to a cube, the solution is not accurate at all. Here I used the fundamental solution to test my result, so it should work for different geometries, I also tested it for ellipsoid and it works. It seems it only work for smooth geometries. Do you know where is the problem? How can I make it work for a cube? Thanks.

import bempp.api

import math

import numpy as np

from numpy import linalg as LA

from scipy import special

from scipy.sparse.linalg import gmres

#grid = bempp.api.shapes.sphere()

grid = bempp.api.shapes.cube()

piecewise_const_space = bempp.api.function_space(grid, "DP", 0)

piecewise_linear_space = bempp.api.function_space(grid, "P", 1)

slp = bempp.api.operators.boundary.laplace.single_layer(piecewise_const_space, piecewise_const_space, piecewise_const_space)

dlp = bempp.api.operators.boundary.laplace.double_layer(piecewise_const_space, piecewise_const_space, piecewise_const_space)

identity = bempp.api.operators.boundary.sparse.identity(piecewise_const_space, piecewise_const_space, piecewise_const_space)

matrix1 = slp.weak_form()

matrix2 = dlp.weak_form()

matrix3 = identity.weak_form()

def u_inc(x,n, domain_index, result):

    rtwo = (x[0])**2 + x[1]**2 + x[2]**2

    z = np.sqrt(rtwo)

    result[0] = 1/z

def u_inc2(x,n, domain_index, result):

    rtwo = (x[0])**2 + x[1]**2 + x[2]**2

    z = np.sqrt(rtwo)

    ip = (x[0])*n[0]+(x[1])*n[1]+(x[2])*n[2]

    result[0] = -ip/(z**3)

f = bempp.api.GridFunction(piecewise_const_space, fun = u_inc2)

p = f.coefficients

result2 = p

k = matrix1 * result2

x,info = gmres(matrix2-0.5*matrix3,k, tol=1E-8)

result = x

f3 = bempp.api.GridFunction(piecewise_const_space, fun = u_inc)

sol = f3.coefficients

result3 = sol

diff =  LA.norm(x-sol)

diff2 =  LA.norm(x-sol)/ LA.norm(sol)

print diff

print diff2

[Matthew Scroggs]

You may need to refine your grid to get better results. Try replacing

grid = bempp.api.shapes.cube()

with

grid = bempp.api.shapes.cube(h=0.1)

The value you give to h will be the smallest edge length in the mesh that will be generated.

Hope this helps,

Matthew

[Sihao Wang]

Hi,

I tried it for grid refinement, but the solution has a big difference with the accurate solution. Thanks.

Best,

Sihao

在 2020年1月21日星期二 UTC-6上午10:23:49，Matthew Scroggs写道：

[Sihao Wang]

Hi,

Sorry I find where the issue is. The cube's center is not at (0,0,0), so I need to find another point inside the cube to make it correct. Thanks.

Best,

Sihao

在 2020年1月21日星期二 UTC-6上午10:23:49，Matthew Scroggs写道：

You may need to refine your grid to get better results. Try replacing