ex2-surface.c (ce18bed930e8f3bfebcf709a18844aba97fe4630) ex2-surface.c (990fdeb6bb8fc9af2da4472bdc0d1f57da5da0e5)
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// libCEED Example 2

--- 34 unchanged lines hidden (view full) ---

43
44#include <ceed.h>
45#include <math.h>
46#include <stdlib.h>
47#include <string.h>
48#include "ex2-surface.h"
49
50// Auxiliary functions.
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// libCEED Example 2

--- 34 unchanged lines hidden (view full) ---

43
44#include <ceed.h>
45#include <math.h>
46#include <stdlib.h>
47#include <string.h>
48#include "ex2-surface.h"
49
50// Auxiliary functions.
51int GetCartesianMeshSize(int dim, int degree, int prob_size, int num_xyz[3]);
52int BuildCartesianRestriction(Ceed ceed, int dim, int num_xyz[3], int degree,
53 int num_comp, CeedInt *size, CeedInt num_qpts,
54 CeedElemRestriction *restr,
51int GetCartesianMeshSize(CeedInt dim, CeedInt degree, CeedInt prob_size,
52 CeedInt num_xyz[3]);
53int BuildCartesianRestriction(Ceed ceed, CeedInt dim, CeedInt num_xyz[3],
54 CeedInt degree, CeedInt num_comp, CeedInt *size,
55 CeedInt num_qpts, CeedElemRestriction *restr,
55 CeedElemRestriction *restr_i);
56 CeedElemRestriction *restr_i);
56int SetCartesianMeshCoords(int dim, int num_xyz[3], int mesh_degree,
57int SetCartesianMeshCoords(CeedInt dim, CeedInt num_xyz[3], CeedInt mesh_degree,
57 CeedVector mesh_coords);
58 CeedVector mesh_coords);
58CeedScalar TransformMeshCoords(int dim, int mesh_size, CeedVector mesh_coords);
59CeedScalar TransformMeshCoords(CeedInt dim, CeedInt mesh_size,
60 CeedVector mesh_coords);
59
60int main(int argc, const char *argv[]) {
61 const char *ceed_spec = "/cpu/self";
61
62int main(int argc, const char *argv[]) {
63 const char *ceed_spec = "/cpu/self";
62 int dim = 3; // dimension of the mesh
63 int num_comp_x = 3; // number of x components
64 int mesh_degree = 4; // polynomial degree for the mesh
65 int sol_degree = 4; // polynomial degree for the solution
66 int num_qpts = sol_degree + 2; // number of 1D quadrature points
67 int prob_size = -1; // approximate problem size
68 int help = 0, test = 0, gallery = 0;
64 CeedInt dim = 3; // dimension of the mesh
65 CeedInt num_comp_x = 3; // number of x components
66 CeedInt mesh_degree = 4; // polynomial degree for the mesh
67 CeedInt sol_degree = 4; // polynomial degree for the solution
68 CeedInt num_qpts = sol_degree + 2; // number of 1D quadrature points
69 CeedInt prob_size = -1; // approximate problem size
70 CeedInt help = 0, test = 0, gallery = 0;
69
70 // Process command line arguments.
71 for (int ia = 1; ia < argc; ia++) {
72 // LCOV_EXCL_START
73 int next_arg = ((ia+1) < argc), parse_error = 0;
74 if (!strcmp(argv[ia],"-h")) {
75 help = 1;
76 } else if (!strcmp(argv[ia],"-c") || !strcmp(argv[ia],"-ceed")) {

--- 26 unchanged lines hidden (view full) ---

103 mesh_degree = fmax(mesh_degree, sol_degree);
104 sol_degree = mesh_degree;
105
106 // Print the values of all options:
107 if (!test || help) {
108 // LCOV_EXCL_START
109 printf("Selected options: [command line option] : <current value>\n");
110 printf(" Ceed specification [-c] : %s\n", ceed_spec);
71
72 // Process command line arguments.
73 for (int ia = 1; ia < argc; ia++) {
74 // LCOV_EXCL_START
75 int next_arg = ((ia+1) < argc), parse_error = 0;
76 if (!strcmp(argv[ia],"-h")) {
77 help = 1;
78 } else if (!strcmp(argv[ia],"-c") || !strcmp(argv[ia],"-ceed")) {

--- 26 unchanged lines hidden (view full) ---

105 mesh_degree = fmax(mesh_degree, sol_degree);
106 sol_degree = mesh_degree;
107
108 // Print the values of all options:
109 if (!test || help) {
110 // LCOV_EXCL_START
111 printf("Selected options: [command line option] : <current value>\n");
112 printf(" Ceed specification [-c] : %s\n", ceed_spec);
111 printf(" Mesh dimension [-d] : %d\n", dim);
112 printf(" Mesh degree [-m] : %d\n", mesh_degree);
113 printf(" Solution degree [-p] : %d\n", sol_degree);
114 printf(" Num. 1D quadr. pts [-q] : %d\n", num_qpts);
115 printf(" Approx. # unknowns [-s] : %d\n", prob_size);
113 printf(" Mesh dimension [-d] : %" CeedInt_FMT "\n", dim);
114 printf(" Mesh degree [-m] : %" CeedInt_FMT "\n", mesh_degree);
115 printf(" Solution degree [-p] : %" CeedInt_FMT "\n", sol_degree);
116 printf(" Num. 1D quadr. pts [-q] : %" CeedInt_FMT "\n", num_qpts);
117 printf(" Approx. # unknowns [-s] : %" CeedInt_FMT "\n", prob_size);
116 printf(" QFunction source [-g] : %s\n", gallery?"gallery":"header");
117 if (help) {
118 printf("Test/quiet mode is %s\n", (test?"ON":"OFF (use -t to enable)"));
119 return 0;
120 }
121 printf("\n");
122 // LCOV_EXCL_STOP
123 }

--- 6 unchanged lines hidden (view full) ---

130 // Construct the mesh and solution bases.
131 CeedBasis mesh_basis, sol_basis;
132 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, mesh_degree + 1,
133 num_qpts, CEED_GAUSS, &mesh_basis);
134 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, sol_degree + 1, num_qpts,
135 CEED_GAUSS, &sol_basis);
136
137 // Determine the mesh size based on the given approximate problem size.
118 printf(" QFunction source [-g] : %s\n", gallery?"gallery":"header");
119 if (help) {
120 printf("Test/quiet mode is %s\n", (test?"ON":"OFF (use -t to enable)"));
121 return 0;
122 }
123 printf("\n");
124 // LCOV_EXCL_STOP
125 }

--- 6 unchanged lines hidden (view full) ---

132 // Construct the mesh and solution bases.
133 CeedBasis mesh_basis, sol_basis;
134 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, mesh_degree + 1,
135 num_qpts, CEED_GAUSS, &mesh_basis);
136 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, sol_degree + 1, num_qpts,
137 CEED_GAUSS, &sol_basis);
138
139 // Determine the mesh size based on the given approximate problem size.
138 int num_xyz[3];
140 CeedInt num_xyz[3];
139 GetCartesianMeshSize(dim, sol_degree, prob_size, num_xyz);
140
141 if (!test) {
142 // LCOV_EXCL_START
141 GetCartesianMeshSize(dim, sol_degree, prob_size, num_xyz);
142
143 if (!test) {
144 // LCOV_EXCL_START
143 printf("Mesh size: nx = %d", num_xyz[0]);
144 if (dim > 1) { printf(", ny = %d", num_xyz[1]); }
145 if (dim > 2) { printf(", nz = %d", num_xyz[2]); }
145 printf("Mesh size: nx = %" CeedInt_FMT, num_xyz[0]);
146 if (dim > 1) { printf(", ny = %" CeedInt_FMT, num_xyz[1]); }
147 if (dim > 2) { printf(", nz = %" CeedInt_FMT, num_xyz[2]); }
146 printf("\n");
147 // LCOV_EXCL_STOP
148 }
149
150 // Build CeedElemRestriction objects describing the mesh and solution discrete
151 // representations.
152 CeedInt mesh_size, sol_size;
153 CeedElemRestriction mesh_restr, sol_restr, q_data_restr_i;
154 BuildCartesianRestriction(ceed, dim, num_xyz, mesh_degree, num_comp_x,
155 &mesh_size, num_qpts, &mesh_restr, NULL);
156 BuildCartesianRestriction(ceed, dim, num_xyz, sol_degree, dim*(dim+1)/2,
157 &sol_size, num_qpts, NULL, &q_data_restr_i);
158 BuildCartesianRestriction(ceed, dim, num_xyz, sol_degree, 1, &sol_size,
159 num_qpts, &sol_restr, NULL);
160 if (!test) {
161 // LCOV_EXCL_START
148 printf("\n");
149 // LCOV_EXCL_STOP
150 }
151
152 // Build CeedElemRestriction objects describing the mesh and solution discrete
153 // representations.
154 CeedInt mesh_size, sol_size;
155 CeedElemRestriction mesh_restr, sol_restr, q_data_restr_i;
156 BuildCartesianRestriction(ceed, dim, num_xyz, mesh_degree, num_comp_x,
157 &mesh_size, num_qpts, &mesh_restr, NULL);
158 BuildCartesianRestriction(ceed, dim, num_xyz, sol_degree, dim*(dim+1)/2,
159 &sol_size, num_qpts, NULL, &q_data_restr_i);
160 BuildCartesianRestriction(ceed, dim, num_xyz, sol_degree, 1, &sol_size,
161 num_qpts, &sol_restr, NULL);
162 if (!test) {
163 // LCOV_EXCL_START
162 printf("Number of mesh nodes : %d\n", mesh_size/dim);
163 printf("Number of solution nodes : %d\n", sol_size);
164 printf("Number of mesh nodes : %" CeedInt_FMT "\n", mesh_size/dim);
165 printf("Number of solution nodes : %" CeedInt_FMT "\n", sol_size);
164 // LCOV_EXCL_STOP
165 }
166
167 // Create a CeedVector with the mesh coordinates.
168 CeedVector mesh_coords;
169 CeedVectorCreate(ceed, mesh_size, &mesh_coords);
170 SetCartesianMeshCoords(dim, num_xyz, mesh_degree, mesh_coords);
171

--- 19 unchanged lines hidden (view full) ---

191 CeedQFunctionAddInput(qf_build, "dx", num_comp_x*dim, CEED_EVAL_GRAD);
192 CeedQFunctionAddInput(qf_build, "weights", 1, CEED_EVAL_WEIGHT);
193 CeedQFunctionAddOutput(qf_build, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE);
194 CeedQFunctionSetContext(qf_build, build_ctx);
195 break;
196 case 1: {
197 // This creates the QFunction via the gallery.
198 char name[16] = "";
166 // LCOV_EXCL_STOP
167 }
168
169 // Create a CeedVector with the mesh coordinates.
170 CeedVector mesh_coords;
171 CeedVectorCreate(ceed, mesh_size, &mesh_coords);
172 SetCartesianMeshCoords(dim, num_xyz, mesh_degree, mesh_coords);
173

--- 19 unchanged lines hidden (view full) ---

193 CeedQFunctionAddInput(qf_build, "dx", num_comp_x*dim, CEED_EVAL_GRAD);
194 CeedQFunctionAddInput(qf_build, "weights", 1, CEED_EVAL_WEIGHT);
195 CeedQFunctionAddOutput(qf_build, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE);
196 CeedQFunctionSetContext(qf_build, build_ctx);
197 break;
198 case 1: {
199 // This creates the QFunction via the gallery.
200 char name[16] = "";
199 snprintf(name, sizeof name, "Poisson%dDBuild", dim);
201 snprintf(name, sizeof name, "Poisson%" CeedInt_FMT "DBuild", dim);
200 CeedQFunctionCreateInteriorByName(ceed, name, &qf_build);
201 break;
202 }
203 }
204
205 // Create the operator that builds the quadrature data for the diffusion
206 // operator.
207 CeedOperator op_build;

--- 5 unchanged lines hidden (view full) ---

213 mesh_basis, CEED_VECTOR_NONE);
214 CeedOperatorSetField(op_build, "qdata", q_data_restr_i,
215 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
216
217 // Compute the quadrature data for the diffusion operator.
218 CeedVector q_data;
219 CeedInt elem_qpts = CeedIntPow(num_qpts, dim);
220 CeedInt num_elem = 1;
202 CeedQFunctionCreateInteriorByName(ceed, name, &qf_build);
203 break;
204 }
205 }
206
207 // Create the operator that builds the quadrature data for the diffusion
208 // operator.
209 CeedOperator op_build;

