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