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