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 3
9 //
10 // This example illustrates a simple usage of libCEED to compute the volume of a
11 // 3D body using matrix-free application of a mass operator. Arbitrary mesh and
12 // solution orders in 1D, 2D and 3D are supported from the same code. This
13 // calculation is executed in triplicate with a 3 component vector system.
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 3
31 // ----------------------------------------------------------------------------
main() -> libceed::Result<()>32 fn main() -> libceed::Result<()> {
33 let options = opt::Opt::parse();
34 example_1_vector(options)
35 }
36
37 #[allow(clippy::erasing_op)]
38 #[allow(clippy::identity_op)]
example_1_vector(options: opt::Opt) -> libceed::Result<()>39 fn example_1_vector(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 8 * 16
60 } else {
61 256 * 1024
62 }
63 } else {
64 problem_size_requested
65 };
66 let ncomp_u = 3;
67
68 // Summary output
69 if !quiet {
70 println!("Selected options: [command line option] : <current value>");
71 println!(" Ceed specification [-c] : {}", ceed_spec);
72 println!(" Mesh dimension [-d] : {}", dim);
73 println!(" Mesh degree [-m] : {}", mesh_degree);
74 println!(" Solution degree [-p] : {}", solution_degree);
75 println!(" Num. 1D quadr. pts [-q] : {}", num_qpts);
76 println!(" Approx. # unknowns [-s] : {}", problem_size);
77 println!(
78 " QFunction source [-g] : {}",
79 if gallery { "gallery" } else { "user closure" }
80 );
81 }
82
83 // Initalize ceed context
84 let ceed = Ceed::init(&ceed_spec);
85
86 // Mesh and solution bases
87 let basis_mesh = ceed.basis_tensor_H1_Lagrange(
88 dim,
89 ncomp_x,
90 mesh_degree + 1,
91 num_qpts,
92 libceed::QuadMode::Gauss,
93 )?;
94 let basis_solution = ceed.basis_tensor_H1_Lagrange(
95 dim,
96 ncomp_u,
97 solution_degree + 1,
98 num_qpts,
99 libceed::QuadMode::Gauss,
100 )?;
101
102 // Determine mesh size from approximate problem size
103 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size);
104 if !quiet {
105 print!("\nMesh size : nx = {}", num_xyz[0]);
106 if dim > 1 {
107 print!(", ny = {}", num_xyz[1]);
108 }
109 if dim > 2 {
110 print!(", nz = {}", num_xyz[2]);
111 }
112 println!();
113 }
114
115 // Build ElemRestriction objects describing the mesh and solution discrete
116 // representations
117 let (rstr_mesh, _) =
118 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?;
119 let (_, rstr_qdata) =
120 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?;
121 let (rstr_solution, _) =
122 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, ncomp_u, num_qpts)?;
123 let mesh_size = rstr_mesh.lvector_size();
124 let solution_size = rstr_solution.lvector_size();
125 if !quiet {
126 println!("Number of mesh nodes : {}", mesh_size / dim);
127 println!("Number of solution nodes : {}", solution_size);
128 }
129
130 // Create a Vector with the mesh coordinates
131 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?;
132
133 // Apply a transformation to the mesh coordinates
134 let exact_volume = transform::transform_mesh_coordinates(dim, mesh_size, &mut mesh_coords)?;
135
136 // QFunction that builds the quadrature data for the mass operator
137 // -- QFunction from user closure
138 let build_mass = move |[jacobian, weights, ..]: QFunctionInputs,
139 [qdata, ..]: QFunctionOutputs| {
140 // Build quadrature data
141 match dim {
142 1 => qdata
143 .iter_mut()
144 .zip(jacobian.iter().zip(weights.iter()))
145 .for_each(|(qdata, (j, weight))| *qdata = j * weight),
146 2 => {
147 let q = qdata.len();
148 qdata.iter_mut().zip(weights.iter()).enumerate().for_each(
149 |(i, (qdata, weight))| {
150 *qdata = (jacobian[i + q * 0] * jacobian[i + q * 3]
151 - jacobian[i + q * 1] * jacobian[i + q * 2])
152 * weight
153 },
154 );
155 }
156 3 => {
157 let q = qdata.len();
158 qdata.iter_mut().zip(weights.iter()).enumerate().for_each(
159 |(i, (qdata, weight))| {
160 *qdata = (jacobian[i + q * 0]
161 * (jacobian[i + q * 4] * jacobian[i + q * 8]
162 - jacobian[i + q * 5] * jacobian[i + q * 7])
163 - jacobian[i + q * 1]
164 * (jacobian[i + q * 3] * jacobian[i + q * 8]
165 - jacobian[i + q * 5] * jacobian[i + q * 6])
166 + jacobian[i + q * 2]
167 * (jacobian[i + q * 3] * jacobian[i + q * 7]
168 - jacobian[i + q * 4] * jacobian[i + q * 6]))
169 * weight
170 },
171 );
172 }
173 _ => unreachable!(),
174 };
175
176 // Return clean error code
177 0
178 };
179 let qf_build_closure = ceed
180 .q_function_interior(1, Box::new(build_mass))?
181 .input("dx", ncomp_x * dim, libceed::EvalMode::Grad)?
182 .input("weights", 1, libceed::EvalMode::Weight)?
183 .output("qdata", 1, libceed::EvalMode::None)?;
184 // -- QFunction from gallery
185 let qf_build_named = {
186 let name = format!("Mass{}DBuild", dim);
187 ceed.q_function_interior_by_name(&name)?
188 };
189 // -- QFunction for use with Operator
190 let qf_build = if gallery {
191 QFunctionOpt::SomeQFunctionByName(&qf_build_named)
192 } else {
193 QFunctionOpt::SomeQFunction(&qf_build_closure)
194 };
195
196 // Operator that build the quadrature data for the mass operator
197 let op_build = ceed
198 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)?
199 .name("build qdata")?
200 .field("dx", &rstr_mesh, &basis_mesh, VectorOpt::Active)?
201 .field(
202 "weights",
203 ElemRestrictionOpt::None,
204 &basis_mesh,
205 VectorOpt::None,
206 )?
207 .field("qdata", &rstr_qdata, BasisOpt::None, VectorOpt::Active)?
208 .check()?;
209
210 // Compute the quadrature data for the mass operator
211 let elem_qpts = num_qpts.pow(dim as u32);
212 let num_elem: usize = num_xyz.iter().take(dim).product();
213 let mut qdata = ceed.vector(num_elem * elem_qpts)?;
214 op_build.apply(&mesh_coords, &mut qdata)?;
215
216 // QFunction that applies the mass operator
217 // -- QFunction from user closure
218 let apply_mass = move |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| {
219 let q = qdata.len();
220 // Apply mass operator
221 for c in 0..ncomp_u {
222 v.iter_mut()
223 .skip(c * q)
224 .zip(u.iter().skip(c * q).zip(qdata.iter()))
225 .for_each(|(v, (u, w))| *v = u * w);
226 }
227 // Return clean error code
228 0
229 };
230 let qf_mass_closure = ceed
231 .q_function_interior(1, Box::new(apply_mass))?
232 .input("u", ncomp_u, libceed::EvalMode::Interp)?
233 .input("qdata", 1, libceed::EvalMode::None)?
234 .output("v", ncomp_u, libceed::EvalMode::Interp)?;
235 // -- QFunction from gallery
236 let qf_mass_named = ceed.q_function_interior_by_name("Vector3MassApply")?;
237 // -- QFunction for use with Operator
238 let qf_mass = if gallery {
239 QFunctionOpt::SomeQFunctionByName(&qf_mass_named)
240 } else {
241 QFunctionOpt::SomeQFunction(&qf_mass_closure)
242 };
243
244 // Mass Operator
245 let op_mass = ceed
246 .operator(qf_mass, QFunctionOpt::None, QFunctionOpt::None)?
247 .name("mass")?
248 .field("u", &rstr_solution, &basis_solution, VectorOpt::Active)?
249 .field("qdata", &rstr_qdata, BasisOpt::None, &qdata)?
250 .field("v", &rstr_solution, &basis_solution, VectorOpt::Active)?
251 .check()?;
252
253 // Solution vectors
254 let mut u = ceed.vector(solution_size)?;
255 let mut v = ceed.vector(solution_size)?;
256
257 // Initialize u with component index
258 u.set_value(0.0)?;
259 for c in 0..ncomp_u {
260 let q = solution_size / ncomp_u;
261 u.view_mut()?.iter_mut().skip(c * q).take(q).for_each(|u| {
262 *u = (c + 1) as libceed::Scalar;
263 });
264 }
265
266 // Apply the mass operator
267 op_mass.apply(&u, &mut v)?;
268
269 // Compute the mesh volume
270 let volume: libceed::Scalar = v.view()?.iter().sum::<libceed::Scalar>()
271 / ((ncomp_u * (ncomp_u + 1)) / 2) as libceed::Scalar;
272
273 // Output results
274 if !quiet {
275 println!("Exact mesh volume : {:.12}", exact_volume);
276 println!("Computed mesh volume : {:.12}", volume);
277 println!(
278 "Volume error : {:.12e}",
279 volume - exact_volume
280 );
281 }
282 let tolerance = match dim {
283 1 => 200.0 * libceed::EPSILON,
284 _ => 1E-5,
285 };
286 let error = (volume - exact_volume).abs();
287 if error > tolerance {
288 println!("Volume error too large: {:.12e}", error);
289 return Err(libceed::Error {
290 message: format!(
291 "Volume error too large - expected: {:.12e}, actual: {:.12e}",
292 tolerance, error
293 ),
294 });
295 }
296 Ok(())
297 }
298
299 // ----------------------------------------------------------------------------
300 // Tests
301 // ----------------------------------------------------------------------------
302 #[cfg(test)]
303 mod tests {
304 use super::*;
305
306 #[test]
example_1_vector_1d()307 fn example_1_vector_1d() {
308 let options = opt::Opt {
309 ceed_spec: "/cpu/self/ref/serial".to_string(),
310 dim: 1,
311 mesh_degree: 4,
312 solution_degree: 4,
313 num_qpts: 6,
314 problem_size_requested: -1,
315 test: true,
316 quiet: true,
317 gallery: false,
318 };
319 assert!(example_1_vector(options).is_ok());
320 }
321
322 #[test]
example_1_vector_2d()323 fn example_1_vector_2d() {
324 let options = opt::Opt {
325 ceed_spec: "/cpu/self/ref/serial".to_string(),
326 dim: 2,
327 mesh_degree: 4,
328 solution_degree: 4,
329 num_qpts: 6,
330 problem_size_requested: -1,
331 test: true,
332 quiet: true,
333 gallery: false,
334 };
335 assert!(example_1_vector(options).is_ok());
336 }
337
338 #[test]
example_1_vector_3d()339 fn example_1_vector_3d() {
340 let options = opt::Opt {
341 ceed_spec: "/cpu/self/ref/serial".to_string(),
342 dim: 3,
343 mesh_degree: 4,
344 solution_degree: 4,
345 num_qpts: 6,
346 problem_size_requested: -1,
347 test: true,
348 quiet: false,
349 gallery: false,
350 };
351 assert!(example_1_vector(options).is_ok());
352 }
353
354 #[test]
example_1_vector_1d_gallery()355 fn example_1_vector_1d_gallery() {
356 let options = opt::Opt {
357 ceed_spec: "/cpu/self/ref/serial".to_string(),
358 dim: 1,
359 mesh_degree: 4,
360 solution_degree: 4,
361 num_qpts: 6,
362 problem_size_requested: -1,
363 test: true,
364 quiet: true,
365 gallery: true,
366 };
367 assert!(example_1_vector(options).is_ok());
368 }
369
370 #[test]
example_1_vector_2d_gallery()371 fn example_1_vector_2d_gallery() {
372 let options = opt::Opt {
373 ceed_spec: "/cpu/self/ref/serial".to_string(),
374 dim: 2,
375 mesh_degree: 4,
376 solution_degree: 4,
377 num_qpts: 6,
378 problem_size_requested: -1,
379 test: true,
380 quiet: true,
381 gallery: true,
382 };
383 assert!(example_1_vector(options).is_ok());
384 }
385
386 #[test]
example_1_vector_3d_gallery()387 fn example_1_vector_3d_gallery() {
388 let options = opt::Opt {
389 ceed_spec: "/cpu/self/ref/serial".to_string(),
390 dim: 3,
391 mesh_degree: 4,
392 solution_degree: 4,
393 num_qpts: 6,
394 problem_size_requested: -1,
395 test: true,
396 quiet: true,
397 gallery: true,
398 };
399 assert!(example_1_vector(options).is_ok());
400 }
401 }
402
403 // ----------------------------------------------------------------------------
404