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