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 1372 : EquationSystem::~EquationSystem()
19 : {
20 1372 : DeleteHBlocks();
21 1372 : DeleteJacobianBlocks();
22 1372 : }
23 :
24 : void
25 3174 : EquationSystem::DeleteHBlocks()
26 : {
27 5018 : for (const auto i : make_range(_h_blocks.NumRows()))
28 3772 : for (const auto j : make_range(_h_blocks.NumCols()))
29 : {
30 1928 : if (_jacobian_blocks.NumRows() && _jacobian_blocks(i, j) == _h_blocks(i, j))
31 0 : _jacobian_blocks(i, j) = nullptr;
32 1928 : delete _h_blocks(i, j);
33 : }
34 3174 : _h_blocks.DeleteAll();
35 3174 : }
36 :
37 : void
38 3552 : EquationSystem::DeleteJacobianBlocks()
39 : {
40 5732 : 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 3552 : _jacobian_blocks.DeleteAll();
45 3552 : }
46 :
47 : bool
48 6895 : EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
49 : const std::string & name) const
50 : {
51 6895 : return std::find(the_vector.begin(), the_vector.end(), name) != the_vector.end();
52 : }
53 :
54 : void
55 1653 : EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
56 : {
57 1653 : if (!VectorContainsName(_coupled_var_names, coupled_var_name))
58 1041 : _coupled_var_names.push_back(coupled_var_name);
59 1653 : }
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 2979 : EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
70 : {
71 2979 : if (!VectorContainsName(_test_var_names, test_var_name))
72 959 : _test_var_names.push_back(test_var_name);
73 2979 : }
74 :
75 : void
76 1372 : EquationSystem::SetTrialVariableNames()
77 : {
78 : // If a coupled variable has an equation associated with it,
79 : // add it to the set of trial variables.
80 2331 : for (const auto & test_var_name : _test_var_names)
81 959 : if (VectorContainsName(_coupled_var_names, test_var_name))
82 959 : _trial_var_names.push_back(test_var_name);
83 :
84 : // Otherwise, add it to the set of eliminated variables.
85 2413 : for (const auto & coupled_var_name : _coupled_var_names)
86 1041 : if (!VectorContainsName(_test_var_names, coupled_var_name))
87 82 : _eliminated_var_names.push_back(coupled_var_name);
88 1372 : }
89 :
90 : void
91 1524 : EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
92 : {
93 1524 : const auto & trial_var_name = kernel->getTrialVariableName();
94 1524 : const auto & test_var_name = kernel->getTestVariableName();
95 1524 : AddCoupledVariableNameIfMissing(trial_var_name);
96 1524 : AddTestVariableNameIfMissing(test_var_name);
97 : // Register new kernels map if not present for the test variable
98 1524 : if (!_kernels_map.Has(test_var_name))
99 : {
100 : auto kernel_field_map =
101 912 : std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
102 912 : _kernels_map.Register(test_var_name, std::move(kernel_field_map));
103 912 : }
104 : // Register new kernels map if not present for the test/trial variable pair
105 1524 : if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
106 : {
107 1018 : auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
108 1018 : _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
109 1018 : }
110 1524 : _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
111 1524 : }
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 1119 : EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
138 : {
139 1119 : const auto & test_var_name = bc->getTestVariableName();
140 1119 : AddTestVariableNameIfMissing(test_var_name);
141 : // Register new essential bc map if not present for the test variable
142 1119 : if (!_essential_bc_map.Has(test_var_name))
143 : {
144 836 : auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
145 836 : _essential_bc_map.Register(test_var_name, std::move(bcs));
146 836 : }
147 1119 : _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
148 1119 : }
149 :
150 : void
151 1310 : EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
152 : Moose::MFEM::ComplexGridFunctions & cmplx_gridfunctions,
153 : mfem::AssemblyLevel assembly_level)
154 : {
155 1310 : _assembly_level = assembly_level;
156 :
157 1310 : 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 1310 : SetTrialVariableNames();
163 :
164 2228 : for (auto & test_var_name : _test_var_names)
165 : {
166 918 : 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 918 : _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
175 : }
176 :
177 2228 : for (auto & trial_var_name : _trial_var_names)
178 : {
179 918 : 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 918 : _var_ess_constraints.emplace_back(
188 1836 : std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(trial_var_name)->ParFESpace()));
189 : }
190 :
191 : // Store pointers to FESpaces of all coupled variables
192 2310 : for (auto & coupled_var_name : _coupled_var_names)
193 1000 : _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 1531 : 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 1310 : _gfuncs = &gridfunctions;
203 1310 : }
204 :
205 : void
206 1839 : EquationSystem::ApplyEssentialBC(const std::string & var_name,
207 : mfem::ParGridFunction & trial_gf,
208 : mfem::Array<int> & global_ess_markers)
209 : {
210 1839 : if (_essential_bc_map.Has(var_name))
211 : {
212 1427 : auto & bcs = _essential_bc_map.GetRef(var_name);
213 3682 : for (auto & bc : bcs)
214 : {
215 : // Set constrained DoFs values on essential boundaries
216 2255 : bc->ApplyBC(trial_gf);
217 : // Fetch marker array labelling essential boundaries of current BC
218 2255 : mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
219 : // Add these boundary markers to the set of markers labelling all essential boundaries
220 11214 : for (const auto i : make_range(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max()))
221 8959 : global_ess_markers[i] = std::max(global_ess_markers[i], ess_bdrs[i]);
222 2255 : }
223 : }
224 1839 : }
225 :
226 : void
227 1797 : EquationSystem::ApplyEssentialBCs()
228 : {
229 1797 : _ess_tdof_lists.resize(_trial_var_names.size());
230 3636 : for (const auto i : index_range(_trial_var_names))
231 : {
232 1839 : const auto & trial_var_name = _trial_var_names.at(i);
233 1839 : 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 1839 : trial_gf.Update();
237 :
238 : // Initial guess for non-linear problems (initial condition or the previous time step solution)
239 1839 : trial_gf = _gfuncs->GetRef(trial_var_name);
240 :
241 1839 : mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
242 1839 : 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 1839 : ApplyEssentialBC(trial_var_name, trial_gf, global_ess_markers);
246 1839 : trial_gf.ParFESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
247 1839 : }
248 1797 : }
249 :
250 : void
251 1823 : EquationSystem::EliminateCoupledVariables()
252 : {
253 3688 : for (const auto & test_var_name : _test_var_names)
254 3043 : 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 1823 : }
262 :
263 : void
264 1838 : 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 1838 : if (_assembly_level == mfem::AssemblyLevel::LEGACY)
272 1802 : 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 1838 : }
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 1761 : EquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
310 : mfem::BlockVector & trueX,
311 : mfem::BlockVector & trueRHS)
312 : {
313 : // Allocate block operator
314 1761 : DeleteHBlocks();
315 1761 : _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
316 1761 : _h_blocks = nullptr;
317 : // Zero out RHS and sync memory
318 1761 : trueRHS = 0.0;
319 1761 : trueRHS.SyncToBlocks();
320 :
321 3564 : for (const auto i : index_range(_test_var_names))
322 : {
323 1803 : auto test_var_name = _test_var_names.at(i);
324 :
325 3690 : for (const auto j : index_range(_trial_var_names))
326 : {
327 1887 : auto trial_var_name = _trial_var_names.at(j);
328 :
329 1887 : mfem::Vector aux_x, aux_rhs;
330 1887 : mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
331 1887 : mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
332 :
333 1887 : if (test_var_name == trial_var_name)
334 : {
335 : mooseAssert(i == j, "Trial and test variables must have the same ordering.");
336 1803 : auto blf = _blfs.Get(test_var_name);
337 1803 : blf->FormLinearSystem(_ess_tdof_lists.at(j),
338 1803 : *_var_ess_constraints.at(j),
339 1803 : *_lfs.Get(test_var_name),
340 : *aux_a,
341 : aux_x,
342 : aux_rhs,
343 : /*copy_interior=*/true);
344 1803 : trueX.GetBlock(j) = aux_x;
345 : }
346 84 : else if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(trial_var_name))
347 : {
348 73 : auto mblf = _mblfs.Get(test_var_name)->Get(trial_var_name);
349 73 : mblf->FormRectangularLinearSystem(_ess_tdof_lists.at(j),
350 73 : _ess_tdof_lists.at(i),
351 73 : *_var_ess_constraints.at(j),
352 73 : aux_lf = 0,
353 : *aux_a,
354 : aux_x,
355 : aux_rhs);
356 : }
357 : else
358 11 : continue;
359 :
360 1876 : trueRHS.GetBlock(i) += aux_rhs;
361 1876 : _h_blocks(i, j) = aux_a;
362 1920 : }
363 1803 : }
364 : // Sync memory
365 1761 : trueX.SyncFromBlocks();
366 1761 : trueRHS.SyncFromBlocks();
367 :
368 : // Create monolithic matrix
369 1761 : op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
370 1761 : }
371 :
372 : void
373 1838 : EquationSystem::FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
374 : {
375 1838 : height = trueX.Size();
376 1838 : width = trueRHS.Size();
377 : // Store block offsets
378 1838 : _block_true_offsets.SetSize(trueX.NumBlocks() + 1);
379 1838 : _block_true_offsets[0] = 0;
380 3718 : for (unsigned i = 0; i < _trial_var_names.size(); i++)
381 1880 : _block_true_offsets[i + 1] = trueX.BlockSize(i);
382 1838 : _block_true_offsets.PartialSum();
383 1838 : FormLinearSystem(_linear_operator, trueX, trueRHS);
384 1838 : }
385 :
386 : void
387 5430 : EquationSystem::Mult(const mfem::Vector & sol, mfem::Vector & residual) const
388 : {
389 : // Update gridfunctions that may be referenced by coefficients within nonlinear integrators
390 5430 : const mfem::BlockVector blockSolution(const_cast<mfem::Vector &>(sol), _block_true_offsets);
391 5430 : SetTrialVariablesFromTrueVectors(blockSolution);
392 :
393 5430 : 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 2906 : residual = 0.0;
408 2906 : _linear_operator->Mult(sol, residual);
409 : }
410 :
411 5430 : sol.HostRead();
412 5430 : residual.HostRead();
413 5430 : }
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 3674 : EquationSystem::GetGradient(const mfem::Vector & u) const
448 : {
449 3674 : if (_non_linear)
450 2180 : const_cast<EquationSystem *>(this)->FormJacobianMatrix(u);
451 : else
452 1494 : _jacobian = _linear_operator;
453 :
454 3674 : return *_jacobian;
455 : }
456 :
457 : void
458 6179 : EquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
459 : {
460 12466 : for (const auto i : index_range(_trial_var_names))
461 : {
462 6287 : auto & trial_var_name = _trial_var_names.at(i);
463 6287 : trueX.GetBlock(i).SyncAliasMemory(trueX);
464 6287 : _gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
465 : }
466 6179 : }
467 :
468 : void
469 1823 : EquationSystem::BuildLinearForms()
470 : {
471 : // Register linear forms
472 3688 : for (const auto i : index_range(_test_var_names))
473 : {
474 1865 : auto test_var_name = _test_var_names.at(i);
475 1865 : _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
476 1865 : _lfs.GetRef(test_var_name) = 0.0;
477 1865 : }
478 :
479 3688 : for (auto & test_var_name : _test_var_names)
480 : {
481 : // Apply kernels
482 1865 : auto lf = _lfs.GetShared(test_var_name);
483 1865 : ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
484 1865 : ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
485 1865 : lf->Assemble();
486 1865 : }
487 :
488 : // Apply essential boundary conditions
489 1823 : ApplyEssentialBCs();
490 :
491 : // Eliminate trivially eliminated variables by subtracting contributions from linear forms
492 1823 : EliminateCoupledVariables();
493 1823 : }
494 :
495 : void
496 775 : EquationSystem::BuildNonlinearForms()
497 : {
498 : // Register non-linear Action forms
499 1574 : for (const auto i : index_range(_test_var_names))
500 : {
501 799 : auto test_var_name = _test_var_names.at(i);
502 799 : _nlfs.Register(test_var_name, std::make_shared<mfem::ParNonlinearForm>(_test_pfespaces.at(i)));
503 : // Apply kernels
504 799 : auto nlf = _nlfs.GetShared(test_var_name);
505 799 : nlf->SetEssentialTrueDofs(_ess_tdof_lists.at(i));
506 799 : ApplyDomainNLFIntegrators(test_var_name, nlf, _kernels_map);
507 799 : ApplyBoundaryNLFIntegrators(test_var_name, nlf, _integrated_bc_map);
508 799 : }
509 775 : }
510 :
511 : void
512 775 : EquationSystem::BuildBilinearForms()
513 : {
514 : // Register bilinear forms
515 1574 : for (const auto i : index_range(_test_var_names))
516 : {
517 799 : auto test_var_name = _test_var_names.at(i);
518 799 : _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
519 :
520 : // Apply kernels
521 799 : auto blf = _blfs.GetShared(test_var_name);
522 799 : blf->SetAssemblyLevel(_assembly_level);
523 1598 : ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
524 799 : test_var_name, test_var_name, blf, _integrated_bc_map);
525 1598 : ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
526 799 : test_var_name, test_var_name, blf, _kernels_map);
527 : // Assemble
528 799 : blf->Assemble();
529 799 : }
530 775 : }
531 :
532 : void
533 775 : 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 1574 : for (const auto i : index_range(_test_var_names))
542 : {
543 799 : auto test_var_name = _test_var_names.at(i);
544 : auto test_mblfs =
545 799 : std::make_shared<Moose::MFEM::NamedFieldsMap<Moose::MFEM::ParMixedBilinearForm>>();
546 1728 : for (const auto j : index_range(_coupled_var_names))
547 : {
548 929 : const auto & coupled_var_name = _coupled_var_names.at(j);
549 1858 : auto mblf = std::make_shared<Moose::MFEM::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
550 929 : _test_pfespaces.at(i));
551 : // Register MixedBilinearForm if kernels exist for it, and assemble kernels
552 929 : if (_kernels_map.Has(test_var_name) &&
553 1834 : _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
554 905 : test_var_name != coupled_var_name)
555 : {
556 119 : mblf->SetAssemblyLevel(_assembly_level);
557 : // Apply all mixed kernels with this test/trial pair
558 238 : ApplyDomainBLFIntegrators<Moose::MFEM::ParMixedBilinearForm>(
559 119 : coupled_var_name, test_var_name, mblf, _kernels_map);
560 : // Assemble mixed bilinear forms
561 119 : mblf->Assemble();
562 : // Register mixed bilinear forms associated with a single trial variable
563 : // for the current test variable
564 119 : test_mblfs->Register(coupled_var_name, mblf);
565 : }
566 929 : }
567 : // Register all mixed bilinear form sets associated with a single test variable
568 799 : _mblfs.Register(test_var_name, test_mblfs);
569 799 : }
570 775 : }
571 :
572 : void
573 1823 : EquationSystem::BuildEquationSystem()
574 : {
575 1823 : BuildBilinearForms();
576 1823 : BuildMixedBilinearForms();
577 1823 : BuildLinearForms();
578 1823 : BuildNonlinearForms();
579 1823 : }
580 :
581 : } // namespace Moose::MFEM
582 :
583 : #endif
|
