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