-
Notifications
You must be signed in to change notification settings - Fork 2
/
element_gradients.c
89 lines (82 loc) · 2.45 KB
/
element_gradients.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#include "element_gradients.h"
#include <assert.h>
#include <string.h>
#include "algebra.h"
#include "jacobian.h"
#include "loop.h"
#include "mesh.h"
#include "tables.h"
#include "tag.h"
LOOP_KERNEL(execute,
unsigned elem_dim,
unsigned const* verts_of_elems,
double const* coords_of_verts,
unsigned ncomps,
double const* comps_of_verts,
unsigned ncomps_out,
double* out,
unsigned verts_per_elem)
unsigned const* verts_of_elem = verts_of_elems + i * verts_per_elem;
double elem_coords[4][3];
double const* elem_comps[4];
for (unsigned j = 0; j < verts_per_elem; ++j) {
unsigned vert = verts_of_elem[j];
elem_comps[j] = comps_of_verts + vert * ncomps;
copy_vector(coords_of_verts + vert * 3, elem_coords[j], 3);
}
double jac[3][3];
element_jacobian(elem_dim, elem_coords, jac);
double jaci[3][3];
invert_jacobian(elem_dim, jac, jaci);
double* grad = out + i * ncomps_out;
for (unsigned j = 0; j < ncomps_out; ++j)
grad[j] = 0;
for (unsigned j = 0; j < elem_dim; ++j)
for (unsigned k = 0; k < ncomps; ++k)
for (unsigned l = 0; l < 3; ++l)
grad[k * 3 + l] += jaci[j][l] *
(elem_comps[j + 1][k] - elem_comps[0][k]);
}
double* element_gradients(
unsigned elem_dim,
unsigned nelems,
unsigned const* verts_of_elems,
double const* coords_of_verts,
unsigned ncomps,
double const* comps_of_verts)
{
assert(elem_dim > 0);
assert(ncomps > 0);
unsigned ncomps_out = ncomps * 3;
double* out = LOOP_MALLOC(double, ncomps_out * nelems);
unsigned verts_per_elem = the_down_degrees[elem_dim][0];
assert(verts_per_elem == elem_dim + 1);
assert(verts_per_elem > 1);
LOOP_EXEC(execute,nelems,
elem_dim,
verts_of_elems,
coords_of_verts,
ncomps,
comps_of_verts,
ncomps_out,
out,
verts_per_elem);
return out;
}
struct const_tag* mesh_element_gradients(
struct mesh* m, char const* name)
{
struct const_tag* t = mesh_find_tag(m, 0, name);
double* data = element_gradients(mesh_dim(m), mesh_count(m, mesh_dim(m)),
mesh_ask_down(m, mesh_dim(m), 0),
mesh_find_tag(m, 0, "coordinates")->d.f64,
t->ncomps, t->d.f64);
static char const* prefix = "grad_";
char* grad_name = LOOP_HOST_MALLOC(char, strlen(t->name) + strlen(prefix) + 1);
strcpy(grad_name, prefix);
strcat(grad_name, t->name);
struct const_tag* out = mesh_add_tag(
m, mesh_dim(m), TAG_F64, grad_name, t->ncomps * 3, data);
loop_host_free(grad_name);
return out;
}