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