--- 5 unchanged lines hidden (view full) ---

215 mesh_basis, CEED_VECTOR_NONE);
216 CeedOperatorSetField(op_build, "qdata", q_data_restr_i,
217 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
218
219 // Compute the quadrature data for the diffusion operator.
220 CeedVector q_data;
221 CeedInt elem_qpts = CeedIntPow(num_qpts, dim);
222 CeedInt num_elem = 1;
221 for (int d = 0; d < dim; d++)
223 for (CeedInt d = 0; d < dim; d++)
222 num_elem *= num_xyz[d];
223 CeedVectorCreate(ceed, num_elem*elem_qpts*dim*(dim+1)/2, &q_data);
224 CeedOperatorApply(op_build, mesh_coords, q_data,
225 CEED_REQUEST_IMMEDIATE);
226
227 // Create the QFunction that defines the action of the diffusion operator.
228 CeedQFunction qf_apply;
229 switch (gallery) {

--- 4 unchanged lines hidden (view full) ---

234 CeedQFunctionAddInput(qf_apply, "du", dim, CEED_EVAL_GRAD);
235 CeedQFunctionAddInput(qf_apply, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE);
236 CeedQFunctionAddOutput(qf_apply, "dv", dim, CEED_EVAL_GRAD);
237 CeedQFunctionSetContext(qf_apply, build_ctx);
238 break;
239 case 1: {
240 // This creates the QFunction via the gallery.
241 char name[16] = "";
224 num_elem *= num_xyz[d];
225 CeedVectorCreate(ceed, num_elem*elem_qpts*dim*(dim+1)/2, &q_data);
226 CeedOperatorApply(op_build, mesh_coords, q_data,
227 CEED_REQUEST_IMMEDIATE);
228
229 // Create the QFunction that defines the action of the diffusion operator.
230 CeedQFunction qf_apply;
231 switch (gallery) {

--- 4 unchanged lines hidden (view full) ---

236 CeedQFunctionAddInput(qf_apply, "du", dim, CEED_EVAL_GRAD);
237 CeedQFunctionAddInput(qf_apply, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE);
238 CeedQFunctionAddOutput(qf_apply, "dv", dim, CEED_EVAL_GRAD);
239 CeedQFunctionSetContext(qf_apply, build_ctx);
240 break;
241 case 1: {
242 // This creates the QFunction via the gallery.
243 char name[16] = "";
242 snprintf(name, sizeof name, "Poisson%dDApply", dim);
244 snprintf(name, sizeof name, "Poisson%" CeedInt_FMT "DApply", dim);
243 CeedQFunctionCreateInteriorByName(ceed, name, &qf_apply);
244 break;
245 }
246 }
247
248 // Create the diffusion operator.
249 CeedOperator op_apply;
250 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE,

--- 62 unchanged lines hidden (view full) ---

313 CeedElemRestrictionDestroy(&mesh_restr);
314 CeedElemRestrictionDestroy(&q_data_restr_i);
315 CeedBasisDestroy(&sol_basis);
316 CeedBasisDestroy(&mesh_basis);
317 CeedDestroy(&ceed);
318 return 0;
319}
320
245 CeedQFunctionCreateInteriorByName(ceed, name, &qf_apply);
246 break;
247 }
248 }
249
250 // Create the diffusion operator.
251 CeedOperator op_apply;
252 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE,

