[petsc-users] dmplex face normals orientation
Matthew Knepley
knepley at gmail.com
Tue Apr 18 13:20:39 CDT 2017
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(ierr);
>> 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_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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170418/1e095279/attachment.html>
More information about the petsc-users
mailing list