dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Monday 07 July 2008 11:28:40 Martin Sandve Alnæs wrote:
> 2008/7/7 Garth N. Wells <gnw20@xxxxxxxxx>:
> > Anders Logg wrote:
> >> On Fri, Jul 04, 2008 at 11:25:52AM +0200, Martin Sandve Alnæs wrote:
> >>> 2008/7/3 Martin Sandve Alnæs <martinal@xxxxxxxxx>:
> >>>> 2008/7/3 Martin Sandve Alnæs <martinal@xxxxxxxxx>:
> >>>>> I've imple
> >>>>>
> >>>>> 2008/7/3 Dag Lindbo <dag@xxxxxxxxxx>:
> >>>>>> Martin Sandve Alnæs wrote:
> >>>>>>> 2008/7/2 Dag Lindbo <dag@xxxxxxxxxx>:
> >>>>>>>> Martin Sandve Alnæs wrote:
> >>>>>>>>> I just added checks for valid Functions to the assembler, but
> >>>>>>>>> this turned out to break some demos so I changed it to warnings.
> >>>>>>>>
> >>>>>>>> This seems like helpful checks to make!
> >>>>>>>>
> >>>>>>>>> The reason is that Function(mesh, 0.0), which creates a
> >>>>>>>>> ConstantFunction, is used for zero vector functions as well,
> >>>>>>>>> which leaves the value shape of such a function undefined. This
> >>>>>>>>> makes error checking impossible, which is not very nice.
> >>>>>>>>>
> >>>>>>>>> In my opinion, the correct next step is to keep these warnings,
> >>>>>>>>> and turn them into errors after "a while". We should either add
> >>>>>>>>> support for ConstantFunctions with tensor shape (would be nice),
> >>>>>>>>> or make the demos and apps use user-define functions.
> >>>>>>>>
> >>>>>>>> Adding tensor shape to the existing concept of ConstantFunction
> >>>>>>>> seems like the way to go. IMHO, the convenience if this over user
> >>>>>>>> defined functions is too great to discard.
> >>>>>>>>
> >>>>>>>> Would it not just amount to having ConstantFunction::rank() return
> >>>>>>>> the proper value as determined by the mutable member size, instead
> >>>>>>>> of 0? (on line 28 of ConstantFunction.cpp)
> >>>>>>>
> >>>>>>> No, because this size isn't known at construction time, only when
> >>>>>>> interpolate is called with a finite_element object. A suitable (set
> >>>>>>> of) constructor(s) needs to be defined, mirrored in Function, and
> >>>>>>> handled in the swig interface for python.
> >>>>>>
> >>>>>> Yes, I see that this is needed if the rank and shape of the function
> >>>>>> is to be set constructor. Clearly, this makes for a nice user
> >>>>>> interface! My suggestion was more along the lines of a quick fix to
> >>>>>> get rid of the warnings: rank() returns zero if size == 1, returns 1
> >>>>>> if size == 2 or 3 etc.
> >>>>>>
> >>>>>> /Dag
> >>>>>>
> >>>>>>> // Scalar constructor
> >>>>>>> Function f(mesh, 0.0);
> >>>>>>>
> >>>>>>> // General tensor constructor:
> >>>>>>> Function f(Mesh& mesh, uint rank, uint* shape, uint* values);
> >>>>>>>
> >>>>>>> // vector
> >>>>>>> uint rank = 1;
> >>>>>>> uint shape[1] = { 2 };
> >>>>>>> real values[2] = { 1., 2. };
> >>>>>>> Function f(mesh, rank, shape, values);
> >>>>>>> // or
> >>>>>>> uint size = 2;
> >>>>>>> real values[2] = { 1., 2. };
> >>>>>>> Function f(mesh, 1, &size, values);
> >>>>>>>
> >>>>>>> // 3x3 tensor:
> >>>>>>> uint rank = 2;
> >>>>>>> uint shape[2] = { 3, 3 };
> >>>>>>> real values[9] = {
> >>>>>>>    1., 2., 3.,
> >>>>>>>    4., 5., 6.,
> >>>>>>>    7., 8., 9. };
> >>>>>>> Function f(mesh, rank, shape, values);
> >>>>>
> >>>>> I just implemented this, but can't reproduce the old behaviour at the
> >>>>> same time, so the demos must be updated and I guess many user
> >>>>> applications as well. I'd like to have consensus to do this change
> >>>>> before I continue with the demos and check in.
> >>>>>
> >>>>> The typical example is a zero function intended to be a vector:
> >>>>>  Function f(mesh, 0.0);
> >>>>>
> >>>>> Which now would become:
> >>>>>  uint fsize = 3;
> >>>>>  real fvalues[3] = { 0.0, 0.0, 0.0 };
> >>>>>  Function f(mesh, 1, &fsize, fvalues);
> >>>>>
> >>>>> I could evt. add a vector constructor:
> >>>>>  real fvalues[3] = { 0.0, 0.0, 0.0 };
> >>>>>  Function f(mesh, 3, fvalues);
> >>>>
> >>>> And:
> >>>>
> >>>>  Function f(mesh, 3, 0.0);
> >>>>
> >>>> which does the same as the above.
> >>>
> >>> I'd rather not wait a month before I merge, so I'll push these changes.
> >>> Somebody who knows the logic should look at dolfin_function_pre.i
> >>
> >> I think Ola knows these typemaps best.
> >>
> >>> and fix the typemaps there, until then the python syntax for constant
> >>> vector functions will be:
> >>>
> >>> f = Function(mesh, size, array)
> >>>
> >>> E.g.
> >>> f = Function(mesh, 3, numpy.array((1.,2.,3.)))
> >>>
> >>> Or simply:
> >>> f = Function(mesh, 3, 0.0)
> >>
> >> I think that it would be nicer to create a vector-valued function in
> >> C++ by
> >>
> >>   Array values(0.0, 0.0, 0.0);
> >>   Function f(mesh, values);
> >>
> >> and in Python by
> >>
> >>   f = Function(mesh, [0.0, 0.0, 0.0])
> >
> > Looks nice to me.
> >
> > Garth
>
> I've done this and added
>
> Array<uint> shape(2,2);
> Array<real> values(1.0, 2.0,
>                                 3.0, 4.0);
> Function f(mesh, shape, values);
>
> as well.
>
> I kept this notation:
>   Function f(mesh, 3, 0.0);
> since it's shorter than defining the Array first, and the size (3) can
> be a variable.
>
> In Python, the typemaps for Array doesn't seem to work,
> so at the moment the syntax there is now:
>
> f = Function(mesh, 3, 0.0)
> f = Function(mesh, ArrayDouble(1.0, 2.0, 3.0))
> f = Function(mesh, ArrayUInt(2,2), ArrayDouble(1.0, 2.0, 3.0, 4.0))
>
> Note that ArrayUInt(2) will give a length-two array with zeros,
> in contrast with how I used it in the commit message...
> This isn't a problem, since one of the first two lines above
> should be fine for vector functions.

