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