1*c4762a1bSJed Brown static char help[] = "Variable-Viscosity Stokes Problem in 2d.\n\ 2*c4762a1bSJed Brown Exact solutions provided by Mirko Velic.\n\n\n"; 3*c4762a1bSJed Brown 4*c4762a1bSJed Brown #include<petsc.h> 5*c4762a1bSJed Brown 6*c4762a1bSJed Brown #include "ex75.h" 7*c4762a1bSJed Brown 8*c4762a1bSJed Brown typedef struct { 9*c4762a1bSJed Brown PetscBool fem; /* Flag for FEM tests */ 10*c4762a1bSJed Brown } AppCtx; 11*c4762a1bSJed Brown 12*c4762a1bSJed Brown PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 13*c4762a1bSJed Brown { 14*c4762a1bSJed Brown PetscErrorCode ierr; 15*c4762a1bSJed Brown 16*c4762a1bSJed Brown PetscFunctionBeginUser; 17*c4762a1bSJed Brown options->fem = PETSC_FALSE; 18*c4762a1bSJed Brown 19*c4762a1bSJed Brown ierr = PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");CHKERRQ(ierr); 20*c4762a1bSJed Brown ierr = PetscOptionsBool("-fem", "Run FEM tests", "ex75.c", options->fem, &options->fem, NULL);CHKERRQ(ierr); 21*c4762a1bSJed Brown ierr = PetscOptionsEnd(); 22*c4762a1bSJed Brown PetscFunctionReturn(0); 23*c4762a1bSJed Brown } 24*c4762a1bSJed Brown 25*c4762a1bSJed Brown /* 26*c4762a1bSJed Brown SolKxSolution - Exact Stokes solutions for exponentially varying viscosity 27*c4762a1bSJed Brown 28*c4762a1bSJed Brown Input Parameters: 29*c4762a1bSJed Brown + x - The x coordinate at which to evaluate the solution 30*c4762a1bSJed Brown . z - The z coordinate at which to evaluate the solution 31*c4762a1bSJed Brown . kn - The constant defining the x-dependence of the forcing function 32*c4762a1bSJed Brown . km - The constant defining the z-dependence of the forcing function 33*c4762a1bSJed Brown - B - The viscosity coefficient 34*c4762a1bSJed Brown 35*c4762a1bSJed Brown Output Parameters: 36*c4762a1bSJed Brown + vx - The x-velocity at (x,z) 37*c4762a1bSJed Brown . vz - The z-velocity at (x,z) 38*c4762a1bSJed Brown . p - The pressure at (x,z) 39*c4762a1bSJed Brown . sxx - The stress sigma_xx at (x,z) 40*c4762a1bSJed Brown . sxz - The stress sigma_xz at (x,z) 41*c4762a1bSJed Brown - szz - The stress sigma_zz at (x,z) 42*c4762a1bSJed Brown 43*c4762a1bSJed Brown Note: 44*c4762a1bSJed Brown $ The domain is the square 0 <= x,z <= 1. We solve the Stokes equation for incompressible flow with free-slip boundary 45*c4762a1bSJed Brown $ conditions everywhere. The forcing term f is given by 46*c4762a1bSJed Brown $ 47*c4762a1bSJed Brown $ fx = 0 48*c4762a1bSJed Brown $ fz = sigma*sin(km*z)*cos(kn*x) 49*c4762a1bSJed Brown $ 50*c4762a1bSJed Brown $ where 51*c4762a1bSJed Brown $ 52*c4762a1bSJed Brown $ km = m*Pi (m may be non-integral) 53*c4762a1bSJed Brown $ kn = n*Pi 54*c4762a1bSJed Brown $ 55*c4762a1bSJed Brown $ meaning that the density rho is -sigma*sin(km*z)*cos(kn*x). The viscosity eta is exp(2*B*x). 56*c4762a1bSJed Brown */ 57*c4762a1bSJed Brown PetscErrorCode SolKxSolution(PetscReal x, PetscReal z, PetscReal kn, PetscReal km, PetscReal B, PetscScalar *vx, PetscScalar *vz, PetscScalar *p, PetscScalar *sxx, PetscScalar *sxz, PetscScalar *szz) 58*c4762a1bSJed Brown { 59*c4762a1bSJed Brown PetscScalar sigma; 60*c4762a1bSJed Brown PetscScalar _C1,_C2,_C3,_C4; 61*c4762a1bSJed Brown PetscScalar Rp, UU, VV; 62*c4762a1bSJed Brown PetscScalar a,b,r,_aa,_bb,AA,BB,Rm; 63*c4762a1bSJed Brown PetscScalar num1,num2,num3,num4,den1; 64*c4762a1bSJed Brown 65*c4762a1bSJed Brown PetscScalar t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; 66*c4762a1bSJed Brown PetscScalar t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21; 67*c4762a1bSJed Brown PetscScalar t22,t23,t24,t25,t26,t28,t29,t30,t31,t32; 68*c4762a1bSJed Brown PetscScalar t33,t34,t35,t36,t37,t38,t39,t40,t41,t42; 69*c4762a1bSJed Brown PetscScalar t44,t45,t46,t47,t48,t49,t51,t52,t53,t54; 70*c4762a1bSJed Brown PetscScalar t56,t58,t61,t62,t63,t64,t65,t66,t67,t68; 71*c4762a1bSJed Brown PetscScalar t69,t70,t71,t72,t73,t74,t75,t76,t77,t78; 72*c4762a1bSJed Brown PetscScalar t79,t80,t81,t82,t83,t84,t85,t86,t87,t88; 73*c4762a1bSJed Brown PetscScalar t89,t90,t91,t92,t93,t94,t95,t96,t97,t98; 74*c4762a1bSJed Brown PetscScalar t99,t100,t101,t103,t104,t105,t106,t107,t108,t109; 75*c4762a1bSJed Brown PetscScalar t110,t111,t112,t113,t114,t115,t116,t117,t118,t119; 76*c4762a1bSJed Brown PetscScalar t120,t121,t123,t125,t127,t128,t130,t131,t132,t133; 77*c4762a1bSJed Brown PetscScalar t135,t136,t138,t140,t141,t142,t143,t152,t160,t162; 78*c4762a1bSJed Brown 79*c4762a1bSJed Brown PetscFunctionBegin; 80*c4762a1bSJed Brown /*************************************************************************/ 81*c4762a1bSJed Brown /*************************************************************************/ 82*c4762a1bSJed Brown /* rho = -sin(km*z)*cos(kn*x) */ 83*c4762a1bSJed Brown /* viscosity Z= exp(2*B*z) */ 84*c4762a1bSJed Brown /* solution valid for km not zero -- should get trivial solution if km=0 */ 85*c4762a1bSJed Brown sigma = 1.0; 86*c4762a1bSJed Brown /*************************************************************************/ 87*c4762a1bSJed Brown /*************************************************************************/ 88*c4762a1bSJed Brown a = B*B + km*km; 89*c4762a1bSJed Brown b = 2.0*km*B; 90*c4762a1bSJed Brown r = sqrt(a*a + b*b); 91*c4762a1bSJed Brown Rp = sqrt( (r+a)/2.0 ); 92*c4762a1bSJed Brown Rm = sqrt( (r-a)/2.0 ); 93*c4762a1bSJed Brown UU = Rp - B; 94*c4762a1bSJed Brown VV = Rp + B; 95*c4762a1bSJed Brown 96*c4762a1bSJed Brown /*******************************************/ 97*c4762a1bSJed Brown /* calculate the constants */ 98*c4762a1bSJed Brown /*******************************************/ 99*c4762a1bSJed Brown t1 = kn * kn; 100*c4762a1bSJed Brown t4 = km * km; 101*c4762a1bSJed Brown t6 = t4 * t4; 102*c4762a1bSJed Brown t7 = B * B; 103*c4762a1bSJed Brown t9 = 0.4e1 * t7 * t4; 104*c4762a1bSJed Brown t12 = 0.8e1 * t7 * kn * km; 105*c4762a1bSJed Brown t14 = 0.4e1 * t7 * t1; 106*c4762a1bSJed Brown t16 = 0.2e1 * t4 * t1; 107*c4762a1bSJed Brown t17 = t1 * t1; 108*c4762a1bSJed Brown _aa = -0.4e1 * B * t1 * sigma * (t4 + t1) / (t6 + t9 + t12 + t14 + t16 + t17) / (t6 + t9 - t12 + t14 + t16 + t17); 109*c4762a1bSJed Brown 110*c4762a1bSJed Brown t2 = kn * kn; 111*c4762a1bSJed Brown t3 = t2 * t2; 112*c4762a1bSJed Brown t4 = B * B; 113*c4762a1bSJed Brown t6 = 0.4e1 * t4 * t2; 114*c4762a1bSJed Brown t7 = km * km; 115*c4762a1bSJed Brown t9 = 0.4e1 * t7 * t4; 116*c4762a1bSJed Brown t10 = t7 * t7; 117*c4762a1bSJed Brown t12 = 0.2e1 * t7 * t2; 118*c4762a1bSJed Brown t16 = 0.8e1 * t4 * kn * km; 119*c4762a1bSJed Brown _bb = sigma * kn * (t3 - t6 + t9 + t10 + t12) / (t10 + t9 + t16 + t6 + t12 + t3) / (t10 + t9 - t16 + t6 + t12 + t3); 120*c4762a1bSJed Brown 121*c4762a1bSJed Brown AA = _aa; 122*c4762a1bSJed Brown BB = _bb; 123*c4762a1bSJed Brown 124*c4762a1bSJed Brown t1 = Rm * Rm; 125*c4762a1bSJed Brown t2 = B - Rp; 126*c4762a1bSJed Brown t4 = Rp + B; 127*c4762a1bSJed Brown t6 = UU * x; 128*c4762a1bSJed Brown t9 = exp(t6 - 0.4e1 * Rp); 129*c4762a1bSJed Brown t13 = kn * kn; 130*c4762a1bSJed Brown t15 = B * B; 131*c4762a1bSJed Brown t18 = Rp * Rp; 132*c4762a1bSJed Brown t19 = t18 * B; 133*c4762a1bSJed Brown t20 = t15 * Rp; 134*c4762a1bSJed Brown t22 = t1 * Rp; 135*c4762a1bSJed Brown t24 = B * t1; 136*c4762a1bSJed Brown t32 = 0.8e1 * t15 * BB * kn * Rp; 137*c4762a1bSJed Brown t34 = 0.2e1 * Rm; 138*c4762a1bSJed Brown t35 = cos(t34); 139*c4762a1bSJed Brown t37 = Rm * Rp; 140*c4762a1bSJed Brown t49 = sin(t34); 141*c4762a1bSJed Brown t63 = exp(t6 - 0.2e1 * Rp); 142*c4762a1bSJed Brown t65 = Rm * t2; 143*c4762a1bSJed Brown t67 = 0.2e1 * B * kn; 144*c4762a1bSJed Brown t68 = B * Rm; 145*c4762a1bSJed Brown t69 = t67 + t68 + t37; 146*c4762a1bSJed Brown t73 = 0.3e1 * t15; 147*c4762a1bSJed Brown t75 = 0.2e1 * B * Rp; 148*c4762a1bSJed Brown t76 = t73 - t75 + t1 - t13 - t18; 149*c4762a1bSJed Brown t78 = t65 * t76 * BB; 150*c4762a1bSJed Brown t80 = Rm - kn; 151*c4762a1bSJed Brown t81 = cos(t80); 152*c4762a1bSJed Brown t83 = t68 - t67 + t37; 153*c4762a1bSJed Brown t88 = Rm + kn; 154*c4762a1bSJed Brown t89 = cos(t88); 155*c4762a1bSJed Brown t92 = t65 * t76 * AA; 156*c4762a1bSJed Brown t97 = sin(t80); 157*c4762a1bSJed Brown t103 = sin(t88); 158*c4762a1bSJed Brown t108 = exp(t6 - 0.3e1 * Rp - B); 159*c4762a1bSJed Brown t110 = Rm * t4; 160*c4762a1bSJed Brown t111 = t67 + t68 - t37; 161*c4762a1bSJed Brown t115 = t73 + t75 + t1 - t13 - t18; 162*c4762a1bSJed Brown t117 = t110 * t115 * BB; 163*c4762a1bSJed Brown t120 = -t67 + t68 - t37; 164*c4762a1bSJed Brown t127 = t110 * t115 * AA; 165*c4762a1bSJed Brown t140 = exp(t6 - Rp - B); 166*c4762a1bSJed Brown num1 = -0.4e1 * t1 * t2 * t4 * AA * t9 + ((0.2e1 * Rp * (-B * t13 + 0.3e1 * t15 * B - t19 - 0.2e1 * t20 - 0.2e1 * t22 - t24) * AA - t32) * t35 + (0.2e1 * t37 * (t1 - t13 + 0.5e1 * t15 - t18) * AA - 0.8e1 * B * BB * kn * Rm * Rp) * t49 - 0.2e1 * B * (0.3e1 * t20 - t18 * Rp - 0.2e1 * t19 - Rp * t13 - t22 - 0.2e1 * t24) * AA + t32) * t63 + ((0.2e1 * t65 * t69 * AA + t78) * t81 + (0.2e1 * t65 * t83 * AA - t78) * t89 + (t92 - 0.2e1 * t65 * t69 * BB) * t97 + (t92 + 0.2e1 * t65 * t83 * BB) * t103) * t108 + ((-0.2e1 * t110 * t111 * AA - t117) * t81 + (-0.2e1 * t110 * t120 * AA + t117) * t89 + (-t127 + 0.2e1 * t110 * t111 * BB) * t97 + (-t127 - 0.2e1 * t110 * t120 * BB) * t103) * t140; 167*c4762a1bSJed Brown 168*c4762a1bSJed Brown t1 = Rp + B; 169*c4762a1bSJed Brown t2 = Rm * t1; 170*c4762a1bSJed Brown t3 = B * B; 171*c4762a1bSJed Brown t4 = 0.3e1 * t3; 172*c4762a1bSJed Brown t5 = B * Rp; 173*c4762a1bSJed Brown t7 = Rm * Rm; 174*c4762a1bSJed Brown t8 = kn * kn; 175*c4762a1bSJed Brown t9 = Rp * Rp; 176*c4762a1bSJed Brown t10 = t4 + 0.2e1 * t5 + t7 - t8 - t9; 177*c4762a1bSJed Brown t12 = t2 * t10 * AA; 178*c4762a1bSJed Brown t14 = B * Rm; 179*c4762a1bSJed Brown t20 = UU * x; 180*c4762a1bSJed Brown t23 = exp(t20 - 0.4e1 * Rp); 181*c4762a1bSJed Brown t25 = Rm * Rp; 182*c4762a1bSJed Brown t32 = Rm * kn; 183*c4762a1bSJed Brown t37 = 0.2e1 * Rm; 184*c4762a1bSJed Brown t38 = cos(t37); 185*c4762a1bSJed Brown t41 = t3 * B; 186*c4762a1bSJed Brown t44 = t3 * Rp; 187*c4762a1bSJed Brown t48 = B * t7; 188*c4762a1bSJed Brown t53 = t3 * BB; 189*c4762a1bSJed Brown t54 = kn * Rp; 190*c4762a1bSJed Brown t58 = sin(t37); 191*c4762a1bSJed Brown t69 = exp(t20 - 0.2e1 * Rp); 192*c4762a1bSJed Brown t71 = t9 * Rp; 193*c4762a1bSJed Brown t72 = Rm * t71; 194*c4762a1bSJed Brown t73 = t3 * Rm; 195*c4762a1bSJed Brown t75 = 0.5e1 * t73 * Rp; 196*c4762a1bSJed Brown t77 = 0.8e1 * t44 * kn; 197*c4762a1bSJed Brown t78 = t25 * t8; 198*c4762a1bSJed Brown t79 = t7 * Rm; 199*c4762a1bSJed Brown t80 = B * t79; 200*c4762a1bSJed Brown t81 = t14 * t8; 201*c4762a1bSJed Brown t82 = t79 * Rp; 202*c4762a1bSJed Brown t84 = 0.3e1 * t41 * Rm; 203*c4762a1bSJed Brown t85 = t14 * t9; 204*c4762a1bSJed Brown t86 = -t72 + t75 + t77 - t78 + t80 - t81 + t82 + t84 + t85; 205*c4762a1bSJed Brown t88 = t7 * t9; 206*c4762a1bSJed Brown t89 = t5 * t8; 207*c4762a1bSJed Brown t90 = t7 * t3; 208*c4762a1bSJed Brown t91 = B * t71; 209*c4762a1bSJed Brown t92 = t48 * Rp; 210*c4762a1bSJed Brown t94 = 0.2e1 * t14 * t54; 211*c4762a1bSJed Brown t96 = 0.3e1 * Rp * t41; 212*c4762a1bSJed Brown t98 = 0.2e1 * t73 * kn; 213*c4762a1bSJed Brown t100 = 0.2e1 * t9 * t3; 214*c4762a1bSJed Brown t101 = -t88 - t89 - t90 - t91 - t92 - t94 + t96 - t98 - t100; 215*c4762a1bSJed Brown t105 = Rm - kn; 216*c4762a1bSJed Brown t106 = cos(t105); 217*c4762a1bSJed Brown t108 = t75 - t77 - t78 + t85 - t72 - t81 + t80 + t84 + t82; 218*c4762a1bSJed Brown t110 = -t100 + t96 - t91 + t94 + t98 - t92 - t89 - t88 - t90; 219*c4762a1bSJed Brown t114 = Rm + kn; 220*c4762a1bSJed Brown t115 = cos(t114); 221*c4762a1bSJed Brown t121 = sin(t105); 222*c4762a1bSJed Brown t127 = sin(t114); 223*c4762a1bSJed Brown t132 = exp(t20 - 0.3e1 * Rp - B); 224*c4762a1bSJed Brown t135 = 0.2e1 * B * kn; 225*c4762a1bSJed Brown t136 = t135 + t14 - t25; 226*c4762a1bSJed Brown t142 = -t135 + t14 - t25; 227*c4762a1bSJed Brown t152 = t2 * t10 * BB; 228*c4762a1bSJed Brown t162 = exp(t20 - Rp - B); 229*c4762a1bSJed Brown num2 = (0.2e1 * t12 - 0.8e1 * t14 * kn * t1 * BB) * t23 + ((-0.2e1 * t25 * (t7 - t8 + 0.5e1 * t3 - t9) * AA + 0.8e1 * B * BB * t32 * Rp) * t38 + (0.2e1 * Rp * (-B * t8 + 0.3e1 * t41 - t9 * B - 0.2e1 * t44 - 0.2e1 * t7 * Rp - t48) * AA - 0.8e1 * t53 * t54) * t58 - 0.2e1 * t14 * (t4 + t9 - t8 + t7) * AA + 0.8e1 * t53 * t32) * t69 + ((-t86 * AA - 0.2e1 * t101 * BB) * t106 + (-t108 * AA + 0.2e1 * t110 * BB) * t115 + (-0.2e1 * t101 * AA + t86 * BB) * t121 + (-0.2e1 * t110 * AA - t108 * BB) * t127) * t132 + ((t12 - 0.2e1 * t2 * t136 * BB) * t106 + (t12 + 0.2e1 * t2 * t142 * BB) * t115 + (-0.2e1 * t2 * t136 * AA - t152) * t121 + (-0.2e1 * t2 * t142 * AA + t152) * t127) * t162; 230*c4762a1bSJed Brown 231*c4762a1bSJed Brown t1 = Rm * Rm; 232*c4762a1bSJed Brown t2 = B - Rp; 233*c4762a1bSJed Brown t4 = Rp + B; 234*c4762a1bSJed Brown t6 = VV * x; 235*c4762a1bSJed Brown t7 = exp(-t6); 236*c4762a1bSJed Brown t11 = B * t1; 237*c4762a1bSJed Brown t12 = Rp * Rp; 238*c4762a1bSJed Brown t13 = t12 * B; 239*c4762a1bSJed Brown t14 = B * B; 240*c4762a1bSJed Brown t15 = t14 * Rp; 241*c4762a1bSJed Brown t19 = kn * kn; 242*c4762a1bSJed Brown t21 = t1 * Rp; 243*c4762a1bSJed Brown t30 = 0.8e1 * t14 * BB * kn * Rp; 244*c4762a1bSJed Brown t32 = 0.2e1 * Rm; 245*c4762a1bSJed Brown t33 = cos(t32); 246*c4762a1bSJed Brown t35 = Rm * Rp; 247*c4762a1bSJed Brown t47 = sin(t32); 248*c4762a1bSJed Brown t61 = exp(-t6 - 0.2e1 * Rp); 249*c4762a1bSJed Brown t63 = Rm * t2; 250*c4762a1bSJed Brown t65 = 0.2e1 * B * kn; 251*c4762a1bSJed Brown t66 = B * Rm; 252*c4762a1bSJed Brown t67 = t65 + t66 + t35; 253*c4762a1bSJed Brown t71 = 0.3e1 * t14; 254*c4762a1bSJed Brown t73 = 0.2e1 * B * Rp; 255*c4762a1bSJed Brown t74 = t71 - t73 + t1 - t19 - t12; 256*c4762a1bSJed Brown t76 = t63 * t74 * BB; 257*c4762a1bSJed Brown t78 = Rm - kn; 258*c4762a1bSJed Brown t79 = cos(t78); 259*c4762a1bSJed Brown t81 = t66 - t65 + t35; 260*c4762a1bSJed Brown t86 = Rm + kn; 261*c4762a1bSJed Brown t87 = cos(t86); 262*c4762a1bSJed Brown t90 = t63 * t74 * AA; 263*c4762a1bSJed Brown t95 = sin(t78); 264*c4762a1bSJed Brown t101 = sin(t86); 265*c4762a1bSJed Brown t106 = exp(-t6 - 0.3e1 * Rp - B); 266*c4762a1bSJed Brown t108 = Rm * t4; 267*c4762a1bSJed Brown t109 = t65 + t66 - t35; 268*c4762a1bSJed Brown t113 = t71 + t73 + t1 - t19 - t12; 269*c4762a1bSJed Brown t115 = t108 * t113 * BB; 270*c4762a1bSJed Brown t118 = -t65 + t66 - t35; 271*c4762a1bSJed Brown t125 = t108 * t113 * AA; 272*c4762a1bSJed Brown t138 = exp(-t6 - Rp - B); 273*c4762a1bSJed Brown num3 = -0.4e1 * t1 * t2 * t4 * AA * t7 + ((-0.2e1 * Rp * (-t11 - t13 + 0.2e1 * t15 + 0.3e1 * t14 * B - B * t19 + 0.2e1 * t21) * AA + t30) * t33 + (-0.2e1 * t35 * (t1 - t19 + 0.5e1 * t14 - t12) * AA + 0.8e1 * B * BB * kn * Rm * Rp) * t47 + 0.2e1 * B * (-t12 * Rp + 0.2e1 * t11 + 0.3e1 * t15 + 0.2e1 * t13 - t21 - Rp * t19) * AA - t30) * t61 + ((-0.2e1 * t63 * t67 * AA - t76) * t79 + (-0.2e1 * t63 * t81 * AA + t76) * t87 + (-t90 + 0.2e1 * t63 * t67 * BB) * t95 + (-t90 - 0.2e1 * t63 * t81 * BB) * t101) * t106 + ((0.2e1 * t108 * t109 * AA + t115) * t79 + (0.2e1 * t108 * t118 * AA - t115) * t87 + (t125 - 0.2e1 * t108 * t109 * BB) * t95 + (t125 + 0.2e1 * t108 * t118 * BB) * t101) * t138; 274*c4762a1bSJed Brown 275*c4762a1bSJed Brown t1 = B - Rp; 276*c4762a1bSJed Brown t2 = Rm * t1; 277*c4762a1bSJed Brown t3 = B * B; 278*c4762a1bSJed Brown t4 = 0.3e1 * t3; 279*c4762a1bSJed Brown t5 = B * Rp; 280*c4762a1bSJed Brown t7 = Rm * Rm; 281*c4762a1bSJed Brown t8 = kn * kn; 282*c4762a1bSJed Brown t9 = Rp * Rp; 283*c4762a1bSJed Brown t10 = t4 - 0.2e1 * t5 + t7 - t8 - t9; 284*c4762a1bSJed Brown t12 = t2 * t10 * AA; 285*c4762a1bSJed Brown t14 = B * Rm; 286*c4762a1bSJed Brown t20 = VV * x; 287*c4762a1bSJed Brown t21 = exp(-t20); 288*c4762a1bSJed Brown t23 = Rm * Rp; 289*c4762a1bSJed Brown t30 = Rm * kn; 290*c4762a1bSJed Brown t35 = 0.2e1 * Rm; 291*c4762a1bSJed Brown t36 = cos(t35); 292*c4762a1bSJed Brown t38 = B * t7; 293*c4762a1bSJed Brown t40 = t3 * Rp; 294*c4762a1bSJed Brown t42 = t3 * B; 295*c4762a1bSJed Brown t51 = t3 * BB; 296*c4762a1bSJed Brown t52 = kn * Rp; 297*c4762a1bSJed Brown t56 = sin(t35); 298*c4762a1bSJed Brown t67 = exp(-t20 - 0.2e1 * Rp); 299*c4762a1bSJed Brown t70 = 0.2e1 * B * kn; 300*c4762a1bSJed Brown t71 = t70 + t14 + t23; 301*c4762a1bSJed Brown t76 = Rm - kn; 302*c4762a1bSJed Brown t77 = cos(t76); 303*c4762a1bSJed Brown t79 = t14 - t70 + t23; 304*c4762a1bSJed Brown t84 = Rm + kn; 305*c4762a1bSJed Brown t85 = cos(t84); 306*c4762a1bSJed Brown t91 = t2 * t10 * BB; 307*c4762a1bSJed Brown t93 = sin(t76); 308*c4762a1bSJed Brown t99 = sin(t84); 309*c4762a1bSJed Brown t104 = exp(-t20 - 0.3e1 * Rp - B); 310*c4762a1bSJed Brown t106 = t9 * Rp; 311*c4762a1bSJed Brown t107 = Rm * t106; 312*c4762a1bSJed Brown t108 = t3 * Rm; 313*c4762a1bSJed Brown t110 = 0.5e1 * t108 * Rp; 314*c4762a1bSJed Brown t112 = 0.8e1 * t40 * kn; 315*c4762a1bSJed Brown t113 = t23 * t8; 316*c4762a1bSJed Brown t114 = t7 * Rm; 317*c4762a1bSJed Brown t115 = B * t114; 318*c4762a1bSJed Brown t116 = t14 * t8; 319*c4762a1bSJed Brown t117 = t114 * Rp; 320*c4762a1bSJed Brown t119 = 0.3e1 * t42 * Rm; 321*c4762a1bSJed Brown t120 = t14 * t9; 322*c4762a1bSJed Brown t121 = t107 - t110 - t112 + t113 + t115 - t116 - t117 + t119 + t120; 323*c4762a1bSJed Brown t123 = t38 * Rp; 324*c4762a1bSJed Brown t125 = 0.2e1 * t14 * t52; 325*c4762a1bSJed Brown t127 = 0.3e1 * Rp * t42; 326*c4762a1bSJed Brown t128 = t7 * t3; 327*c4762a1bSJed Brown t130 = 0.2e1 * t9 * t3; 328*c4762a1bSJed Brown t131 = t7 * t9; 329*c4762a1bSJed Brown t132 = B * t106; 330*c4762a1bSJed Brown t133 = t5 * t8; 331*c4762a1bSJed Brown t135 = 0.2e1 * t108 * kn; 332*c4762a1bSJed Brown t136 = -t123 - t125 + t127 + t128 + t130 + t131 - t132 - t133 + t135; 333*c4762a1bSJed Brown t141 = -t110 + t112 + t113 + t120 + t107 - t116 + t115 + t119 - t117; 334*c4762a1bSJed Brown t143 = t125 - t132 + t130 - t135 + t127 + t131 - t123 + t128 - t133; 335*c4762a1bSJed Brown t160 = exp(-t20 - Rp - B); 336*c4762a1bSJed Brown num4 = (0.2e1 * t12 - 0.8e1 * t14 * kn * t1 * BB) * t21 + ((0.2e1 * t23 * (t7 - t8 + 0.5e1 * t3 - t9) * AA - 0.8e1 * B * BB * t30 * Rp) * t36 + (-0.2e1 * Rp * (-t38 - t9 * B + 0.2e1 * t40 + 0.3e1 * t42 - B * t8 + 0.2e1 * t7 * Rp) * AA + 0.8e1 * t51 * t52) * t56 - 0.2e1 * t14 * (t4 + t9 - t8 + t7) * AA + 0.8e1 * t51 * t30) * t67 + ((t12 - 0.2e1 * t2 * t71 * BB) * t77 + (t12 + 0.2e1 * t2 * t79 * BB) * t85 + (-0.2e1 * t2 * t71 * AA - t91) * t93 + (-0.2e1 * t2 * t79 * AA + t91) * t99) * t104 + ((-t121 * AA + 0.2e1 * t136 * BB) * t77 + (-t141 * AA - 0.2e1 * t143 * BB) * t85 + (0.2e1 * t136 * AA + t121 * BB) * t93 + (0.2e1 * t143 * AA - t141 * BB) * t99) * t160; 337*c4762a1bSJed Brown 338*c4762a1bSJed Brown 339*c4762a1bSJed Brown t1 = Rm * Rm; 340*c4762a1bSJed Brown t2 = Rp * Rp; 341*c4762a1bSJed Brown t3 = t1 * t2; 342*c4762a1bSJed Brown t4 = B * B; 343*c4762a1bSJed Brown t5 = t1 * t4; 344*c4762a1bSJed Brown t9 = exp(-0.4e1 * Rp); 345*c4762a1bSJed Brown t15 = cos(0.2e1 * Rm); 346*c4762a1bSJed Brown t22 = exp(-0.2e1 * Rp); 347*c4762a1bSJed Brown den1 = (-0.4e1 * t3 + 0.4e1 * t5) * t9 + ((0.8e1 * t1 + 0.8e1 * t4) * t2 * t15 - 0.8e1 * t5 - 0.8e1 * t2 * t4) * t22 - 0.4e1 * t3 + 0.4e1 * t5; 348*c4762a1bSJed Brown 349*c4762a1bSJed Brown _C1=num1/den1; _C2=num2/den1; _C3=num3/den1; _C4=num4/den1; 350*c4762a1bSJed Brown 351*c4762a1bSJed Brown /*******************************************/ 352*c4762a1bSJed Brown /* calculate solution */ 353*c4762a1bSJed Brown /*******************************************/ 354*c4762a1bSJed Brown t1 = Rm * x; 355*c4762a1bSJed Brown t2 = cos(t1); 356*c4762a1bSJed Brown t4 = sin(t1); 357*c4762a1bSJed Brown t10 = exp(-0.2e1 * x * B); 358*c4762a1bSJed Brown t12 = kn * x; 359*c4762a1bSJed Brown t13 = cos(t12); 360*c4762a1bSJed Brown t16 = sin(t12); 361*c4762a1bSJed Brown *vx = -km * (_C1 * t2 + _C2 * t4 + _C3 * t2 + _C4 * t4 + t10 * AA * t13 + t10 * BB * t16); 362*c4762a1bSJed Brown 363*c4762a1bSJed Brown t2 = Rm * x; 364*c4762a1bSJed Brown t3 = cos(t2); 365*c4762a1bSJed Brown t6 = sin(t2); 366*c4762a1bSJed Brown t22 = exp(-0.2e1 * x * B); 367*c4762a1bSJed Brown t23 = B * t22; 368*c4762a1bSJed Brown t24 = kn * x; 369*c4762a1bSJed Brown t25 = cos(t24); 370*c4762a1bSJed Brown t29 = sin(t24); 371*c4762a1bSJed Brown *vz = UU * _C1 * t3 + UU * _C2 * t6 - _C1 * t6 * Rm + _C2 * t3 * Rm - VV * _C3 * t3 - VV * _C4 * t6 - _C3 * t6 * Rm + _C4 * t3 * Rm - 0.2e1 * t23 * AA * t25 - 0.2e1 * t23 * BB * t29 - t22 * AA * t29 * kn + t22 * BB * t25 * kn; 372*c4762a1bSJed Brown 373*c4762a1bSJed Brown t3 = exp(0.2e1 * x * B); 374*c4762a1bSJed Brown t4 = t3 * B; 375*c4762a1bSJed Brown t8 = km * km; 376*c4762a1bSJed Brown t9 = t3 * t8; 377*c4762a1bSJed Brown t11 = 0.3e1 * t9 * Rm; 378*c4762a1bSJed Brown t12 = Rm * Rm; 379*c4762a1bSJed Brown t14 = t3 * t12 * Rm; 380*c4762a1bSJed Brown t15 = UU * UU; 381*c4762a1bSJed Brown t19 = 0.4e1 * t4 * UU * Rm - t11 - t14 + 0.3e1 * t3 * t15 * Rm; 382*c4762a1bSJed Brown t20 = Rm * x; 383*c4762a1bSJed Brown t21 = sin(t20); 384*c4762a1bSJed Brown t26 = 0.2e1 * t9 * B; 385*c4762a1bSJed Brown t33 = 0.2e1 * t4 * t12; 386*c4762a1bSJed Brown t36 = -t3 * t15 * UU - t26 + 0.3e1 * t9 * UU + 0.3e1 * t3 * UU * t12 + t33 - 0.2e1 * t4 * t15; 387*c4762a1bSJed Brown t37 = cos(t20); 388*c4762a1bSJed Brown t46 = VV * VV; 389*c4762a1bSJed Brown t53 = -t11 - t14 + 0.3e1 * t3 * t46 * Rm - 0.4e1 * t4 * VV * Rm; 390*c4762a1bSJed Brown t64 = -t26 + t33 + t3 * t46 * VV - 0.3e1 * t9 * VV - 0.2e1 * t4 * t46 - 0.3e1 * t3 * VV * t12; 391*c4762a1bSJed Brown t73 = kn * kn; 392*c4762a1bSJed Brown t74 = t73 * kn; 393*c4762a1bSJed Brown t79 = B * B; 394*c4762a1bSJed Brown t86 = B * t8; 395*c4762a1bSJed Brown t90 = kn * x; 396*c4762a1bSJed Brown t91 = sin(t90); 397*c4762a1bSJed Brown t106 = cos(t90); 398*c4762a1bSJed Brown *sxx = -((t19 * t21 + t36 * t37) * _C1 + (t36 * t21 - t19 * t37) * _C2 + (t53 * t21 + t64 * t37) * _C3 + (t64 * t21 - t53 * t37) * _C4 + (-AA * t74 - 0.4e1 * BB * t73 * B + 0.4e1 * t79 * AA * kn - 0.3e1 * t8 * AA * kn - 0.8e1 * t86 * BB) * t91 + (-0.8e1 * t86 * AA - 0.4e1 * AA * t73 * B - 0.4e1 * t79 * BB * kn + 0.3e1 * t8 * BB * kn + BB * t74) * t106) / km; 399*c4762a1bSJed Brown 400*c4762a1bSJed Brown t3 = exp(0.2e1 * x * B); 401*c4762a1bSJed Brown t4 = km * km; 402*c4762a1bSJed Brown t5 = t3 * t4; 403*c4762a1bSJed Brown t6 = Rm * x; 404*c4762a1bSJed Brown t7 = cos(t6); 405*c4762a1bSJed Brown t8 = _C1 * t7; 406*c4762a1bSJed Brown t10 = sin(t6); 407*c4762a1bSJed Brown t11 = _C2 * t10; 408*c4762a1bSJed Brown t13 = _C3 * t7; 409*c4762a1bSJed Brown t15 = _C4 * t10; 410*c4762a1bSJed Brown t18 = kn * x; 411*c4762a1bSJed Brown t19 = cos(t18); 412*c4762a1bSJed Brown t22 = sin(t18); 413*c4762a1bSJed Brown t24 = UU * UU; 414*c4762a1bSJed Brown t25 = t3 * t24; 415*c4762a1bSJed Brown t28 = t3 * UU; 416*c4762a1bSJed Brown t38 = Rm * Rm; 417*c4762a1bSJed Brown t39 = t7 * t38; 418*c4762a1bSJed Brown t42 = t10 * t38; 419*c4762a1bSJed Brown t44 = t5 * t8 + t5 * t11 + t5 * t13 + t5 * t15 + t4 * AA * t19 + t4 * BB * t22 + t25 * t8 + t25 * t11 - 0.2e1 * t28 * _C1 * t10 * Rm + 0.2e1 * t28 * _C2 * t7 * Rm - t3 * _C1 * t39 - t3 * _C2 * t42; 420*c4762a1bSJed Brown t45 = VV * VV; 421*c4762a1bSJed Brown t46 = t3 * t45; 422*c4762a1bSJed Brown t49 = t3 * VV; 423*c4762a1bSJed Brown t62 = B * B; 424*c4762a1bSJed Brown t78 = kn * kn; 425*c4762a1bSJed Brown t82 = t46 * t13 + t46 * t15 + 0.2e1 * t49 * _C3 * t10 * Rm - 0.2e1 * t49 * _C4 * t7 * Rm - t3 * _C3 * t39 - t3 * _C4 * t42 + 0.4e1 * t62 * AA * t19 + 0.4e1 * t62 * BB * t22 + 0.4e1 * B * AA * t22 * kn - 0.4e1 * B * BB * t19 * kn - AA * t19 * t78 - BB * t22 * t78; 426*c4762a1bSJed Brown *sxz = t44 + t82; 427*c4762a1bSJed Brown 428*c4762a1bSJed Brown t3 = exp(0.2e1 * x * B); 429*c4762a1bSJed Brown t4 = t3 * B; 430*c4762a1bSJed Brown t8 = km * km; 431*c4762a1bSJed Brown t9 = t3 * t8; 432*c4762a1bSJed Brown t10 = t9 * Rm; 433*c4762a1bSJed Brown t11 = Rm * Rm; 434*c4762a1bSJed Brown t13 = t3 * t11 * Rm; 435*c4762a1bSJed Brown t14 = UU * UU; 436*c4762a1bSJed Brown t18 = 0.4e1 * t4 * UU * Rm - t10 - t13 + 0.3e1 * t3 * t14 * Rm; 437*c4762a1bSJed Brown t19 = Rm * x; 438*c4762a1bSJed Brown t20 = sin(t19); 439*c4762a1bSJed Brown t25 = 0.2e1 * t9 * B; 440*c4762a1bSJed Brown t31 = 0.2e1 * t4 * t11; 441*c4762a1bSJed Brown t34 = -t3 * t14 * UU - t25 + t9 * UU + 0.3e1 * t3 * UU * t11 + t31 - 0.2e1 * t4 * t14; 442*c4762a1bSJed Brown t35 = cos(t19); 443*c4762a1bSJed Brown t44 = VV * VV; 444*c4762a1bSJed Brown t51 = -t10 - t13 + 0.3e1 * t3 * t44 * Rm - 0.4e1 * t4 * VV * Rm; 445*c4762a1bSJed Brown t61 = -t25 + t31 + t3 * t44 * VV - t9 * VV - 0.2e1 * t4 * t44 - 0.3e1 * t3 * VV * t11; 446*c4762a1bSJed Brown t70 = kn * kn; 447*c4762a1bSJed Brown t71 = t70 * kn; 448*c4762a1bSJed Brown t76 = B * B; 449*c4762a1bSJed Brown t82 = B * t8; 450*c4762a1bSJed Brown t86 = kn * x; 451*c4762a1bSJed Brown t87 = sin(t86); 452*c4762a1bSJed Brown t101 = cos(t86); 453*c4762a1bSJed Brown *p = ((t18 * t20 + t34 * t35) * _C1 + (t34 * t20 - t18 * t35) * _C2 + (t51 * t20 + t61 * t35) * _C3 + (t61 * t20 - t51 * t35) * _C4 + (-AA * t71 - 0.4e1 * BB * t70 * B + 0.4e1 * t76 * AA * kn - t8 * AA * kn - 0.4e1 * t82 * BB) * t87 + (-0.4e1 * t82 * AA - 0.4e1 * AA * t70 * B - 0.4e1 * t76 * BB * kn + t8 * BB * kn + BB * t71) * t101) / km; 454*c4762a1bSJed Brown 455*c4762a1bSJed Brown t3 = exp(0.2e1 * x * B); 456*c4762a1bSJed Brown t4 = UU * UU; 457*c4762a1bSJed Brown t8 = km * km; 458*c4762a1bSJed Brown t9 = t3 * t8; 459*c4762a1bSJed Brown t10 = t9 * Rm; 460*c4762a1bSJed Brown t11 = Rm * Rm; 461*c4762a1bSJed Brown t13 = t3 * t11 * Rm; 462*c4762a1bSJed Brown t14 = t3 * B; 463*c4762a1bSJed Brown t18 = 0.3e1 * t3 * t4 * Rm + t10 - t13 + 0.4e1 * t14 * UU * Rm; 464*c4762a1bSJed Brown t19 = Rm * x; 465*c4762a1bSJed Brown t20 = sin(t19); 466*c4762a1bSJed Brown t23 = 0.2e1 * t9 * B; 467*c4762a1bSJed Brown t33 = 0.2e1 * t14 * t11; 468*c4762a1bSJed Brown t34 = -t23 + 0.3e1 * t3 * UU * t11 - t9 * UU - t3 * t4 * UU - 0.2e1 * t4 * t14 + t33; 469*c4762a1bSJed Brown t35 = cos(t19); 470*c4762a1bSJed Brown t47 = VV * VV; 471*c4762a1bSJed Brown t51 = t10 - 0.4e1 * t14 * VV * Rm + 0.3e1 * t3 * t47 * Rm - t13; 472*c4762a1bSJed Brown t61 = t9 * VV - t23 + t3 * t47 * VV - 0.2e1 * t14 * t47 + t33 - 0.3e1 * t3 * VV * t11; 473*c4762a1bSJed Brown t70 = B * B; 474*c4762a1bSJed Brown t74 = kn * kn; 475*c4762a1bSJed Brown t75 = t74 * kn; 476*c4762a1bSJed Brown t83 = kn * x; 477*c4762a1bSJed Brown t84 = sin(t83); 478*c4762a1bSJed Brown t96 = cos(t83); 479*c4762a1bSJed Brown *szz = -((t18 * t20 + t34 * t35) * _C1 + (t34 * t20 - t18 * t35) * _C2 + (t51 * t20 + t61 * t35) * _C3 + (t61 * t20 - t51 * t35) * _C4 + (0.4e1 * t70 * AA * kn - AA * t75 - 0.4e1 * BB * t74 * B + t8 * AA * kn) * t84 + (-t8 * BB * kn - 0.4e1 * AA * t74 * B - 0.4e1 * t70 * BB * kn + BB * t75) * t96) / km; 480*c4762a1bSJed Brown 481*c4762a1bSJed Brown /* vx = Vx, vz = Vz, sxx = xx-component of stress tensor, sxz = xz-component of stress tensor, p = pressure, szz = zz-component of stress tensor */ 482*c4762a1bSJed Brown *vx *= cos(km*z); /* Vx */ 483*c4762a1bSJed Brown *vz *= sin(km*z); /* Vz */ 484*c4762a1bSJed Brown *p *= cos(km*z); /* p */ 485*c4762a1bSJed Brown *sxx *= cos(km*z); /* sxx total stress */ 486*c4762a1bSJed Brown *sxz *= sin(km*z); /* tzx stress */ 487*c4762a1bSJed Brown *szz *= cos(km*z); /* szz total stress */ 488*c4762a1bSJed Brown 489*c4762a1bSJed Brown /* rho = -sigma*sin(km*z)*cos(kn*x); */ /* density */ 490*c4762a1bSJed Brown PetscFunctionReturn(0); 491*c4762a1bSJed Brown } 492*c4762a1bSJed Brown 493*c4762a1bSJed Brown PetscErrorCode SolKxWrapperV(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar v[], void *ctx) 494*c4762a1bSJed Brown { 495*c4762a1bSJed Brown PetscReal B = 100.0; 496*c4762a1bSJed Brown PetscReal kn = 100*M_PI; 497*c4762a1bSJed Brown PetscReal km = 100*M_PI; 498*c4762a1bSJed Brown PetscScalar p, sxx, sxz, szz; 499*c4762a1bSJed Brown 500*c4762a1bSJed Brown PetscFunctionBeginUser; 501*c4762a1bSJed Brown SolKxSolution(x[0], x[1], kn, km, B, &v[0], &v[1], &p, &sxx, &sxz, &szz); 502*c4762a1bSJed Brown PetscFunctionReturn(0); 503*c4762a1bSJed Brown } 504*c4762a1bSJed Brown 505*c4762a1bSJed Brown PetscErrorCode SolKxWrapperP(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar v[], void *ctx) 506*c4762a1bSJed Brown { 507*c4762a1bSJed Brown PetscReal B = 100.0; 508*c4762a1bSJed Brown PetscReal kn = 100*M_PI; 509*c4762a1bSJed Brown PetscReal km = 100*M_PI; 510*c4762a1bSJed Brown PetscScalar vx, vz, sxx, sxz, szz; 511*c4762a1bSJed Brown 512*c4762a1bSJed Brown PetscFunctionBeginUser; 513*c4762a1bSJed Brown SolKxSolution(x[0], x[1], kn, km, B, &vx, &vz, &v[0], &sxx, &sxz, &szz); 514*c4762a1bSJed Brown PetscFunctionReturn(0); 515*c4762a1bSJed Brown } 516*c4762a1bSJed Brown 517*c4762a1bSJed Brown /* 518*c4762a1bSJed Brown Compare the C implementation with generated data from Maple 519*c4762a1bSJed Brown */ 520*c4762a1bSJed Brown PetscErrorCode MapleTest(MPI_Comm comm, AppCtx *ctx) 521*c4762a1bSJed Brown { 522*c4762a1bSJed Brown const PetscInt n = 41; 523*c4762a1bSJed Brown PetscScalar vxMaple[41][41], vzMaple[41][41], pMaple[41][41], sxxMaple[41][41], sxzMaple[41][41], szzMaple[41][41]; 524*c4762a1bSJed Brown PetscReal x[41], z[41]; 525*c4762a1bSJed Brown PetscReal kn, km, B; 526*c4762a1bSJed Brown PetscInt i, j; 527*c4762a1bSJed Brown PetscErrorCode ierr; 528*c4762a1bSJed Brown 529*c4762a1bSJed Brown PetscFunctionBegin; 530*c4762a1bSJed Brown ierr = SolKxData5(x, z, &kn, &km, &B, vxMaple, vzMaple, pMaple, sxxMaple, sxzMaple, szzMaple);CHKERRQ(ierr); 531*c4762a1bSJed Brown for (i = 0; i < n; ++i) { 532*c4762a1bSJed Brown for (j = 0; j < n; ++j) { 533*c4762a1bSJed Brown PetscScalar vx, vz, p, sxx, sxz, szz; 534*c4762a1bSJed Brown PetscReal norm; 535*c4762a1bSJed Brown 536*c4762a1bSJed Brown ierr = SolKxSolution(x[i], z[j], kn, km, B, &vx, &vz, &p, &sxx, &sxz, &szz);CHKERRQ(ierr); 537*c4762a1bSJed Brown norm = sqrt(PetscSqr(PetscAbsScalar(vx - vxMaple[i][j])) + PetscSqr(PetscAbsScalar(vz - vzMaple[i][j]))); 538*c4762a1bSJed Brown if (norm > 1.0e-10) { 539*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_SELF, "%0.17e %0.17e %0.17e %0.17e %0.17e %0.17e %0.17e %0.17e %0.17e\n", 540*c4762a1bSJed Brown (double)x[i], (double)z[j], (double)PetscAbsScalar(vx - vxMaple[i][j]), (double)PetscAbsScalar(vz - vzMaple[i][j]), (double)PetscAbsScalar(p - pMaple[i][j]), 541*c4762a1bSJed Brown (double)PetscAbsScalar(sxx - sxxMaple[i][j]), (double)PetscAbsScalar(sxz - sxzMaple[i][j]), (double)PetscAbsScalar(szz - szzMaple[i][j]), (double)norm); 542*c4762a1bSJed Brown SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid solution, error %g", (double)norm); 543*c4762a1bSJed Brown } 544*c4762a1bSJed Brown } 545*c4762a1bSJed Brown } 546*c4762a1bSJed Brown ierr = PetscPrintf(comm, "Verified Maple test 5\n");CHKERRQ(ierr); 547*c4762a1bSJed Brown PetscFunctionReturn(0); 548*c4762a1bSJed Brown } 549*c4762a1bSJed Brown 550*c4762a1bSJed Brown PetscErrorCode FEMTest(MPI_Comm comm, AppCtx *ctx) 551*c4762a1bSJed Brown { 552*c4762a1bSJed Brown DM dm; 553*c4762a1bSJed Brown Vec u; 554*c4762a1bSJed Brown PetscErrorCode (*funcs[2])(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *) = {SolKxWrapperV, SolKxWrapperP}; 555*c4762a1bSJed Brown PetscReal discError; 556*c4762a1bSJed Brown PetscErrorCode ierr; 557*c4762a1bSJed Brown 558*c4762a1bSJed Brown PetscFunctionBegin; 559*c4762a1bSJed Brown if (!ctx->fem) PetscFunctionReturn(0); 560*c4762a1bSJed Brown /* Create DM */ 561*c4762a1bSJed Brown ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_FALSE, &dm);CHKERRQ(ierr); 562*c4762a1bSJed Brown ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 563*c4762a1bSJed Brown /* Project solution into FE space */ 564*c4762a1bSJed Brown ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 565*c4762a1bSJed Brown ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, u);CHKERRQ(ierr); 566*c4762a1bSJed Brown ierr = DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &discError);CHKERRQ(ierr); 567*c4762a1bSJed Brown ierr = VecViewFromOptions(u, NULL, "-vec_view");CHKERRQ(ierr); 568*c4762a1bSJed Brown /* Cleanup */ 569*c4762a1bSJed Brown ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 570*c4762a1bSJed Brown ierr = DMDestroy(&dm);CHKERRQ(ierr); 571*c4762a1bSJed Brown PetscFunctionReturn(0); 572*c4762a1bSJed Brown } 573*c4762a1bSJed Brown 574*c4762a1bSJed Brown int main(int argc, char **argv) 575*c4762a1bSJed Brown { 576*c4762a1bSJed Brown AppCtx user; /* user-defined work context */ 577*c4762a1bSJed Brown PetscErrorCode ierr; 578*c4762a1bSJed Brown 579*c4762a1bSJed Brown ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 580*c4762a1bSJed Brown ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 581*c4762a1bSJed Brown ierr = MapleTest(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 582*c4762a1bSJed Brown ierr = FEMTest(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 583*c4762a1bSJed Brown ierr = PetscFinalize(); 584*c4762a1bSJed Brown return ierr; 585*c4762a1bSJed Brown } 586