Lines Matching refs:n
7 "# libCEED for Python examples\n",
8 "\n",
9 …oject.org/) (CEED) of the [Exascale Computing Project](https://www.exascaleproject.org/) (ECP).\n",
10 "\n",
18 "## Setting up libCEED for Python\n",
19 "\n",
36 "## CeedBasis\n",
37 "\n",
54 "%matplotlib inline\n",
55 "import numpy as np\n",
56 "import matplotlib.pyplot as plt\n",
57 "plt.style.use('ggplot')\n",
58 "\n",
59 "def eval(dim, x):\n",
60 " result, center = 1, 0.1\n",
61 " for d in range(dim):\n",
62 " result *= np.tanh(x[d] - center)\n",
63 " center += 0.1\n",
64 " return result\n",
65 "\n",
66 "def feval(x_1, x_2):\n",
67 " return x_1*x_1 + x_2*x_2 + x_1*x_2 + 1\n",
68 "\n",
69 "def dfeval(x_1, x_2):\n",
77 "## $H^1$ Lagrange bases in 1D\n",
78 "\n",
88 "import libceed\n",
89 "\n",
90 "ceed = libceed.Ceed()\n",
91 "\n",
92 "b = ceed.BasisTensorH1Lagrange(\n",
93 " dim=1, # topological dimension\n",
94 " ncomp=1, # number of components\n",
95 " P=4, # number of basis functions (nodes) per dimension\n",
96 " Q=4, # number of quadrature points per dimension\n",
97 " qmode=libceed.GAUSS_LOBATTO)\n",
114 "P = b.get_num_nodes()\n",
115 "Q_viz = 50\n",
116 "basis_viz = ceed.BasisTensorH1Lagrange(1, 1, P, Q_viz, libceed.GAUSS_LOBATTO)\n",
117 "\n",
118 "# Construct P \"elements\" with one node activated\n",
119 "I = ceed.Vector(P * P)\n",
120 "with I.array_write(P, P) as x:\n",
121 " x[...] = np.eye(P)\n",
122 "\n",
123 "basis_fns = ceed.Vector(P * Q_viz)\n",
124 "basis_viz.apply(4, libceed.EVAL_INTERP, I, basis_fns)\n",
125 "\n",
126 "qpts_viz, _ = ceed.lobatto_quadrature(Q_viz)\n",
127 "with basis_fns.array_read(Q_viz, P) as B_array:\n",
128 " plt.plot(qpts_viz, B_array)\n",
129 "\n",
130 "# Mark tho Lobatto nodes\n",
131 "nodes, _ = ceed.lobatto_quadrature(P)\n",
148 "b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS)\n",
149 "print(b)\n",
150 "\n",
151 "with basis_fns.array_read(Q_viz, P) as B_array:\n",
152 " plt.plot(qpts_viz, B_array)\n",
153 "# Mark tho Gauss quadrature points\n",
154 "qpts, _ = ceed.gauss_quadrature(P)\n",
162 …erlying functions are not an intrinsic property of a `libceed.Basis` in libCEED, the sizes are.\n",
172 "b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO)\n",
173 "\n",
174 "p = b.get_num_nodes()\n",
175 "print('p =', p)\n",
176 "\n",
177 "q = b.get_num_quadrature_points()\n",
194 "for dim in range(1, 4):\n",
195 " Q = 4\n",
196 " Q_dim = Q**dim\n",
197 " X_dim = 2**dim\n",
198 " x = np.empty(X_dim*dim, dtype=\"float64\")\n",
199 " u_array = np.empty(Q_dim, dtype=\"float64\")\n",
200 "\n",
201 " for d in range(dim):\n",
202 " for i in range(X_dim):\n",
203 " x[d*X_dim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n",
204 "\n",
205 " X = ceed.Vector(X_dim*dim)\n",
206 " X.set_array(x, cmode=libceed.USE_POINTER)\n",
207 " X_q = ceed.Vector(Q_dim*dim)\n",
208 " X_q.set_value(0)\n",
209 " U = ceed.Vector(Q_dim)\n",
210 " U.set_value(0)\n",
211 " U_q = ceed.Vector(Q_dim)\n",
212 "\n",
213 " basis_x_lobatto = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO)\n",
214 " basis_u_lobatto = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO)\n",
215 "\n",
216 " basis_x_lobatto.apply(1, libceed.EVAL_INTERP, X, X_q)\n",
217 "\n",
218 " with X_q.array_read() as x_array:\n",
219 " for i in range(Q_dim):\n",
220 " x = np.empty(dim, dtype=\"float64\")\n",
221 " for d in range(dim):\n",
222 " x[d] = x_array[d*Q_dim + i]\n",
223 " u_array[i] = eval(dim, x)\n",
224 "\n",
225 " U_q.set_array(u_array, cmode=libceed.USE_POINTER)\n",
226 "\n",
227 " # This operation is the identity because the quadrature is collocated\n",
228 " basis_u_lobatto.T.apply(1, libceed.EVAL_INTERP, U_q, U)\n",
229 "\n",
230 " basis_x_gauss = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)\n",
231 " basis_u_gauss = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)\n",
232 "\n",
233 " basis_x_gauss.apply(1, libceed.EVAL_INTERP, X, X_q)\n",
234 " basis_u_gauss.apply(1, libceed.EVAL_INTERP, U, U_q)\n",
235 "\n",
236 " with X_q.array_read() as x_array, U_q.array_read() as u_array:\n",
237 " if dim == 2:\n",
238 " # Default ordering is contiguous in x direction, but\n",
239 " # pyplot expects meshgrid convention, which is transposed.\n",
240 " x, y = x_array.reshape(2, Q, Q).transpose(0, 2, 1)\n",
241 " plt.scatter(x, y, c=np.array(u_array).reshape(Q, Q))\n",
242 " plt.xlim(-1, 1)\n",
243 " plt.ylim(-1, 1)\n",
260 "for dim in range (1, 4):\n",
261 " P, Q = 8, 10\n",
262 " P_dim = P**dim\n",
263 " Q_dim = Q**dim\n",
264 " X_dim = 2**dim\n",
265 " sum_1 = sum_2 = 0\n",
266 " x_array = np.empty(X_dim*dim, dtype=\"float64\")\n",
267 " u_array = np.empty(P_dim, dtype=\"float64\")\n",
268 "\n",
269 " for d in range(dim):\n",
270 " for i in range(X_dim):\n",
271 " x_array[d*X_dim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n",
272 "\n",
273 " X = ceed.Vector(X_dim*dim)\n",
274 " X.set_array(x_array, cmode=libceed.USE_POINTER)\n",
275 " X_q = ceed.Vector(P_dim*dim)\n",
276 " X_q.set_value(0)\n",
277 " U = ceed.Vector(P_dim)\n",
278 " U_q = ceed.Vector(Q_dim*dim)\n",
279 " U_q.set_value(0)\n",
280 " Ones = ceed.Vector(Q_dim*dim)\n",
281 " Ones.set_value(1)\n",
282 " G_transpose_ones = ceed.Vector(P_dim)\n",
283 " G_transpose_ones.set_value(0)\n",
284 "\n",
285 " # Get function values at quadrature points\n",
286 " basis_x_lobatto = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO)\n",
287 " basis_x_lobatto.apply(1, libceed.EVAL_INTERP, X, X_q)\n",
288 "\n",
289 " with X_q.array_read() as x_array:\n",
290 " for i in range(P_dim):\n",
291 " x = np.empty(dim, dtype=\"float64\")\n",
292 " for d in range(dim):\n",
293 " x[d] = x_array[d*P_dim + i]\n",
294 " u_array[i] = eval(dim, x)\n",
295 "\n",
296 " U.set_array(u_array, cmode=libceed.USE_POINTER)\n",
297 "\n",
298 " # Calculate G u at quadrature points, G' * 1 at dofs\n",
299 " basis_u_gauss = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS)\n",
300 " basis_u_gauss.apply(1, libceed.EVAL_GRAD, U, U_q)\n",
301 " basis_u_gauss.T.apply(1, libceed.EVAL_GRAD, Ones, G_transpose_ones)\n",
302 "\n",
303 " # Check if 1' * G * u = u' * (G' * 1)\n",
304 " with G_transpose_ones.array_read() as g_array, U_q.array_read() as uq_array:\n",
305 " for i in range(P_dim):\n",
306 " sum_1 += g_array[i]*u_array[i]\n",
307 " for i in range(dim*Q_dim):\n",
308 " sum_2 += uq_array[i]\n",
309 "\n",
310 " # Check that (1' * G * u - u' * (G' * 1)) is numerically zero\n",