defining a discrete function

Andi!

You can just access the vector directly and fill it with what ever values you
want:

  FunctionSpace V(mesh);
  Function f(V);
  Vector& v = V.vector();

  for (int i; i<v.size(); i++)
    v[i] = somevalue(i);

Johan

On Wednesday October 13 2010 11:18:48 andi wrote:
> New question #129296 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/129296
>
> Hello, I would like to ask how someone can define a discrete function using
> a vector which contains the values of the function at the nodes. The demo
> 'advection-diffusion' shows how the values at the nodes can be read from
> an .xml file, but I couldn't find any relevant example on how you can
> define in C++ a vector containing these values. I would be grateful if
> anyone can suggest a solution.
>
> Andi

This is also a problem that I'd also like a solution for. Because, the Vector interface doesn't support the '[]' operator, in the past, I've had to resort to setting the vector via:

void set(const double* block, uint m, const uint* rows)

But, this seems somewhat convoluted.

Therefore, trying what you suggested, Johan, one gets the error:

error: lvalue required as left operand of assignment

Jack

On 14/10/10 11:50, Jack wrote:
> Question #129296 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/129296
>
> Jack requested for more information:
> This is also a problem that I'd also like a solution for. Because, the
> Vector interface doesn't support the '[]' operator, in the past, I've
> had to resort to setting the vector via:
>
> void set(const double* block, uint m, const uint* rows)
>
> But, this seems somewhat convoluted.
>
> Therefore, trying what you suggested, Johan, one gets the error:
>
> error: lvalue required as left operand of assignment
>
> Jack
>
>
Jack, using

void set(const double* block, uint m, const uint* rows)

is the way I got a solution as well. Furthermore, I had to define a
GenericVector, and not a Vector. But can I ask you another question?
To make a proper assignment of the vector values for the function we
want to define, we need to figure out which the dof mapping is. Have you
used this feature? If yes, how you did it?

Andi

On Thursday October 14 2010 06:24:26 andi wrote:
> Question #129296 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/129296
>
> Status: Answered => Open
>
> andi is still having a problem:
>
> On 14/10/10 11:50, Jack wrote:
> > Question #129296 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/129296
> >
> > Jack requested for more information:
> > This is also a problem that I'd also like a solution for. Because, the
> > Vector interface doesn't support the '[]' operator, in the past, I've
> > had to resort to setting the vector via:
> >
> > void set(const double* block, uint m, const uint* rows)
> >
> > But, this seems somewhat convoluted.
> >
> > Therefore, trying what you suggested, Johan, one gets the error:
> >
> > error: lvalue required as left operand of assignment

Yes, I was a little bit too hasty with the suggestion. I am probably poluted
with too much PyDOLFIN coding ;)

Because our Vectors are just thin wrappers around other linear algebra
objects, arent the "set operator[]" implemented. Just the "get operator[]".

You could construct an array and fill this with values too. This can then be
used to set the values of the vector:

  Array<double> values(v.size());

  // Fill values...
  v.set_local(values);

> void set(const double* block, uint m, const uint* rows)
>
> is the way I got a solution as well. Furthermore, I had to define a
> GenericVector, and not a Vector.
>
> But can I ask you another question?
> To make a proper assignment of the vector values for the function we
> want to define, we need to figure out which the dof mapping is. Have you
> used this feature? If yes, how you did it?

You might post another question with a particular example for this question.

Johan

Johan

This works brilliantly. Thanks. 'But', how do you get it to work with MPI. With MPI, you get:

roam146-22:JJ$ mpirun -n 2 ./test
Process 0: Number of global vertices: 10000
Process 0: Number of global cells: 19602
Process 0: Partitioned mesh, edge cut is 102.
Process 1: Partitioned mesh, edge cut is 102.
Assertion failed: (values.size() == local_size), function set_local, file /Users/JJ/FEniCS/src/dolfin-0.9.9/dolfin/la/PETScVector.cpp, line 146.
[roam146-22:05747] *** Process received signal ***
[roam146-22:05747] Signal: Abort trap (6)
[roam146-22:05747] Signal code: (0)
Assertion failed: (values.size() == local_size), function set_local, file /Users/JJ/FEniCS/src/dolfin-0.9.9/dolfin/la/PETScVector.cpp, line 146.

Jack

On Thursday October 14 2010 09:20:00 Jack wrote:
> Question #129296 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/129296
>
> Jack requested for more information:
> Johan
>
> This works brilliantly. Thanks.

:)

> 'But', how do you get it to work with
> MPI. With MPI, you get:

Then you need to reduce the size of values to the local size of the vector,
and associate the correct local range to some global values youhave.

Something like:

  Array<double> values(v.local_size());
  for (int i=0, j=v.local_range.first;
       i<v.local_size(), j<v.local_range.second; i++, j++)
     values[i] = somedata[j];

  v.set_local(values);

Johan

