1using LibCEED, Printf 2 3include("common.jl") 4 5function transform_mesh_coords!(dim, mesh_size, mesh_coords) 6 @witharray coords = mesh_coords begin 7 if dim == 1 8 for i = 1:mesh_size 9 # map [0,1] to [0,1] varying the mesh density 10 coords[i] = 0.5 + 1.0/sqrt(3.0)*sin((2.0/3.0)*pi*(coords[i] - 0.5)) 11 end 12 exact_volume = 1.0 13 else 14 num_nodes = mesh_size÷dim 15 @inbounds @simd for i = 1:num_nodes 16 # map (x,y) from [0,1]x[0,1] to the quarter annulus with polar 17 # coordinates, (r,phi) in [1,2]x[0,pi/2] with area = 3/4*pi 18 u = coords[i] 19 v = coords[i+num_nodes] 20 u = 1.0 + u 21 v = pi/2*v 22 coords[i] = u*cos(v) 23 coords[i+num_nodes] = u*sin(v) 24 end 25 exact_volume = 3.0/4.0*pi 26 end 27 return exact_volume 28 end 29end 30 31function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size) 32 ncompx = dim 33 prob_size < 0 && (prob_size = 256*1024) 34 35 ceed = Ceed(ceed_spec) 36 mesh_basis = 37 create_tensor_h1_lagrange_basis(ceed, dim, ncompx, mesh_order + 1, num_qpts, GAUSS) 38 sol_basis = 39 create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS) 40 41 # Determine the mesh size based on the given approximate problem size. 42 nxyz = get_cartesian_mesh_size(dim, sol_order, prob_size) 43 println("Mesh size: ", nxyz) 44 45 # Build CeedElemRestriction objects describing the mesh and solution discrete 46 # representations. 47 mesh_size, mesh_rstr, _ = 48 build_cartesian_restriction(ceed, dim, nxyz, mesh_order, ncompx, num_qpts) 49 num_q_comp = 1 + div(dim*(dim + 1), 2) 50 sol_size, _, qdata_rstr_i = build_cartesian_restriction( 51 ceed, 52 dim, 53 nxyz, 54 sol_order, 55 num_q_comp, 56 num_qpts, 57 mode=StridedOnly, 58 ) 59 sol_size, sol_rstr, sol_rstr_i = build_cartesian_restriction( 60 ceed, 61 dim, 62 nxyz, 63 sol_order, 64 1, 65 num_qpts, 66 mode=RestrictionAndStrided, 67 ) 68 println("Number of mesh nodes : ", div(mesh_size, dim)) 69 println("Number of solution nodes : ", sol_size) 70 71 # Create a CeedVector with the mesh coordinates. 72 mesh_coords = CeedVector(ceed, mesh_size) 73 set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords) 74 # Apply a transformation to the mesh. 75 exact_vol = transform_mesh_coords!(dim, mesh_size, mesh_coords) 76 77 #Create the Q-function that builds the mass+diffusion operator ( i.e it computes the quadrature data) and set its context data. 78 @interior_qf build_qfunc = ( 79 ceed, 80 dim=dim, 81 (dx, :in, EVAL_GRAD, dim, dim), # ← THIS LINE: dx input 82 (weights, :in, EVAL_WEIGHT), # ← weights input 83 (qdata, :out, EVAL_NONE, num_q_comp), # ← qdata output 84 begin 85 # Compute determinant 86 det_J = det(dx) 87 88 # Store mass component 89 qdata[1] = weights*det_J 90 91 # Store diffusion components (J^T * J) 92 idx = 2 93 for i = 1:dim 94 for j = i:dim 95 qdata[idx] = dx[:, i]'*dx[:, j] 96 idx += 1 97 end 98 end 99 end, 100 ) 101 102 # Create the operator that builds the quadrature data for the mass+diffusion operator. 103 build_oper = Operator( 104 ceed, 105 qf=build_qfunc, 106 fields=[ 107 (:dx, mesh_rstr, mesh_basis, CeedVectorActive()), 108 (:weights, ElemRestrictionNone(), mesh_basis, CeedVectorNone()), 109 (:qdata, qdata_rstr_i, BasisNone(), CeedVectorActive()), 110 ], 111 ) 112 113 # Compute the quadrature data for the mass+diff operator. 114 elem_qpts = num_qpts^dim 115 num_elem = prod(nxyz) 116 qdata = CeedVector(ceed, num_elem*elem_qpts*num_q_comp) 117 print("Computing the quadrature data for the mass+diffusion operator ...") 118 flush(stdout) 119 apply!(build_oper, mesh_coords, qdata) 120 println(" done.") 121 122 # Create the Q-function that defines the action of the mass+diffusion operator. 123 @interior_qf apply_qfunc = ( 124 ceed, 125 dim=dim, 126 (u, :in, EVAL_INTERP), 127 (du, :in, EVAL_GRAD, dim), 128 (qdata, :in, EVAL_NONE, num_q_comp), 129 (v, :out, EVAL_INTERP), 130 (dv, :out, EVAL_GRAD, dim), 131 begin 132 # Apply mass: v = qdata[1] * u 133 v .= qdata[1].*u 134 135 # Apply diffusion: dv = (qdata[2:end]) * du 136 # The qdata contains the symmetric diffusion tensor (J^T*J) 137 # dv_i = sum_j (J^T*J)_{i,j} * du_j 138 139 # For efficiency, rebuild the matrix from stored components 140 idx = 2 141 for i = 1:dim 142 dv_i = 0.0 143 for j = 1:dim 144 # Reconstruct symmetric matrix element 145 if j >= i 146 mat_idx = idx + div((j - 1)*j, 2) + (i - 1) 147 else 148 mat_idx = idx + div((i - 1)*i, 2) + (j - 1) 149 end 150 dv_i += qdata[mat_idx]*du[j] 151 end 152 dv[i] = dv_i 153 end 154 end, 155 ) 156 157 # Create the mass+diffusion operator. 158 oper = Operator( 159 ceed, 160 qf=apply_qfunc, 161 fields=[ 162 (:u, sol_rstr, sol_basis, CeedVectorActive()), 163 (:du, sol_rstr, sol_basis, CeedVectorActive()), 164 (:qdata, qdata_rstr_i, BasisNone(), qdata), 165 (:v, sol_rstr, sol_basis, CeedVectorActive()), 166 (:dv, sol_rstr, sol_basis, CeedVectorActive()), 167 ], 168 ) 169 170 # Compute the mesh volume using the mass+diffusion operator: vol = 1^T \cdot (M + K) \cdot 1 171 print("Computing the mesh volume using the formula: vol = 1^T * (M + K) * 1...") 172 flush(stdout) 173 # Create auxiliary solution-size vectors. 174 u = CeedVector(ceed, sol_size) 175 v = CeedVector(ceed, sol_size) 176 # Initialize 'u' with ones. 177 u[] = 1.0 178 # Apply the mass+diffusion operator: 'u' -> 'v'. 179 apply!(oper, u, v) 180 # Compute and print the sum of the entries of 'v' giving the mesh volume. 181 vol = witharray_read(sum, v, MEM_HOST) 182 183 println(" done.") 184 @printf("Exact mesh volume : % .14g\n", exact_vol) 185 @printf("Computed mesh volume : % .14g\n", vol) 186 @printf("Volume error : % .14g\n", vol - exact_vol) 187end 188 189# Entry point 190run_ex3( 191 ceed_spec="/cpu/self", 192 dim=3, 193 mesh_order=4, 194 sol_order=4, 195 num_qpts=4 + 2, 196 prob_size=-1, 197) 198