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