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