FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
2d_p1_unrolled.hpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6#pragma once
7
8#include <kernel/base_header.hpp>
9#include <kernel/shape.hpp>
10#include <kernel/meshopt/base.hpp>
11#include <kernel/meshopt/rumpf_functional.hpp>
12
13namespace FEAT
14{
15 namespace Meshopt
16 {
18
25 template<typename DataType_>
26 class RumpfFunctionalUnrolled<DataType_,
27 Trafo::Standard::Mapping<Geometry::ConformalMesh<Shape::Simplex<2>, 2, DataType_>>> :
28 public RumpfFunctionalBase<DataType_>
29 {
30 public:
32 typedef RumpfFunctionalBase<DataType_> BaseClass;
33
35 typedef DataType_ DataType;
37 typedef Shape::Simplex<2> ShapeType;
39 typedef Trafo::Standard::Mapping<Geometry::ConformalMesh<ShapeType, 2, DataType_>> TrafoType;
41 static constexpr int world_dim = TrafoType::world_dim;
43 static constexpr int shape_dim = ShapeType::dimension;
45 typedef typename Intern::TrafoFE<TrafoType>::Space TrafoSpace;
46
48 typedef Tiny::Matrix<DataType_, Shape::FaceTraits<ShapeType,0>::count, world_dim> Tx;
50 typedef Tiny::Vector<DataType_, Shape::FaceTraits<ShapeType,0>::count*world_dim> Tgradh;
51
53 typedef Tiny::Matrix<DataType_, shape_dim, world_dim> TgradR;
54
56 typedef typename TrafoType::template Evaluator<ShapeType, DataType>::Type TrafoEvaluator;
58 typedef typename TrafoSpace::template Evaluator<TrafoEvaluator>::Type SpaceEvaluator;
59
60 private:
62 const int _exponent_det;
64 const bool _compute_frobenius;
66 const bool _compute_cof;
68 const bool _compute_det;
69
70 public:
74 explicit RumpfFunctionalUnrolled(
75 const DataType fac_frobenius_,
76 const DataType fac_det_,
77 const DataType fac_cof_,
78 const DataType fac_reg_,
79 const int exponent_det_) :
80 BaseClass( fac_frobenius_,
81 fac_det_,
82 fac_det_*(Math::sqrt( Math::sqr(fac_reg_) + DataType(1) )*Math::pow( DataType(1) + Math::sqrt(Math::sqr(fac_reg_) + DataType(1)), DataType(exponent_det_))),
83 fac_cof_,
84 fac_reg_),
85 _exponent_det(exponent_det_),
86 _compute_frobenius( (fac_frobenius_ > DataType(0)) ),
87 _compute_cof( false ),
88 _compute_det( (fac_det_ > 0) )
89 {
90 XASSERTM(exponent_det_ == 1 || exponent_det_ == 2,"exponent_det must be 1 or 2!");
91 XASSERTM(fac_cof_ == DataType(0), "In 2d, the cofactor term is redundant, so set fac_cof == 0.");
92 }
93
99 static String name()
100 {
101 return "RumpfFunctionalUnrolled<"+ShapeType::name()+">";
102 }
103
107 String info() const
108 {
109 const Index pad_width(30);
110 return name() + ":" + BaseClass::info() + String("\nexponent_det").pad_back(pad_width, '.')
111 + String(": ") + stringify(_exponent_det);
112 }
113
117 void eval_fval_grad(DataType& fval, Tx& grad, const TgradR& DOXY(mat_tensor), const TrafoEvaluator& DOXY(trafo_eval), const SpaceEvaluator& DOXY(space_eval), const Tx& x, const DataType& h)
118 {
119
120 fval = DataType(0);
121 grad.format(DataType(0));
122
123 DataType fval_frobenius(0);
124 DataType fval_det(0);
125 DataType fval_rec_det(0);
126
127 if(_compute_frobenius)
128 {
129 fval_frobenius = this->_fac_frobenius*compute_frobenius_part(x,h);
130 add_grad_frobenius_part(grad, x, h);
131 }
132
133 if(_compute_det)
134 {
135 if(_exponent_det == 1)
136 {
137 fval_det = this->_fac_det*compute_det_1_part(x,h);
138 fval_rec_det = this->_fac_rec_det*compute_rec_det_1_part(x,h);
139 add_grad_det_1_part(grad, x, h);
140 add_grad_rec_det_1_part(grad, x, h);
141 }
142 else
143 {
144 fval_det = this->_fac_det*compute_det_2_part(x,h);
145 fval_rec_det = this->_fac_rec_det*compute_rec_det_2_part(x,h);
146 add_grad_det_2_part(grad, x, h);
147 add_grad_rec_det_2_part(grad, x, h);
148 }
149 }
150
151 fval = fval_frobenius + fval_det + fval_rec_det;
152
153 }
154
155 void eval_fval_cellwise(DataType& fval, const TgradR& DOXY(mat_tensor), const TrafoEvaluator& DOXY(trafo_eval), const SpaceEvaluator& DOXY(space_eval), const Tx& x, const DataType& h, DataType& fval_frobenius, DataType& fval_cof, DataType& fval_det)
156 {
157 fval = DataType(0);
158 fval_frobenius = DataType(0);
159 fval_cof = DataType(0);
160 fval_det = DataType(0);
161
162 if(_compute_frobenius)
163 {
164 fval_frobenius = this->_fac_frobenius*compute_frobenius_part(x,h);
165 }
166
167 if(_compute_det)
168 {
169 if(_exponent_det == 1)
170 {
171 fval_det = this->_fac_det*compute_det_1_part(x,h);
172 fval_det += this->_fac_rec_det*compute_rec_det_1_part(x,h);
173 }
174 else
175 {
176 fval_det = this->_fac_det*compute_det_2_part(x,h);
177 fval_det += this->_fac_rec_det*compute_rec_det_2_part(x,h);
178 }
179 }
180
181 fval = fval_frobenius + fval_det;
182
183 }
184
188 void NOINLINE add_grad_h_part(Tx& grad, const TgradR& DOXY(mat_tensor), const TrafoEvaluator& DOXY(trafo_eval), const SpaceEvaluator& DOXY(space_eval), const Tx& x, const DataType& h, const Tgradh& grad_h)
189 {
190
191 DataType frobenius_der_h(0);
192 DataType det_der_h(0);
193 DataType rec_det_der_h(0);
194
195 if(_compute_frobenius)
196 {
197 frobenius_der_h = this->_fac_frobenius*compute_frobenius_der_h_part(x,h);
198 }
199
200 if(_compute_det)
201 {
202 if(_exponent_det == 1)
203 {
204 det_der_h = this->_fac_det * compute_det_1_der_h_part(x,h);
205 rec_det_der_h = this->_fac_rec_det * compute_rec_det_1_der_h_part(x,h);
206 }
207 else
208 {
209 det_der_h = this->_fac_det * compute_det_2_der_h_part(x,h);
210 rec_det_der_h = this->_fac_rec_det * compute_rec_det_2_der_h_part(x,h);
211 }
212
213 }
214
215 DataType der_h(frobenius_der_h + det_der_h + rec_det_der_h);
216
217 for(int i(0); i < Tx::m; ++i)
218 {
219 for(int d(0); d < Tx::n; ++d)
220 {
221 grad(i,d) += der_h*grad_h(i*Tx::n + d);
222 }
223 }
224 } // add_grad_h_part
225
226
230 DataType NOINLINE compute_frobenius_part(const Tx& x, const DataType& h)
231 {
232 DataType fval_frobenius_part(0);
233 fval_frobenius_part = DataType(4)/DataType(9)*Math::pow(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)),DataType(2))/(h*h*h*h);
234 return fval_frobenius_part;
235 }
236
240 DataType NOINLINE compute_det_1_part( const Tx& x, const DataType& h)
241 {
242 DataType fval_det_1_part(0);
243 fval_det_1_part = DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h);
244 return fval_det_1_part;
245 }
249 DataType NOINLINE compute_det_2_part( const Tx& x, const DataType& h)
250 {
251 DataType fval_det_2_part(0);
252 fval_det_2_part = DataType(4)/DataType(3)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h);
253 return fval_det_2_part;
254 }
255
259 DataType NOINLINE compute_rec_det_1_part( const Tx& x, const DataType& h)
260 {
261 DataType fval_rec_det_1_part(0);
262 fval_rec_det_1_part = DataType(1)/(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3));
263 return fval_rec_det_1_part;
264 }
265
269 DataType NOINLINE compute_rec_det_2_part( const Tx& x, const DataType& h)
270 {
271 DataType fval_rec_det_2_part(0);
272 fval_rec_det_2_part = Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2));
273 return fval_rec_det_2_part;
274 }
278 DataType NOINLINE compute_frobenius_der_h_part(const Tx& x, const DataType& h)
279 {
280 DataType frobenius_der_h_part(0);
281 frobenius_der_h_part = DataType(16)/DataType(3)*(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)))/(h*h*h)-DataType(16)/DataType(9)*Math::pow(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)),DataType(2))/(h*h*h*h*h);
282 return frobenius_der_h_part;
283 }
284
288 DataType NOINLINE compute_det_1_der_h_part( const Tx& x, const DataType& h)
289 {
290 DataType det_1_der_h_part(0);
291 det_1_der_h_part = -DataType(4)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h);
292 return det_1_der_h_part;
293 }
297 DataType NOINLINE compute_det_2_der_h_part( const Tx& x, const DataType& h)
298 {
299 DataType det_2_der_h_part(0);
300 det_2_der_h_part = -DataType(16)/DataType(3)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h*h);
301 return det_2_der_h_part;
302 }
303
307 DataType NOINLINE compute_rec_det_1_der_h_part( const Tx& x, const DataType& h)
308 {
309 DataType rec_det_1_der_h_part(0);
310 rec_det_1_der_h_part = -Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2))*(-DataType(4)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h)-DataType(8)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h*h));
311 return rec_det_1_der_h_part;
312 }
313
317 DataType NOINLINE compute_rec_det_2_der_h_part( const Tx& x, const DataType& h)
318 {
319 DataType rec_det_2_der_h_part(0);
320 rec_det_2_der_h_part = -DataType(2)*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(3))*(-DataType(4)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h)-DataType(8)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h*h));
321 return rec_det_2_der_h_part;
322 }
323
327 void NOINLINE add_grad_frobenius_part(Tx& grad_frobenius_part, const Tx& x, const DataType& h)
328 {
329 grad_frobenius_part(0,0) += DataType(8)/DataType(9)*this->_fac_frobenius*(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)))/(h*h*h*h)*(-DataType(4)*x(0,0)+DataType(2)*x(1,0)+DataType(2)*x(2,0));
330
331 grad_frobenius_part(0,1) += DataType(8)/DataType(9)*this->_fac_frobenius*(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)))/(h*h*h*h)*(-DataType(4)*x(0,1)+DataType(2)*x(1,1)+DataType(2)*x(2,1));
332
333 grad_frobenius_part(1,0) += DataType(8)/DataType(9)*this->_fac_frobenius*(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)))/(h*h*h*h)*(DataType(2)*x(0,0)-DataType(4)*x(1,0)+DataType(2)*x(2,0));
334
335 grad_frobenius_part(1,1) += DataType(8)/DataType(9)*this->_fac_frobenius*(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)))/(h*h*h*h)*(DataType(2)*x(0,1)-DataType(4)*x(1,1)+DataType(2)*x(2,1));
336
337 grad_frobenius_part(2,0) += DataType(8)/DataType(9)*this->_fac_frobenius*(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)))/(h*h*h*h)*(DataType(2)*x(0,0)+DataType(2)*x(1,0)-DataType(4)*x(2,0));
338
339 grad_frobenius_part(2,1) += DataType(8)/DataType(9)*this->_fac_frobenius*(DataType(3)*h*h-DataType(2)*Math::sqr(x(0,0))+DataType(2)*x(0,0)*x(1,0)+DataType(2)*x(0,0)*x(2,0)-DataType(2)*Math::sqr(x(0,1))+DataType(2)*x(0,1)*x(1,1)+DataType(2)*x(0,1)*x(2,1)-DataType(2)*Math::sqr(x(1,0))+DataType(2)*x(1,0)*x(2,0)-DataType(2)*Math::sqr(x(1,1))+DataType(2)*x(1,1)*x(2,1)-DataType(2)*Math::sqr(x(2,0))-DataType(2)*Math::sqr(x(2,1)))/(h*h*h*h)*(DataType(2)*x(0,1)+DataType(2)*x(1,1)-DataType(4)*x(2,1));
340
341 }
342
343
347 void NOINLINE add_grad_det_1_part(Tx& grad_det_1_part, const Tx& x, const DataType& h)
348 {
349 grad_det_1_part(0,0) += DataType(2)/DataType(3)*this->_fac_det*Math::sqrt(DataType(3))*(x(1,1)-x(2,1))/Math::sqr(h);
350
351 grad_det_1_part(0,1) += DataType(2)/DataType(3)*this->_fac_det*Math::sqrt(DataType(3))*(-x(1,0)+x(2,0))/Math::sqr(h);
352
353 grad_det_1_part(1,0) += DataType(2)/DataType(3)*this->_fac_det*Math::sqrt(DataType(3))*(-x(0,1)+x(2,1))/Math::sqr(h);
354
355 grad_det_1_part(1,1) += DataType(2)/DataType(3)*this->_fac_det*Math::sqrt(DataType(3))*(x(0,0)-x(2,0))/Math::sqr(h);
356
357 grad_det_1_part(2,0) += DataType(2)/DataType(3)*this->_fac_det*Math::sqrt(DataType(3))*(x(0,1)-x(1,1))/Math::sqr(h);
358
359 grad_det_1_part(2,1) += DataType(2)/DataType(3)*this->_fac_det*Math::sqrt(DataType(3))*(-x(0,0)+x(1,0))/Math::sqr(h);
360 }
361
365 void NOINLINE add_grad_det_2_part(Tx& grad_det_2_part, const Tx& x, const DataType& h)
366 {
367 grad_det_2_part(0,0) += DataType(8)/DataType(3)*this->_fac_det*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(1,1)-x(2,1));
368
369 grad_det_2_part(0,1) += DataType(8)/DataType(3)*this->_fac_det*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(1,0)+x(2,0));
370
371 grad_det_2_part(1,0) += DataType(8)/DataType(3)*this->_fac_det*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(0,1)+x(2,1));
372
373 grad_det_2_part(1,1) += DataType(8)/DataType(3)*this->_fac_det*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(0,0)-x(2,0));
374
375 grad_det_2_part(2,0) += DataType(8)/DataType(3)*this->_fac_det*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(0,1)-x(1,1));
376
377 grad_det_2_part(2,1) += DataType(8)/DataType(3)*this->_fac_det*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(0,0)+x(1,0));
378 }
379
383 void NOINLINE add_grad_rec_det_1_part(Tx& grad_rec_det_1_part, const Tx& x, const DataType& h)
384 {
385 grad_rec_det_1_part(0,0) += -this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(1,1)-x(2,1))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(1,1)-x(2,1)));
386
387 grad_rec_det_1_part(0,1) += -this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(-x(1,0)+x(2,0))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(1,0)+x(2,0)));
388
389 grad_rec_det_1_part(1,0) += -this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(-x(0,1)+x(2,1))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(0,1)+x(2,1)));
390
391 grad_rec_det_1_part(1,1) += -this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)-x(2,0))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(0,0)-x(2,0)));
392
393 grad_rec_det_1_part(2,0) += -this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,1)-x(1,1))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(0,1)-x(1,1)));
394
395 grad_rec_det_1_part(2,1) += -this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(2))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(-x(0,0)+x(1,0))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(0,0)+x(1,0)));
396 }
397
401 void NOINLINE add_grad_rec_det_2_part(Tx& grad_rec_det_2_part, const Tx& x, const DataType& h)
402 {
403 grad_rec_det_2_part(0,0) += -DataType(2)*this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(3))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(1,1)-x(2,1))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(1,1)-x(2,1)));
404
405 grad_rec_det_2_part(0,1) += -DataType(2)*this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(3))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(-x(1,0)+x(2,0))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(1,0)+x(2,0)));
406
407 grad_rec_det_2_part(1,0) += -DataType(2)*this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(3))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(-x(0,1)+x(2,1))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(0,1)+x(2,1)));
408
409 grad_rec_det_2_part(1,1) += -DataType(2)*this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(3))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)-x(2,0))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(0,0)-x(2,0)));
410
411 grad_rec_det_2_part(2,0) += -DataType(2)*this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(3))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,1)-x(1,1))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(x(0,1)-x(1,1)));
412
413 grad_rec_det_2_part(2,1) += -DataType(2)*this->_fac_rec_det*Math::pow(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/Math::sqr(h)+Math::sqrt(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h))/DataType(3),-DataType(3))*(DataType(2)/DataType(3)*Math::sqrt(DataType(3))*(-x(0,0)+x(1,0))/Math::sqr(h)+DataType(4)*Math::pow(DataType(9)*this->_fac_reg*this->_fac_reg+DataType(12)*Math::pow(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0),DataType(2))/(h*h*h*h),-DataType(1)/DataType(2))*(x(0,0)*x(1,1)-x(0,0)*x(2,1)-x(1,0)*x(0,1)+x(2,0)*x(0,1)+x(1,0)*x(2,1)-x(1,1)*x(2,0))/(h*h*h*h)*(-x(0,0)+x(1,0)));
414 }
415
416
417 }; // class RumpfFunctionalUnrolled
418
419
420 //#ifdef FEAT_EICKT
421 // extern template class RumpfFunctionalUnrolled<double,
422 // Trafo::Standard::Mapping<Geometry::ConformalMesh<Shape::Simplex<2>, 2, double>>>;
423 //#endif
425 } // namespace Meshopt
426} // namespace FEAT
XASSERTM
#define XASSERTM(expr, msg)
Assertion macro definition with custom message.
Definition: assertion.hpp:263
base_header.hpp
FEAT Kernel base header.
FEAT::Math::sqrt
T_ sqrt(T_ x)
Returns the square-root of a value.
Definition: math.hpp:300
FEAT::Math::pow
T_ pow(T_ x, T_ y)
Returns x raised to the power of y.
Definition: math.hpp:643
FEAT::Math::sqr
T_ sqr(T_ x)
Returns the square of a value.
Definition: math.hpp:95
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::stringify
String stringify(const T_ &item)
Converts an item into a String.
Definition: string.hpp:993
FEAT::SpaceTags::grad
@ grad
specifies whether the space should supply basis function gradients
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122