Mercurial > repos > pointsource_pde / file comparison

comparison: src/dolfinx_access.rs

src/dolfinx_access.rs

changeset 12
a0d61ac3719e
parent 3
c3a4f4bb87f7
equal deleted inserted replaced
11:c4f923017349 12:a0d61ac3719e
8 }; 8 };
9 use alg_tools::loc::Loc; 9 use alg_tools::loc::Loc;
10 use anyhow::bail; 10 use anyhow::bail;
11 use pyo3::ffi::PyImport_AppendInittab; 11 use pyo3::ffi::PyImport_AppendInittab;
12 use std::ffi::CString; 12 use std::ffi::CString;
13 use std::mem::MaybeUninit;
14 13
15 /// Import helper required for the linking to the sys crates to actually happen. 14 /// Import helper required for the linking to the sys crates to actually happen.
16 #[allow(unused_imports)] 15 #[allow(unused_imports)]
17 mod dummy_import { 16 mod dummy_import {
18 use dolfinx_sys; 17 use dolfinx_sys;
96 pub unsafe fn minmax_Function_f64_p2( 95 pub unsafe fn minmax_Function_f64_p2(
97 f: *const Function_f64, 96 f: *const Function_f64,
98 max: bool, 97 max: bool,
99 ) -> Result<CoordValuePair>; 98 ) -> Result<CoordValuePair>;
100 99
100 pub unsafe fn cell_dofs_Function_f64_p2(f: *const Function_f64, cell:i32, coords: &mut [f64; 18], vals: &mut [f64; 6]);
101
101 /// Checks if a Dolfinx `Function` is supported (P2 elements on triangular mesh) 102 /// Checks if a Dolfinx `Function` is supported (P2 elements on triangular mesh)
102 pub unsafe fn check_Function_f64(o: *const PyObject) -> bool; 103 pub unsafe fn check_Function_f64(o: *const PyObject) -> bool;
103 104
104 /// Casts a Python object into a Dolfinx `Function` wrapper. 105 /// Casts a Python object into a Dolfinx `Function` wrapper.
105 pub unsafe fn cast_Function_f64(o: *const PyObject) -> Result<*const Function_f64>; 106 pub unsafe fn cast_Function_f64(o: *const PyObject) -> Result<*const Function_f64>;
145 /// Rust-side helper for minimising P2 elements. This bypasses poorly documented basix, 146 /// Rust-side helper for minimising P2 elements. This bypasses poorly documented basix,
146 /// and does things at a low level directly. 147 /// and does things at a low level directly.
147 /// `f` is the function, `simplex_coords` the coorindates of the triangle 148 /// `f` is the function, `simplex_coords` the coorindates of the triangle
148 /// see [`cell_coords_Function_f64_triangle`], `cell` is the cell number, 149 /// see [`cell_coords_Function_f64_triangle`], `cell` is the cell number,
149 /// and `max` indicates whether to maximise or minimise. 150 /// and `max` indicates whether to maximise or minimise.
150 unsafe fn minmax_dolfinx_p2_cell( 151 unsafe fn minmax_dolfinx_p2(
151 f: *const Function_f64, 152 simplex_coords: [f64; 18],
152 simplex_coords: &[f64; 9], 153 vals: [f64; 6],
153 cell: i32,
154 max: bool, 154 max: bool,
155 ) -> CoordValuePair; 155 ) -> CoordValuePair;
156 } 156 }
157 } 157 }
158 158
159 use ffi::CoordValuePair; 159 use ffi::CoordValuePair;
160
161 /// Convert an array of 2D oordinates obtained from Fenics, which hard-codes weird 3D coordinates.
162 #[inline]
163 fn coords_from_fenics(&[x1, x2, _, y1, y2, _, z1, z2, _]: &[f64; 9]) -> [Loc<2, f64>; 3] {
164 [Loc([x1, x2]), Loc([y1, y2]), Loc([z1, z2])]
165 }
166 // fn coords_from_fenics<const N: usize>(coords: &[[f64; 3]; N]) -> [Loc<F, 2>; N] {
167 // coords.map(|[x1, x2, _]| Loc([x1, x2]))
168 // }
169
170 /// Convert an array of 2D oordinates to Fenics, which hard-codes weird 3D coordinates.
171 #[inline]
172 fn coords_to_fenics([Loc([x1, x2]), Loc([y1, y2]), Loc([z1, z2])]: [Loc<2, f64>; 3]) -> [f64; 9] {
173 [x1, x2, 0.0, y1, y2, 0.0, z1, z2, 0.0]
174 }
175 // fn coords_to_fenics<const N: usize>(coords: [Loc<F, 2>; N]) -> [[f64; 3]; N] {
176 // coords.map(|Loc([x1, x2])| [x1, x2, 0.0])
177 // }
178 160
179 /// Given a dolphinx `Function` `f`, as well as 6 coordinates on a simplex in its mesh 161 /// Given a dolphinx `Function` `f`, as well as 6 coordinates on a simplex in its mesh
180 /// (the corners and midpoints), construct a new internal P2 model on `cell`. 162 /// (the corners and midpoints), construct a new internal P2 model on `cell`.
181 /// If `max` is true, the model is negated compared to `f`. 163 /// If `max` is true, the model is negated compared to `f`.
164 #[inline]
182 unsafe fn model_dolfinx_p2_cell( 165 unsafe fn model_dolfinx_p2_cell(
183 f: *const ffi::Function_f64, 166 [p0x, p0y, _, p1x, p1y, _, p2x, p2y, _, p12x, p12y, _, p20x, p20y, _, p01x, p01y, _]: [f64; 18],
184 simplex_coords: &[f64; 9], 167 [v0, v1, v2, v12, v20, v01]: [f64; 6],
185 cell: i32,
186 max: bool, 168 max: bool,
187 ) -> (P2LocalModel<f64, 2, 3>, Simplex<f64, 2, 3>) { 169 ) -> (P2LocalModel<f64, 2, 3>, Simplex<f64, 2, 3>) {
188 // Form simplex 170 // Form simplex
189 let simplex = Simplex(coords_from_fenics(simplex_coords)); 171 let coords = [Loc([p0x, p0y]), Loc([p1x, p1y]), Loc([p2x, p2y])];
190 // Form an array of all coordinates 172 let simplex = Simplex(coords);
191 let midpoints = simplex.midpoints(); 173 let values = if max {
192 let mut all_coords_uninit: [MaybeUninit<f64>; 6 * 3] = [const { MaybeUninit::uninit() }; 6 * 3]; 174 [-v0, -v1, -v2, -v01, -v12, -v20]
193 all_coords_uninit 175 } else {
194 .iter_mut() 176 [v0, v1, v2, v01, v12, v20]
195 .zip( 177 };
196 simplex_coords
197 .iter()
198 .copied()
199 .chain(coords_to_fenics(midpoints).into_iter()),
200 )
201 .for_each(|(a, b)| {
202 a.write(b);
203 });
204 let all_coords = MaybeUninit::array_assume_init(all_coords_uninit);
205 // Evaluate the function, and invert the signs of maximising
206 let mut values = ffi::eval_Function_f64_cell_6(f, &all_coords, cell);
207 if max {
208 values.iter_mut().for_each(|v| *v = -*v);
209 }
210 // Form local model 178 // Form local model
211 let [ref n0, ref n1, ref n2] = simplex.0; 179 let nodes = [Loc([p0x, p0y]), Loc([p1x, p1y]), Loc([p2x, p2y]), Loc([p01x, p01y]), Loc([p12x, p12y]), Loc([p20x, p20y])];
212 let [ref n01, ref n12, ref n20] = midpoints;
213 let nodes = [*n0, *n1, *n2, *n01, *n12, *n20];
214 let model = P2LocalModel::<f64, 2, 3>::new(&nodes, &values); 180 let model = P2LocalModel::<f64, 2, 3>::new(&nodes, &values);
215 181
216 return (model, simplex); 182 return (model, simplex);
217 } 183 }
218 184
219 /// This is a helper routine for `minmax_dolfinx_p2`. It reconstructs a P2 185 /// This is a helper routine for `minmax_dolfinx_p2`. It reconstructs a P2
220 /// model on the cell by evaluating f at 6 points (the simplex corners given 186 /// model on the cell by evaluating f at 6 points (the simplex corners given
221 /// in `simplex_doords`, as well as edge midpoints). Then it minimises this 187 /// in `simplex_doords`, as well as edge midpoints). Then it minimises this
222 /// model. The result should be accurate, given that `f` is assumed to have a 188 /// model. The result should be accurate, given that `f` is assumed to have a
223 /// P2 model on the simplex. 189 /// P2 model on the simplex.
224 pub unsafe fn minmax_dolfinx_p2_cell( 190 pub unsafe fn minmax_dolfinx_p2(
225 f: *const ffi::Function_f64, 191 simplex_coords: [f64; 18],
226 simplex_coords: &[f64; 9], 192 vals: [f64; 6],
227 cell: i32,
228 max: bool, 193 max: bool,
229 ) -> ffi::CoordValuePair { 194 ) -> ffi::CoordValuePair {
230 let (model, simplex) = model_dolfinx_p2_cell(f, simplex_coords, cell, max); 195 let (model, simplex) = model_dolfinx_p2_cell(simplex_coords, vals, max);
231 let (x, v) = model.minimise(&simplex); 196 let (x, v) = model.minimise(&simplex);
232 CoordValuePair { x: x.into(), v } 197 CoordValuePair { x: x.into(), v }
233 } 198 }
234 199
235 pub const NAME: &str = "dolfinx_access"; 200 pub const NAME: &str = "dolfinx_access";

Mercurial Repository: pointsource_pde

mercurial