Problem Reconstructing GridFunction from a (P,1) space.

[Luiz Augusto Alvim]

Hello, I'm solving the neumann problem for exterior acoustic scattering in bempp-cl. Trouble arises when reconstructing the GridFunction from the solution coefficients saved to a numpy array. When using a DP,0 space, I can use my reconstructed GridFunction to evaluate my potential operators and it all works flawlessly. Nonetheless, when using a P,1 space, I get a "dimension mismatch" error:

###

  pScat =  dlp_pot.evaluate(boundData)

  File "C:\Users\gutoa\pypkg\bempp\bempp\api\assembly\potential_operator.py", line 27, in evaluate

    return self._evaluator.evaluate(grid_fun.coefficients)

  File "C:\Users\gutoa\pypkg\bempp\bempp\core\dense_potential_assembler.py", line 104, in evaluate

    self.space.dof_transformation @ coefficients

  File "C:\Users\gutoa\Anaconda3\envs\tcc\lib\site-packages\scipy\sparse\base.py", line 561, in __matmul__

    return self.__mul__(other)

  File "C:\Users\gutoa\Anaconda3\envs\tcc\lib\site-packages\scipy\sparse\base.py", line 499, in __mul__

    raise ValueError('dimension mismatch')

ValueError: dimension mismatch

###

Is there something I am missing in this process? I am using the appropriate "space" variable, the GridFunction coefficients seem to match up exactly before and after reconstruction. I can't figure out why I am getting this error.

Lastly, I would like to compliment Timo for the incredible work in the package, it is really nice to be able to use this tool openly.

Best Regards,

Luiz Augusto

[Matthew Scroggs]

Hi Luiz,

What spaces are you using the define dlp_pot and boundData? Without these it's hard to work out what's going on.

[mrban...@gmail.com]

Thanks for the quick response Matthew. My dlp_pot and space are defined by:

qrid = bempp.api.import_grid('my_mesh.msh')

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

dlp_pot = bempp.api.operators.potential.helmholtz.double_layer(space, pts, k)

My boundData comes from the following formulation.

            identity = bempp.api.operators.boundary.sparse.identity(

                self.space, self.space, self.space)

            dlp = bempp.api.operators.boundary.helmholtz.double_layer(

                self.space, self.space, self.space, k)

            hyp = bempp.api.operators.boundary.helmholtz.hypersingular(self.space, self.space, self.space, k)

                @bempp.api.complex_callable(jit=False)

                def combined_data(r, n, domain_index, result):

                    result[0]=0

                    for i in range(len(self.r0.T)): 

                        pos  = np.linalg.norm(r-self.r0[:,i].reshape(1,3),axis=1)

                        val  = self.q.flat[i]*np.exp(1j*k*pos)/(pos)

                        result[0] += val/(pos**2) * (1j*k*pos-1)* np.dot(r-self.r0[:,i],n) - (ni*val)

                

            monopole_fun = bempp.api.GridFunction(self.space, fun=combined_data)

            lhs =  -hyp + ni*(0.5*identity - dlp)            

            rhs =-(monopole_fun)    

            

            boundData, info = bempp.api.linalg.gmres(lhs, rhs, tol=1E-5)

I then save my boundData to a numpy array (simulation_data is boundData.coefficients) and load it with:

boundData = bempp.api.GridFunction(self.space, coefficients=simulation_data).

Hope there is enough info now.

Thanks once again for your help.

Luiz

[Matthew Scroggs]

Hi Luiz,

I've tried running this on a small cube mesh and it seems to work ok for me. Could you check try making a minimal failing example so I can try running that to see what's up?

Thanks,

Matthew

[Luiz Augusto Alvim]

Hello Mathew, I figured out why it is not working. It had something to do with the space variable and its saving and loading. I wish I could elaborate more on the solution, but quite frankly I did not fully understand why it worked. If someone else in the forum has similar issues, just contact me and we can figure it out.

Thanks you so much for your help,

Best Regards,

Luiz