[petsc-users] dmplex face normals orientation
Ingo Gaertner
ingogaertner.tus at gmail.com
Thu Apr 20 00:41:43 CDT 2017
Thank you, Matt, this answers my question.
Ingo
2017-04-18 22:21 GMT+02:00 Matthew Knepley <knepley at gmail.com>:
> 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_8361884064281082256_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/20170420/df146ada/attachment-0001.html>
More information about the petsc-users
mailing list