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