How to Obtain the Coordinates of the Vertices of a Mesh with P2 and Higher Elements

On Mon, Jan 18, 2010 at 09:31:15AM -0000, Ted Kord wrote:
> New question #97807 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/97807
>
> Hello
>
> How do I get the vertex coordinates for a mesh with P2 elements? For a P1 mesh, this works fine:
>
> for (VertexIterator vertex(mesh); !vertex.end(); ++vertex) {
> xy = new Array<double>(dim, geometry.x(vertex->index()));
> x0 = (*xy)[0]; x1 = (*xy)[1];
> std::cout << i << "(" << x0 << ", " << x1 << ")" << std::endl;
> ++i;
> }
>
> It doesn't work for a P2 mesh though.
>
> Some parts of my code which depend on knowing the vertex coordinates ultimately fail with the following error (P2 and above only):
>
> Assertion failed: (V.dofmap().global_dimension() <= x.size()), function Function, file dolfin/function/Function.cpp, line 62.
>
> Could someone please help with this.
>
> Thank you.
>
> Ted

You should check the function tabulate_coordinates in the DofMap class:

  V.dofmap().tabulate_coordinates(...)

This should give you what you need.

For details on what it does, check the UFC manual.

--
Anders

Thanks Anders.

Could you please show me an example of how to use this. I've tried the following but get a segmentation fault.

  for (CellIterator c(mesh); !c.end(); ++c)
    V.dofmap().tabulate_coordinates(&coordinates, *c);

Thank you.

Ted

Hi

Sorry, problem solved. Just me being a bit dense.

Ted

Thanks Anders Logg, that solved my question.

Hi

Is there a function that allows me to remove the duplicate coordinates after doing a :

V.dofmap().tabulate_coordinates(...) ?

Ted

On Mon, Jan 18, 2010 at 11:54:22PM -0000, Ted Kord wrote:
> Question #97807 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/97807
>
> Status: Solved => Open
>
> Ted Kord is still having a problem:
> Hi
>
> Is there a function that allows me to remove the duplicate coordinates
> after doing a :
>
> V.dofmap().tabulate_coordinates(...) ?
>
> Ted

Why are there duplicates? Is V a vector-valued space?

If so, you could try tabulate_coordinates for a scalar subspace of V.

--
Anders

No, V's not a vector-valued space.

By duplicates, I mean that the vertices that are on a common edge between elements are repeated in the enumeration. For example, for a 4x4 rectangle with 2 P2 elements, i.e., nx = ny = 1

    Rectangle mesh(0.0, 0.0, 4.0, 4.0, nx, ny)

V.dofmap().tabulate_coordinates(...) produces:

0: (0, 0)
1: (4, 0)
2: (4, 4)
3: (4, 2)
4: (2, 2)
5: (2, 0)
6: (0, 0)
7: (0, 4)
8: (4, 4)
9: (2, 4)
10: (2, 2)
11: (0, 2)

Vertices 0 and 6, 2 and 8, 4 and 10 are repeated.

I'd like to eliminate these without sorting so that the vertices are still evaluated in the correct order.

Is there already a function to do this or is there a better way?

Ted

Garth N. Wells wrote:
>
> Ted Kord wrote:
>> Question #97807 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/97807
>>
>> Status: Answered => Open
>>
>> Ted Kord is still having a problem:
>> No, V's not a vector-valued space.
>>
>> By duplicates, I mean that the vertices that are on a common edge
>> between elements are repeated in the enumeration. For example, for a
>> 4x4 rectangle with 2 P2 elements, i.e., nx = ny = 1
>>
>> Rectangle mesh(0.0, 0.0, 4.0, 4.0, nx, ny)
>>
>> V.dofmap().tabulate_coordinates(...) produces:
>>
>> 0: (0, 0)
>> 1: (4, 0)
>> 2: (4, 4)
>> 3: (4, 2)
>> 4: (2, 2)
>> 5: (2, 0)
>> 6: (0, 0)
>> 7: (0, 4)
>> 8: (4, 4)
>> 9: (2, 4)
>> 10: (2, 2)
>> 11: (0, 2)
>>
>> Vertices 0 and 6, 2 and 8, 4 and 10 are repeated.
>>
>> I'd like to eliminate these without sorting so that the vertices are
>> still evaluated in the correct order.
>>
>> Is there already a function to do this or is there a better way?
>>
>
> You can insert the dofs into a dolfin::Set. It's like a std::std, but

I meant 'std::set'.

Garth

> the order follows the order of insertion.
>
> Garth
>
>> Ted
>>
>
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

Using a std::set though means that the values get sorted but I'd like to preserve the order in which dolfin evaluates the nodes (if I can).

If that's not possible, I'll probably have to redesign.

Ted

Ted Kord wrote:
> Question #97807 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/97807
>
> Status: Answered => Open
>
> Ted Kord is still having a problem:
> Using a std::set though means that the values get sorted but I'd like to
> preserve the order in which dolfin evaluates the nodes (if I can).
>

That's why I suggested dolfin::Set. It doesn't re-order.

Garth

> If that's not possible, I'll probably have to redesign.
>
> Ted
>

Ahh... I see. I'll have a go with dolfin::Set.

Thanks

Thanks Garth Wells, that solved my question.

I've just added a new function to DofMap:

    /// Return the set of dof indices
    Set<dolfin::uint> list(const Mesh& mesh, bool sort = false) const;

It can be particularly useful for mixed elements, since a for 'shallow' sub-function it will give the entries in the 'big' function vector associated with a particular field, e.g.

    boost::shared_ptr<FunctionSpace> V_sub = V[1];
    Set<dolfin::uint> sub_dof_list = V_sub->dofmap().list(mesh, true);
    for (dolfin::uint i = 0; i < sub_dof_list.size(); ++i )
        cout << "i: " << i << " " << sub_dof_list[i] << endl;

On Fri, Jan 29, 2010 at 07:21:03PM -0000, Garth Wells wrote:
> Question #97807 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/97807
>
> Garth Wells posted a new comment:
> I've just added a new function to DofMap:
>
> /// Return the set of dof indices
> Set<dolfin::uint> list(const Mesh& mesh, bool sort = false) const;
>
> It can be particularly useful for mixed elements, since a for 'shallow'
> sub-function it will give the entries in the 'big' function vector
> associated with a particular field, e.g.
>
> boost::shared_ptr<FunctionSpace> V_sub = V[1];
> Set<dolfin::uint> sub_dof_list = V_sub->dofmap().list(mesh, true);
> for (dolfin::uint i = 0; i < sub_dof_list.size(); ++i )
> cout << "i: " << i << " " << sub_dof_list[i] << endl;

list does not seem to be an optimal name for this function. Something
like dof_set() or just dofs() would be better.

--
Anders