[petsc-users] DMDAGetElements and global/local element number

[Dave May]

Please always use "reply-all" so that your messages go to the list.
This is standard mailing list etiquette.  It is important to preserve
threading for people who find this discussion later and so that we do
not waste our time re-answering the same questions that have already
been answered in private side-conversations.  You'll likely get an
answer faster that way too.

On Thu, 12 Sep 2019 at 22:26, Emmanuel Ayala <juaneah at gmail.com> wrote:

> Thank you for the answer.
>
> El jue., 12 de sep. de 2019 a la(s) 15:21, Dave May (
> dave.mayhem23 at gmail.com) escribió:
>
>>
>>
>> On Thu, 12 Sep 2019 at 20:21, Emmanuel Ayala via petsc-users <
>> petsc-users at mcs.anl.gov> wrote:
>>
>>> Hi everyone, it would be great if someone can give me a hint for this
>>> issue, i have been trying to figure out how to solve it, but i did not
>>> succeed
>>>
>>> I'm using DMDA to generate a 3D mesh (DMDA_ELEMENT_Q1). I'm trying to
>>> fill a MPI matrix with some values wich are related to the dofs of each
>>> element node, moreover i need to set this values based on the element
>>> number. Something like:
>>>
>>> mpi_A(total_elements X total_dofs)
>>>
>>>                                                  total_dofs
>>> row_0 (element_0)                      a_0 a_1 a_2 ... a_23
>>> row_1 (element_1)                                     a_0 a_1 a_2 ...
>>> a_23
>>> row_2 (element_2)
>>>     a_0 a_1 a_2 ... a_23
>>> .
>>> .
>>> .
>>> row_n (element_n)                            a_0 a_1 a_2 ... a_23
>>>
>>> The element number is related to the row index. And the matrix values
>>> are set depending of the DOFs related to the element.
>>>
>>> With DMDAGetElements i can read the LOCAL nodes connected to the element
>>> and then the DOFs associated to the element. I can handle the local and
>>> global relations with DMGetLocalToGlobalMapping, MatSetLocalToGlobalMapping
>>> and MatSetValuesLocal. BUT i CAN NOT understand how to know the element
>>> number in LOCAL or GLOBAL contex. DMDAGetElements gives the NUMBER OF
>>> ELEMENTS owned in the local process, but there is not any information about
>>> the local or global ELEMENT NUMBER.
>>>
>>> How to know the local or global element number related to the data
>>> provided by DMDAGetElements?
>>>
>>
>> The DMDA defines cells of the same type (quads (2D) or hex (3D), hence
>> every cell defines the same number of vertices.
>>
>
>> DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
>> nel - number of local elements
>> nen - number of element nodes
>> e - the local indices of the elements' vertices
>>
>> e[] defines the ordering of the elements. e[] is an array containing all
>> of the element-vertex maps. Since each element in the DMDA has the same
>> number of vertices, the first nen values in e[] correspond to the vertices
>> (local index) associated with the first element. The next nen values in e[]
>> correspond to the vertices of the second element.  The vertices for any
>> (local) element with the index "cid" can be sought via e[nen*cid + i] where
>> i would range from 0 to nen-1.
>>
>>
> You are right. I can handle the local information, i think the idea is:
>
>     for ( PetscInt i = 0; i < nel; i++ )
>             for (PetscInt j = 0; j < nen; j++)
>                 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"local element %d :
> e[%d] = %d\n", i, j, e[i*nen+j]);
>
> BUT, it does not give information regarding to the ELEMENT identifier
> (number). I need the element number to ordering the elements inside of a
> MPI matrix. I want to access to each element data by means of the matrix
> row . I mean, in the row_0 there is the information (spreading through the
> columns) of the element_0.
>

I think this is a mis-understanding. The element number is not related to a
row in the matrix. The element is associated with vertices (basis
functions), and each vertex (basis) in the DMDA is given a unique index.
The index of that basis corresponds to a row (column) if it's a test
(trial) function. So if you have any element defined by the array e[], you
know how to insert values into a matrix by using the vertex indices.



>
> The element vertices are numbered starting from 0, for each process. It
> does not give information about the element number.
>
> Why would you ever want, or need, the global element number? What is the
>> use case?
>>
>
> I'm performing topology optimization, and it is part of gradient
> computation. I'm already have the analytic gradient.
> I need the global element number to link nodal displacements with the
> element.
>

Any nodal displacement, except those at the corners of your physical domain
are associated with multiple elements. There isn't a one-to-one map between
nodes and elements.



> I can use a local element number just if I have a equivalence between this
> local number and the global element number.
>

I obviously don't understand what you want to do.

However, here is one way to achieve what you are asking for: specifically
relating local element indices to global indices.
Assuming by "global element number" you are referring to what PETSc calls
the "natural" ordering, then the dumbest way to convert from the local
element index to the natural element index is the following:
(i) upon creation, use DMSetUniformCoordinates to define a unit 1 box;
(ii) compute the cell dimensions dx, dy associated with your uniform grid
layout;
(iii) traverse through the e[] array return by DMDAGetElements. For each
element, get the vertices and compute the centroid cx, cy and then compute
  PetscInt J = (PetscInt)(cy/dy);
  PetscInt I = (PetscInt)(cx/dx);
  PetscInt natural_id = I + J * mx;
where mx is the number of elements in the i direction in your domain (not
the sub-domain).
You can determine mx by calling

DMDAGetInfo(dm,NULL,&mx,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
(iv) after you have computed natural_id for every element, layout the
global coordinates of the DMDA however you want. Following that, you must
call the following

  DM dm, cdm;
  Vec coor,lcoor;

  DMGetCoordinateDM(dm,&cdm);
  DMGetCoordinates(dm,&coor);
  DMGetCoordinatesLocal(dm,&lcoor);
  DMGlobalToLocal(cmd,coor,INSERT_VALUES,lcoor);

The above must be executed to ensure the new coordinates values are
propagated to coords associated with your local sub-domain.

I doubt this will solve your _actual_ problem.

Thanks,
Dave


>
>> Thanks,
>> Dave
>>
>>
>>
>>
>>>
>>> Thank you.
>>>
>>
> Thanks!
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190912/af349703/attachment-0001.html>