Re: [bempp] Construct Mass Matrix

[Matthew Scroggs]

Dear Sihao,

You can construct an identity operator using:

id = bempp.api.operators.boundary.sparse.identity(domain_space, range_space, dual_space)

You can then find its weak form (ie the mass matrix) with:

mass = id.weak_form()

Hope this helps,

Matthew

On 23/11/2019 01:34, wang...@gmail.com wrote:

Hi

I need to construct mass matrices to do my calculation, is there any fast way to construct mass matrices in BEM++? Thanks in advance.

Sihao

--
You received this message because you are subscribed to the Google Groups "bempp" group.
To unsubscribe from this group and stop receiving emails from it, send an email to bempp+un...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/bempp/c0deecaa-a639-451e-a34b-c0d2e059ac10%40googlegroups.com.

[wang...@gmail.com]

Hi Scroggs,

Thank you for your reply, I wanted to transform the operator to weak form to solve my problem, but the solution using the following two ways are different.(here we compare the x) Can you tell me what is the difference?  I tried the weak form of slp for the indirect method and the solutions are same. Thanks.

One is

rhs = (-.5 * identity + dlp) * dirichlet_fun

neumann_fun, info = bempp.api.linalg.gmres(slp, rhs, tol=1E-8)

x = neumann_fun.projections(piecewise_const_space)

Another one is  

matrix1 = slp.weak_form()

matrix2 = dlp.weak_form()

matrix3 = identity.weak_form()

cmatrix = -0.5*matrix3 + matrix2

dirichlet_data = dirichlet_fun.projections(piecewise_const_space)

rhs = (cmatrix) * dirichlet_data

#print rhs

x, info = gmres(matrix1, rhs, tol=1E-8)

import bempp.api

import math

import numpy as np

from numpy import linalg as LA

from scipy.sparse.linalg import gmres

from numpy import linalg as LA

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

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

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

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

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)

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

fourpi = (1.0/4.0/np.pi)**1.0

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

    result[0] = fourpi / (rtwo**0.5)

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

rhs = (-.5 * identity + dlp) * dirichlet_fun

neumann_fun, info = bempp.api.linalg.gmres(slp, rhs, tol=1E-8)

x = neumann_fun.projections(piecewise_const_space)

'''

matrix1 = slp.weak_form()

matrix2 = dlp.weak_form()

matrix3 = identity.weak_form()

cmatrix = -0.5*matrix3 + matrix2

dirichlet_data = dirichlet_fun.projections(piecewise_const_space)

rhs = (cmatrix) * dirichlet_data

#print rhs

x, info = gmres(matrix1, rhs, tol=1E-8)

'''

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

fourpi = (1.0/4.0/np.pi)**1.0

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

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

result[0] = -fourpi / (rtwo**1.5) *ip

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

sol = f3.projections(piecewise_const_space)

diff =  LA.norm(x-sol, np.inf)/ LA.norm(sol, np.inf)

print diff

在 2019年11月25日星期一 UTC-6上午3:26:42，Matthew Scroggs写道：

Dear Sihao,

You can construct an identity operator using:

id = bempp.api.operators.boundary.sparse.identity(domain_space, range_space, dual_space)

You can then find its weak form (ie the mass matrix) with:

mass = id.weak_form()

Hope this helps,

Matthew

On 23/11/2019 01:34, wang...@gmail.com wrote:

Hi

I need to construct mass matrices to do my calculation, is there any fast way to construct mass matrices in BEM++? Thanks in advance.

Sihao

--
You received this message because you are subscribed to the Google Groups "bempp" group.

To unsubscribe from this group and stop receiving emails from it, send an email to be...@googlegroups.com.

[Matthew Scroggs]

Hi Sihao,

In the first example, `x = neumann_fun.projections(piecewise_const_space)` should be changed to `x = neumann_fun.coefficients`, and similarly `dirichlet_data` should be changed in the second code. The codes should then give the same result.

`projections` will apply an inverse mass matrix to the coefficients of the grid function, so gives a different result to `coefficients`.

Hope this helps,

Matthew

[Sihao Wang]

Hi Matthew,

Thanks for your reply. I have a quite relative question 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 a cube. Do you know where is the problem? 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

在 2019年11月28日星期四 UTC-6上午7:14:29，Matthew Scroggs写道：

[Matthew Scroggs]

Hello,

Is your solution better is you use a more refined cube, for example:

`grid = bempp.api.shapes.cube(h=0.05)`

All the best,

Matthew