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