Re: [Getfem-users] Getfem python interface questions

[Julien Pommier]

Hi ShunTim,

Well, you found a bug. I have commited a fix in the svn repository

http://svn.gna.org/viewcvs/getfem/trunk/getfem%2B%2B/interface/src/gf_mdbrick_set.cc?r2=2729&rev=2729&r1=2722&dir_pagestart=50

but of course, you can also cast the 1D array into a 2D array with one row, and
then the interface will be happy.

--
Julien

Quoting LUK ShunTim <address@hidden>:


> 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
> --
>
>