FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
refine_factory.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/cubature/rule.hpp>
10
11namespace FEAT
12{
13 namespace Cubature
14 {
16 namespace Intern
17 {
18 template<
19 typename Rule_,
20 typename Shape_ = typename Rule_::ShapeType>
21 class RuleRefinery;
22 } // namespace Intern
24
25 class RefineFactoryCore
26 {
27 public:
28 template<typename Shape_, typename Weight_, typename Coord_, typename Point_>
29 static void create(
30 Rule<Shape_, Weight_, Coord_, Point_>& rule,
31 const Rule<Shape_, Weight_, Coord_, Point_>& rule_in,
32 Index num_refines = 1)
33 {
34 typedef Rule<Shape_, Weight_, Coord_, Point_> RuleType;
35 typedef Intern::RuleRefinery<RuleType> RefineryType;
36
37 if(num_refines == 0)
38 {
39 rule = std::move(rule_in.clone());
40 }
41 else if(num_refines == 1)
42 {
43 rule = Rule<Shape_, Weight_, Coord_, Point_>(rule_in.get_num_points() * RefineryType::count,
44 "refine:" + rule_in.get_name());
45 RefineryType::refine(rule, rule_in);
46 }
47 else
48 {
49 // copy input rule
50 rule = std::move(rule_in.clone());
51 for(Index i(0); i < num_refines; ++i)
52 {
53 RuleType rule_tmp(rule.get_num_points() * RefineryType::count,
54 "refine*" + stringify(i+1) + ":" + rule_in.get_name());
55 RefineryType::refine(rule_tmp, rule);
56 rule = std::move(rule_tmp);
57 }
58 }
59 }
60 };
61
62 template<typename Factory_>
63 class RefineFactoryBase :
64 public RefineFactoryCore
65 {
66 public:
67 typedef Factory_ FactoryType;
68 typedef typename Factory_::ShapeType ShapeType;
69 typedef RefineFactoryCore BaseClass;
70
71/*
72 template<typename Weight_, typename Coord_, typename Point_>
73 static void create(
74 Rule<ShapeType, Weight_, Coord_, Point_>& rule,
75 const Rule<ShapeType, Weight_, Coord_, Point_>& rule_in,
76 Index num_refines = 1)
77 {
78 typedef Rule<ShapeType, Weight_, Coord_, Point_> RuleType;
79 typedef Intern::RuleRefinery<RuleType> RefineryType;
80
81 if(num_refines == 0)
82 {
83 rule = rule_in;
84 }
85 if(num_refines == 1)
86 {
87 RuleType rule(rule_in.get_num_points() * RefineryType::count, "refine:" + rule_in.get_name());
88 RefineryType::refine(rule, rule_in);
89 }
90
91 // copy input rule
92 rule = rule_in;
93 for(Index i(0); i < num_refines; ++i)
94 {
95 RuleType rule_tmp(rule.get_num_points() * RefineryType::count,
96 "refine*" + stringify(i+1) + ":" + rule_in.get_name());
97 RefineryType::refine(rule_tmp, rule);
98 rule = rule_tmp;
99 }
100 }*/
101
102 using BaseClass::create;
103
104 template<typename Weight_, typename Coord_, typename Point_>
105 static bool create(Rule<ShapeType, Weight_, Coord_, Point_>& rule, const String& name)
106 {
107 // try to find a colon within the name string
108 String::size_type k = name.find_first_of(':');
109 if(k == name.npos)
110 return false;
111
112 // extract substrings until the colon
113 String head(name.substr(0, k));
114 String tail(name.substr(k + 1));
115
116 // let's assume we refine 1 time
117 Index num_refines = 1;
118
119 // try to find a star within the head string
120 k = head.find_first_of('*');
121 if(k != head.npos)
122 {
123 // try to parse refine count
124 if(!String(head.substr(k+1)).parse(num_refines))
125 return false;
126
127 // trim head of refine count
128 head = head.substr(0, k);
129 }
130
131 // check head - this is the name of the formula
132 if(head.trim().compare_no_case("refine") != 0)
133 return false;
134
135 // call factory to create the input rule
136 Rule<ShapeType, Weight_, Coord_, Point_> rule_in;
137 if(!FactoryType::create(rule_in, tail.trim()))
138 return false;
139
140 // convert input rule
141 create(rule, rule_in, num_refines);
142 return true;
143 }
144 };
145
146 template<
147 typename Factory_,
148 bool variadic_ = Factory_::variadic>
149 class RefineFactory DOXY({});
150
151 template<typename Factory_>
152 class RefineFactory<Factory_, false> :
153 public RefineFactoryBase<Factory_>
154 {
155 public:
156 typedef Factory_ FactoryType;
157 typedef typename Factory_::ShapeType ShapeType;
158 typedef RefineFactoryBase<FactoryType> BaseClass;
159 static constexpr bool variadic = false;
160 static constexpr int num_points = FactoryType::num_points;
161
162 protected:
163 Index _num_refines;
164
165 public:
166 explicit RefineFactory(Index num_refines = 1) :
167 _num_refines(num_refines)
168 {
169 }
170
171 using BaseClass::create;
172
173 template<typename Weight_, typename Coord_, typename Point_>
174 void create(Rule<ShapeType, Weight_, Coord_, Point_>& rule)
175 {
176 create(rule, _num_refines);
177 }
178
179 template<typename Weight_, typename Coord_, typename Point_>
180 static void create(Rule<ShapeType, Weight_, Coord_, Point_>& rule, Index num_refines)
181 {
182 Rule<ShapeType, Weight_, Coord_, Point_> rule_in;
183 FactoryType::create(rule_in);
184 create(rule, rule_in, num_refines);
185 }
186 };
187
188 template<typename Factory_>
189 class RefineFactory<Factory_, true> :
190 public RefineFactoryBase<Factory_>
191 {
192 public:
193 typedef Factory_ FactoryType;
194 typedef typename Factory_::ShapeType ShapeType;
195 typedef RefineFactoryBase<FactoryType> BaseClass;
196 static constexpr bool variadic = true;
197 static constexpr int min_points = FactoryType::min_points;
198 static constexpr int max_points = FactoryType::max_points;
199
200 protected:
201 int _num_points;
202 Index _num_refines;
203
204 public:
205 explicit RefineFactory(int num_points, Index num_refines = 1) :
206 _num_points(num_points),
207 _num_refines(num_refines)
208 {
209 }
210
211 using BaseClass::create;
212
213 template<typename Weight_, typename Coord_, typename Point_>
214 void create(Rule<ShapeType, Weight_, Coord_, Point_>& rule)
215 {
216 create(rule, _num_points, _num_refines);
217 }
218
219 template<typename Weight_, typename Coord_, typename Point_>
220 static void create(Rule<ShapeType, Weight_, Coord_, Point_>& rule, int num_points, Index num_refines)
221 {
222 Rule<ShapeType, Weight_, Coord_, Point_> rule_in;
223 FactoryType::create(rule_in, num_points);
224 create(rule, rule_in, num_refines);
225 }
226 };
227
229 namespace Intern
230 {
231 template<typename Rule_>
232 class RuleRefinery<Rule_, Shape::Simplex<1> >
233 {
234 public:
235 typedef typename Rule_::WeightType Weight_;
236 typedef typename Rule_::CoordType Coord_;
237 static constexpr int count = 2;
238
239 static void refine(Rule_& rule, const Rule_& rule_in)
240 {
241 int n = rule_in.get_num_points();
242 for(int i(0); i < n; ++i)
243 {
244 rule.get_coord( i, 0) = Coord_(0.5) * rule_in.get_coord(i,0);
245 rule.get_coord(n+i, 0) = Coord_(0.5) * (rule_in.get_coord(i,0) + Coord_(1));
246 rule.get_weight( i) = Weight_(0.5) * rule_in.get_weight(i);
247 rule.get_weight(n+i) = Weight_(0.5) * rule_in.get_weight(i);
248 }
249 }
250 }; //RuleRefinery<Rule_, Shape::Simplex<1> >
251
252 template<typename Rule_>
253 class RuleRefinery<Rule_, Shape::Simplex<2> >
254 {
255 public:
256 typedef typename Rule_::WeightType Weight_;
257 typedef typename Rule_::CoordType Coord_;
258 static constexpr int count = 4;
259
260 static void refine(Rule_& rule, const Rule_& rule_in)
261 {
262 int n = rule_in.get_num_points();
263
264 // points of the child triangles
265 Coord_ v[4][3][2];
266 v[0][0][0] = Coord_(0);
267 v[0][0][1] = Coord_(0);
268 v[0][1][0] = Coord_(1)/Coord_(2);
269 v[0][1][1] = Coord_(0);
270 v[0][2][0] = Coord_(0);
271 v[0][2][1] = Coord_(1)/Coord_(2);
272
273 v[1][0][0] = Coord_(1)/Coord_(2);
274 v[1][0][1] = Coord_(0);
275 v[1][1][0] = Coord_(1);
276 v[1][1][1] = Coord_(0);
277 v[1][2][0] = Coord_(1)/Coord_(2);
278 v[1][2][1] = Coord_(1)/Coord_(2);
279
280 v[2][0][0] = Coord_(0);
281 v[2][0][1] = Coord_(1)/Coord_(2);
282 v[2][1][0] = Coord_(1)/Coord_(2);
283 v[2][1][1] = Coord_(1)/Coord_(2);
284 v[2][2][0] = Coord_(0);
285 v[2][2][1] = Coord_(1);
286
287 v[3][0][0] = Coord_(1)/Coord_(2);
288 v[3][0][1] = Coord_(1)/Coord_(2);
289 v[3][1][0] = Coord_(0);
290 v[3][1][1] = Coord_(1)/Coord_(2);
291 v[3][2][0] = Coord_(1)/Coord_(2);
292 v[3][2][1] = Coord_(0);
293
294 // refining
295 for(int i(0); i <= 3; ++i)
296 {
297 for(int j(0); j < n; ++j)
298 {
299 Coord_ x = rule_in.get_coord(j,0);
300 Coord_ y = rule_in.get_coord(j,1);
301 Coord_ w = Coord_(1) - x - y;
302 rule.get_coord(i*n + j, 0) = w*v[i][0][0] + x*v[i][1][0] + y*v[i][2][0];
303 rule.get_coord(i*n + j, 1) = w*v[i][0][1] + x*v[i][1][1] + y*v[i][2][1];
304 rule.get_weight(i*n + j) = Weight_(0.25) * rule_in.get_weight(j);
305 }
306 }
307 }
308 }; //RuleRefinery<Rule_, Shape::Simplex<2> >
309
310 template<typename Rule_>
311 class RuleRefinery<Rule_, Shape::Simplex<3> >
312 {
313 public:
314 typedef typename Rule_::WeightType Weight_;
315 typedef typename Rule_::CoordType Coord_;
316 static constexpr int count = 12;
317
318 static void refine(Rule_& rule, const Rule_& rule_in)
319 {
320 // points of the child triangles
321 Coord_ v[12][4][3];
322
323 // v_0-child nr.0
324
325 //v_4
326 v[0][0][0] = Coord_(1)/Coord_(2);
327 v[0][0][1] = Coord_(0);
328 v[0][0][2] = Coord_(0);
329
330 //v_6
331 v[0][1][0] = Coord_(0);
332 v[0][1][1] = Coord_(0);
333 v[0][1][2] = Coord_(1)/Coord_(2);
334
335 //v_5
336 v[0][2][0] = Coord_(0);
337 v[0][2][1] = Coord_(1)/Coord_(2);
338 v[0][2][2] = Coord_(0);
339
340 //v_0
341 v[0][3][0] = Coord_(0);
342 v[0][3][1] = Coord_(0);
343 v[0][3][2] = Coord_(0);
344
345
346 // v_0-child nr.1
347
348 //v_4
349 v[1][0][0] = Coord_(1)/Coord_(2);
350 v[1][0][1] = Coord_(0);
351 v[1][0][2] = Coord_(0);
352
353 //v_5
354 v[1][1][0] = Coord_(0);
355 v[1][1][1] = Coord_(1)/Coord_(2);
356 v[1][1][2] = Coord_(0);
357
358 //v_6
359 v[1][2][0] = Coord_(0);
360 v[1][2][1] = Coord_(0);
361 v[1][2][2] = Coord_(1)/Coord_(2);
362
363 //v_10
364 v[1][3][0] = Coord_(1)/Coord_(4);
365 v[1][3][1] = Coord_(1)/Coord_(4);
366 v[1][3][2] = Coord_(1)/Coord_(4);
367
368
369 // v_1-child nr.0
370
371 //v_4
372 v[2][0][0] = Coord_(1)/Coord_(2);
373 v[2][0][1] = Coord_(0);
374 v[2][0][2] = Coord_(0);
375
376 //v_7
377 v[2][1][0] = Coord_(1)/Coord_(2);
378 v[2][1][1] = Coord_(1)/Coord_(2);
379 v[2][1][2] = Coord_(0);
380
381 //v_8
382 v[2][2][0] = Coord_(1)/Coord_(2);
383 v[2][2][1] = Coord_(0);
384 v[2][2][2] = Coord_(1)/Coord_(2);
385
386 //v_1
387 v[2][3][0] = Coord_(1);
388 v[2][3][1] = Coord_(0);
389 v[2][3][2] = Coord_(0);
390
391
392 // v_1-child nr.1
393
394 //v_4
395 v[3][0][0] = Coord_(1)/Coord_(2);
396 v[3][0][1] = Coord_(0);
397 v[3][0][2] = Coord_(0);
398
399 //v_8
400 v[3][1][0] = Coord_(1)/Coord_(2);
401 v[3][1][1] = Coord_(0);
402 v[3][1][2] = Coord_(1)/Coord_(2);
403
404 //v_7
405 v[3][2][0] = Coord_(1)/Coord_(2);
406 v[3][2][1] = Coord_(1)/Coord_(2);
407 v[3][2][2] = Coord_(0);
408
409 //v_10
410 v[3][3][0] = Coord_(1)/Coord_(4);
411 v[3][3][1] = Coord_(1)/Coord_(4);
412 v[3][3][2] = Coord_(1)/Coord_(4);
413
414
415 // v_2-child nr.0
416
417 //v_5
418 v[4][0][0] = Coord_(0);
419 v[4][0][1] = Coord_(1)/Coord_(2);
420 v[4][0][2] = Coord_(0);
421
422 //v_9
423 v[4][1][0] = Coord_(0);
424 v[4][1][1] = Coord_(1)/Coord_(2);
425 v[4][1][2] = Coord_(1)/Coord_(2);
426
427 //v_7
428 v[4][2][0] = Coord_(1)/Coord_(2);
429 v[4][2][1] = Coord_(1)/Coord_(2);
430 v[4][2][2] = Coord_(0);
431
432 //v_2
433 v[4][3][0] = Coord_(0);
434 v[4][3][1] = Coord_(1);
435 v[4][3][2] = Coord_(0);
436
437
438 // v_2-child nr.1
439
440 //v_5
441 v[5][0][0] = Coord_(0);
442 v[5][0][1] = Coord_(1)/Coord_(2);
443 v[5][0][2] = Coord_(0);
444
445 //v_7
446 v[5][1][0] = Coord_(1)/Coord_(2);
447 v[5][1][1] = Coord_(1)/Coord_(2);
448 v[5][1][2] = Coord_(0);
449
450 //v_9
451 v[5][2][0] = Coord_(0);
452 v[5][2][1] = Coord_(1)/Coord_(2);
453 v[5][2][2] = Coord_(1)/Coord_(2);
454
455 //v_10
456 v[5][3][0] = Coord_(1)/Coord_(4);
457 v[5][3][1] = Coord_(1)/Coord_(4);
458 v[5][3][2] = Coord_(1)/Coord_(4);
459
460
461 // v_3-child nr.0
462
463 //v_6
464 v[6][0][0] = Coord_(0);
465 v[6][0][1] = Coord_(0);
466 v[6][0][2] = Coord_(1)/Coord_(2);
467
468 //v_8
469 v[6][1][0] = Coord_(1)/Coord_(2);
470 v[6][1][1] = Coord_(0);
471 v[6][1][2] = Coord_(1)/Coord_(2);
472
473 //v_9
474 v[6][2][0] = Coord_(0);
475 v[6][2][1] = Coord_(1)/Coord_(2);
476 v[6][2][2] = Coord_(1)/Coord_(2);
477
478 //v_3
479 v[6][3][0] = Coord_(0);
480 v[6][3][1] = Coord_(0);
481 v[6][3][2] = Coord_(1);
482
483
484 // v_3-child nr.1
485
486 //v_6
487 v[7][0][0] = Coord_(0);
488 v[7][0][1] = Coord_(0);
489 v[7][0][2] = Coord_(1)/Coord_(2);
490
491 //v_9
492 v[7][1][0] = Coord_(0);
493 v[7][1][1] = Coord_(1)/Coord_(2);
494 v[7][1][2] = Coord_(1)/Coord_(2);
495
496 //v_8
497 v[7][2][0] = Coord_(1)/Coord_(2);
498 v[7][2][1] = Coord_(0);
499 v[7][2][2] = Coord_(1)/Coord_(2);
500
501 //v_10
502 v[7][3][0] = Coord_(1)/Coord_(4);
503 v[7][3][1] = Coord_(1)/Coord_(4);
504 v[7][3][2] = Coord_(1)/Coord_(4);
505
506
507 //t_0-child
508
509 //v_7
510 v[8][0][0] = Coord_(1)/Coord_(2);
511 v[8][0][1] = Coord_(1)/Coord_(2);
512 v[8][0][2] = Coord_(0);
513
514 //v_8
515 v[8][1][0] = Coord_(1)/Coord_(2);
516 v[8][1][1] = Coord_(0);
517 v[8][1][2] = Coord_(1)/Coord_(2);
518
519 //v_9
520 v[8][2][0] = Coord_(0);
521 v[8][2][1] = Coord_(1)/Coord_(2);
522 v[8][2][2] = Coord_(1)/Coord_(2);
523
524 //v_10
525 v[8][3][0] = Coord_(1)/Coord_(4);
526 v[8][3][1] = Coord_(1)/Coord_(4);
527 v[8][3][2] = Coord_(1)/Coord_(4);
528
529
530 //t_1-child
531
532 //v_5
533 v[9][0][0] = Coord_(0);
534 v[9][0][1] = Coord_(1)/Coord_(2);
535 v[9][0][2] = Coord_(0);
536
537 //v_9
538 v[9][1][0] = Coord_(0);
539 v[9][1][1] = Coord_(1)/Coord_(2);
540 v[9][1][2] = Coord_(1)/Coord_(2);
541
542 //v_6
543 v[9][2][0] = Coord_(0);
544 v[9][2][1] = Coord_(0);
545 v[9][2][2] = Coord_(1)/Coord_(2);
546
547 //v_10
548 v[9][3][0] = Coord_(1)/Coord_(4);
549 v[9][3][1] = Coord_(1)/Coord_(4);
550 v[9][3][2] = Coord_(1)/Coord_(4);
551
552
553 // t_2-child
554
555 //v_4
556 v[10][0][0] = Coord_(1)/Coord_(2);
557 v[10][0][1] = Coord_(0);
558 v[10][0][2] = Coord_(0);
559
560 //v_6
561 v[10][1][0] = Coord_(0);
562 v[10][1][1] = Coord_(0);
563 v[10][1][2] = Coord_(1)/Coord_(2);
564
565 //v_8
566 v[10][2][0] = Coord_(1)/Coord_(2);
567 v[10][2][1] = Coord_(0);
568 v[10][2][2] = Coord_(1)/Coord_(2);
569
570 //v_10
571 v[10][3][0] = Coord_(1)/Coord_(4);
572 v[10][3][1] = Coord_(1)/Coord_(4);
573 v[10][3][2] = Coord_(1)/Coord_(4);
574
575
576 //t_3-child
577
578 //v_4
579 v[11][0][0] = Coord_(1)/Coord_(2);
580 v[11][0][1] = Coord_(0);
581 v[11][0][2] = Coord_(0);
582
583 //v_7
584 v[11][1][0] = Coord_(1)/Coord_(2);
585 v[11][1][1] = Coord_(1)/Coord_(2);
586 v[11][1][2] = Coord_(0);
587
588 //v_5
589 v[11][2][0] = Coord_(0);
590 v[11][2][1] = Coord_(1)/Coord_(2);
591 v[11][2][2] = Coord_(0);
592
593 //v_10
594 v[11][3][0] = Coord_(1)/Coord_(4);
595 v[11][3][1] = Coord_(1)/Coord_(4);
596 v[11][3][2] = Coord_(1)/Coord_(4);
597
598
599 int n = rule_in.get_num_points();
600 Weight_ vol;
601
602 // refining
603 for(int i(0); i <= 11; ++i)
604 {
605 for(int j(0); j < n; ++j)
606 {
607 Coord_ x = rule_in.get_coord(j,0);
608 Coord_ y = rule_in.get_coord(j,1);
609 Coord_ z = rule_in.get_coord(j,2);
610 Coord_ w = Coord_(1) - x - y - z;
611 rule.get_coord(i*n + j, 0) = w*v[i][0][0] + x*v[i][1][0] + y*v[i][2][0] + z*v[i][3][0];
612 rule.get_coord(i*n + j, 1) = w*v[i][0][1] + x*v[i][1][1] + y*v[i][2][1] + z*v[i][3][1];
613 rule.get_coord(i*n + j, 2) = w*v[i][0][2] + x*v[i][1][2] + y*v[i][2][2] + z*v[i][3][2];
614
615 // adjust the weights
616 /*
617 vol = (v[i][0][0] - v[i][3][0])*(v[i][1][1] - v[i][3][1])*(v[i][2][2] - v[i][3][2])
618 + (v[i][1][0] - v[i][3][0])*(v[i][2][1] - v[i][3][1])*(v[i][0][2] - v[i][3][2])
619 + (v[i][2][0] - v[i][3][0])*(v[i][0][1] - v[i][3][1])*(v[i][1][2] - v[i][3][2])
620 - (v[i][2][0] - v[i][3][0])*(v[i][1][1] - v[i][3][1])*(v[i][0][2] - v[i][3][2])
621 - (v[i][1][0] - v[i][3][0])*(v[i][0][1] - v[i][3][1])*(v[i][2][2] - v[i][3][2])
622 - (v[i][0][0] - v[i][3][0])*(v[i][2][1] - v[i][3][1])*(v[i][1][2] - v[i][3][2]);
623 */
624 if(i == 0 || i == 2 || i == 4|| i == 6)
625 {
626 vol = Weight_(1)/Weight_(8);
627 }
628 else
629 {
630 vol = Weight_(1)/Weight_(16);
631 }
632
633 rule.get_weight(i*n + j) = rule_in.get_weight(j)*Math::abs(vol);
634 }
635 }
636 }
637 }; //RuleRefinery<Rule_, Shape::Simplex<3> >
638
639 template<typename Rule_>
640 class RuleRefinery<Rule_, Shape::Hypercube<1> >
641 {
642 public:
643 typedef typename Rule_::WeightType Weight_;
644 typedef typename Rule_::CoordType Coord_;
645 static constexpr int count = 2;
646
647 static void refine(Rule_& rule, const Rule_& rule_in)
648 {
649 int n = rule_in.get_num_points();
650 for(int i(0); i < n; ++i)
651 {
652 rule.get_coord( i, 0) = Coord_(0.5) * rule_in.get_coord(i,0) - Coord_(0.5);
653 rule.get_coord(n+i, 0) = Coord_(0.5) * rule_in.get_coord(i,0) + Coord_(0.5);
654 rule.get_weight( i) = Weight_(0.5) * rule_in.get_weight(i);
655 rule.get_weight(n+i) = Weight_(0.5) * rule_in.get_weight(i);
656 }
657 }
658
659 };//RuleRefinery<Rule_, Shape::Hypercube<1> >
660
661 template<typename Rule_>
662 class RuleRefinery<Rule_, Shape::Hypercube<2> >
663 {
664 public:
665 typedef typename Rule_::WeightType Weight_;
666 typedef typename Rule_::CoordType Coord_;
667 static constexpr int count = 4;
668
669 static void refine(Rule_& rule, const Rule_& rule_in)
670 {
671
672 int n = rule_in.get_num_points();
673
674 // refining
675 for(int i(0); i <= 3; ++i)
676 {
677 // shift
678 Coord_ dx = Coord_(0);
679 Coord_ dy = Coord_(0);
680
681 if(i == 0)
682 {
683 dx = -Coord_(1)/Coord_(2);
684 dy = -Coord_(1)/Coord_(2);
685 }
686 else if(i == 1)
687 {
688 dx = Coord_(1)/Coord_(2);
689 dy = -Coord_(1)/Coord_(2);
690 }
691 else if(i == 2)
692 {
693 dx = -Coord_(1)/Coord_(2);
694 dy = Coord_(1)/Coord_(2);
695 }
696 else if(i == 3)
697 {
698 dx = Coord_(1)/Coord_(2);
699 dy = Coord_(1)/Coord_(2);
700 }
701
702 for(int j(0); j < n; ++j)
703 {
704 rule.get_coord(i*n + j, 0) = Coord_(0.5) * rule_in.get_coord(j,0) + dx;
705 rule.get_coord(i*n + j, 1) = Coord_(0.5) * rule_in.get_coord(j,1) + dy;
706 rule.get_weight(i*n + j) = rule_in.get_weight(j)/Weight_(4);
707 }
708 }
709 }
710 }; //RuleRefinery<Rule_, Shape::Hypercube<2> >
711
712 template<typename Rule_>
713 class RuleRefinery<Rule_, Shape::Hypercube<3> >
714 {
715 public:
716 typedef typename Rule_::WeightType Weight_;
717 typedef typename Rule_::CoordType Coord_;
718 static constexpr int count = 8;
719
720 static void refine(Rule_& rule, const Rule_& rule_in)
721 {
722
723 int n = rule_in.get_num_points();
724
725 // refining
726 for(int i(0); i <= 7; ++i)
727 {
728 // shift
729 Coord_ dx = Coord_(0);
730 Coord_ dy = Coord_(0);
731 Coord_ dz = Coord_(0);
732
733 if(i == 0)
734 {
735 dx = -Coord_(1)/Coord_(2);
736 dy = -Coord_(1)/Coord_(2);
737 dz = -Coord_(1)/Coord_(2);
738 }
739 else if(i == 1)
740 {
741 dx = Coord_(1)/Coord_(2);
742 dy = -Coord_(1)/Coord_(2);
743 dz = -Coord_(1)/Coord_(2);
744 }
745 else if(i == 2)
746 {
747 dx = -Coord_(1)/Coord_(2);
748 dy = Coord_(1)/Coord_(2);
749 dz = -Coord_(1)/Coord_(2);
750 }
751 else if(i == 3)
752 {
753 dx = Coord_(1)/Coord_(2);
754 dy = Coord_(1)/Coord_(2);
755 dz = -Coord_(1)/Coord_(2);
756 }
757 else if(i == 4)
758 {
759 dx = -Coord_(1)/Coord_(2);
760 dy = -Coord_(1)/Coord_(2);
761 dz = Coord_(1)/Coord_(2);
762 }
763 else if(i == 5)
764 {
765 dx = Coord_(1)/Coord_(2);
766 dy = -Coord_(1)/Coord_(2);
767 dz = Coord_(1)/Coord_(2);
768 }
769 else if(i == 6)
770 {
771 dx = -Coord_(1)/Coord_(2);
772 dy = Coord_(1)/Coord_(2);
773 dz = Coord_(1)/Coord_(2);
774 }
775 else if(i == 7)
776 {
777 dx = Coord_(1)/Coord_(2);
778 dy = Coord_(1)/Coord_(2);
779 dz = Coord_(1)/Coord_(2);
780 }
781
782 for(int j(0); j < n; ++j)
783 {
784 rule.get_coord(i*n + j, 0) = Coord_(0.5) * rule_in.get_coord(j,0) + dx;
785 rule.get_coord(i*n + j, 1) = Coord_(0.5) * rule_in.get_coord(j,1) + dy;
786 rule.get_coord(i*n + j, 2) = Coord_(0.5) * rule_in.get_coord(j,2) + dz;
787 rule.get_weight(i*n + j) = rule_in.get_weight(j)/Weight_(8);
788 }
789 }
790 }
791 }; //RuleRefinery<Rule_, Shape::Hypercube<3> >
792
793 } // namespace Intern
795 } // namespace Cubature
796} // namespace FEAT
FEAT::Cubature::Rule
Cubature Rule class template.
Definition: rule.hpp:38
FEAT::String
String class implementation.
Definition: string.hpp:47
FEAT::String::parse
bool parse(T_ &t) const
Parses the string and stores its value in the provided variable.
Definition: string.hpp:838
FEAT::String::compare_no_case
int compare_no_case(const String &other) const
Compares two strings without regard to case.
Definition: string.hpp:680
FEAT::String::trim
String trim(const String &charset) const
Trims the string.
Definition: string.hpp:328
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::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122