v0.14.0
Loading...
Searching...
No Matches
dd_boundary_Tensor2_symmetric.hpp
Go to the documentation of this file.
1/* Takes a second derivative of a Tensor2_symmetric, yielding a
2 Ddg. This is mostly useful near boundaries where you might
3 have to take one-sided derivatives. */
4
5#pragma once
6
7namespace FTensor
8{
9 template <class T, int Dim01, int Dim23, char i, char j, char k, char l>
11 {
16
17 public:
19 operator()(const int N1, const int N2, const int N3, const int N4) const
20 {
21 return N3 == N4
22 ? (boundary(N3, 0)
23 ? (*(a.ptr(N1, N2) + 2 * d_ijk(N3))
24 - 2 * *(a.ptr(N1, N2) + d_ijk(N3)) + a(N1, N2))
25 : (boundary(N3, 1)
26 ? (a(N1, N2) - 2 * *(a.ptr(N1, N2) - d_ijk(N3))
27 + *(a.ptr(N1, N2) - 2 * d_ijk(N3)))
28 : (*(a.ptr(N1, N2) + d_ijk(N3)) - 2 * a(N1, N2)
29 + *(a.ptr(N1, N2) - d_ijk(N3)))))
30 * d_xyz(N3) * d_xyz(N3)
31 : (boundary(N3, 0)
32 ? (boundary(N4, 0)
33 ? ((*(a.ptr(N1, N2) + d_ijk(N3) + d_ijk(N4))
34 - *(a.ptr(N1, N2) + d_ijk(N4))
35 - *(a.ptr(N1, N2) + d_ijk(N3)) + a(N1, N2))
36 * d_xyz(N3) * d_xyz(N4))
37 : (boundary(N4, 1)
38 ? ((*(a.ptr(N1, N2) + d_ijk(N3)) - a(N1, N2)
39 - *(a.ptr(N1, N2) + d_ijk(N3) - d_ijk(N4))
40 + *(a.ptr(N1, N2) - d_ijk(N4)))
41 * d_xyz(N3) * d_xyz(N4))
42 : (*(a.ptr(N1, N2) + d_ijk(N3) + d_ijk(N4))
43 - *(a.ptr(N1, N2) + d_ijk(N4))
44 - *(a.ptr(N1, N2) + d_ijk(N3) - d_ijk(N4))
45 + *(a.ptr(N1, N2) - d_ijk(N4)))
46 * d_xyz(N3) * d_xyz(N4) * 0.5))
47 : (boundary(N3, 1)
48 ? (boundary(N4, 0)
49 ? ((*(a.ptr(N1, N2) + d_ijk(N4))
50 - *(a.ptr(N1, N2) - d_ijk(N3) + d_ijk(N4))
51 - a(N1, N2) + *(a.ptr(N1, N2) - d_ijk(N3)))
52 * d_xyz(N3) * d_xyz(N4))
53 : (boundary(N4, 1)
54 ? ((a(N1, N2) - *(a.ptr(N1, N2) - d_ijk(N3))
55 - *(a.ptr(N1, N2) - d_ijk(N4))
56 + *(a.ptr(N1, N2) - d_ijk(N3)
57 - d_ijk(N4)))
58 * d_xyz(N3) * d_xyz(N4))
59 : (*(a.ptr(N1, N2) + d_ijk(N4))
60 - *(a.ptr(N1, N2) - d_ijk(N3) + d_ijk(N4))
61 - *(a.ptr(N1, N2) - d_ijk(N4))
62 + *(a.ptr(N1, N2) - d_ijk(N3)
63 - d_ijk(N4)))
64 * d_xyz(N3) * d_xyz(N4) * 0.5))
65 : (boundary(N4, 0)
66 ? ((*(a.ptr(N1, N2) + d_ijk(N3) + d_ijk(N4))
67 - *(a.ptr(N1, N2) - d_ijk(N3) + d_ijk(N4))
68 - *(a.ptr(N1, N2) + d_ijk(N3))
69 + *(a.ptr(N1, N2) - d_ijk(N3)))
70 * d_xyz(N3) * d_xyz(N4) * 0.5)
71 : (boundary(N4, 1)
72 ? ((*(a.ptr(N1, N2) + d_ijk(N3))
73 - *(a.ptr(N1, N2) - d_ijk(N3))
74 - *(a.ptr(N1, N2) + d_ijk(N3)
75 - d_ijk(N4))
76 + *(a.ptr(N1, N2) - d_ijk(N3)
77 - d_ijk(N4)))
78 * d_xyz(N3) * d_xyz(N4) * 0.5)
79 : ((*(a.ptr(N1, N2) + d_ijk(N3) + d_ijk(N4))
80 - *(a.ptr(N1, N2) - d_ijk(N3)
81 + d_ijk(N4))
82 - *(a.ptr(N1, N2) + d_ijk(N3)
83 - d_ijk(N4))
84 + *(a.ptr(N1, N2) - d_ijk(N3)
85 - d_ijk(N4)))
86 * d_xyz(N3) * d_xyz(N4) * 0.25)))));
87 }
89 const Tensor1<int, Dim23> &D_ijk,
90 const Tensor1<double, Dim23> &D_xyz,
91 const Tensor2<bool, Dim23, 2> &Boundary)
92 : a(A), d_ijk(D_ijk), d_xyz(D_xyz), boundary(Boundary)
93 {}
94 };
95
96 template <class T, int Dim01, int Dim23, char i, char j, char k, char l>
97 const Ddg_Expr<
98 const dd_boundary_Tensor2_symmetric<T, Dim01, Dim23, i, j, k, l>,
99 typename promote<T, double>::V, Dim01, Dim23, i, j, k, l>
101 const Index<i, Dim01> index1, const Index<j, Dim01> index2,
102 const Index<k, Dim23> index3, const Index<l, Dim23> index4,
103 const Tensor1<int, Dim23> &d_ijk,
104 const Tensor1<double, Dim23> &d_xyz,
105 const Tensor2<bool, Dim23, 2> &boundary)
106 {
107 using Tensor_Expr
110 i, j, k, l>(Tensor_Expr(a, d_ijk, d_xyz, boundary));
111 }
112}
static Number< 2 > N2
static Number< 1 > N1
constexpr double a
promote< T, double >::V operator()(const int N1, const int N2, const int N3, const int N4) const
const Tensor2_symmetric< T *, Dim01 > & a
dd_boundary_Tensor2_symmetric(const Tensor2_symmetric< T *, Dim01 > &A, const Tensor1< int, Dim23 > &D_ijk, const Tensor1< double, Dim23 > &D_xyz, const Tensor2< bool, Dim23, 2 > &Boundary)
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
Tensors class implemented by Walter Landry.
Definition FTensor.hpp:51
const Tensor2_symmetric_Expr< const dd_boundary_Tensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd_boundary(const Tensor0< T * > &a, const Index< i, Dim > index3, const Index< j, Dim > index4, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz, const Tensor2< bool, Dim, 2 > &boundary)