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