FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
bilinear_operator_assembler.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/assembly/asm_traits.hpp>
10#include <kernel/util/likwid_marker.hpp>
11
12namespace FEAT
13{
14 namespace Assembly
15 {
23 class BilinearOperatorAssembler
24 {
25 public:
49 template<
50 typename Matrix_,
51 typename Operator_,
52 typename TestSpace_,
53 typename TrialSpace_,
54 typename CubatureFactory_>
55 static void assemble_matrix2(
56 Matrix_& matrix,
57 Operator_& operat,
58 const TestSpace_& test_space,
59 const TrialSpace_& trial_space,
60 const CubatureFactory_& cubature_factory,
61 typename Matrix_::DataType alpha = typename Matrix_::DataType(1))
62 {
63 // validate matrix dimensions
64 XASSERTM(matrix.rows() == test_space.get_num_dofs(), "invalid matrix dimensions");
65 XASSERTM(matrix.columns() == trial_space.get_num_dofs(), "invalid matrix dimensions");
66
67 // matrix type
68 typedef Matrix_ MatrixType;
69 // operator type
70 typedef Operator_ OperatorType;
71 // test-space type
72 typedef TestSpace_ TestSpaceType;
73 // trial-space type
74 typedef TrialSpace_ TrialSpaceType;
75
76 // assembly traits
77 typedef AsmTraits2<
78 typename MatrixType::DataType,
79 TestSpaceType,
80 TrialSpaceType,
81 OperatorType::trafo_config,
82 OperatorType::test_config,
83 OperatorType::trial_config> AsmTraits;
84
85 // fetch the trafo
86 const typename AsmTraits::TrafoType& trafo = test_space.get_trafo();
87
88 // create a trafo evaluator
89 typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
90
91 // create space evaluators
92 typename AsmTraits::TestEvaluator test_eval(test_space);
93 typename AsmTraits::TrialEvaluator trial_eval(trial_space);
94
95 // create dof-mappings
96 typename AsmTraits::TestDofMapping test_dof_mapping(test_space);
97 typename AsmTraits::TrialDofMapping trial_dof_mapping(trial_space);
98
99 // create a operator evaluator
100 typename OperatorType::template Evaluator<AsmTraits> oper_eval(operat);
101
102 // create trafo evaluation data
103 typename AsmTraits::TrafoEvalData trafo_data;
104
105 // create space evaluation data
106 typename AsmTraits::TestEvalData test_data;
107 typename AsmTraits::TrialEvalData trial_data;
108
109 // the value type of the operator
110 typedef typename OperatorType::template Evaluator<AsmTraits>::ValueType ValueType;
111
112 // ensure that the operator and matrix value types are compatible
113 static_assert(std::is_same<ValueType, typename MatrixType::ValueType>::value,
114 "matrix and bilinear operator have different value types!");
115
116 // create local matrix data
117 typename AsmTraits::template TLocalMatrix<ValueType> loc_mat;
118
119 // create cubature rule
120 typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
121
122 // create matrix scatter-axpy
123 typename MatrixType::ScatterAxpy scatter_axpy(matrix);
124 FEAT_APPLICATION_MARKER_START("bilin_op_asm_mat2");
125 // loop over all cells of the mesh
126 for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
127 {
128 // prepare trafo evaluator
129 trafo_eval.prepare(cell);
130
131 // prepare space evaluators
132 test_eval.prepare(trafo_eval);
133 trial_eval.prepare(trafo_eval);
134
135 // prepare operator evaluator
136 oper_eval.prepare(trafo_eval);
137
138 // fetch number of local dofs
139 int num_loc_test_dofs = test_eval.get_num_local_dofs();
140 int num_loc_trial_dofs = trial_eval.get_num_local_dofs();
141
142 // format local matrix
143 loc_mat.format();
144
145 // loop over all quadrature points and integrate
146 for(int k(0); k < cubature_rule.get_num_points(); ++k)
147 {
148 // compute trafo data
149 trafo_eval(trafo_data, cubature_rule.get_point(k));
150
151 // compute basis function data
152 test_eval(test_data, trafo_data);
153 trial_eval(trial_data, trafo_data);
154
155 // prepare bilinear operator
156 oper_eval.set_point(trafo_data);
157
158 // test function loop
159 for(int i(0); i < num_loc_test_dofs; ++i)
160 {
161 // trial function loop
162 for(int j(0); j < num_loc_trial_dofs; ++j)
163 {
164 // evaluate operator and integrate
165 Tiny::axpy(loc_mat(i,j), oper_eval.eval(trial_data.phi[j], test_data.phi[i]),
166 trafo_data.jac_det * cubature_rule.get_weight(k));
167 // continue with next trial function
168 }
169 // continue with next test function
170 }
171 // continue with next cubature point
172 }
173 FEAT_APPLICATION_MARKER_STOP("bilin_op_asm_mat2");
174
175 // finish operator evaluator
176 oper_eval.finish();
177
178 // finish evaluators
179 trial_eval.finish();
180 test_eval.finish();
181 trafo_eval.finish();
182
183 // initialize dof-mappings
184 test_dof_mapping.prepare(cell);
185 trial_dof_mapping.prepare(cell);
186
187 // incorporate local matrix
188 scatter_axpy(loc_mat, test_dof_mapping, trial_dof_mapping, alpha);
189
190 // finish dof mapping
191 trial_dof_mapping.finish();
192 test_dof_mapping.finish();
193
194 // continue with next cell
195 }
196
197 // okay, that's it
198 }
199
220 template<
221 typename Matrix_,
222 typename Operator_,
223 typename Space_,
224 typename CubatureFactory_>
225 static void assemble_matrix1(
226 Matrix_& matrix,
227 Operator_& operat,
228 const Space_& space,
229 const CubatureFactory_& cubature_factory,
230 typename Matrix_::DataType alpha = typename Matrix_::DataType(1))
231 {
232 // validate matrix dimensions
233 XASSERTM(matrix.rows() == space.get_num_dofs(), "invalid matrix dimensions");
234 XASSERTM(matrix.columns() == space.get_num_dofs(), "invalid matrix dimensions");
235
236 // matrix type
237 typedef Matrix_ MatrixType;
238 // operator type
239 typedef Operator_ OperatorType;
240 // space type
241 typedef Space_ SpaceType;
242
243 // assembly traits
244 typedef AsmTraits1<
245 typename MatrixType::DataType,
246 SpaceType,
247 OperatorType::trafo_config,
248 OperatorType::test_config | OperatorType::trial_config> AsmTraits;
249
250 // fetch the trafo
251 const typename AsmTraits::TrafoType& trafo = space.get_trafo();
252
253 // create a trafo evaluator
254 typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
255
256 // create a space evaluator and evaluation data
257 typename AsmTraits::SpaceEvaluator space_eval(space);
258
259 // create a dof-mapping
260 typename AsmTraits::DofMapping dof_mapping(space);
261
262 // create a operator evaluator
263 typename OperatorType::template Evaluator<AsmTraits> oper_eval(operat);
264
265 // create trafo evaluation data
266 typename AsmTraits::TrafoEvalData trafo_data;
267
268 // create space evaluation data
269 typename AsmTraits::SpaceEvalData space_data;
270
271 // the value type of the operator
272 typedef typename OperatorType::template Evaluator<AsmTraits>::ValueType ValueType;
273
274 // ensure that the operator and matrix value types are compatible
275 static_assert(std::is_same<ValueType, typename MatrixType::ValueType>::value,
276 "matrix and bilinear operator have different value types!");
277
278 // create local matrix data
279 typename AsmTraits::template TLocalMatrix<ValueType> loc_mat;
280
281 // create cubature rule
282 typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
283
284 // create matrix scatter-axpy
285 typename MatrixType::ScatterAxpy scatter_axpy(matrix);
286 FEAT_APPLICATION_MARKER_START("bilin_op_asm:matrix1");
287
288 // loop over all cells of the mesh
289 for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
290 {
291 // prepare trafo evaluator
292 trafo_eval.prepare(cell);
293
294 // prepare space evaluator
295 space_eval.prepare(trafo_eval);
296
297 // prepare operator evaluator
298 oper_eval.prepare(trafo_eval);
299
300 // fetch number of local dofs
301 int num_loc_dofs = space_eval.get_num_local_dofs();
302
303 // format local matrix
304 loc_mat.format();
305
306 // loop over all quadrature points and integrate
307 for(int k(0); k < cubature_rule.get_num_points(); ++k)
308 {
309 // compute trafo data
310 trafo_eval(trafo_data, cubature_rule.get_point(k));
311
312 // compute basis function data
313 space_eval(space_data, trafo_data);
314
315 // prepare bilinear operator
316 oper_eval.set_point(trafo_data);
317
318 // test function loop
319 for(int i(0); i < num_loc_dofs; ++i)
320 {
321 // trial function loop
322 for(int j(0); j < num_loc_dofs; ++j)
323 {
324 // evaluate operator and integrate
325 Tiny::axpy(loc_mat(i,j), oper_eval.eval(space_data.phi[j], space_data.phi[i]),
326 trafo_data.jac_det * cubature_rule.get_weight(k));
327 // continue with next trial function
328 }
329 // continue with next test function
330 }
331 // continue with next cubature point
332 }
333 FEAT_APPLICATION_MARKER_STOP("bilin_op_asm:matrix1");
334 // finish operator evaluator
335 oper_eval.finish();
336
337 // finish evaluators
338 space_eval.finish();
339 trafo_eval.finish();
340
341 // initialize dof-mapping
342 dof_mapping.prepare(cell);
343
344 // incorporate local matrix
345 scatter_axpy(loc_mat, dof_mapping, dof_mapping, alpha);
346
347 // finish dof mapping
348 dof_mapping.finish();
349
350 // continue with next cell
351 }
352
353 // okay, that's it
354 }
355
383 template
384 <
385 typename Vector_,
386 typename Operator_,
387 typename Space_,
388 typename CubatureFactory_
389 >
390 static void apply1(
391 Vector_& ret,
392 const Vector_& coeff_vector,
393 const Operator_& operat,
394 const Space_& space,
395 const CubatureFactory_& cubature_factory,
396 typename Vector_::DataType alpha = typename Vector_::DataType(1))
397 {
398 XASSERTM(ret.size()==space.get_num_dofs(), "Return vector size does not match FE space dimension");
399 XASSERTM(coeff_vector.size()==space.get_num_dofs(), "Coefficient vector size does not match FE space dimension");
400 ret.format();
401
402 // Coefficient vector type
403 typedef Vector_ VectorType;
404 // operator type
405 typedef Operator_ OperatorType;
406 // space type
407 typedef Space_ SpaceType;
408
409 // assembly traits
410 typedef AsmTraits1
411 <
412 typename VectorType::DataType,
413 SpaceType,
414 OperatorType::trafo_config,
415 OperatorType::test_config | OperatorType::trial_config
416 > AsmTraits;
417
418 // fetch the trafo
419 const typename AsmTraits::TrafoType& trafo = space.get_trafo();
420
421 // create a trafo evaluator
422 typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
423
424 // create a space evaluator and evaluation data
425 typename AsmTraits::SpaceEvaluator space_eval(space);
426
427 // create a dof-mapping
428 typename AsmTraits::DofMapping dof_mapping(space);
429
430 // create a operator evaluator
431 typename OperatorType::template Evaluator<AsmTraits> oper_eval(operat);
432
433 // create trafo evaluation data
434 typename AsmTraits::TrafoEvalData trafo_data;
435
436 // create space evaluation data
437 typename AsmTraits::SpaceEvalData space_data;
438
439 // get the value type from the vector
440 typedef typename VectorType::ValueType ValueType;
441
442 // create local vector data
443 typename AsmTraits::template TLocalVector<ValueType> loc_coeff, loc_ret;
444
445 // create cubature rule
446 typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
447
448 // create vector scatter axpy for adding to the global return vector
449 typename VectorType::ScatterAxpy scatter_axpy(ret);
450 // create vector gather axpy for picking the local values from the global vector
451 typename VectorType::GatherAxpy gather_axpy(coeff_vector);
452 FEAT_APPLICATION_MARKER_START("bilin_op_asm:apply1");
453 // loop over all cells of the mesh
454 for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
455 {
456 // format local vectors
457 loc_coeff.format();
458 loc_ret.format();
459
460 // initialize dof-mapping
461 dof_mapping.prepare(cell);
462
463 // prepare trafo evaluator
464 trafo_eval.prepare(cell);
465
466 // incorporate local vector
467 gather_axpy(loc_coeff, dof_mapping);
468
469 // prepare space evaluator
470 space_eval.prepare(trafo_eval);
471
472 // prepare operator evaluator
473 oper_eval.prepare(trafo_eval);
474
475 // fetch number of local dofs
476 int num_loc_dofs = space_eval.get_num_local_dofs();
477
478 // loop over all quadrature points and integrate
479 for(int k(0); k < cubature_rule.get_num_points(); ++k)
480 {
481 // compute trafo data
482 trafo_eval(trafo_data, cubature_rule.get_point(k));
483
484 // compute basis function data
485 space_eval(space_data, trafo_data);
486
487 // prepare bilinear operator
488 oper_eval.set_point(trafo_data);
489
490 // test function loop
491 for(int i(0); i < num_loc_dofs; ++i)
492 {
493 // trial function loop
494 for(int j(0); j < num_loc_dofs; ++j)
495 {
496 // Integrate the FE function
497 Tiny::axpy(loc_ret(i),
498 oper_eval.eval(space_data.phi[j], space_data.phi[i]) * loc_coeff(j),
499 trafo_data.jac_det * cubature_rule.get_weight(k));
500 // continue with next test function
501 }
502 // continue with next trial function
503 }
504 // continue with next cubature point
505 }
506 FEAT_APPLICATION_MARKER_STOP("bilin_op_asm:apply1");
507 // finish operator evaluator
508 oper_eval.finish();
509
510 // finish evaluators
511 space_eval.finish();
512 trafo_eval.finish();
513
514 // incorporate local vector
515 scatter_axpy(loc_ret, dof_mapping, alpha);
516
517 // finish dof mapping
518 dof_mapping.finish();
519
520 // continue with next cell
521 }
522
523 // okay, that's it
524 }
525
556 template
557 <
558 typename Vector_,
559 typename Operator_,
560 typename TestSpace_,
561 typename TrialSpace_,
562 typename CubatureFactory_
563 >
564 static void apply2(
565 Vector_& ret,
566 const Vector_& coeff_vector,
567 const Operator_& operat,
568 const TestSpace_& test_space,
569 const TrialSpace_& trial_space,
570 const CubatureFactory_& cubature_factory,
571 typename Vector_::DataType alpha = typename Vector_::DataType(1))
572 {
573 XASSERTM(ret.size()==test_space.get_num_dofs(), "Return vector size does not match FE space dimension");
574 XASSERTM(coeff_vector.size()==trial_space.get_num_dofs(), "Coefficient vector size does not match FE space dimension");
575 ret.format();
576
577 // Coefficient vector type
578 typedef Vector_ VectorType;
579 // operator type
580 typedef Operator_ OperatorType;
581 // test-space type
582 typedef TestSpace_ TestSpaceType;
583 // trial-space type
584 typedef TrialSpace_ TrialSpaceType;
585
586 // assembly traits
587 typedef AsmTraits2
588 <
589 typename VectorType::DataType,
590 TestSpaceType,
591 TrialSpaceType,
592 OperatorType::trafo_config,
593 OperatorType::test_config,
594 OperatorType::trial_config
595 > AsmTraits;
596
597 // fetch the trafo
598 const typename AsmTraits::TrafoType& trafo = test_space.get_trafo();
599
600 // create a trafo evaluator
601 typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
602
603 // create space evaluators
604 typename AsmTraits::TestEvaluator test_eval(test_space);
605 typename AsmTraits::TrialEvaluator trial_eval(trial_space);
606
607 // create dof-mappings
608 typename AsmTraits::TestDofMapping test_dof_mapping(test_space);
609 typename AsmTraits::TrialDofMapping trial_dof_mapping(trial_space);
610
611 // create a operator evaluator
612 typename OperatorType::template Evaluator<AsmTraits> oper_eval(operat);
613
614 // create trafo evaluation data
615 typename AsmTraits::TrafoEvalData trafo_data;
616
617 // create space evaluation data
618 typename AsmTraits::TestEvalData test_data;
619 typename AsmTraits::TrialEvalData trial_data;
620
621 // get the value type from the vector
622 typedef typename VectorType::ValueType ValueType;
623
624 // create local vector data
625 typename AsmTraits::template TLocalTrialVector<ValueType> loc_coeff;
626 typename AsmTraits::template TLocalTestVector<ValueType> loc_ret;
627
628 // create cubature rule
629 typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
630
631 // create vector scatter axpy for adding to the global return vector
632 typename VectorType::ScatterAxpy scatter_axpy(ret);
633 // create vector gather axpy for picking the local values from the global vector
634 typename VectorType::GatherAxpy gather_axpy(coeff_vector);
635
636 FEAT_APPLICATION_MARKER_START("bilin_op_asm:apply2");
637 // loop over all cells of the mesh
638 for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
639 {
640 // format local vectors
641 loc_coeff.format();
642 loc_ret.format();
643
644 // initialize dof-mappings
645 test_dof_mapping.prepare(cell);
646 trial_dof_mapping.prepare(cell);
647
648 // prepare trafo evaluator
649 trafo_eval.prepare(cell);
650 // prepare space evaluator
651 test_eval.prepare(trafo_eval);
652 trial_eval.prepare(trafo_eval);
653
654 // incorporate local vector
655 gather_axpy(loc_coeff, trial_dof_mapping);
656
657 // prepare operator evaluator
658 oper_eval.prepare(trafo_eval);
659
660 // fetch number of local dofs
661 int num_loc_test_dofs = test_eval.get_num_local_dofs();
662 int num_loc_trial_dofs = trial_eval.get_num_local_dofs();
663
664 // loop over all quadrature points and integrate
665 for(int k(0); k < cubature_rule.get_num_points(); ++k)
666 {
667 // compute trafo data
668 trafo_eval(trafo_data, cubature_rule.get_point(k));
669
670 // compute basis function data
671 test_eval(test_data, trafo_data);
672 trial_eval(trial_data, trafo_data);
673
674 // prepare bilinear operator
675 oper_eval.set_point(trafo_data);
676
677 // test function loop
678 for(int i(0); i < num_loc_test_dofs; ++i)
679 {
680 // trial function loop
681 for(int j(0); j < num_loc_trial_dofs; ++j)
682 {
683 // Integrate the FE function
684 Tiny::axpy(loc_ret(i),
685 oper_eval.eval(trial_data.phi[j], test_data.phi[i]) * loc_coeff(j),
686 trafo_data.jac_det * cubature_rule.get_weight(k));
687 // continue with next trial function
688 }
689 // continue with next test function
690 }
691 // continue with next cubature point
692 }
693 FEAT_KERNEL_MARKER_STOP("bilin_op_asm:apply2");
694
695 // finish operator evaluator
696 oper_eval.finish();
697
698 // finish evaluators
699 test_eval.finish();
700 trial_eval.finish();
701 trafo_eval.finish();
702
703 // incorporate local vector
704 scatter_axpy(loc_ret, test_dof_mapping, alpha);
705
706 // finish dof mappings
707 test_dof_mapping.finish();
708 trial_dof_mapping.finish();
709
710 // continue with next cell
711 }
712 // okay, that's it
713 }
714 }; // class BilinearOperatorAssembler<...>
715 } // namespace Assembly
716} // namespace FEAT
XASSERTM
#define XASSERTM(expr, msg)
Assertion macro definition with custom message.
Definition: assertion.hpp:263
FEAT::Assembly::AsmTraits1
Common single-space assembly traits class template.
Definition: asm_traits.hpp:83
FEAT::Assembly::AsmTraits2
Common test-/trial-space assembly traits class template.
Definition: asm_traits.hpp:227
FEAT::Assembly::BilinearOperatorAssembler
Bilinear Operator Assembler class template.
Definition: bilinear_operator_assembler.hpp:24
FEAT::Assembly::BilinearOperatorAssembler::assemble_matrix2
static void assemble_matrix2(Matrix_ &matrix, Operator_ &operat, const TestSpace_ &test_space, const TrialSpace_ &trial_space, const CubatureFactory_ &cubature_factory, typename Matrix_::DataType alpha=typename Matrix_::DataType(1))
Assembles a bilinear operator into a matrix.
Definition: bilinear_operator_assembler.hpp:55
FEAT::Assembly::BilinearOperatorAssembler::apply2
static void apply2(Vector_ &ret, const Vector_ &coeff_vector, const Operator_ &operat, const TestSpace_ &test_space, const TrialSpace_ &trial_space, const CubatureFactory_ &cubature_factory, typename Vector_::DataType alpha=typename Vector_::DataType(1))
Applies the bilinear operator to an FE function.
Definition: bilinear_operator_assembler.hpp:564
FEAT::Assembly::BilinearOperatorAssembler::apply1
static void apply1(Vector_ &ret, const Vector_ &coeff_vector, const Operator_ &operat, const Space_ &space, const CubatureFactory_ &cubature_factory, typename Vector_::DataType alpha=typename Vector_::DataType(1))
Applies the bilinear operator to an FE function.
Definition: bilinear_operator_assembler.hpp:390
FEAT::Assembly::BilinearOperatorAssembler::assemble_matrix1
static void assemble_matrix1(Matrix_ &matrix, Operator_ &operat, const Space_ &space, const CubatureFactory_ &cubature_factory, typename Matrix_::DataType alpha=typename Matrix_::DataType(1))
Assembles a bilinear operator into a matrix.
Definition: bilinear_operator_assembler.hpp:225
likwid_marker.hpp
FEAT::Tiny::axpy
CUDA_HOST_DEVICE void axpy(T_ &y, const T_ &x, const T_ &alpha)
Performs an AXPY of two scalars.
Definition: tiny_algebra.hpp:3072
FEAT
FEAT namespace.
Definition: adjactor.hpp:12