Mercurial > repos > pointsource_pde / file comparison

comparison: src/dolfinx_access/function.rs

src/dolfinx_access/function.rs

changeset 1
a4137aedcb3a
child 3
c3a4f4bb87f7
equal deleted inserted replaced
0:7ec1cfe19a24 1:a4137aedcb3a
1 /*!
2 Python and C++ wrapper for Dolfinx Function<f64>
3 */
4
5 use super::ffi;
6 use super::ffi::Function_f64;
7 use alg_tools::bounds::{Bounds, GlobalAnalysis, MinMaxMapping};
8 use alg_tools::euclidean::Euclidean;
9 use alg_tools::fe_model::base::LocalModel;
10 use alg_tools::instance::Instance;
11 use alg_tools::loc::Loc;
12 use alg_tools::mapping::{BasicDecomposition, DifferentiableImpl, Mapping, Space};
13 use alg_tools::norms::{Dist, HasDual, Norm, Normed, L2};
14 use itertools::izip;
15 use pyo3::ffi::c_str;
16 use pyo3::prelude::*;
17 use pyo3::types::{PyDict, PyFunction};
18 use std::ffi::CStr;
19
20 /// A helper structure of dealing with dolfinx functions.
21 /// `N` is the domain dimension, `O` the order, and `D` is the codomain dimension.
22 #[allow(non_camel_case_types)]
23 pub struct DolfinxPyFunction_f64<'py, const N: u32, const O: u32 = 1, const D: u32 = 1> {
24 #[allow(dead_code)]
25 /// Python object.
26 ///
27 /// This is crucial for maintaining [`cxx`] live, when we obtained the object from Python.
28 ///
29 /// If we needed to copy-on-write, or initialise without obtaining and object from Python,
30 /// this will become `None`.
31 pub u: Option<Bound<'py, PyAny>>,
32 /// We need to maintain a reference to the FunctionSpace Python wrapper, as it's difficult
33 /// to recreate if we loose `py`. Fenics is a one-way hell full of pointless wrappers upon
34 /// wrappers upon wrappers,very difficult to create having with just the unwrapped object.
35 function_space: Bound<'py, PyAny>,
36 /// Do we have the only reference?
37 ///
38 /// This is check on construction only, when we have the GIL.
39 /// If we have the only reference at that time, this will not change during the lifetime.
40 /// We use this to implement copy-on-write optimisations.
41 owned: bool,
42 /// C++ object pointer
43 cxx: *mut Function_f64,
44 }
45
46 struct Helpers {
47 dot: Py<PyFunction>,
48 norm2_squared: Py<PyFunction>,
49 dist2_squared: Py<PyFunction>,
50 }
51
52 use std::sync::OnceLock;
53 static HELPERS: OnceLock<Helpers> = OnceLock::new();
54
55 fn get_helpers(py: Python<'_>) -> PyResult<&Helpers> {
56 HELPERS.get_or_try_init(|| {
57 let extras = PyModule::import(py, "pointsource_pde.dolfinx_extras")?;
58 let g = |s| -> PyResult<Py<PyFunction>> { Ok(extras.getattr(s)?.cast_into()?.unbind()) };
59 Ok(Helpers {
60 dot: g("dot")?,
61 norm2_squared: g("norm2_squared")?,
62 dist2_squared: g("dist2_squared")?,
63 })
64 })
65 }
66
67 impl<'py, const N: u32, const O: u32, const D: u32> Drop for DolfinxPyFunction_f64<'py, N, O, D> {
68 fn drop(&mut self) {
69 // If Python is not managing our C++ object, we need to free it
70 if let None = self.u {
71 if self.cxx != std::ptr::null_mut() {
72 unsafe { ffi::drop_Function_f64(self.cxx) };
73 }
74 }
75 }
76 }
77
78 impl<'py, const N: u32, const O: u32, const D: u32> Clone for DolfinxPyFunction_f64<'py, N, O, D> {
79 fn clone(&self) -> Self {
80 DolfinxPyFunction_f64 {
81 u: None, // It's a different new function, so no Python instance
82 function_space: self.function_space.clone(),
83 owned: self.owned,
84 cxx: unsafe { ffi::clone_Function_f64(self.cxx) },
85 }
86 }
87 }
88
89 impl<'py, const N: u32, const O: u32, const D: u32> Space for DolfinxPyFunction_f64<'py, N, O, D> {
90 type Principal = Self;
91 type Decomp = BasicDecomposition;
92 }
93
94 impl<'py, const N: u32, const O: u32, const D: u32> DolfinxPyFunction_f64<'py, N, O, D> {
95 /// Create a new instance, assuming `cxx` has been correctly extracted from `u`
96 /// and checked to match `N`, `O`, `D`, and mesh restrictions.
97 pub(crate) fn new_prechecked(u: Bound<'py, PyAny>, cxx: *mut Function_f64) -> PyResult<Self> {
98 let function_space = u.getattr("_V")?;
99 let owned = u.get_refcnt() <= 1;
100 // TODO: check dimensionality
101 Ok(DolfinxPyFunction_f64 { u: Some(u), function_space, owned, cxx })
102 }
103
104 /// Ensures that `self.u` exists.
105 /// Panics if there's a failure to create it if needed.
106 pub(crate) fn _make_py(&mut self, py: Python<'py>) {
107 if let None = self.u {
108 // We need to do a lot of work here, due to Fenics one-wayness.
109
110 // This takes ownership of our pointer, so if we fail to create the Python
111 // object, we have a lot of restoring to do.
112 let cpp_object = unsafe {
113 Py::<PyAny>::from_owned_ptr(
114 py,
115 ffi::wrap_Function_f64(self.cxx) as *mut pyo3::ffi::PyObject,
116 )
117 };
118 let locals = PyDict::new(py);
119 match locals.set_item("cpp_object", cpp_object).and_then(|()| {
120 locals.set_item("V", self.function_space.as_borrowed())?;
121 // TODO: precompile
122 py.run(CREATE_FUNCTION_HACK, None, Some(&locals))?;
123 locals.get_item("wrapper")
124 }) {
125 Err(e) => {
126 // We should restore to recover
127 panic!("Failure to pass object to Python: {:?}", e);
128 }
129 Ok(maybe_p) => {
130 self.u = Some(maybe_p.unwrap());
131 }
132 }
133 }
134 }
135
136 /// Returns a (possibly new) Python presentation, without consuming current object.
137 ///
138 /// Please use this functionality through [`pyo3::types::IntoPyObject`].
139 pub(crate) fn _into_py_ref<'a>(
140 &'a mut self,
141 py: Python<'py>,
142 ) -> PyResult<pyo3::Borrowed<'a, 'py, PyAny>> {
143 self._make_py(py);
144 self.owned = false;
145 match self.u {
146 Some(ref p) => Ok(p.as_borrowed()),
147 None => unreachable!("This cannot happen"),
148 }
149 }
150
151 /// Returns a (possibly new) Python presentation, consuming current object.
152 ///
153 /// Please use this functionality through [`pyo3::types::IntoPyObject`].
154 pub(crate) fn _into_py(mut self, py: Python<'py>) -> PyResult<pyo3::Bound<'py, PyAny>> {
155 self._make_py(py);
156 self.cxx = std::ptr::null_mut();
157 self.owned = false;
158 match std::mem::replace(&mut self.u, None) {
159 Some(p) => Ok(p),
160 None => unreachable!("This cannot happen"),
161 }
162 }
163
164 /// Create a new FEM function on the same function space, with unintialised data.
165 pub fn similar(&self) -> Self {
166 DolfinxPyFunction_f64 {
167 u: None,
168 function_space: self.function_space.clone(),
169 owned: true,
170 cxx: unsafe { ffi::similar_Function_f64(self.cxx) },
171 }
172 }
173
174 /// Create a new FEM function on the same function space, with zero data.
175 pub fn similar_zeros(&self) -> Self {
176 let new = self.similar();
177 for x in unsafe { ffi::data_mut_Function_f64(new.cxx) } {
178 *x = 0.0;
179 }
180 new
181 }
182 }
183
184 pub(crate) static CREATE_FUNCTION_HACK: &CStr = c_str!(
185 "\
186 from dolfinx.fem import Function
187 from dolfinx.la import Vector
188 x = Vector(cpp_object.x)
189 # Fenics is very one-way. It does not provide a proper way to create the wrapper,
190 # directly from the unwrapped C++
191 wrapper = Function.__new__(Function, V, x, \"?\", x.array.dtype)
192 # This all just repeats Function.__init__ without creating a new ._cpp_object.
193 super(Function, wrapper).__init__(V.ufl_function_space())
194 wrapper._cpp_object = cpp_object;
195 wrapper._V = V
196 wrapper._x = x
197 wrapper.name = \"?\"
198 "
199 );
200
201 trait ExpandCoordDolphinx {
202 fn to_dolphinx(self) -> [f64; 3];
203 #[allow(dead_code)]
204 fn from_dolphinx(val: [f64; 3]) -> Self;
205 }
206
207 impl ExpandCoordDolphinx for Loc<1, f64> {
208 #[inline]
209 fn to_dolphinx(self) -> [f64; 3] {
210 let Loc([x1]) = self;
211 [x1, 0.0, 0.0]
212 }
213
214 #[inline]
215 fn from_dolphinx([x1, _, _]: [f64; 3]) -> Self {
216 [x1].into()
217 }
218 }
219
220 impl ExpandCoordDolphinx for Loc<2, f64> {
221 #[inline]
222 fn to_dolphinx(self) -> [f64; 3] {
223 let Loc([x1, x2]) = self;
224 [x1, x2, 0.0]
225 }
226
227 #[inline]
228 fn from_dolphinx([x1, x2, _]: [f64; 3]) -> Self {
229 [x1, x2].into()
230 }
231 }
232
233 impl ExpandCoordDolphinx for Loc<3, f64> {
234 #[inline]
235 fn to_dolphinx(self) -> [f64; 3] {
236 let Loc([x1, x2, x3]) = self;
237 [x1, x2, x3]
238 }
239
240 #[inline]
241 fn from_dolphinx([x1, x2, x3]: [f64; 3]) -> Self {
242 [x1, x2, x3].into()
243 }
244 }
245
246 macro_rules! impl_diff {
247 ($n:literal ; $($o:literal)+) => { $(
248 impl<'py> DifferentiableImpl<Loc<$n, f64>> for DolfinxPyFunction_f64<'py, $n, $o, 1> {
249 type Derivative = Loc<$n, f64>;
250
251 fn differential_impl<I: Instance<Loc<$n, f64>>>(&self, x: I) -> Self::Derivative {
252 let x_own = x.own();
253 let f = self.cxx;
254 let expand_x = x_own.to_dolphinx();
255 // Find cell containing `x`
256 let cell = unsafe { ffi::cell_Function_f64(f, &expand_x) };
257 if cell < 0 {
258 Loc::ORIGIN
259 } else {
260 let (model, _) = unsafe {
261 // Find the coordinates of the corresponding triangle. We guarantee a
262 // triangular mesh in the constructor of `DolfinxPyFunction_f64`.
263 let simplex_coords = ffi::cell_coords_Function_f64_triangle(f, cell);
264 // Construct a new P2 model on that cell, as this is easier than dealing
265 // with the awful documentation of Dolphinx and C++, and probably having
266 // to reimplemen our P2 stuff anyway, to do any optimisation with the basix
267 // basis functions.
268 super::model_dolfinx_p2_cell(f, &simplex_coords, cell, false)
269 };
270 model.differential(&x_own)
271 }
272 }
273 }
274 )+ };
275 }
276
277 macro_rules! impl_mapping {
278 ($($n:literal)+) => { $(
279 impl<'py, const O: u32> Mapping<Loc<$n, f64>> for DolfinxPyFunction_f64<'py, $n, O, 1> {
280 type Codomain = f64;
281
282 fn apply<I: Instance<Loc<$n, f64>>>(&self, x: I) -> f64 {
283 let expand_x = x.own().to_dolphinx();
284 let f = self.cxx;
285 unsafe {
286 // Find cell containing `x`
287 return ffi::eval_Function_f64(f, &expand_x);
288 }
289 }
290 }
291 //impl_diff!($n; 2 3 4 5 6 7 8);
292 )+ };
293 }
294
295 impl_mapping!(1 2 3);
296 // Extension to 1D should be easy.
297 //impl_diff!(1; 2 3 4 5 6 7 8);
298 impl_diff!(2; 2 3 4 5 6 7 8);
299
300 impl<'py> GlobalAnalysis<f64, Bounds<f64>> for DolfinxPyFunction_f64<'py, 2, 2, 1> {
301 fn global_analysis(&self) -> Bounds<f64> {
302 let min =
303 unsafe { ffi::min_Function_f64_p2(self.cxx) }.map_or(f64::NEG_INFINITY, |res| res.v);
304 let max = unsafe { ffi::max_Function_f64_p2(self.cxx) }.map_or(f64::INFINITY, |res| res.v);
305 Bounds(min, max)
306 }
307 }
308
309 /// Only second-order polynomials on a triangular mesh are implemented.
310 ///
311 /// May panic if assumptions of the underlying functionn do not hold despite
312 /// verification in the constructors of [`DolfinxPyFunction_f64`].
313 impl<'py> MinMaxMapping<Loc<2>, f64> for DolfinxPyFunction_f64<'py, 2, 2, 1> {
314 fn minimise(&mut self, _tolerance: f64, _max_steps: usize) -> (Loc<2>, f64) {
315 let res = unsafe { ffi::min_Function_f64_p2(self.cxx) }.unwrap();
316 (res.x.into(), res.v)
317 }
318
319 fn maximise(&mut self, _tolerance: f64, _max_steps: usize) -> (Loc<2>, f64) {
320 let res = unsafe { ffi::max_Function_f64_p2(self.cxx) }.unwrap();
321 (res.x.into(), res.v)
322 }
323 }
324
325 /*
326 enum MaybeValid<A> {
327 Valid(A),
328 Invalid(A),
329 }
330
331 use MaybeValid::*;
332
333 impl<A> MaybeValid<A> {
334 #[inline]
335 fn uninitialised_ok(self) -> A {
336 match self {
337 Valid(a) => a,
338 Invalid(a) => a,
339 }
340 }
341 }
342 */
343
344 impl<'py, const N: u32, const O: u32, const D: u32> DolfinxPyFunction_f64<'py, N, O, D> {
345 /// If owned, call `f_valid` with mutable data.
346 /// If not owned, create a new C++ fem.Function, uninitialised, and call `f_invalid` with its
347 /// mutable data, as well as the original immutable data.
348 #[inline]
349 fn with_mut_data_and_orig<F, G, Z>(&mut self, f_valid: F, f_invalid: G) -> Z
350 where
351 F: Fn(&mut [f64]) -> Z,
352 G: Fn(&mut [f64], &[f64]) -> Z,
353 {
354 if self.owned {
355 return f_valid(unsafe { ffi::data_mut_Function_f64(self.cxx) });
356 } else {
357 // Create a completely new Function
358 let old_data = unsafe { ffi::data_mut_Function_f64(self.cxx) };
359 let new = unsafe { ffi::similar_Function_f64(self.cxx) };
360 let new_data = unsafe { ffi::data_mut_Function_f64(new) };
361 let res = f_invalid(new_data, old_data);
362 // Invalidate Python reference; we're managing the object now
363 self.cxx = new;
364 self.u = None;
365 self.owned = true;
366 return res;
367 }
368 }
369
370 /// If owned, call `f` with mutable data.
371 /// If not owned, create a new C++ fem.Function, uninitialised, and call `f` with its
372 /// mutable data.
373 #[inline]
374 fn with_mut_data<F, Z>(&mut self, f: F) -> Z
375 where
376 F: Fn(&mut [f64]) -> Z,
377 {
378 if self.owned {
379 return f(unsafe { ffi::data_mut_Function_f64(self.cxx) });
380 } else {
381 // Create a completely new Function
382 let new = unsafe { ffi::similar_Function_f64(self.cxx) };
383 let new_data = unsafe { ffi::data_mut_Function_f64(new) };
384 let res = f(new_data);
385 // Invalidate Python reference; we're managing the object now
386 self.cxx = new;
387 self.u = None;
388 self.owned = true;
389 return res;
390 }
391 }
392
393 /// Create a similar new function, and call `f` with both the mutable (unininitialised)
394 /// new data, and original data.
395 ///
396 /// Returns object for the new data.
397 #[inline]
398 fn with_new_data_and_orig<F>(&self, f: F) -> Self
399 where
400 F: Fn(&mut [f64], &[f64]),
401 {
402 // Create a completely new Function
403 let new = self.similar();
404 let old_data = unsafe { ffi::data_Function_f64(self.cxx) };
405 let new_data = unsafe { ffi::data_mut_Function_f64(new.cxx) };
406 f(new_data, old_data);
407 new
408 }
409
410 /// Similar to [`with_mut_data`], but check compatibility with `other` and also pass its
411 /// immutable data slice as the final argument.
412 #[inline]
413 fn with_mut_data_and_orig_binop<F, G, Z>(
414 &mut self,
415 other: &DolfinxPyFunction_f64<'py, N, O, D>,
416 f_valid: F,
417 f_invalid: G,
418 ) -> Z
419 where
420 F: Fn(&mut [f64], &[f64]) -> Z,
421 G: Fn(&mut [f64], &[f64], &[f64]) -> Z,
422 {
423 let data_other = unsafe {
424 ffi::check_compat_Function_f64(self.cxx, other.cxx).unwrap();
425 ffi::data_Function_f64(other.cxx)
426 };
427 self.with_mut_data_and_orig(
428 |data_self_valid| {
429 assert_eq!(data_self_valid.len(), data_other.len());
430 f_valid(data_self_valid, data_other)
431 },
432 |data_self_invalid, data_orig| {
433 assert_eq!(data_self_invalid.len(), data_other.len());
434 assert_eq!(data_orig.len(), data_other.len());
435 f_invalid(data_self_invalid, data_orig, data_other)
436 },
437 )
438 }
439
440 /// Similar to [`with_new_data_and_orig`], but check compatibility with `other` and also
441 /// pass its immutable data slice as the final argument.
442 #[inline]
443 fn with_new_data_and_orig_binop<F>(
444 &self,
445 other: &DolfinxPyFunction_f64<'py, N, O, D>,
446 f: F,
447 ) -> Self
448 where
449 F: Fn(&mut [f64], &[f64], &[f64]),
450 {
451 let data_other = unsafe {
452 ffi::check_compat_Function_f64(self.cxx, other.cxx).unwrap();
453 ffi::data_Function_f64(other.cxx)
454 };
455 self.with_new_data_and_orig(|data, data_orig| {
456 assert_eq!(data.len(), data_other.len());
457 assert_eq!(data_orig.len(), data_other.len());
458 f(data, data_orig, data_other)
459 })
460 }
461 }
462
463 macro_rules! impl_unop {
464 ($trait:ident, $fn:ident) => {
465 impl<'py, const N: u32, const O: u32, const D: u32> $trait
466 for DolfinxPyFunction_f64<'py, N, O, D>
467 {
468 type Output = Self;
469
470 fn $fn(mut self) -> Self {
471 self.with_mut_data_and_orig(
472 |a| a.iter_mut().for_each(|x| *x = x.$fn()),
473 |a, a_orig| {
474 izip!(a.iter_mut(), a_orig.iter()).for_each(|(x, x_orig)| {
475 *x = x_orig.$fn();
476 })
477 },
478 );
479 return self;
480 }
481 }
482
483 impl<'py, 'a, const N: u32, const O: u32, const D: u32> $trait
484 for &'a DolfinxPyFunction_f64<'py, N, O, D>
485 {
486 type Output = DolfinxPyFunction_f64<'py, N, O, D>;
487
488 fn $fn(self) -> Self::Output {
489 self.with_new_data_and_orig(|a, a_orig| {
490 izip!(a.iter_mut(), a_orig.iter()).for_each(|(x, x_orig)| {
491 *x = x_orig.$fn();
492 })
493 })
494 }
495 }
496 };
497 }
498
499 use std::ops::Neg;
500
501 impl_unop!(Neg, neg);
502
503 macro_rules! impl_scalarop {
504 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
505 impl<'py, const N: u32, const O: u32, const D: u32> $trait<f64>
506 for DolfinxPyFunction_f64<'py, N, O, D>
507 {
508 type Output = Self;
509
510 fn $fn(mut self, α: f64) -> Self {
511 self.with_mut_data_and_orig(
512 |a| {
513 a.iter_mut().for_each(|x| {
514 x.$fn_assign(α);
515 })
516 },
517 |a, a_orig| {
518 izip!(a.iter_mut(), a_orig.iter()).for_each(|(x, x_orig)| {
519 *x = x_orig.$fn(α);
520 })
521 },
522 );
523 return self;
524 }
525 }
526
527 impl<'py, 'a, const N: u32, const O: u32, const D: u32> $trait<f64>
528 for &'a DolfinxPyFunction_f64<'py, N, O, D>
529 {
530 type Output = DolfinxPyFunction_f64<'py, N, O, D>;
531
532 fn $fn(self, α: f64) -> Self::Output {
533 self.with_new_data_and_orig(|a, a_orig| {
534 izip!(a.iter_mut(), a_orig.iter()).for_each(|(x, x_orig)| {
535 *x = x_orig.$fn(α);
536 })
537 })
538 }
539 }
540
541 impl<'py, const N: u32, const O: u32, const D: u32> $trait_assign<f64>
542 for DolfinxPyFunction_f64<'py, N, O, D>
543 {
544 fn $fn_assign(&mut self, α: f64) {
545 self.with_mut_data_and_orig(
546 |a| {
547 a.iter_mut().for_each(|x| {
548 x.$fn_assign(α);
549 })
550 },
551 |a, a_orig| {
552 izip!(a.iter_mut(), a_orig.iter()).for_each(|(x, x_orig)| {
553 *x = x_orig.$fn(α);
554 })
555 },
556 );
557 }
558 }
559 };
560 }
561
562 use std::ops::{Div, DivAssign, Mul, MulAssign};
563
564 impl_scalarop!(Mul, mul, MulAssign, mul_assign);
565 impl_scalarop!(Div, div, DivAssign, div_assign);
566
567 macro_rules! impl_binop_consume {
568 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
569 impl<'py, const N: u32, const O: u32, const D: u32> $trait<Self>
570 for DolfinxPyFunction_f64<'py, N, O, D>
571 {
572 type Output = Self;
573
574 #[inline]
575 fn $fn(self, other: Self) -> Self {
576 self.$fn(&other)
577 }
578 }
579
580 impl<'py, const N: u32, const O: u32, const D: u32> $trait_assign<Self>
581 for DolfinxPyFunction_f64<'py, N, O, D>
582 {
583 #[inline]
584 fn $fn_assign(&mut self, other: Self) {
585 self.$fn_assign(&other)
586 }
587 }
588
589 impl<'py, 'a, const N: u32, const O: u32, const D: u32> $trait<&'a Self>
590 for DolfinxPyFunction_f64<'py, N, O, D>
591 {
592 type Output = Self;
593
594 fn $fn(mut self, other: &'a Self) -> Self {
595 self.with_mut_data_and_orig_binop(
596 other,
597 |data, data_other| {
598 izip!(data.iter_mut(), data_other.iter()).for_each(|(x, y)| {
599 x.$fn_assign(y);
600 })
601 },
602 |data, data_orig, data_other| {
603 izip!(data.iter_mut(), data_orig.iter(), data_other.iter()).for_each(
604 |(x, x_orig, y)| {
605 *x = x_orig.$fn(y);
606 },
607 )
608 },
609 );
610 return self;
611 }
612 }
613
614 impl<'py, 'a, const N: u32, const O: u32, const D: u32> $trait_assign<&'a Self>
615 for DolfinxPyFunction_f64<'py, N, O, D>
616 {
617 fn $fn_assign(&mut self, other: &'a Self) {
618 self.with_mut_data_and_orig_binop(
619 other,
620 |data, data_other| {
621 izip!(data.iter_mut(), data_other.iter()).for_each(|(x, y)| {
622 x.$fn_assign(y);
623 })
624 },
625 |data, data_orig, data_other| {
626 izip!(data.iter_mut(), data_orig.iter(), data_other.iter()).for_each(
627 |(x, x_orig, y)| {
628 *x = x_orig.$fn(y);
629 },
630 )
631 },
632 );
633 }
634 }
635 };
636 }
637
638 macro_rules! impl_binop_ref {
639 ($trait:ident, $fn:ident) => {
640 impl<'py, 'b, const N: u32, const O: u32, const D: u32> $trait<Self>
641 for &'b DolfinxPyFunction_f64<'py, N, O, D>
642 {
643 type Output = DolfinxPyFunction_f64<'py, N, O, D>;
644
645 fn $fn(self, other: Self) -> Self::Output {
646 self.$fn(&other)
647 }
648 }
649
650 impl<'py, 'a, 'b, const N: u32, const O: u32, const D: u32> $trait<&'a Self>
651 for &'b DolfinxPyFunction_f64<'py, N, O, D>
652 {
653 type Output = DolfinxPyFunction_f64<'py, N, O, D>;
654
655 fn $fn(self, other: &'a Self) -> Self::Output {
656 self.with_new_data_and_orig_binop(other, |data_dest, data, data_other| {
657 izip!(data_dest.iter_mut(), data.iter(), data_other.iter()).for_each(
658 |(x, y, z)| {
659 *x = y.$fn(z);
660 },
661 )
662 })
663 }
664 }
665 };
666 }
667
668 macro_rules! impl_binop {
669 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
670 impl_binop_consume!($trait, $fn, $trait_assign, $fn_assign);
671 impl_binop_ref!($trait, $fn);
672 };
673 }
674
675 use std::ops::{Add, AddAssign, Sub, SubAssign};
676
677 impl_binop!(Add, add, AddAssign, add_assign);
678 impl_binop!(Sub, sub, SubAssign, sub_assign);
679
680 //use alg_tools::euclidean::Euclidean;
681 use alg_tools::linops::{VectorSpace, AXPY};
682 //use alg_tools::norms::{Dist, HasDual, Norm, Normed, L2};
683 use alg_tools::self_ownable;
684
685 self_ownable!(DolfinxPyFunction_f64<'py, N, O, D> where 'py, const N: u32, const O: u32, const D: u32);
686
687 impl<'py, const N: u32, const O: u32, const D: u32> Dist<L2, f64>
688 for DolfinxPyFunction_f64<'py, N, O, D>
689 {
690 fn dist<I: Instance<Self>>(&self, other: I, p: L2) -> f64 {
691 other.eval_ref(|v| (self - v).norm(p))
692 }
693 }
694
695 impl<'py, const N: u32, const O: u32, const D: u32> Norm<L2, f64>
696 for DolfinxPyFunction_f64<'py, N, O, D>
697 {
698 fn norm(&self, _p: L2) -> f64 {
699 self.norm2_squared().sqrt()
700 }
701 }
702
703 impl<'py, const N: u32, const O: u32, const D: u32> VectorSpace
704 for DolfinxPyFunction_f64<'py, N, O, D>
705 {
706 type Field = f64;
707 type PrincipalV = Self;
708
709 /// Return a similar zero as `self`.
710 fn similar_origin(&self) -> Self {
711 let mut new = self.similar();
712 new.set_zero();
713 new
714 }
715 }
716
717 impl<'py, const N: u32, const O: u32, const D: u32> AXPY for DolfinxPyFunction_f64<'py, N, O, D> {
718 /// Computes `y = βy + αx`, where `y` is `Self`.
719 fn axpy<I: Instance<Self>>(&mut self, α: Self::Field, x: I, β: Self::Field) {
720 x.eval(|r| {
721 self.with_mut_data_and_orig_binop(
722 r,
723 |data, data_other| {
724 izip!(data.iter_mut(), data_other.iter()).for_each(|(y, x)| {
725 *y = β * (*y) + α * (*x);
726 })
727 },
728 |data, data_orig, data_other| {
729 izip!(data.iter_mut(), data_orig.iter(), data_other.iter()).for_each(
730 |(y, y_orig, x)| {
731 *y = β * (*y_orig) + α * (*x);
732 },
733 )
734 },
735 );
736 })
737 }
738
739 /// Set self to zero.
740 fn set_zero(&mut self) {
741 self.with_mut_data(|data| {
742 data.iter_mut().for_each(|y| *y = 0.0);
743 })
744 }
745 }
746
747 impl<'py, const N: u32, const O: u32, const D: u32> HasDual
748 for DolfinxPyFunction_f64<'py, N, O, D>
749 {
750 type DualSpace = Self;
751
752 fn dual_origin(&self) -> <Self::DualSpace as VectorSpace>::PrincipalV {
753 self.similar_origin()
754 }
755 }
756
757 impl<'py, const N: u32, const O: u32, const D: u32> Normed for DolfinxPyFunction_f64<'py, N, O, D> {
758 type NormExp = L2;
759
760 fn norm_exponent(&self) -> L2 {
761 // TODO: maybe H1 more approriate for O=2
762 L2
763 }
764 }
765
766 impl<'py, const N: u32, const O: u32, const D: u32> Euclidean
767 for DolfinxPyFunction_f64<'py, N, O, D>
768 {
769 type PrincipalE = Self;
770
771 fn dot<I: Instance<Self>>(&self, other: I) -> f64 {
772 let py = self.function_space.py();
773 get_helpers(py)
774 .unwrap()
775 .dot
776 .call1(py, (self.own(), other.own()))
777 .unwrap()
778 .extract(py)
779 .unwrap()
780 /*
781 let data = unsafe { ffi::data_mut_Function_f64(self.cxx) };
782 other.eval_decompose(|x| {
783 let data_other = unsafe { ffi::data_mut_Function_f64(x.cxx) };
784 izip!(data.iter(), data_other.iter())
785 .map(|(y, x)| y * x)
786 .sum()
787 })*/
788 }
789
790 fn norm2_squared(&self) -> f64 {
791 let py = self.function_space.py();
792 get_helpers(py)
793 .unwrap()
794 .norm2_squared
795 .call1(py, (self.own(),))
796 .unwrap()
797 .extract(py)
798 .unwrap()
799 }
800
801 fn dist2_squared<I: Instance<Self>>(&self, other: I) -> f64 {
802 //other.eval_ref(|v| self.norm2_squared() + 2.0 * self.dot(v) + v.norm2_squared_div2())
803 let py = self.function_space.py();
804 get_helpers(py)
805 .unwrap()
806 .dist2_squared
807 .call1(py, (self.own(), other.own()))
808 .unwrap()
809 .extract(py)
810 .unwrap()
811 }
812 }

Mercurial Repository: pointsource_pde

mercurial