Re: [Getfem-users] Getfem python interface questions

[LUK ShunTim]

Julien Pommier wrote:
> Quoting LUK ShunTim <address@hidden>:
>
> > Can you explain what do the arguments in mfd.eval('[0, x[1]*x[1]/1000]')
> > in this line
> >
> > b1.set('param', 'M', mfd, mfd.eval('[0, x[1]*x[1]/1000]'))
>
> It evaluates the vector field [0, x^2/1000] on the mesh_fem mfd. There is no
> magic in the eval function, its source code is really short:
>
>     def eval(self, expression):
>         """interpolate an expression on the (lagrangian) MeshFem.
>
>         Examples:
>           mf.eval('x[0]*x[1]') interpolates the function 'x*y'
>           mf.eval('[x[0],x[1]]') interpolates the vector field '[x,y]'
>         """
>         P=self.dof_nodes()
>         nbd = P.shape[1];
>
>         if not self.is_lagrangian:
>             raise RuntimeError('cannot eval on a non-Lagragian MeshFem')
>         if self.qdim() != 1:
>             raise RuntimeError('only works (for now) with qdim == 1')
>         x=P[:,0]; r=numarray.array(eval(expression))
>         Z=numarray.zeros(r.shape + (nbd,),'d')
>         for i in range(0,nbd):
>             x=P[:,i]
>             Z[...,i]=eval(expression)
>         return Z
>
>
> > Does MeshFem.eval() requires a numarray array argument?
>
> No, it is just a string expression that is dynamically evaluated for each dof
> node. For vector fields, just enclose the expression in brackets.
>
>
> --
> Julien

Thanks Julien,

I've found out a silly typo which caused my previous error. Here's the 
complete python script converted from demo_laplacian.m

<demo_laplacian.py>
# Converted from demo_laplacian.m
# with corresponding corresponding matlab commands commented out.

import getfem as gf
from numarray import *

#m = gf_mesh('cartesian',[0:.1:1],[0:.1:1]);
m = gf.Mesh('cartesian', arange(0.0, 1.1, 0.1), arange(0.0, 1.1, 0.1))

#mf = gf_mesh_fem(m,1);
mf = gf.MeshFem(m,1)

#gf_mesh_fem_set(mf,'fem',gf_fem('FEM_QK(2,2)'));
mf.set_fem(gf.Fem('FEM_QK(2,2)'))

#mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,10)'));
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS_PARALLELEPIPED(2,10)'))

#border = gf_mesh_get(m,'outer faces');
border = m.outer_faces()

#gf_mesh_set(m, 'boundary', 1, border);
m.set_region(1, border)

#Uexact = gf_mesh_fem_get(mf, 'eval', { 'y.*(y-1).*x.*(x-1)+x.^5' });
Uexact = mf.eval('x[1]*(x[1]-1.0)*x[0]*(x[0]-1.0)+x[0]**5')

#F = gf_mesh_fem_get(mf, 'eval', {'-(2*(x.^2+y.^2)-2*x-2*y+20*x.^3)' });
F = mf.eval('-(2*(x[0]**2 + x[1]**2) - 2*x[0] - 2*x[1] + 20*x**3)')

#b0=gf_mdbrick('generic elliptic',mim,mf);
b0=gf.MdBrick('generic elliptic',mim, mf)

#b1=gf_mdbrick('dirichlet', b0, 1, mf,'penalized');
b1=gf.MdBrick('dirichlet', b0, 1, mf, 'penalized')

#gf_mdbrick_set(b1, 'param', 'R', mf, Uexact);
b1.set_param('R', mf, Uexact)

#b2=gf_mdbrick('source term',b1);
b2=gf.MdBrick('source term', b1)

#gf_mdbrick_set(b2, 'param', 'source_term', mf, F);
b2.set_param('source_term', mf, F)

#mds=gf_mdstate(b1);
mds=gf.MdState(b1);

#gf_mdbrick_get(b2, 'solve', mds)
b2.solve(mds)

#U=gf_mdstate_get(mds, 'state');
U=mds.state()

#disp(sprintf('H1 norm of error: %g', \
#gf_compute(mf,U-Uexact,'H1 norm',mim)));
print 'H1 norm of error: %15.5f' % gf.compute_H1_norm(mf, U-Uexact, mim)
</demo_laplacian.py>

When the line
        b1.set_param('R', mf, Uexact)
is executed, I got this error traceback

"Traceback (most recent call last):
  File "demo_laplacian.py", line 51, in ?
    b1.set_param('R', mf, Uexact)
  File 
"/home/stluk/getfem/lib/python2.4/site-packages/getfem/getfem.py", line 
1247, in set_param
    return self.set("param", name, *args)
  File 
"/home/stluk/getfem/lib/python2.4/site-packages/getfem/getfem.py", line 
1109, in set
    return getfem('mdbrick_set',self.id, *args)
RuntimeError: (Getfem::InterfaceError) -- wrong size for the parameter 
R, expected an array of size 1x441 ( or 1 for a constant field), got an 
array of size 441"

I'm puzzled as my understanding is there should be no distinction 
between "row-vector" and "column vector" as far as a rank-1 array (in 
numarray) is concerned. What did I miss?

Thanks in advance,
ST
--