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, VectorSliceWrapper}, 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_TOPOLOGY_LINE as isize, 119 Triangle = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_TRIANGLE as isize, 120 Quad = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_QUAD as isize, 121 Tet = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_TET as isize, 122 Pyramid = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_PYRAMID as isize, 123 Prism = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_PRISM as isize, 124 Hex = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_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 /// libCEED error messages - returning an Error without painc!ing indicates 157 /// the function call failed but the data is not corrupted 158 #[derive(Debug)] 159 pub struct Error { 160 pub message: String, 161 } 162 163 impl fmt::Display for Error { 164 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 165 write!(f, "{}", self.message) 166 } 167 } 168 169 // ----------------------------------------------------------------------------- 170 // Internal error checker 171 // ----------------------------------------------------------------------------- 172 #[doc(hidden)] 173 pub(crate) fn check_error(ceed_ptr: bind_ceed::Ceed, ierr: i32) -> Result<i32> { 174 // Return early if code is clean 175 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 176 return Ok(ierr); 177 } 178 // Retrieve error message 179 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 180 let c_str = unsafe { 181 bind_ceed::CeedGetErrorMessage(ceed_ptr, &mut ptr); 182 std::ffi::CStr::from_ptr(ptr) 183 }; 184 let message = c_str.to_string_lossy().to_string(); 185 // Panic if negative code, otherwise return error 186 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 187 panic!("{}", message); 188 } 189 Err(Error { message }) 190 } 191 192 // ----------------------------------------------------------------------------- 193 // Ceed error handler 194 // ----------------------------------------------------------------------------- 195 pub enum ErrorHandler { 196 ErrorAbort, 197 ErrorExit, 198 ErrorReturn, 199 ErrorStore, 200 } 201 202 // ----------------------------------------------------------------------------- 203 // Ceed context wrapper 204 // ----------------------------------------------------------------------------- 205 /// A Ceed is a library context representing control of a logical hardware 206 /// resource. 207 #[derive(Debug)] 208 pub struct Ceed { 209 ptr: bind_ceed::Ceed, 210 } 211 212 // ----------------------------------------------------------------------------- 213 // Destructor 214 // ----------------------------------------------------------------------------- 215 impl Drop for Ceed { 216 fn drop(&mut self) { 217 unsafe { 218 bind_ceed::CeedDestroy(&mut self.ptr); 219 } 220 } 221 } 222 223 // ----------------------------------------------------------------------------- 224 // Cloning 225 // ----------------------------------------------------------------------------- 226 impl Clone for Ceed { 227 /// Perform a shallow clone of a Ceed context 228 /// 229 /// ``` 230 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 231 /// let ceed_clone = ceed.clone(); 232 /// 233 /// println!(" original:{} \n clone:{}", ceed, ceed_clone); 234 /// ``` 235 fn clone(&self) -> Self { 236 let mut ptr_clone = std::ptr::null_mut(); 237 let ierr = unsafe { bind_ceed::CeedReferenceCopy(self.ptr, &mut ptr_clone) }; 238 self.check_error(ierr).expect("failed to clone Ceed"); 239 Self { ptr: ptr_clone } 240 } 241 } 242 243 // ----------------------------------------------------------------------------- 244 // Display 245 // ----------------------------------------------------------------------------- 246 impl fmt::Display for Ceed { 247 /// View a Ceed 248 /// 249 /// ``` 250 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 251 /// println!("{}", ceed); 252 /// ``` 253 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 254 let mut ptr = std::ptr::null_mut(); 255 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 256 let cstring = unsafe { 257 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 258 bind_ceed::CeedView(self.ptr, file); 259 bind_ceed::fclose(file); 260 CString::from_raw(ptr) 261 }; 262 cstring.to_string_lossy().fmt(f) 263 } 264 } 265 266 static REGISTER: Once = Once::new(); 267 268 // ----------------------------------------------------------------------------- 269 // Object constructors 270 // ----------------------------------------------------------------------------- 271 impl Ceed { 272 #[cfg_attr(feature = "katexit", katexit::katexit)] 273 /// Returns a Ceed context initialized with the specified resource 274 /// 275 /// # arguments 276 /// 277 /// * `resource` - Resource to use, e.g., "/cpu/self" 278 /// 279 /// ``` 280 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 281 /// ``` 282 pub fn init(resource: &str) -> Self { 283 Ceed::init_with_error_handler(resource, ErrorHandler::ErrorStore) 284 } 285 286 /// Returns a Ceed context initialized with the specified resource 287 /// 288 /// # arguments 289 /// 290 /// * `resource` - Resource to use, e.g., "/cpu/self" 291 /// 292 /// ``` 293 /// let ceed = libceed::Ceed::init_with_error_handler( 294 /// "/cpu/self/ref/serial", 295 /// libceed::ErrorHandler::ErrorAbort, 296 /// ); 297 /// ``` 298 pub fn init_with_error_handler(resource: &str, handler: ErrorHandler) -> Self { 299 REGISTER.call_once(|| unsafe { 300 bind_ceed::CeedRegisterAll(); 301 bind_ceed::CeedQFunctionRegisterAll(); 302 }); 303 304 // Convert to C string 305 let c_resource = CString::new(resource).expect("CString::new failed"); 306 307 // Get error handler pointer 308 let eh = match handler { 309 ErrorHandler::ErrorAbort => bind_ceed::CeedErrorAbort, 310 ErrorHandler::ErrorExit => bind_ceed::CeedErrorExit, 311 ErrorHandler::ErrorReturn => bind_ceed::CeedErrorReturn, 312 ErrorHandler::ErrorStore => bind_ceed::CeedErrorStore, 313 }; 314 315 // Call to libCEED 316 let mut ptr = std::ptr::null_mut(); 317 let mut ierr = unsafe { bind_ceed::CeedInit(c_resource.as_ptr() as *const i8, &mut ptr) }; 318 if ierr != 0 { 319 panic!("Error initializing backend resource: {}", resource) 320 } 321 ierr = unsafe { bind_ceed::CeedSetErrorHandler(ptr, Some(eh)) }; 322 let ceed = Ceed { ptr }; 323 ceed.check_error(ierr).unwrap(); 324 ceed 325 } 326 327 /// Default initializer for testing 328 #[doc(hidden)] 329 pub fn default_init() -> Self { 330 // Convert to C string 331 let resource = "/cpu/self/ref/serial"; 332 crate::Ceed::init(resource) 333 } 334 335 /// Internal error checker 336 #[doc(hidden)] 337 fn check_error(&self, ierr: i32) -> Result<i32> { 338 // Return early if code is clean 339 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 340 return Ok(ierr); 341 } 342 // Retrieve error message 343 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 344 let c_str = unsafe { 345 bind_ceed::CeedGetErrorMessage(self.ptr, &mut ptr); 346 std::ffi::CStr::from_ptr(ptr) 347 }; 348 let message = c_str.to_string_lossy().to_string(); 349 // Panic if negative code, otherwise return error 350 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 351 panic!("{}", message); 352 } 353 Err(Error { message }) 354 } 355 356 /// Returns full resource name for a Ceed context 357 /// 358 /// ``` 359 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 360 /// let resource = ceed.resource(); 361 /// 362 /// assert_eq!(resource, "/cpu/self/ref/serial".to_string()) 363 /// ``` 364 pub fn resource(&self) -> String { 365 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 366 let c_str = unsafe { 367 bind_ceed::CeedGetResource(self.ptr, &mut ptr); 368 std::ffi::CStr::from_ptr(ptr) 369 }; 370 c_str.to_string_lossy().to_string() 371 } 372 373 /// Returns a Vector of the specified length (does not allocate memory) 374 /// 375 /// # arguments 376 /// 377 /// * `n` - Length of vector 378 /// 379 /// ``` 380 /// # use libceed::prelude::*; 381 /// # fn main() -> libceed::Result<()> { 382 /// # let ceed = libceed::Ceed::default_init(); 383 /// let vec = ceed.vector(10)?; 384 /// # Ok(()) 385 /// # } 386 /// ``` 387 pub fn vector<'a>(&self, n: usize) -> Result<Vector<'a>> { 388 Vector::create(self, n) 389 } 390 391 /// Create a Vector initialized with the data (copied) from a slice 392 /// 393 /// # arguments 394 /// 395 /// * `slice` - Slice containing data 396 /// 397 /// ``` 398 /// # use libceed::prelude::*; 399 /// # fn main() -> libceed::Result<()> { 400 /// # let ceed = libceed::Ceed::default_init(); 401 /// let vec = ceed.vector_from_slice(&[1., 2., 3.])?; 402 /// assert_eq!(vec.length(), 3); 403 /// # Ok(()) 404 /// # } 405 /// ``` 406 pub fn vector_from_slice<'a>(&self, slice: &[crate::Scalar]) -> Result<Vector<'a>> { 407 Vector::from_slice(self, slice) 408 } 409 410 /// Returns an ElemRestriction, $\mathcal{E}$, which extracts the degrees of 411 /// freedom for each element from the local vector into the element vector 412 /// or assembles contributions from each element in the element vector to 413 /// the local vector 414 /// 415 /// # arguments 416 /// 417 /// * `nelem` - Number of elements described in the offsets array 418 /// * `elemsize` - Size (number of "nodes") per element 419 /// * `ncomp` - Number of field components per interpolation node (1 420 /// for scalar fields) 421 /// * `compstride` - Stride between components for the same Lvector "node". 422 /// Data for node `i`, component `j`, element `k` can be 423 /// found in the Lvector at index 424 /// `offsets[i + k*elemsize] + j*compstride`. 425 /// * `lsize` - The size of the Lvector. This vector may be larger 426 /// than the elements and fields given by this 427 /// restriction. 428 /// * `mtype` - Memory type of the offsets array, see CeedMemType 429 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 430 /// ordered list of the offsets (into the input CeedVector) 431 /// for the unknowns corresponding to element `i`, where 432 /// `0 <= i < nelem`. All offsets must be in the range 433 /// `[0, lsize - 1]`. 434 /// 435 /// ``` 436 /// # use libceed::prelude::*; 437 /// # fn main() -> libceed::Result<()> { 438 /// # let ceed = libceed::Ceed::default_init(); 439 /// let nelem = 3; 440 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 441 /// for i in 0..nelem { 442 /// ind[2 * i + 0] = i as i32; 443 /// ind[2 * i + 1] = (i + 1) as i32; 444 /// } 445 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 446 /// # Ok(()) 447 /// # } 448 /// ``` 449 pub fn elem_restriction<'a>( 450 &self, 451 nelem: usize, 452 elemsize: usize, 453 ncomp: usize, 454 compstride: usize, 455 lsize: usize, 456 mtype: MemType, 457 offsets: &[i32], 458 ) -> Result<ElemRestriction<'a>> { 459 ElemRestriction::create( 460 self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, 461 ) 462 } 463 464 /// Returns an ElemRestriction, $\mathcal{E}$, from an local vector to 465 /// an element vector where data can be indexed from the `strides` array 466 /// 467 /// # arguments 468 /// 469 /// * `nelem` - Number of elements described in the offsets array 470 /// * `elemsize` - Size (number of "nodes") per element 471 /// * `ncomp` - Number of field components per interpolation node (1 472 /// for scalar fields) 473 /// * `lsize` - The size of the Lvector. This vector may be larger 474 /// than the elements and fields given by this restriction. 475 /// * `strides` - Array for strides between `[nodes, components, elements]`. 476 /// Data for node `i`, component `j`, element `k` can be 477 /// found in the Lvector at index 478 /// `i*strides[0] + j*strides[1] + k*strides[2]`. 479 /// CEED_STRIDES_BACKEND may be used with vectors created 480 /// by a Ceed backend. 481 /// 482 /// ``` 483 /// # use libceed::prelude::*; 484 /// # fn main() -> libceed::Result<()> { 485 /// # let ceed = libceed::Ceed::default_init(); 486 /// let nelem = 3; 487 /// let strides: [i32; 3] = [1, 2, 2]; 488 /// let r = ceed.strided_elem_restriction(nelem, 2, 1, nelem * 2, strides)?; 489 /// # Ok(()) 490 /// # } 491 /// ``` 492 pub fn strided_elem_restriction<'a>( 493 &self, 494 nelem: usize, 495 elemsize: usize, 496 ncomp: usize, 497 lsize: usize, 498 strides: [i32; 3], 499 ) -> Result<ElemRestriction<'a>> { 500 ElemRestriction::create_strided(self, nelem, elemsize, ncomp, lsize, strides) 501 } 502 503 /// Returns an $H^1$ tensor-product Basis 504 /// 505 /// # arguments 506 /// 507 /// * `dim` - Topological dimension of element 508 /// * `ncomp` - Number of field components (1 for scalar fields) 509 /// * `P1d` - Number of Gauss-Lobatto nodes in one dimension. The 510 /// polynomial degree of the resulting `Q_k` element is 511 /// `k=P-1`. 512 /// * `Q1d` - Number of quadrature points in one dimension 513 /// * `interp1d` - Row-major `(Q1d * P1d)` matrix expressing the values of 514 /// nodal basis functions at quadrature points 515 /// * `grad1d` - Row-major `(Q1d * P1d)` matrix expressing derivatives of 516 /// nodal basis functions at quadrature points 517 /// * `qref1d` - Array of length `Q1d` holding the locations of quadrature 518 /// points on the 1D reference element `[-1, 1]` 519 /// * `qweight1d` - Array of length `Q1d` holding the quadrature weights on 520 /// the reference element 521 /// 522 /// ``` 523 /// # use libceed::prelude::*; 524 /// # fn main() -> libceed::Result<()> { 525 /// # let ceed = libceed::Ceed::default_init(); 526 /// let interp1d = [ 0.62994317, 0.47255875, -0.14950343, 0.04700152, 527 /// -0.07069480, 0.97297619, 0.13253993, -0.03482132, 528 /// -0.03482132, 0.13253993, 0.97297619, -0.07069480, 529 /// 0.04700152, -0.14950343, 0.47255875, 0.62994317]; 530 /// let grad1d = [-2.34183742, 2.78794489, -0.63510411, 0.18899664, 531 /// -0.51670214, -0.48795249, 1.33790510, -0.33325047, 532 // 0.33325047, -1.33790510, 0.48795249, 0.51670214, 533 /// -0.18899664, 0.63510411, -2.78794489, 2.34183742]; 534 /// let qref1d = [-0.86113631, -0.33998104, 0.33998104, 0.86113631]; 535 /// let qweight1d = [ 0.34785485, 0.65214515, 0.65214515, 0.34785485]; 536 /// let b = ceed. 537 /// basis_tensor_H1(2, 1, 4, 4, &interp1d, &grad1d, &qref1d, &qweight1d)?; 538 /// # Ok(()) 539 /// # } 540 /// ``` 541 pub fn basis_tensor_H1<'a>( 542 &self, 543 dim: usize, 544 ncomp: usize, 545 P1d: usize, 546 Q1d: usize, 547 interp1d: &[crate::Scalar], 548 grad1d: &[crate::Scalar], 549 qref1d: &[crate::Scalar], 550 qweight1d: &[crate::Scalar], 551 ) -> Result<Basis<'a>> { 552 Basis::create_tensor_H1( 553 self, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, 554 ) 555 } 556 557 /// Returns an $H^1$ Lagrange tensor-product Basis 558 /// 559 /// # arguments 560 /// 561 /// * `dim` - Topological dimension of element 562 /// * `ncomp` - Number of field components (1 for scalar fields) 563 /// * `P` - Number of Gauss-Lobatto nodes in one dimension. The 564 /// polynomial degree of the resulting `Q_k` element is `k=P-1`. 565 /// * `Q` - Number of quadrature points in one dimension 566 /// * `qmode` - Distribution of the `Q` quadrature points (affects order of 567 /// accuracy for the quadrature) 568 /// 569 /// ``` 570 /// # use libceed::prelude::*; 571 /// # fn main() -> libceed::Result<()> { 572 /// # let ceed = libceed::Ceed::default_init(); 573 /// let b = ceed.basis_tensor_H1_Lagrange(2, 1, 3, 4, QuadMode::Gauss)?; 574 /// # Ok(()) 575 /// # } 576 /// ``` 577 pub fn basis_tensor_H1_Lagrange<'a>( 578 &self, 579 dim: usize, 580 ncomp: usize, 581 P: usize, 582 Q: usize, 583 qmode: QuadMode, 584 ) -> Result<Basis<'a>> { 585 Basis::create_tensor_H1_Lagrange(self, dim, ncomp, P, Q, qmode) 586 } 587 588 /// Returns an $H-1$ Basis 589 /// 590 /// # arguments 591 /// 592 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 593 /// * `ncomp` - Number of field components (1 for scalar fields) 594 /// * `nnodes` - Total number of nodes 595 /// * `nqpts` - Total number of quadrature points 596 /// * `interp` - Row-major `(nqpts * nnodes)` matrix expressing the values of 597 /// nodal basis functions at quadrature points 598 /// * `grad` - Row-major `(nqpts * dim * nnodes)` matrix expressing 599 /// derivatives of nodal basis functions at quadrature points 600 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 601 /// points on the reference element `[-1, 1]` 602 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 603 /// the reference element 604 /// 605 /// ``` 606 /// # use libceed::prelude::*; 607 /// # fn main() -> libceed::Result<()> { 608 /// # let ceed = libceed::Ceed::default_init(); 609 /// let interp = [ 610 /// 0.12000000, 611 /// 0.48000000, 612 /// -0.12000000, 613 /// 0.48000000, 614 /// 0.16000000, 615 /// -0.12000000, 616 /// -0.12000000, 617 /// 0.48000000, 618 /// 0.12000000, 619 /// 0.16000000, 620 /// 0.48000000, 621 /// -0.12000000, 622 /// -0.11111111, 623 /// 0.44444444, 624 /// -0.11111111, 625 /// 0.44444444, 626 /// 0.44444444, 627 /// -0.11111111, 628 /// -0.12000000, 629 /// 0.16000000, 630 /// -0.12000000, 631 /// 0.48000000, 632 /// 0.48000000, 633 /// 0.12000000, 634 /// ]; 635 /// let grad = [ 636 /// -1.40000000, 637 /// 1.60000000, 638 /// -0.20000000, 639 /// -0.80000000, 640 /// 0.80000000, 641 /// 0.00000000, 642 /// 0.20000000, 643 /// -1.60000000, 644 /// 1.40000000, 645 /// -0.80000000, 646 /// 0.80000000, 647 /// 0.00000000, 648 /// -0.33333333, 649 /// 0.00000000, 650 /// 0.33333333, 651 /// -1.33333333, 652 /// 1.33333333, 653 /// 0.00000000, 654 /// 0.20000000, 655 /// 0.00000000, 656 /// -0.20000000, 657 /// -2.40000000, 658 /// 2.40000000, 659 /// 0.00000000, 660 /// -1.40000000, 661 /// -0.80000000, 662 /// 0.00000000, 663 /// 1.60000000, 664 /// 0.80000000, 665 /// -0.20000000, 666 /// 0.20000000, 667 /// -2.40000000, 668 /// 0.00000000, 669 /// 0.00000000, 670 /// 2.40000000, 671 /// -0.20000000, 672 /// -0.33333333, 673 /// -1.33333333, 674 /// 0.00000000, 675 /// 0.00000000, 676 /// 1.33333333, 677 /// 0.33333333, 678 /// 0.20000000, 679 /// -0.80000000, 680 /// 0.00000000, 681 /// -1.60000000, 682 /// 0.80000000, 683 /// 1.40000000, 684 /// ]; 685 /// let qref = [ 686 /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, 687 /// 0.60000000, 688 /// ]; 689 /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; 690 /// let b = ceed.basis_H1( 691 /// ElemTopology::Triangle, 692 /// 1, 693 /// 6, 694 /// 4, 695 /// &interp, 696 /// &grad, 697 /// &qref, 698 /// &qweight, 699 /// )?; 700 /// # Ok(()) 701 /// # } 702 /// ``` 703 pub fn basis_H1<'a>( 704 &self, 705 topo: ElemTopology, 706 ncomp: usize, 707 nnodes: usize, 708 nqpts: usize, 709 interp: &[crate::Scalar], 710 grad: &[crate::Scalar], 711 qref: &[crate::Scalar], 712 qweight: &[crate::Scalar], 713 ) -> Result<Basis<'a>> { 714 Basis::create_H1( 715 self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, 716 ) 717 } 718 719 /// Returns a QFunction for evaluating interior (volumetric) terms 720 /// of a weak form corresponding to the $L^2$ inner product 721 /// 722 /// $$ 723 /// \langle v, F(u) \rangle = \int_\Omega v \cdot f_0 \left( u, \nabla u \right) + \left( \nabla v \right) : f_1 \left( u, \nabla u \right), 724 /// $$ 725 /// 726 /// where $v \cdot f_0$ represents contraction over fields and $\nabla v : f_1$ 727 /// represents contraction over both fields and spatial dimensions. 728 /// 729 /// # arguments 730 /// 731 /// * `vlength` - Vector length. Caller must ensure that number of 732 /// quadrature points is a multiple of vlength. 733 /// * `f` - Boxed closure to evaluate weak form at quadrature points. 734 /// 735 /// ``` 736 /// # use libceed::prelude::*; 737 /// # fn main() -> libceed::Result<()> { 738 /// # let ceed = libceed::Ceed::default_init(); 739 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 740 /// // Iterate over quadrature points 741 /// v.iter_mut() 742 /// .zip(u.iter().zip(weights.iter())) 743 /// .for_each(|(v, (u, w))| *v = u * w); 744 /// 745 /// // Return clean error code 746 /// 0 747 /// }; 748 /// 749 /// let qf = ceed.q_function_interior(1, Box::new(user_f))?; 750 /// # Ok(()) 751 /// # } 752 /// ``` 753 pub fn q_function_interior<'a>( 754 &self, 755 vlength: usize, 756 f: Box<qfunction::QFunctionUserClosure>, 757 ) -> Result<QFunction<'a>> { 758 QFunction::create(self, vlength, f) 759 } 760 761 /// Returns a QFunction for evaluating interior (volumetric) terms 762 /// created by name 763 /// 764 /// # arguments 765 /// 766 /// * `name` - name of QFunction from libCEED gallery 767 /// 768 /// ``` 769 /// # use libceed::prelude::*; 770 /// # fn main() -> libceed::Result<()> { 771 /// # let ceed = libceed::Ceed::default_init(); 772 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 773 /// # Ok(()) 774 /// # } 775 /// ``` 776 pub fn q_function_interior_by_name<'a>(&self, name: &str) -> Result<QFunctionByName<'a>> { 777 QFunctionByName::create(self, name) 778 } 779 780 /// Returns an Operator and associate a QFunction. A Basis and 781 /// ElemRestriction can be associated with QFunction fields via 782 /// set_field(). 783 /// 784 /// # arguments 785 /// 786 /// * `qf` - QFunction defining the action of the operator at quadrature 787 /// points 788 /// * `dqf` - QFunction defining the action of the Jacobian of the qf (or 789 /// qfunction_none) 790 /// * `dqfT` - QFunction defining the action of the transpose of the 791 /// Jacobian of the qf (or qfunction_none) 792 /// 793 /// ``` 794 /// # use libceed::prelude::*; 795 /// # fn main() -> libceed::Result<()> { 796 /// # let ceed = libceed::Ceed::default_init(); 797 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 798 /// let op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 799 /// # Ok(()) 800 /// # } 801 /// ``` 802 pub fn operator<'a, 'b>( 803 &self, 804 qf: impl Into<QFunctionOpt<'b>>, 805 dqf: impl Into<QFunctionOpt<'b>>, 806 dqfT: impl Into<QFunctionOpt<'b>>, 807 ) -> Result<Operator<'a>> { 808 Operator::create(self, qf, dqf, dqfT) 809 } 810 811 /// Returns an Operator that composes the action of several Operators 812 /// 813 /// ``` 814 /// # use libceed::prelude::*; 815 /// # fn main() -> libceed::Result<()> { 816 /// # let ceed = libceed::Ceed::default_init(); 817 /// let op = ceed.composite_operator()?; 818 /// # Ok(()) 819 /// # } 820 /// ``` 821 pub fn composite_operator<'a>(&self) -> Result<CompositeOperator<'a>> { 822 CompositeOperator::create(self) 823 } 824 } 825 826 // ----------------------------------------------------------------------------- 827 // Tests 828 // ----------------------------------------------------------------------------- 829 #[cfg(test)] 830 mod tests { 831 use super::*; 832 833 fn ceed_t501() -> Result<()> { 834 let resource = "/cpu/self/ref/blocked"; 835 let ceed = Ceed::init(resource); 836 let nelem = 4; 837 let p = 3; 838 let q = 4; 839 let ndofs = p * nelem - nelem + 1; 840 841 // Vectors 842 let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 843 let mut qdata = ceed.vector(nelem * q)?; 844 qdata.set_value(0.0)?; 845 let mut u = ceed.vector(ndofs)?; 846 u.set_value(1.0)?; 847 let mut v = ceed.vector(ndofs)?; 848 v.set_value(0.0)?; 849 850 // Restrictions 851 let mut indx: Vec<i32> = vec![0; 2 * nelem]; 852 for i in 0..nelem { 853 indx[2 * i + 0] = i as i32; 854 indx[2 * i + 1] = (i + 1) as i32; 855 } 856 let rx = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &indx)?; 857 let mut indu: Vec<i32> = vec![0; p * nelem]; 858 for i in 0..nelem { 859 indu[p * i + 0] = i as i32; 860 indu[p * i + 1] = (i + 1) as i32; 861 indu[p * i + 2] = (i + 2) as i32; 862 } 863 let ru = ceed.elem_restriction(nelem, 3, 1, 1, ndofs, MemType::Host, &indu)?; 864 let strides: [i32; 3] = [1, q as i32, q as i32]; 865 let rq = ceed.strided_elem_restriction(nelem, q, 1, q * nelem, strides)?; 866 867 // Bases 868 let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 869 let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 870 871 // Build quadrature data 872 let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 873 ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 874 .field("dx", &rx, &bx, VectorOpt::Active)? 875 .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 876 .field("qdata", &rq, BasisOpt::Collocated, VectorOpt::Active)? 877 .apply(&x, &mut qdata)?; 878 879 // Mass operator 880 let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 881 let op_mass = ceed 882 .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 883 .field("u", &ru, &bu, VectorOpt::Active)? 884 .field("qdata", &rq, BasisOpt::Collocated, &qdata)? 885 .field("v", &ru, &bu, VectorOpt::Active)? 886 .check()?; 887 888 v.set_value(0.0)?; 889 op_mass.apply(&u, &mut v)?; 890 891 // Check 892 let sum: Scalar = v.view()?.iter().sum(); 893 assert!( 894 (sum - 2.0).abs() < 1e-15, 895 "Incorrect interval length computed" 896 ); 897 Ok(()) 898 } 899 900 #[test] 901 fn test_ceed_t501() { 902 assert!(ceed_t501().is_ok()); 903 } 904 } 905 906 // ----------------------------------------------------------------------------- 907