Thank you, Matt, this answers my question. Ingo
2017-04-18 22:21 GMT+02:00 Matthew Knepley <[email protected]>: > On Tue, Apr 18, 2017 at 2:33 PM, Ingo Gaertner <[email protected] > > 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 <[email protected]>: >> >>> On Tue, Apr 18, 2017 at 12:46 AM, Ingo Gaertner < >>> [email protected]> 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 <[email protected]>: >>>> >>>>> 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 <[email protected]>: >>>>> >>>>>> On Wed, Apr 12, 2017 at 10:52 AM, Ingo Gaertner < >>>>>> [email protected]> 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 >
