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 // ---------------------------------------------------------------------------- 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)] 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] 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] 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] 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] 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] 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] 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