--- 62 unchanged lines hidden (view full) ---

315 CeedElemRestrictionDestroy(&mesh_restr);
316 CeedElemRestrictionDestroy(&q_data_restr_i);
317 CeedBasisDestroy(&sol_basis);
318 CeedBasisDestroy(&mesh_basis);
319 CeedDestroy(&ceed);
320 return 0;
321}
322
321int GetCartesianMeshSize(int dim, int degree, int prob_size, int num_xyz[3]) {
323int GetCartesianMeshSize(CeedInt dim, CeedInt degree, CeedInt prob_size,
324 CeedInt num_xyz[3]) {
322 // Use the approximate formula:
323 // prob_size ~ num_elem * degree^dim
324 CeedInt num_elem = prob_size / CeedIntPow(degree, dim);
325 CeedInt s = 0; // find s: num_elem/2 < 2^s <= num_elem
326 while (num_elem > 1) {
327 num_elem /= 2;
328 s++;
329 }
330 CeedInt r = s%dim;
325 // Use the approximate formula:
326 // prob_size ~ num_elem * degree^dim
327 CeedInt num_elem = prob_size / CeedIntPow(degree, dim);
328 CeedInt s = 0; // find s: num_elem/2 < 2^s <= num_elem
329 while (num_elem > 1) {
330 num_elem /= 2;
331 s++;
332 }
333 CeedInt r = s%dim;
331 for (int d = 0; d < dim; d++) {
332 int sd = s/dim;
334 for (CeedInt d = 0; d < dim; d++) {
335 CeedInt sd = s/dim;
333 if (r > 0) { sd++; r--; }
334 num_xyz[d] = 1 << sd;
335 }
336 return 0;
337}
338
336 if (r > 0) { sd++; r--; }
337 num_xyz[d] = 1 << sd;
338 }
339 return 0;
340}
341
339int BuildCartesianRestriction(Ceed ceed, int dim, int num_xyz[3], int degree,
340 int num_comp, CeedInt *size, CeedInt num_qpts,
341 CeedElemRestriction *restr,
342int BuildCartesianRestriction(Ceed ceed, CeedInt dim, CeedInt num_xyz[3],
343 CeedInt degree, CeedInt num_comp, CeedInt *size,
344 CeedInt num_qpts, CeedElemRestriction *restr,
342 CeedElemRestriction *restr_i) {
343 CeedInt p = degree + 1;
344 CeedInt num_nodes = CeedIntPow(p, dim); // number of scalar nodes per element
345 CeedInt elem_qpts = CeedIntPow(num_qpts, dim); // number of qpts per element
346 CeedInt nd[3], num_elem = 1, scalar_size = 1;
345 CeedElemRestriction *restr_i) {
346 CeedInt p = degree + 1;
347 CeedInt num_nodes = CeedIntPow(p, dim); // number of scalar nodes per element
348 CeedInt elem_qpts = CeedIntPow(num_qpts, dim); // number of qpts per element
349 CeedInt nd[3], num_elem = 1, scalar_size = 1;
347 for (int d = 0; d < dim; d++) {
350 for (CeedInt d = 0; d < dim; d++) {
348 num_elem *= num_xyz[d];
349 nd[d] = num_xyz[d] * (p - 1) + 1;
350 scalar_size *= nd[d];
351 }
352 *size = scalar_size*num_comp;
353 // elem: 0 1 n-1
354 // |---*-...-*---|---*-...-*---|- ... -|--...--|
355 // num_nodes: 0 1 p-1 p p+1 2*p n*p
356 CeedInt *el_nodes = malloc(sizeof(CeedInt)*num_elem*num_nodes);
357 for (CeedInt e = 0; e < num_elem; e++) {
358 CeedInt e_xyz[3] = {1, 1, 1}, re = e;
351 num_elem *= num_xyz[d];
352 nd[d] = num_xyz[d] * (p - 1) + 1;
353 scalar_size *= nd[d];
354 }
355 *size = scalar_size*num_comp;
356 // elem: 0 1 n-1
357 // |---*-...-*---|---*-...-*---|- ... -|--...--|
358 // num_nodes: 0 1 p-1 p p+1 2*p n*p
359 CeedInt *el_nodes = malloc(sizeof(CeedInt)*num_elem*num_nodes);
360 for (CeedInt e = 0; e < num_elem; e++) {
361 CeedInt e_xyz[3] = {1, 1, 1}, re = e;
359 for (int d = 0; d < dim; d++) { e_xyz[d] = re%num_xyz[d]; re /= num_xyz[d]; }
362 for (CeedInt d = 0; d < dim; d++) { e_xyz[d] = re%num_xyz[d]; re /= num_xyz[d]; }
360 CeedInt *loc_el_nodes = el_nodes + e*num_nodes;
363 CeedInt *loc_el_nodes = el_nodes + e*num_nodes;
361 for (int l_nodes = 0; l_nodes < num_nodes; l_nodes++) {
364 for (CeedInt l_nodes = 0; l_nodes < num_nodes; l_nodes++) {
362 CeedInt g_nodes = 0, g_nodes_stride = 1, r_nodes = l_nodes;
365 CeedInt g_nodes = 0, g_nodes_stride = 1, r_nodes = l_nodes;
363 for (int d = 0; d < dim; d++) {
366 for (CeedInt d = 0; d < dim; d++) {
364 g_nodes += (e_xyz[d] * (p - 1) + r_nodes % p) * g_nodes_stride;
365 g_nodes_stride *= nd[d];
366 r_nodes /= p;
367 }
368 loc_el_nodes[l_nodes] = g_nodes;
369 }
370 }
371 if (restr)

--- 6 unchanged lines hidden (view full) ---

378 CeedElemRestrictionCreateStrided(ceed, num_elem, elem_qpts,
379 num_comp, num_comp * elem_qpts * num_elem,
380 CEED_STRIDES_BACKEND, restr_i);
381 }
382
383 return 0;
384}
385
367 g_nodes += (e_xyz[d] * (p - 1) + r_nodes % p) * g_nodes_stride;
368 g_nodes_stride *= nd[d];
369 r_nodes /= p;
370 }
371 loc_el_nodes[l_nodes] = g_nodes;
372 }
373 }
374 if (restr)

