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 127 let (rstr_solution, _) = 128 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?; 129 let mesh_size = rstr_mesh.lvector_size(); 130 let solution_size = rstr_solution.lvector_size(); 131 if !quiet { 132 println!("Number of mesh nodes : {}", mesh_size / dim); 133 println!("Number of solution nodes : {}", solution_size); 134 } 135 136 // Create a Vector with the mesh coordinates 137 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?; 138 139 // Apply a transformation to the mesh coordinates 140 let exact_area = transform::transform_mesh_coordinates(dim, &mut mesh_coords)?; 141 142 // QFunction that builds the quadrature data for the diff operator 143 // -- QFunction from user closure 144 let build_diff = move |[jacobian, weights, ..]: QFunctionInputs, 145 [qdata, ..]: QFunctionOutputs| { 146 // Build quadrature data 147 match dim { 148 1 => qdata 149 .iter_mut() 150 .zip(jacobian.iter().zip(weights.iter())) 151 .for_each(|(qdata, (j, weight))| *qdata = weight / j), 152 2 => { 153 let q = qdata.len() / 3; 154 for i in 0..q { 155 let j11 = jacobian[i + q * 0]; 156 let j21 = jacobian[i + q * 1]; 157 let j12 = jacobian[i + q * 2]; 158 let j22 = jacobian[i + q * 3]; 159 let qw = weights[i] / (j11 * j22 - j21 * j12); 160 qdata[i + q * 0] = qw * (j12 * j12 + j22 * j22); 161 qdata[i + q * 1] = qw * (j11 * j11 + j21 * j21); 162 qdata[i + q * 2] = -qw * (j11 * j12 + j21 * j22); 163 } 164 } 165 3 => { 166 let q = qdata.len() / 6; 167 for i in 0..q { 168 let mut a = [0.0; 9]; 169 for j in 0..3 { 170 for k in 0..3 { 171 a[k * 3 + j] = jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 1) % 3))] 172 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 2) % 3))] 173 - jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 2) % 3))] 174 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 1) % 3))]; 175 } 176 } 177 let qw = weights[i] 178 / (jacobian[i + q * 0] * a[0 * 3 + 0] 179 + jacobian[i + q * 1] * a[0 * 3 + 1] 180 + jacobian[i + q * 2] * a[0 * 3 + 2]); 181 qdata[i + q * 0] = qw 182 * (a[0 * 3 + 0] * a[0 * 3 + 0] 183 + a[0 * 3 + 1] * a[0 * 3 + 1] 184 + a[0 * 3 + 2] * a[0 * 3 + 2]); 185 qdata[i + q * 1] = qw 186 * (a[1 * 3 + 0] * a[1 * 3 + 0] 187 + a[1 * 3 + 1] * a[1 * 3 + 1] 188 + a[1 * 3 + 2] * a[1 * 3 + 2]); 189 qdata[i + q * 2] = qw 190 * (a[2 * 3 + 0] * a[2 * 3 + 0] 191 + a[2 * 3 + 1] * a[2 * 3 + 1] 192 + a[2 * 3 + 2] * a[2 * 3 + 2]); 193 qdata[i + q * 3] = qw 194 * (a[1 * 3 + 0] * a[2 * 3 + 0] 195 + a[1 * 3 + 1] * a[2 * 3 + 1] 196 + a[1 * 3 + 2] * a[2 * 3 + 2]); 197 qdata[i + q * 4] = qw 198 * (a[0 * 3 + 0] * a[2 * 3 + 0] 199 + a[0 * 3 + 1] * a[2 * 3 + 1] 200 + a[0 * 3 + 2] * a[2 * 3 + 2]); 201 qdata[i + q * 5] = qw 202 * (a[0 * 3 + 0] * a[1 * 3 + 0] 203 + a[0 * 3 + 1] * a[1 * 3 + 1] 204 + a[0 * 3 + 2] * a[1 * 3 + 2]); 205 } 206 } 207 _ => unreachable!(), 208 }; 209 210 // Return clean error code 211 0 212 }; 213 let qf_build_closure = ceed 214 .q_function_interior(1, Box::new(build_diff))? 215 .input("dx", ncomp_x * dim, libceed::EvalMode::Grad)? 216 .input("weights", 1, libceed::EvalMode::Weight)? 217 .output("qdata", dim * (dim + 1) / 2, libceed::EvalMode::None)?; 218 // -- QFunction from gallery 219 let qf_build_named = { 220 let name = format!("Poisson{}DBuild", dim); 221 ceed.q_function_interior_by_name(&name)? 222 }; 223 // -- QFunction for use with Operator 224 let qf_build = if gallery { 225 QFunctionOpt::SomeQFunctionByName(&qf_build_named) 226 } else { 227 QFunctionOpt::SomeQFunction(&qf_build_closure) 228 }; 229 230 // Operator that build the quadrature data for the diff operator 231 let op_build = ceed 232 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)? 233 .name("build qdata")? 234 .field("dx", &rstr_mesh, &basis_mesh, VectorOpt::Active)? 235 .field( 236 "weights", 237 ElemRestrictionOpt::None, 238 &basis_mesh, 239 VectorOpt::None, 240 )? 241 .field("qdata", &rstr_qdata, BasisOpt::None, VectorOpt::Active)? 242 .check()?; 243 244 // Compute the quadrature data for the diff operator 245 let elem_qpts = num_qpts.pow(dim as u32); 246 let num_elem: usize = num_xyz.iter().take(dim).product(); 247 let mut qdata = ceed.vector(num_elem * elem_qpts * dim * (dim + 1) / 2)?; 248 op_build.apply(&mesh_coords, &mut qdata)?; 249 250 // QFunction that applies the diff operator 251 // -- QFunction from user closure 252 let apply_diff = move |[ug, qdata, ..]: QFunctionInputs, [vg, ..]: QFunctionOutputs| { 253 // Apply diffusion operator 254 match dim { 255 1 => vg 256 .iter_mut() 257 .zip(ug.iter().zip(qdata.iter())) 258 .for_each(|(vg, (ug, w))| *vg = ug * w), 259 2 => { 260 let q = qdata.len() / 3; 261 for i in 0..q { 262 let du = [ug[i + q * 0], ug[i + q * 1]]; 263 let dxdxdxdx_t = [ 264 [qdata[i + 0 * q], qdata[i + 2 * q]], 265 [qdata[i + 2 * q], qdata[i + 1 * q]], 266 ]; 267 for j in 0..2 { 268 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j]; 269 } 270 } 271 } 272 3 => { 273 let q = qdata.len() / 6; 274 for i in 0..q { 275 let du = [ug[i + q * 0], ug[i + q * 1], ug[i + q * 2]]; 276 let dxdxdxdx_t = [ 277 [qdata[i + 0 * q], qdata[i + 5 * q], qdata[i + 4 * q]], 278 [qdata[i + 5 * q], qdata[i + 1 * q], qdata[i + 3 * q]], 279 [qdata[i + 4 * q], qdata[i + 3 * q], qdata[i + 2 * q]], 280 ]; 281 for j in 0..3 { 282 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] 283 + du[1] * dxdxdxdx_t[1][j] 284 + du[2] * dxdxdxdx_t[2][j]; 285 } 286 } 287 } 288 _ => unreachable!(), 289 }; 290 291 // Return clean error code 292 0 293 }; 294 let qf_diff_closure = ceed 295 .q_function_interior(1, Box::new(apply_diff))? 296 .input("du", dim, libceed::EvalMode::Grad)? 297 .input("qdata", dim * (dim + 1) / 2, libceed::EvalMode::None)? 298 .output("dv", dim, libceed::EvalMode::Grad)?; 299 // -- QFunction from gallery 300 let qf_diff_named = { 301 let name = format!("Poisson{}DApply", dim); 302 ceed.q_function_interior_by_name(&name)? 303 }; 304 // -- QFunction for use with Operator 305 let qf_diff = if gallery { 306 QFunctionOpt::SomeQFunctionByName(&qf_diff_named) 307 } else { 308 QFunctionOpt::SomeQFunction(&qf_diff_closure) 309 }; 310 311 // Diff Operator 312 let op_diff = ceed 313 .operator(qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 314 .name("Poisson")? 315 .field("du", &rstr_solution, &basis_solution, VectorOpt::Active)? 316 .field("qdata", &rstr_qdata, BasisOpt::None, &qdata)? 317 .field("dv", &rstr_solution, &basis_solution, VectorOpt::Active)? 318 .check()?; 319 320 // Solution vectors 321 let mut u = ceed.vector(solution_size)?; 322 let mut v = ceed.vector(solution_size)?; 323 324 // Initialize u with sum of node coordinates 325 let coords = mesh_coords.view()?; 326 u.set_value(0.0)?; 327 for (i, u) in u.view_mut()?.iter_mut().enumerate() { 328 *u = (0..dim).map(|d| coords[i + d * solution_size]).sum(); 329 } 330 331 // Apply the diff operator 332 op_diff.apply(&u, &mut v)?; 333 334 // Compute the mesh surface area 335 let area: libceed::Scalar = v.view()?.iter().map(|v| (*v).abs()).sum(); 336 337 // Output results 338 if !quiet { 339 println!("Exact mesh surface area : {:.12}", exact_area); 340 println!("Computed mesh surface_area : {:.12}", area); 341 println!("Surface area error : {:.12e}", area - exact_area); 342 } 343 let tolerance = match dim { 344 1 => 10000.0 * libceed::EPSILON, 345 _ => 1E-1, 346 }; 347 let error = (area - exact_area).abs(); 348 if error > tolerance { 349 println!("Volume error too large: {:.12e}", error); 350 return Err(libceed::Error { 351 message: format!( 352 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 353 tolerance, error 354 ), 355 }); 356 } 357 Ok(()) 358 } 359 360 // ---------------------------------------------------------------------------- 361 // Tests 362 // ---------------------------------------------------------------------------- 363 #[cfg(test)] 364 mod tests { 365 use super::*; 366 367 #[test] 368 fn example_2_1d() { 369 let options = opt::Opt { 370 ceed_spec: "/cpu/self/ref/serial".to_string(), 371 dim: 1, 372 mesh_degree: 4, 373 solution_degree: 4, 374 num_qpts: 6, 375 problem_size_requested: -1, 376 test: true, 377 quiet: true, 378 gallery: false, 379 }; 380 assert!(example_2(options).is_ok()); 381 } 382 383 #[test] 384 fn example_2_2d() { 385 let options = opt::Opt { 386 ceed_spec: "/cpu/self/ref/serial".to_string(), 387 dim: 2, 388 mesh_degree: 4, 389 solution_degree: 4, 390 num_qpts: 6, 391 problem_size_requested: -1, 392 test: true, 393 quiet: true, 394 gallery: false, 395 }; 396 assert!(example_2(options).is_ok()); 397 } 398 399 #[test] 400 fn example_2_3d() { 401 let options = opt::Opt { 402 ceed_spec: "/cpu/self/ref/serial".to_string(), 403 dim: 3, 404 mesh_degree: 4, 405 solution_degree: 4, 406 num_qpts: 6, 407 problem_size_requested: -1, 408 test: true, 409 quiet: false, 410 gallery: false, 411 }; 412 assert!(example_2(options).is_ok()); 413 } 414 415 #[test] 416 fn example_2_1d_gallery() { 417 let options = opt::Opt { 418 ceed_spec: "/cpu/self/ref/serial".to_string(), 419 dim: 1, 420 mesh_degree: 4, 421 solution_degree: 4, 422 num_qpts: 6, 423 problem_size_requested: -1, 424 test: true, 425 quiet: true, 426 gallery: true, 427 }; 428 assert!(example_2(options).is_ok()); 429 } 430 431 #[test] 432 fn example_2_2d_gallery() { 433 let options = opt::Opt { 434 ceed_spec: "/cpu/self/ref/serial".to_string(), 435 dim: 2, 436 mesh_degree: 4, 437 solution_degree: 4, 438 num_qpts: 6, 439 problem_size_requested: -1, 440 test: true, 441 quiet: true, 442 gallery: true, 443 }; 444 assert!(example_2(options).is_ok()); 445 } 446 447 #[test] 448 fn example_2_3d_gallery() { 449 let options = opt::Opt { 450 ceed_spec: "/cpu/self/ref/serial".to_string(), 451 dim: 3, 452 mesh_degree: 4, 453 solution_degree: 4, 454 num_qpts: 6, 455 problem_size_requested: -1, 456 test: true, 457 quiet: true, 458 gallery: true, 459 }; 460 assert!(example_2(options).is_ok()); 461 } 462 } 463 464 // ---------------------------------------------------------------------------- 465