1 // Copyright (c) 2017-2025, 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 // ---------------------------------------------------------------------------- 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)] 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] 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] 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] 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] 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] 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] 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