1. Vijay Mahadevan
  2. petsc

Commits

Vijay Mahadevan  committed 0add189

Adding an extra argument to box mesh creation in order to specify the domain bounds. Modifying examples accordingly.

  • Participants
  • Parent commits bc27cf6
  • Branches dmmoab

Comments (0)

Files changed (5)

File include/petscdmmoab.h

View file
 PETSC_EXTERN PetscErrorCode DMMoabGetBoundaryMarkers(DM dm,const moab::Range **bdvtx,const moab::Range** bdelems,const moab::Range** bdfaces);
 
 /* DM utility creation interface */
-PETSC_EXTERN PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm,PetscInt,PetscInt,PetscInt,DM*);
+PETSC_EXTERN PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm,PetscInt,const PetscReal*,PetscInt,PetscInt,DM*);
 PETSC_EXTERN PetscErrorCode DMMoabLoadFromFile(MPI_Comm,PetscInt,const char*,const char*,DM*);
 
 #endif

File src/dm/impls/moab/dmmbutil.cxx

View file
   PetscFunctionReturn(0);
 }
 
-static void set_structured_coordinates(PetscInt i, PetscInt j, PetscInt k, PetscReal hxyz, PetscInt vcount, std::vector<double*>& vcoords)
+static void set_structured_coordinates(PetscInt i, PetscInt j, PetscInt k, PetscReal hx, PetscReal hy, PetscReal hz, PetscInt vcount, std::vector<double*>& vcoords)
 {
-  vcoords[0][vcount] = i*hxyz;
-  vcoords[1][vcount] = j*hxyz;
-  vcoords[2][vcount] = k*hxyz;
+  vcoords[0][vcount] = i*hx;
+  vcoords[1][vcount] = j*hy;
+  vcoords[2][vcount] = k*hz;
 }
 
 static void set_element_connectivity(PetscInt dim, moab::EntityType etype, PetscInt offset, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity)
 
 #undef __FUNCT__
 #define __FUNCT__ "DMMoabCreateBoxMesh"
-PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt nele, PetscInt user_nghost, DM *dm)
+PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, const PetscReal* bounds, PetscInt nele, PetscInt user_nghost, DM *dm)
 {
   PetscErrorCode  ierr;
   moab::ErrorCode merr;
   moab::EntityHandle  *connectivity = 0;
   moab::EntityType etype;
   PetscInt    ise[6];
-  PetscReal   xse[6];
+  PetscReal   xse[6],defbounds[6];
   /* TODO: Fix nghost > 0 - now relying on exchange_ghost_cells */
   const PetscInt nghost=0;
 
   /* Compute the co-ordinates of vertices and global IDs */
   std::vector<int>    vgid(locnpts+ghnpts);
   int vcount=0;
-  double hxyz=1.0/nele;
+
+  if (!bounds) { /* default box mesh is defined on a unit-cube */
+    defbounds[0]=0.0; defbounds[1]=1.0;
+    defbounds[2]=0.0; defbounds[3]=1.0;
+    defbounds[4]=0.0; defbounds[5]=1.0;
+    bounds=defbounds;
+  }
+  else {
+    /* validate the bounds data */
+    if(bounds[0] >= bounds[1]) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"X-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[0],bounds[1]);
+    if(dim > 1 && (bounds[2] >= bounds[3])) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Y-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[2],bounds[3]);
+    if(dim > 2 && (bounds[4] >= bounds[5])) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Z-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[4],bounds[5]);
+  }
+
+  const double hx=(bounds[1]-bounds[0])/nele;
+  const double hy=(dim > 1 ? (bounds[3]-bounds[2])/nele : 0.0);
+  const double hz=(dim > 2 ? (bounds[5]-bounds[4])/nele : 0.0);
 
   /* create all the owned vertices */
   for (k = ise[4]; k <= ise[5]; k++) {
     for (j = ise[2]; j <= ise[3]; j++) {
       for (i = ise[0]; i <= ise[1]; i++, vcount++) {
-        set_structured_coordinates(i,j,k,hxyz,vcount,vcoords);
+        set_structured_coordinates(i,j,k,hx,hy,hz,vcount,vcoords);
         vgid[vcount] = (k*npts+j)*npts+i+1;
       }
     }
     for (k = (dim==3?ise[4]-nghost:ise[4]); k <= (dim==3?ise[4]-1:ise[5]); k++) {
       for (j = (dim==2?ise[2]-nghost:ise[2]); j <= (dim==2?ise[2]-1:ise[3]); j++) {
         for (i = (dim>1?ise[0]:ise[0]-nghost); i <= (dim>1?ise[1]:ise[0]-1); i++, vcount++) {
-          set_structured_coordinates(i,j,k,hxyz,vcount,vcoords);
+          set_structured_coordinates(i,j,k,hx,hy,hz,vcount,vcoords);
           vgid[vcount] = (k*npts+j)*npts+i+1;
         }
       }
     for (k = (dim==3?ise[5]+1:ise[4]); k <= (dim==3?ise[5]+nghost:ise[5]); k++) {
       for (j = (dim==2?ise[3]+1:ise[2]); j <= (dim==2?ise[3]+nghost:ise[3]); j++) {
         for (i = (dim>1?ise[0]:ise[1]+1); i <= (dim>1?ise[1]:ise[1]+nghost); i++, vcount++) {
-          set_structured_coordinates(i,j,k,hxyz,vcount,vcoords);
+          set_structured_coordinates(i,j,k,hx,hy,hz,vcount,vcoords);
           vgid[vcount] = (k*npts+j)*npts+i+1;
         }
       }

File src/dm/impls/moab/examples/tests/ex2.c

View file
   moab::Interface*  mbImpl;
   moab::Tag         solndofs;
   moab::Range       ownedvtx;
+  const PetscReal   bounds[2] = {0.0,1.0};
   const char        *fields[2] = {"U","V"};
   PetscScalar       deflt[2]={0.0,0.0};
 
   ierr = Initialize_AppContext(&user);CHKERRQ(ierr);
 
   /* Fill in the user defined work context: */
-  ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, user->n, 1, &dm);CHKERRQ(ierr);
+  ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, bounds, user->n, 1, &dm);CHKERRQ(ierr);
   ierr = DMMoabSetBlockSize(dm, user->nvars);CHKERRQ(ierr);
   ierr = DMMoabSetFields(dm, user->nvars, fields);CHKERRQ(ierr);
   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);

File src/dm/impls/moab/examples/tests/ex4.c

View file
    u = f(x,y) for x = 0, x = 1, y = 0, y = 1
 
 or pure Neumman boundary conditions
+
+Usage:
+    mpiexec -n 2 ./ex4 -bc_type dirichlet -nu .01 -rho .01 -file input/quad_2p.h5m -dmmb_rw_dbg 0 -n 50
+
 */
 
 static char help[] = "\
     ierr = DMMoabLoadFromFile(PETSC_COMM_WORLD, user.dim, user.filename, (np==1 ? "" : ""), &dm);CHKERRQ(ierr);
   }
   else {
-    ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, user.n, 1, &dm);CHKERRQ(ierr);
+    ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, NULL, user.n, 1, &dm);CHKERRQ(ierr);
   }
   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
   ierr = DMMoabSetFields(dm, 1, fields);CHKERRQ(ierr);

File src/dm/impls/moab/examples/tests/ex5.c

View file
   user.delxyz  = 1.0/(user.n);
 
   /* create the DM object for dimension specified */
-  ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, user.n, 1, &dm);CHKERRQ(ierr);
+  ierr = DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, NULL, user.n, 1, &dm);CHKERRQ(ierr);
   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
   ierr = DMMoabSetBlockSize(dm, 2);CHKERRQ(ierr);
   ierr = DMMoabSetFields(dm, 2, fields);CHKERRQ(ierr);