Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #ifndef LIBMESH_FDM_GRADIENT_H
19 : #define LIBMESH_FDM_GRADIENT_H
20 :
21 : // libMesh includes
22 : #include "libmesh/fem_function_base.h"
23 : #include "libmesh/libmesh_common.h"
24 : #include "libmesh/tensor_tools.h"
25 :
26 : namespace libMesh
27 : {
28 :
29 : template <typename GradType>
30 0 : class FDMGradient : public FEMFunctionBase<GradType>
31 : {
32 : public:
33 : typedef typename TensorTools::DecrementRank<GradType>::type ValType;
34 :
35 0 : FDMGradient (FEMFunctionBase<ValType> & value_func,
36 : Real eps) :
37 0 : _val_func(value_func.clone()), _eps(eps)
38 0 : {}
39 :
40 0 : virtual void init_context (const FEMContext & c) override
41 0 : { _val_func->init_context(c); }
42 :
43 0 : virtual std::unique_ptr<FEMFunctionBase<GradType>> clone () const override
44 0 : { return std::make_unique<FDMGradient<GradType>>(*_val_func, _eps); }
45 :
46 0 : virtual GradType operator() (const FEMContext & c,
47 : const Point & p,
48 : const Real time = 0.) override
49 : {
50 0 : GradType g;
51 :
52 0 : auto & val = *_val_func;
53 :
54 0 : Real one_over_dim = Real(0.5) / _eps;
55 :
56 0 : g(0) = (val(c, p+Point(_eps), time) -
57 0 : val(c, p+Point(-_eps), time)) * one_over_dim;
58 : #if LIBMESH_DIM > 1
59 0 : g(1) = (val(c, p+Point(0,_eps), time) -
60 0 : val(c, p+Point(0,-_eps), time)) * one_over_dim;
61 : #endif
62 : #if LIBMESH_DIM > 2
63 0 : g(2) = (val(c, p+Point(0,0,_eps), time) -
64 0 : val(c, p+Point(0,0,-_eps), time)) * one_over_dim;
65 : #endif
66 :
67 0 : return g;
68 : }
69 :
70 0 : virtual void operator() (const FEMContext & c,
71 : const Point & p,
72 : const Real time,
73 : DenseVector<GradType> & output) override
74 : {
75 0 : auto sz = output.size();
76 0 : DenseVector<ValType> v(sz);
77 :
78 0 : auto & val = *_val_func;
79 :
80 0 : val(c, p+Point(_eps), time, v);
81 0 : for (auto i : make_range(sz))
82 0 : output(i)(0) = v(i);
83 :
84 0 : val(c, p+Point(-_eps), time, v);
85 0 : for (auto i : make_range(sz))
86 : {
87 0 : output(i)(0) -= v(i);
88 0 : output(i)(0) /= 2;
89 : }
90 :
91 : #if LIBMESH_DIM > 1
92 0 : val(c, p+Point(0,_eps), time, v);
93 0 : for (auto i : make_range(sz))
94 0 : output(i)(1) = v(i);
95 :
96 0 : val(c, p+Point(0,-_eps), time, v);
97 0 : for (auto i : make_range(sz))
98 : {
99 0 : output(i)(1) -= v(i);
100 0 : output(i)(1) /= 2;
101 : }
102 : #endif
103 : #if LIBMESH_DIM > 2
104 0 : val(c, p+Point(0,0,_eps), time, v);
105 0 : for (auto i : make_range(sz))
106 0 : output(i)(2) = v(i);
107 :
108 0 : val(c, p+Point(0,0,-_eps), time, v);
109 0 : for (auto i : make_range(sz))
110 : {
111 0 : output(i)(2) -= v(i);
112 0 : output(i)(2) /= 2;
113 : }
114 : #endif
115 0 : }
116 :
117 :
118 0 : virtual GradType component (const FEMContext & c,
119 : unsigned int i,
120 : const Point & p,
121 : Real time) override
122 : {
123 0 : GradType g;
124 :
125 0 : auto & val = *_val_func;
126 :
127 0 : Real one_over_dim = Real(0.5) / _eps;
128 :
129 0 : g(0) = (val.component(c, i, p+Point(_eps), time) -
130 0 : val.component(c, i, p+Point(-_eps), time)) * one_over_dim;
131 : #if LIBMESH_DIM > 1
132 0 : g(1) = (val.component(c, i, p+Point(0,_eps), time) -
133 0 : val.component(c, i, p+Point(0,-_eps), time)) * one_over_dim;
134 : #endif
135 : #if LIBMESH_DIM > 2
136 0 : g(2) = (val.component(c, i, p+Point(0,0,_eps), time) -
137 0 : val.component(c, i, p+Point(0,0,-_eps), time)) * one_over_dim;
138 : #endif
139 :
140 0 : return g;
141 : }
142 :
143 : private:
144 :
145 : std::unique_ptr<FEMFunctionBase<ValType>> _val_func;
146 :
147 : Real _eps;
148 : };
149 :
150 :
151 : } // namespace libMesh
152 :
153 : #endif // LIBMESH_FDM_GRADIENT_H
|
