1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7 //
8 // libCEED Example 2
9 //
10 // This example illustrates a simple usage of libCEED to compute the surface
11 // area of a 3D body using matrix-free application of a diffusion operator.
12 // Arbitrary mesh and solution degrees in 1D, 2D and 3D are supported from the
13 // same code.
14 //
15 // The example has no dependencies, and is designed to be self-contained. For
16 // additional examples that use external discretization libraries (MFEM, PETSc,
17 // etc.) see the subdirectories in libceed/examples.
18 //
19 // All libCEED objects use a Ceed device object constructed based on a command
20 // line argument (-ceed).
21
22 use clap::Parser;
23 use libceed::{
24 BasisOpt, Ceed, ElemRestrictionOpt, QFunctionInputs, QFunctionOpt, QFunctionOutputs, VectorOpt,
25 };
26 mod opt;
27 mod transform;
28
29 // ----------------------------------------------------------------------------
30 // Example 2
31 // ----------------------------------------------------------------------------
main() -> libceed::Result<()>32 fn main() -> libceed::Result<()> {
33 let options = opt::Opt::parse();
34 example_2(options)
35 }
36
37 #[allow(clippy::erasing_op)]
38 #[allow(clippy::identity_op)]
example_2(options: opt::Opt) -> libceed::Result<()>39 fn example_2(options: opt::Opt) -> libceed::Result<()> {
40 // Process command line arguments
41 let opt::Opt {
42 ceed_spec,
43 dim,
44 mesh_degree,
45 solution_degree,
46 num_qpts,
47 problem_size_requested,
48 test,
49 quiet,
50 gallery,
51 } = options;
52 assert!((0..=3).contains(&dim));
53 assert!(mesh_degree >= 1);
54 assert!(solution_degree >= 1);
55 assert!(num_qpts >= 1);
56 let ncomp_x = dim;
57 let problem_size: i64 = if problem_size_requested < 0 {
58 if test {
59 16 * 16 * (dim * dim) as i64
60 } else {
61 256 * 1024
62 }
63 } else {
64 problem_size_requested
65 };
66
67 // Summary output
68 if !quiet {
69 println!("Selected options: [command line option] : <current value>");
70 println!(" Ceed specification [-c] : {}", ceed_spec);
71 println!(" Mesh dimension [-d] : {}", dim);
72 println!(" Mesh degree [-m] : {}", mesh_degree);
73 println!(" Solution degree [-p] : {}", solution_degree);
74 println!(" Num. 1D quadr. pts [-q] : {}", num_qpts);
75 println!(" Approx. # unknowns [-s] : {}", problem_size);
76 println!(
77 " QFunction source [-g] : {}",
78 if gallery { "gallery" } else { "user closure" }
79 );
80 }
81
82 // Initalize ceed context
83 let ceed = Ceed::init(&ceed_spec);
84
85 // Mesh and solution bases
86 let basis_mesh = ceed.basis_tensor_H1_Lagrange(
87 dim,
88 ncomp_x,
89 mesh_degree + 1,
90 num_qpts,
91 libceed::QuadMode::Gauss,
92 )?;
93 let basis_solution = ceed.basis_tensor_H1_Lagrange(
94 dim,
95 1,
96 solution_degree + 1,
97 num_qpts,
98 libceed::QuadMode::Gauss,
99 )?;
100
101 // Determine mesh size from approximate problem size
102 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size);
103 if !quiet {
104 print!("\nMesh size : nx = {}", num_xyz[0]);
105 if dim > 1 {
106 print!(", ny = {}", num_xyz[1]);
107 }
108 if dim > 2 {
109 print!(", nz = {}", num_xyz[2]);
110 }
111 println!();
112 }
113
114 // Build ElemRestriction objects describing the mesh and solution discrete
115 // representations
116 let (rstr_mesh, _) =
117 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?;
118 let (_, rstr_qdata) = mesh::build_cartesian_restriction(
119 &ceed,
120 dim,
121 num_xyz,
122 solution_degree,
123 dim * (dim + 1) / 2,
124 num_qpts,
125 )?;
126 let (rstr_solution, _) =
127 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?;
128 let mesh_size = rstr_mesh.lvector_size();
129 let solution_size = rstr_solution.lvector_size();
130 if !quiet {
131 println!("Number of mesh nodes : {}", mesh_size / dim);
132 println!("Number of solution nodes : {}", solution_size);
133 }
134
135 // Create a Vector with the mesh coordinates
136 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?;
137
138 // Apply a transformation to the mesh coordinates
139 let exact_area = transform::transform_mesh_coordinates(dim, &mut mesh_coords)?;
140
141 // QFunction that builds the quadrature data for the diff operator
142 // -- QFunction from user closure
143 let build_diff = move |[jacobian, weights, ..]: QFunctionInputs,
144 [qdata, ..]: QFunctionOutputs| {
145 // Build quadrature data
146 match dim {
147 1 => qdata
148 .iter_mut()
149 .zip(jacobian.iter().zip(weights.iter()))
150 .for_each(|(qdata, (j, weight))| *qdata = weight / j),
151 2 => {
152 let q = qdata.len() / 3;
153 for i in 0..q {
154 let j11 = jacobian[i + q * 0];
155 let j21 = jacobian[i + q * 1];
156 let j12 = jacobian[i + q * 2];
157 let j22 = jacobian[i + q * 3];
158 let qw = weights[i] / (j11 * j22 - j21 * j12);
159 qdata[i + q * 0] = qw * (j12 * j12 + j22 * j22);
160 qdata[i + q * 1] = qw * (j11 * j11 + j21 * j21);
161 qdata[i + q * 2] = -qw * (j11 * j12 + j21 * j22);
162 }
163 }
164 3 => {
165 let q = qdata.len() / 6;
166 for i in 0..q {
167 let mut a = [0.0; 9];
168 for j in 0..3 {
169 for k in 0..3 {
170 a[k * 3 + j] = jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 1) % 3))]
171 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 2) % 3))]
172 - jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 2) % 3))]
173 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 1) % 3))];
174 }
175 }
176 let qw = weights[i]
177 / (jacobian[i + q * 0] * a[0 * 3 + 0]
178 + jacobian[i + q * 1] * a[0 * 3 + 1]
179 + jacobian[i + q * 2] * a[0 * 3 + 2]);
180 qdata[i + q * 0] = qw
181 * (a[0 * 3 + 0] * a[0 * 3 + 0]
182 + a[0 * 3 + 1] * a[0 * 3 + 1]
183 + a[0 * 3 + 2] * a[0 * 3 + 2]);
184 qdata[i + q * 1] = qw
185 * (a[1 * 3 + 0] * a[1 * 3 + 0]
186 + a[1 * 3 + 1] * a[1 * 3 + 1]
187 + a[1 * 3 + 2] * a[1 * 3 + 2]);
188 qdata[i + q * 2] = qw
189 * (a[2 * 3 + 0] * a[2 * 3 + 0]
190 + a[2 * 3 + 1] * a[2 * 3 + 1]
191 + a[2 * 3 + 2] * a[2 * 3 + 2]);
192 qdata[i + q * 3] = qw
193 * (a[1 * 3 + 0] * a[2 * 3 + 0]
194 + a[1 * 3 + 1] * a[2 * 3 + 1]
195 + a[1 * 3 + 2] * a[2 * 3 + 2]);
196 qdata[i + q * 4] = qw
197 * (a[0 * 3 + 0] * a[2 * 3 + 0]
198 + a[0 * 3 + 1] * a[2 * 3 + 1]
199 + a[0 * 3 + 2] * a[2 * 3 + 2]);
200 qdata[i + q * 5] = qw
201 * (a[0 * 3 + 0] * a[1 * 3 + 0]
202 + a[0 * 3 + 1] * a[1 * 3 + 1]
203 + a[0 * 3 + 2] * a[1 * 3 + 2]);
204 }
205 }
206 _ => unreachable!(),
207 };
208
209 // Return clean error code
210 0
211 };
212 let qf_build_closure = ceed
213 .q_function_interior(1, Box::new(build_diff))?
214 .input("dx", ncomp_x * dim, libceed::EvalMode::Grad)?
215 .input("weights", 1, libceed::EvalMode::Weight)?
216 .output("qdata", dim * (dim + 1) / 2, libceed::EvalMode::None)?;
217 // -- QFunction from gallery
218 let qf_build_named = {
219 let name = format!("Poisson{}DBuild", dim);
220 ceed.q_function_interior_by_name(&name)?
221 };
222 // -- QFunction for use with Operator
223 let qf_build = if gallery {
224 QFunctionOpt::SomeQFunctionByName(&qf_build_named)
225 } else {
226 QFunctionOpt::SomeQFunction(&qf_build_closure)
227 };
228
229 // Operator that build the quadrature data for the diff operator
230 let op_build = ceed
231 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)?
232 .name("build qdata")?
233 .field("dx", &rstr_mesh, &basis_mesh, VectorOpt::Active)?
234 .field(
235 "weights",
236 ElemRestrictionOpt::None,
237 &basis_mesh,
238 VectorOpt::None,
239 )?
240 .field("qdata", &rstr_qdata, BasisOpt::None, VectorOpt::Active)?
241 .check()?;
242
243 // Compute the quadrature data for the diff operator
244 let elem_qpts = num_qpts.pow(dim as u32);
245 let num_elem: usize = num_xyz.iter().take(dim).product();
246 let mut qdata = ceed.vector(num_elem * elem_qpts * dim * (dim + 1) / 2)?;
247 op_build.apply(&mesh_coords, &mut qdata)?;
248
249 // QFunction that applies the diff operator
250 // -- QFunction from user closure
251 let apply_diff = move |[ug, qdata, ..]: QFunctionInputs, [vg, ..]: QFunctionOutputs| {
252 // Apply diffusion operator
253 match dim {
254 1 => vg
255 .iter_mut()
256 .zip(ug.iter().zip(qdata.iter()))
257 .for_each(|(vg, (ug, w))| *vg = ug * w),
258 2 => {
259 let q = qdata.len() / 3;
260 for i in 0..q {
261 let du = [ug[i + q * 0], ug[i + q * 1]];
262 let dxdxdxdx_t = [
263 [qdata[i + 0 * q], qdata[i + 2 * q]],
264 [qdata[i + 2 * q], qdata[i + 1 * q]],
265 ];
266 for j in 0..2 {
267 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j];
268 }
269 }
270 }
271 3 => {
272 let q = qdata.len() / 6;
273 for i in 0..q {
274 let du = [ug[i + q * 0], ug[i + q * 1], ug[i + q * 2]];
275 let dxdxdxdx_t = [
276 [qdata[i + 0 * q], qdata[i + 5 * q], qdata[i + 4 * q]],
277 [qdata[i + 5 * q], qdata[i + 1 * q], qdata[i + 3 * q]],
278 [qdata[i + 4 * q], qdata[i + 3 * q], qdata[i + 2 * q]],
279 ];
280 for j in 0..3 {
281 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j]
282 + du[1] * dxdxdxdx_t[1][j]
283 + du[2] * dxdxdxdx_t[2][j];
284 }
285 }
286 }
287 _ => unreachable!(),
288 };
289
290 // Return clean error code
291 0
292 };
293 let qf_diff_closure = ceed
294 .q_function_interior(1, Box::new(apply_diff))?
295 .input("du", dim, libceed::EvalMode::Grad)?
296 .input("qdata", dim * (dim + 1) / 2, libceed::EvalMode::None)?
297 .output("dv", dim, libceed::EvalMode::Grad)?;
298 // -- QFunction from gallery
299 let qf_diff_named = {
300 let name = format!("Poisson{}DApply", dim);
301 ceed.q_function_interior_by_name(&name)?
302 };
303 // -- QFunction for use with Operator
304 let qf_diff = if gallery {
305 QFunctionOpt::SomeQFunctionByName(&qf_diff_named)
306 } else {
307 QFunctionOpt::SomeQFunction(&qf_diff_closure)
308 };
309
310 // Diff Operator
311 let op_diff = ceed
312 .operator(qf_diff, QFunctionOpt::None, QFunctionOpt::None)?
313 .name("Poisson")?
314 .field("du", &rstr_solution, &basis_solution, VectorOpt::Active)?
315 .field("qdata", &rstr_qdata, BasisOpt::None, &qdata)?
316 .field("dv", &rstr_solution, &basis_solution, VectorOpt::Active)?
317 .check()?;
318
319 // Solution vectors
320 let mut u = ceed.vector(solution_size)?;
321 let mut v = ceed.vector(solution_size)?;
322
323 // Initialize u with sum of node coordinates
324 let coords = mesh_coords.view()?;
325 u.set_value(0.0)?;
326 for (i, u) in u.view_mut()?.iter_mut().enumerate() {
327 *u = (0..dim).map(|d| coords[i + d * solution_size]).sum();
328 }
329
330 // Apply the diff operator
331 op_diff.apply(&u, &mut v)?;
332
333 // Compute the mesh surface area
334 let area: libceed::Scalar = v.view()?.iter().map(|v| (*v).abs()).sum();
335
336 // Output results
337 if !quiet {
338 println!("Exact mesh surface area : {:.12}", exact_area);
339 println!("Computed mesh surface_area : {:.12}", area);
340 println!("Surface area error : {:.12e}", area - exact_area);
341 }
342 let tolerance = match dim {
343 1 => 10000.0 * libceed::EPSILON,
344 _ => 1E-1,
345 };
346 let error = (area - exact_area).abs();
347 if error > tolerance {
348 println!("Volume error too large: {:.12e}", error);
349 return Err(libceed::Error {
350 message: format!(
351 "Volume error too large - expected: {:.12e}, actual: {:.12e}",
352 tolerance, error
353 ),
354 });
355 }
356 Ok(())
357 }
358
359 // ----------------------------------------------------------------------------
360 // Tests
361 // ----------------------------------------------------------------------------
362 #[cfg(test)]
363 mod tests {
364 use super::*;
365
366 #[test]
example_2_1d()367 fn example_2_1d() {
368 let options = opt::Opt {
369 ceed_spec: "/cpu/self/ref/serial".to_string(),
370 dim: 1,
371 mesh_degree: 4,
372 solution_degree: 4,
373 num_qpts: 6,
374 problem_size_requested: -1,
375 test: true,
376 quiet: true,
377 gallery: false,
378 };
379 assert!(example_2(options).is_ok());
380 }
381
382 #[test]
example_2_2d()383 fn example_2_2d() {
384 let options = opt::Opt {
385 ceed_spec: "/cpu/self/ref/serial".to_string(),
386 dim: 2,
387 mesh_degree: 4,
388 solution_degree: 4,
389 num_qpts: 6,
390 problem_size_requested: -1,
391 test: true,
392 quiet: true,
393 gallery: false,
394 };
395 assert!(example_2(options).is_ok());
396 }
397
398 #[test]
example_2_3d()399 fn example_2_3d() {
400 let options = opt::Opt {
401 ceed_spec: "/cpu/self/ref/serial".to_string(),
402 dim: 3,
403 mesh_degree: 4,
404 solution_degree: 4,
405 num_qpts: 6,
406 problem_size_requested: -1,
407 test: true,
408 quiet: false,
409 gallery: false,
410 };
411 assert!(example_2(options).is_ok());
412 }
413
414 #[test]
example_2_1d_gallery()415 fn example_2_1d_gallery() {
416 let options = opt::Opt {
417 ceed_spec: "/cpu/self/ref/serial".to_string(),
418 dim: 1,
419 mesh_degree: 4,
420 solution_degree: 4,
421 num_qpts: 6,
422 problem_size_requested: -1,
423 test: true,
424 quiet: true,
425 gallery: true,
426 };
427 assert!(example_2(options).is_ok());
428 }
429
430 #[test]
example_2_2d_gallery()431 fn example_2_2d_gallery() {
432 let options = opt::Opt {
433 ceed_spec: "/cpu/self/ref/serial".to_string(),
434 dim: 2,
435 mesh_degree: 4,
436 solution_degree: 4,
437 num_qpts: 6,
438 problem_size_requested: -1,
439 test: true,
440 quiet: true,
441 gallery: true,
442 };
443 assert!(example_2(options).is_ok());
444 }
445
446 #[test]
example_2_3d_gallery()447 fn example_2_3d_gallery() {
448 let options = opt::Opt {
449 ceed_spec: "/cpu/self/ref/serial".to_string(),
450 dim: 3,
451 mesh_degree: 4,
452 solution_degree: 4,
453 num_qpts: 6,
454 problem_size_requested: -1,
455 test: true,
456 quiet: true,
457 gallery: true,
458 };
459 assert!(example_2(options).is_ok());
460 }
461 }
462
463 // ----------------------------------------------------------------------------
464