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 fn example_2(options: opt::Opt) -> libceed::Result<()> { 38 // Process command line arguments 39 let opt::Opt { 40 ceed_spec, 41 dim, 42 mesh_degree, 43 solution_degree, 44 num_qpts, 45 problem_size_requested, 46 test, 47 quiet, 48 gallery, 49 } = options; 50 assert!((0..=3).contains(&dim)); 51 assert!(mesh_degree >= 1); 52 assert!(solution_degree >= 1); 53 assert!(num_qpts >= 1); 54 let ncomp_x = dim; 55 let problem_size: i64 = if problem_size_requested < 0 { 56 if test { 57 16 * 16 * (dim * dim) as i64 58 } else { 59 256 * 1024 60 } 61 } else { 62 problem_size_requested 63 }; 64 65 // Summary output 66 if !quiet { 67 println!("Selected options: [command line option] : <current value>"); 68 println!(" Ceed specification [-c] : {}", ceed_spec); 69 println!(" Mesh dimension [-d] : {}", dim); 70 println!(" Mesh degree [-m] : {}", mesh_degree); 71 println!(" Solution degree [-p] : {}", solution_degree); 72 println!(" Num. 1D quadr. pts [-q] : {}", num_qpts); 73 println!(" Approx. # unknowns [-s] : {}", problem_size); 74 println!( 75 " QFunction source [-g] : {}", 76 if gallery { "gallery" } else { "user closure" } 77 ); 78 } 79 80 // Initalize ceed context 81 let ceed = Ceed::init(&ceed_spec); 82 83 // Mesh and solution bases 84 let basis_mesh = ceed.basis_tensor_H1_Lagrange( 85 dim, 86 ncomp_x, 87 mesh_degree + 1, 88 num_qpts, 89 libceed::QuadMode::Gauss, 90 )?; 91 let basis_solution = ceed.basis_tensor_H1_Lagrange( 92 dim, 93 1, 94 solution_degree + 1, 95 num_qpts, 96 libceed::QuadMode::Gauss, 97 )?; 98 99 // Determine mesh size from approximate problem size 100 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size); 101 if !quiet { 102 print!("\nMesh size : nx = {}", num_xyz[0]); 103 if dim > 1 { 104 print!(", ny = {}", num_xyz[1]); 105 } 106 if dim > 2 { 107 print!(", nz = {}", num_xyz[2]); 108 } 109 println!(); 110 } 111 112 // Build ElemRestriction objects describing the mesh and solution discrete 113 // representations 114 let (rstr_mesh, _) = 115 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?; 116 let (_, rstr_qdata) = mesh::build_cartesian_restriction( 117 &ceed, 118 dim, 119 num_xyz, 120 solution_degree, 121 dim * (dim + 1) / 2, 122 num_qpts, 123 )?; 124 125 let (rstr_solution, _) = 126 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?; 127 let mesh_size = rstr_mesh.lvector_size(); 128 let solution_size = rstr_solution.lvector_size(); 129 if !quiet { 130 println!("Number of mesh nodes : {}", mesh_size / dim); 131 println!("Number of solution nodes : {}", solution_size); 132 } 133 134 // Create a Vector with the mesh coordinates 135 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?; 136 137 // Apply a transformation to the mesh coordinates 138 let exact_area = transform::transform_mesh_coordinates(dim, &mut mesh_coords)?; 139 140 // QFunction that builds the quadrature data for the diff operator 141 // -- QFunction from user closure 142 let build_diff = move |[jacobian, weights, ..]: QFunctionInputs, 143 [qdata, ..]: QFunctionOutputs| { 144 // Build quadrature data 145 match dim { 146 1 => qdata 147 .iter_mut() 148 .zip(jacobian.iter().zip(weights.iter())) 149 .for_each(|(qdata, (j, weight))| *qdata = weight / j), 150 2 => { 151 let q = qdata.len() / 3; 152 for i in 0..q { 153 let j11 = jacobian[i + q * 0]; 154 let j21 = jacobian[i + q * 1]; 155 let j12 = jacobian[i + q * 2]; 156 let j22 = jacobian[i + q * 3]; 157 let qw = weights[i] / (j11 * j22 - j21 * j12); 158 qdata[i + q * 0] = qw * (j12 * j12 + j22 * j22); 159 qdata[i + q * 1] = qw * (j11 * j11 + j21 * j21); 160 qdata[i + q * 2] = -qw * (j11 * j12 + j21 * j22); 161 } 162 } 163 3 => { 164 let q = qdata.len() / 6; 165 for i in 0..q { 166 let mut a = [0.0; 9]; 167 for j in 0..3 { 168 for k in 0..3 { 169 a[k * 3 + j] = jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 1) % 3))] 170 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 2) % 3))] 171 - jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 2) % 3))] 172 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 1) % 3))]; 173 } 174 } 175 let qw = weights[i] 176 / (jacobian[i + q * 0] * a[0 * 3 + 0] 177 + jacobian[i + q * 1] * a[0 * 3 + 1] 178 + jacobian[i + q * 2] * a[0 * 3 + 2]); 179 qdata[i + q * 0] = qw 180 * (a[0 * 3 + 0] * a[0 * 3 + 0] 181 + a[0 * 3 + 1] * a[0 * 3 + 1] 182 + a[0 * 3 + 2] * a[0 * 3 + 2]); 183 qdata[i + q * 1] = qw 184 * (a[1 * 3 + 0] * a[1 * 3 + 0] 185 + a[1 * 3 + 1] * a[1 * 3 + 1] 186 + a[1 * 3 + 2] * a[1 * 3 + 2]); 187 qdata[i + q * 2] = qw 188 * (a[2 * 3 + 0] * a[2 * 3 + 0] 189 + a[2 * 3 + 1] * a[2 * 3 + 1] 190 + a[2 * 3 + 2] * a[2 * 3 + 2]); 191 qdata[i + q * 3] = qw 192 * (a[1 * 3 + 0] * a[2 * 3 + 0] 193 + a[1 * 3 + 1] * a[2 * 3 + 1] 194 + a[1 * 3 + 2] * a[2 * 3 + 2]); 195 qdata[i + q * 4] = qw 196 * (a[0 * 3 + 0] * a[2 * 3 + 0] 197 + a[0 * 3 + 1] * a[2 * 3 + 1] 198 + a[0 * 3 + 2] * a[2 * 3 + 2]); 199 qdata[i + q * 5] = qw 200 * (a[0 * 3 + 0] * a[1 * 3 + 0] 201 + a[0 * 3 + 1] * a[1 * 3 + 1] 202 + a[0 * 3 + 2] * a[1 * 3 + 2]); 203 } 204 } 205 _ => unreachable!(), 206 }; 207 208 // Return clean error code 209 0 210 }; 211 let qf_build_closure = ceed 212 .q_function_interior(1, Box::new(build_diff))? 213 .input("dx", ncomp_x * dim, libceed::EvalMode::Grad)? 214 .input("weights", 1, libceed::EvalMode::Weight)? 215 .output("qdata", dim * (dim + 1) / 2, libceed::EvalMode::None)?; 216 // -- QFunction from gallery 217 let qf_build_named = { 218 let name = format!("Poisson{}DBuild", dim); 219 ceed.q_function_interior_by_name(&name)? 220 }; 221 // -- QFunction for use with Operator 222 let qf_build = if gallery { 223 QFunctionOpt::SomeQFunctionByName(&qf_build_named) 224 } else { 225 QFunctionOpt::SomeQFunction(&qf_build_closure) 226 }; 227 228 // Operator that build the quadrature data for the diff operator 229 let op_build = ceed 230 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)? 231 .name("build qdata")? 232 .field("dx", &rstr_mesh, &basis_mesh, VectorOpt::Active)? 233 .field( 234 "weights", 235 ElemRestrictionOpt::None, 236 &basis_mesh, 237 VectorOpt::None, 238 )? 239 .field("qdata", &rstr_qdata, BasisOpt::None, VectorOpt::Active)? 240 .check()?; 241 242 // Compute the quadrature data for the diff operator 243 let elem_qpts = num_qpts.pow(dim as u32); 244 let num_elem: usize = num_xyz.iter().take(dim).product(); 245 let mut qdata = ceed.vector(num_elem * elem_qpts * dim * (dim + 1) / 2)?; 246 op_build.apply(&mesh_coords, &mut qdata)?; 247 248 // QFunction that applies the diff operator 249 // -- QFunction from user closure 250 let apply_diff = move |[ug, qdata, ..]: QFunctionInputs, [vg, ..]: QFunctionOutputs| { 251 // Apply diffusion operator 252 match dim { 253 1 => vg 254 .iter_mut() 255 .zip(ug.iter().zip(qdata.iter())) 256 .for_each(|(vg, (ug, w))| *vg = ug * w), 257 2 => { 258 let q = qdata.len() / 3; 259 for i in 0..q { 260 let du = [ug[i + q * 0], ug[i + q * 1]]; 261 let dxdxdxdx_t = [ 262 [qdata[i + 0 * q], qdata[i + 2 * q]], 263 [qdata[i + 2 * q], qdata[i + 1 * q]], 264 ]; 265 for j in 0..2 { 266 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j]; 267 } 268 } 269 } 270 3 => { 271 let q = qdata.len() / 6; 272 for i in 0..q { 273 let du = [ug[i + q * 0], ug[i + q * 1], ug[i + q * 2]]; 274 let dxdxdxdx_t = [ 275 [qdata[i + 0 * q], qdata[i + 5 * q], qdata[i + 4 * q]], 276 [qdata[i + 5 * q], qdata[i + 1 * q], qdata[i + 3 * q]], 277 [qdata[i + 4 * q], qdata[i + 3 * q], qdata[i + 2 * q]], 278 ]; 279 for j in 0..3 { 280 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] 281 + du[1] * dxdxdxdx_t[1][j] 282 + du[2] * dxdxdxdx_t[2][j]; 283 } 284 } 285 } 286 _ => unreachable!(), 287 }; 288 289 // Return clean error code 290 0 291 }; 292 let qf_diff_closure = ceed 293 .q_function_interior(1, Box::new(apply_diff))? 294 .input("du", dim, libceed::EvalMode::Grad)? 295 .input("qdata", dim * (dim + 1) / 2, libceed::EvalMode::None)? 296 .output("dv", dim, libceed::EvalMode::Grad)?; 297 // -- QFunction from gallery 298 let qf_diff_named = { 299 let name = format!("Poisson{}DApply", dim); 300 ceed.q_function_interior_by_name(&name)? 301 }; 302 // -- QFunction for use with Operator 303 let qf_diff = if gallery { 304 QFunctionOpt::SomeQFunctionByName(&qf_diff_named) 305 } else { 306 QFunctionOpt::SomeQFunction(&qf_diff_closure) 307 }; 308 309 // Diff Operator 310 let op_diff = ceed 311 .operator(qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 312 .name("Poisson")? 313 .field("du", &rstr_solution, &basis_solution, VectorOpt::Active)? 314 .field("qdata", &rstr_qdata, BasisOpt::None, &qdata)? 315 .field("dv", &rstr_solution, &basis_solution, VectorOpt::Active)? 316 .check()?; 317 318 // Solution vectors 319 let mut u = ceed.vector(solution_size)?; 320 let mut v = ceed.vector(solution_size)?; 321 322 // Initialize u with sum of node coordinates 323 let coords = mesh_coords.view()?; 324 u.set_value(0.0)?; 325 for (i, u) in u.view_mut()?.iter_mut().enumerate() { 326 *u = (0..dim).map(|d| coords[i + d * solution_size]).sum(); 327 } 328 329 // Apply the diff operator 330 op_diff.apply(&u, &mut v)?; 331 332 // Compute the mesh surface area 333 let area: libceed::Scalar = v.view()?.iter().map(|v| (*v).abs()).sum(); 334 335 // Output results 336 if !quiet { 337 println!("Exact mesh surface area : {:.12}", exact_area); 338 println!("Computed mesh surface_area : {:.12}", area); 339 println!("Surface area error : {:.12e}", area - exact_area); 340 } 341 let tolerance = match dim { 342 1 => 10000.0 * libceed::EPSILON, 343 _ => 1E-1, 344 }; 345 let error = (area - exact_area).abs(); 346 if error > tolerance { 347 println!("Volume error too large: {:.12e}", error); 348 return Err(libceed::Error { 349 message: format!( 350 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 351 tolerance, error 352 ), 353 }); 354 } 355 Ok(()) 356 } 357 358 // ---------------------------------------------------------------------------- 359 // Tests 360 // ---------------------------------------------------------------------------- 361 #[cfg(test)] 362 mod tests { 363 use super::*; 364 365 #[test] 366 fn example_2_1d() { 367 let options = opt::Opt { 368 ceed_spec: "/cpu/self/ref/serial".to_string(), 369 dim: 1, 370 mesh_degree: 4, 371 solution_degree: 4, 372 num_qpts: 6, 373 problem_size_requested: -1, 374 test: true, 375 quiet: true, 376 gallery: false, 377 }; 378 assert!(example_2(options).is_ok()); 379 } 380 381 #[test] 382 fn example_2_2d() { 383 let options = opt::Opt { 384 ceed_spec: "/cpu/self/ref/serial".to_string(), 385 dim: 2, 386 mesh_degree: 4, 387 solution_degree: 4, 388 num_qpts: 6, 389 problem_size_requested: -1, 390 test: true, 391 quiet: true, 392 gallery: false, 393 }; 394 assert!(example_2(options).is_ok()); 395 } 396 397 #[test] 398 fn example_2_3d() { 399 let options = opt::Opt { 400 ceed_spec: "/cpu/self/ref/serial".to_string(), 401 dim: 3, 402 mesh_degree: 4, 403 solution_degree: 4, 404 num_qpts: 6, 405 problem_size_requested: -1, 406 test: true, 407 quiet: false, 408 gallery: false, 409 }; 410 assert!(example_2(options).is_ok()); 411 } 412 413 #[test] 414 fn example_2_1d_gallery() { 415 let options = opt::Opt { 416 ceed_spec: "/cpu/self/ref/serial".to_string(), 417 dim: 1, 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: true, 425 }; 426 assert!(example_2(options).is_ok()); 427 } 428 429 #[test] 430 fn example_2_2d_gallery() { 431 let options = opt::Opt { 432 ceed_spec: "/cpu/self/ref/serial".to_string(), 433 dim: 2, 434 mesh_degree: 4, 435 solution_degree: 4, 436 num_qpts: 6, 437 problem_size_requested: -1, 438 test: true, 439 quiet: true, 440 gallery: true, 441 }; 442 assert!(example_2(options).is_ok()); 443 } 444 445 #[test] 446 fn example_2_3d_gallery() { 447 let options = opt::Opt { 448 ceed_spec: "/cpu/self/ref/serial".to_string(), 449 dim: 3, 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(options).is_ok()); 459 } 460 } 461 462 // ---------------------------------------------------------------------------- 463