--- 6 unchanged lines hidden (view full) ---

381 CeedElemRestrictionCreateStrided(ceed, num_elem, elem_qpts,
382 num_comp, num_comp * elem_qpts * num_elem,
383 CEED_STRIDES_BACKEND, restr_i);
384 }
385
386 return 0;
387}
388
386int SetCartesianMeshCoords(int dim, int num_xyz[3], int mesh_degree,
389int SetCartesianMeshCoords(CeedInt dim, CeedInt num_xyz[3], CeedInt mesh_degree,
387 CeedVector mesh_coords) {
388 CeedInt p = mesh_degree + 1;
389 CeedInt nd[3], num_elem = 1, scalar_size = 1;
390 CeedVector mesh_coords) {
391 CeedInt p = mesh_degree + 1;
392 CeedInt nd[3], num_elem = 1, scalar_size = 1;
390 for (int d = 0; d < dim; d++) {
393 for (CeedInt d = 0; d < dim; d++) {
391 num_elem *= num_xyz[d];
392 nd[d] = num_xyz[d] * (p - 1) + 1;
393 scalar_size *= nd[d];
394 }
395 CeedScalar *coords;
396 CeedVectorGetArrayWrite(mesh_coords, CEED_MEM_HOST, &coords);
397 CeedScalar *nodes = malloc(sizeof(CeedScalar) * p);
398 // The H1 basis uses Lobatto quadrature points as nodes.
399 CeedLobattoQuadrature(p, nodes, NULL); // nodes are in [-1,1]
400 for (CeedInt i = 0; i < p; i++) { nodes[i] = 0.5 + 0.5 * nodes[i]; }
401 for (CeedInt gs_nodes = 0; gs_nodes < scalar_size; gs_nodes++) {
402 CeedInt r_nodes = gs_nodes;
394 num_elem *= num_xyz[d];
395 nd[d] = num_xyz[d] * (p - 1) + 1;
396 scalar_size *= nd[d];
397 }
398 CeedScalar *coords;
399 CeedVectorGetArrayWrite(mesh_coords, CEED_MEM_HOST, &coords);
400 CeedScalar *nodes = malloc(sizeof(CeedScalar) * p);
401 // The H1 basis uses Lobatto quadrature points as nodes.
402 CeedLobattoQuadrature(p, nodes, NULL); // nodes are in [-1,1]
403 for (CeedInt i = 0; i < p; i++) { nodes[i] = 0.5 + 0.5 * nodes[i]; }
404 for (CeedInt gs_nodes = 0; gs_nodes < scalar_size; gs_nodes++) {
405 CeedInt r_nodes = gs_nodes;
403 for (int d = 0; d < dim; d++) {
406 for (CeedInt d = 0; d < dim; d++) {
404 CeedInt d1d = r_nodes % nd[d];
405 coords[gs_nodes + scalar_size * d] = ((d1d / (p - 1)) + nodes[d1d %
406 (p - 1)]) / num_xyz[d];
407 r_nodes /= nd[d];
408 }
409 }
410 free(nodes);
411 CeedVectorRestoreArray(mesh_coords, &coords);
412 return 0;
413}
414
415#ifndef M_PI
416#define M_PI 3.14159265358979323846
417#endif
418
407 CeedInt d1d = r_nodes % nd[d];
408 coords[gs_nodes + scalar_size * d] = ((d1d / (p - 1)) + nodes[d1d %
409 (p - 1)]) / num_xyz[d];
410 r_nodes /= nd[d];
411 }
412 }
413 free(nodes);
414 CeedVectorRestoreArray(mesh_coords, &coords);
415 return 0;
416}
417
418#ifndef M_PI
419#define M_PI 3.14159265358979323846
420#endif
421
419CeedScalar TransformMeshCoords(int dim, int mesh_size, CeedVector mesh_coords) {
422CeedScalar TransformMeshCoords(CeedInt dim, CeedInt mesh_size,
423 CeedVector mesh_coords) {
420 CeedScalar exact_sa = (dim == 1 ? 2 : dim == 2 ? 4 : 6);
421 CeedScalar *coords;
422
423 CeedVectorGetArray(mesh_coords, CEED_MEM_HOST, &coords);
424 for (CeedInt i = 0; i < mesh_size; i++) {
425 // map [0,1] to [0,1] varying the mesh density
426 coords[i] = 0.5 + 1./sqrt(3.) * sin((2./3.) * M_PI * (coords[i] - 0.5));
427 }
428 CeedVectorRestoreArray(mesh_coords, &coords);
429
430 return exact_sa;
431}
424 CeedScalar exact_sa = (dim == 1 ? 2 : dim == 2 ? 4 : 6);
425 CeedScalar *coords;
426
427 CeedVectorGetArray(mesh_coords, CEED_MEM_HOST, &coords);
428 for (CeedInt i = 0; i < mesh_size; i++) {
429 // map [0,1] to [0,1] varying the mesh density
430 coords[i] = 0.5 + 1./sqrt(3.) * sin((2./3.) * M_PI * (coords[i] - 0.5));
431 }
432 CeedVectorRestoreArray(mesh_coords, &coords);
433
434 return exact_sa;
435}