laplacian

In general escript-finley does not allow 2nd order derivatives due its internal representation of
gradients. you could calcute s=Laplace(sol) by solving the PDE

s = - (- grad(sol)_i)_{,i}

which makes X= - grad(sol) and y= grad(sol)_i * n_i in the escript PDE . This then can be
interpolated to FunctionOnBoundary(domain)

The calculation of the tangential Laplace will not work with this approach.

Thanks Lutz,

I found out a way to reformulate the problem using only the gradient and the normal vectors on the boundary.

After performing the calculations, I have a Data object Y, with values for all points on the boundaries but I only need to update the value of the points tagged with a specific value, so I do:

 for i in xrange( BC.getNumberOfDataPoints() ):
    if BC.getTagNumber(i) == 2:
 Yi = Y.getTupleForDataPoint(i)
        tmp = wavenumber*Z0*Yi[0]
        BC.setValueOfDataPoint(i, [[-tmp,-tmp],[tmp,-tmp]])

  pde.setValue(d=BC)

Is this the proper way to update the BC point by point? For a large model, that loop can become expansive, no?

Yes the loop would be expensive and it also uses
setValueOfDataPoint()
which we recommend avoiding _if possible_.

I haven't actually tried this but perhaps something like the following would do what you want:

# create a function for all points in Y
tmp=Y[0]*wavenumber*Z0
update=tmp*[[-1, -1], [1, -1]]

# create a function which is 1 for tag 2 and zero for everything else
mask=Data(0, BC.getFunctionSpace())
mask.setTaggedValue(2,1)

#selectively copy values into BC where the mask is >0
BC.copyWithMask(update, mask)

Thanks Joel,
The copyWithMask method appears to work properly.

I am still having some issues with that boundary condition I am trying to implement.

I am currently doing:
g_bnd = grad(sol, where = FunctionOnBoundary(domain))

to get the gradient on the boundaries but the component of the gradient normal to the boundary is always 0. Is this the intended behavior?

In my case, I need the normal component. Should I calculate the gradient on Function(domain) instead and use a Projector to get the values on the Face elements?

Thanks again for the help.

I am trying to calculate the laplacian as suggested by Lutz in answer #1 above but I am not getting correct results. Here is my script with a simple 2d test case. The equation is:

laplacian(p) + k**2*p = 0

with n_i*grad(p)_i = 0 on all boundaries but one where it is set to 1.

I am looking forward for your comments.

-----

cs = 343. # speed of sound (m/sec )
rho = 1.25 # density (kg/m**3)
frequency = 100. # Hz

L = 0.3

wavenumber = 2*pi*frequency/cs

domain = Rectangle(l0=L,l1=0.02,n0=30,n1=10, order=2)

pde = LinearPDE(domain, numEquations=1)

A = pde.createCoefficient("A")
for i in range(domain.getDim()):
  A[i,i]=-1

pde.setValue(A=A, D=wavenumber**2)

# normal acceleration BC on input plane
anBC = pde.createCoefficient("y")
anBC.setTaggedValue(1, 1)
pde.setValue(y=anBC)

sol = pde.getSolution()

laplacian = LinearPDE(domain, numEquations=1)
laplacian.setValue(D=1, X=-grad(sol) , y=inner(grad(sol,where=FunctionOnBoundary(domain)), domain.getNormal()))

l = laplacian.getSolution()

saveVTK("p.%04d.vtu"%frequency, p=sol, gp=grad(sol), l=l)

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Will take a look at this problem in the next few days.

Ok. Thanks.

What is the result for sol you are expecting?

Do you refer to my original question or to the comment #5 above?

The equation I solved is simply laplacian(p)+k**2*p=0, so when calculating the gradient I should obtain k**2*p as the answer. This is not very useful but I should be able to make it work. This is not what I currently get.

Then, if it works, I'll try to reformulate the method that you proposed in #1 and do a system of three equations to retrieve the 3 spacial components of the vectorial gradient of grad(sol), from which I should be able to calculate the laplacian parallel to the surface and finally go on with implementing the boundary layer impedance model.

Thanks again for your support.

You need to use rich face element when calling "grad(sol,where=FunctionOnBoundary(domain)" and using the component normal to the boundary, see http://esys.esscc.uq.edu.au/esys13/nightly/user/html/node86.html .
Otherwise only the tangential components are calculated properly while the normal component is set to zero - which in you case means that y=0 in the laplacian PDE.

Unfortunately, Rectangle does not use rich face elements by default. To use rich face element you need to call

    Rectangle(l0=L,l1=0.02,n0=60,n1=20, order=2, useElementsOnFace=True)

Hope this helps.