stokes: use -dirichlet LIST for boundary conditions

jedbrown · Apr 24, 2011 · a8ed699 · a8ed699
1 parent ffb6d09
commit a8ed699
Showing 1 changed file with 19 additions and 7 deletions.
diff --git a/src/fs/tests/stokes.c b/src/fs/tests/stokes.c
@@ -22,6 +22,7 @@ static const char help[] = "Solve non-Newtonian Stokes problem using dual order
 
 #include <stokesexact.h>        /* Generated by stokesexact.py */
 
+#define ALEN(a) ((dInt)(sizeof(a)/sizeof(a)[0]))
 static PetscLogEvent LOG_StokesShellMult;
 typedef struct _p_Stokes *Stokes;
 
@@ -60,10 +61,12 @@ struct _p_Stokes {
   IS                     ublock,pblock;
   VecScatter             extractVelocity,extractPressure;
   dInt                   constBDeg,pressureCodim;
-  dBool                  errorview,cardinalMass,neumann300;
+  dBool                  errorview,cardinalMass;
   char                   mattype_A[256],mattype_D[256];
   dQuadratureMethod      function_qmethod,jacobian_qmethod;
   dRulesetIterator       regioniter[EVAL_UB];
+  dInt                   dirichlet[16]; /* Set numbers for Dirichlet conditions, 0 means unused */
+  dBool                  alldirichlet;
 };
 
 static dErr StokesCreate(MPI_Comm comm,Stokes *stokes)
@@ -81,6 +84,10 @@ static dErr StokesCreate(MPI_Comm comm,Stokes *stokes)
   stk->rheo.A        = 1;
   stk->rheo.eps      = 1;
   stk->rheo.p        = 2;
+  stk->dirichlet[0]  = 100;
+  stk->dirichlet[1]  = 200;
+  stk->dirichlet[2]  = 300;
+  stk->alldirichlet  = dTRUE;
   stk->function_qmethod = dQUADRATURE_METHOD_FAST;
   stk->jacobian_qmethod = dQUADRATURE_METHOD_SPARSE;
 
@@ -150,7 +157,14 @@ static dErr StokesSetFromOptions(Stokes stk)
     err = PetscOptionsList("-stokes_D_mat_type","Matrix type for velocity operator","",MatList,stk->mattype_D,stk->mattype_D,sizeof(stk->mattype_D),NULL);dCHK(err);
     err = PetscOptionsEnum("-stokes_f_qmethod","Quadrature method for residual evaluation/matrix-free","",dQuadratureMethods,(PetscEnum)stk->function_qmethod,(PetscEnum*)&stk->function_qmethod,NULL);dCHK(err);
     err = PetscOptionsEnum("-stokes_jac_qmethod","Quadrature to use for Jacobian assembly","",dQuadratureMethods,(PetscEnum)stk->jacobian_qmethod,(PetscEnum*)&stk->jacobian_qmethod,NULL);dCHK(err);
-    err = PetscOptionsBool("-neumann300","Use boundary set 300 as Neumann conditions","",stk->neumann300,&stk->neumann300,NULL);dCHK(err);
+    {
+      dBool flg; dInt n = ALEN(stk->dirichlet);
+      err = PetscOptionsIntArray("-dirichlet","List of boundary sets on which to impose Dirichlet conditions","",stk->dirichlet,&n,&flg);dCHK(err);
+      if (flg) {
+        for (dInt i=n; i<ALEN(stk->dirichlet); i++) stk->dirichlet[i] = 0; /* Clear out any leftover values */
+        if (n < 3) stk->alldirichlet = dFALSE;                             /* @bug More work to determine independent of the mesh whether all the boundaries are Dirichlet */
+      }
+    }
   } err = PetscOptionsEnd();dCHK(err);
 
   switch (exact) {
@@ -191,10 +205,8 @@ static dErr StokesSetFromOptions(Stokes stk)
   err = dFSSetBlockSize(fsu,3);dCHK(err);
   err = dFSSetMesh(fsu,mesh,domain);dCHK(err);
   err = dFSSetDegree(fsu,jac,dtag);dCHK(err);
-  err = dFSRegisterBoundary(fsu,100,dFSBSTATUS_DIRICHLET,NULL,NULL);dCHK(err);
-  err = dFSRegisterBoundary(fsu,200,dFSBSTATUS_DIRICHLET,NULL,NULL);dCHK(err);
-  if (!stk->neumann300) {
-    err = dFSRegisterBoundary(fsu,300,dFSBSTATUS_DIRICHLET,NULL,NULL);dCHK(err);
+  for (dInt i=0; i<ALEN(stk->dirichlet) && stk->dirichlet[i]>0; i++) {
+    err = dFSRegisterBoundary(fsu,stk->dirichlet[i],dFSBSTATUS_DIRICHLET,NULL,NULL);dCHK(err);
   }
   err = PetscObjectSetOptionsPrefix((dObject)fsu,"u");dCHK(err);
   err = dFSSetFromOptions(fsu);dCHK(err);
@@ -1056,7 +1068,7 @@ int main(int argc,char *argv[])
     err = VecDestroy(&b);dCHK(err);
   }
 
-  if (!stk->neumann300) {                             /* Set null space */
+  if (stk->alldirichlet) {                             /* Set null space */
     KSP ksp;
     MatNullSpace matnull;
     err = StokesGetNullSpace(stk,&matnull);dCHK(err);