> roam146-22:JJ$ mpirun -n 2 ./test
> Process 0: Number of global vertices: 10000
> Process 0: Number of global cells: 19602
> Process 0: Partitioned mesh, edge cut is 102.
> Process 1: Partitioned mesh, edge cut is 102.
> Assertion failed: (values.size() == local_size), function set_local, file
> /Users/JJ/FEniCS/src/dolfin-0.9.9/dolfin/la/PETScVector.cpp, line 146.
> [roam146-22:05747] *** Process received signal ***
> [roam146-22:05747] Signal: Abort trap (6)
> [roam146-22:05747] Signal code: (0)
> Assertion failed: (values.size() == local_size), function set_local, file
> /Users/JJ/FEniCS/src/dolfin-0.9.9/dolfin/la/PETScVector.cpp, line 146.
>
>
> Jack

Thanks again, Johan, works brilliantly.

If I may just impose on you again.

The only part of my code that now breaks down with MPI is where I have to allocate to each member of a std::vector, the corresponding values from a Function u, i.e.,

int j = 0;
for (std::vector<foo*>::iterator p = instances.begin(); p != instances.end(); ++p) {
 (*p)->set_value( u.vector()[j] );
 ++j;
}

Is there a way to reduce this to the local size of the u vector and associate the correct local range?

Many thanks

Jack

On Thursday October 14 2010 15:34:09 Jack wrote:
> Question #129296 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/129296
>
> Jack requested for more information:
> Thanks again, Johan, works brilliantly.
>
> If I may just impose on you again.
>
> The only part of my code that now breaks down with MPI is where I have
> to allocate to each member of a std::vector, the corresponding values
> from a Function u, i.e.,
>
> int j = 0;
> for (std::vector<foo*>::iterator p = instances.begin(); p !=
> instances.end(); ++p) { (*p)->set_value( u.vector()[j] );
> ++j;
> }

If instances is a std::vector of size v.local_size() then you should be able
to do something like:

   Array<double> values(v.size());
   v.get_local(values);

   for (std::vector<foo*>::iterator p = instances.begin(); p !=
   instances.end(); ++p) { (*p)->set_value( values[j] );
  ++j;
   }

Johan
> Is there a way to reduce this to the local size of the u vector and
> associate the correct local range?
>
> Many thanks
>
> Jack

Hi

What if instances is a std::vector is of size 'V.dofmap().global.dimension()' ?

I'm trying to split it across the processors but it still keeps crashing.

roam146-22:JJ$ mpirun -n 2 ./test
Process 0: Number of global vertices: 10000
Process 0: Number of global cells: 19602
Process 0: Partitioned mesh, edge cut is 102.
Process 1: Partitioned mesh, edge cut is 102.
Process 0: Solving linear system of size 10000 x 10000 (PETSc Krylov solver).
Process 0: PETSc Krylov solver (cg, bjacobi) converged in 7 iterations.
Process 0: PETSc Krylov solver preconditioner (bjacobi) sub-methods: (preonly, ilu)
Process 1: PETSc Krylov solver (cg, bjacobi) converged in 7 iterations.
Process 1: PETSc Krylov solver preconditioner (bjacobi) sub-methods: (preonly, ilu)
Assertion failed: (i < _size), function operator[], file /Users/JJ/FEniCS/include/dolfin/common/Array.h, line 133.

Many thanks

Jack

It is difficult to advice further as it is not clear what you are doing. But
remember accessing pure local data structures are best done in scaled to local
sized data structures.

If you want to keep a global data structure filled with duplicated data points
you need to start using MPI communication, which might not be what you want,
and might qualify for a new question thread.

Johan

On Friday October 15 2010 04:10:26 Jack wrote:
> Question #129296 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/129296
>
> Jack requested for more information:
> Hi
>
> What if instances is a std::vector is of size
> 'V.dofmap().global.dimension()' ?
>
> I'm trying to split it across the processors but it still keeps
> crashing.
>
> roam146-22:JJ$ mpirun -n 2 ./test
> Process 0: Number of global vertices: 10000
> Process 0: Number of global cells: 19602
> Process 0: Partitioned mesh, edge cut is 102.
> Process 1: Partitioned mesh, edge cut is 102.
> Process 0: Solving linear system of size 10000 x 10000 (PETSc Krylov
> solver). Process 0: PETSc Krylov solver (cg, bjacobi) converged in 7
> iterations. Process 0: PETSc Krylov solver preconditioner (bjacobi)
> sub-methods: (preonly, ilu) Process 1: PETSc Krylov solver (cg, bjacobi)
> converged in 7 iterations. Process 1: PETSc Krylov solver preconditioner
> (bjacobi) sub-methods: (preonly, ilu) Assertion failed: (i < _size),
> function operator[], file /Users/JJ/FEniCS/include/dolfin/common/Array.h,
> line 133.
>
> Many thanks
>
> Jack

Thanks Johan. I'll start a new thread.

Also, Andi, my apology for hijacking your thread.

Jack