1 const char help[] = 2 "A test of H-div conforming discretizations on different cell types.\n"; 3 4 #include <petscdmplex.h> 5 #include <petscds.h> 6 #include <petscsnes.h> 7 #include <petscconvest.h> 8 #include <petscfe.h> 9 #include <petsc/private/petscfeimpl.h> 10 11 /* 12 * We are using the system 13 * 14 * \vec{u} = \vec{\hat{u}} 15 * p = \div{\vec{u}} in low degree approximation space 16 * d = \div{\vec{u}} - p == 0 in higher degree approximation space 17 * 18 * That is, we are using the field d to compute the error between \div{\vec{u}} 19 * computed in a space 1 degree higher than p and the value of p which is 20 * \div{u} computed in the low degree space. If H-div 21 * elements are implemented correctly then this should be identically zero since 22 * the divergence of a function in H(div) should be exactly representable in L_2 23 * by definition. 24 */ 25 static PetscErrorCode zero_func(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar *u,void *ctx) 26 { 27 PetscInt c; 28 for (c = 0; c < Nc; ++c) u[c] = 0; 29 return 0; 30 } 31 /* Linear Exact Functions 32 \vec{u} = \vec{x}; 33 p = dim; 34 */ 35 static PetscErrorCode linear_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar *u,void *ctx) 36 { 37 PetscInt c; 38 for (c = 0; c < Nc; ++c) u[c] = x[c]; 39 return 0; 40 } 41 static PetscErrorCode linear_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar *u,void *ctx) 42 { 43 u[0] = dim; 44 return 0; 45 } 46 47 /* Sinusoidal Exact Functions 48 * u_i = \sin{2*\pi*x_i} 49 * p = \Sum_{i=1}^{dim} 2*\pi*cos{2*\pi*x_i} 50 * */ 51 52 static PetscErrorCode sinusoid_u(PetscInt dim,PetscReal time,const PetscReal 53 x[],PetscInt Nc,PetscScalar *u,void *ctx) 54 { 55 PetscInt c; 56 for (c = 0; c< Nc; ++c) u[c] = PetscSinReal(2*PETSC_PI*x[c]); 57 return 0; 58 } 59 static PetscErrorCode sinusoid_p(PetscInt dim,PetscReal time,const PetscReal 60 x[],PetscInt Nc,PetscScalar *u,void *ctx) 61 { 62 PetscInt d; 63 u[0]=0; 64 for (d=0; d<dim; ++d) u[0] += 2*PETSC_PI*PetscCosReal(2*PETSC_PI*x[d]); 65 return 0; 66 } 67 68 /* Pointwise residual for u = u*. Need one of these for each possible u* */ 69 static void f0_v_linear(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 70 { 71 PetscInt i; 72 PetscScalar *u_rhs; 73 74 PetscCalloc1(dim,&u_rhs); 75 (void) linear_u(dim,t,x,dim,u_rhs,NULL); 76 for (i = 0; i < dim; ++i) f0[i] = u[uOff[0]+i]-u_rhs[i]; 77 PetscFree(u_rhs); 78 } 79 80 static void f0_v_sinusoid(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 81 { 82 PetscInt i; 83 PetscScalar *u_rhs; 84 85 PetscCalloc1(dim,&u_rhs); 86 (void) sinusoid_u(dim,t,x,dim,u_rhs,NULL); 87 for (i = 0; i < dim; ++i) f0[i] = u[uOff[0]+i]-u_rhs[i]; 88 PetscFree(u_rhs); 89 } 90 91 /* Residual function for enforcing p = \div{u}. */ 92 static void f0_q(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 93 { 94 PetscInt i; 95 PetscScalar divu; 96 97 divu = 0.; 98 for (i = 0; i< dim; ++i) divu += u_x[uOff_x[0]+i*dim+i]; 99 f0[0] = u[uOff[1]] - divu; 100 } 101 102 /* Residual function for p_err = \div{u} - p. */ 103 static void f0_w(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 104 { 105 PetscInt i; 106 PetscScalar divu; 107 108 divu = 0.; 109 for (i = 0; i < dim; ++i) divu += u_x[uOff_x[0] + i*dim +i]; 110 f0[0] = u[uOff[2]] - u[uOff[1]] + divu; 111 } 112 113 /* Boundary residual for the embedding system. Need one for each form of 114 * solution. These enforce u = \hat{u} at the boundary. */ 115 static void f0_bd_u_sinusoid(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],const PetscReal n[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 116 { 117 PetscInt d; 118 PetscScalar *u_rhs; 119 PetscCalloc1(dim,&u_rhs); 120 (void) sinusoid_u(dim,t,x,dim,u_rhs,NULL); 121 122 for (d=0; d<dim; ++d) f0[d] = u_rhs[d]; 123 PetscFree(u_rhs); 124 125 } 126 127 static void f0_bd_u_linear(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],const PetscReal n[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 128 { 129 PetscInt d; 130 PetscScalar *u_rhs; 131 PetscCalloc1(dim,&u_rhs); 132 (void) linear_u(dim,t,x,dim,u_rhs,NULL); 133 134 for (d=0; d<dim; ++d) f0[d] = u_rhs[d]; 135 PetscFree(u_rhs); 136 } 137 /* Jacobian functions. For the following, v is the test function associated with 138 * u, q the test function associated with p, and w the test function associated 139 * with d. */ 140 /* <v, u> */ 141 static void g0_vu(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g0[]) 142 { 143 PetscInt c; 144 for (c = 0; c < dim; ++c) g0[c * dim + c] = 1.0; 145 } 146 147 /* <q, p> */ 148 static void g0_qp(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g0[]) 149 { 150 PetscInt d; 151 for (d=0; d< dim; ++d) g0[d * dim + d] = 1.0; 152 } 153 154 /* -<q, \div{u}> For the embedded system. This is different from the method of 155 * manufactured solution because instead of computing <q,\div{u}> - <q,f> we 156 * need <q,p> - <q,\div{u}.*/ 157 static void g1_qu(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g1[]) 158 { 159 PetscInt d; 160 for (d = 0; d < dim; ++d) g1[d * dim + d] = -1.0; 161 } 162 163 /* <w, p> This is only used by the embedded system. Where we need to compute 164 * <w,d> - <w,p> + <w, \div{u}>*/ 165 static void g0_wp(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g0[]) 166 { 167 PetscInt d; 168 for (d=0; d< dim; ++d) g0[d * dim + d] = -1.0; 169 } 170 171 /* <w, d> */ 172 static void g0_wd(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g0[]) 173 { 174 PetscInt c; 175 for (c = 0; c < dim; ++c) g0[c*dim+c] = 1.0; 176 } 177 178 /* <w, \div{u}> for the embedded system. */ 179 static void g1_wu(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g1[]) 180 { 181 PetscInt d; 182 for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; 183 } 184 185 /* Enum and string array for selecting mesh perturbation options */ 186 typedef enum { NONE = 0,PERTURB = 1,SKEW = 2,SKEW_PERTURB = 3 } Transform; 187 const char* const TransformTypes[] = {"none","perturb","skew","skew_perturb","Perturbation","",NULL}; 188 189 /* Enum and string array for selecting the form of the exact solution*/ 190 typedef enum 191 { LINEAR = 0,SINUSOIDAL = 1 } Solution; 192 const char* const SolutionTypes[] = {"linear","sinusoidal","Solution","",NULL}; 193 194 typedef struct 195 { 196 PetscBool simplex; 197 PetscInt dim; 198 Transform mesh_transform; 199 Solution sol_form; 200 } UserCtx; 201 202 /* Process command line options and initialize the UserCtx struct */ 203 static PetscErrorCode ProcessOptions(MPI_Comm comm,UserCtx * user) 204 { 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 /* Default to 2D, unperturbed triangle mesh and Linear solution.*/ 209 user->simplex = PETSC_TRUE; 210 user->dim = 2; 211 user->mesh_transform = NONE; 212 user->sol_form = LINEAR; 213 214 ierr = PetscOptionsBegin(comm,"","H-div Test Options","DMPLEX");CHKERRQ(ierr); 215 ierr = PetscOptionsBool("-simplex","Whether to use simplices (true) or tensor-product (false) cells in " "the mesh","ex39.c",user->simplex,&user->simplex,NULL);CHKERRQ(ierr); 216 ierr = PetscOptionsInt("-dim","Number of solution dimensions","ex39.c",user->dim,&user->dim,NULL);CHKERRQ(ierr); 217 ierr = PetscOptionsEnum("-mesh_transform","Method used to perturb the mesh vertices. Options are skew, perturb, skew_perturb,or none","ex39.c",TransformTypes,(PetscEnum) user->mesh_transform,(PetscEnum*) &user->mesh_transform,NULL);CHKERRQ(ierr); 218 ierr = PetscOptionsEnum("-sol_form","Form of the exact solution. Options are Linear or Sinusoidal","ex39.c",SolutionTypes,(PetscEnum) user->sol_form,(PetscEnum*) &user->sol_form,NULL);CHKERRQ(ierr); 219 ierr = PetscOptionsEnd();CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 223 /* Perturb the position of each mesh vertex by a small amount.*/ 224 static PetscErrorCode PerturbMesh(DM *mesh,PetscScalar *coordVals,PetscInt npoints,PetscInt dim) 225 { 226 PetscInt i,j,k; 227 PetscErrorCode ierr; 228 PetscReal minCoords[3],maxCoords[3],maxPert[3],randVal,amp; 229 PetscRandom ran; 230 231 PetscFunctionBegin; 232 ierr = DMGetCoordinateDim(*mesh,&dim);CHKERRQ(ierr); 233 ierr = DMGetLocalBoundingBox(*mesh,minCoords,maxCoords);CHKERRQ(ierr); 234 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&ran);CHKERRQ(ierr); 235 236 /* Compute something approximately equal to half an edge length. This is the 237 * most we can perturb points and gaurantee that there won't be any topology 238 * issues. */ 239 for (k = 0; k < dim; ++k) maxPert[k] = 0.025 * (maxCoords[k] - minCoords[k]) / (PetscPowReal(npoints,1. / dim) - 1); 240 /* For each mesh vertex */ 241 for (i = 0; i < npoints; ++i) { 242 /* For each coordinate of the vertex */ 243 for (j = 0; j < dim; ++j) { 244 /* Generate a random amplitude in [-0.5*maxPert, 0.5*maxPert] */ 245 ierr = PetscRandomGetValueReal(ran,&randVal);CHKERRQ(ierr); 246 amp = maxPert[j] * (randVal - 0.5); 247 /* Add the perturbation to the vertex*/ 248 coordVals[dim * i + j] += amp; 249 } 250 } 251 252 PetscRandomDestroy(&ran); 253 PetscFunctionReturn(0); 254 } 255 256 /* Apply a global skew transformation to the mesh. */ 257 static PetscErrorCode SkewMesh(DM * mesh,PetscScalar * coordVals,PetscInt npoints,PetscInt dim) 258 { 259 PetscInt i,j,k,l; 260 PetscErrorCode ierr; 261 PetscScalar * transMat; 262 PetscScalar tmpcoord[3]; 263 PetscRandom ran; 264 PetscReal randVal; 265 266 PetscFunctionBegin; 267 ierr = PetscCalloc1(dim * dim,&transMat);CHKERRQ(ierr); 268 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&ran);CHKERRQ(ierr); 269 270 /* Make a matrix representing a skew transformation */ 271 for (i = 0; i < dim; ++i) { 272 for (j = 0; j < dim; ++j) { 273 ierr = PetscRandomGetValueReal(ran,&randVal);CHKERRQ(ierr); 274 if (i == j) transMat[i * dim + j] = 1.; 275 else if (j < i) transMat[i * dim + j] = 2 * (j + i)*randVal; 276 else transMat[i * dim + j] = 0; 277 } 278 } 279 280 /* Multiply each coordinate vector by our tranformation.*/ 281 for (i = 0; i < npoints; ++i) { 282 for (j = 0; j < dim; ++j) { 283 tmpcoord[j] = 0; 284 for (k = 0; k < dim; ++k) tmpcoord[j] += coordVals[dim * i + k] * transMat[dim * k + j]; 285 } 286 for (l = 0; l < dim; ++l) coordVals[dim * i + l] = tmpcoord[l]; 287 } 288 ierr = PetscFree(transMat);CHKERRQ(ierr); 289 ierr = PetscRandomDestroy(&ran);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 /* Accesses the mesh coordinate array and performs the transformation operations 294 * specified by the user options */ 295 static PetscErrorCode TransformMesh(UserCtx * user,DM * mesh) 296 { 297 PetscErrorCode ierr; 298 PetscInt dim,npoints; 299 PetscScalar * coordVals; 300 Vec coords; 301 302 PetscFunctionBegin; 303 ierr = DMGetCoordinates(*mesh,&coords);CHKERRQ(ierr); 304 ierr = VecGetArray(coords,&coordVals);CHKERRQ(ierr); 305 ierr = VecGetLocalSize(coords,&npoints);CHKERRQ(ierr); 306 ierr = DMGetCoordinateDim(*mesh,&dim);CHKERRQ(ierr); 307 npoints = npoints / dim; 308 309 switch (user->mesh_transform) { 310 case PERTURB: 311 ierr = PerturbMesh(mesh,coordVals,npoints,dim);CHKERRQ(ierr); 312 break; 313 case SKEW: 314 ierr = SkewMesh(mesh,coordVals,npoints,dim);CHKERRQ(ierr); 315 break; 316 case SKEW_PERTURB: 317 ierr = SkewMesh(mesh,coordVals,npoints,dim);CHKERRQ(ierr); 318 ierr = PerturbMesh(mesh,coordVals,npoints,dim);CHKERRQ(ierr); 319 break; 320 default: 321 PetscFunctionReturn(-1); 322 } 323 ierr = VecRestoreArray(coords,&coordVals);CHKERRQ(ierr); 324 ierr = DMSetCoordinates(*mesh,coords);CHKERRQ(ierr); 325 PetscFunctionReturn(0); 326 } 327 328 static PetscErrorCode CreateMesh(MPI_Comm comm,UserCtx * user,DM * mesh) 329 { 330 PetscErrorCode ierr; 331 DMLabel label; 332 const char *name = "marker"; 333 DM dmDist = NULL; 334 PetscPartitioner part; 335 336 PetscFunctionBegin; 337 /* Create box mesh from user parameters */ 338 ierr = DMPlexCreateBoxMesh(comm,user->dim,user->simplex,NULL,NULL,NULL,NULL,PETSC_TRUE,mesh);CHKERRQ(ierr); 339 340 /* Make sure the mesh gets properly distributed if running in parallel */ 341 ierr = DMPlexGetPartitioner(*mesh,&part);CHKERRQ(ierr); 342 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 343 ierr = DMPlexDistribute(*mesh,0,NULL,&dmDist);CHKERRQ(ierr); 344 if (dmDist) { 345 ierr = DMDestroy(mesh);CHKERRQ(ierr); 346 *mesh = dmDist; 347 } 348 349 /* Mark the boundaries, we will need this later when setting up the system of 350 * equations */ 351 ierr = DMCreateLabel(*mesh,name);CHKERRQ(ierr); 352 ierr = DMGetLabel(*mesh,name,&label);CHKERRQ(ierr); 353 ierr = DMPlexMarkBoundaryFaces(*mesh,1,label);CHKERRQ(ierr); 354 ierr = DMPlexLabelComplete(*mesh,label);CHKERRQ(ierr); 355 ierr = DMLocalizeCoordinates(*mesh);CHKERRQ(ierr); 356 ierr = PetscObjectSetName((PetscObject) *mesh,"Mesh");CHKERRQ(ierr); 357 358 /* Perform any mesh transformations if specified by user */ 359 if (user->mesh_transform != NONE) { 360 ierr = TransformMesh(user,mesh);CHKERRQ(ierr); 361 } 362 363 /* Get any other mesh options from the command line */ 364 ierr = DMSetApplicationContext(*mesh,user);CHKERRQ(ierr); 365 ierr = DMSetFromOptions(*mesh);CHKERRQ(ierr); 366 ierr = DMViewFromOptions(*mesh,NULL,"-dm_view");CHKERRQ(ierr); 367 368 ierr = DMDestroy(&dmDist);CHKERRQ(ierr); 369 PetscFunctionReturn(0); 370 } 371 372 /* Setup the system of equations that we wish to solve */ 373 static PetscErrorCode SetupProblem(DM dm,UserCtx * user) 374 { 375 PetscDS prob; 376 DMLabel label; 377 PetscErrorCode ierr; 378 const PetscInt id=1; 379 380 PetscFunctionBegin; 381 ierr = DMGetDS(dm,&prob);CHKERRQ(ierr); 382 /* All of these are independent of the user's choice of solution */ 383 ierr = PetscDSSetResidual(prob,1,f0_q,NULL);CHKERRQ(ierr); 384 ierr = PetscDSSetResidual(prob,2,f0_w,NULL);CHKERRQ(ierr); 385 ierr = PetscDSSetJacobian(prob,0,0,g0_vu,NULL,NULL,NULL);CHKERRQ(ierr); 386 ierr = PetscDSSetJacobian(prob,1,0,NULL,g1_qu,NULL,NULL);CHKERRQ(ierr); 387 ierr = PetscDSSetJacobian(prob,1,1,g0_qp,NULL,NULL,NULL);CHKERRQ(ierr); 388 ierr = PetscDSSetJacobian(prob,2,0,NULL,g1_wu,NULL,NULL);CHKERRQ(ierr); 389 ierr = PetscDSSetJacobian(prob,2,1,g0_wp,NULL,NULL,NULL);CHKERRQ(ierr); 390 ierr = PetscDSSetJacobian(prob,2,2,g0_wd,NULL,NULL,NULL);CHKERRQ(ierr); 391 392 /* Field 2 is the error between \div{u} and pressure in a higher dimensional 393 * space. If all is right this should be machine zero. */ 394 ierr = PetscDSSetExactSolution(prob,2,zero_func,NULL);CHKERRQ(ierr); 395 396 switch (user->sol_form) { 397 case LINEAR: 398 ierr = PetscDSSetResidual(prob,0,f0_v_linear,NULL);CHKERRQ(ierr); 399 ierr = PetscDSSetBdResidual(prob,0,f0_bd_u_linear,NULL);CHKERRQ(ierr); 400 ierr = PetscDSSetExactSolution(prob,0,linear_u,NULL);CHKERRQ(ierr); 401 ierr = PetscDSSetExactSolution(prob,1,linear_p,NULL);CHKERRQ(ierr); 402 break; 403 case SINUSOIDAL: 404 ierr = PetscDSSetResidual(prob,0,f0_v_sinusoid,NULL);CHKERRQ(ierr); 405 ierr = PetscDSSetBdResidual(prob,0,f0_bd_u_sinusoid,NULL);CHKERRQ(ierr); 406 ierr = PetscDSSetExactSolution(prob,0,sinusoid_u,NULL);CHKERRQ(ierr); 407 ierr = PetscDSSetExactSolution(prob,1,sinusoid_p,NULL);CHKERRQ(ierr); 408 break; 409 default: 410 PetscFunctionReturn(-1); 411 } 412 413 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 414 ierr = PetscDSAddBoundary(prob,DM_BC_NATURAL,"Boundary Integral",label,1,&id,0,0,NULL,(void (*)(void))NULL,NULL,user,NULL);CHKERRQ(ierr); 415 PetscFunctionReturn(0); 416 } 417 418 /* Create the finite element spaces we will use for this system */ 419 static PetscErrorCode SetupDiscretization(DM mesh,PetscErrorCode (*setup)(DM,UserCtx*),UserCtx *user) 420 { 421 DM cdm = mesh; 422 PetscFE fevel,fepres,fedivErr; 423 const PetscInt dim = user->dim; 424 PetscErrorCode ierr; 425 426 427 PetscFunctionBegin; 428 /* Create FE objects and give them names so that options can be set from 429 * command line */ 430 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) mesh),dim,dim,user->simplex,"velocity_",-1,&fevel);CHKERRQ(ierr); 431 ierr = PetscObjectSetName((PetscObject) fevel,"velocity");CHKERRQ(ierr); 432 433 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) mesh),dim,1,user->simplex,"pressure_",-1,&fepres);CHKERRQ(ierr); 434 ierr = PetscObjectSetName((PetscObject) fepres,"pressure");CHKERRQ(ierr); 435 436 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) 437 mesh),dim,1,user->simplex,"divErr_",-1,&fedivErr);CHKERRQ(ierr); 438 ierr = PetscObjectSetName((PetscObject) fedivErr,"divErr");CHKERRQ(ierr); 439 440 ierr = PetscFECopyQuadrature(fevel,fepres);CHKERRQ(ierr); 441 ierr = PetscFECopyQuadrature(fevel,fedivErr);CHKERRQ(ierr); 442 443 /* Associate the FE objects with the mesh and setup the system */ 444 ierr = DMSetField(mesh,0,NULL,(PetscObject) fevel);CHKERRQ(ierr); 445 ierr = DMSetField(mesh,1,NULL,(PetscObject) fepres);CHKERRQ(ierr); 446 ierr = DMSetField(mesh,2,NULL,(PetscObject) fedivErr);CHKERRQ(ierr); 447 ierr = DMCreateDS(mesh);CHKERRQ(ierr); 448 ierr = (*setup)(mesh,user);CHKERRQ(ierr); 449 450 while (cdm) { 451 ierr = DMCopyDisc(mesh,cdm);CHKERRQ(ierr); 452 ierr = DMGetCoarseDM(cdm,&cdm);CHKERRQ(ierr); 453 } 454 455 /* The Mesh now owns the fields, so we can destroy the FEs created in this 456 * function */ 457 ierr = PetscFEDestroy(&fevel);CHKERRQ(ierr); 458 ierr = PetscFEDestroy(&fepres);CHKERRQ(ierr); 459 ierr = PetscFEDestroy(&fedivErr);CHKERRQ(ierr); 460 ierr = DMDestroy(&cdm);CHKERRQ(ierr); 461 PetscFunctionReturn(0); 462 } 463 464 465 466 int main(int argc,char **argv) 467 { 468 PetscInt i; 469 UserCtx user; 470 DM mesh; 471 SNES snes; 472 Vec computed,divErr; 473 PetscReal divErrNorm; 474 PetscErrorCode ierr; 475 IS * fieldIS; 476 PetscBool exampleSuccess = PETSC_FALSE; 477 const PetscReal errTol = 10. * PETSC_SMALL; 478 479 char stdFormat[] = "L2 Norm of the Divergence Error is: %g\n H(div) elements working correctly: %s\n"; 480 481 /* Initialize PETSc */ 482 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 483 ierr = ProcessOptions(PETSC_COMM_WORLD,&user);CHKERRQ(ierr); 484 485 /* Set up the system, we need to create a solver and a mesh and then assign 486 * the correct spaces into the mesh*/ 487 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 488 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&mesh);CHKERRQ(ierr); 489 ierr = SNESSetDM(snes,mesh);CHKERRQ(ierr); 490 ierr = SetupDiscretization(mesh,SetupProblem,&user);CHKERRQ(ierr); 491 ierr = DMPlexSetSNESLocalFEM(mesh,&user,&user,&user);CHKERRQ(ierr); 492 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 493 494 /* Grab field IS so that we can view the solution by field */ 495 ierr = DMCreateFieldIS(mesh,NULL,NULL,&fieldIS);CHKERRQ(ierr); 496 497 /* Create a vector to store the SNES solution, solve the system and grab the 498 * solution from SNES */ 499 ierr = DMCreateGlobalVector(mesh,&computed);CHKERRQ(ierr); 500 ierr = PetscObjectSetName((PetscObject) computed,"computedSolution");CHKERRQ(ierr); 501 ierr = VecSet(computed,0.0);CHKERRQ(ierr); 502 ierr = SNESSolve(snes,NULL,computed);CHKERRQ(ierr); 503 ierr = SNESGetSolution(snes,&computed);CHKERRQ(ierr); 504 ierr = VecViewFromOptions(computed,NULL,"-computedSolution_view");CHKERRQ(ierr); 505 506 /* Now we pull out the portion of the vector corresponding to the 3rd field 507 * which is the error between \div{u} computed in a higher dimensional space 508 * and p = \div{u} computed in a low dimension space. We report the L2 norm of 509 * this vector which should be zero if the H(div) spaces are implemented 510 * correctly. */ 511 ierr = VecGetSubVector(computed,fieldIS[2],&divErr);CHKERRQ(ierr); 512 ierr = VecNorm(divErr,NORM_2,&divErrNorm);CHKERRQ(ierr); 513 ierr = VecRestoreSubVector(computed,fieldIS[2],&divErr);CHKERRQ(ierr); 514 exampleSuccess = (PetscBool)(divErrNorm <= errTol); 515 516 517 ierr = PetscPrintf(PETSC_COMM_WORLD,stdFormat,divErrNorm,exampleSuccess ? "true" : "false");CHKERRQ(ierr); 518 519 520 /* Tear down */ 521 ierr = VecDestroy(&divErr);CHKERRQ(ierr); 522 ierr = VecDestroy(&computed);CHKERRQ(ierr); 523 for (i = 0; i < 3; ++i) { 524 ierr = ISDestroy(&fieldIS[i]);CHKERRQ(ierr); 525 } 526 ierr = PetscFree(fieldIS);CHKERRQ(ierr); 527 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 528 ierr = DMDestroy(&mesh);CHKERRQ(ierr); 529 ierr = PetscFinalize(); 530 return ierr; 531 } 532 533 /*TEST 534 testset: 535 suffix: 2d_bdm 536 requires: triangle 537 args: -dim 2 \ 538 -simplex true \ 539 -velocity_petscfe_default_quadrature_order 1 \ 540 -velocity_petscspace_degree 1 \ 541 -velocity_petscdualspace_type bdm \ 542 -divErr_petscspace_degree 1 \ 543 -divErr_petscdualspace_lagrange_continuity false \ 544 -dm_refine 0 \ 545 -snes_error_if_not_converged \ 546 -ksp_rtol 1e-10 \ 547 -ksp_error_if_not_converged \ 548 -pc_type fieldsplit\ 549 -pc_fieldsplit_detect_saddle_point\ 550 -pc_fieldsplit_type schur\ 551 -pc_fieldsplit_schur_precondition full 552 test: 553 suffix: linear 554 args: -sol_form linear -mesh_transform none 555 test: 556 suffix: sinusoidal 557 args: -sol_form sinusoidal -mesh_transform none 558 test: 559 suffix: sinusoidal_skew 560 args: -sol_form sinusoidal -mesh_transform skew 561 test: 562 suffix: sinusoidal_perturb 563 args: -sol_form sinusoidal -mesh_transform perturb 564 test: 565 suffix: sinusoidal_skew_perturb 566 args: -sol_form sinusoidal -mesh_transform skew_perturb 567 568 testset: 569 TODO: broken 570 suffix: 2d_bdmq 571 args: -dim 2 \ 572 -simplex false \ 573 -velocity_petscspace_degree 1 \ 574 -velocity_petscdualspace_type bdm \ 575 -velocity_petscdualspace_lagrange_tensor 1 \ 576 -divErr_petscspace_degree 1 \ 577 -divErr_petscdualspace_lagrange_continuity false \ 578 -dm_refine 0 \ 579 -snes_error_if_not_converged \ 580 -ksp_rtol 1e-10 \ 581 -ksp_error_if_not_converged \ 582 -pc_type fieldsplit\ 583 -pc_fieldsplit_detect_saddle_point\ 584 -pc_fieldsplit_type schur\ 585 -pc_fieldsplit_schur_precondition full 586 test: 587 suffix: linear 588 args: -sol_form linear -mesh_transform none 589 test: 590 suffix: sinusoidal 591 args: -sol_form sinusoidal -mesh_transform none 592 test: 593 suffix: sinusoidal_skew 594 args: -sol_form sinusoidal -mesh_transform skew 595 test: 596 suffix: sinusoidal_perturb 597 args: -sol_form sinusoidal -mesh_transform perturb 598 test: 599 suffix: sinusoidal_skew_perturb 600 args: -sol_form sinusoidal -mesh_transform skew_perturb 601 602 testset: 603 suffix: 3d_bdm 604 requires: ctetgen 605 args: -dim 3 \ 606 -simplex true \ 607 -velocity_petscspace_degree 1 \ 608 -velocity_petscdualspace_type bdm \ 609 -divErr_petscspace_degree 1 \ 610 -divErr_petscdualspace_lagrange_continuity false \ 611 -dm_refine 0 \ 612 -snes_error_if_not_converged \ 613 -ksp_rtol 1e-10 \ 614 -ksp_error_if_not_converged \ 615 -pc_type fieldsplit \ 616 -pc_fieldsplit_detect_saddle_point \ 617 -pc_fieldsplit_type schur \ 618 -pc_fieldsplit_schur_precondition full 619 test: 620 suffix: linear 621 args: -sol_form linear -mesh_transform none 622 test: 623 suffix: sinusoidal 624 args: -sol_form sinusoidal -mesh_transform none 625 test: 626 suffix: sinusoidal_skew 627 args: -sol_form sinusoidal -mesh_transform skew 628 test: 629 suffix: sinusoidal_perturb 630 args: -sol_form sinusoidal -mesh_transform perturb 631 test: 632 suffix: sinusoidal_skew_perturb 633 args: -sol_form sinusoidal -mesh_transform skew_perturb 634 635 testset: 636 TODO: broken 637 suffix: 3d_bdmq 638 requires: ctetgen 639 args: -dim 3 \ 640 -simplex false \ 641 -velocity_petscspace_degree 1 \ 642 -velocity_petscdualspace_type bdm \ 643 -velocity_petscdualspace_lagrange_tensor 1 \ 644 -divErr_petscspace_degree 1 \ 645 -divErr_petscdualspace_lagrange_continuity false \ 646 -dm_refine 0 \ 647 -snes_error_if_not_converged \ 648 -ksp_rtol 1e-10 \ 649 -ksp_error_if_not_converged \ 650 -pc_type fieldsplit \ 651 -pc_fieldsplit_detect_saddle_point \ 652 -pc_fieldsplit_type schur \ 653 -pc_fieldsplit_schur_precondition full 654 test: 655 suffix: linear 656 args: -sol_form linear -mesh_transform none 657 test: 658 suffix: sinusoidal 659 args: -sol_form sinusoidal -mesh_transform none 660 test: 661 suffix: sinusoidal_skew 662 args: -sol_form sinusoidal -mesh_transform skew 663 test: 664 suffix: sinusoidal_perturb 665 args: -sol_form sinusoidal -mesh_transform perturb 666 test: 667 suffix: sinusoidal_skew_perturb 668 args: -sol_form sinusoidal -mesh_transform skew_perturb 669 670 test: 671 suffix: quad_rt_0 672 args: -dim 2 -simplex false -mesh_transform skew \ 673 -divErr_petscspace_degree 1 \ 674 -divErr_petscdualspace_lagrange_continuity false \ 675 -dm_refine 0 \ 676 -snes_error_if_not_converged \ 677 -ksp_rtol 1e-10 \ 678 -ksp_error_if_not_converged \ 679 -pc_type fieldsplit\ 680 -pc_fieldsplit_detect_saddle_point\ 681 -pc_fieldsplit_type schur\ 682 -pc_fieldsplit_schur_precondition full \ 683 -velocity_petscfe_default_quadrature_order 1 \ 684 -velocity_petscspace_type sum \ 685 -velocity_petscspace_variables 2 \ 686 -velocity_petscspace_components 2 \ 687 -velocity_petscspace_sum_spaces 2 \ 688 -velocity_petscspace_sum_concatenate true \ 689 -velocity_subspace0_petscspace_variables 2 \ 690 -velocity_subspace0_petscspace_type tensor \ 691 -velocity_subspace0_petscspace_tensor_spaces 2 \ 692 -velocity_subspace0_petscspace_tensor_uniform false \ 693 -velocity_subspace0_subspace_0_petscspace_degree 1 \ 694 -velocity_subspace0_subspace_1_petscspace_degree 0 \ 695 -velocity_subspace1_petscspace_variables 2 \ 696 -velocity_subspace1_petscspace_type tensor \ 697 -velocity_subspace1_petscspace_tensor_spaces 2 \ 698 -velocity_subspace1_petscspace_tensor_uniform false \ 699 -velocity_subspace1_subspace_0_petscspace_degree 0 \ 700 -velocity_subspace1_subspace_1_petscspace_degree 1 \ 701 -velocity_petscdualspace_form_degree -1 \ 702 -velocity_petscdualspace_order 1 \ 703 -velocity_petscdualspace_lagrange_trimmed true 704 TEST*/ 705