Newer
Older
#include "types.h"
#include "entity_data.h"
/* Field derivatives kernel, case for elements with planar faces */
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
void compute_field_derivatives(const entity_data_cpu<1>& ed,
const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz)
{
/* The elements must be ordered by orientation -> entity::sort_by_orientation() */
size_t orient_elem_base = 0;
for (size_t iO = 0; iO < ed.num_orientations; iO++)
{
size_t num_elems = ed.num_elems[iO];
size_t orient_base = orient_elem_base * ed.num_bf;
for (size_t iT = 0; iT < num_elems; iT++)
{
size_t ent_iT = orient_elem_base + iT;
size_t dofofs = ed.dofs_base + orient_base + iT*ed.num_bf;
for (size_t odof = 0; odof < ed.num_bf; odof++)
{
double acc_df_dx = 0.0;
double acc_df_dy = 0.0;
double acc_df_dz = 0.0;
for (size_t refsys_dir = 0; refsys_dir < 3; refsys_dir++)
{
for (size_t idof = 0; idof < ed.num_bf; idof++)
{
auto dm_row = iO*ed.num_bf + odof;
auto dm_col = refsys_dir*ed.num_bf + idof;
double D = ed.differentiation_matrix(dm_row, dm_col);
acc_df_dx += ed.jacobians(3*ent_iT + 0, refsys_dir) * D * f(dofofs+idof);
acc_df_dy += ed.jacobians(3*ent_iT + 1, refsys_dir) * D * f(dofofs+idof);
acc_df_dz += ed.jacobians(3*ent_iT + 2, refsys_dir) * D * f(dofofs+idof);
}
}
df_dx(dofofs+odof) = acc_df_dx;
df_dy(dofofs+odof) = acc_df_dy;
df_dz(dofofs+odof) = acc_df_dz;
}
}
orient_elem_base += num_elems;
}
}