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