Added random vertices distribution

Author: rezgarshakeri

examples/Hdiv-mixed/conv_test_result.csv

@@ -1,10 +1,7 @@
 run,mesh_res,error_u,error_p
-0,2,23.34221,0.01392
-1,3,14.88356,0.00793
-2,4,11.04429,0.00485
-3,5,8.90114,0.00337
-4,6,6.49004,0.00260
-5,7,4.87041,0.00207
-6,8,3.99474,0.00165
-7,9,3.24332,0.00136
-8,10,2.65430,0.00114
+0,2,26.96934,0.01516
+1,3,18.26207,0.00982
+2,4,14.36768,0.00532
+3,5,11.16971,0.00438
+4,6,9.12436,0.00348
+5,7,7.85373,0.00257

examples/Hdiv-mixed/convrate_mixed.png

@@ -12,5 +12,6 @@
 // ---------------------------------------------------------------------------
 PetscErrorCode CreateDM(MPI_Comm comm, VecType vec_type, DM *dm);
 PetscErrorCode PerturbVerticesSmooth(DM dm);
+PetscErrorCode PerturbVerticesRandom(DM dm);
 
 #endif // setupdm_h

examples/Hdiv-mixed/include/setup-dm.h

@@ -30,6 +30,7 @@
 // (boundary is not working)
 //   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1
 //   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1,2,3,4
+//   ./main -pc_type svd -problem darcy3d -dm_plex_dim 3 -dm_plex_box_faces 4,4,4 -view_solution -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,0.1,1
 const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n";
 
 #include "main.h"
@@ -83,8 +84,8 @@ int main(int argc, char **argv) {
   // ---------------------------------------------------------------------------
   // -- Initialize backend
   Ceed ceed;
-  //CeedInit("/cpu/self/ref/serial", &ceed);
-  CeedInit(app_ctx->ceed_resource, &ceed);
+  CeedInit("/cpu/self/ref/serial", &ceed);
+  //CeedInit(app_ctx->ceed_resource, &ceed);
   CeedMemType mem_type_backend;
   CeedGetPreferredMemType(ceed, &mem_type_backend);
 
@@ -118,9 +119,11 @@ int main(int argc, char **argv) {
   PetscCall( CreateDM(comm, vec_type, &dm_H1) );
   // TODO: add mesh option
   // perturb to have smooth random mesh
-  PetscCall( PerturbVerticesSmooth(dm) );
-  PetscCall( PerturbVerticesSmooth(dm_H1) );
+  //PetscCall( PerturbVerticesSmooth(dm) );
+  //PetscCall( PerturbVerticesSmooth(dm_H1) );
 
+  PetscCall(PerturbVerticesRandom(dm) );
+  PetscCall(PerturbVerticesRandom(dm_H1) );
   // ---------------------------------------------------------------------------
   // Process command line options
   // ---------------------------------------------------------------------------

examples/Hdiv-mixed/main.c

@@ -91,7 +91,6 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, DM dm,
   PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
   for (PetscInt i=0; i<2; i++) domain_size[i] = domain_max[i] - domain_min[i];
 
-  printf("lx %f, ly %f \n",domain_size[0], domain_size[1]);
   darcy_ctx->kappa = kappa;
   darcy_ctx->rho_a0 = rho_a0;
   darcy_ctx->g = g;

examples/Hdiv-mixed/problems/darcy2d.c

@@ -68,3 +68,59 @@ PetscErrorCode PerturbVerticesSmooth(DM dm) {
   PetscCall( VecRestoreArray(coordinates,&coords) );
   PetscFunctionReturn(0);
 }
+
+
+PetscErrorCode PerturbVerticesRandom(DM dm) {
+
+  PetscFunctionBegin;
+  Vec          coordinates;
+  PetscSection coordSection;
+  PetscScalar *coords;
+  PetscInt     v,vStart,vEnd,offset, dim;
+  PetscReal    x, y, z;
+
+  PetscCall( DMGetDimension(dm,&dim) );
+  PetscCall( DMGetCoordinateSection(dm, &coordSection) );
+  PetscCall( DMGetCoordinatesLocal(dm, &coordinates) );
+  PetscCall( DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd) );
+  PetscCall( VecGetArray(coordinates,&coords) );
+  PetscInt c_end, c_start, num_elem;
+  PetscCall( DMPlexGetHeightStratum(dm, 0, &c_start, &c_end) );
+  num_elem = c_end - c_start;
+
+  for(v=vStart; v<vEnd; v++) {
+    PetscCall( PetscSectionGetOffset(coordSection,v,&offset) );
+    if(dim==2) {
+      PetscScalar nx = sqrt(num_elem);
+      PetscReal domain_min[2], domain_max[2], domain_size[2];
+      PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
+      for (PetscInt i=0; i<2; i++) domain_size[i] = domain_max[i] - domain_min[i];
+      PetscReal hx = domain_size[0]/nx, hy = domain_size[1]/nx;
+      x = coords[offset]; y = coords[offset+1];
+      // perturb randomly O(h*sqrt(2)/3)
+      PetscReal rx = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hx*0.471404);
+      PetscReal ry = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hy*0.471404);
+      PetscReal t = ((PetscReal)rand())/((PetscReal)RAND_MAX)*PETSC_PI;
+      coords[offset  ] = x + rx*PetscCosReal(t);
+      coords[offset+1] = y + ry*PetscSinReal(t);
+    } else {
+      PetscScalar nx = cbrt(num_elem);
+      PetscReal domain_min[3], domain_max[3], domain_size[3];
+      PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
+      for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i];
+      PetscReal hx = domain_size[0]/nx, hy = domain_size[1]/nx,
+                hz = domain_size[2]/nx;
+      x = coords[offset]; y = coords[offset+1], z = coords[offset+2];
+      // This is because 'boundary' is broken in 3D
+      PetscReal rx = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hx*0.471404);
+      PetscReal ry = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hy*0.471404);
+      PetscReal rz = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hz*0.471404);
+      PetscReal t = ((PetscReal)rand())/((PetscReal)RAND_MAX)*PETSC_PI;
+      coords[offset  ] = x + rx*PetscCosReal(t);
+      coords[offset+1] = y + ry*PetscCosReal(t);
+      coords[offset+2] = z + rz*PetscSinReal(t);
+    }
+  }
+  PetscCall( VecRestoreArray(coordinates,&coords) );
+  PetscFunctionReturn(0);
+}