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