1 // Copyright (c) 2017-2024, 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 //! A Ceed Operator defines the finite/spectral element operator associated to a 9 //! Ceed QFunction. A Ceed Operator connects Ceed ElemRestrictions, 10 //! Ceed Bases, and Ceed QFunctions. 11 12 use crate::prelude::*; 13 14 // ----------------------------------------------------------------------------- 15 // Operator Field context wrapper 16 // ----------------------------------------------------------------------------- 17 #[derive(Debug)] 18 pub struct OperatorField<'a> { 19 pub(crate) ptr: bind_ceed::CeedOperatorField, 20 pub(crate) vector: crate::Vector<'a>, 21 pub(crate) elem_restriction: crate::ElemRestriction<'a>, 22 pub(crate) basis: crate::Basis<'a>, 23 _lifeline: PhantomData<&'a ()>, 24 } 25 26 // ----------------------------------------------------------------------------- 27 // Implementations 28 // ----------------------------------------------------------------------------- 29 impl<'a> OperatorField<'a> { 30 pub(crate) fn from_raw( 31 ptr: bind_ceed::CeedOperatorField, 32 ceed: crate::Ceed, 33 ) -> crate::Result<Self> { 34 let vector = { 35 let mut vector_ptr = std::ptr::null_mut(); 36 ceed.check_error(unsafe { 37 bind_ceed::CeedOperatorFieldGetVector(ptr, &mut vector_ptr) 38 })?; 39 crate::Vector::from_raw(vector_ptr)? 40 }; 41 let elem_restriction = { 42 let mut elem_restriction_ptr = std::ptr::null_mut(); 43 ceed.check_error(unsafe { 44 bind_ceed::CeedOperatorFieldGetElemRestriction(ptr, &mut elem_restriction_ptr) 45 })?; 46 crate::ElemRestriction::from_raw(elem_restriction_ptr)? 47 }; 48 let basis = { 49 let mut basis_ptr = std::ptr::null_mut(); 50 ceed.check_error(unsafe { bind_ceed::CeedOperatorFieldGetBasis(ptr, &mut basis_ptr) })?; 51 crate::Basis::from_raw(basis_ptr)? 52 }; 53 Ok(Self { 54 ptr, 55 vector, 56 elem_restriction, 57 basis, 58 _lifeline: PhantomData, 59 }) 60 } 61 62 /// Get the name of an OperatorField 63 /// 64 /// ``` 65 /// # use libceed::prelude::*; 66 /// # fn main() -> libceed::Result<()> { 67 /// # let ceed = libceed::Ceed::default_init(); 68 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 69 /// 70 /// // Operator field arguments 71 /// let ne = 3; 72 /// let q = 4 as usize; 73 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 74 /// for i in 0..ne { 75 /// ind[2 * i + 0] = i as i32; 76 /// ind[2 * i + 1] = (i + 1) as i32; 77 /// } 78 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 79 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 80 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 81 /// 82 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 83 /// 84 /// // Operator fields 85 /// let op = ceed 86 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 87 /// .field("dx", &r, &b, VectorOpt::Active)? 88 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 89 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 90 /// 91 /// let inputs = op.inputs()?; 92 /// 93 /// assert_eq!(inputs[0].name(), "dx", "Incorrect input name"); 94 /// assert_eq!(inputs[1].name(), "weights", "Incorrect input name"); 95 /// # Ok(()) 96 /// # } 97 /// ``` 98 pub fn name(&self) -> &str { 99 let mut name_ptr: *mut std::os::raw::c_char = std::ptr::null_mut(); 100 unsafe { 101 bind_ceed::CeedOperatorFieldGetName( 102 self.ptr, 103 &mut name_ptr as *const _ as *mut *const _, 104 ); 105 } 106 unsafe { CStr::from_ptr(name_ptr) }.to_str().unwrap() 107 } 108 109 /// Get the ElemRestriction of an OperatorField 110 /// 111 /// ``` 112 /// # use libceed::prelude::*; 113 /// # fn main() -> libceed::Result<()> { 114 /// # let ceed = libceed::Ceed::default_init(); 115 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 116 /// 117 /// // Operator field arguments 118 /// let ne = 3; 119 /// let q = 4 as usize; 120 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 121 /// for i in 0..ne { 122 /// ind[2 * i + 0] = i as i32; 123 /// ind[2 * i + 1] = (i + 1) as i32; 124 /// } 125 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 126 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 127 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 128 /// 129 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 130 /// 131 /// // Operator fields 132 /// let op = ceed 133 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 134 /// .field("dx", &r, &b, VectorOpt::Active)? 135 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 136 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 137 /// 138 /// let inputs = op.inputs()?; 139 /// 140 /// assert!( 141 /// inputs[0].elem_restriction().is_some(), 142 /// "Incorrect field ElemRestriction" 143 /// ); 144 /// if let ElemRestrictionOpt::Some(r) = inputs[0].elem_restriction() { 145 /// assert_eq!( 146 /// r.num_elements(), 147 /// ne, 148 /// "Incorrect field ElemRestriction number of elements" 149 /// ); 150 /// } 151 /// 152 /// assert!( 153 /// inputs[1].elem_restriction().is_none(), 154 /// "Incorrect field ElemRestriction" 155 /// ); 156 /// 157 /// let outputs = op.outputs()?; 158 /// 159 /// assert!( 160 /// outputs[0].elem_restriction().is_some(), 161 /// "Incorrect field ElemRestriction" 162 /// ); 163 /// if let ElemRestrictionOpt::Some(r) = outputs[0].elem_restriction() { 164 /// assert_eq!( 165 /// r.num_elements(), 166 /// ne, 167 /// "Incorrect field ElemRestriction number of elements" 168 /// ); 169 /// } 170 /// # Ok(()) 171 /// # } 172 /// ``` 173 pub fn elem_restriction(&self) -> ElemRestrictionOpt { 174 if self.elem_restriction.ptr == unsafe { bind_ceed::CEED_ELEMRESTRICTION_NONE } { 175 ElemRestrictionOpt::None 176 } else { 177 ElemRestrictionOpt::Some(&self.elem_restriction) 178 } 179 } 180 181 /// Get the Basis of an OperatorField 182 /// 183 /// ``` 184 /// # use libceed::prelude::*; 185 /// # fn main() -> libceed::Result<()> { 186 /// # let ceed = libceed::Ceed::default_init(); 187 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 188 /// 189 /// // Operator field arguments 190 /// let ne = 3; 191 /// let q = 4 as usize; 192 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 193 /// for i in 0..ne { 194 /// ind[2 * i + 0] = i as i32; 195 /// ind[2 * i + 1] = (i + 1) as i32; 196 /// } 197 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 198 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 199 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 200 /// 201 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 202 /// 203 /// // Operator fields 204 /// let op = ceed 205 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 206 /// .field("dx", &r, &b, VectorOpt::Active)? 207 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 208 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 209 /// 210 /// let inputs = op.inputs()?; 211 /// 212 /// assert!(inputs[0].basis().is_some(), "Incorrect field Basis"); 213 /// if let BasisOpt::Some(b) = inputs[0].basis() { 214 /// assert_eq!( 215 /// b.num_quadrature_points(), 216 /// q, 217 /// "Incorrect field Basis number of quadrature points" 218 /// ); 219 /// } 220 /// assert!(inputs[1].basis().is_some(), "Incorrect field Basis"); 221 /// if let BasisOpt::Some(b) = inputs[1].basis() { 222 /// assert_eq!( 223 /// b.num_quadrature_points(), 224 /// q, 225 /// "Incorrect field Basis number of quadrature points" 226 /// ); 227 /// } 228 /// 229 /// let outputs = op.outputs()?; 230 /// 231 /// assert!(outputs[0].basis().is_none(), "Incorrect field Basis"); 232 /// # Ok(()) 233 /// # } 234 /// ``` 235 pub fn basis(&self) -> BasisOpt { 236 if self.basis.ptr == unsafe { bind_ceed::CEED_BASIS_NONE } { 237 BasisOpt::None 238 } else { 239 BasisOpt::Some(&self.basis) 240 } 241 } 242 243 /// Get the Vector of an OperatorField 244 /// 245 /// ``` 246 /// # use libceed::prelude::*; 247 /// # fn main() -> libceed::Result<()> { 248 /// # let ceed = libceed::Ceed::default_init(); 249 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 250 /// 251 /// // Operator field arguments 252 /// let ne = 3; 253 /// let q = 4 as usize; 254 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 255 /// for i in 0..ne { 256 /// ind[2 * i + 0] = i as i32; 257 /// ind[2 * i + 1] = (i + 1) as i32; 258 /// } 259 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 260 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 261 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 262 /// 263 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 264 /// 265 /// // Operator fields 266 /// let op = ceed 267 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 268 /// .field("dx", &r, &b, VectorOpt::Active)? 269 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 270 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 271 /// 272 /// let inputs = op.inputs()?; 273 /// 274 /// assert!(inputs[0].vector().is_active(), "Incorrect field Vector"); 275 /// assert!(inputs[1].vector().is_none(), "Incorrect field Vector"); 276 /// 277 /// let outputs = op.outputs()?; 278 /// 279 /// assert!(outputs[0].vector().is_active(), "Incorrect field Vector"); 280 /// # Ok(()) 281 /// # } 282 /// ``` 283 pub fn vector(&self) -> VectorOpt { 284 if self.vector.ptr == unsafe { bind_ceed::CEED_VECTOR_ACTIVE } { 285 VectorOpt::Active 286 } else if self.vector.ptr == unsafe { bind_ceed::CEED_VECTOR_NONE } { 287 VectorOpt::None 288 } else { 289 VectorOpt::Some(&self.vector) 290 } 291 } 292 } 293 294 // ----------------------------------------------------------------------------- 295 // Operator context wrapper 296 // ----------------------------------------------------------------------------- 297 #[derive(Debug)] 298 pub(crate) struct OperatorCore<'a> { 299 ptr: bind_ceed::CeedOperator, 300 _lifeline: PhantomData<&'a ()>, 301 } 302 303 #[derive(Debug)] 304 pub struct Operator<'a> { 305 op_core: OperatorCore<'a>, 306 } 307 308 #[derive(Debug)] 309 pub struct CompositeOperator<'a> { 310 op_core: OperatorCore<'a>, 311 } 312 313 // ----------------------------------------------------------------------------- 314 // Destructor 315 // ----------------------------------------------------------------------------- 316 impl<'a> Drop for OperatorCore<'a> { 317 fn drop(&mut self) { 318 unsafe { 319 bind_ceed::CeedOperatorDestroy(&mut self.ptr); 320 } 321 } 322 } 323 324 // ----------------------------------------------------------------------------- 325 // Display 326 // ----------------------------------------------------------------------------- 327 impl<'a> fmt::Display for OperatorCore<'a> { 328 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 329 let mut ptr = std::ptr::null_mut(); 330 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 331 let cstring = unsafe { 332 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 333 bind_ceed::CeedOperatorView(self.ptr, file); 334 bind_ceed::fclose(file); 335 CString::from_raw(ptr) 336 }; 337 cstring.to_string_lossy().fmt(f) 338 } 339 } 340 341 /// View an Operator 342 /// 343 /// ``` 344 /// # use libceed::prelude::*; 345 /// # fn main() -> libceed::Result<()> { 346 /// # let ceed = libceed::Ceed::default_init(); 347 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 348 /// 349 /// // Operator field arguments 350 /// let ne = 3; 351 /// let q = 4 as usize; 352 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 353 /// for i in 0..ne { 354 /// ind[2 * i + 0] = i as i32; 355 /// ind[2 * i + 1] = (i + 1) as i32; 356 /// } 357 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 358 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 359 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 360 /// 361 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 362 /// 363 /// // Operator fields 364 /// let op = ceed 365 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 366 /// .name("mass")? 367 /// .field("dx", &r, &b, VectorOpt::Active)? 368 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 369 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 370 /// 371 /// println!("{}", op); 372 /// # Ok(()) 373 /// # } 374 /// ``` 375 impl<'a> fmt::Display for Operator<'a> { 376 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 377 self.op_core.fmt(f) 378 } 379 } 380 381 /// View a composite Operator 382 /// 383 /// ``` 384 /// # use libceed::prelude::*; 385 /// # fn main() -> libceed::Result<()> { 386 /// # let ceed = libceed::Ceed::default_init(); 387 /// 388 /// // Sub operator field arguments 389 /// let ne = 3; 390 /// let q = 4 as usize; 391 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 392 /// for i in 0..ne { 393 /// ind[2 * i + 0] = i as i32; 394 /// ind[2 * i + 1] = (i + 1) as i32; 395 /// } 396 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 397 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 398 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 399 /// 400 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 401 /// 402 /// let qdata_mass = ceed.vector(q * ne)?; 403 /// let qdata_diff = ceed.vector(q * ne)?; 404 /// 405 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 406 /// let op_mass = ceed 407 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 408 /// .name("Mass term")? 409 /// .field("u", &r, &b, VectorOpt::Active)? 410 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 411 /// .field("v", &r, &b, VectorOpt::Active)?; 412 /// 413 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 414 /// let op_diff = ceed 415 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 416 /// .name("Poisson term")? 417 /// .field("du", &r, &b, VectorOpt::Active)? 418 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 419 /// .field("dv", &r, &b, VectorOpt::Active)?; 420 /// 421 /// let op = ceed 422 /// .composite_operator()? 423 /// .name("Screened Poisson")? 424 /// .sub_operator(&op_mass)? 425 /// .sub_operator(&op_diff)?; 426 /// 427 /// println!("{}", op); 428 /// # Ok(()) 429 /// # } 430 /// ``` 431 impl<'a> fmt::Display for CompositeOperator<'a> { 432 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 433 self.op_core.fmt(f) 434 } 435 } 436 437 // ----------------------------------------------------------------------------- 438 // Core functionality 439 // ----------------------------------------------------------------------------- 440 impl<'a> OperatorCore<'a> { 441 // Raw Ceed for error handling 442 #[doc(hidden)] 443 fn ceed(&self) -> bind_ceed::Ceed { 444 unsafe { bind_ceed::CeedOperatorReturnCeed(self.ptr) } 445 } 446 447 // Error handling 448 #[doc(hidden)] 449 fn check_error(&self, ierr: i32) -> crate::Result<i32> { 450 crate::check_error(|| self.ceed(), ierr) 451 } 452 453 // Common implementations 454 pub fn check(&self) -> crate::Result<i32> { 455 self.check_error(unsafe { bind_ceed::CeedOperatorCheckReady(self.ptr) }) 456 } 457 458 pub fn name(&self, name: &str) -> crate::Result<i32> { 459 let name_c = CString::new(name).expect("CString::new failed"); 460 self.check_error(unsafe { bind_ceed::CeedOperatorSetName(self.ptr, name_c.as_ptr()) }) 461 } 462 463 pub fn apply(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 464 self.check_error(unsafe { 465 bind_ceed::CeedOperatorApply( 466 self.ptr, 467 input.ptr, 468 output.ptr, 469 bind_ceed::CEED_REQUEST_IMMEDIATE, 470 ) 471 }) 472 } 473 474 pub fn apply_add(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 475 self.check_error(unsafe { 476 bind_ceed::CeedOperatorApplyAdd( 477 self.ptr, 478 input.ptr, 479 output.ptr, 480 bind_ceed::CEED_REQUEST_IMMEDIATE, 481 ) 482 }) 483 } 484 485 pub fn linear_assemble_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 486 self.check_error(unsafe { 487 bind_ceed::CeedOperatorLinearAssembleDiagonal( 488 self.ptr, 489 assembled.ptr, 490 bind_ceed::CEED_REQUEST_IMMEDIATE, 491 ) 492 }) 493 } 494 495 pub fn linear_assemble_add_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 496 self.check_error(unsafe { 497 bind_ceed::CeedOperatorLinearAssembleAddDiagonal( 498 self.ptr, 499 assembled.ptr, 500 bind_ceed::CEED_REQUEST_IMMEDIATE, 501 ) 502 }) 503 } 504 505 pub fn linear_assemble_point_block_diagonal( 506 &self, 507 assembled: &mut Vector, 508 ) -> crate::Result<i32> { 509 self.check_error(unsafe { 510 bind_ceed::CeedOperatorLinearAssemblePointBlockDiagonal( 511 self.ptr, 512 assembled.ptr, 513 bind_ceed::CEED_REQUEST_IMMEDIATE, 514 ) 515 }) 516 } 517 518 pub fn linear_assemble_add_point_block_diagonal( 519 &self, 520 assembled: &mut Vector, 521 ) -> crate::Result<i32> { 522 self.check_error(unsafe { 523 bind_ceed::CeedOperatorLinearAssembleAddPointBlockDiagonal( 524 self.ptr, 525 assembled.ptr, 526 bind_ceed::CEED_REQUEST_IMMEDIATE, 527 ) 528 }) 529 } 530 } 531 532 // ----------------------------------------------------------------------------- 533 // Operator 534 // ----------------------------------------------------------------------------- 535 impl<'a> Operator<'a> { 536 // Constructor 537 pub fn create<'b>( 538 ceed: &crate::Ceed, 539 qf: impl Into<QFunctionOpt<'b>>, 540 dqf: impl Into<QFunctionOpt<'b>>, 541 dqfT: impl Into<QFunctionOpt<'b>>, 542 ) -> crate::Result<Self> { 543 let mut ptr = std::ptr::null_mut(); 544 ceed.check_error(unsafe { 545 bind_ceed::CeedOperatorCreate( 546 ceed.ptr, 547 qf.into().to_raw(), 548 dqf.into().to_raw(), 549 dqfT.into().to_raw(), 550 &mut ptr, 551 ) 552 })?; 553 Ok(Self { 554 op_core: OperatorCore { 555 ptr, 556 _lifeline: PhantomData, 557 }, 558 }) 559 } 560 561 fn from_raw(ptr: bind_ceed::CeedOperator) -> crate::Result<Self> { 562 Ok(Self { 563 op_core: OperatorCore { 564 ptr, 565 _lifeline: PhantomData, 566 }, 567 }) 568 } 569 570 /// Set name for Operator printing 571 /// 572 /// * 'name' - Name to set 573 /// 574 /// ``` 575 /// # use libceed::prelude::*; 576 /// # fn main() -> libceed::Result<()> { 577 /// # let ceed = libceed::Ceed::default_init(); 578 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 579 /// 580 /// // Operator field arguments 581 /// let ne = 3; 582 /// let q = 4 as usize; 583 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 584 /// for i in 0..ne { 585 /// ind[2 * i + 0] = i as i32; 586 /// ind[2 * i + 1] = (i + 1) as i32; 587 /// } 588 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 589 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 590 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 591 /// 592 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 593 /// 594 /// // Operator fields 595 /// let op = ceed 596 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 597 /// .name("mass")? 598 /// .field("dx", &r, &b, VectorOpt::Active)? 599 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 600 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 601 /// # Ok(()) 602 /// # } 603 /// ``` 604 #[allow(unused_mut)] 605 pub fn name(mut self, name: &str) -> crate::Result<Self> { 606 self.op_core.name(name)?; 607 Ok(self) 608 } 609 610 /// Apply Operator to a vector 611 /// 612 /// * `input` - Input Vector 613 /// * `output` - Output Vector 614 /// 615 /// ``` 616 /// # use libceed::prelude::*; 617 /// # fn main() -> libceed::Result<()> { 618 /// # let ceed = libceed::Ceed::default_init(); 619 /// let ne = 4; 620 /// let p = 3; 621 /// let q = 4; 622 /// let ndofs = p * ne - ne + 1; 623 /// 624 /// // Vectors 625 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 626 /// let mut qdata = ceed.vector(ne * q)?; 627 /// qdata.set_value(0.0); 628 /// let u = ceed.vector_from_slice(&vec![1.0; ndofs])?; 629 /// let mut v = ceed.vector(ndofs)?; 630 /// v.set_value(0.0); 631 /// 632 /// // Restrictions 633 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 634 /// for i in 0..ne { 635 /// indx[2 * i + 0] = i as i32; 636 /// indx[2 * i + 1] = (i + 1) as i32; 637 /// } 638 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 639 /// let mut indu: Vec<i32> = vec![0; p * ne]; 640 /// for i in 0..ne { 641 /// indu[p * i + 0] = i as i32; 642 /// indu[p * i + 1] = (i + 1) as i32; 643 /// indu[p * i + 2] = (i + 2) as i32; 644 /// } 645 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 646 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 647 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 648 /// 649 /// // Bases 650 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 651 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 652 /// 653 /// // Build quadrature data 654 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 655 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 656 /// .field("dx", &rx, &bx, VectorOpt::Active)? 657 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 658 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 659 /// .apply(&x, &mut qdata)?; 660 /// 661 /// // Mass operator 662 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 663 /// let op_mass = ceed 664 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 665 /// .field("u", &ru, &bu, VectorOpt::Active)? 666 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 667 /// .field("v", &ru, &bu, VectorOpt::Active)?; 668 /// 669 /// v.set_value(0.0); 670 /// op_mass.apply(&u, &mut v)?; 671 /// 672 /// // Check 673 /// let sum: Scalar = v.view()?.iter().sum(); 674 /// let error: Scalar = (sum - 2.0).abs(); 675 /// assert!( 676 /// error < 50.0 * libceed::EPSILON, 677 /// "Incorrect interval length computed. Expected: 2.0, Found: {}, Error: {:.12e}", 678 /// sum, 679 /// error 680 /// ); 681 /// # Ok(()) 682 /// # } 683 /// ``` 684 pub fn apply(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 685 self.op_core.apply(input, output) 686 } 687 688 /// Apply Operator to a vector and add result to output vector 689 /// 690 /// * `input` - Input Vector 691 /// * `output` - Output Vector 692 /// 693 /// ``` 694 /// # use libceed::prelude::*; 695 /// # fn main() -> libceed::Result<()> { 696 /// # let ceed = libceed::Ceed::default_init(); 697 /// let ne = 4; 698 /// let p = 3; 699 /// let q = 4; 700 /// let ndofs = p * ne - ne + 1; 701 /// 702 /// // Vectors 703 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 704 /// let mut qdata = ceed.vector(ne * q)?; 705 /// qdata.set_value(0.0); 706 /// let u = ceed.vector_from_slice(&vec![1.0; ndofs])?; 707 /// let mut v = ceed.vector(ndofs)?; 708 /// 709 /// // Restrictions 710 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 711 /// for i in 0..ne { 712 /// indx[2 * i + 0] = i as i32; 713 /// indx[2 * i + 1] = (i + 1) as i32; 714 /// } 715 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 716 /// let mut indu: Vec<i32> = vec![0; p * ne]; 717 /// for i in 0..ne { 718 /// indu[p * i + 0] = i as i32; 719 /// indu[p * i + 1] = (i + 1) as i32; 720 /// indu[p * i + 2] = (i + 2) as i32; 721 /// } 722 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 723 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 724 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 725 /// 726 /// // Bases 727 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 728 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 729 /// 730 /// // Build quadrature data 731 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 732 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 733 /// .field("dx", &rx, &bx, VectorOpt::Active)? 734 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 735 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 736 /// .apply(&x, &mut qdata)?; 737 /// 738 /// // Mass operator 739 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 740 /// let op_mass = ceed 741 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 742 /// .field("u", &ru, &bu, VectorOpt::Active)? 743 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 744 /// .field("v", &ru, &bu, VectorOpt::Active)?; 745 /// 746 /// v.set_value(1.0); 747 /// op_mass.apply_add(&u, &mut v)?; 748 /// 749 /// // Check 750 /// let sum: Scalar = v.view()?.iter().sum(); 751 /// assert!( 752 /// (sum - (2.0 + ndofs as Scalar)).abs() < 10.0 * libceed::EPSILON, 753 /// "Incorrect interval length computed and added" 754 /// ); 755 /// # Ok(()) 756 /// # } 757 /// ``` 758 pub fn apply_add(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 759 self.op_core.apply_add(input, output) 760 } 761 762 /// Provide a field to a Operator for use by its QFunction 763 /// 764 /// * `fieldname` - Name of the field (to be matched with the name used by 765 /// the QFunction) 766 /// * `r` - ElemRestriction 767 /// * `b` - Basis in which the field resides or `BasisOpt::None` 768 /// * `v` - Vector to be used by Operator or `VectorOpt::Active` if field 769 /// is active or `VectorOpt::None` if using `Weight` with the 770 /// QFunction 771 /// 772 /// 773 /// ``` 774 /// # use libceed::prelude::*; 775 /// # fn main() -> libceed::Result<()> { 776 /// # let ceed = libceed::Ceed::default_init(); 777 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 778 /// let mut op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 779 /// 780 /// // Operator field arguments 781 /// let ne = 3; 782 /// let q = 4; 783 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 784 /// for i in 0..ne { 785 /// ind[2 * i + 0] = i as i32; 786 /// ind[2 * i + 1] = (i + 1) as i32; 787 /// } 788 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 789 /// 790 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 791 /// 792 /// // Operator field 793 /// op = op.field("dx", &r, &b, VectorOpt::Active)?; 794 /// # Ok(()) 795 /// # } 796 /// ``` 797 #[allow(unused_mut)] 798 pub fn field<'b>( 799 mut self, 800 fieldname: &str, 801 r: impl Into<ElemRestrictionOpt<'b>>, 802 b: impl Into<BasisOpt<'b>>, 803 v: impl Into<VectorOpt<'b>>, 804 ) -> crate::Result<Self> { 805 let fieldname = CString::new(fieldname).expect("CString::new failed"); 806 let fieldname = fieldname.as_ptr() as *const i8; 807 self.op_core.check_error(unsafe { 808 bind_ceed::CeedOperatorSetField( 809 self.op_core.ptr, 810 fieldname, 811 r.into().to_raw(), 812 b.into().to_raw(), 813 v.into().to_raw(), 814 ) 815 })?; 816 Ok(self) 817 } 818 819 /// Get a slice of Operator inputs 820 /// 821 /// ``` 822 /// # use libceed::prelude::*; 823 /// # fn main() -> libceed::Result<()> { 824 /// # let ceed = libceed::Ceed::default_init(); 825 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 826 /// 827 /// // Operator field arguments 828 /// let ne = 3; 829 /// let q = 4 as usize; 830 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 831 /// for i in 0..ne { 832 /// ind[2 * i + 0] = i as i32; 833 /// ind[2 * i + 1] = (i + 1) as i32; 834 /// } 835 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 836 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 837 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 838 /// 839 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 840 /// 841 /// // Operator fields 842 /// let op = ceed 843 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 844 /// .field("dx", &r, &b, VectorOpt::Active)? 845 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 846 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 847 /// 848 /// let inputs = op.inputs()?; 849 /// 850 /// assert_eq!(inputs.len(), 2, "Incorrect inputs array"); 851 /// # Ok(()) 852 /// # } 853 /// ``` 854 pub fn inputs(&self) -> crate::Result<Vec<crate::OperatorField>> { 855 // Get array of raw C pointers for inputs 856 let mut num_inputs = 0; 857 let mut inputs_ptr = std::ptr::null_mut(); 858 self.op_core.check_error(unsafe { 859 bind_ceed::CeedOperatorGetFields( 860 self.op_core.ptr, 861 &mut num_inputs, 862 &mut inputs_ptr, 863 std::ptr::null_mut() as *mut bind_ceed::CeedInt, 864 std::ptr::null_mut() as *mut *mut bind_ceed::CeedOperatorField, 865 ) 866 })?; 867 // Convert raw C pointers to fixed length slice 868 let inputs_slice = unsafe { 869 std::slice::from_raw_parts( 870 inputs_ptr as *mut bind_ceed::CeedOperatorField, 871 num_inputs as usize, 872 ) 873 }; 874 // And finally build vec 875 let ceed = { 876 let ceed_raw = self.op_core.ceed(); 877 let mut ptr = std::ptr::null_mut(); 878 unsafe { 879 bind_ceed::CeedReferenceCopy(ceed_raw, &mut ptr); // refcount 880 } 881 crate::Ceed { ptr } 882 }; 883 let inputs = (0..num_inputs as usize) 884 .map(|i| crate::OperatorField::from_raw(inputs_slice[i], ceed.clone())) 885 .collect::<crate::Result<Vec<_>>>()?; 886 Ok(inputs) 887 } 888 889 /// Get a slice of Operator outputs 890 /// 891 /// ``` 892 /// # use libceed::prelude::*; 893 /// # fn main() -> libceed::Result<()> { 894 /// # let ceed = libceed::Ceed::default_init(); 895 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 896 /// 897 /// // Operator field arguments 898 /// let ne = 3; 899 /// let q = 4 as usize; 900 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 901 /// for i in 0..ne { 902 /// ind[2 * i + 0] = i as i32; 903 /// ind[2 * i + 1] = (i + 1) as i32; 904 /// } 905 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 906 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 907 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 908 /// 909 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 910 /// 911 /// // Operator fields 912 /// let op = ceed 913 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 914 /// .field("dx", &r, &b, VectorOpt::Active)? 915 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 916 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 917 /// 918 /// let outputs = op.outputs()?; 919 /// 920 /// assert_eq!(outputs.len(), 1, "Incorrect outputs array"); 921 /// # Ok(()) 922 /// # } 923 /// ``` 924 pub fn outputs(&self) -> crate::Result<Vec<crate::OperatorField>> { 925 // Get array of raw C pointers for outputs 926 let mut num_outputs = 0; 927 let mut outputs_ptr = std::ptr::null_mut(); 928 self.op_core.check_error(unsafe { 929 bind_ceed::CeedOperatorGetFields( 930 self.op_core.ptr, 931 std::ptr::null_mut() as *mut bind_ceed::CeedInt, 932 std::ptr::null_mut() as *mut *mut bind_ceed::CeedOperatorField, 933 &mut num_outputs, 934 &mut outputs_ptr, 935 ) 936 })?; 937 // Convert raw C pointers to fixed length slice 938 let outputs_slice = unsafe { 939 std::slice::from_raw_parts( 940 outputs_ptr as *mut bind_ceed::CeedOperatorField, 941 num_outputs as usize, 942 ) 943 }; 944 // And finally build vec 945 let ceed = { 946 let ceed_raw = self.op_core.ceed(); 947 let mut ptr = std::ptr::null_mut(); 948 unsafe { 949 bind_ceed::CeedReferenceCopy(ceed_raw, &mut ptr); // refcount 950 } 951 crate::Ceed { ptr } 952 }; 953 let outputs = (0..num_outputs as usize) 954 .map(|i| crate::OperatorField::from_raw(outputs_slice[i], ceed.clone())) 955 .collect::<crate::Result<Vec<_>>>()?; 956 Ok(outputs) 957 } 958 959 /// Check if Operator is setup correctly 960 /// 961 /// ``` 962 /// # use libceed::prelude::*; 963 /// # fn main() -> libceed::Result<()> { 964 /// # let ceed = libceed::Ceed::default_init(); 965 /// let ne = 4; 966 /// let p = 3; 967 /// let q = 4; 968 /// let ndofs = p * ne - ne + 1; 969 /// 970 /// // Restrictions 971 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 972 /// for i in 0..ne { 973 /// indx[2 * i + 0] = i as i32; 974 /// indx[2 * i + 1] = (i + 1) as i32; 975 /// } 976 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 977 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 978 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 979 /// 980 /// // Bases 981 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 982 /// 983 /// // Build quadrature data 984 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 985 /// let op_build = ceed 986 /// .operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 987 /// .field("dx", &rx, &bx, VectorOpt::Active)? 988 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 989 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 990 /// 991 /// // Check 992 /// assert!(op_build.check().is_ok(), "Operator setup incorrectly"); 993 /// # Ok(()) 994 /// # } 995 /// ``` 996 pub fn check(self) -> crate::Result<Self> { 997 self.op_core.check()?; 998 Ok(self) 999 } 1000 1001 /// Get the number of elements associated with an Operator 1002 /// 1003 /// 1004 /// ``` 1005 /// # use libceed::prelude::*; 1006 /// # fn main() -> libceed::Result<()> { 1007 /// # let ceed = libceed::Ceed::default_init(); 1008 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1009 /// let mut op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 1010 /// 1011 /// // Operator field arguments 1012 /// let ne = 3; 1013 /// let q = 4; 1014 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 1015 /// for i in 0..ne { 1016 /// ind[2 * i + 0] = i as i32; 1017 /// ind[2 * i + 1] = (i + 1) as i32; 1018 /// } 1019 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 1020 /// 1021 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1022 /// 1023 /// // Operator field 1024 /// op = op.field("dx", &r, &b, VectorOpt::Active)?; 1025 /// 1026 /// // Check number of elements 1027 /// assert_eq!(op.num_elements(), ne, "Incorrect number of elements"); 1028 /// # Ok(()) 1029 /// # } 1030 /// ``` 1031 pub fn num_elements(&self) -> usize { 1032 let mut nelem = 0; 1033 unsafe { bind_ceed::CeedOperatorGetNumElements(self.op_core.ptr, &mut nelem) }; 1034 usize::try_from(nelem).unwrap() 1035 } 1036 1037 /// Get the number of quadrature points associated with each element of 1038 /// an Operator 1039 /// 1040 /// 1041 /// ``` 1042 /// # use libceed::prelude::*; 1043 /// # fn main() -> libceed::Result<()> { 1044 /// # let ceed = libceed::Ceed::default_init(); 1045 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1046 /// let mut op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 1047 /// 1048 /// // Operator field arguments 1049 /// let ne = 3; 1050 /// let q = 4; 1051 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 1052 /// for i in 0..ne { 1053 /// ind[2 * i + 0] = i as i32; 1054 /// ind[2 * i + 1] = (i + 1) as i32; 1055 /// } 1056 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 1057 /// 1058 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1059 /// 1060 /// // Operator field 1061 /// op = op.field("dx", &r, &b, VectorOpt::Active)?; 1062 /// 1063 /// // Check number of quadrature points 1064 /// assert_eq!( 1065 /// op.num_quadrature_points(), 1066 /// q, 1067 /// "Incorrect number of quadrature points" 1068 /// ); 1069 /// # Ok(()) 1070 /// # } 1071 /// ``` 1072 pub fn num_quadrature_points(&self) -> usize { 1073 let mut Q = 0; 1074 unsafe { bind_ceed::CeedOperatorGetNumQuadraturePoints(self.op_core.ptr, &mut Q) }; 1075 usize::try_from(Q).unwrap() 1076 } 1077 1078 /// Assemble the diagonal of a square linear Operator 1079 /// 1080 /// This overwrites a Vector with the diagonal of a linear Operator. 1081 /// 1082 /// Note: Currently only non-composite Operators with a single field and 1083 /// composite Operators with single field sub-operators are supported. 1084 /// 1085 /// * `op` - Operator to assemble QFunction 1086 /// * `assembled` - Vector to store assembled Operator diagonal 1087 /// 1088 /// ``` 1089 /// # use libceed::prelude::*; 1090 /// # fn main() -> libceed::Result<()> { 1091 /// # let ceed = libceed::Ceed::default_init(); 1092 /// let ne = 4; 1093 /// let p = 3; 1094 /// let q = 4; 1095 /// let ndofs = p * ne - ne + 1; 1096 /// 1097 /// // Vectors 1098 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1099 /// let mut qdata = ceed.vector(ne * q)?; 1100 /// qdata.set_value(0.0); 1101 /// let mut u = ceed.vector(ndofs)?; 1102 /// u.set_value(1.0); 1103 /// let mut v = ceed.vector(ndofs)?; 1104 /// v.set_value(0.0); 1105 /// 1106 /// // Restrictions 1107 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1108 /// for i in 0..ne { 1109 /// indx[2 * i + 0] = i as i32; 1110 /// indx[2 * i + 1] = (i + 1) as i32; 1111 /// } 1112 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1113 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1114 /// for i in 0..ne { 1115 /// indu[p * i + 0] = (2 * i) as i32; 1116 /// indu[p * i + 1] = (2 * i + 1) as i32; 1117 /// indu[p * i + 2] = (2 * i + 2) as i32; 1118 /// } 1119 /// let ru = ceed.elem_restriction(ne, p, 1, 1, ndofs, MemType::Host, &indu)?; 1120 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1121 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1122 /// 1123 /// // Bases 1124 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1125 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 1126 /// 1127 /// // Build quadrature data 1128 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1129 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1130 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1131 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1132 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1133 /// .apply(&x, &mut qdata)?; 1134 /// 1135 /// // Mass operator 1136 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1137 /// let op_mass = ceed 1138 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1139 /// .field("u", &ru, &bu, VectorOpt::Active)? 1140 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1141 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1142 /// 1143 /// // Diagonal 1144 /// let mut diag = ceed.vector(ndofs)?; 1145 /// diag.set_value(0.0); 1146 /// op_mass.linear_assemble_diagonal(&mut diag)?; 1147 /// 1148 /// // Manual diagonal computation 1149 /// let mut true_diag = ceed.vector(ndofs)?; 1150 /// true_diag.set_value(0.0)?; 1151 /// for i in 0..ndofs { 1152 /// u.set_value(0.0); 1153 /// { 1154 /// let mut u_array = u.view_mut()?; 1155 /// u_array[i] = 1.; 1156 /// } 1157 /// 1158 /// op_mass.apply(&u, &mut v)?; 1159 /// 1160 /// { 1161 /// let v_array = v.view()?; 1162 /// let mut true_array = true_diag.view_mut()?; 1163 /// true_array[i] = v_array[i]; 1164 /// } 1165 /// } 1166 /// 1167 /// // Check 1168 /// diag.view()? 1169 /// .iter() 1170 /// .zip(true_diag.view()?.iter()) 1171 /// .for_each(|(computed, actual)| { 1172 /// assert!( 1173 /// (*computed - *actual).abs() < 10.0 * libceed::EPSILON, 1174 /// "Diagonal entry incorrect" 1175 /// ); 1176 /// }); 1177 /// # Ok(()) 1178 /// # } 1179 /// ``` 1180 pub fn linear_assemble_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 1181 self.op_core.linear_assemble_diagonal(assembled) 1182 } 1183 1184 /// Assemble the diagonal of a square linear Operator 1185 /// 1186 /// This sums into a Vector with the diagonal of a linear Operator. 1187 /// 1188 /// Note: Currently only non-composite Operators with a single field and 1189 /// composite Operators with single field sub-operators are supported. 1190 /// 1191 /// * `op` - Operator to assemble QFunction 1192 /// * `assembled` - Vector to store assembled Operator diagonal 1193 /// 1194 /// 1195 /// ``` 1196 /// # use libceed::prelude::*; 1197 /// # fn main() -> libceed::Result<()> { 1198 /// # let ceed = libceed::Ceed::default_init(); 1199 /// let ne = 4; 1200 /// let p = 3; 1201 /// let q = 4; 1202 /// let ndofs = p * ne - ne + 1; 1203 /// 1204 /// // Vectors 1205 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1206 /// let mut qdata = ceed.vector(ne * q)?; 1207 /// qdata.set_value(0.0); 1208 /// let mut u = ceed.vector(ndofs)?; 1209 /// u.set_value(1.0); 1210 /// let mut v = ceed.vector(ndofs)?; 1211 /// v.set_value(0.0); 1212 /// 1213 /// // Restrictions 1214 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1215 /// for i in 0..ne { 1216 /// indx[2 * i + 0] = i as i32; 1217 /// indx[2 * i + 1] = (i + 1) as i32; 1218 /// } 1219 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1220 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1221 /// for i in 0..ne { 1222 /// indu[p * i + 0] = (2 * i) as i32; 1223 /// indu[p * i + 1] = (2 * i + 1) as i32; 1224 /// indu[p * i + 2] = (2 * i + 2) as i32; 1225 /// } 1226 /// let ru = ceed.elem_restriction(ne, p, 1, 1, ndofs, MemType::Host, &indu)?; 1227 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1228 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1229 /// 1230 /// // Bases 1231 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1232 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 1233 /// 1234 /// // Build quadrature data 1235 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1236 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1237 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1238 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1239 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1240 /// .apply(&x, &mut qdata)?; 1241 /// 1242 /// // Mass operator 1243 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1244 /// let op_mass = ceed 1245 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1246 /// .field("u", &ru, &bu, VectorOpt::Active)? 1247 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1248 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1249 /// 1250 /// // Diagonal 1251 /// let mut diag = ceed.vector(ndofs)?; 1252 /// diag.set_value(1.0); 1253 /// op_mass.linear_assemble_add_diagonal(&mut diag)?; 1254 /// 1255 /// // Manual diagonal computation 1256 /// let mut true_diag = ceed.vector(ndofs)?; 1257 /// true_diag.set_value(0.0)?; 1258 /// for i in 0..ndofs { 1259 /// u.set_value(0.0); 1260 /// { 1261 /// let mut u_array = u.view_mut()?; 1262 /// u_array[i] = 1.; 1263 /// } 1264 /// 1265 /// op_mass.apply(&u, &mut v)?; 1266 /// 1267 /// { 1268 /// let v_array = v.view()?; 1269 /// let mut true_array = true_diag.view_mut()?; 1270 /// true_array[i] = v_array[i] + 1.0; 1271 /// } 1272 /// } 1273 /// 1274 /// // Check 1275 /// diag.view()? 1276 /// .iter() 1277 /// .zip(true_diag.view()?.iter()) 1278 /// .for_each(|(computed, actual)| { 1279 /// assert!( 1280 /// (*computed - *actual).abs() < 10.0 * libceed::EPSILON, 1281 /// "Diagonal entry incorrect" 1282 /// ); 1283 /// }); 1284 /// # Ok(()) 1285 /// # } 1286 /// ``` 1287 pub fn linear_assemble_add_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 1288 self.op_core.linear_assemble_add_diagonal(assembled) 1289 } 1290 1291 /// Assemble the point block diagonal of a square linear Operator 1292 /// 1293 /// This overwrites a Vector with the point block diagonal of a linear 1294 /// Operator. 1295 /// 1296 /// Note: Currently only non-composite Operators with a single field and 1297 /// composite Operators with single field sub-operators are supported. 1298 /// 1299 /// * `op` - Operator to assemble QFunction 1300 /// * `assembled` - Vector to store assembled CeedOperator point block 1301 /// diagonal, provided in row-major form with an 1302 /// `ncomp * ncomp` block at each node. The dimensions of 1303 /// this vector are derived from the active vector for 1304 /// the CeedOperator. The array has shape 1305 /// `[nodes, component out, component in]`. 1306 /// 1307 /// ``` 1308 /// # use libceed::prelude::*; 1309 /// # fn main() -> libceed::Result<()> { 1310 /// # let ceed = libceed::Ceed::default_init(); 1311 /// let ne = 4; 1312 /// let p = 3; 1313 /// let q = 4; 1314 /// let ncomp = 2; 1315 /// let ndofs = p * ne - ne + 1; 1316 /// 1317 /// // Vectors 1318 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1319 /// let mut qdata = ceed.vector(ne * q)?; 1320 /// qdata.set_value(0.0); 1321 /// let mut u = ceed.vector(ncomp * ndofs)?; 1322 /// u.set_value(1.0); 1323 /// let mut v = ceed.vector(ncomp * ndofs)?; 1324 /// v.set_value(0.0); 1325 /// 1326 /// // Restrictions 1327 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1328 /// for i in 0..ne { 1329 /// indx[2 * i + 0] = i as i32; 1330 /// indx[2 * i + 1] = (i + 1) as i32; 1331 /// } 1332 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1333 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1334 /// for i in 0..ne { 1335 /// indu[p * i + 0] = (2 * i) as i32; 1336 /// indu[p * i + 1] = (2 * i + 1) as i32; 1337 /// indu[p * i + 2] = (2 * i + 2) as i32; 1338 /// } 1339 /// let ru = ceed.elem_restriction(ne, p, ncomp, ndofs, ncomp * ndofs, MemType::Host, &indu)?; 1340 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1341 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1342 /// 1343 /// // Bases 1344 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1345 /// let bu = ceed.basis_tensor_H1_Lagrange(1, ncomp, p, q, QuadMode::Gauss)?; 1346 /// 1347 /// // Build quadrature data 1348 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1349 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1350 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1351 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1352 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1353 /// .apply(&x, &mut qdata)?; 1354 /// 1355 /// // Mass operator 1356 /// let mut mass_2_comp = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 1357 /// // Number of quadrature points 1358 /// let q = qdata.len(); 1359 /// 1360 /// // Iterate over quadrature points 1361 /// for i in 0..q { 1362 /// v[i + 0 * q] = u[i + 1 * q] * qdata[i]; 1363 /// v[i + 1 * q] = u[i + 0 * q] * qdata[i]; 1364 /// } 1365 /// 1366 /// // Return clean error code 1367 /// 0 1368 /// }; 1369 /// 1370 /// let qf_mass = ceed 1371 /// .q_function_interior(1, Box::new(mass_2_comp))? 1372 /// .input("u", 2, EvalMode::Interp)? 1373 /// .input("qdata", 1, EvalMode::None)? 1374 /// .output("v", 2, EvalMode::Interp)?; 1375 /// 1376 /// let op_mass = ceed 1377 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1378 /// .field("u", &ru, &bu, VectorOpt::Active)? 1379 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1380 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1381 /// 1382 /// // Diagonal 1383 /// let mut diag = ceed.vector(ncomp * ncomp * ndofs)?; 1384 /// diag.set_value(0.0); 1385 /// op_mass.linear_assemble_point_block_diagonal(&mut diag)?; 1386 /// 1387 /// // Manual diagonal computation 1388 /// let mut true_diag = ceed.vector(ncomp * ncomp * ndofs)?; 1389 /// true_diag.set_value(0.0)?; 1390 /// for i in 0..ndofs { 1391 /// for j in 0..ncomp { 1392 /// u.set_value(0.0); 1393 /// { 1394 /// let mut u_array = u.view_mut()?; 1395 /// u_array[i + j * ndofs] = 1.; 1396 /// } 1397 /// 1398 /// op_mass.apply(&u, &mut v)?; 1399 /// 1400 /// { 1401 /// let v_array = v.view()?; 1402 /// let mut true_array = true_diag.view_mut()?; 1403 /// for k in 0..ncomp { 1404 /// true_array[i * ncomp * ncomp + k * ncomp + j] = v_array[i + k * ndofs]; 1405 /// } 1406 /// } 1407 /// } 1408 /// } 1409 /// 1410 /// // Check 1411 /// diag.view()? 1412 /// .iter() 1413 /// .zip(true_diag.view()?.iter()) 1414 /// .for_each(|(computed, actual)| { 1415 /// assert!( 1416 /// (*computed - *actual).abs() < 10.0 * libceed::EPSILON, 1417 /// "Diagonal entry incorrect" 1418 /// ); 1419 /// }); 1420 /// # Ok(()) 1421 /// # } 1422 /// ``` 1423 pub fn linear_assemble_point_block_diagonal( 1424 &self, 1425 assembled: &mut Vector, 1426 ) -> crate::Result<i32> { 1427 self.op_core.linear_assemble_point_block_diagonal(assembled) 1428 } 1429 1430 /// Assemble the point block diagonal of a square linear Operator 1431 /// 1432 /// This sums into a Vector with the point block diagonal of a linear 1433 /// Operator. 1434 /// 1435 /// Note: Currently only non-composite Operators with a single field and 1436 /// composite Operators with single field sub-operators are supported. 1437 /// 1438 /// * `op` - Operator to assemble QFunction 1439 /// * `assembled` - Vector to store assembled CeedOperator point block 1440 /// diagonal, provided in row-major form with an 1441 /// `ncomp * ncomp` block at each node. The dimensions of 1442 /// this vector are derived from the active vector for 1443 /// the CeedOperator. The array has shape 1444 /// `[nodes, component out, component in]`. 1445 /// 1446 /// ``` 1447 /// # use libceed::prelude::*; 1448 /// # fn main() -> libceed::Result<()> { 1449 /// # let ceed = libceed::Ceed::default_init(); 1450 /// let ne = 4; 1451 /// let p = 3; 1452 /// let q = 4; 1453 /// let ncomp = 2; 1454 /// let ndofs = p * ne - ne + 1; 1455 /// 1456 /// // Vectors 1457 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1458 /// let mut qdata = ceed.vector(ne * q)?; 1459 /// qdata.set_value(0.0); 1460 /// let mut u = ceed.vector(ncomp * ndofs)?; 1461 /// u.set_value(1.0); 1462 /// let mut v = ceed.vector(ncomp * ndofs)?; 1463 /// v.set_value(0.0); 1464 /// 1465 /// // Restrictions 1466 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1467 /// for i in 0..ne { 1468 /// indx[2 * i + 0] = i as i32; 1469 /// indx[2 * i + 1] = (i + 1) as i32; 1470 /// } 1471 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1472 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1473 /// for i in 0..ne { 1474 /// indu[p * i + 0] = (2 * i) as i32; 1475 /// indu[p * i + 1] = (2 * i + 1) as i32; 1476 /// indu[p * i + 2] = (2 * i + 2) as i32; 1477 /// } 1478 /// let ru = ceed.elem_restriction(ne, p, ncomp, ndofs, ncomp * ndofs, MemType::Host, &indu)?; 1479 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1480 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1481 /// 1482 /// // Bases 1483 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1484 /// let bu = ceed.basis_tensor_H1_Lagrange(1, ncomp, p, q, QuadMode::Gauss)?; 1485 /// 1486 /// // Build quadrature data 1487 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1488 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1489 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1490 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1491 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1492 /// .apply(&x, &mut qdata)?; 1493 /// 1494 /// // Mass operator 1495 /// let mut mass_2_comp = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 1496 /// // Number of quadrature points 1497 /// let q = qdata.len(); 1498 /// 1499 /// // Iterate over quadrature points 1500 /// for i in 0..q { 1501 /// v[i + 0 * q] = u[i + 1 * q] * qdata[i]; 1502 /// v[i + 1 * q] = u[i + 0 * q] * qdata[i]; 1503 /// } 1504 /// 1505 /// // Return clean error code 1506 /// 0 1507 /// }; 1508 /// 1509 /// let qf_mass = ceed 1510 /// .q_function_interior(1, Box::new(mass_2_comp))? 1511 /// .input("u", 2, EvalMode::Interp)? 1512 /// .input("qdata", 1, EvalMode::None)? 1513 /// .output("v", 2, EvalMode::Interp)?; 1514 /// 1515 /// let op_mass = ceed 1516 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1517 /// .field("u", &ru, &bu, VectorOpt::Active)? 1518 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1519 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1520 /// 1521 /// // Diagonal 1522 /// let mut diag = ceed.vector(ncomp * ncomp * ndofs)?; 1523 /// diag.set_value(1.0); 1524 /// op_mass.linear_assemble_add_point_block_diagonal(&mut diag)?; 1525 /// 1526 /// // Manual diagonal computation 1527 /// let mut true_diag = ceed.vector(ncomp * ncomp * ndofs)?; 1528 /// true_diag.set_value(0.0)?; 1529 /// for i in 0..ndofs { 1530 /// for j in 0..ncomp { 1531 /// u.set_value(0.0); 1532 /// { 1533 /// let mut u_array = u.view_mut()?; 1534 /// u_array[i + j * ndofs] = 1.; 1535 /// } 1536 /// 1537 /// op_mass.apply(&u, &mut v)?; 1538 /// 1539 /// { 1540 /// let v_array = v.view()?; 1541 /// let mut true_array = true_diag.view_mut()?; 1542 /// for k in 0..ncomp { 1543 /// true_array[i * ncomp * ncomp + k * ncomp + j] = v_array[i + k * ndofs]; 1544 /// } 1545 /// } 1546 /// } 1547 /// } 1548 /// 1549 /// // Check 1550 /// diag.view()? 1551 /// .iter() 1552 /// .zip(true_diag.view()?.iter()) 1553 /// .for_each(|(computed, actual)| { 1554 /// assert!( 1555 /// (*computed - 1.0 - *actual).abs() < 10.0 * libceed::EPSILON, 1556 /// "Diagonal entry incorrect" 1557 /// ); 1558 /// }); 1559 /// # Ok(()) 1560 /// # } 1561 /// ``` 1562 pub fn linear_assemble_add_point_block_diagonal( 1563 &self, 1564 assembled: &mut Vector, 1565 ) -> crate::Result<i32> { 1566 self.op_core 1567 .linear_assemble_add_point_block_diagonal(assembled) 1568 } 1569 1570 /// Create a multigrid coarse Operator and level transfer Operators for a 1571 /// given Operator, creating the prolongation basis from the fine and 1572 /// coarse grid interpolation 1573 /// 1574 /// * `p_mult_fine` - Lvector multiplicity in parallel gather/scatter 1575 /// * `rstr_coarse` - Coarse grid restriction 1576 /// * `basis_coarse` - Coarse grid active vector basis 1577 /// 1578 /// ``` 1579 /// # use libceed::prelude::*; 1580 /// # fn main() -> libceed::Result<()> { 1581 /// # let ceed = libceed::Ceed::default_init(); 1582 /// let ne = 15; 1583 /// let p_coarse = 3; 1584 /// let p_fine = 5; 1585 /// let q = 6; 1586 /// let ndofs_coarse = p_coarse * ne - ne + 1; 1587 /// let ndofs_fine = p_fine * ne - ne + 1; 1588 /// 1589 /// // Vectors 1590 /// let x_array = (0..ne + 1) 1591 /// .map(|i| 2.0 * i as Scalar / ne as Scalar - 1.0) 1592 /// .collect::<Vec<Scalar>>(); 1593 /// let x = ceed.vector_from_slice(&x_array)?; 1594 /// let mut qdata = ceed.vector(ne * q)?; 1595 /// qdata.set_value(0.0); 1596 /// let mut u_coarse = ceed.vector(ndofs_coarse)?; 1597 /// u_coarse.set_value(1.0); 1598 /// let mut u_fine = ceed.vector(ndofs_fine)?; 1599 /// u_fine.set_value(1.0); 1600 /// let mut v_coarse = ceed.vector(ndofs_coarse)?; 1601 /// v_coarse.set_value(0.0); 1602 /// let mut v_fine = ceed.vector(ndofs_fine)?; 1603 /// v_fine.set_value(0.0); 1604 /// let mut multiplicity = ceed.vector(ndofs_fine)?; 1605 /// multiplicity.set_value(1.0); 1606 /// 1607 /// // Restrictions 1608 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1609 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1610 /// 1611 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1612 /// for i in 0..ne { 1613 /// indx[2 * i + 0] = i as i32; 1614 /// indx[2 * i + 1] = (i + 1) as i32; 1615 /// } 1616 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1617 /// 1618 /// let mut indu_coarse: Vec<i32> = vec![0; p_coarse * ne]; 1619 /// for i in 0..ne { 1620 /// for j in 0..p_coarse { 1621 /// indu_coarse[p_coarse * i + j] = (i + j) as i32; 1622 /// } 1623 /// } 1624 /// let ru_coarse = ceed.elem_restriction( 1625 /// ne, 1626 /// p_coarse, 1627 /// 1, 1628 /// 1, 1629 /// ndofs_coarse, 1630 /// MemType::Host, 1631 /// &indu_coarse, 1632 /// )?; 1633 /// 1634 /// let mut indu_fine: Vec<i32> = vec![0; p_fine * ne]; 1635 /// for i in 0..ne { 1636 /// for j in 0..p_fine { 1637 /// indu_fine[p_fine * i + j] = (i + j) as i32; 1638 /// } 1639 /// } 1640 /// let ru_fine = ceed.elem_restriction(ne, p_fine, 1, 1, ndofs_fine, MemType::Host, &indu_fine)?; 1641 /// 1642 /// // Bases 1643 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1644 /// let bu_coarse = ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, q, QuadMode::Gauss)?; 1645 /// let bu_fine = ceed.basis_tensor_H1_Lagrange(1, 1, p_fine, q, QuadMode::Gauss)?; 1646 /// 1647 /// // Build quadrature data 1648 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1649 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1650 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1651 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1652 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1653 /// .apply(&x, &mut qdata)?; 1654 /// 1655 /// // Mass operator 1656 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1657 /// let op_mass_fine = ceed 1658 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1659 /// .field("u", &ru_fine, &bu_fine, VectorOpt::Active)? 1660 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1661 /// .field("v", &ru_fine, &bu_fine, VectorOpt::Active)?; 1662 /// 1663 /// // Multigrid setup 1664 /// let (op_mass_coarse, op_prolong, op_restrict) = 1665 /// op_mass_fine.create_multigrid_level(&multiplicity, &ru_coarse, &bu_coarse)?; 1666 /// 1667 /// // Coarse problem 1668 /// u_coarse.set_value(1.0); 1669 /// op_mass_coarse.apply(&u_coarse, &mut v_coarse)?; 1670 /// 1671 /// // Check 1672 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1673 /// assert!( 1674 /// (sum - 2.0).abs() < 50.0 * libceed::EPSILON, 1675 /// "Incorrect interval length computed" 1676 /// ); 1677 /// 1678 /// // Prolong 1679 /// op_prolong.apply(&u_coarse, &mut u_fine)?; 1680 /// 1681 /// // Fine problem 1682 /// op_mass_fine.apply(&u_fine, &mut v_fine)?; 1683 /// 1684 /// // Check 1685 /// let sum: Scalar = v_fine.view()?.iter().sum(); 1686 /// assert!( 1687 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1688 /// "Incorrect interval length computed" 1689 /// ); 1690 /// 1691 /// // Restrict 1692 /// op_restrict.apply(&v_fine, &mut v_coarse)?; 1693 /// 1694 /// // Check 1695 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1696 /// assert!( 1697 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1698 /// "Incorrect interval length computed" 1699 /// ); 1700 /// # Ok(()) 1701 /// # } 1702 /// ``` 1703 pub fn create_multigrid_level<'b>( 1704 &self, 1705 p_mult_fine: &Vector, 1706 rstr_coarse: &ElemRestriction, 1707 basis_coarse: &Basis, 1708 ) -> crate::Result<(Operator<'b>, Operator<'b>, Operator<'b>)> { 1709 let mut ptr_coarse = std::ptr::null_mut(); 1710 let mut ptr_prolong = std::ptr::null_mut(); 1711 let mut ptr_restrict = std::ptr::null_mut(); 1712 self.op_core.check_error(unsafe { 1713 bind_ceed::CeedOperatorMultigridLevelCreate( 1714 self.op_core.ptr, 1715 p_mult_fine.ptr, 1716 rstr_coarse.ptr, 1717 basis_coarse.ptr, 1718 &mut ptr_coarse, 1719 &mut ptr_prolong, 1720 &mut ptr_restrict, 1721 ) 1722 })?; 1723 let op_coarse = Operator::from_raw(ptr_coarse)?; 1724 let op_prolong = Operator::from_raw(ptr_prolong)?; 1725 let op_restrict = Operator::from_raw(ptr_restrict)?; 1726 Ok((op_coarse, op_prolong, op_restrict)) 1727 } 1728 1729 /// Create a multigrid coarse Operator and level transfer Operators for a 1730 /// given Operator with a tensor basis for the active basis 1731 /// 1732 /// * `p_mult_fine` - Lvector multiplicity in parallel gather/scatter 1733 /// * `rstr_coarse` - Coarse grid restriction 1734 /// * `basis_coarse` - Coarse grid active vector basis 1735 /// * `interp_c_to_f` - Matrix for coarse to fine 1736 /// 1737 /// ``` 1738 /// # use libceed::prelude::*; 1739 /// # fn main() -> libceed::Result<()> { 1740 /// # let ceed = libceed::Ceed::default_init(); 1741 /// let ne = 15; 1742 /// let p_coarse = 3; 1743 /// let p_fine = 5; 1744 /// let q = 6; 1745 /// let ndofs_coarse = p_coarse * ne - ne + 1; 1746 /// let ndofs_fine = p_fine * ne - ne + 1; 1747 /// 1748 /// // Vectors 1749 /// let x_array = (0..ne + 1) 1750 /// .map(|i| 2.0 * i as Scalar / ne as Scalar - 1.0) 1751 /// .collect::<Vec<Scalar>>(); 1752 /// let x = ceed.vector_from_slice(&x_array)?; 1753 /// let mut qdata = ceed.vector(ne * q)?; 1754 /// qdata.set_value(0.0); 1755 /// let mut u_coarse = ceed.vector(ndofs_coarse)?; 1756 /// u_coarse.set_value(1.0); 1757 /// let mut u_fine = ceed.vector(ndofs_fine)?; 1758 /// u_fine.set_value(1.0); 1759 /// let mut v_coarse = ceed.vector(ndofs_coarse)?; 1760 /// v_coarse.set_value(0.0); 1761 /// let mut v_fine = ceed.vector(ndofs_fine)?; 1762 /// v_fine.set_value(0.0); 1763 /// let mut multiplicity = ceed.vector(ndofs_fine)?; 1764 /// multiplicity.set_value(1.0); 1765 /// 1766 /// // Restrictions 1767 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1768 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1769 /// 1770 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1771 /// for i in 0..ne { 1772 /// indx[2 * i + 0] = i as i32; 1773 /// indx[2 * i + 1] = (i + 1) as i32; 1774 /// } 1775 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1776 /// 1777 /// let mut indu_coarse: Vec<i32> = vec![0; p_coarse * ne]; 1778 /// for i in 0..ne { 1779 /// for j in 0..p_coarse { 1780 /// indu_coarse[p_coarse * i + j] = (i + j) as i32; 1781 /// } 1782 /// } 1783 /// let ru_coarse = ceed.elem_restriction( 1784 /// ne, 1785 /// p_coarse, 1786 /// 1, 1787 /// 1, 1788 /// ndofs_coarse, 1789 /// MemType::Host, 1790 /// &indu_coarse, 1791 /// )?; 1792 /// 1793 /// let mut indu_fine: Vec<i32> = vec![0; p_fine * ne]; 1794 /// for i in 0..ne { 1795 /// for j in 0..p_fine { 1796 /// indu_fine[p_fine * i + j] = (i + j) as i32; 1797 /// } 1798 /// } 1799 /// let ru_fine = ceed.elem_restriction(ne, p_fine, 1, 1, ndofs_fine, MemType::Host, &indu_fine)?; 1800 /// 1801 /// // Bases 1802 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1803 /// let bu_coarse = ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, q, QuadMode::Gauss)?; 1804 /// let bu_fine = ceed.basis_tensor_H1_Lagrange(1, 1, p_fine, q, QuadMode::Gauss)?; 1805 /// 1806 /// // Build quadrature data 1807 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1808 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1809 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1810 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1811 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1812 /// .apply(&x, &mut qdata)?; 1813 /// 1814 /// // Mass operator 1815 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1816 /// let op_mass_fine = ceed 1817 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1818 /// .field("u", &ru_fine, &bu_fine, VectorOpt::Active)? 1819 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1820 /// .field("v", &ru_fine, &bu_fine, VectorOpt::Active)?; 1821 /// 1822 /// // Multigrid setup 1823 /// let mut interp_c_to_f: Vec<Scalar> = vec![0.; p_coarse * p_fine]; 1824 /// { 1825 /// let mut coarse = ceed.vector(p_coarse)?; 1826 /// let mut fine = ceed.vector(p_fine)?; 1827 /// let basis_c_to_f = 1828 /// ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, p_fine, QuadMode::GaussLobatto)?; 1829 /// for i in 0..p_coarse { 1830 /// coarse.set_value(0.0); 1831 /// { 1832 /// let mut array = coarse.view_mut()?; 1833 /// array[i] = 1.; 1834 /// } 1835 /// basis_c_to_f.apply( 1836 /// 1, 1837 /// TransposeMode::NoTranspose, 1838 /// EvalMode::Interp, 1839 /// &coarse, 1840 /// &mut fine, 1841 /// )?; 1842 /// let array = fine.view()?; 1843 /// for j in 0..p_fine { 1844 /// interp_c_to_f[j * p_coarse + i] = array[j]; 1845 /// } 1846 /// } 1847 /// } 1848 /// let (op_mass_coarse, op_prolong, op_restrict) = op_mass_fine.create_multigrid_level_tensor_H1( 1849 /// &multiplicity, 1850 /// &ru_coarse, 1851 /// &bu_coarse, 1852 /// &interp_c_to_f, 1853 /// )?; 1854 /// 1855 /// // Coarse problem 1856 /// u_coarse.set_value(1.0); 1857 /// op_mass_coarse.apply(&u_coarse, &mut v_coarse)?; 1858 /// 1859 /// // Check 1860 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1861 /// assert!( 1862 /// (sum - 2.0).abs() < 10.0 * libceed::EPSILON, 1863 /// "Incorrect interval length computed" 1864 /// ); 1865 /// 1866 /// // Prolong 1867 /// op_prolong.apply(&u_coarse, &mut u_fine)?; 1868 /// 1869 /// // Fine problem 1870 /// op_mass_fine.apply(&u_fine, &mut v_fine)?; 1871 /// 1872 /// // Check 1873 /// let sum: Scalar = v_fine.view()?.iter().sum(); 1874 /// assert!( 1875 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1876 /// "Incorrect interval length computed" 1877 /// ); 1878 /// 1879 /// // Restrict 1880 /// op_restrict.apply(&v_fine, &mut v_coarse)?; 1881 /// 1882 /// // Check 1883 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1884 /// assert!( 1885 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1886 /// "Incorrect interval length computed" 1887 /// ); 1888 /// # Ok(()) 1889 /// # } 1890 /// ``` 1891 pub fn create_multigrid_level_tensor_H1<'b>( 1892 &self, 1893 p_mult_fine: &Vector, 1894 rstr_coarse: &ElemRestriction, 1895 basis_coarse: &Basis, 1896 interpCtoF: &Vec<Scalar>, 1897 ) -> crate::Result<(Operator<'b>, Operator<'b>, Operator<'b>)> { 1898 let mut ptr_coarse = std::ptr::null_mut(); 1899 let mut ptr_prolong = std::ptr::null_mut(); 1900 let mut ptr_restrict = std::ptr::null_mut(); 1901 self.op_core.check_error(unsafe { 1902 bind_ceed::CeedOperatorMultigridLevelCreateTensorH1( 1903 self.op_core.ptr, 1904 p_mult_fine.ptr, 1905 rstr_coarse.ptr, 1906 basis_coarse.ptr, 1907 interpCtoF.as_ptr(), 1908 &mut ptr_coarse, 1909 &mut ptr_prolong, 1910 &mut ptr_restrict, 1911 ) 1912 })?; 1913 let op_coarse = Operator::from_raw(ptr_coarse)?; 1914 let op_prolong = Operator::from_raw(ptr_prolong)?; 1915 let op_restrict = Operator::from_raw(ptr_restrict)?; 1916 Ok((op_coarse, op_prolong, op_restrict)) 1917 } 1918 1919 /// Create a multigrid coarse Operator and level transfer Operators for a 1920 /// given Operator with a non-tensor basis for the active basis 1921 /// 1922 /// * `p_mult_fine` - Lvector multiplicity in parallel gather/scatter 1923 /// * `rstr_coarse` - Coarse grid restriction 1924 /// * `basis_coarse` - Coarse grid active vector basis 1925 /// * `interp_c_to_f` - Matrix for coarse to fine 1926 /// 1927 /// ``` 1928 /// # use libceed::prelude::*; 1929 /// # fn main() -> libceed::Result<()> { 1930 /// # let ceed = libceed::Ceed::default_init(); 1931 /// let ne = 15; 1932 /// let p_coarse = 3; 1933 /// let p_fine = 5; 1934 /// let q = 6; 1935 /// let ndofs_coarse = p_coarse * ne - ne + 1; 1936 /// let ndofs_fine = p_fine * ne - ne + 1; 1937 /// 1938 /// // Vectors 1939 /// let x_array = (0..ne + 1) 1940 /// .map(|i| 2.0 * i as Scalar / ne as Scalar - 1.0) 1941 /// .collect::<Vec<Scalar>>(); 1942 /// let x = ceed.vector_from_slice(&x_array)?; 1943 /// let mut qdata = ceed.vector(ne * q)?; 1944 /// qdata.set_value(0.0); 1945 /// let mut u_coarse = ceed.vector(ndofs_coarse)?; 1946 /// u_coarse.set_value(1.0); 1947 /// let mut u_fine = ceed.vector(ndofs_fine)?; 1948 /// u_fine.set_value(1.0); 1949 /// let mut v_coarse = ceed.vector(ndofs_coarse)?; 1950 /// v_coarse.set_value(0.0); 1951 /// let mut v_fine = ceed.vector(ndofs_fine)?; 1952 /// v_fine.set_value(0.0); 1953 /// let mut multiplicity = ceed.vector(ndofs_fine)?; 1954 /// multiplicity.set_value(1.0); 1955 /// 1956 /// // Restrictions 1957 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1958 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1959 /// 1960 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1961 /// for i in 0..ne { 1962 /// indx[2 * i + 0] = i as i32; 1963 /// indx[2 * i + 1] = (i + 1) as i32; 1964 /// } 1965 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1966 /// 1967 /// let mut indu_coarse: Vec<i32> = vec![0; p_coarse * ne]; 1968 /// for i in 0..ne { 1969 /// for j in 0..p_coarse { 1970 /// indu_coarse[p_coarse * i + j] = (i + j) as i32; 1971 /// } 1972 /// } 1973 /// let ru_coarse = ceed.elem_restriction( 1974 /// ne, 1975 /// p_coarse, 1976 /// 1, 1977 /// 1, 1978 /// ndofs_coarse, 1979 /// MemType::Host, 1980 /// &indu_coarse, 1981 /// )?; 1982 /// 1983 /// let mut indu_fine: Vec<i32> = vec![0; p_fine * ne]; 1984 /// for i in 0..ne { 1985 /// for j in 0..p_fine { 1986 /// indu_fine[p_fine * i + j] = (i + j) as i32; 1987 /// } 1988 /// } 1989 /// let ru_fine = ceed.elem_restriction(ne, p_fine, 1, 1, ndofs_fine, MemType::Host, &indu_fine)?; 1990 /// 1991 /// // Bases 1992 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1993 /// let bu_coarse = ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, q, QuadMode::Gauss)?; 1994 /// let bu_fine = ceed.basis_tensor_H1_Lagrange(1, 1, p_fine, q, QuadMode::Gauss)?; 1995 /// 1996 /// // Build quadrature data 1997 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1998 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1999 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2000 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2001 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2002 /// .apply(&x, &mut qdata)?; 2003 /// 2004 /// // Mass operator 2005 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2006 /// let op_mass_fine = ceed 2007 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2008 /// .field("u", &ru_fine, &bu_fine, VectorOpt::Active)? 2009 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 2010 /// .field("v", &ru_fine, &bu_fine, VectorOpt::Active)?; 2011 /// 2012 /// // Multigrid setup 2013 /// let mut interp_c_to_f: Vec<Scalar> = vec![0.; p_coarse * p_fine]; 2014 /// { 2015 /// let mut coarse = ceed.vector(p_coarse)?; 2016 /// let mut fine = ceed.vector(p_fine)?; 2017 /// let basis_c_to_f = 2018 /// ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, p_fine, QuadMode::GaussLobatto)?; 2019 /// for i in 0..p_coarse { 2020 /// coarse.set_value(0.0); 2021 /// { 2022 /// let mut array = coarse.view_mut()?; 2023 /// array[i] = 1.; 2024 /// } 2025 /// basis_c_to_f.apply( 2026 /// 1, 2027 /// TransposeMode::NoTranspose, 2028 /// EvalMode::Interp, 2029 /// &coarse, 2030 /// &mut fine, 2031 /// )?; 2032 /// let array = fine.view()?; 2033 /// for j in 0..p_fine { 2034 /// interp_c_to_f[j * p_coarse + i] = array[j]; 2035 /// } 2036 /// } 2037 /// } 2038 /// let (op_mass_coarse, op_prolong, op_restrict) = op_mass_fine.create_multigrid_level_H1( 2039 /// &multiplicity, 2040 /// &ru_coarse, 2041 /// &bu_coarse, 2042 /// &interp_c_to_f, 2043 /// )?; 2044 /// 2045 /// // Coarse problem 2046 /// u_coarse.set_value(1.0); 2047 /// op_mass_coarse.apply(&u_coarse, &mut v_coarse)?; 2048 /// 2049 /// // Check 2050 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 2051 /// assert!( 2052 /// (sum - 2.0).abs() < 10.0 * libceed::EPSILON, 2053 /// "Incorrect interval length computed" 2054 /// ); 2055 /// 2056 /// // Prolong 2057 /// op_prolong.apply(&u_coarse, &mut u_fine)?; 2058 /// 2059 /// // Fine problem 2060 /// op_mass_fine.apply(&u_fine, &mut v_fine)?; 2061 /// 2062 /// // Check 2063 /// let sum: Scalar = v_fine.view()?.iter().sum(); 2064 /// assert!( 2065 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 2066 /// "Incorrect interval length computed" 2067 /// ); 2068 /// 2069 /// // Restrict 2070 /// op_restrict.apply(&v_fine, &mut v_coarse)?; 2071 /// 2072 /// // Check 2073 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 2074 /// assert!( 2075 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 2076 /// "Incorrect interval length computed" 2077 /// ); 2078 /// # Ok(()) 2079 /// # } 2080 /// ``` 2081 pub fn create_multigrid_level_H1<'b>( 2082 &self, 2083 p_mult_fine: &Vector, 2084 rstr_coarse: &ElemRestriction, 2085 basis_coarse: &Basis, 2086 interpCtoF: &[Scalar], 2087 ) -> crate::Result<(Operator<'b>, Operator<'b>, Operator<'b>)> { 2088 let mut ptr_coarse = std::ptr::null_mut(); 2089 let mut ptr_prolong = std::ptr::null_mut(); 2090 let mut ptr_restrict = std::ptr::null_mut(); 2091 self.op_core.check_error(unsafe { 2092 bind_ceed::CeedOperatorMultigridLevelCreateH1( 2093 self.op_core.ptr, 2094 p_mult_fine.ptr, 2095 rstr_coarse.ptr, 2096 basis_coarse.ptr, 2097 interpCtoF.as_ptr(), 2098 &mut ptr_coarse, 2099 &mut ptr_prolong, 2100 &mut ptr_restrict, 2101 ) 2102 })?; 2103 let op_coarse = Operator::from_raw(ptr_coarse)?; 2104 let op_prolong = Operator::from_raw(ptr_prolong)?; 2105 let op_restrict = Operator::from_raw(ptr_restrict)?; 2106 Ok((op_coarse, op_prolong, op_restrict)) 2107 } 2108 } 2109 2110 // ----------------------------------------------------------------------------- 2111 // Composite Operator 2112 // ----------------------------------------------------------------------------- 2113 impl<'a> CompositeOperator<'a> { 2114 // Constructor 2115 pub fn create(ceed: &crate::Ceed) -> crate::Result<Self> { 2116 let mut ptr = std::ptr::null_mut(); 2117 ceed.check_error(unsafe { bind_ceed::CeedCompositeOperatorCreate(ceed.ptr, &mut ptr) })?; 2118 Ok(Self { 2119 op_core: OperatorCore { 2120 ptr, 2121 _lifeline: PhantomData, 2122 }, 2123 }) 2124 } 2125 2126 /// Set name for CompositeOperator printing 2127 /// 2128 /// * 'name' - Name to set 2129 /// 2130 /// ``` 2131 /// # use libceed::prelude::*; 2132 /// # fn main() -> libceed::Result<()> { 2133 /// # let ceed = libceed::Ceed::default_init(); 2134 /// 2135 /// // Sub operator field arguments 2136 /// let ne = 3; 2137 /// let q = 4 as usize; 2138 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 2139 /// for i in 0..ne { 2140 /// ind[2 * i + 0] = i as i32; 2141 /// ind[2 * i + 1] = (i + 1) as i32; 2142 /// } 2143 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 2144 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2145 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2146 /// 2147 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2148 /// 2149 /// let qdata_mass = ceed.vector(q * ne)?; 2150 /// let qdata_diff = ceed.vector(q * ne)?; 2151 /// 2152 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2153 /// let op_mass = ceed 2154 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2155 /// .name("Mass term")? 2156 /// .field("u", &r, &b, VectorOpt::Active)? 2157 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 2158 /// .field("v", &r, &b, VectorOpt::Active)?; 2159 /// 2160 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2161 /// let op_diff = ceed 2162 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 2163 /// .name("Poisson term")? 2164 /// .field("du", &r, &b, VectorOpt::Active)? 2165 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 2166 /// .field("dv", &r, &b, VectorOpt::Active)?; 2167 /// 2168 /// let op = ceed 2169 /// .composite_operator()? 2170 /// .name("Screened Poisson")? 2171 /// .sub_operator(&op_mass)? 2172 /// .sub_operator(&op_diff)?; 2173 /// # Ok(()) 2174 /// # } 2175 /// ``` 2176 #[allow(unused_mut)] 2177 pub fn name(mut self, name: &str) -> crate::Result<Self> { 2178 self.op_core.name(name)?; 2179 Ok(self) 2180 } 2181 2182 /// Apply Operator to a vector 2183 /// 2184 /// * `input` - Inpuht Vector 2185 /// * `output` - Output Vector 2186 /// 2187 /// ``` 2188 /// # use libceed::prelude::*; 2189 /// # fn main() -> libceed::Result<()> { 2190 /// # let ceed = libceed::Ceed::default_init(); 2191 /// let ne = 4; 2192 /// let p = 3; 2193 /// let q = 4; 2194 /// let ndofs = p * ne - ne + 1; 2195 /// 2196 /// // Vectors 2197 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 2198 /// let mut qdata_mass = ceed.vector(ne * q)?; 2199 /// qdata_mass.set_value(0.0); 2200 /// let mut qdata_diff = ceed.vector(ne * q)?; 2201 /// qdata_diff.set_value(0.0); 2202 /// let mut u = ceed.vector(ndofs)?; 2203 /// u.set_value(1.0); 2204 /// let mut v = ceed.vector(ndofs)?; 2205 /// v.set_value(0.0); 2206 /// 2207 /// // Restrictions 2208 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 2209 /// for i in 0..ne { 2210 /// indx[2 * i + 0] = i as i32; 2211 /// indx[2 * i + 1] = (i + 1) as i32; 2212 /// } 2213 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 2214 /// let mut indu: Vec<i32> = vec![0; p * ne]; 2215 /// for i in 0..ne { 2216 /// indu[p * i + 0] = i as i32; 2217 /// indu[p * i + 1] = (i + 1) as i32; 2218 /// indu[p * i + 2] = (i + 2) as i32; 2219 /// } 2220 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 2221 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2222 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2223 /// 2224 /// // Bases 2225 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2226 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 2227 /// 2228 /// // Build quadrature data 2229 /// let qf_build_mass = ceed.q_function_interior_by_name("Mass1DBuild")?; 2230 /// ceed.operator(&qf_build_mass, QFunctionOpt::None, QFunctionOpt::None)? 2231 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2232 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2233 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2234 /// .apply(&x, &mut qdata_mass)?; 2235 /// 2236 /// let qf_build_diff = ceed.q_function_interior_by_name("Poisson1DBuild")?; 2237 /// ceed.operator(&qf_build_diff, QFunctionOpt::None, QFunctionOpt::None)? 2238 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2239 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2240 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2241 /// .apply(&x, &mut qdata_diff)?; 2242 /// 2243 /// // Application operator 2244 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2245 /// let op_mass = ceed 2246 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2247 /// .field("u", &ru, &bu, VectorOpt::Active)? 2248 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 2249 /// .field("v", &ru, &bu, VectorOpt::Active)?; 2250 /// 2251 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2252 /// let op_diff = ceed 2253 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 2254 /// .field("du", &ru, &bu, VectorOpt::Active)? 2255 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 2256 /// .field("dv", &ru, &bu, VectorOpt::Active)?; 2257 /// 2258 /// let op_composite = ceed 2259 /// .composite_operator()? 2260 /// .sub_operator(&op_mass)? 2261 /// .sub_operator(&op_diff)?; 2262 /// 2263 /// v.set_value(0.0); 2264 /// op_composite.apply(&u, &mut v)?; 2265 /// 2266 /// // Check 2267 /// let sum: Scalar = v.view()?.iter().sum(); 2268 /// assert!( 2269 /// (sum - 2.0).abs() < 10.0 * libceed::EPSILON, 2270 /// "Incorrect interval length computed" 2271 /// ); 2272 /// # Ok(()) 2273 /// # } 2274 /// ``` 2275 pub fn apply(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 2276 self.op_core.apply(input, output) 2277 } 2278 2279 /// Apply Operator to a vector and add result to output vector 2280 /// 2281 /// * `input` - Input Vector 2282 /// * `output` - Output Vector 2283 /// 2284 /// ``` 2285 /// # use libceed::prelude::*; 2286 /// # fn main() -> libceed::Result<()> { 2287 /// # let ceed = libceed::Ceed::default_init(); 2288 /// let ne = 4; 2289 /// let p = 3; 2290 /// let q = 4; 2291 /// let ndofs = p * ne - ne + 1; 2292 /// 2293 /// // Vectors 2294 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 2295 /// let mut qdata_mass = ceed.vector(ne * q)?; 2296 /// qdata_mass.set_value(0.0); 2297 /// let mut qdata_diff = ceed.vector(ne * q)?; 2298 /// qdata_diff.set_value(0.0); 2299 /// let mut u = ceed.vector(ndofs)?; 2300 /// u.set_value(1.0); 2301 /// let mut v = ceed.vector(ndofs)?; 2302 /// v.set_value(0.0); 2303 /// 2304 /// // Restrictions 2305 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 2306 /// for i in 0..ne { 2307 /// indx[2 * i + 0] = i as i32; 2308 /// indx[2 * i + 1] = (i + 1) as i32; 2309 /// } 2310 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 2311 /// let mut indu: Vec<i32> = vec![0; p * ne]; 2312 /// for i in 0..ne { 2313 /// indu[p * i + 0] = i as i32; 2314 /// indu[p * i + 1] = (i + 1) as i32; 2315 /// indu[p * i + 2] = (i + 2) as i32; 2316 /// } 2317 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 2318 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2319 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2320 /// 2321 /// // Bases 2322 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2323 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 2324 /// 2325 /// // Build quadrature data 2326 /// let qf_build_mass = ceed.q_function_interior_by_name("Mass1DBuild")?; 2327 /// ceed.operator(&qf_build_mass, QFunctionOpt::None, QFunctionOpt::None)? 2328 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2329 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2330 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2331 /// .apply(&x, &mut qdata_mass)?; 2332 /// 2333 /// let qf_build_diff = ceed.q_function_interior_by_name("Poisson1DBuild")?; 2334 /// ceed.operator(&qf_build_diff, QFunctionOpt::None, QFunctionOpt::None)? 2335 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2336 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2337 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2338 /// .apply(&x, &mut qdata_diff)?; 2339 /// 2340 /// // Application operator 2341 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2342 /// let op_mass = ceed 2343 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2344 /// .field("u", &ru, &bu, VectorOpt::Active)? 2345 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 2346 /// .field("v", &ru, &bu, VectorOpt::Active)?; 2347 /// 2348 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2349 /// let op_diff = ceed 2350 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 2351 /// .field("du", &ru, &bu, VectorOpt::Active)? 2352 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 2353 /// .field("dv", &ru, &bu, VectorOpt::Active)?; 2354 /// 2355 /// let op_composite = ceed 2356 /// .composite_operator()? 2357 /// .sub_operator(&op_mass)? 2358 /// .sub_operator(&op_diff)?; 2359 /// 2360 /// v.set_value(1.0); 2361 /// op_composite.apply_add(&u, &mut v)?; 2362 /// 2363 /// // Check 2364 /// let sum: Scalar = v.view()?.iter().sum(); 2365 /// assert!( 2366 /// (sum - (2.0 + ndofs as Scalar)).abs() < 10.0 * libceed::EPSILON, 2367 /// "Incorrect interval length computed" 2368 /// ); 2369 /// # Ok(()) 2370 /// # } 2371 /// ``` 2372 pub fn apply_add(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 2373 self.op_core.apply_add(input, output) 2374 } 2375 2376 /// Add a sub-Operator to a Composite Operator 2377 /// 2378 /// * `subop` - Sub-Operator 2379 /// 2380 /// ``` 2381 /// # use libceed::prelude::*; 2382 /// # fn main() -> libceed::Result<()> { 2383 /// # let ceed = libceed::Ceed::default_init(); 2384 /// let mut op = ceed.composite_operator()?; 2385 /// 2386 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2387 /// let op_mass = ceed.operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)?; 2388 /// op = op.sub_operator(&op_mass)?; 2389 /// 2390 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2391 /// let op_diff = ceed.operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)?; 2392 /// op = op.sub_operator(&op_diff)?; 2393 /// # Ok(()) 2394 /// # } 2395 /// ``` 2396 #[allow(unused_mut)] 2397 pub fn sub_operator(mut self, subop: &Operator) -> crate::Result<Self> { 2398 self.op_core.check_error(unsafe { 2399 bind_ceed::CeedCompositeOperatorAddSub(self.op_core.ptr, subop.op_core.ptr) 2400 })?; 2401 Ok(self) 2402 } 2403 2404 /// Check if CompositeOperator is setup correctly 2405 /// 2406 /// ``` 2407 /// # use libceed::prelude::*; 2408 /// # fn main() -> libceed::Result<()> { 2409 /// # let ceed = libceed::Ceed::default_init(); 2410 /// let ne = 4; 2411 /// let p = 3; 2412 /// let q = 4; 2413 /// let ndofs = p * ne - ne + 1; 2414 /// 2415 /// // Restrictions 2416 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 2417 /// for i in 0..ne { 2418 /// indx[2 * i + 0] = i as i32; 2419 /// indx[2 * i + 1] = (i + 1) as i32; 2420 /// } 2421 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 2422 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2423 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2424 /// 2425 /// // Bases 2426 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2427 /// 2428 /// // Build quadrature data 2429 /// let qf_build_mass = ceed.q_function_interior_by_name("Mass1DBuild")?; 2430 /// let op_build_mass = ceed 2431 /// .operator(&qf_build_mass, QFunctionOpt::None, QFunctionOpt::None)? 2432 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2433 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2434 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 2435 /// 2436 /// let qf_build_diff = ceed.q_function_interior_by_name("Poisson1DBuild")?; 2437 /// let op_build_diff = ceed 2438 /// .operator(&qf_build_diff, QFunctionOpt::None, QFunctionOpt::None)? 2439 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2440 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2441 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 2442 /// 2443 /// let op_build = ceed 2444 /// .composite_operator()? 2445 /// .sub_operator(&op_build_mass)? 2446 /// .sub_operator(&op_build_diff)?; 2447 /// 2448 /// // Check 2449 /// assert!(op_build.check().is_ok(), "Operator setup incorrectly"); 2450 /// # Ok(()) 2451 /// # } 2452 /// ``` 2453 pub fn check(self) -> crate::Result<Self> { 2454 self.op_core.check()?; 2455 Ok(self) 2456 } 2457 2458 /// Assemble the diagonal of a square linear Operator 2459 /// 2460 /// This overwrites a Vector with the diagonal of a linear Operator. 2461 /// 2462 /// Note: Currently only non-composite Operators with a single field and 2463 /// composite Operators with single field sub-operators are supported. 2464 /// 2465 /// * `op` - Operator to assemble QFunction 2466 /// * `assembled` - Vector to store assembled Operator diagonal 2467 pub fn linear_assemble_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 2468 self.op_core.linear_assemble_add_diagonal(assembled) 2469 } 2470 2471 /// Assemble the point block diagonal of a square linear Operator 2472 /// 2473 /// This overwrites a Vector with the point block diagonal of a linear 2474 /// Operator. 2475 /// 2476 /// Note: Currently only non-composite Operators with a single field and 2477 /// composite Operators with single field sub-operators are supported. 2478 /// 2479 /// * `op` - Operator to assemble QFunction 2480 /// * `assembled` - Vector to store assembled Operator diagonal 2481 pub fn linear_assemble_add_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 2482 self.op_core.linear_assemble_add_diagonal(assembled) 2483 } 2484 2485 /// Assemble the diagonal of a square linear Operator 2486 /// 2487 /// This overwrites a Vector with the diagonal of a linear Operator. 2488 /// 2489 /// Note: Currently only non-composite Operators with a single field and 2490 /// composite Operators with single field sub-operators are supported. 2491 /// 2492 /// * `op` - Operator to assemble QFunction 2493 /// * `assembled` - Vector to store assembled CeedOperator point block 2494 /// diagonal, provided in row-major form with an 2495 /// `ncomp * ncomp` block at each node. The dimensions of 2496 /// this vector are derived from the active vector for 2497 /// the CeedOperator. The array has shape 2498 /// `[nodes, component out, component in]`. 2499 pub fn linear_assemble_point_block_diagonal( 2500 &self, 2501 assembled: &mut Vector, 2502 ) -> crate::Result<i32> { 2503 self.op_core.linear_assemble_add_diagonal(assembled) 2504 } 2505 2506 /// Assemble the diagonal of a square linear Operator 2507 /// 2508 /// This sums into a Vector with the diagonal of a linear Operator. 2509 /// 2510 /// Note: Currently only non-composite Operators with a single field and 2511 /// composite Operators with single field sub-operators are supported. 2512 /// 2513 /// * `op` - Operator to assemble QFunction 2514 /// * `assembled` - Vector to store assembled CeedOperator point block 2515 /// diagonal, provided in row-major form with an 2516 /// `ncomp * ncomp` block at each node. The dimensions of 2517 /// this vector are derived from the active vector for 2518 /// the CeedOperator. The array has shape 2519 /// `[nodes, component out, component in]`. 2520 pub fn linear_assemble_add_point_block_diagonal( 2521 &self, 2522 assembled: &mut Vector, 2523 ) -> crate::Result<i32> { 2524 self.op_core.linear_assemble_add_diagonal(assembled) 2525 } 2526 } 2527 2528 // ----------------------------------------------------------------------------- 2529