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 PetscErrorCode ierr; 377 const PetscInt id=1; 378 379 PetscFunctionBegin; 380 ierr = DMGetDS(dm,&prob);CHKERRQ(ierr); 381 /* All of these are independent of the user's choice of solution */ 382 ierr = PetscDSSetResidual(prob,1,f0_q,NULL);CHKERRQ(ierr); 383 ierr = PetscDSSetResidual(prob,2,f0_w,NULL);CHKERRQ(ierr); 384 ierr = PetscDSSetJacobian(prob,0,0,g0_vu,NULL,NULL,NULL);CHKERRQ(ierr); 385 ierr = PetscDSSetJacobian(prob,1,0,NULL,g1_qu,NULL,NULL);CHKERRQ(ierr); 386 ierr = PetscDSSetJacobian(prob,1,1,g0_qp,NULL,NULL,NULL);CHKERRQ(ierr); 387 ierr = PetscDSSetJacobian(prob,2,0,NULL,g1_wu,NULL,NULL);CHKERRQ(ierr); 388 ierr = PetscDSSetJacobian(prob,2,1,g0_wp,NULL,NULL,NULL);CHKERRQ(ierr); 389 ierr = PetscDSSetJacobian(prob,2,2,g0_wd,NULL,NULL,NULL);CHKERRQ(ierr); 390 391 /* Field 2 is the error between \div{u} and pressure in a higher dimensional 392 * space. If all is right this should be machine zero. */ 393 ierr = PetscDSSetExactSolution(prob,2,zero_func,NULL);CHKERRQ(ierr); 394 395 switch (user->sol_form) { 396 case LINEAR: 397 ierr = PetscDSSetResidual(prob,0,f0_v_linear,NULL);CHKERRQ(ierr); 398 ierr = PetscDSSetBdResidual(prob,0,f0_bd_u_linear,NULL);CHKERRQ(ierr); 399 ierr = PetscDSSetExactSolution(prob,0,linear_u,NULL);CHKERRQ(ierr); 400 ierr = PetscDSSetExactSolution(prob,1,linear_p,NULL);CHKERRQ(ierr); 401 break; 402 case SINUSOIDAL: 403 ierr = PetscDSSetResidual(prob,0,f0_v_sinusoid,NULL);CHKERRQ(ierr); 404 ierr = PetscDSSetBdResidual(prob,0,f0_bd_u_sinusoid,NULL);CHKERRQ(ierr); 405 ierr = PetscDSSetExactSolution(prob,0,sinusoid_u,NULL);CHKERRQ(ierr); 406 ierr = PetscDSSetExactSolution(prob,1,sinusoid_p,NULL);CHKERRQ(ierr); 407 break; 408 default: 409 PetscFunctionReturn(-1); 410 } 411 412 ierr = PetscDSAddBoundary(prob,DM_BC_NATURAL,"Boundary Integral","marker",0,0,NULL,(void (*)(void))NULL,NULL,1,&id,user);CHKERRQ(ierr); 413 PetscFunctionReturn(0); 414 } 415 416 /* Create the finite element spaces we will use for this system */ 417 static PetscErrorCode SetupDiscretization(DM mesh,PetscErrorCode (*setup)(DM,UserCtx*),UserCtx *user) 418 { 419 DM cdm = mesh; 420 PetscFE fevel,fepres,fedivErr; 421 const PetscInt dim = user->dim; 422 PetscErrorCode ierr; 423 424 425 PetscFunctionBegin; 426 /* Create FE objects and give them names so that options can be set from 427 * command line */ 428 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) mesh),dim,dim,user->simplex,"velocity_",-1,&fevel);CHKERRQ(ierr); 429 ierr = PetscObjectSetName((PetscObject) fevel,"velocity");CHKERRQ(ierr); 430 431 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) mesh),dim,1,user->simplex,"pressure_",-1,&fepres);CHKERRQ(ierr); 432 ierr = PetscObjectSetName((PetscObject) fepres,"pressure");CHKERRQ(ierr); 433 434 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) 435 mesh),dim,1,user->simplex,"divErr_",-1,&fedivErr);CHKERRQ(ierr); 436 ierr = PetscObjectSetName((PetscObject) fedivErr,"divErr");CHKERRQ(ierr); 437 438 ierr = PetscFECopyQuadrature(fevel,fepres);CHKERRQ(ierr); 439 ierr = PetscFECopyQuadrature(fevel,fedivErr);CHKERRQ(ierr); 440 441 /* Associate the FE objects with the mesh and setup the system */ 442 ierr = DMSetField(mesh,0,NULL,(PetscObject) fevel);CHKERRQ(ierr); 443 ierr = DMSetField(mesh,1,NULL,(PetscObject) fepres);CHKERRQ(ierr); 444 ierr = DMSetField(mesh,2,NULL,(PetscObject) fedivErr);CHKERRQ(ierr); 445 ierr = DMCreateDS(mesh);CHKERRQ(ierr); 446 ierr = (*setup)(mesh,user);CHKERRQ(ierr); 447 448 while (cdm) { 449 ierr = DMCopyDisc(mesh,cdm);CHKERRQ(ierr); 450 ierr = DMGetCoarseDM(cdm,&cdm);CHKERRQ(ierr); 451 } 452 453 /* The Mesh now owns the fields, so we can destroy the FEs created in this 454 * function */ 455 ierr = PetscFEDestroy(&fevel);CHKERRQ(ierr); 456 ierr = PetscFEDestroy(&fepres);CHKERRQ(ierr); 457 ierr = PetscFEDestroy(&fedivErr);CHKERRQ(ierr); 458 ierr = DMDestroy(&cdm);CHKERRQ(ierr); 459 PetscFunctionReturn(0); 460 } 461 462 463 464 int main(int argc,char **argv) 465 { 466 PetscInt i; 467 UserCtx user; 468 DM mesh; 469 SNES snes; 470 Vec computed,divErr; 471 PetscReal divErrNorm; 472 PetscErrorCode ierr; 473 IS * fieldIS; 474 PetscBool exampleSuccess = PETSC_FALSE; 475 const PetscReal errTol = 10. * PETSC_SMALL; 476 477 char stdFormat[] = "L2 Norm of the Divergence Error is: %g\n H(div) elements working correctly: %s\n"; 478 479 /* Initialize PETSc */ 480 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 481 ierr = ProcessOptions(PETSC_COMM_WORLD,&user);CHKERRQ(ierr); 482 483 /* Set up the system, we need to create a solver and a mesh and then assign 484 * the correct spaces into the mesh*/ 485 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 486 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&mesh);CHKERRQ(ierr); 487 ierr = SNESSetDM(snes,mesh);CHKERRQ(ierr); 488 ierr = SetupDiscretization(mesh,SetupProblem,&user);CHKERRQ(ierr); 489 ierr = DMPlexSetSNESLocalFEM(mesh,&user,&user,&user);CHKERRQ(ierr); 490 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 491 492 /* Grab field IS so that we can view the solution by field */ 493 ierr = DMCreateFieldIS(mesh,NULL,NULL,&fieldIS);CHKERRQ(ierr); 494 495 /* Create a vector to store the SNES solution, solve the system and grab the 496 * solution from SNES */ 497 ierr = DMCreateGlobalVector(mesh,&computed);CHKERRQ(ierr); 498 ierr = PetscObjectSetName((PetscObject) computed,"computedSolution");CHKERRQ(ierr); 499 ierr = VecSet(computed,0.0);CHKERRQ(ierr); 500 ierr = SNESSolve(snes,NULL,computed);CHKERRQ(ierr); 501 ierr = SNESGetSolution(snes,&computed);CHKERRQ(ierr); 502 ierr = VecViewFromOptions(computed,NULL,"-computedSolution_view");CHKERRQ(ierr); 503 504 /* Now we pull out the portion of the vector corresponding to the 3rd field 505 * which is the error between \div{u} computed in a higher dimensional space 506 * and p = \div{u} computed in a low dimension space. We report the L2 norm of 507 * this vector which should be zero if the H(div) spaces are implemented 508 * correctly. */ 509 ierr = VecGetSubVector(computed,fieldIS[2],&divErr);CHKERRQ(ierr); 510 ierr = VecNorm(divErr,NORM_2,&divErrNorm); CHKERRQ(ierr); 511 ierr = VecRestoreSubVector(computed,fieldIS[2],&divErr);CHKERRQ(ierr); 512 exampleSuccess = (PetscBool)(divErrNorm <= errTol); 513 514 515 ierr = PetscPrintf(PETSC_COMM_WORLD,stdFormat,divErrNorm,exampleSuccess ? "true" : "false");CHKERRQ(ierr); 516 517 518 /* Tear down */ 519 ierr = VecDestroy(&divErr);CHKERRQ(ierr); 520 ierr = VecDestroy(&computed);CHKERRQ(ierr); 521 for (i = 0; i < 3; ++i) { 522 ierr = ISDestroy(&fieldIS[i]);CHKERRQ(ierr); 523 } 524 ierr = PetscFree(fieldIS);CHKERRQ(ierr); 525 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 526 ierr = DMDestroy(&mesh);CHKERRQ(ierr); 527 ierr = PetscFinalize(); 528 return ierr; 529 } 530 531 /*TEST 532 testset: 533 suffix: 2d_bdm 534 requires: triangle 535 args: -dim 2 \ 536 -simplex true \ 537 -velocity_petscfe_default_quadrature_order 1 \ 538 -velocity_petscspace_degree 1 \ 539 -velocity_petscdualspace_type bdm \ 540 -divErr_petscspace_degree 1 \ 541 -divErr_petscdualspace_lagrange_continuity false \ 542 -dm_refine 0 \ 543 -snes_error_if_not_converged \ 544 -ksp_rtol 1e-10 \ 545 -ksp_error_if_not_converged \ 546 -pc_type fieldsplit\ 547 -pc_fieldsplit_detect_saddle_point\ 548 -pc_fieldsplit_type schur\ 549 -pc_fieldsplit_schur_precondition full 550 test: 551 suffix: linear 552 args: -sol_form linear -mesh_transform none 553 test: 554 suffix: sinusoidal 555 args: -sol_form sinusoidal -mesh_transform none 556 test: 557 suffix: sinusoidal_skew 558 args: -sol_form sinusoidal -mesh_transform skew 559 test: 560 suffix: sinusoidal_perturb 561 args: -sol_form sinusoidal -mesh_transform perturb 562 test: 563 suffix: sinusoidal_skew_perturb 564 args: -sol_form sinusoidal -mesh_transform skew_perturb 565 566 testset: 567 TODO: broken 568 suffix: 2d_bdmq 569 args: -dim 2 \ 570 -simplex false \ 571 -velocity_petscspace_degree 1 \ 572 -velocity_petscdualspace_type bdm \ 573 -velocity_petscdualspace_lagrange_tensor 1 \ 574 -divErr_petscspace_degree 1 \ 575 -divErr_petscdualspace_lagrange_continuity false \ 576 -dm_refine 0 \ 577 -snes_error_if_not_converged \ 578 -ksp_rtol 1e-10 \ 579 -ksp_error_if_not_converged \ 580 -pc_type fieldsplit\ 581 -pc_fieldsplit_detect_saddle_point\ 582 -pc_fieldsplit_type schur\ 583 -pc_fieldsplit_schur_precondition full 584 test: 585 suffix: linear 586 args: -sol_form linear -mesh_transform none 587 test: 588 suffix: sinusoidal 589 args: -sol_form sinusoidal -mesh_transform none 590 test: 591 suffix: sinusoidal_skew 592 args: -sol_form sinusoidal -mesh_transform skew 593 test: 594 suffix: sinusoidal_perturb 595 args: -sol_form sinusoidal -mesh_transform perturb 596 test: 597 suffix: sinusoidal_skew_perturb 598 args: -sol_form sinusoidal -mesh_transform skew_perturb 599 600 testset: 601 suffix: 3d_bdm 602 requires: ctetgen 603 args: -dim 3 \ 604 -simplex true \ 605 -velocity_petscspace_degree 1 \ 606 -velocity_petscdualspace_type bdm \ 607 -divErr_petscspace_degree 1 \ 608 -divErr_petscdualspace_lagrange_continuity false \ 609 -dm_refine 0 \ 610 -snes_error_if_not_converged \ 611 -ksp_rtol 1e-10 \ 612 -ksp_error_if_not_converged \ 613 -pc_type fieldsplit \ 614 -pc_fieldsplit_detect_saddle_point \ 615 -pc_fieldsplit_type schur \ 616 -pc_fieldsplit_schur_precondition full 617 test: 618 suffix: linear 619 args: -sol_form linear -mesh_transform none 620 test: 621 suffix: sinusoidal 622 args: -sol_form sinusoidal -mesh_transform none 623 test: 624 suffix: sinusoidal_skew 625 args: -sol_form sinusoidal -mesh_transform skew 626 test: 627 suffix: sinusoidal_perturb 628 args: -sol_form sinusoidal -mesh_transform perturb 629 test: 630 suffix: sinusoidal_skew_perturb 631 args: -sol_form sinusoidal -mesh_transform skew_perturb 632 633 testset: 634 TODO: broken 635 suffix: 3d_bdmq 636 requires: ctetgen 637 args: -dim 3 \ 638 -simplex false \ 639 -velocity_petscspace_degree 1 \ 640 -velocity_petscdualspace_type bdm \ 641 -velocity_petscdualspace_lagrange_tensor 1 \ 642 -divErr_petscspace_degree 1 \ 643 -divErr_petscdualspace_lagrange_continuity false \ 644 -dm_refine 0 \ 645 -snes_error_if_not_converged \ 646 -ksp_rtol 1e-10 \ 647 -ksp_error_if_not_converged \ 648 -pc_type fieldsplit \ 649 -pc_fieldsplit_detect_saddle_point \ 650 -pc_fieldsplit_type schur \ 651 -pc_fieldsplit_schur_precondition full 652 test: 653 suffix: linear 654 args: -sol_form linear -mesh_transform none 655 test: 656 suffix: sinusoidal 657 args: -sol_form sinusoidal -mesh_transform none 658 test: 659 suffix: sinusoidal_skew 660 args: -sol_form sinusoidal -mesh_transform skew 661 test: 662 suffix: sinusoidal_perturb 663 args: -sol_form sinusoidal -mesh_transform perturb 664 test: 665 suffix: sinusoidal_skew_perturb 666 args: -sol_form sinusoidal -mesh_transform skew_perturb 667 668 test: 669 suffix: quad_rt_0 670 args: -dim 2 -simplex false -mesh_transform skew \ 671 -divErr_petscspace_degree 1 \ 672 -divErr_petscdualspace_lagrange_continuity false \ 673 -dm_refine 0 \ 674 -snes_error_if_not_converged \ 675 -ksp_rtol 1e-10 \ 676 -ksp_error_if_not_converged \ 677 -pc_type fieldsplit\ 678 -pc_fieldsplit_detect_saddle_point\ 679 -pc_fieldsplit_type schur\ 680 -pc_fieldsplit_schur_precondition full \ 681 -velocity_petscfe_default_quadrature_order 1 \ 682 -velocity_petscspace_type sum \ 683 -velocity_petscspace_variables 2 \ 684 -velocity_petscspace_components 2 \ 685 -velocity_petscspace_sum_spaces 2 \ 686 -velocity_petscspace_sum_concatenate true \ 687 -velocity_subspace0_petscspace_variables 2 \ 688 -velocity_subspace0_petscspace_type tensor \ 689 -velocity_subspace0_petscspace_tensor_spaces 2 \ 690 -velocity_subspace0_petscspace_tensor_uniform false \ 691 -velocity_subspace0_subspace_0_petscspace_degree 1 \ 692 -velocity_subspace0_subspace_1_petscspace_degree 0 \ 693 -velocity_subspace1_petscspace_variables 2 \ 694 -velocity_subspace1_petscspace_type tensor \ 695 -velocity_subspace1_petscspace_tensor_spaces 2 \ 696 -velocity_subspace1_petscspace_tensor_uniform false \ 697 -velocity_subspace1_subspace_0_petscspace_degree 0 \ 698 -velocity_subspace1_subspace_1_petscspace_degree 1 \ 699 -velocity_petscdualspace_form_degree -1 \ 700 -velocity_petscdualspace_order 1 \ 701 -velocity_petscdualspace_lagrange_trimmed true 702 TEST*/ 703