Add random mesh distribution

Author: rezgarshakeri

examples/Hdiv-mixed/conv_test.sh

@@ -44,7 +44,7 @@ declare -A run_flags
 declare -A test_flags
     test_flags[res_start]=2
     test_flags[res_stride]=1
-    test_flags[res_end]=6
+    test_flags[res_end]=12
 
 file_name=conv_test_result.csv
 

examples/Hdiv-mixed/conv_test_result.csv

@@ -1,6 +1,12 @@
 run,mesh_res,error_u,error_p
-0,2,0.99085,0.05737
-1,3,0.51742,0.03392
-2,4,0.30793,0.02117
-3,5,0.20227,0.01422
-4,6,0.14246,0.01014
+0,2,0.85636,0.09544
+1,3,0.43484,0.05265
+2,4,0.26359,0.03280
+3,5,0.17618,0.02214
+4,6,0.12591,0.01588
+5,7,0.09423,0.01191
+6,8,0.07307,0.00925
+7,9,0.05826,0.00738
+8,10,0.04750,0.00602
+9,11,0.03946,0.00501
+10,12,0.03328,0.00422

examples/Hdiv-mixed/convrate_mixed.png

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

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

@@ -107,6 +107,9 @@ int main(int argc, char **argv) {
   // ---------------------------------------------------------------------------
   DM             dm;
   PetscCall( CreateDM(comm, vec_type, &dm) );
+  // TODO: add mesh option
+  // perturb to have smooth random mesh
+  PetscCall( PerturbVerticesSmooth(dm) );
 
   // ---------------------------------------------------------------------------
   // Setup FE

examples/Hdiv-mixed/main.c

@@ -1,4 +1,5 @@
 #include "../include/setup-dm.h"
+#include "petscerror.h"
 
 // ---------------------------------------------------------------------------
 // Setup DM
@@ -18,4 +19,41 @@ PetscErrorCode CreateDM(MPI_Comm comm, VecType vec_type, DM *dm) {
   PetscCall( DMViewFromOptions(*dm, NULL, "-dm_view") );
 
   PetscFunctionReturn(0);
-};
+};
+
+PetscErrorCode PerturbVerticesSmooth(DM dm) {
+
+  Vec          coordinates;
+  PetscSection coordSection;
+  PetscScalar *coords;
+  PetscInt     v,vStart,vEnd,offset,dim;
+  PetscReal    x,y,z;
+
+  PetscFunctionBeginUser;
+
+  PetscCall( DMGetDimension(dm, &dim) );
+  PetscCall( DMGetCoordinateSection(dm, &coordSection) );
+  PetscCall( DMGetCoordinatesLocal(dm, &coordinates) );
+  PetscCall( DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd) );
+  PetscCall( VecGetArray(coordinates,&coords) );
+  for(v=vStart; v<vEnd; v++) {
+    PetscCall( PetscSectionGetOffset(coordSection,v,&offset) );
+    if(dim==2) {
+      x = coords[offset]; y = coords[offset+1];
+      coords[offset]   = x + 0.06*PetscSinReal(2.0*PETSC_PI*x)*PetscSinReal(
+                           2.0*PETSC_PI*y);
+      coords[offset+1] = y - 0.05*PetscSinReal(2.0*PETSC_PI*x)*PetscSinReal(
+                           2.0*PETSC_PI*y);
+    } else {
+      x = coords[offset]; y = coords[offset+1]; z = coords[offset+2];
+      coords[offset]   = x + 0.03*PetscSinReal(3*PETSC_PI*x)*PetscCosReal(
+                           3*PETSC_PI*y)*PetscCosReal(3*PETSC_PI*z);
+      coords[offset+1] = y - 0.04*PetscCosReal(3*PETSC_PI*x)*PetscSinReal(
+                           3*PETSC_PI*y)*PetscCosReal(3*PETSC_PI*z);
+      coords[offset+2] = z + 0.05*PetscCosReal(3*PETSC_PI*x)*PetscCosReal(
+                           3*PETSC_PI*y)*PetscSinReal(3*PETSC_PI*z);
+    }
+  }
+  PetscCall( VecRestoreArray(coordinates,&coords) );
+  PetscFunctionReturn(0);
+}