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