1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative 16 17 //! # libCEED Rust Interface 18 //! 19 //! This is the documentation for the high level libCEED Rust interface. 20 //! See the full libCEED user manual [here](https://libceed.readthedocs.io). 21 //! 22 //! libCEED is a low-level API for for the efficient high-order discretization methods 23 //! developed by the ECP co-design Center for Efficient Exascale Discretizations (CEED). 24 //! While our focus is on high-order finite elements, the approach is mostly algebraic 25 //! and thus applicable to other discretizations in factored form. 26 //! 27 //! ## Usage 28 //! 29 //! To call libCEED from a Rust package, the following `Cargo.toml` can be used. 30 //! ```toml 31 //! [dependencies] 32 //! libceed = { git = "https://github.com/CEED/libCEED", branch = "main" } 33 //! ``` 34 //! 35 //! Supported features: 36 //! * `static` (default): link to static libceed.a 37 //! * `system`: use libceed from a system directory (otherwise, install from source) 38 //! 39 //! ``` 40 //! extern crate libceed; 41 //! 42 //! fn main() { 43 //! let ceed = libceed::Ceed::init("/cpu/self/ref"); 44 //! let xc = ceed.vector_from_slice(&[0., 0.5, 1.0]).unwrap(); 45 //! let xs = xc.view(); 46 //! assert_eq!(xs[..], [0., 0.5, 1.0]); 47 //! } 48 //! ``` 49 //! 50 //! ## Examples 51 //! 52 //! Examples of libCEED can be found in the libCEED GitHub repository under the 53 //! `examples/rust` folder. 54 //! 55 //! ## Development 56 //! 57 //! To develop libCEED, use `cargo build` in the `rust/` directory to install a 58 //! local copy and build the bindings. If you need custom flags for the C 59 //! project, we recommend using `make configure` to cache arguments. If you 60 //! disable the `static` feature, then you'll need to set `LD_LIBRARY_PATH` for 61 //! doctests to be able to find it. You can do this in `$CEED_DIR/lib` and set 62 //! `PKG_CONFIG_PATH`. 63 //! 64 //! Note: the `LD_LIBRARY_PATH` workarounds will become unnecessary if [this 65 //! issue](https://github.com/rust-lang/cargo/issues/1592) is resolved -- it's 66 //! currently closed, but the problem still exists. 67 68 // ----------------------------------------------------------------------------- 69 // Exceptions 70 // ----------------------------------------------------------------------------- 71 #![allow(non_snake_case)] 72 73 // ----------------------------------------------------------------------------- 74 // Crate prelude 75 // ----------------------------------------------------------------------------- 76 use crate::prelude::*; 77 use std::sync::Once; 78 79 pub mod prelude { 80 pub(crate) use libceed_sys::bind_ceed; 81 pub use crate::{ 82 basis::{self, Basis, BasisOpt}, 83 elem_restriction::{self, ElemRestriction, ElemRestrictionOpt}, 84 operator::{self, CompositeOperator, Operator}, 85 qfunction::{ 86 self, QFunction, QFunctionByName, QFunctionInputs, QFunctionOpt, QFunctionOutputs, 87 }, 88 vector::{self, Vector, VectorOpt}, 89 ElemTopology, EvalMode, MemType, NormType, QuadMode, TransposeMode, CEED_STRIDES_BACKEND, 90 MAX_QFUNCTION_FIELDS, 91 }; 92 pub(crate) use std::convert::TryFrom; 93 pub(crate) use std::ffi::CString; 94 pub(crate) use std::fmt; 95 } 96 97 // ----------------------------------------------------------------------------- 98 // Modules 99 // ----------------------------------------------------------------------------- 100 pub mod basis; 101 pub mod elem_restriction; 102 pub mod operator; 103 pub mod qfunction; 104 pub mod vector; 105 106 // ----------------------------------------------------------------------------- 107 // Constants for library 108 // ----------------------------------------------------------------------------- 109 const MAX_BUFFER_LENGTH: u64 = 4096; 110 pub const MAX_QFUNCTION_FIELDS: usize = 16; 111 pub const CEED_STRIDES_BACKEND: [i32; 3] = [0; 3]; 112 113 // ----------------------------------------------------------------------------- 114 // Enums for libCEED 115 // ----------------------------------------------------------------------------- 116 #[derive(Clone, Copy, PartialEq, Eq)] 117 /// Many Ceed interfaces take or return pointers to memory. This enum is used to 118 /// specify where the memory being provided or requested must reside. 119 pub enum MemType { 120 Host = bind_ceed::CeedMemType_CEED_MEM_HOST as isize, 121 Device = bind_ceed::CeedMemType_CEED_MEM_DEVICE as isize, 122 } 123 124 #[derive(Clone, Copy, PartialEq, Eq)] 125 // OwnPointer will not be used by user but is included for internal use 126 #[allow(dead_code)] 127 /// Conveys ownership status of arrays passed to Ceed interfaces. 128 pub(crate) enum CopyMode { 129 CopyValues = bind_ceed::CeedCopyMode_CEED_COPY_VALUES as isize, 130 UsePointer = bind_ceed::CeedCopyMode_CEED_USE_POINTER as isize, 131 OwnPointer = bind_ceed::CeedCopyMode_CEED_OWN_POINTER as isize, 132 } 133 134 #[derive(Clone, Copy, PartialEq, Eq)] 135 /// Denotes type of vector norm to be computed 136 pub enum NormType { 137 One = bind_ceed::CeedNormType_CEED_NORM_1 as isize, 138 Two = bind_ceed::CeedNormType_CEED_NORM_2 as isize, 139 Max = bind_ceed::CeedNormType_CEED_NORM_MAX as isize, 140 } 141 142 #[derive(Clone, Copy, PartialEq, Eq)] 143 /// Denotes whether a linear transformation or its transpose should be applied 144 pub enum TransposeMode { 145 NoTranspose = bind_ceed::CeedTransposeMode_CEED_NOTRANSPOSE as isize, 146 Transpose = bind_ceed::CeedTransposeMode_CEED_TRANSPOSE as isize, 147 } 148 149 #[derive(Clone, Copy, PartialEq, Eq)] 150 /// Type of quadrature; also used for location of nodes 151 pub enum QuadMode { 152 Gauss = bind_ceed::CeedQuadMode_CEED_GAUSS as isize, 153 GaussLobatto = bind_ceed::CeedQuadMode_CEED_GAUSS_LOBATTO as isize, 154 } 155 156 #[derive(Clone, Copy, PartialEq, Eq)] 157 /// Type of basis shape to create non-tensor H1 element basis 158 pub enum ElemTopology { 159 Line = bind_ceed::CeedElemTopology_CEED_LINE as isize, 160 Triangle = bind_ceed::CeedElemTopology_CEED_TRIANGLE as isize, 161 Quad = bind_ceed::CeedElemTopology_CEED_QUAD as isize, 162 Tet = bind_ceed::CeedElemTopology_CEED_TET as isize, 163 Pyramid = bind_ceed::CeedElemTopology_CEED_PYRAMID as isize, 164 Prism = bind_ceed::CeedElemTopology_CEED_PRISM as isize, 165 Hex = bind_ceed::CeedElemTopology_CEED_HEX as isize, 166 } 167 168 #[derive(Clone, Copy, PartialEq, Eq)] 169 /// Basis evaluation mode 170 pub enum EvalMode { 171 None = bind_ceed::CeedEvalMode_CEED_EVAL_NONE as isize, 172 Interp = bind_ceed::CeedEvalMode_CEED_EVAL_INTERP as isize, 173 Grad = bind_ceed::CeedEvalMode_CEED_EVAL_GRAD as isize, 174 Div = bind_ceed::CeedEvalMode_CEED_EVAL_DIV as isize, 175 Curl = bind_ceed::CeedEvalMode_CEED_EVAL_CURL as isize, 176 Weight = bind_ceed::CeedEvalMode_CEED_EVAL_WEIGHT as isize, 177 } 178 179 // ----------------------------------------------------------------------------- 180 // Ceed error 181 // ----------------------------------------------------------------------------- 182 type Result<T> = std::result::Result<T, CeedError>; 183 184 #[derive(Debug)] 185 pub struct CeedError { 186 pub message: String, 187 } 188 189 impl fmt::Display for CeedError { 190 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 191 write!(f, "{}", self.message) 192 } 193 } 194 195 // ----------------------------------------------------------------------------- 196 // Ceed error handler 197 // ----------------------------------------------------------------------------- 198 pub enum CeedErrorHandler { 199 ErrorAbort, 200 ErrorExit, 201 ErrorReturn, 202 ErrorStore, 203 } 204 205 // ----------------------------------------------------------------------------- 206 // Ceed context wrapper 207 // ----------------------------------------------------------------------------- 208 /// A Ceed is a library context representing control of a logical hardware 209 /// resource. 210 #[derive(Debug)] 211 pub struct Ceed { 212 ptr: bind_ceed::Ceed, 213 } 214 215 // ----------------------------------------------------------------------------- 216 // Destructor 217 // ----------------------------------------------------------------------------- 218 impl Drop for Ceed { 219 fn drop(&mut self) { 220 unsafe { 221 bind_ceed::CeedDestroy(&mut self.ptr); 222 } 223 } 224 } 225 226 // ----------------------------------------------------------------------------- 227 // Display 228 // ----------------------------------------------------------------------------- 229 impl fmt::Display for Ceed { 230 /// View a Ceed 231 /// 232 /// ``` 233 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 234 /// println!("{}", ceed); 235 /// ``` 236 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 237 let mut ptr = std::ptr::null_mut(); 238 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 239 let cstring = unsafe { 240 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 241 bind_ceed::CeedView(self.ptr, file); 242 bind_ceed::fclose(file); 243 CString::from_raw(ptr) 244 }; 245 cstring.to_string_lossy().fmt(f) 246 } 247 } 248 249 static REGISTER: Once = Once::new(); 250 251 // ----------------------------------------------------------------------------- 252 // Object constructors 253 // ----------------------------------------------------------------------------- 254 impl Ceed { 255 /// Returns a Ceed context initialized with the specified resource 256 /// 257 /// # arguments 258 /// 259 /// * `resource` - Resource to use, e.g., "/cpu/self" 260 /// 261 /// ``` 262 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 263 /// ``` 264 pub fn init(resource: &str) -> Self { 265 Ceed::init_with_error_handler(resource, CeedErrorHandler::ErrorStore) 266 } 267 268 /// Returns a Ceed context initialized with the specified resource 269 /// 270 /// # arguments 271 /// 272 /// * `resource` - Resource to use, e.g., "/cpu/self" 273 /// 274 /// ``` 275 /// let ceed = libceed::Ceed::init_with_error_handler( 276 /// "/cpu/self/ref/serial", 277 /// libceed::CeedErrorHandler::ErrorAbort, 278 /// ); 279 /// ``` 280 pub fn init_with_error_handler(resource: &str, handler: CeedErrorHandler) -> Self { 281 REGISTER.call_once(|| unsafe { 282 bind_ceed::CeedRegisterAll(); 283 bind_ceed::CeedQFunctionRegisterAll(); 284 }); 285 286 // Convert to C string 287 let c_resource = CString::new(resource).expect("CString::new failed"); 288 289 // Get error handler pointer 290 let eh = match handler { 291 CeedErrorHandler::ErrorAbort => bind_ceed::CeedErrorAbort, 292 CeedErrorHandler::ErrorExit => bind_ceed::CeedErrorExit, 293 CeedErrorHandler::ErrorReturn => bind_ceed::CeedErrorReturn, 294 CeedErrorHandler::ErrorStore => bind_ceed::CeedErrorStore, 295 }; 296 297 // Call to libCEED 298 let mut ptr = std::ptr::null_mut(); 299 let mut ierr = unsafe { bind_ceed::CeedInit(c_resource.as_ptr() as *const i8, &mut ptr) }; 300 if ierr != 0 { 301 panic!("Error initializing backend resource: {}", resource) 302 } 303 ierr = unsafe { bind_ceed::CeedSetErrorHandler(ptr, Some(eh)) }; 304 let ceed = Ceed { ptr }; 305 ceed.check_error(ierr).unwrap(); 306 ceed 307 } 308 309 /// Default initializer for testing 310 #[doc(hidden)] 311 pub fn default_init() -> Self { 312 // Convert to C string 313 let resource = "/cpu/self/ref/serial"; 314 crate::Ceed::init(resource) 315 } 316 317 /// Internal error checker 318 #[doc(hidden)] 319 fn check_error(&self, ierr: i32) -> Result<i32> { 320 // Return early if code is clean 321 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 322 return Ok(ierr); 323 } 324 // Retrieve error message 325 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 326 let c_str = unsafe { 327 bind_ceed::CeedGetErrorMessage(self.ptr, &mut ptr); 328 std::ffi::CStr::from_ptr(ptr) 329 }; 330 let message = c_str.to_string_lossy().to_string(); 331 // Panic if negative code, otherwise return error 332 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 333 panic!("{}", message); 334 } 335 Err(CeedError { message }) 336 } 337 338 /// Returns full resource name for a Ceed context 339 /// 340 /// ``` 341 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 342 /// let resource = ceed.resource(); 343 /// 344 /// assert_eq!(resource, "/cpu/self/ref/serial".to_string()) 345 /// ``` 346 pub fn resource(&self) -> String { 347 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 348 let c_str = unsafe { 349 bind_ceed::CeedGetResource(self.ptr, &mut ptr); 350 std::ffi::CStr::from_ptr(ptr) 351 }; 352 c_str.to_string_lossy().to_string() 353 } 354 355 /// Returns a CeedVector of the specified length (does not allocate memory) 356 /// 357 /// # arguments 358 /// 359 /// * `n` - Length of vector 360 /// 361 /// ``` 362 /// # use libceed::prelude::*; 363 /// # let ceed = libceed::Ceed::default_init(); 364 /// let vec = ceed.vector(10).unwrap(); 365 /// ``` 366 pub fn vector(&self, n: usize) -> Result<Vector> { 367 Vector::create(self, n) 368 } 369 370 /// Create a Vector initialized with the data (copied) from a slice 371 /// 372 /// # arguments 373 /// 374 /// * `slice` - Slice containing data 375 /// 376 /// ``` 377 /// # use libceed::prelude::*; 378 /// # let ceed = libceed::Ceed::default_init(); 379 /// let vec = ceed.vector_from_slice(&[1., 2., 3.]).unwrap(); 380 /// assert_eq!(vec.length(), 3); 381 /// ``` 382 pub fn vector_from_slice(&self, slice: &[f64]) -> Result<Vector> { 383 Vector::from_slice(self, slice) 384 } 385 386 /// Returns a ElemRestriction 387 /// 388 /// # arguments 389 /// 390 /// * `nelem` - Number of elements described in the offsets array 391 /// * `elemsize` - Size (number of "nodes") per element 392 /// * `ncomp` - Number of field components per interpolation node (1 393 /// for scalar fields) 394 /// * `compstride` - Stride between components for the same Lvector "node". 395 /// Data for node `i`, component `j`, element `k` can be 396 /// found in the Lvector at index 397 /// `offsets[i + k*elemsize] + j*compstride`. 398 /// * `lsize` - The size of the Lvector. This vector may be larger 399 /// than the elements and fields given by this 400 /// restriction. 401 /// * `mtype` - Memory type of the offsets array, see CeedMemType 402 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 403 /// ordered list of the offsets (into the input CeedVector) 404 /// for the unknowns corresponding to element `i`, where 405 /// `0 <= i < nelem`. All offsets must be in the range 406 /// `[0, lsize - 1]`. 407 /// 408 /// ``` 409 /// # use libceed::prelude::*; 410 /// # let ceed = libceed::Ceed::default_init(); 411 /// let nelem = 3; 412 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 413 /// for i in 0..nelem { 414 /// ind[2 * i + 0] = i as i32; 415 /// ind[2 * i + 1] = (i + 1) as i32; 416 /// } 417 /// let r = ceed 418 /// .elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind) 419 /// .unwrap(); 420 /// ``` 421 pub fn elem_restriction( 422 &self, 423 nelem: usize, 424 elemsize: usize, 425 ncomp: usize, 426 compstride: usize, 427 lsize: usize, 428 mtype: MemType, 429 offsets: &[i32], 430 ) -> Result<ElemRestriction> { 431 ElemRestriction::create( 432 self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, 433 ) 434 } 435 436 /// Returns a ElemRestriction 437 /// 438 /// # arguments 439 /// 440 /// * `nelem` - Number of elements described in the offsets array 441 /// * `elemsize` - Size (number of "nodes") per element 442 /// * `ncomp` - Number of field components per interpolation node (1 443 /// for scalar fields) 444 /// * `compstride` - Stride between components for the same Lvector "node". 445 /// Data for node `i`, component `j`, element `k` can be 446 /// found in the Lvector at index 447 /// `offsets[i + k*elemsize] + j*compstride`. 448 /// * `lsize` - The size of the Lvector. This vector may be larger 449 /// than the elements and fields given by this restriction. 450 /// * `strides` - Array for strides between `[nodes, components, elements]`. 451 /// Data for node `i`, component `j`, element `k` can be 452 /// found in the Lvector at index 453 /// `i*strides[0] + j*strides[1] + k*strides[2]`. 454 /// CEED_STRIDES_BACKEND may be used with vectors created 455 /// by a Ceed backend. 456 /// 457 /// ``` 458 /// # use libceed::prelude::*; 459 /// # let ceed = libceed::Ceed::default_init(); 460 /// let nelem = 3; 461 /// let strides: [i32; 3] = [1, 2, 2]; 462 /// let r = ceed 463 /// .strided_elem_restriction(nelem, 2, 1, nelem * 2, strides) 464 /// .unwrap(); 465 /// ``` 466 pub fn strided_elem_restriction( 467 &self, 468 nelem: usize, 469 elemsize: usize, 470 ncomp: usize, 471 lsize: usize, 472 strides: [i32; 3], 473 ) -> Result<ElemRestriction> { 474 ElemRestriction::create_strided(self, nelem, elemsize, ncomp, lsize, strides) 475 } 476 477 /// Returns a tensor-product basis 478 /// 479 /// # arguments 480 /// 481 /// * `dim` - Topological dimension of element 482 /// * `ncomp` - Number of field components (1 for scalar fields) 483 /// * `P1d` - Number of Gauss-Lobatto nodes in one dimension. The 484 /// polynomial degree of the resulting `Q_k` element is 485 /// `k=P-1`. 486 /// * `Q1d` - Number of quadrature points in one dimension 487 /// * `interp1d` - Row-major `(Q1d * P1d)` matrix expressing the values of 488 /// nodal basis functions at quadrature points 489 /// * `grad1d` - Row-major `(Q1d * P1d)` matrix expressing derivatives of 490 /// nodal basis functions at quadrature points 491 /// * `qref1d` - Array of length `Q1d` holding the locations of quadrature 492 /// points on the 1D reference element `[-1, 1]` 493 /// * `qweight1d` - Array of length `Q1d` holding the quadrature weights on 494 /// the reference element 495 /// 496 /// ``` 497 /// # use libceed::prelude::*; 498 /// # let ceed = libceed::Ceed::default_init(); 499 /// let interp1d = [ 0.62994317, 0.47255875, -0.14950343, 0.04700152, 500 /// -0.07069480, 0.97297619, 0.13253993, -0.03482132, 501 /// -0.03482132, 0.13253993, 0.97297619, -0.07069480, 502 /// 0.04700152, -0.14950343, 0.47255875, 0.62994317]; 503 /// let grad1d = [-2.34183742, 2.78794489, -0.63510411, 0.18899664, 504 /// -0.51670214, -0.48795249, 1.33790510, -0.33325047, 505 // 0.33325047, -1.33790510, 0.48795249, 0.51670214, 506 /// -0.18899664, 0.63510411, -2.78794489, 2.34183742]; 507 /// let qref1d = [-0.86113631, -0.33998104, 0.33998104, 0.86113631]; 508 /// let qweight1d = [ 0.34785485, 0.65214515, 0.65214515, 0.34785485]; 509 /// let b = ceed. 510 /// basis_tensor_H1(2, 1, 4, 4, &interp1d, &grad1d, &qref1d, &qweight1d).unwrap(); 511 /// ``` 512 pub fn basis_tensor_H1( 513 &self, 514 dim: usize, 515 ncomp: usize, 516 P1d: usize, 517 Q1d: usize, 518 interp1d: &[f64], 519 grad1d: &[f64], 520 qref1d: &[f64], 521 qweight1d: &[f64], 522 ) -> Result<Basis> { 523 Basis::create_tensor_H1( 524 self, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, 525 ) 526 } 527 528 /// Returns a tensor-product Lagrange basis 529 /// 530 /// # arguments 531 /// 532 /// * `dim` - Topological dimension of element 533 /// * `ncomp` - Number of field components (1 for scalar fields) 534 /// * `P` - Number of Gauss-Lobatto nodes in one dimension. The 535 /// polynomial degree of the resulting `Q_k` element is `k=P-1`. 536 /// * `Q` - Number of quadrature points in one dimension 537 /// * `qmode` - Distribution of the `Q` quadrature points (affects order of 538 /// accuracy for the quadrature) 539 /// 540 /// ``` 541 /// # use libceed::prelude::*; 542 /// # let ceed = libceed::Ceed::default_init(); 543 /// let b = ceed 544 /// .basis_tensor_H1_Lagrange(2, 1, 3, 4, QuadMode::Gauss) 545 /// .unwrap(); 546 /// ``` 547 pub fn basis_tensor_H1_Lagrange( 548 &self, 549 dim: usize, 550 ncomp: usize, 551 P: usize, 552 Q: usize, 553 qmode: QuadMode, 554 ) -> Result<Basis> { 555 Basis::create_tensor_H1_Lagrange(self, dim, ncomp, P, Q, qmode) 556 } 557 558 /// Returns a tensor-product basis 559 /// 560 /// # arguments 561 /// 562 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 563 /// * `ncomp` - Number of field components (1 for scalar fields) 564 /// * `nnodes` - Total number of nodes 565 /// * `nqpts` - Total number of quadrature points 566 /// * `interp` - Row-major `(nqpts * nnodes)` matrix expressing the values of 567 /// nodal basis functions at quadrature points 568 /// * `grad` - Row-major `(nqpts * dim * nnodes)` matrix expressing 569 /// derivatives of nodal basis functions at quadrature points 570 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 571 /// points on the reference element `[-1, 1]` 572 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 573 /// the reference element 574 /// 575 /// ``` 576 /// # use libceed::prelude::*; 577 /// # let ceed = libceed::Ceed::default_init(); 578 /// let interp = [ 579 /// 0.12000000, 580 /// 0.48000000, 581 /// -0.12000000, 582 /// 0.48000000, 583 /// 0.16000000, 584 /// -0.12000000, 585 /// -0.12000000, 586 /// 0.48000000, 587 /// 0.12000000, 588 /// 0.16000000, 589 /// 0.48000000, 590 /// -0.12000000, 591 /// -0.11111111, 592 /// 0.44444444, 593 /// -0.11111111, 594 /// 0.44444444, 595 /// 0.44444444, 596 /// -0.11111111, 597 /// -0.12000000, 598 /// 0.16000000, 599 /// -0.12000000, 600 /// 0.48000000, 601 /// 0.48000000, 602 /// 0.12000000, 603 /// ]; 604 /// let grad = [ 605 /// -1.40000000, 606 /// 1.60000000, 607 /// -0.20000000, 608 /// -0.80000000, 609 /// 0.80000000, 610 /// 0.00000000, 611 /// 0.20000000, 612 /// -1.60000000, 613 /// 1.40000000, 614 /// -0.80000000, 615 /// 0.80000000, 616 /// 0.00000000, 617 /// -0.33333333, 618 /// 0.00000000, 619 /// 0.33333333, 620 /// -1.33333333, 621 /// 1.33333333, 622 /// 0.00000000, 623 /// 0.20000000, 624 /// 0.00000000, 625 /// -0.20000000, 626 /// -2.40000000, 627 /// 2.40000000, 628 /// 0.00000000, 629 /// -1.40000000, 630 /// -0.80000000, 631 /// 0.00000000, 632 /// 1.60000000, 633 /// 0.80000000, 634 /// -0.20000000, 635 /// 0.20000000, 636 /// -2.40000000, 637 /// 0.00000000, 638 /// 0.00000000, 639 /// 2.40000000, 640 /// -0.20000000, 641 /// -0.33333333, 642 /// -1.33333333, 643 /// 0.00000000, 644 /// 0.00000000, 645 /// 1.33333333, 646 /// 0.33333333, 647 /// 0.20000000, 648 /// -0.80000000, 649 /// 0.00000000, 650 /// -1.60000000, 651 /// 0.80000000, 652 /// 1.40000000, 653 /// ]; 654 /// let qref = [ 655 /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, 656 /// 0.60000000, 657 /// ]; 658 /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; 659 /// let b = ceed 660 /// .basis_H1( 661 /// ElemTopology::Triangle, 662 /// 1, 663 /// 6, 664 /// 4, 665 /// &interp, 666 /// &grad, 667 /// &qref, 668 /// &qweight, 669 /// ) 670 /// .unwrap(); 671 /// ``` 672 pub fn basis_H1( 673 &self, 674 topo: ElemTopology, 675 ncomp: usize, 676 nnodes: usize, 677 nqpts: usize, 678 interp: &[f64], 679 grad: &[f64], 680 qref: &[f64], 681 qweight: &[f64], 682 ) -> Result<Basis> { 683 Basis::create_H1( 684 self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, 685 ) 686 } 687 688 /// Returns a CeedQFunction for evaluating interior (volumetric) terms 689 /// 690 /// # arguments 691 /// 692 /// * `vlength` - Vector length. Caller must ensure that number of 693 /// quadrature points is a multiple of vlength. 694 /// * `f` - Boxed closure to evaluate action at quadrature points. 695 /// 696 /// ``` 697 /// # use libceed::prelude::*; 698 /// # let ceed = libceed::Ceed::default_init(); 699 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 700 /// // Iterate over quadrature points 701 /// v.iter_mut() 702 /// .zip(u.iter().zip(weights.iter())) 703 /// .for_each(|(v, (u, w))| *v = u * w); 704 /// 705 /// // Return clean error code 706 /// 0 707 /// }; 708 /// 709 /// let qf = ceed.q_function_interior(1, Box::new(user_f)).unwrap(); 710 /// ``` 711 pub fn q_function_interior( 712 &self, 713 vlength: usize, 714 f: Box<qfunction::QFunctionUserClosure>, 715 ) -> Result<QFunction> { 716 QFunction::create(self, vlength, f) 717 } 718 719 /// Returns a CeedQFunction for evaluating interior (volumetric) terms 720 /// created by name 721 /// 722 /// ``` 723 /// # use libceed::prelude::*; 724 /// # let ceed = libceed::Ceed::default_init(); 725 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild").unwrap(); 726 /// ``` 727 pub fn q_function_interior_by_name(&self, name: &str) -> Result<QFunctionByName> { 728 QFunctionByName::create(self, name) 729 } 730 731 /// Returns a Operator and associate a QFunction. A Basis and 732 /// ElemRestriction can be associated with QFunction fields with 733 /// set_field(). 734 /// 735 /// * `qf` - QFunction defining the action of the operator at quadrature 736 /// points 737 /// * `dqf` - QFunction defining the action of the Jacobian of the qf (or 738 /// qfunction_none) 739 /// * `dqfT` - QFunction defining the action of the transpose of the 740 /// Jacobian of the qf (or qfunction_none) 741 /// 742 /// ``` 743 /// # use libceed::prelude::*; 744 /// # let ceed = libceed::Ceed::default_init(); 745 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild").unwrap(); 746 /// let op = ceed 747 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None) 748 /// .unwrap(); 749 /// ``` 750 pub fn operator<'b>( 751 &self, 752 qf: impl Into<QFunctionOpt<'b>>, 753 dqf: impl Into<QFunctionOpt<'b>>, 754 dqfT: impl Into<QFunctionOpt<'b>>, 755 ) -> Result<Operator> { 756 Operator::create(self, qf, dqf, dqfT) 757 } 758 759 /// Returns an Operator that composes the action of several Operators 760 /// 761 /// ``` 762 /// # use libceed::prelude::*; 763 /// # let ceed = libceed::Ceed::default_init(); 764 /// let op = ceed.composite_operator().unwrap(); 765 /// ``` 766 pub fn composite_operator(&self) -> Result<CompositeOperator> { 767 CompositeOperator::create(self) 768 } 769 } 770 771 // ----------------------------------------------------------------------------- 772 // Tests 773 // ----------------------------------------------------------------------------- 774 #[cfg(test)] 775 mod tests { 776 use super::*; 777 778 fn ceed_t501() -> Result<i32> { 779 let resource = "/cpu/self/ref/blocked"; 780 let ceed = Ceed::init(resource); 781 let nelem = 4; 782 let p = 3; 783 let q = 4; 784 let ndofs = p * nelem - nelem + 1; 785 786 // Vectors 787 let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 788 let mut qdata = ceed.vector(nelem * q)?; 789 qdata.set_value(0.0)?; 790 let mut u = ceed.vector(ndofs)?; 791 u.set_value(1.0)?; 792 let mut v = ceed.vector(ndofs)?; 793 v.set_value(0.0)?; 794 795 // Restrictions 796 let mut indx: Vec<i32> = vec![0; 2 * nelem]; 797 for i in 0..nelem { 798 indx[2 * i + 0] = i as i32; 799 indx[2 * i + 1] = (i + 1) as i32; 800 } 801 let rx = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &indx)?; 802 let mut indu: Vec<i32> = vec![0; p * nelem]; 803 for i in 0..nelem { 804 indu[p * i + 0] = i as i32; 805 indu[p * i + 1] = (i + 1) as i32; 806 indu[p * i + 2] = (i + 2) as i32; 807 } 808 let ru = ceed.elem_restriction(nelem, 3, 1, 1, ndofs, MemType::Host, &indu)?; 809 let strides: [i32; 3] = [1, q as i32, q as i32]; 810 let rq = ceed.strided_elem_restriction(nelem, q, 1, q * nelem, strides)?; 811 812 // Bases 813 let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 814 let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 815 816 // Build quadrature data 817 let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 818 ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 819 .field("dx", &rx, &bx, VectorOpt::Active)? 820 .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 821 .field("qdata", &rq, BasisOpt::Collocated, VectorOpt::Active)? 822 .apply(&x, &mut qdata)?; 823 824 // Mass operator 825 let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 826 let op_mass = ceed 827 .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 828 .field("u", &ru, &bu, VectorOpt::Active)? 829 .field("qdata", &rq, BasisOpt::Collocated, &qdata)? 830 .field("v", &ru, &bu, VectorOpt::Active)?; 831 832 v.set_value(0.0)?; 833 op_mass.apply(&u, &mut v)?; 834 835 // Check 836 let sum: f64 = v.view().iter().sum(); 837 assert!( 838 (sum - 2.0).abs() < 1e-15, 839 "Incorrect interval length computed" 840 ); 841 Ok(0) 842 } 843 844 #[test] 845 fn test_ceed_t501() { 846 assert!(ceed_t501().is_ok()); 847 } 848 } 849 850 // ----------------------------------------------------------------------------- 851