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