Commits

Matt Knepley  committed 7ff4606

DMPlex: Reorganized geometry testing

  • Participants
  • Parent commits 4017564

Comments (0)

Files changed (3)

File config/builder.py

                                                                #{'numProcs': 3, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
                                                                #{'numProcs': 5, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'}
 ],
-                        'src/ts/examples/tutorials/ex11':      [{'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo'},
+                        'src/ts/examples/tutorials/ex11':      [# 2D 0-6
+                                                                {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo'},
                                                                 {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside-quad-15.exo'},
                                                                 {'numProcs': 2, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo'},
                                                                 {'numProcs': 2, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside-quad-15.exo'},
                                                                 {'numProcs': 8, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside-quad.exo'},
                                                                 {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo -ts_type rosw'},
                                                                 {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/squaremotor-30.exo -ufv_split_faces'},
+                                                                # 3D 7
                                                                 {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/blockcylinder-50.exo -bc_inflow 100,101,200 -bc_outflow 201'}],
                         }
 

File src/dm/impls/plex/examples/tests/ex8.c

 #include <petscdmplex.h>
 
 #undef __FUNCT__
+#define __FUNCT__ "ChangeCoordinates"
+PetscErrorCode ChangeCoordinates(DM dm, PetscInt spaceDim, PetscScalar vertexCoords[])
+{
+  PetscSection   coordSection;
+  Vec            coordinates;
+  PetscScalar   *coords;
+  PetscInt       vStart, vEnd, v, d, coordSize;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
+  ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr);
+  for (v = vStart; v < vEnd; ++v) {
+    ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
+  ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
+  ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
+  for (v = vStart; v < vEnd; ++v) {
+    PetscInt off;
+
+    ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
+    for (d = 0; d < spaceDim; ++d) {
+      coords[off+d] = vertexCoords[(v-vStart)*spaceDim+d];
+    }
+  }
+  ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
+  ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
+  ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
+  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CheckFEMGeometry"
+PetscErrorCode CheckFEMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx)
+{
+  PetscReal      v0[3], J[9], invJ[9], detJ;
+  PetscInt       d, i, j;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexComputeCellGeometry(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr);
+  for (d = 0; d < spaceDim; ++d) {
+    if (v0[d] != v0Ex[d]) {
+      switch (spaceDim) {
+      case 2: SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g) != (%g, %g)", v0[0], v0[1], v0Ex[0], v0Ex[1]);break;
+      case 3: SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g, %g) != (%g, %g, %g)", v0[0], v0[1], v0[2], v0Ex[0], v0Ex[1], v0Ex[2]);break;
+      default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid space dimension %d", spaceDim);
+      }
+    }
+  }
+  for (i = 0; i < spaceDim; ++i) {
+    for (j = 0; j < spaceDim; ++j) {
+      if (fabs(J[i*spaceDim+j]    - JEx[i*spaceDim+j])    > 1.0e-9) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J[%d,%d]: %g != %g", i, j, J[i*spaceDim+j], JEx[i*spaceDim+j]);
+      if (fabs(invJ[i*spaceDim+j] - invJEx[i*spaceDim+j]) > 1.0e-9) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ[%d,%d]: %g != %g", i, j, invJ[i*spaceDim+j], invJEx[i*spaceDim+j]);
+    }
+  }
+  if (fabs(detJ - detJEx) > 1.0e-9) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g != %g", detJ, detJEx);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CheckFVMGeometry"
+PetscErrorCode CheckFVMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx)
+{
+  PetscReal      centroid[3], normal[3], vol;
+  PetscInt       d;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexComputeCellGeometryFVM(dm, cell, &vol, centroid, normal);CHKERRQ(ierr);
+  for (d = 0; d < spaceDim; ++d) {
+    if (fabs(centroid[d] - centroidEx[d]) > 1.0e-9) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid centroid[%d]: %g != %g", d, centroid[d], centroidEx[d]);
+    if (fabs(normal[d] - normalEx[d]) > 1.0e-9) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid normal[%d]: %g != %g", d, normal[d], normalEx[d]);
+  }
+  if (fabs(volEx - vol) > 1.0e-9) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid volume = %g != %g", vol, volEx);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "TestTriangle"
 PetscErrorCode TestTriangle(MPI_Comm comm, PetscBool interpolate)
 {
   DM             dm;
-  PetscReal      v0[3], J[9], invJ[9], detJ;
-  PetscReal      centroid[3], normal[3], vol;
-  PetscInt       dim, i, j;
+  PetscInt       dim;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
     ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
   }
   /* Check reference geometry: determinant is scaled by reference volume (2.0) */
-  ierr = DMPlexComputeCellGeometry(dm, 0, v0, J, invJ, &detJ);CHKERRQ(ierr);
-  if ((v0[0] != -1.0) || (v0[1] != -1.0)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g)", v0[0], v0[1]);
-  for (i = 0; i < dim; ++i) {
-    for (j = 0; j < dim; ++j) {
-      if (fabs(J[i*dim+j] - (i == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J");
-      if (fabs(invJ[i*dim+j] - (i == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ");
-    }
-  }
-  if (fabs(detJ - 1.0) > 1.0e-9) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g should be 1.0", detJ);
-  if (interpolate) {
-    ierr = DMPlexComputeCellGeometryFVM(dm, 0, &vol, centroid, normal);CHKERRQ(ierr);
-    for (i = 0; i < dim; ++i) {
-      if (fabs(centroid[i] + 0.333333) > 1.0e-6) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid centroid[%d], %g != -0.333333 (%g)", i, centroid[i]);
-      if (fabs(normal[i]) > 1.0e-9) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid normal[%d], %g != 0.0", i, normal[i]);
-    }
-    if (fabs(detJ*2.0 - vol) > 1.0e-9) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid volume = %g should be 2.0", vol);
+  {
+    PetscReal v0Ex[2]       = {-1.0, -1.0};
+    PetscReal JEx[4]        = {1.0, 0.0, 0.0, 1.0};
+    PetscReal invJEx[4]     = {1.0, 0.0, 0.0, 1.0};
+    PetscReal detJEx        = 1.0;
+    PetscReal centroidEx[2] = {-0.333333333333, -0.333333333333};
+    PetscReal normalEx[2]   = {0.0, 0.0};
+    PetscReal volEx         = 2.0;
+
+    ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+    if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
   }
   /* Check random triangles: rotate and translate */
 
-  /* Move to 3D */
+  /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (2.0) */
   dim = 3;
   {
-    PetscSection coordSection;
-    Vec          coordinates;
-    PetscScalar  vertexCoords[9] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0};
-    PetscScalar *coords;
-    PetscInt     vStart, vEnd, v, d, coordSize, spaceDim = 3;
-
-    ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
-    ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
-    ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr);
-    for (v = vStart; v < vEnd; ++v) {
-      ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr);
-      ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr);
-    }
-    ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
-    ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
-    ierr = VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);CHKERRQ(ierr);
-    ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
-    ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
-    ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
-    ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
-    for (v = vStart; v < vEnd; ++v) {
-      PetscInt off;
-
-      ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
-      for (d = 0; d < spaceDim; ++d) {
-        coords[off+d] = vertexCoords[(v-vStart)*spaceDim+d];
-      }
-    }
-    ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
-    ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
-    ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
-    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
-  }
-  /* Check reference geometry: determinant is scaled by reference volume (2.0) */
-  ierr = DMPlexComputeCellGeometry(dm, 0, v0, J, invJ, &detJ);CHKERRQ(ierr);
-  if ((v0[0] != -1.0) || (v0[1] != -1.0) || (v0[2] != 0.0)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0(%g, %g, %g)", v0[0], v0[1], v0[2]);
-  for (i = 0; i < dim; ++i) {
-    for (j = 0; j < dim; ++j) {
-      if (fabs(J[i*dim+j] - (i == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J[%d,%d]", i, j);
-      if (fabs(invJ[i*dim+j] - (i == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ[%d,%d]", i, j);
-    }
-  }
-  if (fabs(detJ - 1.0) > 1.0e-9) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g should be 1.0", detJ);
-  if (interpolate) {
-    ierr = DMPlexComputeCellGeometryFVM(dm, 0, &vol, centroid, normal);CHKERRQ(ierr);
-    for (i = 0; i < dim; ++i) {
-      if (fabs(centroid[i] - (i < 2 ? -0.333333 : 0.0)) > 1.0e-6) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid centroid[%d], %g != %g", i, centroid[i], (i < 2 ? -0.333333 : 0.0));
-      if (fabs(normal[i] - (i < 2 ? 0.0 : 1.0)) > 1.0e-9) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid normal[%d], %g != %g", i, normal[i], (i < 2 ? 0.0 : 1.0));
-    }
-    if (fabs(detJ*2.0 - vol) > 1.0e-9) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid volume = %g should be 2.0", vol);
+    PetscScalar vertexCoords[9] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0};
+    PetscReal v0Ex[3]       = {-1.0, -1.0, 0.0};
+    PetscReal JEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+    PetscReal invJEx[9]     = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+    PetscReal detJEx        = 1.0;
+    PetscReal centroidEx[3] = {-0.333333333333, -0.333333333333, 0.0};
+    PetscReal normalEx[3]   = {0.0, 0.0, 1.0};
+    PetscReal volEx         = 2.0;
+
+    ierr = ChangeCoordinates(dm, dim, vertexCoords);CHKERRQ(ierr);
+    ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+    if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
   }
   /* Rotated reference element */
   {
-    PetscSection coordSection;
-    Vec          coordinates;
-    PetscScalar  vertexCoords[9] = {0.0, -1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0};
-    PetscScalar *coords;
-    PetscInt     vStart, vEnd, v, d, spaceDim = 3;
-
-    ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
-    ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
-    ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
-    ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
-    for (v = vStart; v < vEnd; ++v) {
-      PetscInt off;
-
-      ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
-      for (d = 0; d < spaceDim; ++d) {
-        coords[off+d] = vertexCoords[(v-vStart)*spaceDim+d];
-      }
-    }
-    ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
-    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
-  }
-  ierr = DMPlexComputeCellGeometry(dm, 0, v0, J, invJ, &detJ);CHKERRQ(ierr);
-  if ((v0[0] != 0.0) || (v0[1] != -1.0) || (v0[2] != -1.0)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0(%g, %g, %g)", v0[0], v0[1], v0[2]);
-  for (i = 0; i < dim; ++i) {
-    for (j = 0; j < dim; ++j) {
-      if (fabs(J[i*dim+j] - ((i+2)%3 == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J");
-      if (fabs(invJ[i*dim+j] - ((i+1)%3 == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ");
-    }
-  }
-  if (fabs(detJ - 1.0) > 1.0e-9) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g should be 1.0", detJ);
-  if (interpolate) {
-    ierr = DMPlexComputeCellGeometryFVM(dm, 0, &vol, centroid, normal);CHKERRQ(ierr);
-    for (i = 0; i < dim; ++i) {
-      if (fabs(centroid[i] - (i > 0 ? -0.333333 : 0.0)) > 1.0e-6) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid centroid[%d], %g != %g", i, centroid[i], (i > 0 ? -0.333333 : 0.0));
-      if (fabs(normal[i] - (i > 0 ? 0.0 : 1.0)) > 1.0e-9) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid normal[%d], %g != %g", i, normal[i], (i > 0 ? 0.0 : 1.0));
-    }
-    if (fabs(detJ*2.0 - vol) > 1.0e-9) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid volume = %g should be 2.0", vol);
+    PetscScalar vertexCoords[9] = {0.0, -1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0};
+    PetscReal v0Ex[3]   = {0.0, -1.0, -1.0};
+    PetscReal JEx[9]    = {0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0};
+    PetscReal invJEx[9] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0};
+    PetscReal detJEx    = 1.0;
+    PetscReal centroidEx[3] = {0.0, -0.333333333333, -0.333333333333};
+    PetscReal normalEx[3]   = {1.0, 0.0, 0.0};
+    PetscReal volEx         = 2.0;
+
+    ierr = ChangeCoordinates(dm, dim, vertexCoords);CHKERRQ(ierr);
+    ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+    if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
   }
   /* Cleanup */
   ierr = DMDestroy(&dm);CHKERRQ(ierr);

File src/dm/impls/plex/examples/tests/output/ex8_0.out

 [0]: 0 <---- 2 (0)
 [0]: 0 <---- 3 (0)
 coordinates with 1 fields
-  field 0 with 2 components
+  field 0 with 3 components
 Process 0:
   (   1) dim  3 offset   0 -1 -1 0
   (   2) dim  3 offset   3 1 -1 0
 [0]: 0 <---- 2 (0)
 [0]: 0 <---- 3 (0)
 coordinates with 1 fields
-  field 0 with 2 components
+  field 0 with 3 components
 Process 0:
   (   1) dim  3 offset   0 0 -1 -1
   (   2) dim  3 offset   3 0 1 -1
 [0]: 6 <---- 3 (0)
 [0]: 6 <---- 1 (0)
 coordinates with 1 fields
-  field 0 with 2 components
+  field 0 with 3 components
 Process 0:
   (   1) dim  3 offset   0 -1 -1 0
   (   2) dim  3 offset   3 1 -1 0
 [0]: 6 <---- 3 (0)
 [0]: 6 <---- 1 (0)
 coordinates with 1 fields
-  field 0 with 2 components
+  field 0 with 3 components
 Process 0:
   (   1) dim  3 offset   0 0 -1 -1
   (   2) dim  3 offset   3 0 1 -1