[petsc-users] dmplex face normals orientation
Matthew Knepley
knepley at gmail.com
Tue Apr 18 15:21:28 CDT 2017
On Tue, Apr 18, 2017 at 2:33 PM, Ingo Gaertner <ingogaertner.tus at gmail.com>
wrote:
> The part that does not make sense is:
> The code calculates that
> face 11 (or edge 11 if you call the boundary of a 2D cell an edge) has the
> centroid c=(0.750000 1.000000 0.000000) and the normal n=(0.000000 0.500000
> 0.000000) and that
> face 12 (or edge 12) has the centroid c=(0.000000 0.500000 0.000000) and
> the normal n=(-1.000000 -0.000000 0.000000).
> I understood your previous answer ("The normal should be outward, unless
> the face has orientation < 0") such that I have to reverse these normals to
> get an outward normal, because faces 11 and 12 have orientation=-2. But the
> normals n=(0.000000 0.500000 0.000000) and n=(-1.000000 -0.000000
> 0.000000) do already point outward. If I reverse them, they point inward.
>
> I need a CONSISTENT rule how to make use of the normals obtained from
> DMPlexComputeGeometryFVM(dm, &cellgeom,&facegeom) to calculate OUTWARD
> pointing normals with respect to EACH INDIVIDUAL cell.
> Or do I have to iterate over each face of each cell and use a geometric
> check to see if the calculated normals are pointing inward or outward?
>
I apologize. I did not understand the question before. The convention I
have chosen might not be the best one,
but it seems appropriate for FVM.
The orientation of a face normal is chosen such that
n . r > 0
where r is the vector from the centroid of the left face to
the centroid of the right face. If we do GetSupport() for
a face we get back {l, r}, where r could be empty, meaning
the cell at infinity.
This convention means that I never have to check orientations
when I am using the facegeom[] stuff in FVM, I just need to have
{l, r} which I generally have in the loop.
Thanks,
Matt
> Thanks
> Ingo
>
>
>
> 2017-04-18 20:20 GMT+02:00 Matthew Knepley <knepley at gmail.com>:
>
>> On Tue, Apr 18, 2017 at 12:46 AM, Ingo Gaertner <
>> ingogaertner.tus at gmail.com> wrote:
>>
>>> Dear Matt,
>>> please explain your previous answer ("The normal should be outward,
>>> unless the face has orientation < 0") with respect to my example.
>>> As I think, the example shows that the face normals for faces 11 and 12
>>> are pointing outward, although they have orientation=-2. For all other
>>> faces the normal direction and their orientation sign agree with what you
>>> said.
>>>
>>
>> 1) You should destroy the objects at the end
>>
>> ierr = VecDestroy(&cellgeom);CHKERRQ(ierr);
>> ierr = VecDestroy(&facegeom);CHKERRQ(ierr);
>> ierr = DMDestroy(&dm);CHKERRQ(ierr);
>>
>> 2) You should call
>>
>> ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
>> ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
>>
>> after DM creation to make it easier to customize and debug. Then we can
>> use
>>
>> -dm_view ::ascii_info_detail
>>
>> to look at the DM.
>>
>> 3) Lets look at Cell 0
>>
>> [0]: 0 <---- 8 (0)
>> [0]: 0 <---- 13 (0)
>> [0]: 0 <---- 10 (-2)
>> [0]: 0 <---- 12 (-2)
>>
>> There are two edges reversed. The edges themselves {8, 13, 10, 12}
>> should proceed counter-clockwise from the bottom, and have vertices
>>
>> [0]: 8 <---- 2 (0)
>> [0]: 8 <---- 3 (0)
>> [0]: 13 <---- 3 (0)
>> [0]: 13 <---- 6 (0)
>> [0]: 10 <---- 5 (0)
>> [0]: 10 <---- 6 (0)
>> [0]: 12 <---- 2 (0)
>> [0]: 12 <---- 5 (0)
>>
>> so we get as we expect
>>
>> 2 --> 3 --> 6 --> 5
>>
>> which agrees with the coordinates
>>
>> ( 2) dim 2 offset 0 0. 0.
>> ( 3) dim 2 offset 2 0.5 0.
>> ( 4) dim 2 offset 4 1. 0.
>> ( 5) dim 2 offset 6 0. 1.
>> ( 6) dim 2 offset 8 0.5 1.
>> ( 7) dim 2 offset 10 1. 1.
>>
>> Which part does not make sense?
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> Thanks
>>> Ingo
>>>
>>> 2017-04-14 11:28 GMT+02:00 Ingo Gaertner <ingogaertner.tus at gmail.com>:
>>>
>>>> Thank you, Matt,
>>>> as you say, the face orientations do change sign when switching between
>>>> the two adjacent cells. (I confused my program output. You are right.) But
>>>> it seems not to be always correct to keep the normal direction for
>>>> orientation>=0 and invert it for orientation<0:
>>>>
>>>> I include sample code below to make my question more clear. The program
>>>> creates a HexBoxMesh split horizontally in two cells (at x=0.5) It produces
>>>> this output:
>>>>
>>>> "Face centroids (c) and normals(n):
>>>> face #008 c=(0.250000 0.000000 0.000000) n=(-0.000000 -0.500000
>>>> 0.000000)
>>>> face #009 c=(0.750000 0.000000 0.000000) n=(-0.000000 -0.500000
>>>> 0.000000)
>>>> face #010 c=(0.250000 1.000000 0.000000) n=(0.000000 0.500000 0.000000)
>>>> face #011 c=(0.750000 1.000000 0.000000) n=(0.000000 0.500000 0.000000)
>>>> face #012 c=(0.000000 0.500000 0.000000) n=(-1.000000 -0.000000
>>>> 0.000000)
>>>> face #013 c=(0.500000 0.500000 0.000000) n=(1.000000 0.000000 0.000000)
>>>> face #014 c=(1.000000 0.500000 0.000000) n=(1.000000 0.000000 0.000000)
>>>> Cell faces orientations:
>>>> cell #0, faces:[8 13 10 12] orientations:[0 0 -2 -2]
>>>> cell #1, faces:[9 14 11 13] orientations:[0 0 -2 -2]"
>>>>
>>>> Looking at the face normals, all boundary normals point outside (good).
>>>> The normal of face #013 points outside with respect to the left cell
>>>> #0, but inside w.r.t. the right cell #1.
>>>>
>>>> Face 13 is shared between both cells. It has orientation 0 for cell #0,
>>>> but orientation -2 for cell #1 (good).
>>>> What I don't understand is the orientation of face 12 (cell 0) and of
>>>> face 11 (cell 1). These are negative, which would make them point into the
>>>> cell. Have I done some other stupid mistake?
>>>>
>>>> Thanks
>>>> Ingo
>>>>
>>>> Here is the code:
>>>>
>>>> static char help[] = "Check face normals orientations.\n\n";
>>>> #include <petscdmplex.h>
>>>> #undef __FUNCT__
>>>> #define __FUNCT__ "main"
>>>>
>>>> int main(int argc, char **argv)
>>>> {
>>>> DM dm;
>>>> PetscErrorCode ierr;
>>>> PetscFVCellGeom *cgeom;
>>>> PetscFVFaceGeom *fgeom;
>>>> Vec cellgeom,facegeom;
>>>> int dim=2;
>>>> int cells[]={2,1};
>>>> int cStart,cEnd,fStart,fEnd;
>>>> int coneSize,supportSize;
>>>> const int *cone,*coneOrientation;
>>>>
>>>> ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
>>>> ierr = PetscOptionsGetInt(NULL, NULL, "-dim", &dim,
>>>> NULL);CHKERRQ(ierr);
>>>> ierr = DMPlexCreateHexBoxMesh(PETSC_COMM_WORLD, dim,
>>>> cells,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
>>>> &dm);CHKERRQ(ierr);
>>>>
>>>> ierr = DMPlexComputeGeometryFVM(dm, &cellgeom,&facegeom);CHKERRQ(i
>>>> err);
>>>> ierr = VecGetArray(cellgeom, (PetscScalar**)&cgeom);CHKERRQ(ierr);
>>>> ierr = VecGetArray(facegeom, (PetscScalar**)&fgeom);CHKERRQ(ierr);
>>>> ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
>>>> ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
>>>>
>>>> fprintf(stderr,"Face centroids (c) and normals(n):\n");
>>>> for (int f=fStart;f<fEnd;f++){
>>>> fprintf(stderr,"face #%03d c=(%03f %03f %03f) n=(%03f %03f
>>>> %03f)\n",f,
>>>> fgeom[f-fStart].centroid[0],fgeom[f-fStart].centroid[1],fgeo
>>>> m[f-fStart].centroid[2],
>>>> fgeom[f-fStart].normal[0],fgeom[f-fStart].normal[1],fgeom[f-
>>>> fStart].normal[2]);
>>>> }
>>>>
>>>> fprintf(stderr,"Cell faces orientations:\n");
>>>> for (int c=cStart;c<cEnd;c++){
>>>> ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr);
>>>> ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
>>>> ierr = DMPlexGetConeOrientation(dm, c,
>>>> &coneOrientation);CHKERRQ(ierr);
>>>> if (dim==2){
>>>> if (coneSize!=4){
>>>> fprintf(stderr,"Expected coneSize 4, got %d.\n",coneSize);
>>>> exit(1);
>>>> }
>>>> fprintf(stderr,"cell #%d, faces:[%d %d %d %d] orientations:[%d
>>>> %d %d %d]\n",c,
>>>> cone[0],cone[1],cone[2],cone[3],
>>>> coneOrientation[0],coneOrientation[1],coneOrientation[2],con
>>>> eOrientation[3]
>>>> );
>>>> } else if (dim==3){
>>>> if (coneSize!=6){
>>>> fprintf(stderr,"Expected coneSize 6, got %d.\n",coneSize);
>>>> exit(1);
>>>> }
>>>> fprintf(stderr,"cell #%d, faces:[%d %d %d %d %d %d]
>>>> orientations:[%d %d %d %d %d %d]\n",c,
>>>> cone[0],cone[1],cone[2],cone[3],cone[4],cone[5],
>>>> coneOrientation[0],coneOrientation[1],coneOrientation[2],con
>>>> eOrientation[3],coneOrientation[4],coneOrientation[5]
>>>> );
>>>> } else {
>>>> fprintf(stderr,"Dimension %d not implemented.\n",dim);
>>>> exit(1);
>>>> }
>>>> }
>>>> ierr = PetscFinalize();
>>>>
>>>> }
>>>>
>>>>
>>>> 2017-04-14 11:00 GMT+02:00 Matthew Knepley <knepley at gmail.com>:
>>>>
>>>>> On Wed, Apr 12, 2017 at 10:52 AM, Ingo Gaertner <
>>>>> ingogaertner.tus at gmail.com> wrote:
>>>>>
>>>>>> Hello,
>>>>>> I have problems determining the orientation of the face normals of a
>>>>>> DMPlex.
>>>>>>
>>>>>> I create a DMPlex, for example with DMPlexCreateHexBoxMesh().
>>>>>> Next, I get the face normals using DMPlexComputeGeometryFVM(DM dm,
>>>>>> Vec *cellgeom, Vec *facegeom). facegeom gives the correct normals, but I
>>>>>> don't know how the inside/outside is defined with respect to the adjacant
>>>>>> cells?
>>>>>>
>>>>>
>>>>> The normal should be outward, unless the face has orientation < 0.
>>>>>
>>>>>
>>>>>> Finally, I iterate over all cells. For each cell I iterate over the
>>>>>> bounding faces (obtained from DMPlexGetCone) and try to obtain their
>>>>>> orientation with respect to the current cell using
>>>>>> DMPlexGetConeOrientation(). However, the six integers for the orientation
>>>>>> are the same for each cell. I expect them to flip between neighbour cells,
>>>>>> because if a face normal is pointing outside for any cell, the same normal
>>>>>> is pointing inside for its neighbour. Apparently I have a misunderstanding
>>>>>> here.
>>>>>>
>>>>>
>>>>> I see the orientations changing sign for adjacent cells. Want to send
>>>>> a simple code? You should see this
>>>>> for examples. You can run SNES ex12 with -dm_view ::ascii_info_detail
>>>>> to see the change in sign.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Matt
>>>>>
>>>>>
>>>>>> How can I make use of the face normals in facegeom and the
>>>>>> orientation values from DMPlexGetConeOrientation() to get the outside face
>>>>>> normals for each cell?
>>>>>>
>>>>>> Thank you
>>>>>> Ingo
>>>>>>
>>>>>>
>>>>>> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail> Virenfrei.
>>>>>> www.avast.com
>>>>>> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail>
>>>>>> <#m_1492622701141221377_m_-6171374370410339176_m_8179380476569161657_m_2360740100696784902_m_-1897645241821352774_m_-6432344443275881332_DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> What most experimenters take for granted before they begin their
>>>>> experiments is infinitely more interesting than any results to which their
>>>>> experiments lead.
>>>>> -- Norbert Wiener
>>>>>
>>>>
>>>>
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170418/0592cafe/attachment-0001.html>
More information about the petsc-users
mailing list