1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 // Fenced `rust` code blocks included from README.md are executed as part of doctests. 9 #![doc = include_str!("../README.md")] 10 // ----------------------------------------------------------------------------- 11 // Exceptions 12 // ----------------------------------------------------------------------------- 13 #![allow(non_snake_case)] 14 15 // ----------------------------------------------------------------------------- 16 // Crate prelude 17 // ----------------------------------------------------------------------------- 18 use crate::prelude::*; 19 use std::sync::Once; 20 21 pub mod prelude { 22 pub use crate::{ 23 basis::{self, Basis, BasisOpt}, 24 elem_restriction::{self, ElemRestriction, ElemRestrictionOpt}, 25 operator::{self, CompositeOperator, Operator, OperatorField}, 26 qfunction::{ 27 self, QFunction, QFunctionByName, QFunctionField, QFunctionInputs, QFunctionOpt, 28 QFunctionOutputs, 29 }, 30 vector::{self, Vector, VectorOpt, VectorSliceWrapper}, 31 ElemTopology, EvalMode, MemType, NormType, QuadMode, Scalar, TransposeMode, 32 CEED_STRIDES_BACKEND, EPSILON, MAX_QFUNCTION_FIELDS, 33 }; 34 pub(crate) use libceed_sys::bind_ceed; 35 pub(crate) use std::convert::TryFrom; 36 pub(crate) use std::ffi::{CStr, CString}; 37 pub(crate) use std::fmt; 38 pub(crate) use std::marker::PhantomData; 39 } 40 41 // ----------------------------------------------------------------------------- 42 // Modules 43 // ----------------------------------------------------------------------------- 44 pub mod basis; 45 pub mod elem_restriction; 46 pub mod operator; 47 pub mod qfunction; 48 pub mod vector; 49 50 // ----------------------------------------------------------------------------- 51 // Typedef for scalar 52 // ----------------------------------------------------------------------------- 53 pub type Scalar = bind_ceed::CeedScalar; 54 55 // ----------------------------------------------------------------------------- 56 // Constants for library 57 // ----------------------------------------------------------------------------- 58 const MAX_BUFFER_LENGTH: usize = 4096; 59 pub const MAX_QFUNCTION_FIELDS: usize = 16; 60 pub const CEED_STRIDES_BACKEND: [i32; 3] = [0; 3]; 61 pub const EPSILON: crate::Scalar = bind_ceed::CEED_EPSILON as crate::Scalar; 62 63 // ----------------------------------------------------------------------------- 64 // Enums for libCEED 65 // ----------------------------------------------------------------------------- 66 #[derive(Clone, Copy, PartialEq, Eq)] 67 /// Many Ceed interfaces take or return pointers to memory. This enum is used to 68 /// specify where the memory being provided or requested must reside. 69 pub enum MemType { 70 Host = bind_ceed::CeedMemType_CEED_MEM_HOST as isize, 71 Device = bind_ceed::CeedMemType_CEED_MEM_DEVICE as isize, 72 } 73 74 #[derive(Clone, Copy, PartialEq, Eq)] 75 // OwnPointer will not be used by user but is included for internal use 76 #[allow(dead_code)] 77 /// Conveys ownership status of arrays passed to Ceed interfaces. 78 pub(crate) enum CopyMode { 79 CopyValues = bind_ceed::CeedCopyMode_CEED_COPY_VALUES as isize, 80 UsePointer = bind_ceed::CeedCopyMode_CEED_USE_POINTER as isize, 81 OwnPointer = bind_ceed::CeedCopyMode_CEED_OWN_POINTER as isize, 82 } 83 84 #[derive(Clone, Copy, PartialEq, Eq)] 85 /// Denotes type of vector norm to be computed 86 pub enum NormType { 87 One = bind_ceed::CeedNormType_CEED_NORM_1 as isize, 88 Two = bind_ceed::CeedNormType_CEED_NORM_2 as isize, 89 Max = bind_ceed::CeedNormType_CEED_NORM_MAX as isize, 90 } 91 92 #[derive(Clone, Copy, PartialEq, Eq)] 93 /// Denotes whether a linear transformation or its transpose should be applied 94 pub enum TransposeMode { 95 NoTranspose = bind_ceed::CeedTransposeMode_CEED_NOTRANSPOSE as isize, 96 Transpose = bind_ceed::CeedTransposeMode_CEED_TRANSPOSE as isize, 97 } 98 99 #[derive(Clone, Copy, PartialEq, Eq)] 100 /// Type of quadrature; also used for location of nodes 101 pub enum QuadMode { 102 Gauss = bind_ceed::CeedQuadMode_CEED_GAUSS as isize, 103 GaussLobatto = bind_ceed::CeedQuadMode_CEED_GAUSS_LOBATTO as isize, 104 } 105 106 #[derive(Clone, Copy, PartialEq, Eq)] 107 /// Type of basis shape to create non-tensor H1 element basis 108 pub enum ElemTopology { 109 Line = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_LINE as isize, 110 Triangle = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_TRIANGLE as isize, 111 Quad = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_QUAD as isize, 112 Tet = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_TET as isize, 113 Pyramid = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_PYRAMID as isize, 114 Prism = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_PRISM as isize, 115 Hex = bind_ceed::CeedElemTopology_CEED_TOPOLOGY_HEX as isize, 116 } 117 118 #[derive(Clone, Copy, Debug, PartialEq, Eq)] 119 /// Basis evaluation mode 120 pub enum EvalMode { 121 None = bind_ceed::CeedEvalMode_CEED_EVAL_NONE as isize, 122 Interp = bind_ceed::CeedEvalMode_CEED_EVAL_INTERP as isize, 123 Grad = bind_ceed::CeedEvalMode_CEED_EVAL_GRAD as isize, 124 Div = bind_ceed::CeedEvalMode_CEED_EVAL_DIV as isize, 125 Curl = bind_ceed::CeedEvalMode_CEED_EVAL_CURL as isize, 126 Weight = bind_ceed::CeedEvalMode_CEED_EVAL_WEIGHT as isize, 127 } 128 impl EvalMode { 129 pub(crate) fn from_u32(value: u32) -> EvalMode { 130 match value { 131 bind_ceed::CeedEvalMode_CEED_EVAL_NONE => EvalMode::None, 132 bind_ceed::CeedEvalMode_CEED_EVAL_INTERP => EvalMode::Interp, 133 bind_ceed::CeedEvalMode_CEED_EVAL_GRAD => EvalMode::Grad, 134 bind_ceed::CeedEvalMode_CEED_EVAL_DIV => EvalMode::Div, 135 bind_ceed::CeedEvalMode_CEED_EVAL_CURL => EvalMode::Curl, 136 bind_ceed::CeedEvalMode_CEED_EVAL_WEIGHT => EvalMode::Weight, 137 _ => panic!("Unknown value: {}", value), 138 } 139 } 140 } 141 142 // ----------------------------------------------------------------------------- 143 // Ceed error 144 // ----------------------------------------------------------------------------- 145 pub type Result<T> = std::result::Result<T, Error>; 146 147 /// libCEED error messages - returning an Error without painc!ing indicates 148 /// the function call failed but the data is not corrupted 149 #[derive(Debug)] 150 pub struct Error { 151 pub message: String, 152 } 153 154 impl fmt::Display for Error { 155 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 156 write!(f, "{}", self.message) 157 } 158 } 159 160 // ----------------------------------------------------------------------------- 161 // Internal error checker 162 // ----------------------------------------------------------------------------- 163 #[doc(hidden)] 164 pub(crate) fn check_error(ceed_ptr: bind_ceed::Ceed, ierr: i32) -> Result<i32> { 165 // Return early if code is clean 166 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 167 return Ok(ierr); 168 } 169 // Retrieve error message 170 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 171 let c_str = unsafe { 172 bind_ceed::CeedGetErrorMessage(ceed_ptr, &mut ptr); 173 std::ffi::CStr::from_ptr(ptr) 174 }; 175 let message = c_str.to_string_lossy().to_string(); 176 // Panic if negative code, otherwise return error 177 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 178 panic!("{}", message); 179 } 180 Err(Error { message }) 181 } 182 183 // ----------------------------------------------------------------------------- 184 // Ceed error handler 185 // ----------------------------------------------------------------------------- 186 pub enum ErrorHandler { 187 ErrorAbort, 188 ErrorExit, 189 ErrorReturn, 190 ErrorStore, 191 } 192 193 // ----------------------------------------------------------------------------- 194 // Ceed context wrapper 195 // ----------------------------------------------------------------------------- 196 /// A Ceed is a library context representing control of a logical hardware 197 /// resource. 198 #[derive(Debug)] 199 pub struct Ceed { 200 ptr: bind_ceed::Ceed, 201 } 202 203 // ----------------------------------------------------------------------------- 204 // Destructor 205 // ----------------------------------------------------------------------------- 206 impl Drop for Ceed { 207 fn drop(&mut self) { 208 unsafe { 209 bind_ceed::CeedDestroy(&mut self.ptr); 210 } 211 } 212 } 213 214 // ----------------------------------------------------------------------------- 215 // Cloning 216 // ----------------------------------------------------------------------------- 217 impl Clone for Ceed { 218 /// Perform a shallow clone of a Ceed context 219 /// 220 /// ``` 221 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 222 /// let ceed_clone = ceed.clone(); 223 /// 224 /// println!(" original:{} \n clone:{}", ceed, ceed_clone); 225 /// ``` 226 fn clone(&self) -> Self { 227 let mut ptr_clone = std::ptr::null_mut(); 228 let ierr = unsafe { bind_ceed::CeedReferenceCopy(self.ptr, &mut ptr_clone) }; 229 self.check_error(ierr).expect("failed to clone Ceed"); 230 Self { ptr: ptr_clone } 231 } 232 } 233 234 // ----------------------------------------------------------------------------- 235 // Display 236 // ----------------------------------------------------------------------------- 237 impl fmt::Display for Ceed { 238 /// View a Ceed 239 /// 240 /// ``` 241 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 242 /// println!("{}", ceed); 243 /// ``` 244 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 245 let mut ptr = std::ptr::null_mut(); 246 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 247 let cstring = unsafe { 248 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 249 bind_ceed::CeedView(self.ptr, file); 250 bind_ceed::fclose(file); 251 CString::from_raw(ptr) 252 }; 253 cstring.to_string_lossy().fmt(f) 254 } 255 } 256 257 static REGISTER: Once = Once::new(); 258 259 // ----------------------------------------------------------------------------- 260 // Object constructors 261 // ----------------------------------------------------------------------------- 262 impl Ceed { 263 #[cfg_attr(feature = "katexit", katexit::katexit)] 264 /// Returns a Ceed context initialized with the specified resource 265 /// 266 /// # arguments 267 /// 268 /// * `resource` - Resource to use, e.g., "/cpu/self" 269 /// 270 /// ``` 271 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 272 /// ``` 273 pub fn init(resource: &str) -> Self { 274 Ceed::init_with_error_handler(resource, ErrorHandler::ErrorStore) 275 } 276 277 /// Returns a Ceed context initialized with the specified resource 278 /// 279 /// # arguments 280 /// 281 /// * `resource` - Resource to use, e.g., "/cpu/self" 282 /// 283 /// ``` 284 /// let ceed = libceed::Ceed::init_with_error_handler( 285 /// "/cpu/self/ref/serial", 286 /// libceed::ErrorHandler::ErrorAbort, 287 /// ); 288 /// ``` 289 pub fn init_with_error_handler(resource: &str, handler: ErrorHandler) -> Self { 290 REGISTER.call_once(|| unsafe { 291 bind_ceed::CeedRegisterAll(); 292 bind_ceed::CeedQFunctionRegisterAll(); 293 }); 294 295 // Convert to C string 296 let c_resource = CString::new(resource).expect("CString::new failed"); 297 298 // Get error handler pointer 299 let eh = match handler { 300 ErrorHandler::ErrorAbort => bind_ceed::CeedErrorAbort, 301 ErrorHandler::ErrorExit => bind_ceed::CeedErrorExit, 302 ErrorHandler::ErrorReturn => bind_ceed::CeedErrorReturn, 303 ErrorHandler::ErrorStore => bind_ceed::CeedErrorStore, 304 }; 305 306 // Call to libCEED 307 let mut ptr = std::ptr::null_mut(); 308 let mut ierr = unsafe { bind_ceed::CeedInit(c_resource.as_ptr() as *const i8, &mut ptr) }; 309 if ierr != 0 { 310 panic!("Error initializing backend resource: {}", resource) 311 } 312 ierr = unsafe { bind_ceed::CeedSetErrorHandler(ptr, Some(eh)) }; 313 let ceed = Ceed { ptr }; 314 ceed.check_error(ierr).unwrap(); 315 ceed 316 } 317 318 /// Default initializer for testing 319 #[doc(hidden)] 320 pub fn default_init() -> Self { 321 // Convert to C string 322 let resource = "/cpu/self/ref/serial"; 323 crate::Ceed::init(resource) 324 } 325 326 /// Internal error checker 327 #[doc(hidden)] 328 fn check_error(&self, ierr: i32) -> Result<i32> { 329 // Return early if code is clean 330 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 331 return Ok(ierr); 332 } 333 // Retrieve error message 334 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 335 let c_str = unsafe { 336 bind_ceed::CeedGetErrorMessage(self.ptr, &mut ptr); 337 std::ffi::CStr::from_ptr(ptr) 338 }; 339 let message = c_str.to_string_lossy().to_string(); 340 // Panic if negative code, otherwise return error 341 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 342 panic!("{}", message); 343 } 344 Err(Error { message }) 345 } 346 347 /// Returns full resource name for a Ceed context 348 /// 349 /// ``` 350 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 351 /// let resource = ceed.resource(); 352 /// 353 /// assert_eq!(resource, "/cpu/self/ref/serial".to_string()) 354 /// ``` 355 pub fn resource(&self) -> String { 356 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 357 let c_str = unsafe { 358 bind_ceed::CeedGetResource(self.ptr, &mut ptr); 359 std::ffi::CStr::from_ptr(ptr) 360 }; 361 c_str.to_string_lossy().to_string() 362 } 363 364 /// Returns a Vector of the specified length (does not allocate memory) 365 /// 366 /// # arguments 367 /// 368 /// * `n` - Length of vector 369 /// 370 /// ``` 371 /// # use libceed::prelude::*; 372 /// # fn main() -> libceed::Result<()> { 373 /// # let ceed = libceed::Ceed::default_init(); 374 /// let vec = ceed.vector(10)?; 375 /// # Ok(()) 376 /// # } 377 /// ``` 378 pub fn vector<'a>(&self, n: usize) -> Result<Vector<'a>> { 379 Vector::create(self, n) 380 } 381 382 /// Create a Vector initialized with the data (copied) from a slice 383 /// 384 /// # arguments 385 /// 386 /// * `slice` - Slice containing data 387 /// 388 /// ``` 389 /// # use libceed::prelude::*; 390 /// # fn main() -> libceed::Result<()> { 391 /// # let ceed = libceed::Ceed::default_init(); 392 /// let vec = ceed.vector_from_slice(&[1., 2., 3.])?; 393 /// assert_eq!(vec.length(), 3); 394 /// # Ok(()) 395 /// # } 396 /// ``` 397 pub fn vector_from_slice<'a>(&self, slice: &[crate::Scalar]) -> Result<Vector<'a>> { 398 Vector::from_slice(self, slice) 399 } 400 401 /// Returns an ElemRestriction, $\mathcal{E}$, which extracts the degrees of 402 /// freedom for each element from the local vector into the element vector 403 /// or assembles contributions from each element in the element vector to 404 /// the local vector 405 /// 406 /// # arguments 407 /// 408 /// * `nelem` - Number of elements described in the offsets array 409 /// * `elemsize` - Size (number of "nodes") per element 410 /// * `ncomp` - Number of field components per interpolation node (1 411 /// for scalar fields) 412 /// * `compstride` - Stride between components for the same Lvector "node". 413 /// Data for node `i`, component `j`, element `k` can be 414 /// found in the Lvector at index 415 /// `offsets[i + k*elemsize] + j*compstride`. 416 /// * `lsize` - The size of the Lvector. This vector may be larger 417 /// than the elements and fields given by this 418 /// restriction. 419 /// * `mtype` - Memory type of the offsets array, see CeedMemType 420 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 421 /// ordered list of the offsets (into the input CeedVector) 422 /// for the unknowns corresponding to element `i`, where 423 /// `0 <= i < nelem`. All offsets must be in the range 424 /// `[0, lsize - 1]`. 425 /// 426 /// ``` 427 /// # use libceed::prelude::*; 428 /// # fn main() -> libceed::Result<()> { 429 /// # let ceed = libceed::Ceed::default_init(); 430 /// let nelem = 3; 431 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 432 /// for i in 0..nelem { 433 /// ind[2 * i + 0] = i as i32; 434 /// ind[2 * i + 1] = (i + 1) as i32; 435 /// } 436 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 437 /// # Ok(()) 438 /// # } 439 /// ``` 440 pub fn elem_restriction<'a>( 441 &self, 442 nelem: usize, 443 elemsize: usize, 444 ncomp: usize, 445 compstride: usize, 446 lsize: usize, 447 mtype: MemType, 448 offsets: &[i32], 449 ) -> Result<ElemRestriction<'a>> { 450 ElemRestriction::create( 451 self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, 452 ) 453 } 454 455 /// Returns an oriented ElemRestriction, $\mathcal{E}$, which extracts the 456 /// degrees of freedom for each element from the local vector into the 457 /// element vector or assembles contributions from each element in the 458 /// element vector to the local vector 459 /// 460 /// # arguments 461 /// 462 /// * `nelem` - Number of elements described in the offsets array 463 /// * `elemsize` - Size (number of "nodes") per element 464 /// * `ncomp` - Number of field components per interpolation node (1 465 /// for scalar fields) 466 /// * `compstride` - Stride between components for the same Lvector "node". 467 /// Data for node `i`, component `j`, element `k` can be 468 /// found in the Lvector at index 469 /// `offsets[i + k*elemsize] + j*compstride`. 470 /// * `lsize` - The size of the Lvector. This vector may be larger 471 /// than the elements and fields given by this 472 /// restriction. 473 /// * `mtype` - Memory type of the offsets array, see CeedMemType 474 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 475 /// ordered list of the offsets (into the input CeedVector) 476 /// for the unknowns corresponding to element `i`, where 477 /// `0 <= i < nelem`. All offsets must be in the range 478 /// `[0, lsize - 1]`. 479 /// * `orients` - Array of shape `[nelem, elemsize]`. Row `i` holds the 480 /// ordered list of the orientations for the unknowns 481 /// corresponding to element `i`, with bool `false` used 482 /// for positively oriented and `true` to flip the 483 /// orientation. 484 /// 485 /// ``` 486 /// # use libceed::prelude::*; 487 /// # fn main() -> libceed::Result<()> { 488 /// # let ceed = libceed::Ceed::default_init(); 489 /// let nelem = 3; 490 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 491 /// let mut orients: Vec<bool> = vec![false; 2 * nelem]; 492 /// for i in 0..nelem { 493 /// ind[2 * i + 0] = i as i32; 494 /// ind[2 * i + 1] = (i + 1) as i32; 495 /// orients[2 * i + 0] = (i % 2) > 0; // flip the dofs on element 1, 3, ... 496 /// orients[2 * i + 1] = (i % 2) > 0; 497 /// } 498 /// let r = 499 /// ceed.oriented_elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind, &orients)?; 500 /// # Ok(()) 501 /// # } 502 /// ``` 503 pub fn oriented_elem_restriction<'a>( 504 &self, 505 nelem: usize, 506 elemsize: usize, 507 ncomp: usize, 508 compstride: usize, 509 lsize: usize, 510 mtype: MemType, 511 offsets: &[i32], 512 orients: &[bool], 513 ) -> Result<ElemRestriction<'a>> { 514 ElemRestriction::create_oriented( 515 self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, orients, 516 ) 517 } 518 519 /// Returns a curl-oriented ElemRestriction, $\mathcal{E}$, which extracts 520 /// the degrees of freedom for each element from the local vector into 521 /// the element vector or assembles contributions from each element in 522 /// the element vector to the local vector 523 /// 524 /// # arguments 525 /// 526 /// * `nelem` - Number of elements described in the offsets array 527 /// * `elemsize` - Size (number of "nodes") per element 528 /// * `ncomp` - Number of field components per interpolation node (1 529 /// for scalar fields) 530 /// * `compstride` - Stride between components for the same Lvector "node". 531 /// Data for node `i`, component `j`, element `k` can be 532 /// found in the Lvector at index 533 /// `offsets[i + k*elemsize] + j*compstride`. 534 /// * `lsize` - The size of the Lvector. This vector may be larger 535 /// than the elements and fields given by this 536 /// restriction. 537 /// * `mtype` - Memory type of the offsets array, see CeedMemType 538 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 539 /// ordered list of the offsets (into the input CeedVector) 540 /// for the unknowns corresponding to element `i`, where 541 /// `0 <= i < nelem`. All offsets must be in the range 542 /// `[0, lsize - 1]`. 543 /// * `curlorients` - Array of shape `[nelem, 3 * elemsize]`. Row `i` holds 544 /// a row-major tridiagonal matrix (`curlorients[i, 0] = 545 /// curlorients[i, 3 * elemsize - 1] = 0`, where 546 /// `0 <= i < nelem`) which is applied to the element 547 /// unknowns upon restriction. 548 /// 549 /// ``` 550 /// # use libceed::prelude::*; 551 /// # fn main() -> libceed::Result<()> { 552 /// # let ceed = libceed::Ceed::default_init(); 553 /// let nelem = 3; 554 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 555 /// let mut curlorients: Vec<i8> = vec![0; 3 * 2 * nelem]; 556 /// for i in 0..nelem { 557 /// ind[2 * i + 0] = i as i32; 558 /// ind[2 * i + 1] = (i + 1) as i32; 559 /// curlorients[3 * 2 * i] = 0 as i8; 560 /// curlorients[3 * 2 * (i + 1) - 1] = 0 as i8; 561 /// if (i % 2 > 0) { 562 /// // T = [0 -1] 563 /// // [-1 0] 564 /// curlorients[3 * 2 * i + 1] = 0 as i8; 565 /// curlorients[3 * 2 * i + 2] = -1 as i8; 566 /// curlorients[3 * 2 * i + 3] = -1 as i8; 567 /// curlorients[3 * 2 * i + 4] = 0 as i8; 568 /// } else { 569 /// // T = I 570 /// curlorients[3 * 2 * i + 1] = 1 as i8; 571 /// curlorients[3 * 2 * i + 2] = 0 as i8; 572 /// curlorients[3 * 2 * i + 3] = 0 as i8; 573 /// curlorients[3 * 2 * i + 4] = 1 as i8; 574 /// } 575 /// } 576 /// let r = ceed.curl_oriented_elem_restriction( 577 /// nelem, 578 /// 2, 579 /// 1, 580 /// 1, 581 /// nelem + 1, 582 /// MemType::Host, 583 /// &ind, 584 /// &curlorients, 585 /// )?; 586 /// # Ok(()) 587 /// # } 588 /// ``` 589 pub fn curl_oriented_elem_restriction<'a>( 590 &self, 591 nelem: usize, 592 elemsize: usize, 593 ncomp: usize, 594 compstride: usize, 595 lsize: usize, 596 mtype: MemType, 597 offsets: &[i32], 598 curlorients: &[i8], 599 ) -> Result<ElemRestriction<'a>> { 600 ElemRestriction::create_curl_oriented( 601 self, 602 nelem, 603 elemsize, 604 ncomp, 605 compstride, 606 lsize, 607 mtype, 608 offsets, 609 curlorients, 610 ) 611 } 612 613 /// Returns an ElemRestriction, $\mathcal{E}$, from an local vector to 614 /// an element vector where data can be indexed from the `strides` array 615 /// 616 /// # arguments 617 /// 618 /// * `nelem` - Number of elements described in the offsets array 619 /// * `elemsize` - Size (number of "nodes") per element 620 /// * `ncomp` - Number of field components per interpolation node (1 621 /// for scalar fields) 622 /// * `lsize` - The size of the Lvector. This vector may be larger 623 /// than the elements and fields given by this restriction. 624 /// * `strides` - Array for strides between `[nodes, components, elements]`. 625 /// Data for node `i`, component `j`, element `k` can be 626 /// found in the Lvector at index 627 /// `i*strides[0] + j*strides[1] + k*strides[2]`. 628 /// CEED_STRIDES_BACKEND may be used with vectors created 629 /// by a Ceed backend. 630 /// 631 /// ``` 632 /// # use libceed::prelude::*; 633 /// # fn main() -> libceed::Result<()> { 634 /// # let ceed = libceed::Ceed::default_init(); 635 /// let nelem = 3; 636 /// let strides: [i32; 3] = [1, 2, 2]; 637 /// let r = ceed.strided_elem_restriction(nelem, 2, 1, nelem * 2, strides)?; 638 /// # Ok(()) 639 /// # } 640 /// ``` 641 pub fn strided_elem_restriction<'a>( 642 &self, 643 nelem: usize, 644 elemsize: usize, 645 ncomp: usize, 646 lsize: usize, 647 strides: [i32; 3], 648 ) -> Result<ElemRestriction<'a>> { 649 ElemRestriction::create_strided(self, nelem, elemsize, ncomp, lsize, strides) 650 } 651 652 /// Returns an $H^1$ tensor-product Basis 653 /// 654 /// # arguments 655 /// 656 /// * `dim` - Topological dimension of element 657 /// * `ncomp` - Number of field components (1 for scalar fields) 658 /// * `P1d` - Number of Gauss-Lobatto nodes in one dimension. The 659 /// polynomial degree of the resulting `Q_k` element is 660 /// `k=P-1`. 661 /// * `Q1d` - Number of quadrature points in one dimension 662 /// * `interp1d` - Row-major `(Q1d * P1d)` matrix expressing the values of 663 /// nodal basis functions at quadrature points 664 /// * `grad1d` - Row-major `(Q1d * P1d)` matrix expressing derivatives of 665 /// nodal basis functions at quadrature points 666 /// * `qref1d` - Array of length `Q1d` holding the locations of quadrature 667 /// points on the 1D reference element `[-1, 1]` 668 /// * `qweight1d` - Array of length `Q1d` holding the quadrature weights on 669 /// the reference element 670 /// 671 /// ``` 672 /// # use libceed::prelude::*; 673 /// # fn main() -> libceed::Result<()> { 674 /// # let ceed = libceed::Ceed::default_init(); 675 /// let interp1d = [ 0.62994317, 0.47255875, -0.14950343, 0.04700152, 676 /// -0.07069480, 0.97297619, 0.13253993, -0.03482132, 677 /// -0.03482132, 0.13253993, 0.97297619, -0.07069480, 678 /// 0.04700152, -0.14950343, 0.47255875, 0.62994317]; 679 /// let grad1d = [-2.34183742, 2.78794489, -0.63510411, 0.18899664, 680 /// -0.51670214, -0.48795249, 1.33790510, -0.33325047, 681 // 0.33325047, -1.33790510, 0.48795249, 0.51670214, 682 /// -0.18899664, 0.63510411, -2.78794489, 2.34183742]; 683 /// let qref1d = [-0.86113631, -0.33998104, 0.33998104, 0.86113631]; 684 /// let qweight1d = [ 0.34785485, 0.65214515, 0.65214515, 0.34785485]; 685 /// let b = ceed. 686 /// basis_tensor_H1(2, 1, 4, 4, &interp1d, &grad1d, &qref1d, &qweight1d)?; 687 /// # Ok(()) 688 /// # } 689 /// ``` 690 pub fn basis_tensor_H1<'a>( 691 &self, 692 dim: usize, 693 ncomp: usize, 694 P1d: usize, 695 Q1d: usize, 696 interp1d: &[crate::Scalar], 697 grad1d: &[crate::Scalar], 698 qref1d: &[crate::Scalar], 699 qweight1d: &[crate::Scalar], 700 ) -> Result<Basis<'a>> { 701 Basis::create_tensor_H1( 702 self, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, 703 ) 704 } 705 706 /// Returns an $H^1$ Lagrange tensor-product Basis 707 /// 708 /// # arguments 709 /// 710 /// * `dim` - Topological dimension of element 711 /// * `ncomp` - Number of field components (1 for scalar fields) 712 /// * `P` - Number of Gauss-Lobatto nodes in one dimension. The 713 /// polynomial degree of the resulting `Q_k` element is `k=P-1`. 714 /// * `Q` - Number of quadrature points in one dimension 715 /// * `qmode` - Distribution of the `Q` quadrature points (affects order of 716 /// accuracy for the quadrature) 717 /// 718 /// ``` 719 /// # use libceed::prelude::*; 720 /// # fn main() -> libceed::Result<()> { 721 /// # let ceed = libceed::Ceed::default_init(); 722 /// let b = ceed.basis_tensor_H1_Lagrange(2, 1, 3, 4, QuadMode::Gauss)?; 723 /// # Ok(()) 724 /// # } 725 /// ``` 726 pub fn basis_tensor_H1_Lagrange<'a>( 727 &self, 728 dim: usize, 729 ncomp: usize, 730 P: usize, 731 Q: usize, 732 qmode: QuadMode, 733 ) -> Result<Basis<'a>> { 734 Basis::create_tensor_H1_Lagrange(self, dim, ncomp, P, Q, qmode) 735 } 736 737 /// Returns an $H-1$ Basis 738 /// 739 /// # arguments 740 /// 741 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 742 /// * `ncomp` - Number of field components (1 for scalar fields) 743 /// * `nnodes` - Total number of nodes 744 /// * `nqpts` - Total number of quadrature points 745 /// * `interp` - Row-major `(nqpts * nnodes)` matrix expressing the values of 746 /// nodal basis functions at quadrature points 747 /// * `grad` - Row-major `(dim * nqpts * nnodes)` matrix expressing 748 /// derivatives of nodal basis functions at quadrature points 749 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 750 /// points on the reference element `[-1, 1]` 751 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 752 /// the reference element 753 /// 754 /// ``` 755 /// # use libceed::prelude::*; 756 /// # fn main() -> libceed::Result<()> { 757 /// # let ceed = libceed::Ceed::default_init(); 758 /// let interp = [ 759 /// 0.12000000, 760 /// 0.48000000, 761 /// -0.12000000, 762 /// 0.48000000, 763 /// 0.16000000, 764 /// -0.12000000, 765 /// -0.12000000, 766 /// 0.48000000, 767 /// 0.12000000, 768 /// 0.16000000, 769 /// 0.48000000, 770 /// -0.12000000, 771 /// -0.11111111, 772 /// 0.44444444, 773 /// -0.11111111, 774 /// 0.44444444, 775 /// 0.44444444, 776 /// -0.11111111, 777 /// -0.12000000, 778 /// 0.16000000, 779 /// -0.12000000, 780 /// 0.48000000, 781 /// 0.48000000, 782 /// 0.12000000, 783 /// ]; 784 /// let grad = [ 785 /// -1.40000000, 786 /// 1.60000000, 787 /// -0.20000000, 788 /// -0.80000000, 789 /// 0.80000000, 790 /// 0.00000000, 791 /// 0.20000000, 792 /// -1.60000000, 793 /// 1.40000000, 794 /// -0.80000000, 795 /// 0.80000000, 796 /// 0.00000000, 797 /// -0.33333333, 798 /// 0.00000000, 799 /// 0.33333333, 800 /// -1.33333333, 801 /// 1.33333333, 802 /// 0.00000000, 803 /// 0.20000000, 804 /// 0.00000000, 805 /// -0.20000000, 806 /// -2.40000000, 807 /// 2.40000000, 808 /// 0.00000000, 809 /// -1.40000000, 810 /// -0.80000000, 811 /// 0.00000000, 812 /// 1.60000000, 813 /// 0.80000000, 814 /// -0.20000000, 815 /// 0.20000000, 816 /// -2.40000000, 817 /// 0.00000000, 818 /// 0.00000000, 819 /// 2.40000000, 820 /// -0.20000000, 821 /// -0.33333333, 822 /// -1.33333333, 823 /// 0.00000000, 824 /// 0.00000000, 825 /// 1.33333333, 826 /// 0.33333333, 827 /// 0.20000000, 828 /// -0.80000000, 829 /// 0.00000000, 830 /// -1.60000000, 831 /// 0.80000000, 832 /// 1.40000000, 833 /// ]; 834 /// let qref = [ 835 /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, 836 /// 0.60000000, 837 /// ]; 838 /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; 839 /// let b = ceed.basis_H1( 840 /// ElemTopology::Triangle, 841 /// 1, 842 /// 6, 843 /// 4, 844 /// &interp, 845 /// &grad, 846 /// &qref, 847 /// &qweight, 848 /// )?; 849 /// # Ok(()) 850 /// # } 851 /// ``` 852 pub fn basis_H1<'a>( 853 &self, 854 topo: ElemTopology, 855 ncomp: usize, 856 nnodes: usize, 857 nqpts: usize, 858 interp: &[crate::Scalar], 859 grad: &[crate::Scalar], 860 qref: &[crate::Scalar], 861 qweight: &[crate::Scalar], 862 ) -> Result<Basis<'a>> { 863 Basis::create_H1( 864 self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, 865 ) 866 } 867 868 /// Returns an $H(div)$ Basis 869 /// 870 /// # arguments 871 /// 872 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 873 /// * `ncomp` - Number of field components (1 for scalar fields) 874 /// * `nnodes` - Total number of nodes 875 /// * `nqpts` - Total number of quadrature points 876 /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the 877 /// values of basis functions at quadrature points 878 /// * `div` - Row-major `(nqpts * nnodes)` matrix expressing the 879 /// divergence of basis functions at quadrature points 880 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 881 /// points on the reference element `[-1, 1]` 882 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 883 /// the reference element 884 /// 885 /// ``` 886 /// # use libceed::prelude::*; 887 /// # fn main() -> libceed::Result<()> { 888 /// # let ceed = libceed::Ceed::default_init(); 889 /// let interp = [ 890 /// 0.00000000, 891 /// -1.57735027, 892 /// 0.57735027, 893 /// 0.00000000, 894 /// 0.00000000, 895 /// -1.57735027, 896 /// 0.57735027, 897 /// 0.00000000, 898 /// 0.00000000, 899 /// -1.57735027, 900 /// 0.57735027, 901 /// 0.00000000, 902 /// 0.00000000, 903 /// -0.42264973, 904 /// -0.57735027, 905 /// 0.00000000, 906 /// 0.42264973, 907 /// 0.00000000, 908 /// 0.00000000, 909 /// 0.57735027, 910 /// 0.42264973, 911 /// 0.00000000, 912 /// 0.00000000, 913 /// 0.57735027, 914 /// 1.57735027, 915 /// 0.00000000, 916 /// 0.00000000, 917 /// -0.57735027, 918 /// 1.57735027, 919 /// 0.00000000, 920 /// 0.00000000, 921 /// -0.57735027, 922 /// ]; 923 /// let div = [ 924 /// -1.00000000, 925 /// 1.00000000, 926 /// -1.00000000, 927 /// 1.00000000, 928 /// -1.00000000, 929 /// 1.00000000, 930 /// -1.00000000, 931 /// 1.00000000, 932 /// -1.00000000, 933 /// 1.00000000, 934 /// -1.00000000, 935 /// 1.00000000, 936 /// -1.00000000, 937 /// 1.00000000, 938 /// -1.00000000, 939 /// 1.00000000, 940 /// ]; 941 /// let qref = [ 942 /// 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 943 /// 0.57735026, 944 /// ]; 945 /// let qweight = [1.00000000, 1.00000000, 1.00000000, 1.00000000]; 946 /// let b = ceed.basis_Hdiv(ElemTopology::Quad, 1, 4, 4, &interp, &div, &qref, &qweight)?; 947 /// # Ok(()) 948 /// # } 949 /// ``` 950 pub fn basis_Hdiv<'a>( 951 &self, 952 topo: ElemTopology, 953 ncomp: usize, 954 nnodes: usize, 955 nqpts: usize, 956 interp: &[crate::Scalar], 957 div: &[crate::Scalar], 958 qref: &[crate::Scalar], 959 qweight: &[crate::Scalar], 960 ) -> Result<Basis<'a>> { 961 Basis::create_Hdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight) 962 } 963 964 /// Returns an $H(curl)$ Basis 965 /// 966 /// # arguments 967 /// 968 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 969 /// * `ncomp` - Number of field components (1 for scalar fields) 970 /// * `nnodes` - Total number of nodes 971 /// * `nqpts` - Total number of quadrature points 972 /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the 973 /// values of basis functions at quadrature points 974 /// * `curl` - Row-major `(curl_comp * nqpts * nnodes)`, `curl_comp = 1 if 975 /// dim < 3 else dim` matrix expressing the curl of basis 976 /// functions at quadrature points 977 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 978 /// points on the reference element `[-1, 1]` 979 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 980 /// the reference element 981 /// 982 /// ``` 983 /// # use libceed::prelude::*; 984 /// # fn main() -> libceed::Result<()> { 985 /// # let ceed = libceed::Ceed::default_init(); 986 /// let interp = [ 987 /// -0.20000000, 988 /// 0.20000000, 989 /// 0.80000000, 990 /// -0.20000000, 991 /// 0.20000000, 992 /// 0.80000000, 993 /// -0.33333333, 994 /// 0.33333333, 995 /// 0.66666667, 996 /// -0.60000000, 997 /// 0.60000000, 998 /// 0.40000000, 999 /// 0.20000000, 1000 /// 0.80000000, 1001 /// 0.20000000, 1002 /// 0.60000000, 1003 /// 0.40000000, 1004 /// 0.60000000, 1005 /// 0.33333333, 1006 /// 0.66666667, 1007 /// 0.33333333, 1008 /// 0.20000000, 1009 /// 0.80000000, 1010 /// 0.20000000, 1011 /// ]; 1012 /// let curl = [ 1013 /// 2.00000000, 1014 /// -2.00000000, 1015 /// -2.00000000, 1016 /// 2.00000000, 1017 /// -2.00000000, 1018 /// -2.00000000, 1019 /// 2.00000000, 1020 /// -2.00000000, 1021 /// -2.00000000, 1022 /// 2.00000000, 1023 /// -2.00000000, 1024 /// -2.00000000, 1025 /// ]; 1026 /// let qref = [ 1027 /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, 1028 /// 0.60000000, 1029 /// ]; 1030 /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; 1031 /// let b = ceed.basis_Hcurl( 1032 /// ElemTopology::Triangle, 1033 /// 1, 1034 /// 3, 1035 /// 4, 1036 /// &interp, 1037 /// &curl, 1038 /// &qref, 1039 /// &qweight, 1040 /// )?; 1041 /// # Ok(()) 1042 /// # } 1043 /// ``` 1044 pub fn basis_Hcurl<'a>( 1045 &self, 1046 topo: ElemTopology, 1047 ncomp: usize, 1048 nnodes: usize, 1049 nqpts: usize, 1050 interp: &[crate::Scalar], 1051 curl: &[crate::Scalar], 1052 qref: &[crate::Scalar], 1053 qweight: &[crate::Scalar], 1054 ) -> Result<Basis<'a>> { 1055 Basis::create_Hcurl( 1056 self, topo, ncomp, nnodes, nqpts, interp, curl, qref, qweight, 1057 ) 1058 } 1059 1060 /// Returns a QFunction for evaluating interior (volumetric) terms 1061 /// of a weak form corresponding to the $L^2$ inner product 1062 /// 1063 /// $$ 1064 /// \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), 1065 /// $$ 1066 /// 1067 /// where $v \cdot f_0$ represents contraction over fields and $\nabla v : f_1$ 1068 /// represents contraction over both fields and spatial dimensions. 1069 /// 1070 /// # arguments 1071 /// 1072 /// * `vlength` - Vector length. Caller must ensure that number of 1073 /// quadrature points is a multiple of vlength. 1074 /// * `f` - Boxed closure to evaluate weak form at quadrature points. 1075 /// 1076 /// ``` 1077 /// # use libceed::prelude::*; 1078 /// # fn main() -> libceed::Result<()> { 1079 /// # let ceed = libceed::Ceed::default_init(); 1080 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 1081 /// // Iterate over quadrature points 1082 /// v.iter_mut() 1083 /// .zip(u.iter().zip(weights.iter())) 1084 /// .for_each(|(v, (u, w))| *v = u * w); 1085 /// 1086 /// // Return clean error code 1087 /// 0 1088 /// }; 1089 /// 1090 /// let qf = ceed.q_function_interior(1, Box::new(user_f))?; 1091 /// # Ok(()) 1092 /// # } 1093 /// ``` 1094 pub fn q_function_interior<'a>( 1095 &self, 1096 vlength: usize, 1097 f: Box<qfunction::QFunctionUserClosure>, 1098 ) -> Result<QFunction<'a>> { 1099 QFunction::create(self, vlength, f) 1100 } 1101 1102 /// Returns a QFunction for evaluating interior (volumetric) terms 1103 /// created by name 1104 /// 1105 /// # arguments 1106 /// 1107 /// * `name` - name of QFunction from libCEED gallery 1108 /// 1109 /// ``` 1110 /// # use libceed::prelude::*; 1111 /// # fn main() -> libceed::Result<()> { 1112 /// # let ceed = libceed::Ceed::default_init(); 1113 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1114 /// # Ok(()) 1115 /// # } 1116 /// ``` 1117 pub fn q_function_interior_by_name<'a>(&self, name: &str) -> Result<QFunctionByName<'a>> { 1118 QFunctionByName::create(self, name) 1119 } 1120 1121 /// Returns an Operator and associate a QFunction. A Basis and 1122 /// ElemRestriction can be associated with QFunction fields via 1123 /// set_field(). 1124 /// 1125 /// # arguments 1126 /// 1127 /// * `qf` - QFunction defining the action of the operator at quadrature 1128 /// points 1129 /// * `dqf` - QFunction defining the action of the Jacobian of the qf (or 1130 /// qfunction_none) 1131 /// * `dqfT` - QFunction defining the action of the transpose of the 1132 /// Jacobian of the qf (or qfunction_none) 1133 /// 1134 /// ``` 1135 /// # use libceed::prelude::*; 1136 /// # fn main() -> libceed::Result<()> { 1137 /// # let ceed = libceed::Ceed::default_init(); 1138 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1139 /// let op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 1140 /// # Ok(()) 1141 /// # } 1142 /// ``` 1143 pub fn operator<'a, 'b>( 1144 &self, 1145 qf: impl Into<QFunctionOpt<'b>>, 1146 dqf: impl Into<QFunctionOpt<'b>>, 1147 dqfT: impl Into<QFunctionOpt<'b>>, 1148 ) -> Result<Operator<'a>> { 1149 Operator::create(self, qf, dqf, dqfT) 1150 } 1151 1152 /// Returns an Operator that composes the action of several Operators 1153 /// 1154 /// ``` 1155 /// # use libceed::prelude::*; 1156 /// # fn main() -> libceed::Result<()> { 1157 /// # let ceed = libceed::Ceed::default_init(); 1158 /// let op = ceed.composite_operator()?; 1159 /// # Ok(()) 1160 /// # } 1161 /// ``` 1162 pub fn composite_operator<'a>(&self) -> Result<CompositeOperator<'a>> { 1163 CompositeOperator::create(self) 1164 } 1165 } 1166 1167 // ----------------------------------------------------------------------------- 1168 // Tests 1169 // ----------------------------------------------------------------------------- 1170 #[cfg(test)] 1171 mod tests { 1172 use super::*; 1173 1174 fn ceed_t501() -> Result<()> { 1175 let resource = "/cpu/self/ref/blocked"; 1176 let ceed = Ceed::init(resource); 1177 let nelem = 4; 1178 let p = 3; 1179 let q = 4; 1180 let ndofs = p * nelem - nelem + 1; 1181 1182 // Vectors 1183 let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1184 let mut qdata = ceed.vector(nelem * q)?; 1185 qdata.set_value(0.0)?; 1186 let mut u = ceed.vector(ndofs)?; 1187 u.set_value(1.0)?; 1188 let mut v = ceed.vector(ndofs)?; 1189 v.set_value(0.0)?; 1190 1191 // Restrictions 1192 let mut indx: Vec<i32> = vec![0; 2 * nelem]; 1193 for i in 0..nelem { 1194 indx[2 * i + 0] = i as i32; 1195 indx[2 * i + 1] = (i + 1) as i32; 1196 } 1197 let rx = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &indx)?; 1198 let mut indu: Vec<i32> = vec![0; p * nelem]; 1199 for i in 0..nelem { 1200 indu[p * i + 0] = i as i32; 1201 indu[p * i + 1] = (i + 1) as i32; 1202 indu[p * i + 2] = (i + 2) as i32; 1203 } 1204 let ru = ceed.elem_restriction(nelem, 3, 1, 1, ndofs, MemType::Host, &indu)?; 1205 let strides: [i32; 3] = [1, q as i32, q as i32]; 1206 let rq = ceed.strided_elem_restriction(nelem, q, 1, q * nelem, strides)?; 1207 1208 // Bases 1209 let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1210 let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 1211 1212 // Build quadrature data 1213 let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1214 ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1215 .field("dx", &rx, &bx, VectorOpt::Active)? 1216 .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1217 .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1218 .apply(&x, &mut qdata)?; 1219 1220 // Mass operator 1221 let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1222 let op_mass = ceed 1223 .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1224 .field("u", &ru, &bu, VectorOpt::Active)? 1225 .field("qdata", &rq, BasisOpt::None, &qdata)? 1226 .field("v", &ru, &bu, VectorOpt::Active)? 1227 .check()?; 1228 1229 v.set_value(0.0)?; 1230 op_mass.apply(&u, &mut v)?; 1231 1232 // Check 1233 let sum: Scalar = v.view()?.iter().sum(); 1234 assert!( 1235 (sum - 2.0).abs() < 1e-15, 1236 "Incorrect interval length computed" 1237 ); 1238 Ok(()) 1239 } 1240 1241 #[test] 1242 fn test_ceed_t501() { 1243 assert!(ceed_t501().is_ok()); 1244 } 1245 } 1246 1247 // ----------------------------------------------------------------------------- 1248