I had a general python-sequence to std::vector<number> typemap laying waiting 
to be used. I have added an in-typemap for any python sequence for Array<uin> 
and Array<real>. The typemap will take any python sequence and copying it to 
either an Array<uint> or an Array<real>. An exception is raised if any of
the elements in the input sequence is not a number.

The following is now possible:
  f = Function(mesh, (1.0, 2.0, 3.0))
  f = Function(mesh, [2,2], numpy.array([1.0, 2.0, 3.0, 4.0]))

and if one do:
  f = Function(mesh, (2,2), (1.0, 2.0, 3.0, "jada"))

one will get 
  <type 'exceptions.ValueError'>: *** Error: Sequence elements must be numbers

I have also out-typemaps for the same structures, then returning a numpy 
array. But as there could potentially be quite large such arrays, from e.g. 
some of the mesh data, and the data is copied, I was not sure of adding this. 
On the other hand an out-typemap will sure make the python interface a bit 
smoother and coping the underlaying data we ensure that no data is either 
lost or leaked.

Johan


> --
> Martin
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev


# HG changeset patch
# User "Johan Hake <hake@xxxxxxxxx>"
# Date 1215547833 -7200
# Node ID 99e2390c5c862f411a04f374d907b1c4b303eaaf
# Parent  509da519281b7cb04f36c1d60fdd6fe1bd608e5a
Added in-typemaps for Array<uint> and Array<real>

