1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 use libceed::{prelude::*, Ceed}; 18 19 // ---------------------------------------------------------------------------- 20 // Determine problem size in each dimension from size and dimenison 21 // ---------------------------------------------------------------------------- 22 pub fn cartesian_mesh_size(dim: usize, solution_degree: usize, problem_size: i64) -> [usize; 3] { 23 // Use the approximate formula: 24 // prob_size ~ num_elem * degree^dim 25 let mut num_elem = problem_size / solution_degree.pow(dim as u32) as i64; 26 let mut s = 0; // find s: num_elem / 2 < 2^s <= num_elem 27 while num_elem > 1 { 28 num_elem /= 2; 29 s += 1; 30 } 31 32 // Size per dimension 33 let mut r = s % dim; 34 let mut num_xyz = [0; 3]; 35 for d in 0..dim { 36 let mut sd = s / dim; 37 if r > 0 { 38 sd += 1; 39 r -= 1; 40 } 41 num_xyz[d] = 1 << sd; 42 } 43 num_xyz 44 } 45 46 // ---------------------------------------------------------------------------- 47 // Build element restriction objects for the mesh 48 // ---------------------------------------------------------------------------- 49 pub fn build_cartesian_restriction( 50 ceed: &Ceed, 51 dim: usize, 52 num_xyz: [usize; 3], 53 degree: usize, 54 num_comp: usize, 55 num_qpts: usize, 56 ) -> libceed::Result<(ElemRestriction, ElemRestriction)> { 57 let p = degree + 1; 58 let num_nodes = p.pow(dim as u32); // number of nodes per element 59 let elem_qpts = num_qpts.pow(dim as u32); // number of quadrature pts per element 60 61 // Problem dimensions 62 let mut num_d = [0; 3]; 63 let mut num_elem = 1; 64 let mut scalar_size = 1; 65 for d in 0..dim { 66 num_elem *= num_xyz[d]; 67 num_d[d] = num_xyz[d] * (p - 1) + 1; 68 scalar_size *= num_d[d]; 69 } 70 71 // elem: 0 1 n-1 72 // |---*-...-*---|---*-...-*---|- ... -|--...--| 73 // nodes: 0 1 p-1 p p+1 2*p n*p 74 let mut elem_nodes = vec![0; num_elem * num_nodes]; 75 for e in 0..num_elem { 76 let mut e_xyz = [1; 3]; 77 let mut re = e; 78 for d in 0..dim { 79 e_xyz[d] = re % num_xyz[d]; 80 re /= num_xyz[d]; 81 } 82 let loc_offset = e * num_nodes; 83 for loc_nodes in 0..num_nodes { 84 let mut global_nodes = 0; 85 let mut global_nodes_stride = 1; 86 let mut r_nodes = loc_nodes; 87 for d in 0..dim { 88 global_nodes += (e_xyz[d] * (p - 1) + r_nodes % p) * global_nodes_stride; 89 global_nodes_stride *= num_d[d]; 90 r_nodes /= p; 91 } 92 elem_nodes[loc_offset + loc_nodes] = global_nodes as i32; 93 } 94 } 95 96 // Mesh/solution data restriction 97 let restr = ceed.elem_restriction( 98 num_elem, 99 num_nodes, 100 num_comp, 101 scalar_size, 102 num_comp * scalar_size, 103 MemType::Host, 104 &elem_nodes, 105 )?; 106 107 // Quadratue data restriction 108 let restr_qdata = ceed.strided_elem_restriction( 109 num_elem, 110 elem_qpts, 111 num_comp, 112 num_comp * elem_qpts * num_elem, 113 CEED_STRIDES_BACKEND, 114 )?; 115 Ok((restr, restr_qdata)) 116 } 117 118 // ---------------------------------------------------------------------------- 119 // Set mesh coordinates 120 // ---------------------------------------------------------------------------- 121 pub fn cartesian_mesh_coords( 122 ceed: &Ceed, 123 dim: usize, 124 num_xyz: [usize; 3], 125 mesh_degree: usize, 126 mesh_size: usize, 127 ) -> libceed::Result<Vector> { 128 let p = mesh_degree + 1; 129 let mut num_d = [0; 3]; 130 let mut scalar_size = 1; 131 for d in 0..dim { 132 num_d[d] = num_xyz[d] * (p - 1) + 1; 133 scalar_size *= num_d[d]; 134 } 135 136 // Lobatto points 137 let lobatto_basis = ceed.basis_tensor_H1_Lagrange(1, 1, 2, p, QuadMode::GaussLobatto)?; 138 let nodes_corners = ceed.vector_from_slice(&[0.0, 1.0])?; 139 let mut nodes_full = ceed.vector(p)?; 140 lobatto_basis.apply( 141 1, 142 TransposeMode::NoTranspose, 143 EvalMode::Interp, 144 &nodes_corners, 145 &mut nodes_full, 146 )?; 147 148 // Coordinates for mesh 149 let mut mesh_coords = ceed.vector(mesh_size)?; 150 { 151 let mut coords = mesh_coords.view_mut()?; 152 let nodes = nodes_full.view()?; 153 for gs_nodes in 0..scalar_size { 154 let mut r_nodes = gs_nodes; 155 for d in 0..dim { 156 let d_1d = r_nodes % num_d[d]; 157 coords[gs_nodes + scalar_size * d] = 158 ((d_1d / (p - 1)) as Scalar + nodes[d_1d % (p - 1)]) / num_xyz[d] as Scalar; 159 r_nodes /= num_d[d]; 160 } 161 } 162 } 163 Ok(mesh_coords) 164 } 165 166 // ---------------------------------------------------------------------------- 167