Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #ifdef MOOSE_MFEM_ENABLED
11 :
12 : #include "EquationSystem.h"
13 : #include "libmesh/int_range.h"
14 :
15 : namespace Moose::MFEM
16 : {
17 :
18 1361 : EquationSystem::~EquationSystem()
19 : {
20 1361 : DeleteHBlocks();
21 1361 : DeleteJacobianBlocks();
22 1361 : }
23 :
24 : void
25 3152 : EquationSystem::DeleteHBlocks()
26 : {
27 4974 : for (const auto i : make_range(_h_blocks.NumRows()))
28 3706 : for (const auto j : make_range(_h_blocks.NumCols()))
29 : {
30 1884 : if (_jacobian_blocks.NumRows() && _jacobian_blocks(i, j) == _h_blocks(i, j))
31 0 : _jacobian_blocks(i, j) = nullptr;
32 1884 : delete _h_blocks(i, j);
33 : }
34 3152 : _h_blocks.DeleteAll();
35 3152 : }
36 :
37 : void
38 3541 : EquationSystem::DeleteJacobianBlocks()
39 : {
40 5721 : for (const auto i : make_range(_jacobian_blocks.NumRows()))
41 4360 : for (const auto j : make_range(_jacobian_blocks.NumCols()))
42 2180 : if (!_h_blocks.NumRows() || _jacobian_blocks(i, j) != _h_blocks(i, j))
43 2180 : delete _jacobian_blocks(i, j);
44 3541 : _jacobian_blocks.DeleteAll();
45 3541 : }
46 :
47 : bool
48 6752 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
49 : const std::string & name) const
50 : {
51 6752 : return std::find(the_vector.begin(), the_vector.end(), name) != the_vector.end();
52 : }
53 :
54 : void
55 1620 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
56 : {
57 1620 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
58 1019 : _coupled_var_names.push_back(coupled_var_name);
59 1620 : }
60 :
61 : void
62 145 : EquationSystem::AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name)
63 : {
64 145 : if (!VectorContainsName(_eliminated_var_names, eliminated_var_name))
65 139 : _eliminated_var_names.push_back(eliminated_var_name);
66 145 : }
67 :
68 : void
69 2913 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
70 : {
71 2913 : if (!VectorContainsName(_test_var_names, test_var_name))
72 937 : _test_var_names.push_back(test_var_name);
73 2913 : }
74 :
75 : void
76 1361 : EquationSystem::SetTrialVariableNames()
77 : {
78 : // If a coupled variable has an equation associated with it,
79 : // add it to the set of trial variables.
80 2298 : for (const auto & test_var_name : _test_var_names)
81 937 : if (VectorContainsName(_coupled_var_names, test_var_name))
82 937 : _trial_var_names.push_back(test_var_name);
83 :
84 : // Otherwise, add it to the set of eliminated variables.
85 2380 : for (const auto & coupled_var_name : _coupled_var_names)
86 1019 : if (!VectorContainsName(_test_var_names, coupled_var_name))
87 82 : _eliminated_var_names.push_back(coupled_var_name);
88 1361 : }
89 :
90 : void
91 1491 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
92 : {
93 1491 : const auto & trial_var_name = kernel->getTrialVariableName();
94 1491 : const auto & test_var_name = kernel->getTestVariableName();
95 1491 : AddCoupledVariableNameIfMissing(trial_var_name);
96 1491 : AddTestVariableNameIfMissing(test_var_name);
97 : // Register new kernels map if not present for the test variable
98 1491 : if (!_kernels_map.Has(test_var_name))
99 : {
100 : auto kernel_field_map =
101 890 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
102 890 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
103 890 : }
104 : // Register new kernels map if not present for the test/trial variable pair
105 1491 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
106 : {
107 985 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
108 985 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
109 985 : }
110 1491 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
111 1491 : }
112 :
113 : void
114 56 : EquationSystem::AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> bc)
115 : {
116 56 : const auto & trial_var_name = bc->getTrialVariableName();
117 56 : const auto & test_var_name = bc->getTestVariableName();
118 56 : AddCoupledVariableNameIfMissing(trial_var_name);
119 56 : AddTestVariableNameIfMissing(test_var_name);
120 : // Register new integrated bc map if not present for the test variable
121 56 : if (!_integrated_bc_map.Has(test_var_name))
122 : {
123 : auto integrated_bc_field_map = std::make_shared<
124 50 : Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>();
125 50 : _integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
126 50 : }
127 : // Register new integrated bc map if not present for the test/trial variable pair
128 56 : if (!_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
129 : {
130 50 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMIntegratedBC>>>();
131 50 : _integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
132 50 : }
133 56 : _integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
134 56 : }
135 :
136 : void
137 1086 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
138 : {
139 1086 : const auto & test_var_name = bc->getTestVariableName();
140 1086 : AddTestVariableNameIfMissing(test_var_name);
141 : // Register new essential bc map if not present for the test variable
142 1086 : if (!_essential_bc_map.Has(test_var_name))
143 : {
144 814 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
145 814 : _essential_bc_map.Register(test_var_name, std::move(bcs));
146 814 : }
147 1086 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
148 1086 : }
149 :
150 : void
151 1299 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
152 : Moose::MFEM::ComplexGridFunctions & cmplx_gridfunctions,
153 : mfem::AssemblyLevel assembly_level)
154 : {
155 1299 : _assembly_level = assembly_level;
156 :
157 1299 : if (cmplx_gridfunctions.size())
158 0 : mooseError("Complex variables have been created but the executioner numeric type has not been "
159 : "set to complex. Please set Executioner/numeric_type = complex.");
160 :
161 : // Extract which coupled variables are to be trivially eliminated and which are trial variables
162 1299 : SetTrialVariableNames();
163 :
164 2195 : for (auto & test_var_name : _test_var_names)
165 : {
166 896 : if (!gridfunctions.Has(test_var_name))
167 : {
168 0 : mooseError("MFEM variable ",
169 : test_var_name,
170 : " requested by equation system during initialization was "
171 : "not found in gridfunctions");
172 : }
173 : // Store pointers to test FESpaces
174 896 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
175 : }
176 :
177 2195 : for (auto & trial_var_name : _trial_var_names)
178 : {
179 896 : if (!gridfunctions.Has(trial_var_name))
180 : {
181 0 : mooseError("MFEM variable ",
182 : trial_var_name,
183 : " requested by equation system during initialization was "
184 : "not found in gridfunctions");
185 : }
186 : // Create auxiliary gridfunctions for storing essential constraints from Dirichlet conditions
187 896 : _var_ess_constraints.emplace_back(
188 1792 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(trial_var_name)->ParFESpace()));
189 : }
190 :
191 : // Store pointers to FESpaces of all coupled variables
192 2277 : for (auto & coupled_var_name : _coupled_var_names)
193 978 : _coupled_pfespaces.push_back(gridfunctions.Get(coupled_var_name)->ParFESpace());
194 :
195 : // Store pointers to coupled variable GridFunctions that are to be eliminated prior to forming the
196 : // jacobian
197 1520 : for (auto & eliminated_var_name : _eliminated_var_names)
198 221 : _eliminated_variables.Register(eliminated_var_name,
199 442 : gridfunctions.GetShared(eliminated_var_name));
200 :
201 : // Get a reference to the GridFunctions
202 1299 : _gfuncs = &gridfunctions;
203 1299 : }
204 :
205 : void
206 1817 : EquationSystem::ApplyEssentialBC(const std::string & var_name,
207 : mfem::ParGridFunction & trial_gf,
208 : mfem::Array<int> & global_ess_markers)
209 : {
210 1817 : if (_essential_bc_map.Has(var_name))
211 : {
212 1405 : auto & bcs = _essential_bc_map.GetRef(var_name);
213 3627 : for (auto & bc : bcs)
214 : {
215 : // Set constrained DoFs values on essential boundaries
216 2222 : bc->ApplyBC(trial_gf);
217 : // Fetch marker array labelling essential boundaries of current BC
218 2222 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
219 : // Add these boundary markers to the set of markers labelling all essential boundaries
220 11060 : for (const auto i : make_range(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max()))
221 8838 : global_ess_markers[i] = std::max(global_ess_markers[i], ess_bdrs[i]);
222 2222 : }
223 : }
224 1817 : }
225 :
226 : void
227 1786 : EquationSystem::ApplyEssentialBCs()
228 : {
229 1786 : _ess_tdof_lists.resize(_trial_var_names.size());
230 3603 : for (const auto i : index_range(_trial_var_names))
231 : {
232 1817 : const auto & trial_var_name = _trial_var_names.at(i);
233 1817 : mfem::ParGridFunction & trial_gf = *_var_ess_constraints.at(i);
234 :
235 : // Make sure we update the size, if this mesh has changed recently for instance
236 1817 : trial_gf.Update();
237 :
238 : // Initial guess for non-linear problems (initial condition or the previous time step solution)
239 1817 : trial_gf = _gfuncs->GetRef(trial_var_name);
240 :
241 1817 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
242 1817 : global_ess_markers = 0;
243 : // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
244 : // essential boundaries to the global_ess_markers array
245 1817 : ApplyEssentialBC(trial_var_name, trial_gf, global_ess_markers);
246 1817 : trial_gf.ParFESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
247 1817 : }
248 1786 : }
249 :
250 : void
251 1812 : EquationSystem::EliminateCoupledVariables()
252 : {
253 3655 : for (const auto & test_var_name : _test_var_names)
254 3021 : for (const auto & eliminated_var_name : _eliminated_var_names)
255 1296 : if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(eliminated_var_name) &&
256 118 : !VectorContainsName(_test_var_names, eliminated_var_name))
257 : {
258 82 : auto & mblf = *_mblfs.Get(test_var_name)->Get(eliminated_var_name);
259 82 : mblf.AddMult(*_eliminated_variables.Get(eliminated_var_name), *_lfs.Get(test_var_name), -1);
260 : }
261 1812 : }
262 :
263 : void
264 1827 : EquationSystem::FormLinearSystem(mfem::OperatorHandle & op,
265 : mfem::BlockVector & trueX,
266 : mfem::BlockVector & trueRHS)
267 : {
268 : mooseAssert(_test_var_names.size() == _trial_var_names.size(),
269 : "Number of test and trial variables must be the same for block matrix assembly.");
270 :
271 1827 : if (_assembly_level == mfem::AssemblyLevel::LEGACY)
272 1791 : FormSystemMatrix(op, trueX, trueRHS);
273 : else
274 : {
275 : mooseAssert(_test_var_names.size() == 1 && _test_var_names.size() == _trial_var_names.size(),
276 : "Non-legacy assembly is only supported for single test and trial variable systems");
277 36 : FormSystemOperator(op, trueX, trueRHS);
278 : }
279 1827 : }
280 :
281 : void
282 36 : EquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
283 : mfem::BlockVector & trueX,
284 : mfem::BlockVector & trueRHS)
285 : {
286 36 : auto & test_var_name = _test_var_names.at(0);
287 36 : mfem::Vector aux_x, aux_rhs;
288 36 : mfem::OperatorPtr aux_a;
289 :
290 36 : auto blf = _blfs.Get(test_var_name);
291 36 : blf->FormLinearSystem(_ess_tdof_lists.at(0),
292 36 : *_var_ess_constraints.at(0),
293 36 : *_lfs.Get(test_var_name),
294 : aux_a,
295 : aux_x,
296 : aux_rhs,
297 : /*copy_interior=*/true);
298 :
299 36 : trueX.GetBlock(0) = aux_x;
300 36 : trueRHS.GetBlock(0) = aux_rhs;
301 36 : trueX.SyncFromBlocks();
302 36 : trueRHS.SyncFromBlocks();
303 :
304 36 : op.Reset(aux_a.Ptr());
305 36 : aux_a.SetOperatorOwner(false);
306 36 : }
307 :
308 : void
309 1750 : EquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
310 : mfem::BlockVector & trueX,
311 : mfem::BlockVector & trueRHS)
312 : {
313 : // Allocate block operator
314 1750 : DeleteHBlocks();
315 1750 : _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
316 1750 : _h_blocks = nullptr;
317 : // Zero out RHS and sync memory
318 1750 : trueRHS = 0.0;
319 1750 : trueRHS.SyncToBlocks();
320 :
321 3531 : for (const auto i : index_range(_test_var_names))
322 : {
323 1781 : auto test_var_name = _test_var_names.at(i);
324 :
325 3624 : for (const auto j : index_range(_trial_var_names))
326 : {
327 1843 : auto trial_var_name = _trial_var_names.at(j);
328 :
329 1843 : mfem::Vector aux_x, aux_rhs;
330 1843 : mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
331 1843 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
332 :
333 1843 : if (test_var_name == trial_var_name)
334 : {
335 : mooseAssert(i == j, "Trial and test variables must have the same ordering.");
336 1781 : auto blf = _blfs.Get(test_var_name);
337 1781 : blf->FormLinearSystem(_ess_tdof_lists.at(j),
338 1781 : *_var_ess_constraints.at(j),
339 1781 : *_lfs.Get(test_var_name),
340 : *aux_a,
341 : aux_x,
342 : aux_rhs,
343 : /*copy_interior=*/true);
344 1781 : trueX.GetBlock(j) = aux_x;
345 : }
346 62 : else if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(trial_var_name))
347 : {
348 62 : auto mblf = _mblfs.Get(test_var_name)->Get(trial_var_name);
349 62 : mblf->FormRectangularLinearSystem(_ess_tdof_lists.at(j),
350 62 : _ess_tdof_lists.at(i),
351 62 : *_var_ess_constraints.at(j),
352 62 : aux_lf = 0,
353 : *aux_a,
354 : aux_x,
355 : aux_rhs);
356 : }
357 : else
358 0 : continue;
359 :
360 1843 : trueRHS.GetBlock(i) += aux_rhs;
361 1843 : _h_blocks(i, j) = aux_a;
362 1843 : }
363 1781 : }
364 : // Sync memory
365 1750 : trueX.SyncFromBlocks();
366 1750 : trueRHS.SyncFromBlocks();
367 :
368 : // Create monolithic matrix
369 1750 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
370 1750 : }
371 :
372 : void
373 1827 : EquationSystem::FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
374 : {
375 1827 : height = trueX.Size();
376 1827 : width = trueRHS.Size();
377 : // Store block offsets
378 1827 : _block_true_offsets.SetSize(trueX.NumBlocks() + 1);
379 1827 : _block_true_offsets[0] = 0;
380 3685 : for (unsigned i = 0; i < _trial_var_names.size(); i++)
381 1858 : _block_true_offsets[i + 1] = trueX.BlockSize(i);
382 1827 : _block_true_offsets.PartialSum();
383 1827 : FormLinearSystem(_linear_operator, trueX, trueRHS);
384 1827 : }
385 :
386 : void
387 5408 : EquationSystem::Mult(const mfem::Vector & sol, mfem::Vector & residual) const
388 : {
389 : // Update gridfunctions that may be referenced by coefficients within nonlinear integrators
390 5408 : const mfem::BlockVector blockSolution(const_cast<mfem::Vector &>(sol), _block_true_offsets);
391 5408 : SetTrialVariablesFromTrueVectors(blockSolution);
392 :
393 5408 : if (_non_linear)
394 : {
395 2524 : mfem::BlockVector blockResidual(residual, _block_true_offsets);
396 5048 : for (unsigned int i = 0; i < _test_var_names.size(); i++)
397 : {
398 2524 : auto & test_var_name = _test_var_names.at(i);
399 2524 : auto nlf = _nlfs.GetShared(test_var_name);
400 2524 : nlf->Mult(blockSolution.GetBlock(i), blockResidual.GetBlock(i));
401 2524 : blockResidual.GetBlock(i).SyncAliasMemory(blockResidual);
402 2524 : }
403 2524 : _linear_operator->AddMult(sol, residual);
404 2524 : }
405 : else
406 : {
407 2884 : residual = 0.0;
408 2884 : _linear_operator->Mult(sol, residual);
409 : }
410 :
411 5408 : sol.HostRead();
412 5408 : residual.HostRead();
413 5408 : }
414 :
415 : void
416 2180 : EquationSystem::FormJacobianMatrix(const mfem::Vector & u)
417 : {
418 2180 : DeleteJacobianBlocks();
419 2180 : _jacobian_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
420 2180 : _jacobian_blocks = nullptr;
421 :
422 2180 : const mfem::BlockVector update_vector(const_cast<mfem::Vector &>(u), _block_true_offsets);
423 4360 : for (const auto i : index_range(_test_var_names))
424 : {
425 2180 : auto test_var_name = _test_var_names.at(i);
426 2180 : if (_nlfs.Has(test_var_name))
427 : {
428 2180 : auto nlf = _nlfs.Get(test_var_name);
429 : mfem::HypreParMatrix * nlf_jac =
430 2180 : dynamic_cast<mfem::HypreParMatrix *>(&nlf->GetGradient(update_vector.GetBlock(i)));
431 : mooseAssert(nlf_jac,
432 : "Jacobian contribution of nonlinear form associated with " + test_var_name +
433 : " is not castable into a HypreParMatrix");
434 2180 : _jacobian_blocks(i, i) = mfem::ParAdd(_h_blocks(i, i), nlf_jac);
435 : }
436 : else
437 0 : _jacobian_blocks(i, i) = _h_blocks(i, i);
438 4360 : for (const auto j : index_range(_trial_var_names))
439 2180 : if (i != j) // nlf->GetGradient only contributes to on-diagonal blocks
440 0 : _jacobian_blocks(i, j) = _h_blocks(i, j);
441 2180 : }
442 : // Create monolithic matrix
443 2180 : _jacobian.Reset(mfem::HypreParMatrixFromBlocks(_jacobian_blocks));
444 2180 : }
445 :
446 : mfem::Operator &
447 3663 : EquationSystem::GetGradient(const mfem::Vector & u) const
448 : {
449 3663 : if (_non_linear)
450 2180 : const_cast<EquationSystem *>(this)->FormJacobianMatrix(u);
451 : else
452 1483 : _jacobian = _linear_operator;
453 :
454 3663 : return *_jacobian;
455 : }
456 :
457 : void
458 6146 : EquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
459 : {
460 12367 : for (const auto i : index_range(_trial_var_names))
461 : {
462 6221 : auto & trial_var_name = _trial_var_names.at(i);
463 6221 : trueX.GetBlock(i).SyncAliasMemory(trueX);
464 6221 : _gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
465 : }
466 6146 : }
467 :
468 : void
469 1812 : EquationSystem::BuildLinearForms()
470 : {
471 : // Register linear forms
472 3655 : for (const auto i : index_range(_test_var_names))
473 : {
474 1843 : auto test_var_name = _test_var_names.at(i);
475 1843 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
476 1843 : _lfs.GetRef(test_var_name) = 0.0;
477 1843 : }
478 :
479 3655 : for (auto & test_var_name : _test_var_names)
480 : {
481 : // Apply kernels
482 1843 : auto lf = _lfs.GetShared(test_var_name);
483 1843 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
484 1843 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
485 1843 : lf->Assemble();
486 1843 : }
487 :
488 : // Apply essential boundary conditions
489 1812 : ApplyEssentialBCs();
490 :
491 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
492 1812 : EliminateCoupledVariables();
493 1812 : }
494 :
495 : void
496 764 : EquationSystem::BuildNonlinearForms()
497 : {
498 : // Register non-linear Action forms
499 1541 : for (const auto i : index_range(_test_var_names))
500 : {
501 777 : auto test_var_name = _test_var_names.at(i);
502 777 : _nlfs.Register(test_var_name, std::make_shared<mfem::ParNonlinearForm>(_test_pfespaces.at(i)));
503 : // Apply kernels
504 777 : auto nlf = _nlfs.GetShared(test_var_name);
505 777 : nlf->SetEssentialTrueDofs(_ess_tdof_lists.at(i));
506 777 : ApplyDomainNLFIntegrators(test_var_name, nlf, _kernels_map);
507 777 : ApplyBoundaryNLFIntegrators(test_var_name, nlf, _integrated_bc_map);
508 777 : }
509 764 : }
510 :
511 : void
512 764 : EquationSystem::BuildBilinearForms()
513 : {
514 : // Register bilinear forms
515 1541 : for (const auto i : index_range(_test_var_names))
516 : {
517 777 : auto test_var_name = _test_var_names.at(i);
518 777 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
519 :
520 : // Apply kernels
521 777 : auto blf = _blfs.GetShared(test_var_name);
522 777 : blf->SetAssemblyLevel(_assembly_level);
523 1554 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
524 777 : test_var_name, test_var_name, blf, _integrated_bc_map);
525 1554 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
526 777 : test_var_name, test_var_name, blf, _kernels_map);
527 : // Assemble
528 777 : blf->Assemble();
529 777 : }
530 764 : }
531 :
532 : void
533 764 : EquationSystem::BuildMixedBilinearForms()
534 : {
535 : // Register mixed bilinear forms. Note that not all combinations may
536 : // have a kernel.
537 :
538 : // Create mblf for each test/coupled variable pair with an added kernel.
539 : // Mixed bilinear forms with coupled variables that are not trial variables are
540 : // associated with contributions from eliminated variables.
541 1541 : for (const auto i : index_range(_test_var_names))
542 : {
543 777 : auto test_var_name = _test_var_names.at(i);
544 777 : auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
545 1662 : for (const auto j : index_range(_coupled_var_names))
546 : {
547 885 : const auto & coupled_var_name = _coupled_var_names.at(j);
548 1770 : auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
549 885 : _test_pfespaces.at(i));
550 : // Register MixedBilinearForm if kernels exist for it, and assemble kernels
551 885 : if (_kernels_map.Has(test_var_name) &&
552 1757 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
553 872 : test_var_name != coupled_var_name)
554 : {
555 108 : mblf->SetAssemblyLevel(_assembly_level);
556 : // Apply all mixed kernels with this test/trial pair
557 216 : ApplyDomainBLFIntegrators<mfem::ParMixedBilinearForm>(
558 108 : coupled_var_name, test_var_name, mblf, _kernels_map);
559 : // Assemble mixed bilinear forms
560 108 : mblf->Assemble();
561 : // Register mixed bilinear forms associated with a single trial variable
562 : // for the current test variable
563 108 : test_mblfs->Register(coupled_var_name, mblf);
564 : }
565 885 : }
566 : // Register all mixed bilinear form sets associated with a single test variable
567 777 : _mblfs.Register(test_var_name, test_mblfs);
568 777 : }
569 764 : }
570 :
571 : void
572 1812 : EquationSystem::BuildEquationSystem()
573 : {
574 1812 : BuildBilinearForms();
575 1812 : BuildMixedBilinearForms();
576 1812 : BuildLinearForms();
577 1812 : BuildNonlinearForms();
578 1812 : }
579 :
580 : } // namespace Moose::MFEM
581 :
582 : #endif
|
