FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
node_functional.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// includes, FEAT
9#include <kernel/base_header.hpp>
10#include <kernel/space/node_functional_base.hpp>
11#include <kernel/cubature/dynamic_factory.hpp>
12#include <kernel/space/bernstein2/scalar_basis.hpp>
13#include <kernel/util/tiny_algebra.hpp>
14
15
16namespace FEAT
17{
18 namespace Space
19 {
20 namespace Bernstein2
21 {
22 template<
23 typename Space_,
24 typename ShapeType_,
25 typename DataType_>
26 class NodeFunctional :
27 public NodeFunctionalNull<Space_, DataType_>
28 {
29 public:
30 explicit NodeFunctional(const Space_& space) :
31 NodeFunctionalNull<Space_, DataType_>(space)
32 {
33 }
34 };
35
36 template<
37 typename Space_,
38 typename DataType_>
39 class NodeFunctional<Space_, Shape::Vertex, DataType_> :
40 public NodeFunctionalBase<Space_, DataType_>
41 {
42 public:
43 typedef NodeFunctionalBase<Space_, DataType_> BaseClass;
44 static constexpr int max_assigned_dofs = 1;
45
46 template<typename Function_>
47 struct Value
48 {
49 typedef typename Analytic::EvalTraits<DataType_, Function_>::ValueType Type;
50 };
51
52 protected:
53 typedef typename Space_::TrafoType TrafoType;
54 typedef typename TrafoType::template Evaluator<Shape::Vertex, DataType_>::Type TrafoEvalType;
55 typedef typename TrafoEvalType::EvalTraits TrafoEvalTraits;
56 typedef typename TrafoEvalTraits::DataType DataType;
57 typedef typename TrafoEvalTraits::DomainPointType DomainPointType;
58
59 // declare trafo evaluation data
60 typedef typename TrafoEvalType::template ConfigTraits<TrafoTags::img_point>::EvalDataType TrafoEvalData;
61
62 TrafoEvalType _trafo_eval;
63
64 public:
65 explicit NodeFunctional(const Space_& space) :
66 BaseClass(space),
67 _trafo_eval(space.get_trafo())
68 {
69 }
70
71 void prepare(Index cell_index)
72 {
73 BaseClass::prepare(cell_index);
74 _trafo_eval.prepare(cell_index);
75 }
76
77 void finish()
78 {
79 _trafo_eval.finish();
80 BaseClass::finish();
81 }
82
83 int get_num_assigned_dofs() const
84 {
85 return max_assigned_dofs;
86 }
87
88 template<typename NodeData_, typename Function_>
89 void operator()(NodeData_& node_data, const Function_& function) const
90 {
91 static_assert(std::is_base_of<Analytic::Function, Function_>::value, "invalid function object");
92 static_assert(Function_::can_value, "function cannot compute values");
93
94 // declare our evaluation traits
95 typedef Analytic::EvalTraits<DataType_, Function_> FuncEvalTraits;
96
97 // declare function evaluator
98 typename Function_::template Evaluator<FuncEvalTraits> func_eval(function);
99
100 // compute trafo data
101 DomainPointType dom_point(DataType(0));
102 TrafoEvalData trafo_data;
103 _trafo_eval(trafo_data, dom_point);
104
105 // evaluate function
106 node_data[0] = func_eval.value(trafo_data.img_point);
107 }
108 };
109
110 template<
111 typename Space_,
112 typename DataType_>
113 class NodeFunctional<Space_, Shape::Hypercube<1>, DataType_> :
114 public NodeFunctionalBase<Space_, DataType_>
115 {
116 public:
117 typedef NodeFunctionalBase<Space_, DataType_> BaseClass;
118 static constexpr int max_assigned_dofs = 1;
119
120
121 template<typename Function_>
122 struct Value
123 {
124 typedef typename Analytic::EvalTraits<DataType_, Function_>::ValueType Type;
125 };
126
127 protected:
128 typedef typename Space_::TrafoType TrafoType;
129 typedef Shape::Hypercube<1> ShapeType;
130 typedef typename TrafoType::template Evaluator<ShapeType, DataType_>::Type TrafoEvalType;
131 typedef typename TrafoEvalType::EvalTraits TrafoEvalTraits;
132 typedef typename TrafoEvalTraits::DataType DataType;
133 typedef typename TrafoEvalTraits::DomainPointType DomainPointType;
134
135 static constexpr int image_dim = TrafoEvalTraits::image_dim;
136
137 // The coef Matrix gets calculatet in prepare.
138 // Every row represents the coefficients of a basis function.
139 // The i-th basis functions can be evaluated in (x,y) via Intern::base_q(x,y,coef[i])
140 Tiny::Matrix<DataType_, 3, 3> coef;
141
142 static constexpr TrafoTags trafo_tags = TrafoTags::img_point | TrafoTags::jac_det;
143 typedef typename TrafoEvalType::template ConfigTraits<trafo_tags>::EvalDataType TrafoEvalData;
144
145 typedef Cubature::Rule<ShapeType, DataType_, DataType_, DomainPointType> CubRuleType;
146
147 TrafoEvalType _trafo_eval;
148 CubRuleType _cub_rule;
149
150 public:
151 explicit NodeFunctional(const Space_& space) :
152 BaseClass(space),
153 _trafo_eval(space.get_trafo()),
154 _cub_rule(Cubature::ctor_factory, Cubature::DynamicFactory("gauss-legendre:3"))
155 {
156 coef.format();
157 }
158
159 // prepare() is called for every FE.
160 // In prepare the coefficientmatrix is calculated for every FE because otherwise the isoparametic trafo doesnt work.
161 void prepare(Index cell_index)
162 {
163 coef.format();
164 BaseClass::prepare(cell_index);
165 _trafo_eval.prepare(cell_index);
166
167 // declare trafo evaluation data
168 TrafoEvalData trafo_data;
169
170
171 Tiny::Matrix<DataType_, 3, 3> temp(DataType_(0));
172 DataType mean(DataType(0));
173
174 for (int i(0); i < _cub_rule.get_num_points(); ++i)
175 {
176 // get current cubature point
177 _trafo_eval(trafo_data, _cub_rule.get_point(i));
178
179 // calculate p0,p1,p2,q0,q1,q2 at the i-th cubature point[x] and temporary store them in temp like this:
180 //temp=[p0(x),q0(x),0;
181 // p1(x),q1(x),0;
182 // p2(x),q2(x),0;]
183 for (int j = 0; j < 3; ++j)
184 {
185 temp(j, 0) = Intern::p(j, trafo_data.dom_point[0]);
186 temp(j, 1) = Intern::q(j, trafo_data.dom_point[0]);
187 }
188
189 DataType_ mult = _cub_rule.get_weight(i) * trafo_data.jac_det;
190 mean += mult;
191 for (int px = 0; px < 3; ++px)
192 {
193 for (int qx = 0; qx < 3; ++qx)
194 {
195 coef(px, qx) +=
196 (mult * temp(px, 0) * temp(qx, 1));
197 } // end for qx
198 } // end for px
199 temp.format();
200 } // end cubature
201 // scale coefficient matrix
202 coef *= DataType(2) / mean;
203 temp.set_inverse(coef);
204 coef.set_transpose(temp);
205 } // end prepare
206
207 void finish()
208 {
209 _trafo_eval.finish();
210 BaseClass::finish();
211 }
212
213 int get_num_assigned_dofs() const
214 {
215 return max_assigned_dofs;
216 }
217
218 template<typename NodeData_, typename Function_>
219 void operator()(NodeData_& node_data, const Function_& function) const
220 {
221 static_assert(std::is_base_of<Analytic::Function, Function_>::value, "invalid function object");
222 static_assert(Function_::can_value, "function cannot compute values");
223
224 // declare our evaluation traits
225 typedef Analytic::EvalTraits<DataType_, Function_> FuncEvalTraits;
226
227 // declare function evaluator
228 typename Function_::template Evaluator<FuncEvalTraits> func_eval(function);
229
230 // compute trafo data
231 typename FuncEvalTraits::ValueType value(DataType_(0));
232 DataType_ mean(DataType_(0));
233 TrafoEvalData trafo_data;
234
235 // integrate over facet with cubature factory
236 for (int i(0); i < _cub_rule.get_num_points(); ++i)
237 {
238 _trafo_eval(trafo_data, _cub_rule.get_point(i));
239 value += (_cub_rule.get_weight(i) * trafo_data.jac_det) * func_eval.value(trafo_data.img_point)
240 *Intern::base_q(trafo_data.dom_point[0], coef[2]);
241 mean += _cub_rule.get_weight(i) *trafo_data.jac_det;
242 }
243
244 // set node_data to the integral mean
245 node_data[0] = (DataType_(2) / mean)* value;
246 }
247 };
248
249 template<
250 typename Space_,
251 typename DataType_>
252 class NodeFunctional<Space_, Shape::Hypercube<2>, DataType_> :
253 public NodeFunctionalBase<Space_, DataType_>
254 {
255 public:
256 typedef NodeFunctionalBase<Space_, DataType_> BaseClass;
257 static constexpr int max_assigned_dofs = 1;
258
259 template<typename Function_>
260 struct Value
261 {
262 typedef typename Analytic::EvalTraits<DataType_, Function_>::ValueType Type;
263 };
264
265 protected:
266 typedef typename Space_::TrafoType TrafoType;
267 typedef Shape::Hypercube<2> ShapeType;
268 typedef typename TrafoType::template Evaluator<ShapeType, DataType_>::Type TrafoEvalType;
269 typedef typename TrafoEvalType::EvalTraits TrafoEvalTraits;
270 typedef typename TrafoEvalTraits::DataType DataType;
271 typedef typename TrafoEvalTraits::DomainPointType DomainPointType;
272
273 static constexpr int image_dim = TrafoEvalTraits::image_dim;
274
275 // The coef Matrix gets calculatet in prepare.
276 // Every row represents the coefficients of a basis function.
277 // The i-th basis functions can be evaluated in (x,y) via Intern::base_q(x,y,coef[i])
278 Tiny::Matrix<DataType_, 9, 9> coef;
279
280
281 // declare trafo evaluation data
282 static constexpr TrafoTags trafo_tags = TrafoTags::img_point | TrafoTags::jac_det;
283 typedef typename TrafoEvalType::template ConfigTraits<trafo_tags>::EvalDataType TrafoEvalData;
284
285 typedef Cubature::Rule<ShapeType, DataType_, DataType_, DomainPointType> CubRuleType;
286
287 TrafoEvalType _trafo_eval;
288 CubRuleType _cub_rule;
289
290 public:
291 explicit NodeFunctional(const Space_& space) :
292 BaseClass(space),
293 _trafo_eval(space.get_trafo()),
294 _cub_rule(Cubature::ctor_factory, Cubature::DynamicFactory("gauss-legendre:3"))
295 {
296 coef.format();
297 }
298
299 // prepare() is called for every FE.
300 // In prepare the coefficientmatrix is calculated for every FE because of the changing jacobi determinant.
301 void prepare(Index cell_index)
302 {
303 coef.format();
304 BaseClass::prepare(cell_index);
305 _trafo_eval.prepare(cell_index);
306
307 // declare trafo evaluation data
308 TrafoEvalData trafo_data;
309
310 Tiny::Matrix<DataType_, 9, 9> temp(DataType_(0));
311 DataType mean(DataType(0));
312 for (int i(0); i < _cub_rule.get_num_points(); ++i)
313 {
314 // get current cubature point
315 _trafo_eval(trafo_data, _cub_rule.get_point(i));
316
317 // calculate p0,p1,p2,q0,q1,q2 at the i-th cubature point[x,y] and temporary store them in temp like this:
318 //temp=[p0(x),p0(y),q0(x),q0(y),0,..0;
319 // p1(x),p1(y),q1(x),q1(y),0,..0;
320 // p2(x),p2(y),q2(x),q2(y),0,..0;
321 // 0,......,0]
322 for (int j = 0; j < 3; ++j)
323 {
324 temp(j, 0) = Intern::p(j, trafo_data.dom_point[0]);
325 temp(j, 1) = Intern::p(j, trafo_data.dom_point[1]);
326 temp(j, 2) = Intern::q(j, trafo_data.dom_point[0]);
327 temp(j, 3) = Intern::q(j, trafo_data.dom_point[1]);
328 }
329
330 DataType_ mult = _cub_rule.get_weight(i) * trafo_data.jac_det;
331 mean += mult;
332 for (int px = 0; px < 3; ++px)
333 {
334 for (int py = 0; py < 3; ++py)
335 {
336 for (int qx = 0; qx < 3; ++qx)
337 {
338 for (int qy = 0; qy < 3; ++qy)
339 {
340 coef(px * 3 + py, qx * 3 + qy) +=
341 (mult * temp(px, 0) * temp(py, 1) * temp(qx, 2) * temp(qy, 3));
342 } // end for qy
343 } // end for qx
344 } // end for p_y
345 } // end for px
346 temp.format();
347 } // end cubature
348 // scale coefficient matrix
349 coef *= DataType(4) / mean;
350 temp.set_inverse(coef);
351 coef.set_transpose(temp);
352 } // end prepare
353
354 void finish()
355 {
356 _trafo_eval.finish();
357 BaseClass::finish();
358 }
359
360 int get_num_assigned_dofs() const
361 {
362 return max_assigned_dofs;
363 }
364
365 template<typename NodeData_, typename Function_>
366 void operator()(NodeData_& node_data, const Function_& function) const
367 {
368 static_assert(std::is_base_of<Analytic::Function, Function_>::value, "invalid function object");
369 static_assert(Function_::can_value, "function cannot compute values");
370
371 // declare our evaluation traits
372 typedef Analytic::EvalTraits<DataType_, Function_> FuncEvalTraits;
373
374 // declare function evaluator
375 typename Function_::template Evaluator<FuncEvalTraits> func_eval(function);
376
377 // compute trafo data
378 typename FuncEvalTraits::ValueType value(DataType_(0));
379 DataType_ mean(DataType_(0));
380 TrafoEvalData trafo_data;
381
382
383 // integrate over facet
384 for (int i(0); i < _cub_rule.get_num_points(); ++i)
385 {
386 _trafo_eval(trafo_data, _cub_rule.get_point(i));
387 value += (_cub_rule.get_weight(i) * trafo_data.jac_det) * func_eval.value(trafo_data.img_point)
388 * Intern::base_q(trafo_data.dom_point[0], trafo_data.dom_point[1], coef[8]);
389 mean += _cub_rule.get_weight(i) * trafo_data.jac_det;
390 }
391
392 // set node_date to integral mean
393 node_data[0] = (DataType_(4) / mean) * value;
394 }
395 };
396
397 template<
398 typename Space_,
399 typename DataType_>
400 class NodeFunctional<Space_, Shape::Hypercube<3>, DataType_> :
401 public NodeFunctionalBase<Space_, DataType_>
402 {
403 public:
404 typedef NodeFunctionalBase<Space_, DataType_> BaseClass;
405 static constexpr int max_assigned_dofs = 1;
406
407 template<typename Function_>
408 struct Value
409 {
410 typedef typename Analytic::EvalTraits<DataType_, Function_>::ValueType Type;
411 };
412
413 protected:
414 typedef typename Space_::TrafoType TrafoType;
415 typedef Shape::Hypercube<3> ShapeType;
416 typedef typename TrafoType::template Evaluator<ShapeType, DataType_>::Type TrafoEvalType;
417 typedef typename TrafoEvalType::EvalTraits TrafoEvalTraits;
418 typedef typename TrafoEvalTraits::DataType DataType;
419 typedef typename TrafoEvalTraits::DomainPointType DomainPointType;
420
421 static constexpr int image_dim = TrafoEvalTraits::image_dim;
422
423 // The coef Matrix gets calculatet in prepare.
424 // Every row represents the coefficients of a basis function.
425 // The i-th basis functions can be evaluated in (x,y) via Intern::base_q(x,y,coef[i])
426 Tiny::Matrix<DataType_, 27, 27> coef;
427
428
429 // declare trafo evaluation data
430 static constexpr TrafoTags trafo_tags = TrafoTags::img_point | TrafoTags::jac_det;
431 typedef typename TrafoEvalType::template ConfigTraits<trafo_tags>::EvalDataType TrafoEvalData;
432
433 typedef Cubature::Rule<ShapeType, DataType_, DataType_, DomainPointType> CubRuleType;
434
435 TrafoEvalType _trafo_eval;
436 CubRuleType _cub_rule;
437
438 public:
439 explicit NodeFunctional(const Space_& space) :
440 BaseClass(space),
441 _trafo_eval(space.get_trafo()),
442 _cub_rule(Cubature::ctor_factory, Cubature::DynamicFactory("gauss-legendre:3"))
443 {
444 coef.format();
445 }
446
447 // prepare() is called for every FE.
448 // In prepare the coefficientmatrix is calculated for every FE because of the changing jacobi determinant.
449 void prepare(Index cell_index)
450 {
451 coef.format();
452 BaseClass::prepare(cell_index);
453 _trafo_eval.prepare(cell_index);
454
455 // declare trafo evaluation data
456 TrafoEvalData trafo_data;
457
458 Tiny::Matrix<DataType_, 27, 27> temp(DataType_(0));
459 DataType mean(DataType(0));
460 for (int i(0); i < _cub_rule.get_num_points(); ++i)
461 {
462 // get current cubature point
463 _trafo_eval(trafo_data, _cub_rule.get_point(i));
464
465 // calculate p0,p1,p2,q0,q1,q2 at the i-th cubature point[x,y,z] and temporary store them in temp like this:
466 //temp=[p0(x),p0(y),p0(z),q0(x),q0(y),q0(z),0,..0;
467 // p1(x),p1(y),p1(z),q1(x),q1(y),q1(z),0,..0;
468 // p2(x),p2(y),p2(z),q2(x),q2(y),q2(z),0,..0;
469 // 0,......,0]
470 for (int j = 0; j < 3; ++j)
471 {
472 temp(j, 0) = Intern::p(j, trafo_data.dom_point[0]);
473 temp(j, 1) = Intern::p(j, trafo_data.dom_point[1]);
474 temp(j, 2) = Intern::p(j, trafo_data.dom_point[2]);
475 temp(j, 3) = Intern::q(j, trafo_data.dom_point[0]);
476 temp(j, 4) = Intern::q(j, trafo_data.dom_point[1]);
477 temp(j, 5) = Intern::q(j, trafo_data.dom_point[2]);
478 }
479
480 DataType_ mult = _cub_rule.get_weight(i) * trafo_data.jac_det;
481 mean += mult;
482 for (int pz = 0; pz < 3; ++pz)
483 {
484 for (int px = 0; px < 3; ++px)
485 {
486 for (int py = 0; py < 3; ++py)
487 {
488 for (int qz = 0; qz < 3; ++qz)
489 {
490 for (int qx = 0; qx < 3; ++qx)
491 {
492 for (int qy = 0; qy < 3; ++qy)
493 {
494 coef(9*pz+px * 3 + py, 9*qz+qx * 3 + qy) +=
495 (mult * temp(px, 0) * temp(py, 1)*temp(pz,2) * temp(qx, 3) * temp(qy, 4)*temp(qz,5));
496 } // end for qy
497 } // end for qx
498 } // end for qz
499 } // end for py
500 } // end for px
501 } //end for pz
502 temp.format();
503 } // end cubature
504 // scale coefficient matrix
505 coef *= DataType(8) / mean;
506 temp.set_inverse(coef);
507 coef.set_transpose(temp);
508 } // end prepare
509
510 void finish()
511 {
512 _trafo_eval.finish();
513 BaseClass::finish();
514 }
515
516 int get_num_assigned_dofs() const
517 {
518 return max_assigned_dofs;
519 }
520
521 template<typename NodeData_, typename Function_>
522 void operator()(NodeData_& node_data, const Function_& function) const
523 {
524 static_assert(std::is_base_of<Analytic::Function, Function_>::value, "invalid function object");
525 static_assert(Function_::can_value, "function cannot compute values");
526
527 // declare our evaluation traits
528 typedef Analytic::EvalTraits<DataType_, Function_> FuncEvalTraits;
529
530 // declare function evaluator
531 typename Function_::template Evaluator<FuncEvalTraits> func_eval(function);
532
533 // compute trafo data
534 typename FuncEvalTraits::ValueType value(DataType_(0));
535 DataType_ mean(DataType_(0));
536 TrafoEvalData trafo_data;
537
538
539 // integrate over facet
540 for (int i(0); i < _cub_rule.get_num_points(); ++i)
541 {
542 _trafo_eval(trafo_data, _cub_rule.get_point(i));
543 value += (_cub_rule.get_weight(i) * trafo_data.jac_det) * func_eval.value(trafo_data.img_point)
544 * Intern::base_q(trafo_data.dom_point[0], trafo_data.dom_point[1], trafo_data.dom_point[2], coef[26]);
545 mean += _cub_rule.get_weight(i) * trafo_data.jac_det;
546 }
547
548 // set node_date to integral mean
549 node_data[0] = (DataType_(8) / mean) * value;
550 }
551 };
552 } // namespace Bernstein2
553 } // namespace Space
554} // namespace FEAT
base_header.hpp
FEAT Kernel base header.
FEAT::Cubature::Rule
Cubature Rule class template.
Definition: rule.hpp:38
FEAT::Space::Bernstein2::Evaluator
Bernstein-2 Element Evaluator class template declaration.
Definition: evaluator.hpp:36
FEAT::Space::NodeFunctionalBase
Node-functional base class template.
Definition: node_functional_base.hpp:30
FEAT::Space::NodeFunctionalBase::finish
void finish()
Releases the node-functional from the current cell.
Definition: node_functional_base.hpp:70
FEAT::Space::NodeFunctionalBase::prepare
void prepare(Index cell_index)
Prepares the node-functional for a given cell.
Definition: node_functional_base.hpp:62
FEAT::Space::NodeFunctionalNull
Null-Node-Functional class template.
Definition: node_functional_base.hpp:122
FEAT::Space::NodeFunctionalNull::get_num_assigned_dofs
int get_num_assigned_dofs() const
Returns the number of assigned dofs on the current cell.
Definition: node_functional_base.hpp:138
FEAT::Space::NodeFunctionalNull::BaseClass
NodeFunctionalBase< Space_, DataType_ > BaseClass
base-class typedef
Definition: node_functional_base.hpp:128
FEAT::Space::NodeFunctionalNull::operator()
void operator()(NodeData_ &node_data, const Function_ &function) const
Evaluation operator.
Definition: node_functional_base.hpp:147
FEAT::Tiny::Matrix
Tiny Matrix class template.
Definition: tiny_algebra.hpp:84
FEAT::Tiny::Matrix::set_transpose
CUDA_HOST_DEVICE Matrix & set_transpose(const Matrix< T_, n_, m_, sma_, sna_ > &a)
Sets this matrix to the transpose of another matrix.
Definition: tiny_algebra.hpp:1879
FEAT::Tiny::Matrix::format
CUDA_HOST_DEVICE void format(DataType alpha=DataType(0))
Formats the matrix.
Definition: tiny_algebra.hpp:1638
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::SpaceTags::value
@ value
specifies whether the space should supply basis function values
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122
FEAT::TrafoTags
TrafoTags
Trafo configuration tags enum.
Definition: eval_tags.hpp:22
FEAT::TrafoTags::img_point
@ img_point
specifies whether the trafo should supply image point coordinates
FEAT::TrafoTags::dom_point
@ dom_point
specifies whether the trafo should supply domain point coordinates
FEAT::TrafoTags::jac_det
@ jac_det
specifies whether the trafo should supply jacobian determinants
FEAT::Shape::Hypercube
Hypercube shape tag struct template.
Definition: shape.hpp:64