FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
stokes_3field.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
9#include <kernel/lafem/dense_vector.hpp>
10#include <kernel/lafem/dense_vector_blocked.hpp>
11#include <kernel/lafem/filter_chain.hpp>
12#include <kernel/lafem/sparse_matrix_csr.hpp>
13#include <kernel/lafem/sparse_matrix_bcsr.hpp>
14#include <kernel/lafem/sparse_matrix_bwrappedcsr.hpp>
15#include <kernel/lafem/null_matrix.hpp>
16#include <kernel/lafem/tuple_matrix.hpp>
17#include <kernel/lafem/unit_filter.hpp>
18#include <kernel/lafem/unit_filter_blocked.hpp>
19#include <kernel/lafem/mean_filter.hpp>
20#include <kernel/lafem/none_filter.hpp>
21#include <kernel/lafem/slip_filter.hpp>
22#include <kernel/lafem/tuple_filter.hpp>
23#include <kernel/lafem/tuple_mirror.hpp>
24#include <kernel/lafem/tuple_diag_matrix.hpp>
25#include <kernel/lafem/transfer.hpp>
26#include <kernel/assembly/mirror_assembler.hpp>
27#include <kernel/assembly/symbolic_assembler.hpp>
28#include <kernel/assembly/bilinear_operator_assembler.hpp>
29#include <kernel/assembly/domain_assembler_basic_jobs.hpp>
30#include <kernel/assembly/common_operators.hpp>
31#include <kernel/assembly/gpdv_assembler.hpp>
32#include <kernel/assembly/grid_transfer.hpp>
33#include <kernel/assembly/mean_filter_assembler.hpp>
34#include <kernel/assembly/unit_filter_assembler.hpp>
35#include <kernel/global/gate.hpp>
36#include <kernel/global/muxer.hpp>
37#include <kernel/global/vector.hpp>
38#include <kernel/global/matrix.hpp>
39#include <kernel/global/filter.hpp>
40#include <kernel/global/mean_filter.hpp>
41#include <kernel/global/transfer.hpp>
42
43#include <control/domain/domain_control.hpp>
44#include <control/asm/gate_asm.hpp>
45#include <control/asm/muxer_asm.hpp>
46#include <control/asm/splitter_asm.hpp>
47#include <control/asm/transfer_asm.hpp>
48#include <control/asm/transfer_voxel_asm.hpp>
49#include <control/asm/unit_filter_asm.hpp>
50#include <control/asm/mean_filter_asm.hpp>
51#include <control/asm/slip_filter_asm.hpp>
52
53namespace FEAT
54{
55 namespace Control
56 {
64 template
65 <
66 int dim_,
67 int nsc_ = (dim_*(dim_+1))/2,
68 typename DataType_ = Real,
69 typename IndexType_ = Index,
70 typename MatrixBlockA_ = LAFEM::SparseMatrixBCSR<DataType_, IndexType_, dim_, dim_>,
71 typename MatrixBlockB_ = LAFEM::SparseMatrixBCSR<DataType_, IndexType_, dim_, 1>,
72 typename MatrixBlockD_ = LAFEM::SparseMatrixBCSR<DataType_, IndexType_, 1, dim_>,
73 typename MatrixBlockM_ = LAFEM::SparseMatrixBCSR<DataType_, IndexType_, nsc_, nsc_>,
74 typename MatrixBlockK_ = LAFEM::SparseMatrixBCSR<DataType_, IndexType_, nsc_, dim_>,
75 typename MatrixBlockL_ = LAFEM::SparseMatrixBCSR<DataType_, IndexType_, dim_, nsc_>,
76 typename ScalarMatrix_ = LAFEM::SparseMatrixCSR<DataType_, IndexType_>,
77 typename TransferMatrixV_ = LAFEM::SparseMatrixBWrappedCSR<DataType_, IndexType_, dim_>,
78 typename TransferMatrixP_ = LAFEM::SparseMatrixCSR<DataType_, IndexType_>,
79 typename TransferMatrixS_ = LAFEM::SparseMatrixBWrappedCSR<DataType_, IndexType_, nsc_>
80 >
82 {
83 public:
84 // basic types
85 typedef DataType_ DataType;
86 typedef IndexType_ IndexType;
87 static constexpr int dim = dim_;
88 static constexpr int nsc = nsc_;
89
90 // scalar types
91 typedef ScalarMatrix_ LocalScalarMatrix;
92 typedef typename LocalScalarMatrix::VectorTypeL LocalScalarVector;
93
94 // define local blocked matrix types
95 typedef MatrixBlockA_ LocalMatrixBlockA;
96 typedef MatrixBlockB_ LocalMatrixBlockB;
97 typedef MatrixBlockD_ LocalMatrixBlockD;
98 typedef MatrixBlockM_ LocalMatrixBlockM;
99 typedef MatrixBlockK_ LocalMatrixBlockK;
100 typedef MatrixBlockL_ LocalMatrixBlockL;
104 typedef LAFEM::TupleMatrix<
109
110 // define local vector types
111 typedef typename LocalMatrixBlockB::VectorTypeL LocalVeloVector;
112 typedef typename LocalMatrixBlockD::VectorTypeL LocalPresVector;
113 typedef typename LocalMatrixBlockM::VectorTypeL LocalStressVector;
115
116 // define local transfer matrix types
117 typedef TransferMatrixV_ LocalVeloTransferMatrix;
118 typedef TransferMatrixP_ LocalPresTransferMatrix;
119 typedef TransferMatrixS_ LocalStressTransferMatrix;
121
122 // define local transfer operators
127
128 // define mirror types
130 typedef ScalarMirror VeloMirror;
131 typedef ScalarMirror PresMirror;
134
135 // define gates
141
142 // define muxers
148
149 // define global vector types
154
155 // define global matrix types
163
164 // define global transfer types
169
170 /* ***************************************************************************************** */
171
174 PresGate gate_pres;
175 StressGate gate_stress;
176 SystemGate gate_sys;
177 ScalarGate gate_scalar_velo;
178 ScalarGate gate_scalar_stress;
179
182 PresMuxer coarse_muxer_pres;
183 StressMuxer coarse_muxer_stress;
184 SystemMuxer coarse_muxer_sys;
185 ScalarMuxer coarse_muxer_scalar_velo;
186 ScalarMuxer coarse_muxer_scalar_stress;
187
190 GlobalMatrixBlockA matrix_a;
191 GlobalMatrixBlockB matrix_b;
192 GlobalMatrixBlockD matrix_d;
193 GlobalMatrixBlockM matrix_m;
194 GlobalMatrixBlockK matrix_k;
195 GlobalMatrixBlockL matrix_l;
196
199 NullMatrixBlockSP null_matrix_sp;
200 NullMatrixBlockPS null_matrix_ps;
201
204 GlobalPresTransfer transfer_pres;
205 GlobalStressTransfer transfer_stress;
206 GlobalSystemTransfer transfer_sys;
207
210 matrix_sys(&gate_sys, &gate_sys),
211 matrix_a(&gate_velo, &gate_velo),
212 matrix_b(&gate_velo, &gate_pres),
213 matrix_d(&gate_pres, &gate_velo),
214 matrix_m(&gate_stress, &gate_stress),
215 matrix_k(&gate_stress, &gate_velo),
216 matrix_l(&gate_velo, &gate_stress),
218 transfer_pres(&coarse_muxer_pres),
219 transfer_stress(&coarse_muxer_stress),
220 transfer_sys(&coarse_muxer_sys)
221 {
222 }
223
224 // no copies, no problems
226 Stokes3FieldSystemLevel& operator=(const Stokes3FieldSystemLevel&) = delete;
227
229 {
230 }
231
233 std::size_t bytes() const
234 {
235 return this->matrix_sys.local().bytes () + transfer_sys.bytes();
236 }
237
238 void compile_system_matrix()
239 {
240 matrix_sys.local().template at<0,0>() = matrix_a.local().clone(LAFEM::CloneMode::Shallow);
241 matrix_sys.local().template at<0,1>() = matrix_b.local().clone(LAFEM::CloneMode::Shallow);
242 matrix_sys.local().template at<0,2>() = matrix_l.local().clone(LAFEM::CloneMode::Shallow);
243 matrix_sys.local().template at<1,0>() = matrix_d.local().clone(LAFEM::CloneMode::Shallow);
245 matrix_sys.local().template at<1,2>() = null_matrix_ps.clone(LAFEM::CloneMode::Shallow);
246 matrix_sys.local().template at<2,0>() = matrix_k.local().clone(LAFEM::CloneMode::Shallow);
247 matrix_sys.local().template at<2,1>() = null_matrix_sp.clone(LAFEM::CloneMode::Shallow);
248 matrix_sys.local().template at<2,2>() = matrix_m.local().clone(LAFEM::CloneMode::Shallow);
249 }
250
251 void compile_system_transfer()
252 {
253 transfer_sys.get_mat_prol().template at<0,0>() = transfer_velo.get_mat_prol().clone(LAFEM::CloneMode::Shallow);
254 transfer_sys.get_mat_rest().template at<0,0>() = transfer_velo.get_mat_rest().clone(LAFEM::CloneMode::Shallow);
255 transfer_sys.get_mat_trunc().template at<0,0>() = transfer_velo.get_mat_trunc().clone(LAFEM::CloneMode::Shallow);
256 transfer_sys.get_mat_prol().template at<1,1>() = transfer_pres.get_mat_prol().clone(LAFEM::CloneMode::Shallow);
257 transfer_sys.get_mat_rest().template at<1,1>() = transfer_pres.get_mat_rest().clone(LAFEM::CloneMode::Shallow);
258 transfer_sys.get_mat_trunc().template at<1,1>() = transfer_pres.get_mat_trunc().clone(LAFEM::CloneMode::Shallow);
259 transfer_sys.get_mat_prol().template at<2,2>() = transfer_stress.get_mat_prol().clone(LAFEM::CloneMode::Shallow);
260 transfer_sys.get_mat_rest().template at<2,2>() = transfer_stress.get_mat_rest().clone(LAFEM::CloneMode::Shallow);
261 transfer_sys.get_mat_trunc().template at<2,2>() = transfer_stress.get_mat_trunc().clone(LAFEM::CloneMode::Shallow);
262 transfer_sys.compile();
263 }
264
265 template<typename DomainLevel_>
266 void assemble_gates(const Domain::VirtualLevel<DomainLevel_>& virt_dom_lvl)
267 {
268 Asm::asm_gate(virt_dom_lvl, virt_dom_lvl->space_velo, this->gate_velo, true);
269 Asm::asm_gate(virt_dom_lvl, virt_dom_lvl->space_pres, this->gate_pres, true);
270 Asm::asm_gate(virt_dom_lvl, virt_dom_lvl->space_stress, this->gate_stress, true);
271 Asm::build_gate_tuple(this->gate_sys, this->gate_velo, this->gate_pres, this->gate_stress);
272 this->gate_scalar_velo.convert(this->gate_velo, LocalScalarVector(virt_dom_lvl->space_velo.get_num_dofs()));
273 this->gate_scalar_stress.convert(this->gate_stress, LocalScalarVector(virt_dom_lvl->space_stress.get_num_dofs()));
274 }
275
276 template<typename DomainLevel_>
277 void assemble_coarse_muxers(const Domain::VirtualLevel<DomainLevel_>& virt_lvl_coarse)
278 {
279 Asm::asm_muxer(virt_lvl_coarse, [](const DomainLevel_& dl){return &dl.space_velo;}, this->coarse_muxer_velo);
280 Asm::asm_muxer(virt_lvl_coarse, [](const DomainLevel_& dl){return &dl.space_pres;}, this->coarse_muxer_pres);
281 Asm::asm_muxer(virt_lvl_coarse, [](const DomainLevel_& dl){return &dl.space_stress;}, this->coarse_muxer_stress);
282 Asm::build_muxer_tuple(this->coarse_muxer_sys, this->gate_sys.get_freqs(), this->coarse_muxer_velo, this->coarse_muxer_pres, this->coarse_muxer_stress);
283 this->coarse_muxer_scalar_velo.convert(this->coarse_muxer_velo, this->gate_scalar_velo.get_freqs().clone(LAFEM::CloneMode::Shallow));
284 this->coarse_muxer_scalar_stress.convert(this->coarse_muxer_stress, this->gate_scalar_stress.get_freqs().clone(LAFEM::CloneMode::Shallow));
285 }
286
287 template<typename DomainLevel_>
288 void assemble_transfers(
289 const Stokes3FieldSystemLevel& sys_lvl_coarse,
290 const Domain::VirtualLevel<DomainLevel_>& virt_lvl_fine,
291 const Domain::VirtualLevel<DomainLevel_>& virt_lvl_coarse,
292 const String& cubature, bool trunc_v = false, bool trunc_p = false, bool trunc_s = false, bool shrink = true)
293 {
294 Asm::asm_transfer_blocked(virt_lvl_fine, virt_lvl_coarse, cubature, trunc_v, shrink,
295 [](const DomainLevel_& dl) {return &dl.space_velo;},
296 this->transfer_velo.local(), this->coarse_muxer_velo, this->gate_velo, sys_lvl_coarse.gate_velo);
297 Asm::asm_transfer_scalar(virt_lvl_fine, virt_lvl_coarse, cubature, trunc_p, shrink,
298 [](const DomainLevel_& dl) {return &dl.space_pres;},
299 this->transfer_pres.local(), this->coarse_muxer_pres, this->gate_pres, sys_lvl_coarse.gate_pres);
300 Asm::asm_transfer_blocked(virt_lvl_fine, virt_lvl_coarse, cubature, trunc_s, shrink,
301 [](const DomainLevel_& dl) {return &dl.space_stress;},
302 this->transfer_stress.local(), this->coarse_muxer_stress, this->gate_stress, sys_lvl_coarse.gate_stress);
303
304 this->transfer_velo.compile();
305 this->transfer_pres.compile();
306 this->transfer_stress.compile();
307 this->compile_system_transfer();
308 }
309
310 template<typename DomainLevel_>
311 void assemble_transfers(
312 const Domain::VirtualLevel<DomainLevel_>& virt_lvl_fine,
313 const Domain::VirtualLevel<DomainLevel_>& virt_lvl_coarse,
314 const String& cubature, bool trunc_v = false, bool trunc_p = false, bool trunc_s = false, bool shrink = true)
315 {
316 // if the coarse level is a parent, then we need the coarse system level
317 XASSERT(!virt_lvl_coarse.is_parent());
318
319 Asm::asm_transfer_blocked(virt_lvl_fine, virt_lvl_coarse, cubature, trunc_v, shrink,
320 [](const DomainLevel_& dl) {return &dl.space_velo;},
321 this->transfer_velo.local(), this->coarse_muxer_velo, this->gate_velo, this->gate_velo);
322 Asm::asm_transfer_scalar(virt_lvl_fine, virt_lvl_coarse, cubature, trunc_p, shrink,
323 [](const DomainLevel_& dl) {return &dl.space_pres;},
324 this->transfer_pres.local(), this->coarse_muxer_pres, this->gate_pres, this->gate_pres);
325 Asm::asm_transfer_blocked(virt_lvl_fine, virt_lvl_coarse, cubature, trunc_s, shrink,
326 [](const DomainLevel_& dl) {return &dl.space_stress;},
327 this->transfer_stress.local(), this->coarse_muxer_stress, this->gate_stress, this->gate_stress);
328
329 this->transfer_velo.compile();
330 this->transfer_pres.compile();
331 this->transfer_stress.compile();
332 this->compile_system_transfer();
333 }
334
335 template<typename SpaceVelo_, typename SpacePres_, typename Cubature_>
336 void assemble_grad_div_matrices(const SpaceVelo_& space_velo, const SpacePres_& space_pres, const Cubature_& cubature)
337 {
338 Assembly::GradPresDivVeloAssembler::assemble(this->matrix_b.local(), this->matrix_d.local(), space_velo, space_pres, cubature);
339 }
340
341 template<typename SpaceVelo_, typename SpacePres_, typename SpaceStress_>
342 void assemble_structs(const SpaceVelo_& space_velo, const SpacePres_& space_pres, const SpaceStress_& space_stress)
343 {
344 Assembly::SymbolicAssembler::assemble_matrix_std1(this->matrix_a.local(), space_velo);
345 Assembly::SymbolicAssembler::assemble_matrix_std2(this->matrix_b.local(), space_velo, space_pres);
346 Assembly::SymbolicAssembler::assemble_matrix_std2(this->matrix_d.local(), space_pres, space_velo);
347 Assembly::SymbolicAssembler::assemble_matrix_std1(this->matrix_m.local(), space_stress);
348 Assembly::SymbolicAssembler::assemble_matrix_std2(this->matrix_l.local(), space_velo, space_stress);
349 Assembly::SymbolicAssembler::assemble_matrix_std2(this->matrix_k.local(), space_stress, space_velo);
350 Assembly::SymbolicAssembler::assemble_matrix_std1(this->null_matrix_pp, space_pres);
351 Assembly::SymbolicAssembler::assemble_matrix_std2(this->null_matrix_ps, space_pres, space_stress);
352 Assembly::SymbolicAssembler::assemble_matrix_std2(this->null_matrix_sp, space_stress, space_pres);
353 }
354 }; // class Stokes3FieldSystemLevel<...>
355 } // namespace Control
356} // namespace FEAT
#define XASSERT(expr)
Assertion macro definition.
Definition: assertion.hpp:262
FEAT Kernel base header.
static void assemble(LAFEM::SparseMatrixBCSR< DataType_, IndexType_, dim_, 1 > &matrix_b, LAFEM::SparseMatrixBCSR< DataType_, IndexType_, 1, dim_ > &matrix_d, const SpaceVelo_ &space_velo, const SpacePres_ &space_pres, const String &cubature_name, const DataType_ scale_b=-DataType_(1), const DataType_ scale_d=-DataType_(1))
Assembles the B and D matrices.
static void assemble_matrix_std1(MatrixType_ &matrix, const Space_ &space)
Assembles a standard matrix structure from a single space.
static void assemble_matrix_std2(MatrixType_ &matrix, const TestSpace_ &test_space, const TrialSpace_ &trial_space)
Assembles a standard matrix structure from a test-trial-space pair.
GlobalVeloTransfer transfer_velo
our global transfer operator
std::size_t bytes() const
Returns the total amount of bytes allocated.
NullMatrixBlockPP null_matrix_pp
null-matrix blocks
GlobalSystemMatrix matrix_sys
our global system matrix
VeloMuxer coarse_muxer_velo
our coarse-level system muxer
Global gate implementation.
Definition: gate.hpp:51
void convert(const Gate< LVT2_, MT2_ > &other)
Conversion function for same vector container type but with different MDI-Type.
Definition: gate.hpp:191
const LocalVector_ & get_freqs() const
Returns a const reference to the frequencies vector.
Definition: gate.hpp:179
Global Matrix wrapper class template.
Definition: matrix.hpp:40
LocalMatrix_ & local()
Returns a reference to the internal local LAFEM matrix object.
Definition: matrix.hpp:126
Global multiplexer/demultiplexer implementation.
Definition: muxer.hpp:68
void convert(const Muxer< LVT2_, MT2_ > &other)
Conversion function for same vector container type but with different MDI-Type.
Definition: muxer.hpp:155
Global grid-transfer operator class template.
Definition: transfer.hpp:23
Transfer clone(LAFEM::CloneMode clone_mode=LAFEM::CloneMode::Weak) const
Creates a clone of this object.
Definition: transfer.hpp:125
std::size_t bytes() const
Definition: transfer.hpp:143
LocalTransfer_ & local()
Definition: transfer.hpp:131
Global vector wrapper class template.
Definition: vector.hpp:68
Null Matrix implementation.
Definition: null_matrix.hpp:56
NullMatrix clone(CloneMode clone_mode=CloneMode::Weak) const
Clone operation.
Grid-Transfer operator class template.
Definition: transfer.hpp:27
Tuple-Diag-Matrix meta class template.
Variadic TupleMatrix class template.
std::size_t bytes() const
Returns the total amount of bytes allocated.
TupleMatrix row helper class template.
TupleVector meta-mirror class template.
Variadic TupleVector class template.
Handles vector prolongation, restriction and serialization.
FEAT namespace.
Definition: adjactor.hpp:12
double Real
Real data type.
std::uint64_t Index
Index data type.