The typemap will take any python sequence and copying it to either
an Array<uint> or an Array<real>. An exception is raised if any of
the elements in the input sequence is not a number.

diff -r 509da519281b -r 99e2390c5c86 dolfin/swig/ignores.i
--- a/dolfin/swig/ignores.i	Tue Jul 08 20:45:23 2008 +0100
+++ b/dolfin/swig/ignores.i	Tue Jul 08 22:10:33 2008 +0200
@@ -24,3 +24,4 @@
 %ignore dolfin::cout;
 %ignore dolfin::endl;
 %ignore *::operator<<(unsigned int);
+%ignore dolfin::MeshConnectivity::set(uint entity, uint* connections);
diff -r 509da519281b -r 99e2390c5c86 dolfin/swig/typemaps.i
--- a/dolfin/swig/typemaps.i	Tue Jul 08 20:45:23 2008 +0100
+++ b/dolfin/swig/typemaps.i	Tue Jul 08 22:10:33 2008 +0200
@@ -25,6 +25,69 @@
 %typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY) dolfin::uint* {
     // General typemap
     $1 = PyArray_Check($input) ? 1 : 0;
+}
+
+// Typemap check
+%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY) dolfin::Array<dolfin::real>& {
+    // General typemap
+    $1 = PySequence_Check($input) ? 1 : 0;
+}
+
+// Typemap check
+%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY) dolfin::Array<dolfin::uint>& {
+    // General typemap
+    $1 = PySequence_Check($input) ? 1 : 0;
+}
+
+// Typemap for sending any sequence as input to functions expecting an Array of real
+%typemap(in) const dolfin::Array<dolfin::real>& (dolfin::Array<dolfin::real> tmp) {
+  int i;
+  if (!PySequence_Check($input)) {
+    PyErr_SetString(PyExc_ValueError,"*** Error: Expected a sequence");
+    return NULL;
+  }
+  int pyseq_length = PySequence_Size($input);
+  if (!pyseq_length > 0){
+    PyErr_SetString(PyExc_RuntimeError,"*** Error: Supply a sequence with length > 0");	
+    return NULL;
+  }
+  tmp.reserve(pyseq_length);
+  for (i = 0; i < pyseq_length; i++) {
+    PyObject *o = PySequence_GetItem($input,i);
+    if (PyNumber_Check(o)) {
+      tmp.push_back(static_cast<dolfin::real>(PyFloat_AsDouble(o)));
+    } else {
+      PyErr_SetString(PyExc_ValueError,"*** Error: Sequence elements must be numbers");
+      return NULL;
+    }
+
+  }
+  $1 = &tmp;
+}
+
+// Typemap for sending any sequence as input to functions expecting an Array of uint
+%typemap(in) const dolfin::Array<dolfin::uint>& (dolfin::Array<dolfin::uint> tmp) {
+  int i;
+  if (!PySequence_Check($input)) {
+    PyErr_SetString(PyExc_ValueError,"*** Error: Expected a sequence");
+    return NULL;
+  }
+  int pyseq_length = PySequence_Size($input);
+  if (!pyseq_length > 0){
+    PyErr_SetString(PyExc_RuntimeError,"*** Error: Supply a sequence with length > 0");	
+    return NULL;
+  }
+  tmp.reserve(pyseq_length);
+  for (i = 0; i < pyseq_length; i++) {
+    PyObject *o = PySequence_GetItem($input,i);
+    if (PyNumber_Check(o)) {
+      tmp.push_back(static_cast<dolfin::uint>(PyInt_AsLong(o)));
+    } else {
+      PyErr_SetString(PyExc_ValueError,"*** Error: Sequence elements must be integers");
+      return NULL;
+    }
+  }
+  $1 = &tmp;
 }
 
 // Typemap for sending numpy arrays as input to functions expecting a C array of real